1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
|
---------------------------------------------------------------------
-- 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 gt.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(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
|
