PostGIS Raster integration #112

mngr777 · 2022-12-06T15:09:53Z

This new functionality is intended to simplify projecting raster data onto H3 grid.

Two main methods of processing raster data are implemented:

Calculate a centroid for each pixel, determine which H3 cell corresponds to the centroid, then group by H3 index.
ST_Clip the raster by H3 cell boundary for each cell intersecting the raster, process each partial raster using ST_SummaryStats or ST_ValueCount.

The second method becomes more efficient as number of pixels per cell increases. There are also fallback functions for cases when H3 cell is smaller than a pixel, but in practice it would porbably be better to increase raster resolution instead.

Continuous data example:

SELECT
    h3,
    (h3_raster_summary_stats_agg(stats)).*
FROM
    rasters,
    h3_raster_summary(rast, 8)
GROUP BY 1;

       h3        | count | sum  |        mean        |       stddev       | min | max 
-----------------+-------+------+--------------------+--------------------+-----+-----
 88194ad4b7fffff |   200 | 1156 |               5.78 | 2.6384844134464758 |   1 |  10
 88194e6d33fffff |   201 | 1162 |  5.781094527363184 |  2.639508622297715 |   1 |  10
 88194e6d5dfffff |   201 | 1066 |  5.303482587064677 | 2.6263945001831512 |   1 |  10
 <...>

Discrete data example (produces one JSON object per H3 index):

WITH
    summary AS (
        SELECT h3, val, h3_raster_class_summary_item_agg(summary) AS item
        FROM
            test_rasters_bm,
            h3_raster_class_summary(rast, 8)
        GROUP BY 1, 2),
    summary_total AS (
        SELECT h3, val, item, sum((item).count) OVER (PARTITION BY h3) AS total
        FROM summary)
SELECT
    h3,
    jsonb_object_agg(
        concat('class_', val::text),
        h3_raster_class_summary_item_to_jsonb(item) -- val, count, area
            || jsonb_build_object('fraction', (item).count / total)
        ORDER BY val
    ) AS summary
FROM summary_total
GROUP BY 1;

       h3        |                                summary
-----------------+--------------------------------------------------------------------------
 88194e6f3bfffff | {"class_1": {"area": 9913.41706341505, "count": 3, "value": 1, "fraction": 0.058823529411764705}, "class_2": {"area": 13217.889417886734, "count": 4, "value": 2, "fraction": 0.0784313725490196}, <...>}
 88194e6f37fffff | {"class_1": {"area": 46262.61296260357, "count": 14, "value": 1, "fraction": 0.09090909090909091}, "class_2": {"area": 82611.80886179209, "count": 25, "value": 2, "fraction": 0.16233766233766234}, <...>}
<...>

mngr777 · 2023-01-09T22:46:38Z

Some points that may be worth discussing:

Approximating pixel/cell area:
Area is not calculated for each pixel or cell separately.
- Pixel area is always the area of (1, 1) pixel for given raster, h3_raster_class_summary_clip/centroids calculate area by multiplying this value by # of pixels.
- Cell area is approximated as an area of H3 cell at the center of the raster. h3_raster_class_summary_subpixel uses this value for area.
Getting the list of cells intersecting the raster:
polygonToCells returns only cells with centroids inside the polygon. These workarounds seem to work but are not based on some exact calculations.
- Polygon is buffered by h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3
- If raster SRID is not 4326, then the polygon is additionally buffered by max(pixel_height, pixel_width) to account for difference between SRSs.
Selecting the most efficient method:
ST_Clip method is preferred when pixels-per-cell > 250 OR (cells-per-raster > 10000 / (pixels-per-cell - 80)).
This condition was determined experimentally, the graph below shows how much faster (in %) is clip method than centroids method for given pixels-per-cell/cells-per-pixel.
I still don't know why increasing cells-per-raster makes clip method preferable.

mngr777 · 2023-01-10T10:27:20Z

pgxn test fails because tests for released version don't create postgis_raster extension.

mngr777 · 2023-01-10T11:17:15Z

Scripts used to make the graph:
bm.tar.gz

cd batch
./run.sh
./compile.sh > data.csv
./graph.py data.csv

docs/api.md

Komzpa · 2023-01-10T11:42:39Z

docs/api.md

+
+
+## Discrete raster data
+Combining summary from multiple rasters into a single JSON object for each H3 index,


Need a paragraph on why would you want this, and how you define a difference between Discrete/Continuous, and when to pick which.

I've added some examples for types of data and sample output.

Komzpa · 2023-01-10T11:44:31Z

docs/api.md

+Combining summary from multiple rasters into a single JSON object for each H3 index,
+adding `fraction` value (fraction of H3 cell area for each value):
+```
+WITH


This query frightens me as documentation, likely because it lacks sample output to grasp what it does.

Added sample output and a description.

Interested and innocent bystander: It would be good to have an example of this but where the return value is simpler, e.g. just the value at the H3 cell centre, or the return value is the mode value, rather than these complex objects.

@alpha-beta-soup
Thank you, that's a good idea. I've added an example of getting the most frequent value in each cell.
The other example shows an intended use case, and I also tried to show how to calculate a total # of pixels in each cell, since it's not immediately obvious.