From fc659b26bb9921dcaceeeba7734c9dfc9a0ca44d Mon Sep 17 00:00:00 2001 From: Pierre Racine Date: Tue, 29 Nov 2011 16:15:25 +0000 Subject: [PATCH] Added some plpgsql functions git-svn-id: http://svn.osgeo.org/postgis/trunk@8253 b70326c6-7e19-0410-871a-916f4a2858ee --- .../plpgsql/st_areaweightedsummarystats.sql | 230 ++++++++++++++++++ .../scripts/plpgsql/st_createindexraster.sql | 73 ++++++ raster/scripts/plpgsql/st_splittable.sql | 47 ++++ raster/scripts/plpgsql/st_summarystatsagg.sql | 128 ++++++++++ 4 files changed, 478 insertions(+) create mode 100644 raster/scripts/plpgsql/st_areaweightedsummarystats.sql create mode 100644 raster/scripts/plpgsql/st_createindexraster.sql create mode 100644 raster/scripts/plpgsql/st_splittable.sql create mode 100644 raster/scripts/plpgsql/st_summarystatsagg.sql diff --git a/raster/scripts/plpgsql/st_areaweightedsummarystats.sql b/raster/scripts/plpgsql/st_areaweightedsummarystats.sql new file mode 100644 index 000000000..32e5e3ec2 --- /dev/null +++ b/raster/scripts/plpgsql/st_areaweightedsummarystats.sql @@ -0,0 +1,230 @@ +--------------------------------------------------------------------- +-- ST_AreaWeightedSummaryStats AGGREGATE +-- Compute statistics of a value weighted by the area of the corresponding geometry. +-- Specially written to be used with ST_Intersection(raster, geometry) +-- +-- Exemple +-- SELECT (aws).count, +-- (aws).distinctcount, +-- (aws).geom, +-- (aws).totalarea, +-- (aws).meanarea, +-- (aws).totalperimeter, +-- (aws).meanperimeter, +-- (aws).weightedsum, +-- (aws).weightedmean, +-- (aws).maxareavalue, +-- (aws).minareavalue, +-- (aws).maxcombinedareavalue, +-- (aws).mincombinedareavalue, +-- (aws).sum, +-- (aws).mean, +-- (aws).max, +-- (aws).min +-- FROM (SELECT ST_AreaWeightedSummaryStats(gv) aws +-- FROM (SELECT ST_Intersection(rt.rast, gt.geom) gv +-- FROM rasttable rt, geomtable gt +-- WHERE ST_Intersects(rt.rast, gt.geom) +-- ) foo +-- GROUP BY gt.id +-- ) foo2 +--------------------------------------------------------------------- + +--DROP TYPE arealweightedstats CASCADE; +CREATE TYPE arealweightedstats AS ( + count int, + distinctcount int, + geom geometry, + totalarea double precision, + meanarea double precision, + totalperimeter double precision, + meanperimeter double precision, + weightedsum double precision, + weightedmean double precision, + maxareavalue double precision, + minareavalue double precision, + maxcombinedareavalue double precision, + mincombinedareavalue double precision, + sum double precision, + mean double precision, + max double precision, + min double precision +); + +-- DROP TYPE arealweightedstatsstate CASCADE; +CREATE TYPE arealweightedstatsstate AS ( + count int, + distinctvalues double precision[], + unionedgeom geometry, + totalarea double precision, + totalperimeter double precision, + weightedsum double precision, + maxareavalue double precision[], + minareavalue double precision[], + combinedweightedareas double precision[], + sum double precision, + max double precision, + min double precision +); + +--------------------------------------------------------------------- +-- geomval_arealweightedstate +-- State function used by the ST_AreaWeightedSummaryStats aggregate +CREATE OR REPLACE FUNCTION geomval_arealweightedstate(aws arealweightedstatsstate, gv geomval) + RETURNS arealweightedstatsstate + AS $$ + DECLARE + i int; + ret arealweightedstatsstate; + newcombinedweightedareas double precision[] := ($1).combinedweightedareas; + newgeom geometry := ($2).geom; + geomtype text := GeometryType(($2).geom); + BEGIN + IF geomtype = 'GEOMETRYCOLLECTION' THEN + newgeom := ST_CollectionExtract(newgeom, 3); + END IF; + IF newgeom IS NULL OR ST_IsEmpty(newgeom) OR geomtype = 'POINT' OR geomtype = 'LINESTRING' OR geomtype = 'MULTIPOINT' OR geomtype = 'MULTILINESTRING' THEN + ret := aws; + ELSEIF $1 IS NULL THEN + ret := (1, + ARRAY[($2).val], + newgeom, + ST_Area(newgeom), + ST_Perimeter(newgeom), + ($2).val * ST_Area(newgeom), + ARRAY[ST_Area(newgeom), ($2).val], + ARRAY[ST_Area(newgeom), ($2).val], + ARRAY[ST_Area(newgeom)], + ($2).val, + ($2).val, + ($2).val + )::arealweightedstatsstate; + ELSE + -- Search for the new value in the array of distinct values + SELECT n FROM generate_series(1, array_length(($1).distinctvalues, 1)) n WHERE (($1).distinctvalues)[n] = ($2).val INTO i; +RAISE NOTICE 'i=% ',i; + -- If the value already exists, increment the corresponding area with the new area + IF NOT i IS NULL THEN + newcombinedweightedareas[i] := newcombinedweightedareas[i] + ST_Area(newgeom); + END IF; + ret := (($1).count + 1, + CASE WHEN i IS NULL THEN array_append(($1).distinctvalues, ($2).val) ELSE ($1).distinctvalues END, + ST_Union(($1).unionedgeom, newgeom), + ($1).totalarea + ST_Area(newgeom), + ($1).totalperimeter + ST_Perimeter(newgeom), + ($1).weightedsum + ($2).val * ST_Area(newgeom), + CASE WHEN ST_Area(newgeom) > (($1).maxareavalue)[1] THEN ARRAY[ST_Area(newgeom), ($2).val] ELSE ($1).maxareavalue END, + CASE WHEN ST_Area(newgeom) < (($1).minareavalue)[1] THEN ARRAY[ST_Area(newgeom), ($2).val] ELSE ($1).minareavalue END, + CASE WHEN i IS NULL THEN array_append(($1).combinedweightedareas, ST_Area(newgeom)) ELSE ($1).combinedweightedareas END, + ($1).sum + ($2).val, + greatest(($1).max, ($2).val), + least(($1).min, ($2).val) + )::arealweightedstatsstate; + END IF; + RETURN ret; + END; + $$ + LANGUAGE 'plpgsql'; + +CREATE OR REPLACE FUNCTION geomval_arealweightedstate(aws arealweightedstatsstate, geom geometry, val double precision) + RETURNS arealweightedstatsstate + AS $$ + SELECT geomval_arealweightedstate($1, ($2, $3)::geomval); + $$ LANGUAGE 'SQL'; + +--------------------------------------------------------------------- +-- geomval_arealweightedfinal +-- Final function used by the ST_AreaWeightedSummaryStats aggregate +CREATE OR REPLACE FUNCTION geomval_arealweightedfinal(aws arealweightedstatsstate) + RETURNS arealweightedstats + AS $$ + DECLARE + a RECORD; + maxarea double precision = 0.0; + minarea double precision = (($1).combinedweightedareas)[1]; + imax int := 1; + imin int := 1; + ret arealweightedstats; + BEGIN + -- Search for the max and the min areas in the array of all distinct values + FOR a IN SELECT n, (($1).combinedweightedareas)[n] warea FROM generate_series(1, array_length(($1).combinedweightedareas, 1)) n LOOP + IF a.warea > maxarea THEN + imax := a.n; + maxarea = a.warea; + END IF; + IF a.warea < minarea THEN + imin := a.n; + minarea = a.warea; + END IF; + END LOOP; + + ret := (($1).count, + array_length(($1).distinctvalues, 1), + ($1).unionedgeom, + ($1).totalarea, + ($1).totalarea / ($1).count, + ($1).totalperimeter, + ($1).totalperimeter / ($1).count, + ($1).weightedsum, + ($1).weightedsum / ($1).totalarea, + (($1).maxareavalue)[2], + (($1).minareavalue)[2], + (($1).distinctvalues)[imax], + (($1).distinctvalues)[imin], + ($1).sum, + ($1).sum / ($1).count, + ($1).max, + ($1).min + )::arealweightedstats; + RETURN ret; + END; + $$ + LANGUAGE 'plpgsql'; + +--------------------------------------------------------------------- +-- ST_AreaWeightedSummaryStats AGGREGATE +--------------------------------------------------------------------- +CREATE AGGREGATE ST_AreaWeightedSummaryStats(geomval) ( + SFUNC=geomval_arealweightedstate, + STYPE=arealweightedstatsstate, + FINALFUNC=geomval_arealweightedfinal +); + +--------------------------------------------------------------------- +-- ST_AreaWeightedSummaryStats AGGREGATE +--------------------------------------------------------------------- +CREATE AGGREGATE ST_AreaWeightedSummaryStats(geometry, double precision) ( + SFUNC=geomval_arealweightedstate, + STYPE=arealweightedstatsstate, + FINALFUNC=geomval_arealweightedfinal +); + + +SELECT id, + (aws).count, + (aws).distinctcount, + (aws).geom, + (aws).totalarea, + (aws).meanarea, + (aws).totalperimeter, + (aws).meanperimeter, + (aws).weightedsum, + (aws).weightedmean, + (aws).maxareavalue, + (aws).minareavalue, + (aws).maxcombinedareavalue, + (aws).mincombinedareavalue, + (aws).sum, + (aws).mean, + (aws).max, + (aws).min +FROM (SELECT ST_AreaWeightedSummaryStats((geom, weight)::geomval) as aws, id + FROM (SELECT ST_GeomFromEWKT('SRID=4269;POLYGON((0 0,0 10, 10 10, 10 0, 0 0))') as geom, 'a' as id, 100 as weight + UNION ALL + SELECT ST_GeomFromEWKT('SRID=4269;POLYGON((12 0,12 1, 13 1, 13 0, 12 0))') as geom, 'a' as id, 1 as weight + UNION ALL + SELECT ST_GeomFromEWKT('SRID=4269;POLYGON((10 0, 10 2, 12 2, 12 0, 10 0))') as geom, 'b' as id, 4 as weight + ) foo + GROUP BY id + ) foo2 + diff --git a/raster/scripts/plpgsql/st_createindexraster.sql b/raster/scripts/plpgsql/st_createindexraster.sql new file mode 100644 index 000000000..43ab1c6b7 --- /dev/null +++ b/raster/scripts/plpgsql/st_createindexraster.sql @@ -0,0 +1,73 @@ +---------------------------------------------------------------------------------------------------------------------- +-- Create an index raster. Georeference is borrowed from the provided raster. +-- pixeltype - The pixeltype of the resulting raster +-- startvalue - The first index value assigned to the raster. Default to 0. +-- incwithx - If true, the index increments when the x position of the pixel increments. +-- The index decrement on x otherwise. Default to true. +-- incwithy - If true, the index increments when the y position of the pixel increments. +-- The index decrement on y otherwise. Default to true. +-- columnfirst - If true, columns are traversed first. +-- Rows are traversed first otherwise. Default to true. +-- rowscanorder - If true, the raster is traversed in "row scan" order. +-- Row prime order (Boustrophedon) is used otherwise. Default to true. +-- falsecolinc - Overwrite the column index increment which is normally equal to ST_Height() +-- falserowinc - Overwrite the row index increment which is normally equal to ST_Width() +---------------------------------------------------------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION ST_CreateIndexRaster(rast raster, + pixeltype text DEFAULT '32BUI', + startvalue int DEFAULT 0, + incwithx boolean DEFAULT true, + incwithy boolean DEFAULT true, + columnfirst boolean DEFAULT true, + rowscanorder boolean DEFAULT true, + falsecolinc int DEFAULT NULL, + falserowinc int DEFAULT NULL) +RETURNS raster AS +$BODY$ +DECLARE + newraster raster := ST_AddBand(ST_MakeEmptyRaster(rast), pixeltype); + x int; + y int; + w int := ST_Width(newraster); + h int := ST_Height(newraster); + rowincx int := Coalesce(falserowinc, w); + colincx int := Coalesce(falsecolinc, h); + rowincy int := Coalesce(falserowinc, 1); + colincy int := Coalesce(falsecolinc, 1); + xdir int := CASE WHEN Coalesce(incwithx, true) THEN 1 ELSE w END; + ydir int := CASE WHEN Coalesce(incwithy, true) THEN 1 ELSE h END; + xdflag int := Coalesce(incwithx::int, 1); + ydflag int := Coalesce(incwithy::int, 1); + rsflag int := Coalesce(rowscanorder::int, 1); + newstartvalue int := Coalesce(startvalue, 0); + newcolumnfirst boolean := Coalesce(columnfirst, true); +BEGIN + IF newcolumnfirst THEN + IF colincx <= (h - 1) * rowincy THEN + RAISE EXCEPTION 'Column increment (now %) must be greater than the number of index on one column (now % pixel x % = %)...', colincx, h - 1, rowincy, (h - 1) * rowincy; + END IF; + FOR x IN 1..w LOOP + FOR y IN 1..h LOOP + newraster := ST_SetValue(newraster, x, y, abs(x - xdir) * colincx + abs(y - (h ^ ((abs(x - xdir + 1) % 2) | rsflag # ydflag))::int) * rowincy + newstartvalue); + END LOOP; + END LOOP; + ELSE + IF rowincx <= (w - 1) * colincy THEN + RAISE EXCEPTION 'Row increment (now %) must be greater than the number of index on one row (now % pixel x % = %)...', rowincx, w - 1, colincy, (w - 1) * colincy; + END IF; + FOR x IN 1..w LOOP + FOR y IN 1..h LOOP + newraster := ST_SetValue(newraster, x, y, abs(x - (w ^ ((abs(y - ydir + 1) % 2) | rsflag # xdflag))::int) * colincy + abs(y - ydir) * rowincx + newstartvalue); + END LOOP; + END LOOP; + END IF; + RETURN newraster; +END; +$BODY$ +LANGUAGE plpgsql VOLATILE; + +-- test +SELECT ST_AsBinary((gvxy).geom), (gvxy).val, (gvxy).x, (gvxy).y +FROM (SELECT ST_PixelAsPolygons(ST_CreateIndexRaster(ST_MakeEmptyRaster(2, 2, 0, 0, 1), '32BSI', null, null, null, null, null, 3, null)) gvxy) foo + +SELECT Coalesce(null::int, 2); \ No newline at end of file diff --git a/raster/scripts/plpgsql/st_splittable.sql b/raster/scripts/plpgsql/st_splittable.sql new file mode 100644 index 000000000..3c8cb7b35 --- /dev/null +++ b/raster/scripts/plpgsql/st_splittable.sql @@ -0,0 +1,47 @@ +---------------------------------------------------------------------------------------------------------------------- +-- SplitTable +-- Split a table into a series of table which names are composed of the concatenation of a prefix +-- and the value of a column. This function is usefull when loading many raster in one operation but +-- still wanting to split them in different tables afterward. They must have been loaded with the -F +-- raster2pgsql option so that different rasters are identifiable by a column. +-- +-- sourcetablename - The name of the table to split into multiple table +-- targettableschema - The schema in which to create the new set of table +-- targettableprefix - The prefix of the set of table names to create. +-- suffixcolumnname - The name of the column providing the suffix to each table name. +-- +-- Example to split the table 'test' into a set of table starting with 't_' and +-- ending with the value of the column 'rid' to be created in the 'public' schema. +-- +-- SELECT SplitTable('test', 'public', 't_', 'rid');; +---------------------------------------------------------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION SplitTable(sourcetablename text, targettableschema text, targettableprefix text, suffixcolumnname text) +RETURNS int AS +$BODY$ +DECLARE + newtablename text; + uniqueid RECORD; +BEGIN + FOR uniqueid IN EXECUTE 'SELECT DISTINCT ' || quote_ident(suffixcolumnname) || ' AS xyz123 FROM ' || sourcetablename LOOP + newtablename := targettableschema || '.' || targettableprefix || uniqueid.xyz123; + EXECUTE 'CREATE TABLE ' || quote_ident(newtablename) || ' AS SELECT * FROM ' || sourcetablename || ' WHERE ' || suffixcolumnname || ' = ' || uniqueid.xyz123; + END LOOP; + RETURN 1; +END; +$BODY$ +LANGUAGE plpgsql VOLATILE STRICT; + +--------------------------------------- +-- test +CREATE TABLE test AS +SELECT 1 AS rid, ST_MakeEmptyRaster(2,2,0,0,1,1,1,1,4326) AS rast +UNION ALL +SELECT 2 AS rid, ST_MakeEmptyRaster(2,2,0,0,1,1,1,1,4326) AS rast +UNION ALL +SELECT 1 AS rid, ST_MakeEmptyRaster(2,2,0,0,1,1,1,1,4326) AS rast +UNION ALL +SELECT 2 AS rid, ST_MakeEmptyRaster(2,2,0,0,1,1,1,1,4326) AS rast + +SELECT * FROM test; + +SELECT SplitTable('test', 'public', 't_', 'rid'); diff --git a/raster/scripts/plpgsql/st_summarystatsagg.sql b/raster/scripts/plpgsql/st_summarystatsagg.sql new file mode 100644 index 000000000..4e59cb9f9 --- /dev/null +++ b/raster/scripts/plpgsql/st_summarystatsagg.sql @@ -0,0 +1,128 @@ +--------------------------------------------------------------------- +-- ST_SummaryStatsAgg AGGREGATE +-- Compute summary statistics for an aggregation of raster. +-- +-- Exemple +-- SELECT (aws).count, +-- (aws).sum, +-- (aws).mean, +-- (aws).min, +-- (aws).max +-- FROM (SELECT ST_SummaryStatsAgg(gv) aws +-- FROM (SELECT ST_Clip(rt.rast, gt.geom) gv +-- FROM rasttable rt, geomtable gt +-- WHERE ST_Intersects(rt.rast, gt.geom) +-- ) foo +-- GROUP BY gt.id +-- ) foo2 +--------------------------------------------------------------------- + +-- DROP TYPE summarystatsstate CASCADE; +CREATE TYPE summarystatsstate AS ( + count int, + sum double precision, + min double precision, + max double precision +); + +--------------------------------------------------------------------- +-- raster_summarystatsstate +-- State function used by the ST_SummaryStatsAgg aggregate +CREATE OR REPLACE FUNCTION raster_summarystatsstate(sss summarystatsstate, rast raster, nband int DEFAULT 1, exclude_nodata_value boolean DEFAULT TRUE, sample_percent double precision DEFAULT 1) + RETURNS summarystatsstate + AS $$ + DECLARE + newstats summarystats; + ret summarystatsstate; + BEGIN + IF rast IS NULL THEN + RETURN sss; + END IF; + newstats := _ST_SummaryStats(rast, nband, exclude_nodata_value, sample_percent); + IF $1 IS NULL THEN + ret := (newstats.count, + newstats.sum, + newstats.min, + newstats.max)::summarystatsstate; + ELSE + ret := (sss.count + newstats.count, + sss.sum + newstats.sum, + least(sss.min, newstats.min), + greatest(sss.max, newstats.max))::summarystatsstate; + END IF; +RAISE NOTICE 'min=% ',ret.min; + RETURN ret; + END; + $$ + LANGUAGE 'plpgsql'; + +CREATE OR REPLACE FUNCTION raster_summarystatsstate(sss summarystatsstate, rast raster) + RETURNS summarystatsstate + AS $$ + SELECT raster_summarystatsstate($1, $2, 1, true, 1); + $$ LANGUAGE 'SQL'; + +--------------------------------------------------------------------- +-- raster_summarystatsfinal +-- Final function used by the ST_SummaryStatsAgg aggregate +CREATE OR REPLACE FUNCTION raster_summarystatsfinal(sss summarystatsstate) + RETURNS summarystats + AS $$ + DECLARE + ret summarystats; + BEGIN + ret := (($1).count, + ($1).sum, + ($1).sum / ($1).count, + null, + ($1).min, + ($1).max + )::summarystats; + RETURN ret; + END; + $$ + LANGUAGE 'plpgsql'; + +--------------------------------------------------------------------- +-- ST_SummaryStatsAgg AGGREGATE +--------------------------------------------------------------------- +CREATE AGGREGATE ST_SummaryStatsAgg(raster, int, boolean, double precision) ( + SFUNC=raster_summarystatsstate, + STYPE=summarystatsstate, + FINALFUNC=raster_summarystatsfinal +); + +CREATE AGGREGATE ST_SummaryStatsAgg(raster) ( + SFUNC=raster_summarystatsstate, + STYPE=summarystatsstate, + FINALFUNC=raster_summarystatsfinal +); + +-- Test +CREATE OR REPLACE FUNCTION ST_TestRaster(h integer, w integer, val float8) + RETURNS raster AS + $$ + DECLARE + BEGIN + RETURN ST_AddBand(ST_MakeEmptyRaster(h, w, 0, 0, 1, 1, 0, 0, -1), '32BF', val, -1); + END; + $$ + LANGUAGE 'plpgsql'; + +SELECT id, + (sss).count, + (sss).sum, + (sss).mean, + (sss).stddev, + (sss).min, + (sss).max +FROM (SELECT ST_SummaryStatsAgg(rast) as sss, id + FROM (SELECT 1 id, ST_TestRaster(2, 2, 2) rast + UNION ALL + SELECT 1 id, ST_TestRaster(2, 2, 4) rast + UNION ALL + SELECT 2 id, ST_TestRaster(2, 2, 4) rast + ) foo + GROUP BY id + ) foo2 + -- 2.40.0