PostGIS Raster integration

.github/workflows/pgxn.yml

@@ -24,7 +24,7 @@ jobs:
           sudo systemctl start postgresql.service
           sudo -u postgres -i createuser --superuser runner
           sudo -u postgres -i createdb runner
-          psql -c "CREATE EXTENSION postgis;"
+          psql -c "CREATE EXTENSION postgis; CREATE EXTENSION postgis_raster;"
 
       - name: Bundle
         run: scripts/bundle

.github/workflows/test-linux.yml

@@ -35,7 +35,7 @@ jobs:
           sudo apt-get -y install postgresql-${{ matrix.pg }}-postgis-3
           sudo -u postgres -i createuser --superuser runner
           sudo -u postgres -i createdb runner
-          psql -c "CREATE EXTENSION postgis;"
+          psql -c "CREATE EXTENSION postgis; CREATE EXTENSION postgis_raster;"
 
       - name: Generate
         run: cmake -B build -DCMAKE_BUILD_TYPE=${{ matrix.config }}

docs/api.md

@@ -618,3 +618,166 @@ Create a LinkedGeoPolygon describing the outline(s) of a set of hexagons, conver
 Splits polygons when crossing 180th meridian.
 
 
+# Raster processing functions
+
+## Continuous raster data
+For rasters with pixel values representing continuous data (temperature, humidity,
+elevation), the data inside H3 cells can be summarized by calculating number of
+pixels, sum, mean, standard deviation, min and max for each cell inside a raster
+and grouping these stats across multiple rasters by H3 index.
+```
+SELECT
+    (summary).h3 AS h3,
+    (h3_raster_summary_stats_agg((summary).stats)).*
+FROM (
+    SELECT h3_raster_summary(rast, 8) AS summary
+    FROM rasters
+) t
+GROUP BY 1;
+       h3        | count |        sum         |        mean         |       stddev       |  min  |       max
+-----------------+-------+--------------------+---------------------+--------------------+-------+------------------
+ 882d638189fffff |    10 |  4.607657432556152 | 0.46076574325561526 | 1.3822972297668457 |     0 | 4.607657432556152
+ 882d64c4d1fffff |    10 | 3.6940908953547478 |  0.3694090895354748 |  1.099336879464068 |     0 | 3.667332887649536
+ 882d607431fffff |    11 |  6.219290263950825 |  0.5653900239955295 | 1.7624673707119065 |     0 | 6.13831996917724
+<...>
+```
+
+
+*Since vunreleased*
+
+
+### h3_raster_summary_stats_agg(setof `h3_raster_summary_stats`)
+*Since vunreleased*
+
+
+### h3_raster_summary_clip(rast `raster`, resolution `integer`, [nband `integer` = 1]) ⇒ TABLE (h3 `h3index`, stats `h3_raster_summary_stats`)
+*Since vunreleased*
+
+
+Returns `h3_raster_summary_stats` for each H3 cell in raster for a given band. Clips the raster by H3 cell geometries and processes each part separately.
+
+
+### h3_raster_summary_centroids(rast `raster`, resolution `integer`, [nband `integer` = 1]) ⇒ TABLE (h3 `h3index`, stats `h3_raster_summary_stats`)
+*Since vunreleased*
+
+
+Returns `h3_raster_summary_stats` for each H3 cell in raster for a given band. Finds corresponding H3 cell for each pixel, then groups values by H3 index.
+
+
+### h3_raster_summary_subpixel(rast `raster`, resolution `integer`, [nband `integer` = 1]) ⇒ TABLE (h3 `h3index`, stats `h3_raster_summary_stats`)
+*Since vunreleased*
+
+
+Returns `h3_raster_summary_stats` for each H3 cell in raster for a given band. Assumes H3 cell is smaller than a pixel. Finds corresponding pixel for each H3 cell in raster.
+
+
+### h3_raster_summary(rast `raster`, resolution `integer`, [nband `integer` = 1]) ⇒ TABLE (h3 `h3index`, stats `h3_raster_summary_stats`)
+*Since vunreleased*
+
+
+Returns `h3_raster_summary_stats` for each H3 cell in raster for a given band. Attempts to select an appropriate method based on number of pixels per H3 cell.
+
+
+## Discrete raster data
+For rasters where pixels have discrete values corresponding to different classes
+of land cover or land use, H3 cell data summary can be represented by a JSON object
+with separate fields for each class. First, value, number of pixels and approximate
+area are calculated for each H3 cell and value in a raster, then the stats are
+grouped across multiple rasters by H3 index and value, and after that stats for
+different values in a cell are combined into a single JSON object.
+The following example query additionally calculates a fraction of H3 cell pixels
+for each value, using a window function to get a total number of pixels:
+```
+WITH
+    summary AS (
+        -- get aggregated summary for each H3 index/value pair
+        SELECT h3, val, h3_raster_class_summary_item_agg(summary) AS item
+        FROM
+            rasters,
+            h3_raster_class_summary(rast, 8)
+        GROUP BY 1, 2),
+    summary_total AS (
+        -- add total number of pixels per H3 cell
+        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) -- add fraction value
+        ORDER BY val
+    ) AS summary
+FROM summary_total
+GROUP BY 1;
+      h3        |                                                                            summary
+----------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------
+88194e6f3bfffff | {"class_1": {"area": 75855.5748, "count": 46, "value": 1, "fraction": 0.4509}, "class_2": {"area": 92345.9171, "count": 56, "value": 2, "fraction": 0.5490}}
+88194e6f37fffff | {"class_1": {"area": 255600.3064, "count": 155, "value": 1, "fraction": 0.5}, "class_2": {"area": 255600.3064, "count": 155, "value": 2, "fraction": 0.5}}
+88194e6f33fffff | {"class_1": {"area": 336402.9840, "count": 204, "value": 1, "fraction": 0.5125}, "class_2": {"area": 319912.6416, "count": 194, "value": 2, "fraction": 0.4874}}
+<...>
+```
+Area covered by pixels with the most frequent value in each cell:
+```
+SELECT DISTINCT ON (h3)
+    h3, val, (item).area
+FROM (
+    SELECT
+        h3, val, h3_raster_class_summary_item_agg(summary) AS item
+    FROM
+        rasters,
+        h3_raster_class_summary(rast, 8)
+    GROUP BY 1, 2
+) t
+ORDER BY h3, (item).count DESC;
+       h3        | val |        area
+-----------------+-----+--------------------
+ 88194e6f3bfffff |   5 | 23238.699360251427
+ 88194e6f37fffff |   9 |  60863.26022922993
+ 88194e6f33fffff |   8 |  76355.72646939754
+<...>
+```
+
+
+*Since vunreleased*
+
+
+### h3_raster_class_summary_item_to_jsonb(item `h3_raster_class_summary_item`) ⇒ `jsonb`
+*Since vunreleased*
+
+
+Convert raster summary to JSONB, example: `{"count": 10, "value": 2, "area": 16490.3423}`
+
+
+### h3_raster_class_summary_item_agg(setof `h3_raster_class_summary_item`)
+*Since vunreleased*
+
+
+### h3_raster_class_summary_clip(rast `raster`, resolution `integer`, [nband `integer` = 1]) ⇒ TABLE (h3 `h3index`, val `integer`, summary `h3_raster_class_summary_item`)
+*Since vunreleased*
+
+
+Returns `h3_raster_class_summary_item` for each H3 cell and value for a given band. Clips the raster by H3 cell geometries and processes each part separately.
+
+
+### h3_raster_class_summary_centroids(rast `raster`, resolution `integer`, [nband `integer` = 1]) ⇒ TABLE (h3 `h3index`, val `integer`, summary `h3_raster_class_summary_item`)
+*Since vunreleased*
+
+
+Returns `h3_raster_class_summary_item` for each H3 cell and value for a given band. Finds corresponding H3 cell for each pixel, then groups by H3 and value.
+
+
+### h3_raster_class_summary_subpixel(rast `raster`, resolution `integer`, [nband `integer` = 1]) ⇒ TABLE (h3 `h3index`, val `integer`, summary `h3_raster_class_summary_item`)
+*Since vunreleased*
+
+
+Returns `h3_raster_class_summary_item` for each H3 cell and value for a given band. Assumes H3 cell is smaller than a pixel. Finds corresponding pixel for each H3 cell in raster.
+
+
+### h3_raster_class_summary(rast `raster`, resolution `integer`, [nband `integer` = 1]) ⇒ TABLE (h3 `h3index`, val `integer`, summary `h3_raster_class_summary_item`)
+*Since vunreleased*
+
+
+Returns `h3_raster_class_summary_item` for each H3 cell and value for a given band. Attempts to select an appropriate method based on number of pixels per H3 cell.
+
+

h3_postgis/CMakeLists.txt

@@ -7,6 +7,7 @@ PostgreSQL_add_extension(postgresql_h3_postgis
   REQUIRES
     h3
     postgis
+    postgis_raster
   SOURCES
     src/init.c
     src/wkb_bbox3.c
@@ -23,6 +24,7 @@ PostgreSQL_add_extension(postgresql_h3_postgis
     sql/install/05-regions.sql
     sql/install/20-casts.sql
     sql/install/30-wkb.sql
+    sql/install/40-rasters.sql
     sql/install/99-deprecated.sql
   UPDATES
     sql/updates/h3_postgis--4.0.0.sql