From fa6a6c477bbb0433012ac70f2cbd96b66bdd0373 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Tue, 27 Dec 2022 22:41:27 +0300 Subject: [PATCH 01/22] PostGIS Raster integration --- h3_postgis/CMakeLists.txt | 2 + h3_postgis/sql/install/40-rasters.sql | 576 ++++++++++++++++++ .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 507 +++++++++++++++ h3_postgis/test/CMakeLists.txt | 2 + h3_postgis/test/expected/rasters.out | 172 ++++++ h3_postgis/test/sql/rasters.sql | 176 ++++++ 6 files changed, 1435 insertions(+) create mode 100644 h3_postgis/sql/install/40-rasters.sql create mode 100644 h3_postgis/test/expected/rasters.out create mode 100644 h3_postgis/test/sql/rasters.sql diff --git a/h3_postgis/CMakeLists.txt b/h3_postgis/CMakeLists.txt index d8c168ba..8469c22c 100644 --- a/h3_postgis/CMakeLists.txt +++ b/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 diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql new file mode 100644 index 00000000..04d293f2 --- /dev/null +++ b/h3_postgis/sql/install/40-rasters.sql @@ -0,0 +1,576 @@ +/* + * Copyright 2022 Bytes & Brains + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +--| # Raster processing functions + +-- Get nodata value for ST_Clip function +CREATE OR REPLACE FUNCTION __h3_raster_band_nodata( + rast raster, + nband integer) +RETURNS double precision +AS $$ + SELECT coalesce( + ST_BandNoDataValue(rast, nband), + ST_MinPossibleValue(ST_BandPixelType(rast, nband))); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_pixel_area(rast raster) +RETURNS double precision +AS $$ + SELECT ST_Area(ST_PixelAsPolygon(rast, 1, 1)::geography); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +-- Get area of H3 cell close to the center of the raster +CREATE OR REPLACE FUNCTION __h3_raster_cell_area(rast raster, resolution integer) +RETURNS double precision +AS $$ +DECLARE + rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); +BEGIN + SELECT ST_Area( + h3_cell_to_boundary_geography( + ST_Transform(ST_Centroid(rast_geom), 4326), + resolution)); +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_to_polygon(rast raster, nband integer) +RETURNS geometry +AS $$ + SELECT ST_MinConvexHull(rast, nband); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +-- Get all H3 cells potentially intersecting the polygon, +-- result may contain cells outside of polygon +CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cells( + poly geometry, + resolution integer) +RETURNS SETOF h3index +AS $$ + SELECT h3_polygon_to_cells( + ST_Buffer( + ST_Transform(poly, 4326)::geography, + h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3), + resolution); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +-- Get H3 cell geometries intersecting the raster +CREATE OR REPLACE FUNCTION __h3_raster_to_cell_boundaries( + rast raster, + resolution integer, + nband integer) +RETURNS TABLE(h3 h3index, geom geometry) +AS $$ +DECLARE + rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); +BEGIN + RETURN QUERY + WITH + geoms AS ( + SELECT + c.h3, + ST_Transform(h3_cell_to_boundary_geometry(c.h3), ST_SRID(rast)) AS geom + FROM ( + SELECT __h3_raster_polygon_to_cells(rast_geom, resolution) AS h3 + ) c) + SELECT g.* + FROM geoms g + WHERE ST_Intersects(g.geom, rast_geom); +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + +-- Get H3 cell centroids within the raster +CREATE OR REPLACE FUNCTION __h3_raster_to_cell_centroids( + rast raster, + resolution integer, + nband integer) +RETURNS TABLE(h3 h3index, geom geometry) +AS $$ +DECLARE + rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); +BEGIN + RETURN QUERY + SELECT + c.h3, + ST_Transform(h3_cell_to_geometry(c.h3), ST_SRID(rast)) AS geom + FROM ( + SELECT h3_polygon_to_cells(rast_geom, resolution) AS h3 + ) c; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + +--| ## Continuous raster data +--| +--| Combining summary stats from multiple rasters: +--| ``` +--| 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; +--| ``` + +-- NOTE: `count` can be < 1 when cell area is less than pixel area +--@ availability: unreleased +CREATE TYPE h3_raster_summary_stats AS ( + count double precision, + sum double precision, + mean double precision, + stddev double precision, + min double precision, + max double precision +); + +-- ST_SummaryStats result type to h3_raster_summary_stats +CREATE OR REPLACE FUNCTION __h3_raster_to_summary_stats(stats summarystats) +RETURNS h3_raster_summary_stats +AS $$ + SELECT ROW( + (stats).count, + (stats).sum, + (stats).mean, + (stats).stddev, + (stats).min, + (stats).max) +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_summary_stats_agg_transfn( + s1 h3_raster_summary_stats, + s2 h3_raster_summary_stats) +RETURNS h3_raster_summary_stats +AS $$ + WITH total AS ( + SELECT + (s1).count + (s2).count AS count, + (s1).sum + (s2).sum AS sum) + SELECT ROW( + t.count, + t.sum, + t.sum / t.count, + sqrt( + ( + -- sum of squared values: (variance + mean squared) * count + (((s1).stddev * (s1).stddev + (s1).mean * (s1).mean)) * (s1).count + + (((s2).stddev * (s2).stddev + (s2).mean * (s2).mean)) * (s2).count + ) + / t.count + - ((t.sum * t.sum) / (t.count * t.count)) -- mean squared + ), + least((s1).min, (s2).min), + greatest((s1).max, (s2).max) + ) + FROM total AS t +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +--@ availability: unreleased +CREATE OR REPLACE AGGREGATE h3_raster_summary_stats_agg(h3_raster_summary_stats) ( + sfunc = __h3_raster_summary_stats_agg_transfn, + stype = h3_raster_summary_stats, + parallel = safe +); + +--@ availability: unreleased +CREATE OR REPLACE FUNCTION h3_raster_summary_clip( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ +DECLARE + nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); +BEGIN + RETURN QUERY + WITH parts AS ( + SELECT + g.h3, + ST_Clip(rast, nband, g.geom, nodata) AS part + FROM ( + -- h3, geom + SELECT (__h3_raster_to_cell_boundaries(rast, resolution, nband)).* + ) g) + SELECT + p.h3, + __h3_raster_to_summary_stats(ST_SummaryStats(part, nband, TRUE)) AS stats + FROM parts AS p + WHERE NOT ST_BandIsNoData(part, 1); +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_summary_clip(raster, integer, integer) +IS '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.'; + +--@ availability: unreleased +CREATE OR REPLACE FUNCTION h3_raster_summary_centroids( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ + WITH pixels AS ( + -- x, y, val, geom + SELECT (ST_PixelAsCentroids(rast, nband)).*) + SELECT + h3_lat_lng_to_cell(geom, resolution) AS h3, + ROW( + count(val), + sum(val), + avg(val), + stddev_pop(val), + min(val), + max(val) + ) AS stats + FROM pixels + GROUP BY 1 +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_summary_centroids(raster, integer, integer) +IS '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.'; + +CREATE OR REPLACE FUNCTION __h3_raster_summary_subpixel( + rast raster, + resolution integer, + nband integer, + pixel_area double precision, + cell_area double precision) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ + WITH + vals AS ( + SELECT + h3, + ST_Value( + rast, + nband, + ST_WorldToRasterCoordX(rast, geom), + ST_WorldToRasterCoordY(rast, geom) + ) AS val + FROM ( + -- h3, geom + SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* + ) t) + SELECT + h3, + ROW( + cell_area / pixel_area, -- count + val, -- sum + val, -- mean + 0, -- stddev + val, -- min + val -- max + ) + FROM vals; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +--@ availability: unreleased +CREATE OR REPLACE FUNCTION h3_raster_summary_subpixel( + rast raster, + resolution integer, + nband integer DEFAUlT 1) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ +DECLARE + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); +BEGIN + RETURN QUERY SELECT (__h3_raster_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_summary_subpixel(raster, integer, integer) +IS '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.'; + +CREATE OR REPLACE FUNCTION h3_raster_summary( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ +DECLARE + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + pixels_per_cell CONSTANT double precision := cell_area / pixel_area; +BEGIN + IF pixels_per_cell > 350 THEN + RETURN QUERY SELECT (h3_raster_summary_clip(rast, resolution, nband)).*; + ELSIF pixels_per_cell > 1 THEN + RETURN QUERY SELECT (h3_raster_summary_centroids(rast, resolution, nband)).*; + ELSE + RETURN QUERY SELECT (__h3_raster_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + END IF; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_summary(raster, integer, integer) +IS '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 +--| +--| 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 +--| summary AS ( +--| SELECT h3, val, h3_raster_class_summary_item_agg(summary) AS item +--| FROM ( +--| -- h3, val, summary +--| SELECT (h3_raster_class_summary(rast, 8)).* AS summary +--| FROM rasters +--| ) t +--| 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; +--| ``` + +-- NOTE: `count` can be < 1 when cell area is less than pixel area +--@ availability: unreleased +CREATE TYPE h3_raster_class_summary_item AS ( + val integer, + count double precision, + area double precision +); + +--@ availability: unreleased +CREATE OR REPLACE FUNCTION h3_raster_class_summary_item_to_jsonb( + item h3_raster_class_summary_item) +RETURNS jsonb +AS $$ + SELECT jsonb_build_object( + 'value', (item).val, + 'count', (item).count, + 'area', (item).area + ); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_class_summary_item_to_jsonb(h3_raster_class_summary_item) +IS 'Convert raster summary to binary JSON.'; + +CREATE OR REPLACE FUNCTION __h3_raster_class_summary_item_agg_transfn( + s1 h3_raster_class_summary_item, + s2 h3_raster_class_summary_item) +RETURNS h3_raster_class_summary_item +AS $$ + SELECT + s1.val, + s1.count + s2.count, + s1.area + s2.area +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +--@ availability: unreleased +CREATE OR REPLACE AGGREGATE h3_raster_class_summary_item_agg(h3_raster_class_summary_item) ( + stype = h3_raster_class_summary_item, + sfunc = __h3_raster_class_summary_item_agg_transfn, + parallel = safe +); + +-- Get summary items for a raster clipped by H3 cell geometry +CREATE OR REPLACE FUNCTION __h3_raster_class_summary_part( + rast raster, + nband integer, + pixel_area double precision) +RETURNS SETOF h3_raster_class_summary_item +AS $$ + WITH + vals AS (SELECT unnest(ST_DumpValues(rast, nband)) AS val) + SELECT + vals.val::integer, + count(*)::double precision, + count(*) * pixel_area + FROM vals + WHERE val IS NOT NULL + GROUP BY 1 +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_class_summary_clip( + rast raster, + resolution integer, + nband integer, + pixel_area double precision) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ +DECLARE + nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); +BEGIN + RETURN QUERY + WITH + parts AS ( + SELECT + g.h3, + ST_Clip(rast, nband, g.geom, nodata, TRUE) AS part + FROM ( + -- h3, geom + SELECT (__h3_raster_to_cell_boundaries(rast, resolution, nband)).* + ) g), + summary AS ( + SELECT + p.h3, + __h3_raster_class_summary_part(part, nband, pixel_area) AS summary + FROM parts p) + SELECT s.h3, (s.summary).val, s.summary + FROM summary s; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + +--@ availability: unreleased +CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ +DECLARE + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); +BEGIN + RETURN QUERY SELECT (__h3_raster_class_summary_clip(rast, resolution, nband, pixel_area)).*; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_class_summary_clip(raster, integer, integer) +IS '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.'; + + +CREATE OR REPLACE FUNCTION __h3_raster_class_summary_centroids( + rast raster, + resolution integer, + nband integer, + pixel_area double precision) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ + SELECT + h3_lat_lng_to_cell(geom, resolution) AS h3, + val::integer, + ROW( + val::integer, + count(*)::double precision, + count(*) * pixel_area) + FROM ( + -- x, y, val, geom + SELECT (ST_PixelAsCentroids(rast, nband)).* + ) c + GROUP BY 1, 2; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + + +-- For each pixel determine which H3 cell it belongs to then group by H3 index and value. +--@ availability: unreleased +CREATE OR REPLACE FUNCTION h3_raster_class_summary_centroids( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ +DECLARE + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); +BEGIN + RETURN QUERY SELECT (__h3_raster_class_summary_centroids(rast, resolution, nband, pixel_area)).*; +END +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_class_summary_centroids(raster, integer, integer) +IS '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.'; + +--@ availability: unreleased +CREATE OR REPLACE FUNCTION __h3_raster_class_summary_subpixel( + rast raster, + resolution integer, + nband integer, + pixel_area double precision, + cell_area double precision) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ + WITH + vals AS ( + SELECT + h3, + ST_Value( + rast, + nband, + ST_WorldToRasterCoordX(rast, geom), + ST_WorldToRasterCoordY(rast, geom) + ) AS val + FROM ( + -- h3, geom + SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* + ) c) + SELECT + h3, + val::integer, + ROW( + val::integer, + cell_area / pixel_area, + cell_area) + FROM vals v; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +-- Get summary items for each H3 index and value. +-- For each H3 cell centroid determine which pixel it belongs to. +--@ availability: unreleased +CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ +DECLARE + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); +BEGIN + RETURN QUERY SELECT (__h3_raster_class_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_class_summary_subpixel(raster, integer, integer) +IS '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.'; + + +-- Get summary items for each H3 index and value. +-- Select appropriate method based on number of pixels per H3 cell. +--@ availability: unreleased +CREATE OR REPLACE FUNCTION h3_raster_class_summary( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ +DECLARE + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + pixels_per_cell CONSTANT double precision := cell_area / pixel_area; +BEGIN + IF pixels_per_cell > 350 THEN + RETURN QUERY SELECT (__h3_raster_class_summary_clip(rast, resolution, nband, pixel_area)).*; + ELSIF pixels_per_cell > 1 THEN + RETURN QUERY SELECT (__h3_raster_class_summary_centroids(rast, resolution, nband, pixel_area)).*; + ELSE + RETURN QUERY SELECT (__h3_raster_class_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + END IF; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION h3_raster_class_summary(raster, integer, integer) +IS '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.'; diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index 2109e42d..7b03be1a 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -34,3 +34,510 @@ AS 'h3' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION IS 'Create a LinkedGeoPolygon describing the outline(s) of a set of hexagons, converts to EWKB. Splits polygons when crossing 180th meridian.'; + +-- Raster processing + +-- Get nodata value for ST_Clip function +CREATE OR REPLACE FUNCTION __h3_raster_band_nodata( + rast raster, + nband integer) +RETURNS double precision +AS $$ + SELECT coalesce( + ST_BandNoDataValue(rast, nband), + ST_MinPossibleValue(ST_BandPixelType(rast, nband))); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_pixel_area(rast raster) +RETURNS double precision +AS $$ + SELECT ST_Area(ST_PixelAsPolygon(rast, 1, 1)::geography); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +-- Get area of H3 cell close to the center of the raster +CREATE OR REPLACE FUNCTION __h3_raster_cell_area(rast raster, resolution integer) +RETURNS double precision +AS $$ +DECLARE + rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); +BEGIN + SELECT ST_Area( + h3_cell_to_boundary_geography( + ST_Transform(ST_Centroid(rast_geom), 4326), + resolution)); +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_to_polygon(rast raster, nband integer) +RETURNS geometry +AS $$ + SELECT ST_MinConvexHull(rast, nband); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +-- Get all H3 cells potentially intersecting the polygon, +-- result may contain cells outside of polygon +CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cells( + poly geometry, + resolution integer) +RETURNS SETOF h3index +AS $$ + SELECT h3_polygon_to_cells( + ST_Buffer( + ST_Transform(poly, 4326)::geography, + h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3), + resolution); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +-- Get H3 cell geometries intersecting the raster +CREATE OR REPLACE FUNCTION __h3_raster_to_cell_boundaries( + rast raster, + resolution integer, + nband integer) +RETURNS TABLE(h3 h3index, geom geometry) +AS $$ +DECLARE + rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); +BEGIN + RETURN QUERY + WITH + geoms AS ( + SELECT + c.h3, + ST_Transform(h3_cell_to_boundary_geometry(c.h3), ST_SRID(rast)) AS geom + FROM ( + SELECT __h3_raster_polygon_to_cells(rast_geom, resolution) AS h3 + ) c) + SELECT g.* + FROM geoms g + WHERE ST_Intersects(g.geom, rast_geom); +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + +-- Get H3 cell centroids within the raster +CREATE OR REPLACE FUNCTION __h3_raster_to_cell_centroids( + rast raster, + resolution integer, + nband integer) +RETURNS TABLE(h3 h3index, geom geometry) +AS $$ +DECLARE + rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); +BEGIN + RETURN QUERY + SELECT + c.h3, + ST_Transform(h3_cell_to_geometry(c.h3), ST_SRID(rast)) AS geom + FROM ( + SELECT h3_polygon_to_cells(rast_geom, resolution) AS h3 + ) c; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + +-- Raster processing: continuous data + +CREATE TYPE h3_raster_summary_stats AS ( + count double precision, + sum double precision, + mean double precision, + stddev double precision, + min double precision, + max double precision +); + +-- ST_SummaryStats result type to h3_raster_summary_stats +CREATE OR REPLACE FUNCTION __h3_raster_to_summary_stats(stats summarystats) +RETURNS h3_raster_summary_stats +AS $$ + SELECT ROW( + (stats).count, + (stats).sum, + (stats).mean, + (stats).stddev, + (stats).min, + (stats).max) +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_summary_stats_agg_transfn( + s1 h3_raster_summary_stats, + s2 h3_raster_summary_stats) +RETURNS h3_raster_summary_stats +AS $$ + WITH total AS ( + SELECT + (s1).count + (s2).count AS count, + (s1).sum + (s2).sum AS sum) + SELECT ROW( + t.count, + t.sum, + t.sum / t.count, + sqrt( + ( + -- sum of squared values: (variance + mean squared) * count + (((s1).stddev * (s1).stddev + (s1).mean * (s1).mean)) * (s1).count + + (((s2).stddev * (s2).stddev + (s2).mean * (s2).mean)) * (s2).count + ) + / t.count + - ((t.sum * t.sum) / (t.count * t.count)) -- mean squared + ), + least((s1).min, (s2).min), + greatest((s1).max, (s2).max) + ) + FROM total AS t +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE AGGREGATE h3_raster_summary_stats_agg(h3_raster_summary_stats) ( + sfunc = __h3_raster_summary_stats_agg_transfn, + stype = h3_raster_summary_stats, + parallel = safe +); + +CREATE OR REPLACE FUNCTION h3_raster_summary_clip( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ +DECLARE + nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); +BEGIN + RETURN QUERY + WITH parts AS ( + SELECT + g.h3, + ST_Clip(rast, nband, g.geom, nodata) AS part + FROM ( + -- h3, geom + SELECT (__h3_raster_to_cell_boundaries(rast, resolution, nband)).* + ) g) + SELECT + p.h3, + __h3_raster_to_summary_stats(ST_SummaryStats(part, nband, TRUE)) AS stats + FROM parts AS p + WHERE NOT ST_BandIsNoData(part, 1); +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_summary_clip(raster, integer, integer) +IS '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.'; + +CREATE OR REPLACE FUNCTION h3_raster_summary_centroids( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ + WITH pixels AS ( + -- x, y, val, geom + SELECT (ST_PixelAsCentroids(rast, nband)).*) + SELECT + h3_lat_lng_to_cell(geom, resolution) AS h3, + ROW( + count(val), + sum(val), + avg(val), + stddev_pop(val), + min(val), + max(val) + ) AS stats + FROM pixels + GROUP BY 1 +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_summary_centroids(raster, integer, integer) +IS '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.'; + +CREATE OR REPLACE FUNCTION __h3_raster_summary_subpixel( + rast raster, + resolution integer, + nband integer, + pixel_area double precision, + cell_area double precision) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ + WITH + vals AS ( + SELECT + h3, + ST_Value( + rast, + nband, + ST_WorldToRasterCoordX(rast, geom), + ST_WorldToRasterCoordY(rast, geom) + ) AS val + FROM ( + -- h3, geom + SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* + ) t) + SELECT + h3, + ROW( + cell_area / pixel_area, -- count + val, -- sum + val, -- mean + 0, -- stddev + val, -- min + val -- max + ) + FROM vals; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION h3_raster_summary_subpixel( + rast raster, + resolution integer, + nband integer DEFAUlT 1) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ +DECLARE + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); +BEGIN + RETURN QUERY SELECT (__h3_raster_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_summary_subpixel(raster, integer, integer) +IS '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.'; + +CREATE OR REPLACE FUNCTION h3_raster_summary( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ +DECLARE + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + pixels_per_cell CONSTANT double precision := cell_area / pixel_area; +BEGIN + IF pixels_per_cell > 350 THEN + RETURN QUERY SELECT (h3_raster_summary_clip(rast, resolution, nband)).*; + ELSIF pixels_per_cell > 1 THEN + RETURN QUERY SELECT (h3_raster_summary_centroids(rast, resolution, nband)).*; + ELSE + RETURN QUERY SELECT (__h3_raster_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + END IF; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_summary(raster, integer, integer) +IS '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.'; + +-- Raster processing: discrete data + +CREATE TYPE h3_raster_class_summary_item AS ( + val integer, + count double precision, + area double precision +); + +CREATE OR REPLACE FUNCTION h3_raster_class_summary_item_to_jsonb( + item h3_raster_class_summary_item) +RETURNS jsonb +AS $$ + SELECT jsonb_build_object( + 'value', (item).val, + 'count', (item).count, + 'area', (item).area + ); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_class_summary_item_to_jsonb(h3_raster_class_summary_item) +IS 'Convert raster summary to binary JSON.'; + +CREATE OR REPLACE FUNCTION __h3_raster_class_summary_item_agg_transfn( + s1 h3_raster_class_summary_item, + s2 h3_raster_class_summary_item) +RETURNS h3_raster_class_summary_item +AS $$ + SELECT + s1.val, + s1.count + s2.count, + s1.area + s2.area +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE AGGREGATE h3_raster_class_summary_item_agg(h3_raster_class_summary_item) ( + stype = h3_raster_class_summary_item, + sfunc = __h3_raster_class_summary_item_agg_transfn, + parallel = safe +); + +-- Get summary items for a raster clipped by H3 cell geometry +CREATE OR REPLACE FUNCTION __h3_raster_class_summary_part( + rast raster, + nband integer, + pixel_area double precision) +RETURNS SETOF h3_raster_class_summary_item +AS $$ + WITH + vals AS (SELECT unnest(ST_DumpValues(rast, nband)) AS val) + SELECT + vals.val::integer, + count(*)::double precision, + count(*) * pixel_area + FROM vals + WHERE val IS NOT NULL + GROUP BY 1 +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_class_summary_clip( + rast raster, + resolution integer, + nband integer, + pixel_area double precision) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ +DECLARE + nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); +BEGIN + RETURN QUERY + WITH + parts AS ( + SELECT + g.h3, + ST_Clip(rast, nband, g.geom, nodata, TRUE) AS part + FROM ( + -- h3, geom + SELECT (__h3_raster_to_cell_boundaries(rast, resolution, nband)).* + ) g), + summary AS ( + SELECT + p.h3, + __h3_raster_class_summary_part(part, nband, pixel_area) AS summary + FROM parts p) + SELECT s.h3, (s.summary).val, s.summary + FROM summary s; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ +DECLARE + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); +BEGIN + RETURN QUERY SELECT (__h3_raster_class_summary_clip(rast, resolution, nband, pixel_area)).*; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_class_summary_clip(raster, integer, integer) +IS '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.'; + + +CREATE OR REPLACE FUNCTION __h3_raster_class_summary_centroids( + rast raster, + resolution integer, + nband integer, + pixel_area double precision) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ + SELECT + h3_lat_lng_to_cell(geom, resolution) AS h3, + val::integer, + ROW( + val::integer, + count(*)::double precision, + count(*) * pixel_area) + FROM ( + -- x, y, val, geom + SELECT (ST_PixelAsCentroids(rast, nband)).* + ) c + GROUP BY 1, 2; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + + +-- For each pixel determine which H3 cell it belongs to then group by H3 index and value. +CREATE OR REPLACE FUNCTION h3_raster_class_summary_centroids( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ +DECLARE + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); +BEGIN + RETURN QUERY SELECT (__h3_raster_class_summary_centroids(rast, resolution, nband, pixel_area)).*; +END +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_class_summary_centroids(raster, integer, integer) +IS '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.'; + +CREATE OR REPLACE FUNCTION __h3_raster_class_summary_subpixel( + rast raster, + resolution integer, + nband integer, + pixel_area double precision, + cell_area double precision) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ + WITH + vals AS ( + SELECT + h3, + ST_Value( + rast, + nband, + ST_WorldToRasterCoordX(rast, geom), + ST_WorldToRasterCoordY(rast, geom) + ) AS val + FROM ( + -- h3, geom + SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* + ) c) + SELECT + h3, + val::integer, + ROW( + val::integer, + cell_area / pixel_area, + cell_area) + FROM vals v; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +-- Get summary items for each H3 index and value. +-- For each H3 cell centroid determine which pixel it belongs to. +CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ +DECLARE + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); +BEGIN + RETURN QUERY SELECT (__h3_raster_class_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION + h3_raster_class_summary_subpixel(raster, integer, integer) +IS '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.'; + + +-- Get summary items for each H3 index and value. +-- Select appropriate method based on number of pixels per H3 cell. +CREATE OR REPLACE FUNCTION h3_raster_class_summary( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) +AS $$ +DECLARE + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + pixels_per_cell CONSTANT double precision := cell_area / pixel_area; +BEGIN + IF pixels_per_cell > 350 THEN + RETURN QUERY SELECT (__h3_raster_class_summary_clip(rast, resolution, nband, pixel_area)).*; + ELSIF pixels_per_cell > 1 THEN + RETURN QUERY SELECT (__h3_raster_class_summary_centroids(rast, resolution, nband, pixel_area)).*; + ELSE + RETURN QUERY SELECT (__h3_raster_class_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + END IF; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION h3_raster_class_summary(raster, integer, integer) +IS '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.'; diff --git a/h3_postgis/test/CMakeLists.txt b/h3_postgis/test/CMakeLists.txt index 5a3a26b8..48cec9ee 100644 --- a/h3_postgis/test/CMakeLists.txt +++ b/h3_postgis/test/CMakeLists.txt @@ -1,5 +1,6 @@ set(TESTS postgis + rasters ) if(PostgreSQL_REGRESS) @@ -12,6 +13,7 @@ if(PostgreSQL_REGRESS) --outputdir=${CMAKE_CURRENT_BINARY_DIR} --load-extension h3 --load-extension postgis + --load-extension postgis_raster --load-extension h3_postgis ${TESTS} ) diff --git a/h3_postgis/test/expected/rasters.out b/h3_postgis/test/expected/rasters.out new file mode 100644 index 00000000..4314c1a1 --- /dev/null +++ b/h3_postgis/test/expected/rasters.out @@ -0,0 +1,172 @@ +\pset tuples_only on +\set resolution 9 +\set coverage_size 2 +\set raster_size 25 +\set pixel_size 0.0005 +\set value_num 5 +\set lat 51.5 +\set lng -0.025 +CREATE TABLE h3_test_rasters (id SERIAL, rast raster); +INSERT INTO h3_test_rasters (rast) ( + WITH + vals AS ( + SELECT array_agg(row) AS vals + FROM ( + SELECT array_agg((x + y) % :value_num + 1) AS row + FROM + generate_series(1, :raster_size) AS x, + generate_series(1, :raster_size) AS y + GROUP BY y + ) t), + rasts AS ( + SELECT + ST_AddBand( + ST_MakeEmptyCoverage( + :raster_size, :raster_size, + :raster_size * :coverage_size, :raster_size * :coverage_size, + :lng, :lat, + :pixel_size, -(:pixel_size), + 0, 0, + 4326), + ARRAY[ROW(1, '8BUI', 1, 0)]::addbandarg[] + ) AS rast) + SELECT ST_SetValues(r.rast, 1, 1, 1, v.vals) + FROM rasts r, vals v +); +CREATE FUNCTION h3_test_equal( + v1 double precision, + v2 double precision) +RETURNS boolean +AS $$ + SELECT ABS(v1 - v2) < 1e-12; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; +CREATE FUNCTION h3_test_raster_summary_stats_equal( + s1 h3_raster_summary_stats, + s2 h3_raster_summary_stats) +RETURNS boolean +AS $$ + SELECT s1 IS NOT NULL AND s2 IS NOT NULL + AND h3_test_equal((s1).count, (s2).count) + AND h3_test_equal((s1).sum, (s2).sum) + AND h3_test_equal((s1).mean, (s2).mean) + AND h3_test_equal((s1).stddev, (s2).stddev) + AND h3_test_equal((s1).min, (s2).min) + AND h3_test_equal((s1).max, (s2).max); +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; +CREATE FUNCTION h3_test_raster_class_summary_item_equal( + i1 h3_raster_class_summary_item, + i2 h3_raster_class_summary_item) +RETURNS boolean +AS $$ + SELECT i1 IS NOT NULL AND i2 IS NOT NULL + AND h3_test_equal((i1).val, (i2).val) + AND h3_test_equal((i1).count, (i2).count) + AND h3_test_equal((i1).area, (i2).area); +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; +-- Results of `h3_raster_summary_clip` and `h3_raster_summary_centroids` +-- should be identical +WITH + clip AS ( + SELECT + h3, + h3_raster_summary_stats_agg(stats) AS stats + FROM ( + -- h3, stats + SELECT (h3_raster_summary_clip(rast, :resolution)).* + FROM h3_test_rasters + ) t + GROUP BY 1), + centroids AS ( + SELECT + h3, + h3_raster_summary_stats_agg(stats) AS stats + FROM ( + -- h3, stats + SELECT (h3_raster_summary_centroids(rast, :resolution)).* + FROM h3_test_rasters + ) t + GROUP BY 1) +SELECT COUNT(*) +FROM clip a FULL OUTER JOIN centroids b ON a.h3 = b.h3 +WHERE NOT h3_test_raster_summary_stats_equal(a.stats, b.stats); + 0 + +-- Results of `h3_raster_class_summary_clip` and `h3_raster_class_summary_clip` +-- should be identical +WITH + clip AS ( + SELECT + h3, + val, + h3_raster_class_summary_item_agg(summary) AS summary + FROM ( + -- h3, val, summary + SELECT (h3_raster_class_summary_clip(rast, :resolution)).* + FROM h3_test_rasters + ) t + GROUP BY 1, 2), + centroids AS ( + SELECT + h3, + val, + h3_raster_class_summary_item_agg(summary) AS summary + FROM ( + -- h3, val, summary + SELECT (h3_raster_class_summary_centroids(rast, :resolution)).* + FROM h3_test_rasters + ) t + GROUP BY 1, 2) +SELECT COUNT(*) +FROM clip a FULL OUTER JOIN centroids b ON a.h3 = b.h3 AND a.val = b.val +WHERE NOT h3_test_raster_class_summary_item_equal(a.summary, b.summary); + 0 + +-- Stats aggregation check: +-- stats for a cell intersecting multiple rasters (with aggregation) should be +-- the same when calculated on a union of rasters (without aggregation). +WITH + rast AS ( + -- Union all test rasters + SELECT ST_Union(rast) AS rast FROM h3_test_rasters), + middle AS ( + -- Find an H3 cell in a bottom-right corner of a first raster + -- (intersecting 4 rasters) + SELECT + h3_lat_lng_to_cell( + ST_MakePoint( + ST_RasterToWorldCoordX(rast, :raster_size), + ST_RasterToWorldCoordY(rast, :raster_size)), + :resolution + ) AS h3 + FROM rast), + summary1 AS ( + -- Get summary from combined raster + SELECT t.stats + FROM ( + -- h3, stats + SELECT (h3_raster_summary_clip(rast, :resolution)).* + FROM rast + ) t, middle m + WHERE t.h3 = m.h3), + summary2 AS ( + -- Get aggregates summary from separate rasters + SELECT h3_raster_summary_stats_agg(t.stats) AS stats + FROM ( + -- h3, stats + SELECT (h3_raster_summary_clip(rast, :resolution)).* + FROM h3_test_rasters + ) t, middle m + WHERE t.h3 = m.h3 + GROUP BY t.h3) +SELECT h3_test_raster_summary_stats_equal(s1.stats, s2.stats) +FROM summary1 s1, summary2 s2; + t + +DROP FUNCTION h3_test_raster_class_summary_item_equal( + h3_raster_class_summary_item, + h3_raster_class_summary_item); +DROP FUNCTION h3_test_raster_summary_stats_equal( + h3_raster_summary_stats, + h3_raster_summary_stats); +DROP FUNCTION h3_test_equal(double precision, double precision); +DROP TABLE h3_test_rasters; diff --git a/h3_postgis/test/sql/rasters.sql b/h3_postgis/test/sql/rasters.sql new file mode 100644 index 00000000..b79601a0 --- /dev/null +++ b/h3_postgis/test/sql/rasters.sql @@ -0,0 +1,176 @@ +\pset tuples_only on +\set resolution 9 +\set coverage_size 2 +\set raster_size 25 +\set pixel_size 0.0005 +\set value_num 5 +\set lat 51.5 +\set lng -0.025 + +CREATE TABLE h3_test_rasters (id SERIAL, rast raster); + +INSERT INTO h3_test_rasters (rast) ( + WITH + vals AS ( + SELECT array_agg(row) AS vals + FROM ( + SELECT array_agg((x + y) % :value_num + 1) AS row + FROM + generate_series(1, :raster_size) AS x, + generate_series(1, :raster_size) AS y + GROUP BY y + ) t), + rasts AS ( + SELECT + ST_AddBand( + ST_MakeEmptyCoverage( + :raster_size, :raster_size, + :raster_size * :coverage_size, :raster_size * :coverage_size, + :lng, :lat, + :pixel_size, -(:pixel_size), + 0, 0, + 4326), + ARRAY[ROW(1, '8BUI', 1, 0)]::addbandarg[] + ) AS rast) + SELECT ST_SetValues(r.rast, 1, 1, 1, v.vals) + FROM rasts r, vals v +); + +CREATE FUNCTION h3_test_equal( + v1 double precision, + v2 double precision) +RETURNS boolean +AS $$ + SELECT ABS(v1 - v2) < 1e-12; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; + +CREATE FUNCTION h3_test_raster_summary_stats_equal( + s1 h3_raster_summary_stats, + s2 h3_raster_summary_stats) +RETURNS boolean +AS $$ + SELECT s1 IS NOT NULL AND s2 IS NOT NULL + AND h3_test_equal((s1).count, (s2).count) + AND h3_test_equal((s1).sum, (s2).sum) + AND h3_test_equal((s1).mean, (s2).mean) + AND h3_test_equal((s1).stddev, (s2).stddev) + AND h3_test_equal((s1).min, (s2).min) + AND h3_test_equal((s1).max, (s2).max); +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; + +CREATE FUNCTION h3_test_raster_class_summary_item_equal( + i1 h3_raster_class_summary_item, + i2 h3_raster_class_summary_item) +RETURNS boolean +AS $$ + SELECT i1 IS NOT NULL AND i2 IS NOT NULL + AND h3_test_equal((i1).val, (i2).val) + AND h3_test_equal((i1).count, (i2).count) + AND h3_test_equal((i1).area, (i2).area); +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; + +-- Results of `h3_raster_summary_clip` and `h3_raster_summary_centroids` +-- should be identical +WITH + clip AS ( + SELECT + h3, + h3_raster_summary_stats_agg(stats) AS stats + FROM ( + -- h3, stats + SELECT (h3_raster_summary_clip(rast, :resolution)).* + FROM h3_test_rasters + ) t + GROUP BY 1), + centroids AS ( + SELECT + h3, + h3_raster_summary_stats_agg(stats) AS stats + FROM ( + -- h3, stats + SELECT (h3_raster_summary_centroids(rast, :resolution)).* + FROM h3_test_rasters + ) t + GROUP BY 1) +SELECT COUNT(*) +FROM clip a FULL OUTER JOIN centroids b ON a.h3 = b.h3 +WHERE NOT h3_test_raster_summary_stats_equal(a.stats, b.stats); + +-- Results of `h3_raster_class_summary_clip` and `h3_raster_class_summary_clip` +-- should be identical +WITH + clip AS ( + SELECT + h3, + val, + h3_raster_class_summary_item_agg(summary) AS summary + FROM ( + -- h3, val, summary + SELECT (h3_raster_class_summary_clip(rast, :resolution)).* + FROM h3_test_rasters + ) t + GROUP BY 1, 2), + centroids AS ( + SELECT + h3, + val, + h3_raster_class_summary_item_agg(summary) AS summary + FROM ( + -- h3, val, summary + SELECT (h3_raster_class_summary_centroids(rast, :resolution)).* + FROM h3_test_rasters + ) t + GROUP BY 1, 2) +SELECT COUNT(*) +FROM clip a FULL OUTER JOIN centroids b ON a.h3 = b.h3 AND a.val = b.val +WHERE NOT h3_test_raster_class_summary_item_equal(a.summary, b.summary); + +-- Stats aggregation check: +-- stats for a cell intersecting multiple rasters (with aggregation) should be +-- the same when calculated on a union of rasters (without aggregation). +WITH + rast AS ( + -- Union all test rasters + SELECT ST_Union(rast) AS rast FROM h3_test_rasters), + middle AS ( + -- Find an H3 cell in a bottom-right corner of a first raster + -- (intersecting 4 rasters) + SELECT + h3_lat_lng_to_cell( + ST_MakePoint( + ST_RasterToWorldCoordX(rast, :raster_size), + ST_RasterToWorldCoordY(rast, :raster_size)), + :resolution + ) AS h3 + FROM rast), + summary1 AS ( + -- Get summary from combined raster + SELECT t.stats + FROM ( + -- h3, stats + SELECT (h3_raster_summary_clip(rast, :resolution)).* + FROM rast + ) t, middle m + WHERE t.h3 = m.h3), + summary2 AS ( + -- Get aggregates summary from separate rasters + SELECT h3_raster_summary_stats_agg(t.stats) AS stats + FROM ( + -- h3, stats + SELECT (h3_raster_summary_clip(rast, :resolution)).* + FROM h3_test_rasters + ) t, middle m + WHERE t.h3 = m.h3 + GROUP BY t.h3) +SELECT h3_test_raster_summary_stats_equal(s1.stats, s2.stats) +FROM summary1 s1, summary2 s2; + +DROP FUNCTION h3_test_raster_class_summary_item_equal( + h3_raster_class_summary_item, + h3_raster_class_summary_item); +DROP FUNCTION h3_test_raster_summary_stats_equal( + h3_raster_summary_stats, + h3_raster_summary_stats); +DROP FUNCTION h3_test_equal(double precision, double precision); + +DROP TABLE h3_test_rasters; From 48a65303845dac168492f5bbff1fe23b6c7cb596 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Wed, 28 Dec 2022 15:02:24 +0300 Subject: [PATCH 02/22] PostGIS Raster integration: creating postgis_raster extension in tests --- .github/workflows/pgxn.yml | 2 +- .github/workflows/test-linux.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pgxn.yml b/.github/workflows/pgxn.yml index fc24f733..aa1c5ea6 100644 --- a/.github/workflows/pgxn.yml +++ b/.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 diff --git a/.github/workflows/test-linux.yml b/.github/workflows/test-linux.yml index 67cdb946..ada69b8d 100644 --- a/.github/workflows/test-linux.yml +++ b/.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 }} From a86913434f35aaea27638d07f3bec8e2a13a5731 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Wed, 28 Dec 2022 15:35:35 +0300 Subject: [PATCH 03/22] PostGIS Raster integration: explicit conversion of records to return type --- h3_postgis/sql/install/40-rasters.sql | 23 +++++++++++-------- .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 23 +++++++++++-------- 2 files changed, 26 insertions(+), 20 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 04d293f2..7070851c 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -72,7 +72,7 @@ CREATE OR REPLACE FUNCTION __h3_raster_to_cell_boundaries( rast raster, resolution integer, nband integer) -RETURNS TABLE(h3 h3index, geom geometry) +RETURNS TABLE (h3 h3index, geom geometry) AS $$ DECLARE rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); @@ -97,7 +97,7 @@ CREATE OR REPLACE FUNCTION __h3_raster_to_cell_centroids( rast raster, resolution integer, nband integer) -RETURNS TABLE(h3 h3index, geom geometry) +RETURNS TABLE (h3 h3index, geom geometry) AS $$ DECLARE rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); @@ -147,7 +147,8 @@ AS $$ (stats).mean, (stats).stddev, (stats).min, - (stats).max) + (stats).max + )::h3_raster_summary_stats $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_summary_stats_agg_transfn( @@ -174,7 +175,7 @@ AS $$ ), least((s1).min, (s2).min), greatest((s1).max, (s2).max) - ) + )::h3_raster_summary_stats FROM total AS t $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -234,7 +235,7 @@ AS $$ stddev_pop(val), min(val), max(val) - ) AS stats + )::h3_raster_summary_stats AS stats FROM pixels GROUP BY 1 $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -273,7 +274,7 @@ AS $$ 0, -- stddev val, -- min val -- max - ) + )::h3_rasters_summary_stats AS stats FROM vals; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -464,11 +465,12 @@ RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ SELECT h3_lat_lng_to_cell(geom, resolution) AS h3, - val::integer, + val::integer AS val, ROW( val::integer, count(*)::double precision, - count(*) * pixel_area) + count(*) * pixel_area + )::h3_raster_class_summary_item AS summary FROM ( -- x, y, val, geom SELECT (ST_PixelAsCentroids(rast, nband)).* @@ -520,11 +522,12 @@ AS $$ ) c) SELECT h3, - val::integer, + val::integer AS val, ROW( val::integer, cell_area / pixel_area, - cell_area) + cell_area + )::h3_raster_class_summary_item AS summary FROM vals v; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index 7b03be1a..25c6324d 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -93,7 +93,7 @@ CREATE OR REPLACE FUNCTION __h3_raster_to_cell_boundaries( rast raster, resolution integer, nband integer) -RETURNS TABLE(h3 h3index, geom geometry) +RETURNS TABLE (h3 h3index, geom geometry) AS $$ DECLARE rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); @@ -118,7 +118,7 @@ CREATE OR REPLACE FUNCTION __h3_raster_to_cell_centroids( rast raster, resolution integer, nband integer) -RETURNS TABLE(h3 h3index, geom geometry) +RETURNS TABLE (h3 h3index, geom geometry) AS $$ DECLARE rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); @@ -154,7 +154,8 @@ AS $$ (stats).mean, (stats).stddev, (stats).min, - (stats).max) + (stats).max + )::h3_raster_summary_stats $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_summary_stats_agg_transfn( @@ -181,7 +182,7 @@ AS $$ ), least((s1).min, (s2).min), greatest((s1).max, (s2).max) - ) + )::h3_raster_summary_stats FROM total AS t $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -238,7 +239,7 @@ AS $$ stddev_pop(val), min(val), max(val) - ) AS stats + )::h3_raster_summary_stats AS stats FROM pixels GROUP BY 1 $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -277,7 +278,7 @@ AS $$ 0, -- stddev val, -- min val -- max - ) + )::h3_rasters_summary_stats AS stats FROM vals; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -435,11 +436,12 @@ RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ SELECT h3_lat_lng_to_cell(geom, resolution) AS h3, - val::integer, + val::integer AS val, ROW( val::integer, count(*)::double precision, - count(*) * pixel_area) + count(*) * pixel_area + )::h3_raster_class_summary_item AS summary FROM ( -- x, y, val, geom SELECT (ST_PixelAsCentroids(rast, nband)).* @@ -489,11 +491,12 @@ AS $$ ) c) SELECT h3, - val::integer, + val::integer AS val, ROW( val::integer, cell_area / pixel_area, - cell_area) + cell_area + )::h3_raster_class_summary_item AS summary FROM vals v; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; From b9acb570bbc1b3a84f0e0154be0dc0a1bdffe20f Mon Sep 17 00:00:00 2001 From: sergei sh Date: Wed, 28 Dec 2022 15:50:43 +0300 Subject: [PATCH 04/22] PostGIS Raster integration: "CREATE OR REPLACE AGGREGATE" removed for compatibility with Postgres 11 --- h3_postgis/sql/install/40-rasters.sql | 4 ++-- h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 7070851c..a56c7f21 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -180,7 +180,7 @@ AS $$ $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --@ availability: unreleased -CREATE OR REPLACE AGGREGATE h3_raster_summary_stats_agg(h3_raster_summary_stats) ( +CREATE AGGREGATE h3_raster_summary_stats_agg(h3_raster_summary_stats) ( sfunc = __h3_raster_summary_stats_agg_transfn, stype = h3_raster_summary_stats, parallel = safe @@ -384,7 +384,7 @@ AS $$ $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --@ availability: unreleased -CREATE OR REPLACE AGGREGATE h3_raster_class_summary_item_agg(h3_raster_class_summary_item) ( +CREATE AGGREGATE h3_raster_class_summary_item_agg(h3_raster_class_summary_item) ( stype = h3_raster_class_summary_item, sfunc = __h3_raster_class_summary_item_agg_transfn, parallel = safe diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index 25c6324d..8424b02c 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -186,7 +186,7 @@ AS $$ FROM total AS t $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -CREATE OR REPLACE AGGREGATE h3_raster_summary_stats_agg(h3_raster_summary_stats) ( +CREATE AGGREGATE h3_raster_summary_stats_agg(h3_raster_summary_stats) ( sfunc = __h3_raster_summary_stats_agg_transfn, stype = h3_raster_summary_stats, parallel = safe @@ -356,7 +356,7 @@ AS $$ s1.area + s2.area $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -CREATE OR REPLACE AGGREGATE h3_raster_class_summary_item_agg(h3_raster_class_summary_item) ( +CREATE AGGREGATE h3_raster_class_summary_item_agg(h3_raster_class_summary_item) ( stype = h3_raster_class_summary_item, sfunc = __h3_raster_class_summary_item_agg_transfn, parallel = safe From 6e0f7dca36616280d9afe51015b762c126099ef3 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Thu, 29 Dec 2022 13:16:50 +0300 Subject: [PATCH 05/22] PostGIS Raster integration: documentation generation: "RETURNS TABLE", code blocks --- docs/api.md | 121 ++++++++++++++++++++++++++ h3_postgis/sql/install/40-rasters.sql | 6 +- scripts/sql2doc/generate.py | 68 ++++++++++++--- scripts/sql2doc/sql.lark | 15 +++- 4 files changed, 192 insertions(+), 18 deletions(-) diff --git a/docs/api.md b/docs/api.md index b38ac02c..9311cd1d 100644 --- a/docs/api.md +++ b/docs/api.md @@ -618,3 +618,124 @@ 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 +Combining summary stats from multiple rasters: +``` +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; +``` + + +*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`) + + +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 +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 + summary AS ( + SELECT h3, val, h3_raster_class_summary_item_agg(summary) AS item + FROM ( + -- h3, val, summary + SELECT (h3_raster_class_summary(rast, 8)).* AS summary + FROM rasters + ) t + 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; +``` + + +*Since vunreleased* + + +### h3_raster_class_summary_item_to_jsonb(item `h3_raster_class_summary_item`) ⇒ `jsonb` +*Since vunreleased* + + +Convert raster summary to binary JSON. + + +### 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. + + diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index a56c7f21..78a8d576 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -128,7 +128,7 @@ $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; -- NOTE: `count` can be < 1 when cell area is less than pixel area --@ availability: unreleased -CREATE TYPE h3_raster_summary_stats AS ( +CREATE TYPE h3_raster_summary_stats ( count double precision, sum double precision, mean double precision, @@ -282,7 +282,7 @@ $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_summary_subpixel( rast raster, resolution integer, - nband integer DEFAUlT 1) + nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE @@ -351,7 +351,7 @@ IS 'Returns `h3_raster_summary_stats` for each H3 cell in raster for a given ban -- NOTE: `count` can be < 1 when cell area is less than pixel area --@ availability: unreleased -CREATE TYPE h3_raster_class_summary_item AS ( +CREATE TYPE h3_raster_class_summary_item ( val integer, count double precision, area double precision diff --git a/scripts/sql2doc/generate.py b/scripts/sql2doc/generate.py index a4f4ce48..17c10fe5 100644 --- a/scripts/sql2doc/generate.py +++ b/scripts/sql2doc/generate.py @@ -26,9 +26,9 @@ import argparse import glob import re -import sys from pathlib import Path from lark import Lark, Transformer, v_args, visitors +import sys #TEST # see PGXN::API::Indexer::_clean_html_body def to_anchor(text): @@ -148,6 +148,39 @@ def __str__(self): return s +class Column: + def __init__(self, name, coltype): + self.name = name; + self.coltype = coltype; + + def __str__(self): + return "{} `{}`".format(self.name, self.coltype) + +class FunctionReturnsBase: + pass + + +class FunctionReturns(FunctionReturnsBase): + def __init__(self, rettype): + self.rettype = rettype + + def __str__(self): + return "`{}`".format(self.rettype) + + +class FunctionReturnsSet(FunctionReturns): + def __str__(self): + return "SETOF " + super().__str__() + + +class FunctionReturnsTable(FunctionReturnsBase): + def __init__(self, columns): + self.columns = columns + + def __str__(self): + return "TABLE (" + ", ".join([str(col) for col in self.columns]) + ")"; + + class CustomMd(StmtBase): def __init__(self, line): super().__init__() @@ -171,12 +204,11 @@ def __str__(self): class CreateFunctionStmt(StmtBase): - def __init__(self, name: str, arguments, returntype: str, returns_set: bool): + def __init__(self, name: str, arguments, returns: FunctionReturnsBase): super().__init__(2) self.name = name self.arguments = arguments or [] - self.returntype = returntype - self.returns_set = returns_set + self.returns = returns def get_ref_text(self): return "{}({})".format( @@ -184,11 +216,10 @@ def get_ref_text(self): ", ".join([arg.get_typestr() for arg in self.arguments])) def __str__(self): - return "{}({}) ⇒ {}`{}`".format( + return "{}({}) ⇒ {}".format( self.name, ", ".join([str(arg) for arg in self.arguments]), - "SETOF " if self.returns_set else "", - self.returntype) + self.returns) class CreateAggregateStmt(StmtBase): @@ -276,6 +307,7 @@ def custom_md_statement(self, children): # -- MARKDOWN -------------------------------------------------------------- @v_args(inline=True) def custom_markdown(self, line=""): + line = re.sub(r"^--\|[ \t]?", '', line) return CustomMd(line) # -- CREATE TYPE ----------------------------------------------------------- @@ -300,16 +332,29 @@ def create_oper_opt(self, option, value=None): return [str(option), value] # -- CREATE FUNCTION ------------------------------------------------------- + @v_args(inline=True) + def create_fun_ret_table_columns(self, columns): + return FunctionReturnsTable(columns) + def create_fun_rettype(self, children): - # (returntype, returns_set) - return (children[1], children[0] is not None) + rettype = children[1] + if children[0] is not None: + return FunctionReturnsSet(rettype) + return FunctionReturns(rettype) + + def column_list(self, children): + return children @v_args(inline=True) - def create_func_stmt(self, name: str, arguments, returntype, *opts): + def column(self, name: str, coltype: str): + return Column(name, coltype) + + @v_args(inline=True) + def create_func_stmt(self, name: str, arguments, returns, *opts): # skip internal functions if name.startswith("__"): raise visitors.Discard() - return CreateFunctionStmt(name, arguments, returntype[0], returntype[1]) + return CreateFunctionStmt(name, arguments, returns) # -- CREATE AGGREGATE ------------------------------------------------------ @v_args(inline=True) @@ -357,7 +402,6 @@ def CNAME(self, cname): def OPERATOR(self, name): return str(name) - def parse_args(): parser = argparse.ArgumentParser() parser.add_argument( diff --git a/scripts/sql2doc/sql.lark b/scripts/sql2doc/sql.lark index 3c34dd56..c242c08f 100644 --- a/scripts/sql2doc/sql.lark +++ b/scripts/sql2doc/sql.lark @@ -17,8 +17,8 @@ custom_decorated_statement: [custom_decorators] statement custom_decorators: ("--@" /([^\n])+/)+ -custom_markdown: "--|" /([^\n])+/ - | "--|\n" +custom_markdown: /--\|[ \t]?([^\n])+/ + | "--|\n" // ----------------------------------------------------------------------------- // --------------- Basic SQL below --------------------------------------------- @@ -83,12 +83,15 @@ create_oper_opts: create_oper_opt ("," create_oper_opt)* // } ... // [ WITH ( attribute [, ...] ) ] create_func_stmt: "CREATE" ("OR" "REPLACE")? "FUNCTION" fun_name "(" [argument_list] ")" [create_fun_rets] create_fun_opts* -?create_fun_rets: "RETURNS" create_fun_rettype +?create_fun_rets: ("RETURNS" "TABLE" "(" create_fun_ret_table_columns ")") + | ("RETURNS" create_fun_rettype) create_fun_opts: "LANGUAGE" (CNAME|"'" CNAME "'") | ("IMMUTABLE" | "STABLE" | "VOLATILE" | ("NOT"? "LEAKPROOF")) | (("CALLED" "ON" "NULL" "INPUT") | ("RETURNS" "NULL" "ON" "NULL" "INPUT") | "STRICT") | ("PARALLEL" ("UNSAFE" | "RESTRICTED" | "SAFE")) | "AS" string ("," string)? +create_fun_ret_table_columns: column_list +column_list: column ("," column)* argument_list: argument ("," argument)* !create_fun_rettype: ["SETOF"] DATATYPE @@ -141,9 +144,15 @@ comment_on_type: "CAST" "(" DATATYPE "AS" DATATYPE ")" -> comment_on_cast // ----------------------------------------------------------------------------- // SIMPLE RULES +column: CNAME DATATYPE argument: [ARGMODE] [CNAME] DATATYPE ("DEFAULT" expr)? ARGMODE.2: "IN" | "OUT" | "INOUT" DATATYPE_SCALAR: "h3index" + | "raster" + | "summarystats" + | "h3_raster_summary_stats" + | "h3_raster_class_summary_item" + | "jsonb" | "bigint" | "boolean" | "cstring" From ee3373c4c2ea1026ad2a2dfa00f32624407872f3 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Thu, 29 Dec 2022 13:48:36 +0300 Subject: [PATCH 06/22] PostGIS Raster integration: type creation fix --- h3_postgis/sql/install/40-rasters.sql | 4 ++-- scripts/sql2doc/sql.lark | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 78a8d576..3f6b561f 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -128,7 +128,7 @@ $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; -- NOTE: `count` can be < 1 when cell area is less than pixel area --@ availability: unreleased -CREATE TYPE h3_raster_summary_stats ( +CREATE TYPE h3_raster_summary_stats AS ( count double precision, sum double precision, mean double precision, @@ -351,7 +351,7 @@ IS 'Returns `h3_raster_summary_stats` for each H3 cell in raster for a given ban -- NOTE: `count` can be < 1 when cell area is less than pixel area --@ availability: unreleased -CREATE TYPE h3_raster_class_summary_item ( +CREATE TYPE h3_raster_class_summary_item AS ( val integer, count double precision, area double precision diff --git a/scripts/sql2doc/sql.lark b/scripts/sql2doc/sql.lark index c242c08f..81972016 100644 --- a/scripts/sql2doc/sql.lark +++ b/scripts/sql2doc/sql.lark @@ -26,7 +26,7 @@ custom_markdown: /--\|[ \t]?([^\n])+/ // ----------------------------------------------------------------------------- // CREATE TYPE name ... -create_type_stmt: "CREATE" "TYPE" CNAME ("(" /([^\)])+/ ")")? +create_type_stmt: "CREATE" "TYPE" CNAME "AS"? ("(" /([^\)])+/ ")")? // ----------------------------------------------------------------------------- // CREATE CAST ... From 3d004891c4e96e142215d18eb8c630e5dd3aff87 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Fri, 30 Dec 2022 14:17:12 +0300 Subject: [PATCH 07/22] PostGIS Raster integration: raster geometry transform fix --- h3_postgis/sql/install/40-rasters.sql | 101 +++++++++------- .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 114 +++++++++--------- 2 files changed, 112 insertions(+), 103 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 3f6b561f..75b5b6cb 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -17,6 +17,9 @@ --| # Raster processing functions -- Get nodata value for ST_Clip function +-- ST_Clip sets nodata pixel values to minimum value by default, but it won't +-- set band nodata value in this case, which we need later for filtering dumped +-- values. CREATE OR REPLACE FUNCTION __h3_raster_band_nodata( rast raster, nband integer) @@ -27,30 +30,30 @@ AS $$ ST_MinPossibleValue(ST_BandPixelType(rast, nband))); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +CREATE OR REPLACE FUNCTION __h3_raster_to_polygon(rast raster, nband integer) +RETURNS geometry +AS $$ + SELECT ST_MinConvexHull(rast, nband); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + CREATE OR REPLACE FUNCTION __h3_raster_pixel_area(rast raster) RETURNS double precision AS $$ - SELECT ST_Area(ST_PixelAsPolygon(rast, 1, 1)::geography); + SELECT ST_Area(ST_Transform(ST_PixelAsPolygon(rast, 1, 1), 4326)::geography); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -- Get area of H3 cell close to the center of the raster -CREATE OR REPLACE FUNCTION __h3_raster_cell_area(rast raster, resolution integer) +CREATE OR REPLACE FUNCTION __h3_raster_cell_area( + rast raster, + resolution integer, + nband integer) RETURNS double precision AS $$ -DECLARE - rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); -BEGIN SELECT ST_Area( h3_cell_to_boundary_geography( - ST_Transform(ST_Centroid(rast_geom), 4326), - resolution)); -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; - -CREATE OR REPLACE FUNCTION __h3_raster_to_polygon(rast raster, nband integer) -RETURNS geometry -AS $$ - SELECT ST_MinConvexHull(rast, nband); + h3_lat_lng_to_cell( + ST_Transform(ST_Centroid(__h3_raster_to_polygon(rast, nband)), 4326), + resolution))); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -- Get all H3 cells potentially intersecting the polygon, @@ -103,12 +106,15 @@ DECLARE rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); BEGIN RETURN QUERY - SELECT - c.h3, - ST_Transform(h3_cell_to_geometry(c.h3), ST_SRID(rast)) AS geom - FROM ( - SELECT h3_polygon_to_cells(rast_geom, resolution) AS h3 - ) c; + WITH geoms AS ( + SELECT + c.h3, + ST_Transform(h3_cell_to_geometry(c.h3), ST_SRID(rast)) AS geom + FROM ( + SELECT h3_polygon_to_cells(ST_Transform(rast_geom, 4326), resolution) AS h3 + ) c) + SELECT g.* + FROM geoms g; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; @@ -227,7 +233,7 @@ AS $$ -- x, y, val, geom SELECT (ST_PixelAsCentroids(rast, nband)).*) SELECT - h3_lat_lng_to_cell(geom, resolution) AS h3, + h3_lat_lng_to_cell(ST_Transform(geom, 4326), resolution) AS h3, ROW( count(val), sum(val), @@ -252,19 +258,23 @@ CREATE OR REPLACE FUNCTION __h3_raster_summary_subpixel( RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ WITH - vals AS ( + coords AS ( SELECT h3, - ST_Value( - rast, - nband, - ST_WorldToRasterCoordX(rast, geom), - ST_WorldToRasterCoordY(rast, geom) - ) AS val + ST_WorldToRasterCoordX(rast, geom) AS x, + ST_WorldToRasterCoordX(rast, geom) AS y FROM ( -- h3, geom SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* - ) t) + ) t), + vals AS ( + SELECT + h3, + ST_Value(rast, nband, x, y) AS val + FROM coords + WHERE + x BETWEEN 1 AND ST_Width(rast) + AND y BETWEEN 1 AND ST_Height(rast)) SELECT h3, ROW( @@ -274,7 +284,7 @@ AS $$ 0, -- stddev val, -- min val -- max - )::h3_rasters_summary_stats AS stats + )::h3_raster_summary_stats AS stats FROM vals; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -287,7 +297,7 @@ RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); BEGIN RETURN QUERY SELECT (__h3_raster_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; END; @@ -464,7 +474,7 @@ CREATE OR REPLACE FUNCTION __h3_raster_class_summary_centroids( RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ SELECT - h3_lat_lng_to_cell(geom, resolution) AS h3, + h3_lat_lng_to_cell(ST_Transform(geom, 4326), resolution) AS h3, val::integer AS val, ROW( val::integer, @@ -507,19 +517,23 @@ CREATE OR REPLACE FUNCTION __h3_raster_class_summary_subpixel( RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ WITH - vals AS ( + coords AS ( SELECT h3, - ST_Value( - rast, - nband, - ST_WorldToRasterCoordX(rast, geom), - ST_WorldToRasterCoordY(rast, geom) - ) AS val + ST_WorldToRasterCoordX(rast, geom) AS x, + ST_WorldToRasterCoordX(rast, geom) AS y FROM ( -- h3, geom SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* - ) c) + ) t), + vals AS ( + SELECT + h3, + ST_Value(rast, nband, x, y) AS val + FROM coords + WHERE + x BETWEEN 1 AND ST_Width(rast) + AND y BETWEEN 1 AND ST_Height(rast)) SELECT h3, val::integer AS val, @@ -528,7 +542,7 @@ AS $$ cell_area / pixel_area, cell_area )::h3_raster_class_summary_item AS summary - FROM vals v; + FROM vals; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -- Get summary items for each H3 index and value. @@ -541,7 +555,7 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ DECLARE - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); BEGIN RETURN QUERY SELECT (__h3_raster_class_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; @@ -551,7 +565,6 @@ COMMENT ON FUNCTION h3_raster_class_summary_subpixel(raster, integer, integer) IS '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.'; - -- Get summary items for each H3 index and value. -- Select appropriate method based on number of pixels per H3 cell. --@ availability: unreleased @@ -563,7 +576,7 @@ RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 350 THEN diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index 8424b02c..addaf187 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -37,7 +37,6 @@ Splits polygons when crossing 180th meridian.'; -- Raster processing --- Get nodata value for ST_Clip function CREATE OR REPLACE FUNCTION __h3_raster_band_nodata( rast raster, nband integer) @@ -48,34 +47,31 @@ AS $$ ST_MinPossibleValue(ST_BandPixelType(rast, nband))); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +CREATE OR REPLACE FUNCTION __h3_raster_to_polygon(rast raster, nband integer) +RETURNS geometry +AS $$ + SELECT ST_MinConvexHull(rast, nband); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + CREATE OR REPLACE FUNCTION __h3_raster_pixel_area(rast raster) RETURNS double precision AS $$ - SELECT ST_Area(ST_PixelAsPolygon(rast, 1, 1)::geography); + SELECT ST_Area(ST_Transform(ST_PixelAsPolygon(rast, 1, 1), 4326)::geography); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --- Get area of H3 cell close to the center of the raster -CREATE OR REPLACE FUNCTION __h3_raster_cell_area(rast raster, resolution integer) +CREATE OR REPLACE FUNCTION __h3_raster_cell_area( + rast raster, + resolution integer, + nband integer) RETURNS double precision AS $$ -DECLARE - rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); -BEGIN SELECT ST_Area( h3_cell_to_boundary_geography( - ST_Transform(ST_Centroid(rast_geom), 4326), - resolution)); -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; - -CREATE OR REPLACE FUNCTION __h3_raster_to_polygon(rast raster, nband integer) -RETURNS geometry -AS $$ - SELECT ST_MinConvexHull(rast, nband); + h3_lat_lng_to_cell( + ST_Transform(ST_Centroid(__h3_raster_to_polygon(rast, nband)), 4326), + resolution))); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --- Get all H3 cells potentially intersecting the polygon, --- result may contain cells outside of polygon CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cells( poly geometry, resolution integer) @@ -88,7 +84,6 @@ AS $$ resolution); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --- Get H3 cell geometries intersecting the raster CREATE OR REPLACE FUNCTION __h3_raster_to_cell_boundaries( rast raster, resolution integer, @@ -113,7 +108,6 @@ BEGIN END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; --- Get H3 cell centroids within the raster CREATE OR REPLACE FUNCTION __h3_raster_to_cell_centroids( rast raster, resolution integer, @@ -124,12 +118,15 @@ DECLARE rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); BEGIN RETURN QUERY - SELECT - c.h3, - ST_Transform(h3_cell_to_geometry(c.h3), ST_SRID(rast)) AS geom - FROM ( - SELECT h3_polygon_to_cells(rast_geom, resolution) AS h3 - ) c; + WITH geoms AS ( + SELECT + c.h3, + ST_Transform(h3_cell_to_geometry(c.h3), ST_SRID(rast)) AS geom + FROM ( + SELECT h3_polygon_to_cells(ST_Transform(rast_geom, 4326), resolution) AS h3 + ) c) + SELECT g.* + FROM geoms g; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; @@ -144,7 +141,6 @@ CREATE TYPE h3_raster_summary_stats AS ( max double precision ); --- ST_SummaryStats result type to h3_raster_summary_stats CREATE OR REPLACE FUNCTION __h3_raster_to_summary_stats(stats summarystats) RETURNS h3_raster_summary_stats AS $$ @@ -231,7 +227,7 @@ AS $$ -- x, y, val, geom SELECT (ST_PixelAsCentroids(rast, nband)).*) SELECT - h3_lat_lng_to_cell(geom, resolution) AS h3, + h3_lat_lng_to_cell(ST_Transform(geom, 4326), resolution) AS h3, ROW( count(val), sum(val), @@ -256,19 +252,23 @@ CREATE OR REPLACE FUNCTION __h3_raster_summary_subpixel( RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ WITH - vals AS ( + coords AS ( SELECT h3, - ST_Value( - rast, - nband, - ST_WorldToRasterCoordX(rast, geom), - ST_WorldToRasterCoordY(rast, geom) - ) AS val + ST_WorldToRasterCoordX(rast, geom) AS x, + ST_WorldToRasterCoordX(rast, geom) AS y FROM ( -- h3, geom SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* - ) t) + ) t), + vals AS ( + SELECT + h3, + ST_Value(rast, nband, x, y) AS val + FROM coords + WHERE + x BETWEEN 1 AND ST_Width(rast) + AND y BETWEEN 1 AND ST_Height(rast)) SELECT h3, ROW( @@ -278,19 +278,19 @@ AS $$ 0, -- stddev val, -- min val -- max - )::h3_rasters_summary_stats AS stats + )::h3_raster_summary_stats AS stats FROM vals; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_summary_subpixel( rast raster, resolution integer, - nband integer DEFAUlT 1) + nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); BEGIN RETURN QUERY SELECT (__h3_raster_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; END; @@ -362,7 +362,6 @@ CREATE AGGREGATE h3_raster_class_summary_item_agg(h3_raster_class_summary_item) parallel = safe ); --- Get summary items for a raster clipped by H3 cell geometry CREATE OR REPLACE FUNCTION __h3_raster_class_summary_part( rast raster, nband integer, @@ -435,7 +434,7 @@ CREATE OR REPLACE FUNCTION __h3_raster_class_summary_centroids( RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ SELECT - h3_lat_lng_to_cell(geom, resolution) AS h3, + h3_lat_lng_to_cell(ST_Transform(geom, 4326), resolution) AS h3, val::integer AS val, ROW( val::integer, @@ -449,8 +448,6 @@ AS $$ GROUP BY 1, 2; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; - --- For each pixel determine which H3 cell it belongs to then group by H3 index and value. CREATE OR REPLACE FUNCTION h3_raster_class_summary_centroids( rast raster, resolution integer, @@ -476,19 +473,23 @@ CREATE OR REPLACE FUNCTION __h3_raster_class_summary_subpixel( RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ WITH - vals AS ( + coords AS ( SELECT h3, - ST_Value( - rast, - nband, - ST_WorldToRasterCoordX(rast, geom), - ST_WorldToRasterCoordY(rast, geom) - ) AS val + ST_WorldToRasterCoordX(rast, geom) AS x, + ST_WorldToRasterCoordX(rast, geom) AS y FROM ( -- h3, geom SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* - ) c) + ) t), + vals AS ( + SELECT + h3, + ST_Value(rast, nband, x, y) AS val + FROM coords + WHERE + x BETWEEN 1 AND ST_Width(rast) + AND y BETWEEN 1 AND ST_Height(rast)) SELECT h3, val::integer AS val, @@ -497,11 +498,9 @@ AS $$ cell_area / pixel_area, cell_area )::h3_raster_class_summary_item AS summary - FROM vals v; + FROM vals; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --- Get summary items for each H3 index and value. --- For each H3 cell centroid determine which pixel it belongs to. CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( rast raster, resolution integer, @@ -509,7 +508,7 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ DECLARE - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); BEGIN RETURN QUERY SELECT (__h3_raster_class_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; @@ -519,9 +518,6 @@ COMMENT ON FUNCTION h3_raster_class_summary_subpixel(raster, integer, integer) IS '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.'; - --- Get summary items for each H3 index and value. --- Select appropriate method based on number of pixels per H3 cell. CREATE OR REPLACE FUNCTION h3_raster_class_summary( rast raster, resolution integer, @@ -530,7 +526,7 @@ RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); + cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 350 THEN From 423470f239c2351096b631ed61e0f27f9ba52e1c Mon Sep 17 00:00:00 2001 From: sergei sh Date: Mon, 2 Jan 2023 23:19:09 +0300 Subject: [PATCH 08/22] PostGIS Raster integration: rewritten to re-use raster polygon --- h3_postgis/sql/install/40-rasters.sql | 430 +++++++++++------- .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 408 ++++++++++------- 2 files changed, 501 insertions(+), 337 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 75b5b6cb..be2e541d 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -30,93 +30,164 @@ AS $$ ST_MinPossibleValue(ST_BandPixelType(rast, nband))); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -CREATE OR REPLACE FUNCTION __h3_raster_to_polygon(rast raster, nband integer) +CREATE OR REPLACE FUNCTION __h3_raster_to_polygon( + rast raster, + nband integer) RETURNS geometry AS $$ SELECT ST_MinConvexHull(rast, nband); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +-- Area of a single pixel in meters CREATE OR REPLACE FUNCTION __h3_raster_pixel_area(rast raster) RETURNS double precision AS $$ SELECT ST_Area(ST_Transform(ST_PixelAsPolygon(rast, 1, 1), 4326)::geography); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --- Get area of H3 cell close to the center of the raster -CREATE OR REPLACE FUNCTION __h3_raster_cell_area( - rast raster, - resolution integer, - nband integer) +-- Area of a cell close to the center of raster polygon, in meters +CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( + poly geometry, + resolution integer) +RETURNS h3index +AS $$ + SELECT h3_lat_lng_to_cell(ST_Transform(ST_Centroid(poly), 4326), resolution); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell_area( + poly geometry, + resolution integer) RETURNS double precision AS $$ SELECT ST_Area( h3_cell_to_boundary_geography( - h3_lat_lng_to_cell( - ST_Transform(ST_Centroid(__h3_raster_to_polygon(rast, nband)), 4326), - resolution))); + __h3_raster_polygon_centroid_cell(poly, resolution))); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --- Get all H3 cells potentially intersecting the polygon, --- result may contain cells outside of polygon +-- Get list of cells inside of the raster polygon, +-- buffered by `buffer` value (in meters). +-- If SRID != 4326 then additionally buffer by 1 pixel to account for transformation. CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cells( + rast raster, poly geometry, - resolution integer) + resolution integer, + buffer double precision) RETURNS SETOF h3index AS $$ - SELECT h3_polygon_to_cells( - ST_Buffer( - ST_Transform(poly, 4326)::geography, - h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3), - resolution); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +DECLARE + buffered geometry := poly; +BEGIN + IF ST_SRID(rast) != 4326 THEN + buffered := ST_Transform( + ST_Buffer(poly, greatest(ST_PixelWidth(rast), ST_PixelHeight(rast))), + 4326); + END IF; + IF buffer > 0.0 THEN + RETURN QUERY + SELECT h3_polygon_to_cells( + ST_Buffer(buffered::geography, buffer * 1.3), + resolution); + ELSE + RETURN QUERY + SELECT h3_polygon_to_cells(buffered, resolution); + END IF; +END +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; --- Get H3 cell geometries intersecting the raster -CREATE OR REPLACE FUNCTION __h3_raster_to_cell_boundaries( +-- Get geometries of H3 cells interesecting raster polygon. +CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_boundaries_intersects( rast raster, - resolution integer, - nband integer) + poly geometry, + resolution integer) RETURNS TABLE (h3 h3index, geom geometry) AS $$ -DECLARE - rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); -BEGIN - RETURN QUERY WITH geoms AS ( SELECT c.h3, ST_Transform(h3_cell_to_boundary_geometry(c.h3), ST_SRID(rast)) AS geom FROM ( - SELECT __h3_raster_polygon_to_cells(rast_geom, resolution) AS h3 + SELECT __h3_raster_polygon_to_cells( + rast, + poly, + resolution, + h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3) + AS h3 ) c) SELECT g.* FROM geoms g - WHERE ST_Intersects(g.geom, rast_geom); -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + WHERE ST_Intersects(g.geom, poly); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +-- Get raster coordinates of H3 cells with centroids inside the raster polygon +CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_coords_centroid( + rast raster, + poly geometry, + resolution integer) +RETURNS TABLE (h3 h3index, x integer, y integer) +AS $$ + WITH + geoms AS ( + SELECT + h3, + ST_Transform( + h3_cell_to_geometry(h3), + ST_SRID(poly) + ) AS geom + FROM ( + SELECT __h3_raster_polygon_to_cells(rast, poly, resolution, 0.0) AS h3 + ) t), + coords AS ( + SELECT + h3, + ST_WorldToRasterCoordX(rast, geom) AS x, + ST_WorldToRasterCoordY(rast, geom) AS y + FROM geoms) + SELECT h3, x, y + FROM coords + WHERE + x BETWEEN 1 AND ST_Width(rast) + AND y BETWEEN 1 AND ST_Height(rast); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --- Get H3 cell centroids within the raster -CREATE OR REPLACE FUNCTION __h3_raster_to_cell_centroids( +CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_parts( rast raster, + poly geometry, resolution integer, nband integer) -RETURNS TABLE (h3 h3index, geom geometry) +RETURNS TABLE (h3 h3index, part raster) AS $$ -DECLARE - rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); -BEGIN - RETURN QUERY - WITH geoms AS ( - SELECT - c.h3, - ST_Transform(h3_cell_to_geometry(c.h3), ST_SRID(rast)) AS geom - FROM ( - SELECT h3_polygon_to_cells(ST_Transform(rast_geom, 4326), resolution) AS h3 - ) c) - SELECT g.* - FROM geoms g; -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + WITH + parts AS ( + SELECT + h3, + ST_Clip(rast, nband, geom, __h3_raster_band_nodata(rast, nband)) AS part + FROM ( + -- h3, geom + SELECT (__h3_raster_polygon_to_cell_boundaries_intersects(rast, poly, resolution)).* + ) t) + SELECT h3, part + FROM parts + WHERE NOT ST_BandIsNoData(part, nband); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +-- Get values corresponding to all H3 cells with centroids inside the +-- raster polygon. Assumes cell area is less than pixel area. +CREATE OR REPLACE FUNCTION __h3_raster_polygon_subpixel_cell_values( + rast raster, + poly geometry, + resolution integer, + nband integer) +RETURNS TABLE (h3 h3index, val double precision) +AS $$ + SELECT + h3, + ST_Value(rast, nband, x, y) AS val + FROM ( + -- h3, x, y + SELECT (__h3_raster_polygon_to_cell_coords_centroid(rast, poly, resolution)).* + ) t; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --| ## Continuous raster data --| @@ -192,6 +263,22 @@ CREATE AGGREGATE h3_raster_summary_stats_agg(h3_raster_summary_stats) ( parallel = safe ); +CREATE OR REPLACE FUNCTION __h3_raster_polygon_summary_clip( + rast raster, + poly geometry, + resolution integer, + nband integer) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ + SELECT + h3, + __h3_raster_to_summary_stats(ST_SummaryStats(part, nband, TRUE)) AS stats + FROM ( + -- h3, part + SELECT (__h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband)).* + ) t; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + --@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_summary_clip( rast raster, @@ -199,25 +286,13 @@ CREATE OR REPLACE FUNCTION h3_raster_summary_clip( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ -DECLARE - nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); -BEGIN - RETURN QUERY - WITH parts AS ( - SELECT - g.h3, - ST_Clip(rast, nband, g.geom, nodata) AS part - FROM ( - -- h3, geom - SELECT (__h3_raster_to_cell_boundaries(rast, resolution, nband)).* - ) g) - SELECT - p.h3, - __h3_raster_to_summary_stats(ST_SummaryStats(part, nband, TRUE)) AS stats - FROM parts AS p - WHERE NOT ST_BandIsNoData(part, 1); -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + SELECT (__h3_raster_polygon_summary_clip( + rast, + __h3_raster_to_polygon(rast, nband), + resolution, + nband + )).*; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_summary_clip(raster, integer, integer) IS '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.'; @@ -229,9 +304,6 @@ CREATE OR REPLACE FUNCTION h3_raster_summary_centroids( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ - WITH pixels AS ( - -- x, y, val, geom - SELECT (ST_PixelAsCentroids(rast, nband)).*) SELECT h3_lat_lng_to_cell(ST_Transform(geom, 4326), resolution) AS h3, ROW( @@ -242,50 +314,38 @@ AS $$ min(val), max(val) )::h3_raster_summary_stats AS stats - FROM pixels - GROUP BY 1 + FROM ( + -- x, y, val, geom + SELECT (ST_PixelAsCentroids(rast, nband)).* + ) t + GROUP BY 1; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_summary_centroids(raster, integer, integer) IS '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.'; -CREATE OR REPLACE FUNCTION __h3_raster_summary_subpixel( +CREATE OR REPLACE FUNCTION __h3_raster_polygon_summary_subpixel( rast raster, + poly geometry, resolution integer, nband integer, - pixel_area double precision, - cell_area double precision) + pixels_per_cell double precision) RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ - WITH - coords AS ( - SELECT - h3, - ST_WorldToRasterCoordX(rast, geom) AS x, - ST_WorldToRasterCoordX(rast, geom) AS y - FROM ( - -- h3, geom - SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* - ) t), - vals AS ( - SELECT - h3, - ST_Value(rast, nband, x, y) AS val - FROM coords - WHERE - x BETWEEN 1 AND ST_Width(rast) - AND y BETWEEN 1 AND ST_Height(rast)) SELECT h3, ROW( - cell_area / pixel_area, -- count + pixels_per_cell, -- count val, -- sum val, -- mean - 0, -- stddev + 0.0, -- stddev val, -- min val -- max )::h3_raster_summary_stats AS stats - FROM vals; + FROM ( + -- h3, val + SELECT (__h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband)).* + ) t; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --@ availability: unreleased @@ -296,16 +356,24 @@ CREATE OR REPLACE FUNCTION h3_raster_summary_subpixel( RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE + poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); + cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); BEGIN - RETURN QUERY SELECT (__h3_raster_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + RETURN QUERY SELECT (__h3_raster_polygon_summary_subpixel( + rast, + poly, + resolution, + nband, + cell_area / pixel_area) + ).*; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_summary_subpixel(raster, integer, integer) IS '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.'; +--@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_summary( rast raster, resolution integer, @@ -313,16 +381,32 @@ CREATE OR REPLACE FUNCTION h3_raster_summary( RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); - pixels_per_cell CONSTANT double precision := cell_area / pixel_area; + poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); + pixels_per_cell CONSTANT double precision := + __h3_raster_polygon_centroid_cell_area(poly, resolution) + / __h3_raster_pixel_area(rast); BEGIN IF pixels_per_cell > 350 THEN - RETURN QUERY SELECT (h3_raster_summary_clip(rast, resolution, nband)).*; + RETURN QUERY SELECT (__h3_raster_polygon_summary_clip( + rast, + poly, + resolution, + nband + )).*; ELSIF pixels_per_cell > 1 THEN - RETURN QUERY SELECT (h3_raster_summary_centroids(rast, resolution, nband)).*; + RETURN QUERY SELECT (h3_raster_summary_centroids( + rast, + resolution, + nband + )).*; ELSE - RETURN QUERY SELECT (__h3_raster_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + RETURN QUERY SELECT (__h3_raster_polygon_summary_subpixel( + rast, + poly, + resolution, + nband, + pixels_per_cell + )).*; END IF; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; @@ -400,7 +484,6 @@ CREATE AGGREGATE h3_raster_class_summary_item_agg(h3_raster_class_summary_item) parallel = safe ); --- Get summary items for a raster clipped by H3 cell geometry CREATE OR REPLACE FUNCTION __h3_raster_class_summary_part( rast raster, nband integer, @@ -418,35 +501,26 @@ AS $$ GROUP BY 1 $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -CREATE OR REPLACE FUNCTION __h3_raster_class_summary_clip( +CREATE OR REPLACE FUNCTION __h3_raster_class_polygon_summary_clip( rast raster, + poly geometry, resolution integer, nband integer, pixel_area double precision) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ -DECLARE - nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); -BEGIN - RETURN QUERY WITH - parts AS ( - SELECT - g.h3, - ST_Clip(rast, nband, g.geom, nodata, TRUE) AS part - FROM ( - -- h3, geom - SELECT (__h3_raster_to_cell_boundaries(rast, resolution, nband)).* - ) g), summary AS ( SELECT - p.h3, + h3, __h3_raster_class_summary_part(part, nband, pixel_area) AS summary - FROM parts p) - SELECT s.h3, (s.summary).val, s.summary - FROM summary s; -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + FROM ( + -- h3, part + SELECT (__h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband)).* + ) t) + SELECT h3, (summary).val, summary + FROM summary; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( @@ -455,17 +529,18 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ -DECLARE - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); -BEGIN - RETURN QUERY SELECT (__h3_raster_class_summary_clip(rast, resolution, nband, pixel_area)).*; -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + SELECT (__h3_raster_class_polygon_summary_clip( + rast, + __h3_raster_to_polygon(rast, nband), + resolution, + nband, + __h3_raster_pixel_area(rast) + )).* +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_clip(raster, integer, integer) IS '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.'; - CREATE OR REPLACE FUNCTION __h3_raster_class_summary_centroids( rast raster, resolution integer, @@ -488,8 +563,6 @@ AS $$ GROUP BY 1, 2; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; - --- For each pixel determine which H3 cell it belongs to then group by H3 index and value. --@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_class_summary_centroids( rast raster, @@ -497,43 +570,26 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_centroids( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ -DECLARE - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); -BEGIN - RETURN QUERY SELECT (__h3_raster_class_summary_centroids(rast, resolution, nband, pixel_area)).*; -END -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + SELECT (__h3_raster_class_summary_centroids( + rast, + resolution, + nband, + __h3_raster_pixel_area(rast) + )).* +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_centroids(raster, integer, integer) IS '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.'; ---@ availability: unreleased -CREATE OR REPLACE FUNCTION __h3_raster_class_summary_subpixel( +CREATE OR REPLACE FUNCTION __h3_raster_class_polygon_summary_subpixel( rast raster, + poly geometry, resolution integer, nband integer, - pixel_area double precision, - cell_area double precision) + cell_area double precision, + pixel_area double precision) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ - WITH - coords AS ( - SELECT - h3, - ST_WorldToRasterCoordX(rast, geom) AS x, - ST_WorldToRasterCoordX(rast, geom) AS y - FROM ( - -- h3, geom - SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* - ) t), - vals AS ( - SELECT - h3, - ST_Value(rast, nband, x, y) AS val - FROM coords - WHERE - x BETWEEN 1 AND ST_Width(rast) - AND y BETWEEN 1 AND ST_Height(rast)) SELECT h3, val::integer AS val, @@ -542,11 +598,12 @@ AS $$ cell_area / pixel_area, cell_area )::h3_raster_class_summary_item AS summary - FROM vals; + FROM ( + -- h3, val + SELECT (__h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband)).* + ) t; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --- Get summary items for each H3 index and value. --- For each H3 cell centroid determine which pixel it belongs to. --@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( rast raster, @@ -555,18 +612,22 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ DECLARE - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); BEGIN - RETURN QUERY SELECT (__h3_raster_class_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + RETURN QUERY SELECT (__h3_raster_class_polygon_summary_subpixel( + rast, + poly, + resolution, + nband, + __h3_raster_polygon_centroid_cell_area(poly, resolution), + __h3_raster_pixel_area(rast) + )).*; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_subpixel(raster, integer, integer) IS '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.'; --- Get summary items for each H3 index and value. --- Select appropriate method based on number of pixels per H3 cell. --@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_class_summary( rast raster, @@ -575,16 +636,35 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary( RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ DECLARE - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); - pixels_per_cell CONSTANT double precision := cell_area / pixel_area; + poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); + cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 350 THEN - RETURN QUERY SELECT (__h3_raster_class_summary_clip(rast, resolution, nband, pixel_area)).*; + RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( + rast, + poly, + resolution, + nband, + pixel_area + )).*; ELSIF pixels_per_cell > 1 THEN - RETURN QUERY SELECT (__h3_raster_class_summary_centroids(rast, resolution, nband, pixel_area)).*; + RETURN QUERY SELECT (__h3_raster_class_summary_centroids( + rast, + resolution, + nband, + pixel_area + )).*; ELSE - RETURN QUERY SELECT (__h3_raster_class_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + RETURN QUERY SELECT (__h3_raster_class_polygon_summary_subpixel( + rast, + poly, + resolution, + nband, + cell_area, + pixel_area + )).*; END IF; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index addaf187..7069e228 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -47,7 +47,9 @@ AS $$ ST_MinPossibleValue(ST_BandPixelType(rast, nband))); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -CREATE OR REPLACE FUNCTION __h3_raster_to_polygon(rast raster, nband integer) +CREATE OR REPLACE FUNCTION __h3_raster_to_polygon( + rast raster, + nband integer) RETURNS geometry AS $$ SELECT ST_MinConvexHull(rast, nband); @@ -59,76 +61,141 @@ AS $$ SELECT ST_Area(ST_Transform(ST_PixelAsPolygon(rast, 1, 1), 4326)::geography); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -CREATE OR REPLACE FUNCTION __h3_raster_cell_area( - rast raster, - resolution integer, - nband integer) +CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( + poly geometry, + resolution integer) +RETURNS h3index +AS $$ + SELECT h3_lat_lng_to_cell(ST_Transform(ST_Centroid(poly), 4326), resolution); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell_area( + poly geometry, + resolution integer) RETURNS double precision AS $$ SELECT ST_Area( h3_cell_to_boundary_geography( - h3_lat_lng_to_cell( - ST_Transform(ST_Centroid(__h3_raster_to_polygon(rast, nband)), 4326), - resolution))); + __h3_raster_polygon_centroid_cell(poly, resolution))); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cells( + rast raster, poly geometry, - resolution integer) + resolution integer, + buffer double precision) RETURNS SETOF h3index AS $$ - SELECT h3_polygon_to_cells( - ST_Buffer( - ST_Transform(poly, 4326)::geography, - h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3), - resolution); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +DECLARE + buffered geometry := poly; +BEGIN + IF ST_SRID(rast) != 4326 THEN + buffered := ST_Transform( + ST_Buffer(poly, greatest(ST_PixelWidth(rast), ST_PixelHeight(rast))), + 4326); + END IF; + IF buffer > 0.0 THEN + RETURN QUERY + SELECT h3_polygon_to_cells( + ST_Buffer(buffered::geography, buffer * 1.3), + resolution); + ELSE + RETURN QUERY + SELECT h3_polygon_to_cells(buffered, resolution); + END IF; +END +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; -CREATE OR REPLACE FUNCTION __h3_raster_to_cell_boundaries( +CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_boundaries_intersects( rast raster, - resolution integer, - nband integer) + poly geometry, + resolution integer) RETURNS TABLE (h3 h3index, geom geometry) AS $$ -DECLARE - rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); -BEGIN - RETURN QUERY WITH geoms AS ( SELECT c.h3, ST_Transform(h3_cell_to_boundary_geometry(c.h3), ST_SRID(rast)) AS geom FROM ( - SELECT __h3_raster_polygon_to_cells(rast_geom, resolution) AS h3 + SELECT __h3_raster_polygon_to_cells( + rast, + poly, + resolution, + h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3) + AS h3 ) c) SELECT g.* FROM geoms g - WHERE ST_Intersects(g.geom, rast_geom); -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + WHERE ST_Intersects(g.geom, poly); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -CREATE OR REPLACE FUNCTION __h3_raster_to_cell_centroids( +CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_coords_centroid( rast raster, + poly geometry, + resolution integer) +RETURNS TABLE (h3 h3index, x integer, y integer) +AS $$ + WITH + geoms AS ( + SELECT + h3, + ST_Transform( + h3_cell_to_geometry(h3), + ST_SRID(poly) + ) AS geom + FROM ( + SELECT __h3_raster_polygon_to_cells(rast, poly, resolution, 0.0) AS h3 + ) t), + coords AS ( + SELECT + h3, + ST_WorldToRasterCoordX(rast, geom) AS x, + ST_WorldToRasterCoordY(rast, geom) AS y + FROM geoms) + SELECT h3, x, y + FROM coords + WHERE + x BETWEEN 1 AND ST_Width(rast) + AND y BETWEEN 1 AND ST_Height(rast); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_parts( + rast raster, + poly geometry, resolution integer, nband integer) -RETURNS TABLE (h3 h3index, geom geometry) +RETURNS TABLE (h3 h3index, part raster) AS $$ -DECLARE - rast_geom CONSTANT geometry := __h3_raster_to_polygon(rast, nband); -BEGIN - RETURN QUERY - WITH geoms AS ( - SELECT - c.h3, - ST_Transform(h3_cell_to_geometry(c.h3), ST_SRID(rast)) AS geom - FROM ( - SELECT h3_polygon_to_cells(ST_Transform(rast_geom, 4326), resolution) AS h3 - ) c) - SELECT g.* - FROM geoms g; -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + WITH + parts AS ( + SELECT + h3, + ST_Clip(rast, nband, geom, __h3_raster_band_nodata(rast, nband)) AS part + FROM ( + -- h3, geom + SELECT (__h3_raster_polygon_to_cell_boundaries_intersects(rast, poly, resolution)).* + ) t) + SELECT h3, part + FROM parts + WHERE NOT ST_BandIsNoData(part, nband); +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION __h3_raster_polygon_subpixel_cell_values( + rast raster, + poly geometry, + resolution integer, + nband integer) +RETURNS TABLE (h3 h3index, val double precision) +AS $$ + SELECT + h3, + ST_Value(rast, nband, x, y) AS val + FROM ( + -- h3, x, y + SELECT (__h3_raster_polygon_to_cell_coords_centroid(rast, poly, resolution)).* + ) t; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -- Raster processing: continuous data @@ -188,31 +255,35 @@ CREATE AGGREGATE h3_raster_summary_stats_agg(h3_raster_summary_stats) ( parallel = safe ); -CREATE OR REPLACE FUNCTION h3_raster_summary_clip( +CREATE OR REPLACE FUNCTION __h3_raster_polygon_summary_clip( rast raster, + poly geometry, resolution integer, - nband integer DEFAULT 1) + nband integer) RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ -DECLARE - nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); -BEGIN - RETURN QUERY - WITH parts AS ( - SELECT - g.h3, - ST_Clip(rast, nband, g.geom, nodata) AS part - FROM ( - -- h3, geom - SELECT (__h3_raster_to_cell_boundaries(rast, resolution, nband)).* - ) g) SELECT - p.h3, + h3, __h3_raster_to_summary_stats(ST_SummaryStats(part, nband, TRUE)) AS stats - FROM parts AS p - WHERE NOT ST_BandIsNoData(part, 1); -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + FROM ( + -- h3, part + SELECT (__h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband)).* + ) t; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OR REPLACE FUNCTION h3_raster_summary_clip( + rast raster, + resolution integer, + nband integer DEFAULT 1) +RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) +AS $$ + SELECT (__h3_raster_polygon_summary_clip( + rast, + __h3_raster_to_polygon(rast, nband), + resolution, + nband + )).*; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_summary_clip(raster, integer, integer) IS '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.'; @@ -223,9 +294,6 @@ CREATE OR REPLACE FUNCTION h3_raster_summary_centroids( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ - WITH pixels AS ( - -- x, y, val, geom - SELECT (ST_PixelAsCentroids(rast, nband)).*) SELECT h3_lat_lng_to_cell(ST_Transform(geom, 4326), resolution) AS h3, ROW( @@ -236,50 +304,39 @@ AS $$ min(val), max(val) )::h3_raster_summary_stats AS stats - FROM pixels - GROUP BY 1 + FROM ( + -- x, y, val, geom + SELECT (ST_PixelAsCentroids(rast, nband)).* + ) t + GROUP BY 1; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_summary_centroids(raster, integer, integer) IS '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.'; -CREATE OR REPLACE FUNCTION __h3_raster_summary_subpixel( + +CREATE OR REPLACE FUNCTION __h3_raster_polygon_summary_subpixel( rast raster, + poly geometry, resolution integer, nband integer, - pixel_area double precision, - cell_area double precision) + pixels_per_cell double precision) RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ - WITH - coords AS ( - SELECT - h3, - ST_WorldToRasterCoordX(rast, geom) AS x, - ST_WorldToRasterCoordX(rast, geom) AS y - FROM ( - -- h3, geom - SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* - ) t), - vals AS ( - SELECT - h3, - ST_Value(rast, nband, x, y) AS val - FROM coords - WHERE - x BETWEEN 1 AND ST_Width(rast) - AND y BETWEEN 1 AND ST_Height(rast)) SELECT h3, ROW( - cell_area / pixel_area, -- count + pixels_per_cell, -- count val, -- sum val, -- mean - 0, -- stddev + 0.0, -- stddev val, -- min val -- max )::h3_raster_summary_stats AS stats - FROM vals; + FROM ( + -- h3, val + SELECT (__h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband)).* + ) t; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_summary_subpixel( @@ -289,10 +346,17 @@ CREATE OR REPLACE FUNCTION h3_raster_summary_subpixel( RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE + poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); + cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); BEGIN - RETURN QUERY SELECT (__h3_raster_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + RETURN QUERY SELECT (__h3_raster_polygon_summary_subpixel( + rast, + poly, + resolution, + nband, + cell_area / pixel_area) + ).*; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION @@ -306,16 +370,32 @@ CREATE OR REPLACE FUNCTION h3_raster_summary( RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution); - pixels_per_cell CONSTANT double precision := cell_area / pixel_area; + poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); + pixels_per_cell CONSTANT double precision := + __h3_raster_polygon_centroid_cell_area(poly, resolution) + / __h3_raster_pixel_area(rast); BEGIN IF pixels_per_cell > 350 THEN - RETURN QUERY SELECT (h3_raster_summary_clip(rast, resolution, nband)).*; + RETURN QUERY SELECT (__h3_raster_polygon_summary_clip( + rast, + poly, + resolution, + nband + )).*; ELSIF pixels_per_cell > 1 THEN - RETURN QUERY SELECT (h3_raster_summary_centroids(rast, resolution, nband)).*; + RETURN QUERY SELECT (h3_raster_summary_centroids( + rast, + resolution, + nband + )).*; ELSE - RETURN QUERY SELECT (__h3_raster_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + RETURN QUERY SELECT (__h3_raster_polygon_summary_subpixel( + rast, + poly, + resolution, + nband, + pixels_per_cell + )).*; END IF; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; @@ -379,35 +459,26 @@ AS $$ GROUP BY 1 $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -CREATE OR REPLACE FUNCTION __h3_raster_class_summary_clip( +CREATE OR REPLACE FUNCTION __h3_raster_class_polygon_summary_clip( rast raster, + poly geometry, resolution integer, nband integer, pixel_area double precision) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ -DECLARE - nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); -BEGIN - RETURN QUERY WITH - parts AS ( - SELECT - g.h3, - ST_Clip(rast, nband, g.geom, nodata, TRUE) AS part - FROM ( - -- h3, geom - SELECT (__h3_raster_to_cell_boundaries(rast, resolution, nband)).* - ) g), summary AS ( SELECT - p.h3, + h3, __h3_raster_class_summary_part(part, nband, pixel_area) AS summary - FROM parts p) - SELECT s.h3, (s.summary).val, s.summary - FROM summary s; -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + FROM ( + -- h3, part + SELECT (__h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband)).* + ) t) + SELECT h3, (summary).val, summary + FROM summary; +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( rast raster, @@ -415,17 +486,18 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ -DECLARE - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); -BEGIN - RETURN QUERY SELECT (__h3_raster_class_summary_clip(rast, resolution, nband, pixel_area)).*; -END; -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + SELECT (__h3_raster_class_polygon_summary_clip( + rast, + __h3_raster_to_polygon(rast, nband), + resolution, + nband, + __h3_raster_pixel_area(rast) + )).* +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_clip(raster, integer, integer) IS '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.'; - CREATE OR REPLACE FUNCTION __h3_raster_class_summary_centroids( rast raster, resolution integer, @@ -454,42 +526,26 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_centroids( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ -DECLARE - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); -BEGIN - RETURN QUERY SELECT (__h3_raster_class_summary_centroids(rast, resolution, nband, pixel_area)).*; -END -$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; + SELECT (__h3_raster_class_summary_centroids( + rast, + resolution, + nband, + __h3_raster_pixel_area(rast) + )).* +$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_centroids(raster, integer, integer) IS '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.'; -CREATE OR REPLACE FUNCTION __h3_raster_class_summary_subpixel( +CREATE OR REPLACE FUNCTION __h3_raster_class_polygon_summary_subpixel( rast raster, + poly geometry, resolution integer, nband integer, - pixel_area double precision, - cell_area double precision) + cell_area double precision, + pixel_area double precision) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ - WITH - coords AS ( - SELECT - h3, - ST_WorldToRasterCoordX(rast, geom) AS x, - ST_WorldToRasterCoordX(rast, geom) AS y - FROM ( - -- h3, geom - SELECT (__h3_raster_to_cell_centroids(rast, resolution, nband)).* - ) t), - vals AS ( - SELECT - h3, - ST_Value(rast, nband, x, y) AS val - FROM coords - WHERE - x BETWEEN 1 AND ST_Width(rast) - AND y BETWEEN 1 AND ST_Height(rast)) SELECT h3, val::integer AS val, @@ -498,7 +554,10 @@ AS $$ cell_area / pixel_area, cell_area )::h3_raster_class_summary_item AS summary - FROM vals; + FROM ( + -- h3, val + SELECT (__h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband)).* + ) t; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( @@ -508,10 +567,16 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ DECLARE - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); BEGIN - RETURN QUERY SELECT (__h3_raster_class_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + RETURN QUERY SELECT (__h3_raster_class_polygon_summary_subpixel( + rast, + poly, + resolution, + nband, + __h3_raster_polygon_centroid_cell_area(poly, resolution), + __h3_raster_pixel_area(rast) + )).*; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION @@ -525,16 +590,35 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary( RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ DECLARE - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); - cell_area CONSTANT double precision := __h3_raster_cell_area(rast, resolution, nband); - pixels_per_cell CONSTANT double precision := cell_area / pixel_area; + poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); + cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 350 THEN - RETURN QUERY SELECT (__h3_raster_class_summary_clip(rast, resolution, nband, pixel_area)).*; + RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( + rast, + poly, + resolution, + nband, + pixel_area + )).*; ELSIF pixels_per_cell > 1 THEN - RETURN QUERY SELECT (__h3_raster_class_summary_centroids(rast, resolution, nband, pixel_area)).*; + RETURN QUERY SELECT (__h3_raster_class_summary_centroids( + rast, + resolution, + nband, + pixel_area + )).*; ELSE - RETURN QUERY SELECT (__h3_raster_class_summary_subpixel(rast, resolution, nband, pixel_area, cell_area)).*; + RETURN QUERY SELECT (__h3_raster_class_polygon_summary_subpixel( + rast, + poly, + resolution, + nband, + cell_area, + pixel_area + )).*; END IF; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; From cc81a6a196dc4162482ea9f7358d8a2509081016 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Mon, 2 Jan 2023 23:27:09 +0300 Subject: [PATCH 09/22] PostGIS Raster integration: minor documentation update --- docs/api.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/api.md b/docs/api.md index 9311cd1d..ddafd071 100644 --- a/docs/api.md +++ b/docs/api.md @@ -663,6 +663,7 @@ Returns `h3_raster_summary_stats` for each H3 cell in raster for a given band. A ### 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. From a95b8ce72f3dd461d416795873ba3e824b765a19 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Thu, 5 Jan 2023 01:21:59 +0300 Subject: [PATCH 10/22] PostGIS Raster integration: redundant subqueries removed --- h3_postgis/sql/install/40-rasters.sql | 111 +++++++----------- .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 104 ++++++---------- 2 files changed, 81 insertions(+), 134 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index be2e541d..3bdbd387 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -85,7 +85,7 @@ BEGIN IF buffer > 0.0 THEN RETURN QUERY SELECT h3_polygon_to_cells( - ST_Buffer(buffered::geography, buffer * 1.3), + ST_Buffer(buffered::geography, buffer), resolution); ELSE RETURN QUERY @@ -111,8 +111,8 @@ AS $$ rast, poly, resolution, - h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3) - AS h3 + h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3 + ) AS h3 ) c) SELECT g.* FROM geoms g @@ -157,19 +157,21 @@ CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_parts( nband integer) RETURNS TABLE (h3 h3index, part raster) AS $$ +DECLARE + nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); +BEGIN + RETURN QUERY WITH parts AS ( SELECT - h3, - ST_Clip(rast, nband, geom, __h3_raster_band_nodata(rast, nband)) AS part - FROM ( - -- h3, geom - SELECT (__h3_raster_polygon_to_cell_boundaries_intersects(rast, poly, resolution)).* - ) t) - SELECT h3, part - FROM parts - WHERE NOT ST_BandIsNoData(part, nband); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + g.h3, + ST_Clip(rast, nband, g.geom, nodata, TRUE) AS part + FROM __h3_raster_polygon_to_cell_boundaries_intersects(rast, poly, resolution) g) + SELECT p.h3, p.part + FROM parts p + WHERE NOT ST_BandIsNoData(p.part, nband); +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; -- Get values corresponding to all H3 cells with centroids inside the -- raster polygon. Assumes cell area is less than pixel area. @@ -183,10 +185,7 @@ AS $$ SELECT h3, ST_Value(rast, nband, x, y) AS val - FROM ( - -- h3, x, y - SELECT (__h3_raster_polygon_to_cell_coords_centroid(rast, poly, resolution)).* - ) t; + FROM __h3_raster_polygon_to_cell_coords_centroid(rast, poly, resolution); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --| ## Continuous raster data @@ -273,10 +272,7 @@ AS $$ SELECT h3, __h3_raster_to_summary_stats(ST_SummaryStats(part, nband, TRUE)) AS stats - FROM ( - -- h3, part - SELECT (__h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband)).* - ) t; + FROM __h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --@ availability: unreleased @@ -286,12 +282,11 @@ CREATE OR REPLACE FUNCTION h3_raster_summary_clip( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ - SELECT (__h3_raster_polygon_summary_clip( + SELECT __h3_raster_polygon_summary_clip( rast, __h3_raster_to_polygon(rast, nband), resolution, - nband - )).*; + nband); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_summary_clip(raster, integer, integer) @@ -314,10 +309,7 @@ AS $$ min(val), max(val) )::h3_raster_summary_stats AS stats - FROM ( - -- x, y, val, geom - SELECT (ST_PixelAsCentroids(rast, nband)).* - ) t + FROM ST_PixelAsCentroids(rast, nband) GROUP BY 1; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION @@ -342,10 +334,7 @@ AS $$ val, -- min val -- max )::h3_raster_summary_stats AS stats - FROM ( - -- h3, val - SELECT (__h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband)).* - ) t; + FROM __h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband) t; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --@ availability: unreleased @@ -365,8 +354,7 @@ BEGIN poly, resolution, nband, - cell_area / pixel_area) - ).*; + cell_area / pixel_area)).*; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION @@ -386,7 +374,7 @@ DECLARE __h3_raster_polygon_centroid_cell_area(poly, resolution) / __h3_raster_pixel_area(rast); BEGIN - IF pixels_per_cell > 350 THEN + IF pixels_per_cell > 100 THEN RETURN QUERY SELECT (__h3_raster_polygon_summary_clip( rast, poly, @@ -406,7 +394,7 @@ BEGIN resolution, nband, pixels_per_cell - )).*; + )).*; END IF; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; @@ -422,11 +410,9 @@ IS 'Returns `h3_raster_summary_stats` for each H3 cell in raster for a given ban --| WITH --| summary AS ( --| SELECT h3, val, h3_raster_class_summary_item_agg(summary) AS item ---| FROM ( ---| -- h3, val, summary ---| SELECT (h3_raster_class_summary(rast, 8)).* AS summary ---| FROM rasters ---| ) t +--| FROM +--| rasters, +--| 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 @@ -474,7 +460,7 @@ AS $$ SELECT s1.val, s1.count + s2.count, - s1.area + s2.area + s1.area + s2.area; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --@ availability: unreleased @@ -490,15 +476,12 @@ CREATE OR REPLACE FUNCTION __h3_raster_class_summary_part( pixel_area double precision) RETURNS SETOF h3_raster_class_summary_item AS $$ - WITH - vals AS (SELECT unnest(ST_DumpValues(rast, nband)) AS val) - SELECT - vals.val::integer, - count(*)::double precision, - count(*) * pixel_area - FROM vals - WHERE val IS NOT NULL - GROUP BY 1 + SELECT ROW( + value::integer, + count::double precision, + count * pixel_area + )::h3_raster_class_summary_item + FROM ST_ValueCount(rast, nband) t; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_class_polygon_summary_clip( @@ -514,10 +497,7 @@ AS $$ SELECT h3, __h3_raster_class_summary_part(part, nband, pixel_area) AS summary - FROM ( - -- h3, part - SELECT (__h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband)).* - ) t) + FROM __h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband) t) SELECT h3, (summary).val, summary FROM summary; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -529,13 +509,13 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ - SELECT (__h3_raster_class_polygon_summary_clip( + SELECT __h3_raster_class_polygon_summary_clip( rast, __h3_raster_to_polygon(rast, nband), resolution, nband, __h3_raster_pixel_area(rast) - )).* + ); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_clip(raster, integer, integer) @@ -556,10 +536,7 @@ AS $$ count(*)::double precision, count(*) * pixel_area )::h3_raster_class_summary_item AS summary - FROM ( - -- x, y, val, geom - SELECT (ST_PixelAsCentroids(rast, nband)).* - ) c + FROM ST_PixelAsCentroids(rast, nband) GROUP BY 1, 2; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -570,12 +547,11 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_centroids( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ - SELECT (__h3_raster_class_summary_centroids( + SELECT __h3_raster_class_summary_centroids( rast, resolution, nband, - __h3_raster_pixel_area(rast) - )).* + __h3_raster_pixel_area(rast)); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_centroids(raster, integer, integer) @@ -598,10 +574,7 @@ AS $$ cell_area / pixel_area, cell_area )::h3_raster_class_summary_item AS summary - FROM ( - -- h3, val - SELECT (__h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband)).* - ) t; + FROM __h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; --@ availability: unreleased @@ -641,7 +614,7 @@ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN - IF pixels_per_cell > 350 THEN + IF pixels_per_cell > 100 THEN RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( rast, poly, diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index 7069e228..42a708af 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -97,7 +97,7 @@ BEGIN IF buffer > 0.0 THEN RETURN QUERY SELECT h3_polygon_to_cells( - ST_Buffer(buffered::geography, buffer * 1.3), + ST_Buffer(buffered::geography, buffer), resolution); ELSE RETURN QUERY @@ -122,8 +122,8 @@ AS $$ rast, poly, resolution, - h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3) - AS h3 + h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3 + ) AS h3 ) c) SELECT g.* FROM geoms g @@ -167,19 +167,21 @@ CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_parts( nband integer) RETURNS TABLE (h3 h3index, part raster) AS $$ +DECLARE + nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); +BEGIN + RETURN QUERY WITH parts AS ( SELECT - h3, - ST_Clip(rast, nband, geom, __h3_raster_band_nodata(rast, nband)) AS part - FROM ( - -- h3, geom - SELECT (__h3_raster_polygon_to_cell_boundaries_intersects(rast, poly, resolution)).* - ) t) - SELECT h3, part - FROM parts - WHERE NOT ST_BandIsNoData(part, nband); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + g.h3, + ST_Clip(rast, nband, g.geom, nodata, TRUE) AS part + FROM __h3_raster_polygon_to_cell_boundaries_intersects(rast, poly, resolution) g) + SELECT p.h3, p.part + FROM parts p + WHERE NOT ST_BandIsNoData(p.part, nband); +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_subpixel_cell_values( rast raster, @@ -191,10 +193,7 @@ AS $$ SELECT h3, ST_Value(rast, nband, x, y) AS val - FROM ( - -- h3, x, y - SELECT (__h3_raster_polygon_to_cell_coords_centroid(rast, poly, resolution)).* - ) t; + FROM __h3_raster_polygon_to_cell_coords_centroid(rast, poly, resolution); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; -- Raster processing: continuous data @@ -265,10 +264,7 @@ AS $$ SELECT h3, __h3_raster_to_summary_stats(ST_SummaryStats(part, nband, TRUE)) AS stats - FROM ( - -- h3, part - SELECT (__h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband)).* - ) t; + FROM __h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_summary_clip( @@ -277,12 +273,11 @@ CREATE OR REPLACE FUNCTION h3_raster_summary_clip( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ - SELECT (__h3_raster_polygon_summary_clip( + SELECT __h3_raster_polygon_summary_clip( rast, __h3_raster_to_polygon(rast, nband), resolution, - nband - )).*; + nband); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_summary_clip(raster, integer, integer) @@ -304,17 +299,13 @@ AS $$ min(val), max(val) )::h3_raster_summary_stats AS stats - FROM ( - -- x, y, val, geom - SELECT (ST_PixelAsCentroids(rast, nband)).* - ) t + FROM ST_PixelAsCentroids(rast, nband) GROUP BY 1; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_summary_centroids(raster, integer, integer) IS '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.'; - CREATE OR REPLACE FUNCTION __h3_raster_polygon_summary_subpixel( rast raster, poly geometry, @@ -333,10 +324,7 @@ AS $$ val, -- min val -- max )::h3_raster_summary_stats AS stats - FROM ( - -- h3, val - SELECT (__h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband)).* - ) t; + FROM __h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband) t; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_summary_subpixel( @@ -355,8 +343,7 @@ BEGIN poly, resolution, nband, - cell_area / pixel_area) - ).*; + cell_area / pixel_area)).*; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION @@ -375,7 +362,7 @@ DECLARE __h3_raster_polygon_centroid_cell_area(poly, resolution) / __h3_raster_pixel_area(rast); BEGIN - IF pixels_per_cell > 350 THEN + IF pixels_per_cell > 100 THEN RETURN QUERY SELECT (__h3_raster_polygon_summary_clip( rast, poly, @@ -395,7 +382,7 @@ BEGIN resolution, nband, pixels_per_cell - )).*; + )).*; END IF; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; @@ -433,7 +420,7 @@ AS $$ SELECT s1.val, s1.count + s2.count, - s1.area + s2.area + s1.area + s2.area; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE AGGREGATE h3_raster_class_summary_item_agg(h3_raster_class_summary_item) ( @@ -448,15 +435,12 @@ CREATE OR REPLACE FUNCTION __h3_raster_class_summary_part( pixel_area double precision) RETURNS SETOF h3_raster_class_summary_item AS $$ - WITH - vals AS (SELECT unnest(ST_DumpValues(rast, nband)) AS val) - SELECT - vals.val::integer, - count(*)::double precision, - count(*) * pixel_area - FROM vals - WHERE val IS NOT NULL - GROUP BY 1 + SELECT ROW( + value::integer, + count::double precision, + count * pixel_area + )::h3_raster_class_summary_item + FROM ST_ValueCount(rast, nband) t; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_class_polygon_summary_clip( @@ -472,10 +456,7 @@ AS $$ SELECT h3, __h3_raster_class_summary_part(part, nband, pixel_area) AS summary - FROM ( - -- h3, part - SELECT (__h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband)).* - ) t) + FROM __h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband) t) SELECT h3, (summary).val, summary FROM summary; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -486,13 +467,13 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ - SELECT (__h3_raster_class_polygon_summary_clip( + SELECT __h3_raster_class_polygon_summary_clip( rast, __h3_raster_to_polygon(rast, nband), resolution, nband, __h3_raster_pixel_area(rast) - )).* + ); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_clip(raster, integer, integer) @@ -513,10 +494,7 @@ AS $$ count(*)::double precision, count(*) * pixel_area )::h3_raster_class_summary_item AS summary - FROM ( - -- x, y, val, geom - SELECT (ST_PixelAsCentroids(rast, nband)).* - ) c + FROM ST_PixelAsCentroids(rast, nband) GROUP BY 1, 2; $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; @@ -526,12 +504,11 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_centroids( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ - SELECT (__h3_raster_class_summary_centroids( + SELECT __h3_raster_class_summary_centroids( rast, resolution, nband, - __h3_raster_pixel_area(rast) - )).* + __h3_raster_pixel_area(rast)); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_centroids(raster, integer, integer) @@ -554,10 +531,7 @@ AS $$ cell_area / pixel_area, cell_area )::h3_raster_class_summary_item AS summary - FROM ( - -- h3, val - SELECT (__h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband)).* - ) t; + FROM __h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( @@ -595,7 +569,7 @@ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN - IF pixels_per_cell > 350 THEN + IF pixels_per_cell > 100 THEN RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( rast, poly, From b7e2fb4058172136b717738f1c26c513b2d72715 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Thu, 5 Jan 2023 01:23:11 +0300 Subject: [PATCH 11/22] PostGIS Raster integration: documentation update --- docs/api.md | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/docs/api.md b/docs/api.md index ddafd071..06b2a9f5 100644 --- a/docs/api.md +++ b/docs/api.md @@ -676,11 +676,9 @@ adding `fraction` value (fraction of H3 cell area for each value): WITH summary AS ( SELECT h3, val, h3_raster_class_summary_item_agg(summary) AS item - FROM ( - -- h3, val, summary - SELECT (h3_raster_class_summary(rast, 8)).* AS summary - FROM rasters - ) t + FROM + rasters, + 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 ff8bd856e63f2515b454642d9dd6bcfbf5a1cdf4 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Tue, 10 Jan 2023 00:34:59 +0300 Subject: [PATCH 12/22] PostGIS Raster integration: method selection based on pixels per cell and cells per raster --- h3_postgis/sql/install/40-rasters.sql | 18 +++++++++++++----- .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 18 +++++++++++++----- 2 files changed, 26 insertions(+), 10 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 3bdbd387..ff685c7a 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -370,11 +370,15 @@ RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); - pixels_per_cell CONSTANT double precision := - __h3_raster_polygon_centroid_cell_area(poly, resolution) - / __h3_raster_pixel_area(rast); + cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN - IF pixels_per_cell > 100 THEN + IF pixels_per_cell > 250 + OR pixels_per_cell > 80 + -- cells_per_raster > 10000 / (pixels_per_cell - 80) + AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 80) + THEN RETURN QUERY SELECT (__h3_raster_polygon_summary_clip( rast, poly, @@ -614,7 +618,11 @@ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN - IF pixels_per_cell > 100 THEN + IF pixels_per_cell > 250 + OR pixels_per_cell > 80 + -- cells_per_raster > 10000 / (pixels_per_cell - 80) + AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 80) + THEN RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( rast, poly, diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index 42a708af..c657d655 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -358,11 +358,15 @@ RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); - pixels_per_cell CONSTANT double precision := - __h3_raster_polygon_centroid_cell_area(poly, resolution) - / __h3_raster_pixel_area(rast); + cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); + pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN - IF pixels_per_cell > 100 THEN + IF pixels_per_cell > 250 + OR pixels_per_cell > 80 + -- cells_per_raster > 10000 / (pixels_per_cell - 80) + AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 80) + THEN RETURN QUERY SELECT (__h3_raster_polygon_summary_clip( rast, poly, @@ -569,7 +573,11 @@ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN - IF pixels_per_cell > 100 THEN + IF pixels_per_cell > 250 + OR pixels_per_cell > 80 + -- cells_per_raster > 10000 / (pixels_per_cell - 80) + AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 80) + THEN RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( rast, poly, From de88f47f35a1b01d264dba977b9dcaddcc2da377 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Thu, 12 Jan 2023 14:09:50 +0300 Subject: [PATCH 13/22] PostGIS Raster integration: STRICT removed for helper functions that are potentially inlinable --- h3_postgis/sql/install/40-rasters.sql | 30 +++++++++---------- .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 30 +++++++++---------- 2 files changed, 30 insertions(+), 30 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index ff685c7a..86ab61a4 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -28,7 +28,7 @@ AS $$ SELECT coalesce( ST_BandNoDataValue(rast, nband), ST_MinPossibleValue(ST_BandPixelType(rast, nband))); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_to_polygon( rast raster, @@ -36,14 +36,14 @@ CREATE OR REPLACE FUNCTION __h3_raster_to_polygon( RETURNS geometry AS $$ SELECT ST_MinConvexHull(rast, nband); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; -- Area of a single pixel in meters CREATE OR REPLACE FUNCTION __h3_raster_pixel_area(rast raster) RETURNS double precision AS $$ SELECT ST_Area(ST_Transform(ST_PixelAsPolygon(rast, 1, 1), 4326)::geography); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; -- Area of a cell close to the center of raster polygon, in meters CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( @@ -52,7 +52,7 @@ CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( RETURNS h3index AS $$ SELECT h3_lat_lng_to_cell(ST_Transform(ST_Centroid(poly), 4326), resolution); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell_area( poly geometry, @@ -62,7 +62,7 @@ AS $$ SELECT ST_Area( h3_cell_to_boundary_geography( __h3_raster_polygon_centroid_cell(poly, resolution))); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; -- Get list of cells inside of the raster polygon, -- buffered by `buffer` value (in meters). @@ -117,7 +117,7 @@ AS $$ SELECT g.* FROM geoms g WHERE ST_Intersects(g.geom, poly); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; -- Get raster coordinates of H3 cells with centroids inside the raster polygon CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_coords_centroid( @@ -148,7 +148,7 @@ AS $$ WHERE x BETWEEN 1 AND ST_Width(rast) AND y BETWEEN 1 AND ST_Height(rast); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_parts( rast raster, @@ -186,7 +186,7 @@ AS $$ h3, ST_Value(rast, nband, x, y) AS val FROM __h3_raster_polygon_to_cell_coords_centroid(rast, poly, resolution); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; --| ## Continuous raster data --| @@ -225,7 +225,7 @@ AS $$ (stats).min, (stats).max )::h3_raster_summary_stats -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_summary_stats_agg_transfn( s1 h3_raster_summary_stats, @@ -273,7 +273,7 @@ AS $$ h3, __h3_raster_to_summary_stats(ST_SummaryStats(part, nband, TRUE)) AS stats FROM __h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; --@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_summary_clip( @@ -335,7 +335,7 @@ AS $$ val -- max )::h3_raster_summary_stats AS stats FROM __h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband) t; -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; --@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_summary_subpixel( @@ -486,7 +486,7 @@ AS $$ count * pixel_area )::h3_raster_class_summary_item FROM ST_ValueCount(rast, nband) t; -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_class_polygon_summary_clip( rast raster, @@ -504,7 +504,7 @@ AS $$ FROM __h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband) t) SELECT h3, (summary).val, summary FROM summary; -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; --@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( @@ -542,7 +542,7 @@ AS $$ )::h3_raster_class_summary_item AS summary FROM ST_PixelAsCentroids(rast, nband) GROUP BY 1, 2; -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; --@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_class_summary_centroids( @@ -579,7 +579,7 @@ AS $$ cell_area )::h3_raster_class_summary_item AS summary FROM __h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; --@ availability: unreleased CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index c657d655..0490ad98 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -45,7 +45,7 @@ AS $$ SELECT coalesce( ST_BandNoDataValue(rast, nband), ST_MinPossibleValue(ST_BandPixelType(rast, nband))); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_to_polygon( rast raster, @@ -53,13 +53,13 @@ CREATE OR REPLACE FUNCTION __h3_raster_to_polygon( RETURNS geometry AS $$ SELECT ST_MinConvexHull(rast, nband); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_pixel_area(rast raster) RETURNS double precision AS $$ SELECT ST_Area(ST_Transform(ST_PixelAsPolygon(rast, 1, 1), 4326)::geography); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( poly geometry, @@ -67,7 +67,7 @@ CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( RETURNS h3index AS $$ SELECT h3_lat_lng_to_cell(ST_Transform(ST_Centroid(poly), 4326), resolution); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell_area( poly geometry, @@ -77,7 +77,7 @@ AS $$ SELECT ST_Area( h3_cell_to_boundary_geography( __h3_raster_polygon_centroid_cell(poly, resolution))); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cells( rast raster, @@ -128,7 +128,7 @@ AS $$ SELECT g.* FROM geoms g WHERE ST_Intersects(g.geom, poly); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_coords_centroid( rast raster, @@ -158,7 +158,7 @@ AS $$ WHERE x BETWEEN 1 AND ST_Width(rast) AND y BETWEEN 1 AND ST_Height(rast); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_parts( rast raster, @@ -194,7 +194,7 @@ AS $$ h3, ST_Value(rast, nband, x, y) AS val FROM __h3_raster_polygon_to_cell_coords_centroid(rast, poly, resolution); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; -- Raster processing: continuous data @@ -218,7 +218,7 @@ AS $$ (stats).min, (stats).max )::h3_raster_summary_stats -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_summary_stats_agg_transfn( s1 h3_raster_summary_stats, @@ -265,7 +265,7 @@ AS $$ h3, __h3_raster_to_summary_stats(ST_SummaryStats(part, nband, TRUE)) AS stats FROM __h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_summary_clip( rast raster, @@ -325,7 +325,7 @@ AS $$ val -- max )::h3_raster_summary_stats AS stats FROM __h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband) t; -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_summary_subpixel( rast raster, @@ -445,7 +445,7 @@ AS $$ count * pixel_area )::h3_raster_class_summary_item FROM ST_ValueCount(rast, nband) t; -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_class_polygon_summary_clip( rast raster, @@ -463,7 +463,7 @@ AS $$ FROM __h3_raster_polygon_to_cell_parts(rast, poly, resolution, nband) t) SELECT h3, (summary).val, summary FROM summary; -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( rast raster, @@ -500,7 +500,7 @@ AS $$ )::h3_raster_class_summary_item AS summary FROM ST_PixelAsCentroids(rast, nband) GROUP BY 1, 2; -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_class_summary_centroids( rast raster, @@ -536,7 +536,7 @@ AS $$ cell_area )::h3_raster_class_summary_item AS summary FROM __h3_raster_polygon_subpixel_cell_values(rast, poly, resolution, nband); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; +$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION h3_raster_class_summary_subpixel( rast raster, From 51142e7585e34f33bca543ff862929b2f299c9df Mon Sep 17 00:00:00 2001 From: sergei sh Date: Thu, 12 Jan 2023 15:50:34 +0300 Subject: [PATCH 14/22] PostGIS Raster integration: ST_Buffer "join=mitre" parameter added for simpler result geometry --- h3_postgis/sql/install/40-rasters.sql | 10 ++++++++-- h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql | 10 ++++++++-- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 86ab61a4..edb2d522 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -79,13 +79,19 @@ DECLARE BEGIN IF ST_SRID(rast) != 4326 THEN buffered := ST_Transform( - ST_Buffer(poly, greatest(ST_PixelWidth(rast), ST_PixelHeight(rast))), + ST_Buffer( + poly, + greatest(ST_PixelWidth(rast), ST_PixelHeight(rast)), + 'join=mitre'), 4326); END IF; IF buffer > 0.0 THEN RETURN QUERY SELECT h3_polygon_to_cells( - ST_Buffer(buffered::geography, buffer), + ST_Buffer( + buffered::geography, + buffer, + 'join=mitre'), resolution); ELSE RETURN QUERY diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index 0490ad98..fef51c05 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -91,13 +91,19 @@ DECLARE BEGIN IF ST_SRID(rast) != 4326 THEN buffered := ST_Transform( - ST_Buffer(poly, greatest(ST_PixelWidth(rast), ST_PixelHeight(rast))), + ST_Buffer( + poly, + greatest(ST_PixelWidth(rast), ST_PixelHeight(rast)), + 'join=mitre'), 4326); END IF; IF buffer > 0.0 THEN RETURN QUERY SELECT h3_polygon_to_cells( - ST_Buffer(buffered::geography, buffer), + ST_Buffer( + buffered::geography, + buffer, + 'join=mitre'), resolution); ELSE RETURN QUERY From 7c753f71367174bccc2156e0f85b69f15c9c83c5 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Fri, 13 Jan 2023 12:18:29 +0300 Subject: [PATCH 15/22] PostGIS Raster integration: cell geometry filtering fix --- h3_postgis/sql/install/40-rasters.sql | 40 +++++++------------ .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 40 +++++++------------ 2 files changed, 30 insertions(+), 50 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index edb2d522..92a72759 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -107,22 +107,16 @@ CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_boundaries_intersects( resolution integer) RETURNS TABLE (h3 h3index, geom geometry) AS $$ - WITH - geoms AS ( - SELECT - c.h3, - ST_Transform(h3_cell_to_boundary_geometry(c.h3), ST_SRID(rast)) AS geom - FROM ( - SELECT __h3_raster_polygon_to_cells( - rast, - poly, - resolution, - h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3 - ) AS h3 - ) c) - SELECT g.* - FROM geoms g - WHERE ST_Intersects(g.geom, poly); + SELECT h3, geom + FROM + __h3_raster_polygon_to_cells( + rast, + poly, + resolution, + h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3 + ) AS h3, + ST_Transform(h3_cell_to_boundary_geometry(h3), ST_SRID(rast)) AS geom + WHERE ST_Intersects(geom, poly); $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; -- Get raster coordinates of H3 cells with centroids inside the raster polygon @@ -167,15 +161,11 @@ DECLARE nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); BEGIN RETURN QUERY - WITH - parts AS ( - SELECT - g.h3, - ST_Clip(rast, nband, g.geom, nodata, TRUE) AS part - FROM __h3_raster_polygon_to_cell_boundaries_intersects(rast, poly, resolution) g) - SELECT p.h3, p.part - FROM parts p - WHERE NOT ST_BandIsNoData(p.part, nband); + SELECT c.h3, p AS part + FROM + __h3_raster_polygon_to_cell_boundaries_intersects(rast, poly, resolution) AS c, + ST_Clip(rast, nband, c.geom, nodata, TRUE) AS p + WHERE NOT ST_BandIsNoData(p, nband); END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index fef51c05..d8f6d44e 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -118,22 +118,16 @@ CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_boundaries_intersects( resolution integer) RETURNS TABLE (h3 h3index, geom geometry) AS $$ - WITH - geoms AS ( - SELECT - c.h3, - ST_Transform(h3_cell_to_boundary_geometry(c.h3), ST_SRID(rast)) AS geom - FROM ( - SELECT __h3_raster_polygon_to_cells( - rast, - poly, - resolution, - h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3 - ) AS h3 - ) c) - SELECT g.* - FROM geoms g - WHERE ST_Intersects(g.geom, poly); + SELECT h3, geom + FROM + __h3_raster_polygon_to_cells( + rast, + poly, + resolution, + h3_get_hexagon_edge_length_avg(resolution, 'm') * 1.3 + ) AS h3, + ST_Transform(h3_cell_to_boundary_geometry(h3), ST_SRID(rast)) AS geom + WHERE ST_Intersects(geom, poly); $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_to_cell_coords_centroid( @@ -177,15 +171,11 @@ DECLARE nodata CONSTANT double precision := __h3_raster_band_nodata(rast, nband); BEGIN RETURN QUERY - WITH - parts AS ( - SELECT - g.h3, - ST_Clip(rast, nband, g.geom, nodata, TRUE) AS part - FROM __h3_raster_polygon_to_cell_boundaries_intersects(rast, poly, resolution) g) - SELECT p.h3, p.part - FROM parts p - WHERE NOT ST_BandIsNoData(p.part, nband); + SELECT c.h3, p AS part + FROM + __h3_raster_polygon_to_cell_boundaries_intersects(rast, poly, resolution) AS c, + ST_Clip(rast, nband, c.geom, nodata, TRUE) AS p + WHERE NOT ST_BandIsNoData(p, nband); END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; From 4f4bfb364520e5a6369ba6a3928ea9b1b478d8a1 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Fri, 13 Jan 2023 13:00:55 +0300 Subject: [PATCH 16/22] PostGIS Raster integration: method selection condition adjusted --- h3_postgis/sql/install/40-rasters.sql | 12 ++++-------- h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql | 12 ++++-------- 2 files changed, 8 insertions(+), 16 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 92a72759..4f19a68c 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -370,10 +370,8 @@ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN - IF pixels_per_cell > 250 - OR pixels_per_cell > 80 - -- cells_per_raster > 10000 / (pixels_per_cell - 80) - AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 80) + IF pixels_per_cell > 70 + AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 70) THEN RETURN QUERY SELECT (__h3_raster_polygon_summary_clip( rast, @@ -614,10 +612,8 @@ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN - IF pixels_per_cell > 250 - OR pixels_per_cell > 80 - -- cells_per_raster > 10000 / (pixels_per_cell - 80) - AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 80) + IF pixels_per_cell > 70 + AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 70) THEN RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( rast, diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index d8f6d44e..bc50a72b 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -358,10 +358,8 @@ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN - IF pixels_per_cell > 250 - OR pixels_per_cell > 80 - -- cells_per_raster > 10000 / (pixels_per_cell - 80) - AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 80) + IF pixels_per_cell > 70 + AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 70) THEN RETURN QUERY SELECT (__h3_raster_polygon_summary_clip( rast, @@ -569,10 +567,8 @@ DECLARE pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN - IF pixels_per_cell > 250 - OR pixels_per_cell > 80 - -- cells_per_raster > 10000 / (pixels_per_cell - 80) - AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 80) + IF pixels_per_cell > 70 + AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 70) THEN RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( rast, From 42625ac70d5ee980822d6a03259f94936adc69ad Mon Sep 17 00:00:00 2001 From: sergei sh Date: Fri, 13 Jan 2023 14:15:24 +0300 Subject: [PATCH 17/22] PostGIS Raster integration: always using hexagons for cell area estimation --- h3_postgis/sql/install/40-rasters.sql | 13 ++++++++++--- h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql | 11 +++++++++-- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 4f19a68c..4f1cd13e 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -45,15 +45,22 @@ AS $$ SELECT ST_Area(ST_Transform(ST_PixelAsPolygon(rast, 1, 1), 4326)::geography); $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; --- Area of a cell close to the center of raster polygon, in meters CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( poly geometry, resolution integer) RETURNS h3index AS $$ - SELECT h3_lat_lng_to_cell(ST_Transform(ST_Centroid(poly), 4326), resolution); -$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; +DECLARE + cell h3index := h3_lat_lng_to_cell(ST_Transform(ST_Centroid(poly), 4326), resolution); +BEGIN + IF h3_is_pentagon(cell) THEN + SELECT h3 INTO cell FROM h3_grid_disk(cell) AS h3 WHERE h3 != cell LIMIT 1; + END IF; + RETURN cell; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; +-- Area of a cell close to the center of raster polygon, in meters CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell_area( poly geometry, resolution integer) diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index bc50a72b..5e876e25 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -66,8 +66,15 @@ CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( resolution integer) RETURNS h3index AS $$ - SELECT h3_lat_lng_to_cell(ST_Transform(ST_Centroid(poly), 4326), resolution); -$$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; +DECLARE + cell h3index := h3_lat_lng_to_cell(ST_Transform(ST_Centroid(poly), 4326), resolution); +BEGIN + IF h3_is_pentagon(cell) THEN + SELECT h3 INTO cell FROM h3_grid_disk(cell) AS h3 WHERE h3 != cell LIMIT 1; + END IF; + RETURN cell; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell_area( poly geometry, From 18ea23b7127011ff19e87eeaf0ef6a0cb822dd8a Mon Sep 17 00:00:00 2001 From: sergei sh Date: Tue, 17 Jan 2023 18:59:08 +0300 Subject: [PATCH 18/22] PostGIS Raster integration: documentation update --- docs/api.md | 34 +++++++++++++--- h3_postgis/sql/install/40-rasters.sql | 40 +++++++++++++++---- .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 4 +- 3 files changed, 63 insertions(+), 15 deletions(-) diff --git a/docs/api.md b/docs/api.md index 06b2a9f5..583bdbef 100644 --- a/docs/api.md +++ b/docs/api.md @@ -621,7 +621,10 @@ Splits polygons when crossing 180th meridian. # Raster processing functions ## Continuous raster data -Combining summary stats from multiple rasters: +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, @@ -631,6 +634,12 @@ FROM ( 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 +<...> ``` @@ -670,8 +679,14 @@ Returns `h3_raster_summary_stats` for each H3 cell in raster for a given band. A ## Discrete raster data -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): +For rasters where pixels have discrete values corresponding to different classes +of land cover or land use, H3 cell 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 ( @@ -681,18 +696,25 @@ WITH 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) + 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}} +<...> ``` @@ -703,7 +725,7 @@ GROUP BY 1; *Since vunreleased* -Convert raster summary to binary JSON. +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`) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 4f1cd13e..02ea9891 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -1,5 +1,5 @@ /* - * Copyright 2022 Bytes & Brains + * Copyright 2022-2023 Kontur Inc. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. @@ -193,7 +193,11 @@ $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; --| ## Continuous raster data --| ---| Combining summary stats from multiple rasters: +--| 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, @@ -203,6 +207,13 @@ $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; --| 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 +--| <...> --| ``` -- NOTE: `count` can be < 1 when cell area is less than pixel area @@ -409,8 +420,15 @@ IS 'Returns `h3_raster_summary_stats` for each H3 cell in raster for a given ban --| ## Discrete raster data --| ---| 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): +--| For rasters where pixels have discrete values corresponding to different classes +--| of land cover or land use, H3 cell 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 ( @@ -420,18 +438,26 @@ IS 'Returns `h3_raster_summary_stats` for each H3 cell in raster for a given ban --| 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) +--| 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}} +--| <...> --| ``` -- NOTE: `count` can be < 1 when cell area is less than pixel area @@ -455,7 +481,7 @@ AS $$ $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_item_to_jsonb(h3_raster_class_summary_item) -IS 'Convert raster summary to binary JSON.'; +IS 'Convert raster summary to JSONB, example: `{"count": 10, "value": 2, "area": 16490.3423}`'; CREATE OR REPLACE FUNCTION __h3_raster_class_summary_item_agg_transfn( s1 h3_raster_class_summary_item, diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index 5e876e25..4905e8e9 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -1,5 +1,5 @@ /* - * Copyright 2022 Bytes & Brains + * Copyright 2022-2023 Kontur Inc. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. @@ -415,7 +415,7 @@ AS $$ $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_item_to_jsonb(h3_raster_class_summary_item) -IS 'Convert raster summary to binary JSON.'; +IS 'Convert raster summary to JSONB, example: `{"count": 10, "value": 2, "area": 16490.3423}`'; CREATE OR REPLACE FUNCTION __h3_raster_class_summary_item_agg_transfn( s1 h3_raster_class_summary_item, From 711d80cad790df3788051d89be0ab894e1910cfc Mon Sep 17 00:00:00 2001 From: sergei sh Date: Tue, 17 Jan 2023 22:03:22 +0300 Subject: [PATCH 19/22] PostGIS Raster integration: pixel area estimation fix --- h3_postgis/sql/install/40-rasters.sql | 38 ++++++++++++------- .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 37 ++++++++++++------ 2 files changed, 50 insertions(+), 25 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 02ea9891..414cb02d 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -38,11 +38,18 @@ AS $$ SELECT ST_MinConvexHull(rast, nband); $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; --- Area of a single pixel in meters -CREATE OR REPLACE FUNCTION __h3_raster_pixel_area(rast raster) +-- Area of a pixel close to the center of raster polygon, in meters +CREATE OR REPLACE FUNCTION __h3_raster_polygon_pixel_area( + rast raster, + poly geometry) RETURNS double precision AS $$ - SELECT ST_Area(ST_Transform(ST_PixelAsPolygon(rast, 1, 1), 4326)::geography); + SELECT ST_Area( + ST_PixelAsPolygon( + rast, + ST_WorldToRasterCoordX(rast, c), + ST_WorldToRasterCoordY(rast, c))::geography) + FROM ST_Transform(ST_Centroid(poly), 4326) AS c $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( @@ -360,7 +367,7 @@ RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + pixel_area CONSTANT double precision := __h3_raster_polygon_pixel_area(rast, poly); cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); BEGIN RETURN QUERY SELECT (__h3_raster_polygon_summary_subpixel( @@ -385,7 +392,7 @@ AS $$ DECLARE poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + pixel_area CONSTANT double precision := __h3_raster_polygon_pixel_area(rast, poly); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 70 @@ -540,14 +547,18 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ - SELECT __h3_raster_class_polygon_summary_clip( +DECLARE + poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); +BEGIN + RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( rast, - __h3_raster_to_polygon(rast, nband), + poly, resolution, nband, - __h3_raster_pixel_area(rast) - ); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + __h3_raster_polygon_pixel_area(rast, poly) + )).*; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_clip(raster, integer, integer) IS '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.'; @@ -582,7 +593,8 @@ AS $$ rast, resolution, nband, - __h3_raster_pixel_area(rast)); + __h3_raster_polygon_pixel_area(rast, __h3_raster_to_polygon(rast, nband)) + ); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_centroids(raster, integer, integer) @@ -624,7 +636,7 @@ BEGIN resolution, nband, __h3_raster_polygon_centroid_cell_area(poly, resolution), - __h3_raster_pixel_area(rast) + __h3_raster_polygon_pixel_area(rast, poly) )).*; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; @@ -642,7 +654,7 @@ AS $$ DECLARE poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + pixel_area CONSTANT double precision := __h3_raster_polygon_pixel_area(rast, poly); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 70 diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index 4905e8e9..726b0a70 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -55,10 +55,17 @@ AS $$ SELECT ST_MinConvexHull(rast, nband); $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; -CREATE OR REPLACE FUNCTION __h3_raster_pixel_area(rast raster) +CREATE OR REPLACE FUNCTION __h3_raster_polygon_pixel_area( + rast raster, + poly geometry) RETURNS double precision AS $$ - SELECT ST_Area(ST_Transform(ST_PixelAsPolygon(rast, 1, 1), 4326)::geography); + SELECT ST_Area( + ST_PixelAsPolygon( + rast, + ST_WorldToRasterCoordX(rast, c), + ST_WorldToRasterCoordY(rast, c))::geography) + FROM ST_Transform(ST_Centroid(poly), 4326) AS c $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( @@ -338,7 +345,7 @@ RETURNS TABLE (h3 h3index, stats h3_raster_summary_stats) AS $$ DECLARE poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + pixel_area CONSTANT double precision := __h3_raster_polygon_pixel_area(rast, poly); cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); BEGIN RETURN QUERY SELECT (__h3_raster_polygon_summary_subpixel( @@ -362,7 +369,7 @@ AS $$ DECLARE poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + pixel_area CONSTANT double precision := __h3_raster_polygon_pixel_area(rast, poly); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 70 @@ -472,14 +479,18 @@ CREATE OR REPLACE FUNCTION h3_raster_class_summary_clip( nband integer DEFAULT 1) RETURNS TABLE (h3 h3index, val integer, summary h3_raster_class_summary_item) AS $$ - SELECT __h3_raster_class_polygon_summary_clip( +DECLARE + poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); +BEGIN + RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( rast, - __h3_raster_to_polygon(rast, nband), + poly, resolution, nband, - __h3_raster_pixel_area(rast) - ); -$$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; + __h3_raster_polygon_pixel_area(rast, poly) + )).*; +END; +$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_clip(raster, integer, integer) IS '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.'; @@ -513,7 +524,8 @@ AS $$ rast, resolution, nband, - __h3_raster_pixel_area(rast)); + __h3_raster_polygon_pixel_area(rast, __h3_raster_to_polygon(rast, nband)) + ); $$ LANGUAGE SQL IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary_centroids(raster, integer, integer) @@ -554,7 +566,7 @@ BEGIN resolution, nband, __h3_raster_polygon_centroid_cell_area(poly, resolution), - __h3_raster_pixel_area(rast) + __h3_raster_polygon_pixel_area(rast, poly) )).*; END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; @@ -571,7 +583,7 @@ AS $$ DECLARE poly CONSTANT geometry := __h3_raster_to_polygon(rast, nband); cell_area CONSTANT double precision := __h3_raster_polygon_centroid_cell_area(poly, resolution); - pixel_area CONSTANT double precision := __h3_raster_pixel_area(rast); + pixel_area CONSTANT double precision := __h3_raster_polygon_pixel_area(rast, poly); pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 70 @@ -605,3 +617,4 @@ END; $$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE; COMMENT ON FUNCTION h3_raster_class_summary(raster, integer, integer) IS '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.'; + From 1c29e8b5827ad893f9365c35e1fe48ae18ed0713 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Wed, 18 Jan 2023 14:59:26 +0300 Subject: [PATCH 20/22] PostGIS Raster integration: documentation update --- docs/api.md | 3 ++- h3_postgis/sql/install/40-rasters.sql | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/api.md b/docs/api.md index 583bdbef..a9c4788a 100644 --- a/docs/api.md +++ b/docs/api.md @@ -680,7 +680,7 @@ Returns `h3_raster_summary_stats` for each H3 cell in raster for a given band. A ## Discrete raster data For rasters where pixels have discrete values corresponding to different classes -of land cover or land use, H3 cell summary can be represented by a JSON object +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 @@ -690,6 +690,7 @@ 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, diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index 414cb02d..aa0a994c 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -428,7 +428,7 @@ IS 'Returns `h3_raster_summary_stats` for each H3 cell in raster for a given ban --| ## Discrete raster data --| --| For rasters where pixels have discrete values corresponding to different classes ---| of land cover or land use, H3 cell summary can be represented by a JSON object +--| 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 @@ -439,6 +439,7 @@ IS 'Returns `h3_raster_summary_stats` for each H3 cell in raster for a given ban --| ``` --| 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, From 9a08ab71b0cc34c1d327c1ba7311b9c8489480f3 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Wed, 18 Jan 2023 16:46:22 +0300 Subject: [PATCH 21/22] PostGIS Raster integration: documentation update --- docs/api.md | 20 ++++++++++++++++++++ h3_postgis/sql/install/40-rasters.sql | 22 ++++++++++++++++++++++ 2 files changed, 42 insertions(+) diff --git a/docs/api.md b/docs/api.md index a9c4788a..9245f522 100644 --- a/docs/api.md +++ b/docs/api.md @@ -717,6 +717,26 @@ GROUP BY 1; 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* diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index aa0a994c..a1695698 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -467,6 +467,28 @@ IS 'Returns `h3_raster_summary_stats` for each H3 cell in raster for a given ban --| 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 +--| <...> +--| ``` -- NOTE: `count` can be < 1 when cell area is less than pixel area --@ availability: unreleased From 94c4dfc3603317d1c754722f12aedcc7f55d39a1 Mon Sep 17 00:00:00 2001 From: sergei sh Date: Wed, 18 Jan 2023 17:30:21 +0300 Subject: [PATCH 22/22] PostGIS Raster integration: geometry to geography conversion fix --- h3_postgis/sql/install/40-rasters.sql | 16 +++++++++------- .../sql/updates/h3_postgis--4.0.3--4.1.0.sql | 16 +++++++++------- 2 files changed, 18 insertions(+), 14 deletions(-) diff --git a/h3_postgis/sql/install/40-rasters.sql b/h3_postgis/sql/install/40-rasters.sql index a1695698..5931a544 100644 --- a/h3_postgis/sql/install/40-rasters.sql +++ b/h3_postgis/sql/install/40-rasters.sql @@ -45,11 +45,13 @@ CREATE OR REPLACE FUNCTION __h3_raster_polygon_pixel_area( RETURNS double precision AS $$ SELECT ST_Area( - ST_PixelAsPolygon( - rast, - ST_WorldToRasterCoordX(rast, c), - ST_WorldToRasterCoordY(rast, c))::geography) - FROM ST_Transform(ST_Centroid(poly), 4326) AS c + ST_Transform( + ST_PixelAsPolygon( + rast, + ST_WorldToRasterCoordX(rast, c), + ST_WorldToRasterCoordY(rast, c)), + 4326)::geography) + FROM ST_Centroid(poly) AS c $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( @@ -396,7 +398,7 @@ DECLARE pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 70 - AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 70) + AND (ST_Area(ST_Transform(poly, 4326)::geography) / cell_area) > 10000 / (pixels_per_cell - 70) THEN RETURN QUERY SELECT (__h3_raster_polygon_summary_clip( rast, @@ -681,7 +683,7 @@ DECLARE pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 70 - AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 70) + AND (ST_Area(ST_Transform(poly, 4326)::geography) / cell_area) > 10000 / (pixels_per_cell - 70) THEN RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( rast, diff --git a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql index 726b0a70..3ee2330d 100644 --- a/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql +++ b/h3_postgis/sql/updates/h3_postgis--4.0.3--4.1.0.sql @@ -61,11 +61,13 @@ CREATE OR REPLACE FUNCTION __h3_raster_polygon_pixel_area( RETURNS double precision AS $$ SELECT ST_Area( - ST_PixelAsPolygon( - rast, - ST_WorldToRasterCoordX(rast, c), - ST_WorldToRasterCoordY(rast, c))::geography) - FROM ST_Transform(ST_Centroid(poly), 4326) AS c + ST_Transform( + ST_PixelAsPolygon( + rast, + ST_WorldToRasterCoordX(rast, c), + ST_WorldToRasterCoordY(rast, c)), + 4326)::geography) + FROM ST_Centroid(poly) AS c $$ LANGUAGE SQL IMMUTABLE PARALLEL SAFE; CREATE OR REPLACE FUNCTION __h3_raster_polygon_centroid_cell( @@ -373,7 +375,7 @@ DECLARE pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 70 - AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 70) + AND (ST_Area(ST_Transform(poly, 4326)::geography) / cell_area) > 10000 / (pixels_per_cell - 70) THEN RETURN QUERY SELECT (__h3_raster_polygon_summary_clip( rast, @@ -587,7 +589,7 @@ DECLARE pixels_per_cell CONSTANT double precision := cell_area / pixel_area; BEGIN IF pixels_per_cell > 70 - AND (ST_Area(poly::geography) / cell_area) > 10000 / (pixels_per_cell - 70) + AND (ST_Area(ST_Transform(poly, 4326)::geography) / cell_area) > 10000 / (pixels_per_cell - 70) THEN RETURN QUERY SELECT (__h3_raster_class_polygon_summary_clip( rast,