self-hosted/postgis
self-hosted/postgis
1
0
Fork 0
mirror of https://git.osgeo.org/gitea/postgis/postgis synced 2024-10-24 09:02:37 +00:00
Code Issues Packages Projects Releases Wiki Activity

Added some plpgsql functions

Browse source 
git-svn-id: http://svn.osgeo.org/postgis/trunk@8253 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Pierre Racine 2011-11-29 16:15:25 +00:00
parent 576b575b13
commit fc659b26bb
4 changed files with 478 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

230
raster/scripts/plpgsql/st_areaweightedsummarystats.sql Normal file
View file

@ -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

73
raster/scripts/plpgsql/st_createindexraster.sql Normal file
View file

@ -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);

47
raster/scripts/plpgsql/st_splittable.sql Normal file
View file

@ -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');

128
raster/scripts/plpgsql/st_summarystatsagg.sql Normal file
View file

@ -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