From: Pierre Racine Date: Thu, 21 Apr 2011 15:03:55 +0000 (+0000) Subject: -First version of ST_Histogram X-Git-Tag: 2.0.0alpha1~1748 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=d5407a8030fd2a38202379f24760d81084d97f1e;p=postgis -First version of ST_Histogram git-svn-id: http://svn.osgeo.org/postgis/trunk@7055 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/scripts/plpgsql/st_histogram.sql b/raster/scripts/plpgsql/st_histogram.sql new file mode 100644 index 000000000..ff760a081 --- /dev/null +++ b/raster/scripts/plpgsql/st_histogram.sql @@ -0,0 +1,348 @@ +---------------------------------------------------------------------- +-- _ST_Values(rast raster, band int) +-- Return all rast pixels values which center are in a geometry +-- Values are returned as groups of identical adjacent values (value, count) +-- in order to reduce the number of row returned. +---------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION _ST_Values(rast raster, band int, geom geometry, OUT val float8, OUT count int) + RETURNS SETOF record AS + $$ + DECLARE + geomintersect geometry; + m int[]; + x integer := 0; + y integer := 0; + curval float8; + BEGIN + m := ST_GeomExtent2RasterCoord(rast, geom); + -- Get the intersection between with the geometry. + geomintersect := ST_Intersection(geom, ST_ConvexHull(rast)); + + -- If the intersection is empty, return false + IF m[1] IS NULL AND m[2] IS NULL AND m[3] IS NULL AND m[4] IS NULL THEN + RETURN; + END IF; + + count := 0; + val := ST_Value(rast, band, m[1], m[2]); + FOR x IN m[1]..m[3] LOOP + FOR y IN m[2]..m[4] LOOP + -- Check first if the pixel intersects with the geometry. Many won't. + IF ST_Intersects(geom, ST_Centroid(ST_PixelAsPolygon(rast, x, y))) THEN + curval = ST_Value(rast, band, x, y); + IF NOT curval IS NULL THEN + IF curval = val THEN + count := count + 1; + ELSE + RETURN NEXT; + val := curval; + count := 1; + END IF; + END IF; + END IF; + END LOOP; + END LOOP; + RETURN NEXT; + RETURN; + END; + $$ + LANGUAGE 'plpgsql' IMMUTABLE STRICT; + +---------------------------------------------------------------------- +-- _ST_Values(rast raster, band int) +-- Return all rast pixels values +-- Values are returned as groups of identical adjacent values (value, count) +-- in order to reduce the number of row returned. +---------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION _ST_Values(rast raster, band int, OUT val float8, OUT count int) + RETURNS SETOF record AS + $$ + DECLARE + x int; + y int; + width int := ST_Width(rast); + height int := ST_Height(rast); + curval float8; + BEGIN + count := 0; + val := ST_Value(rast, band, 1, 1); + FOR x IN 1..width LOOP + FOR y IN 1..height LOOP + curval = ST_Value(rast, band, x, y); + IF NOT curval IS NULL THEN + IF curval = val THEN + count := count + 1; + ELSE + RETURN NEXT; + val := curval; + count := 1; + END IF; + END IF; + END LOOP; + END LOOP; + RETURN NEXT; + RETURN; + END; + $$ + LANGUAGE 'plpgsql'; + +---------------------------------------------------------------------- +-- ST_Histogram(rast raster, band int) group +-- Return a set of (val, count) rows forming the value histogram for a raster +---------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION ST_Histogram(rast raster, band int, OUT val double precision, OUT count bigint) +RETURNS SETOF record + AS $$ + SELECT (vc).val val, sum((vc).count)::bigint count + FROM (SELECT _ST_Values($1, $2) vc) foo GROUP BY (vc).val; + $$ + LANGUAGE SQL; + +CREATE OR REPLACE FUNCTION ST_Histogram(rast raster, OUT val double precision, OUT count bigint) +RETURNS SETOF record + AS $$ + SELECT (vc).val val, sum((vc).count)::bigint count + FROM (SELECT _ST_Values($1, 1) vc) foo GROUP BY (vc).val; + $$ + LANGUAGE SQL; + +---------------------------------------------------------------------- +-- ST_Histogram(rast raster, band int, geom geometry) group +-- Return a set of (val, count) rows forming the value histogram for the area of a raster covered by a polygon geometry. +-- Pixels are selected only when their center intersects the polygon +---------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION ST_Histogram(rast raster, band int, geom geometry, OUT val double precision, OUT count bigint) +RETURNS SETOF record + AS $$ + SELECT (vc).val val, sum((vc).count)::bigint count + FROM (SELECT _ST_Values($1, $2, $3) vc) foo GROUP BY (vc).val; + $$ + LANGUAGE SQL; + +CREATE OR REPLACE FUNCTION ST_Histogram(rast raster, geom geometry, OUT val double precision, OUT count bigint) +RETURNS SETOF record + AS $$ + SELECT (vc).val val, sum((vc).count)::bigint count + FROM (SELECT _ST_Values($1, 1, $2) vc) foo GROUP BY (vc).val; + $$ + LANGUAGE SQL; + +---------------------------------------------------------------------- +-- This variant might be faster (not using an intermediate _ST_Values function) +---------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION ST_Histogram2(rast raster, band int, OUT val double precision, OUT count bigint) +RETURNS SETOF record + AS $$ + SELECT val, count(*) count + FROM (SELECT ST_Value($1, $2, x, y) val FROM generate_series(1, ST_Width($1)) x , generate_series(1, ST_Height($1)) y) foo + GROUP BY val; + $$ + LANGUAGE SQL IMMUTABLE; + +SELECT (hist).val val, sum((hist).count) count +FROM (SELECT ST_Histogram2(rast, 1) hist FROM srtm_22_03_tiled_10x10) foo +GROUP BY val +ORDER BY count DESC + +---------------------------------------------------------------------- +-- Other variant (not grouping in the function) (not using an intermediate _ST_Values function) +---------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION ST_Histogram3(rast raster, band int, OUT val double precision) +RETURNS SETOF double precision + AS $$ + SELECT ST_Value($1, $2, x, y) val FROM generate_series(1, ST_Width($1)) x , generate_series(1, ST_Height($1)) y; + $$ + LANGUAGE SQL IMMUTABLE; + +SELECT val, count(*) count +FROM (SELECT ST_Histogram3(rast, 1) val FROM srtm_22_03_tiled_10x10) foo +GROUP BY val +ORDER BY count DESC + +SELECT ST_Histogram3(rast, 1) val, count(*) count FROM srtm_22_03_tiled_10x10 +GROUP BY val +ORDER BY count DESC + +---------------------------------------------------------------------- +-- This might actually be the fasters query to get the histogram +---------------------------------------------------------------------- +SELECT val, count(*) count +FROM (SELECT ST_Value(rast, 1, x, y) val + FROM generate_series(1, 10) x, generate_series(1, 10) y, srtm_22_03_tiled_10x10 + ) foo +GROUP BY val +ORDER BY count DESC + +---------------------------------------------------------------------- +-- Example 1: Query returning the histogram for a raster tile (one at a time) +---------------------------------------------------------------------- + +------------------------------- +-- The subquery gets the histogram for each tile +-- The main query split the resulting records in their two components (val & count) +------------------------------- +SELECT rid, + (hist).val, + (hist).count +FROM (SELECT rid, + ST_Histogram(rast) hist + FROM srtm_22_03_tiled_25x25 + WHERE rid = 234 + ) foo +ORDER BY (hist).count DESC + +---------------------------------------------------------------------- +-- Example 2: Query returning the complete histogram for a tiled raster coverage +- (might be very long) +---------------------------------------------------------------------- + +------------------------------- +-- The subquery gets the histogram for each tile +-- The main query split the resulting records in their two +-- components (val & count) and sum the count over all the tiles +------------------------------- +SELECT (hist).val, + SUM((hist).count) count +FROM (SELECT rid, + ST_Histogram(rast) hist + FROM srtm_22_03_tiled_25x25 + ) foo +GROUP BY (hist).val +ORDER BY count DESC + +---------------------------------------------------------------------- +-- Example 3: Query returning the mean pixel value for each tile of a +-- tiled raster (might be very long) +---------------------------------------------------------------------- + +------------------------------- +-- The subquery gets the histogram for each tile +-- The main query split the resulting records in their two +-- components (val & count) computing a mean value per tile at the same time +------------------------------- +SELECT rid, + geom, + round(((SUM((hist).val * (hist).count)) / SUM((hist).count))::numeric, 2) meanval +FROM (SELECT rid, + rast::geometry geom, + ST_Histogram(rast) hist + FROM srtm_22_03_tiled_25x25 + ) foo +GROUP BY rid, geom +ORDER BY rid; + +---------------------------------------------------------------------- +-- Example 4: Query returning the most frequent pixel value for each tile +-- of a tiled raster (might be very long) +-- This example requires an aggregate function tracking the value +-- associated with the maximum count +---------------------------------------------------------------------- +CREATE TYPE dblIntSet AS ( + maxval double precision, + val int +); + +CREATE OR REPLACE FUNCTION maxFromDblIntSet(dblIntSet, dblIntSet) RETURNS dblIntSet AS +$$ SELECT CASE WHEN $1.maxval>$2.maxval THEN $1 ELSE $2 END $$ +LANGUAGE sql; + +CREATE OR REPLACE FUNCTION valFromDblIntSet(dblIntSet) RETURNS int AS +$$ SELECT $1.val $$ +LANGUAGE sql; + +CREATE AGGREGATE maxFromDblIntSet ( + BASETYPE = dblIntSet, + SFUNC = maxFromDblIntSet, + STYPE = dblIntSet, + INITCOND = '(0,0)', + FINALFUNC = valFromDblIntSet +); + +------------------------------- +-- Actual query +-- The subquery gets the histogram for each tile +-- The main query split the resulting records in their two +-- components (val & count) and compute the maximum count and its associated value +------------------------------- +SELECT rid, + geom, + maxFromDblIntSet(ROW((hist).count, (hist).val::int)) mostfreqval, + MAX((hist).count) count +FROM (SELECT rid, + rast::geometry geom, + ST_Histogram(rast) hist + FROM srtm_22_03_tiled_25x25 + ) foo +GROUP BY rid, geom +ORDER BY rid + +---------------------------------------------------------------------- +-- Example 5: Query returning the most frequent pixel value per polygon from a raster +-- Do not use when the raster is big, in this case it should be tiled and +-- the next example (6) should be used instead +---------------------------------------------------------------------- +SELECT polyid, + geom, + maxFromDblIntSet(ROW((hist).count, (hist).val::int)) mostfreqval, + MAX((hist).count) count +FROM ( + SELECT polyid, + geom, + ST_Histogram(rast, geom) hist + FROM srtm_22_03, mypolygons + WHERE ST_Intersects(rast, geom) + ) foo +GROUP BY polyid, geom + +---------------------------------------------------------------------- +-- Example 6: Query returning the most frequent pixel value per polygon on a tiled raster coverage +---------------------------------------------------------------------- + +------------------------------- +-- The first subquery gets the histogram for each tile +-- The second subquery split the resulting records in their two +-- components (val & count) and sum the count for each polygon-value couple +-- The main query compute the maximum count and its associated value +------------------------------- +SELECT polyid, + geom, + maxFromDblIntSet(ROW(count, val)) mostfreqval, + MAX(count) count +FROM (SELECT polyid, + geom, (hist).val::int val, + SUM((hist).count) count + FROM (SELECT polyid, + geom, + ST_Histogram(rast, geom) hist + FROM srtm_22_03_tiled_25x25, mypolygons + WHERE ST_Intersects(rast, geom) + ) foo + GROUP BY polyid, geom, (hist).val + ) bar +GROUP BY polyid, geom + +---------------------------------------------------------------------- +-- Example 7: Query returning the mean pixel value per polygon on a tiled raster coverage +---------------------------------------------------------------------- + +------------------------------- +-- The first subquery gets the histogram for each tile +-- The second subquery split the resulting records in their two +-- components (val & count) and sum the count for each polygon-value couple +-- The main query compute the mean pixel value +------------------------------- +SELECT polyid, + geom, + round((SUM(val * count) / SUM(count))::numeric, 2) meanval +FROM (SELECT polyid, + geom, (hist).val::int val, + SUM((hist).count) count + FROM (SELECT polyid, + geom, + ST_Histogram(rast, geom) hist + FROM srtm_22_03_tiled_25x25, mypolygons + WHERE ST_Intersects(rast, geom) + ) foo + GROUP BY polyid, geom, (hist).val + ) bar +GROUP BY polyid, geom \ No newline at end of file