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
|
----------------------------------------------------------------------
--
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- NOTE: The ST_Clip function found in this file still many improvements before being implemented in C.
CREATE OR REPLACE FUNCTION ST_Clip(rast raster, x int, y int, width int, height int)
RETURNS raster AS
$$
DECLARE
newrast raster := ST_MakeEmptyRaster(width, height, ST_UpperLeftX(rast), ST_UpperLeftY(rast),
ST_ScaleX(rast), ST_ScaleY(rast), ST_SkewX(rast), ST_SkewY(rast), ST_SRID(rast));
numband int := ST_Numbands(rast);
band int;
cx int;
cy int;
newwidth int := CASE WHEN x + width > ST_Width(rast) THEN ST_Width(rast) - x ELSE width END;
newheight int := CASE WHEN y + height > ST_Height(rast) THEN ST_Height(rast) - y ELSE height END;
BEGIN
FOR b IN 1..numband LOOP
newrast := ST_AddBand(newrast, ST_PixelType(rast, band), ST_NodataValue(rast, band), ST_NodataValue(rast, band));
FOR cx IN 1..newwidth LOOP
FOR cy IN 1..newheight LOOP
newrast := ST_SetValue(newrast, band, cx, cy, ST_Value(rast, band, cx + x - 1, cy + y - 1));
END LOOP;
END LOOP;
END LOOP;
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';
-------------------------------------------------------------------
-- ST_Clip
-- Clip the values of a raster to the shape of a polygon.
--
-- rast - raster to be clipped
-- band - limit the result to only one band
-- geom - geometry defining the she to clip the raster
-- nodata - define (if there is none defined) or replace the raster nodata value with this value
-- trimraster - limit the extent of the result to the extent of the geometry
-- Todo:
-- test point
-- test line
-- test polygon smaller than pixel
-- test and optimize raster totally included in polygon
CREATE OR REPLACE FUNCTION ST_Clip(rast raster, band int, geom geometry, nodata float8 DEFAULT null, trimraster boolean DEFAULT false)
RETURNS raster AS
$$
DECLARE
sourceraster raster := rast;
newrast raster;
geomrast raster;
numband int := ST_Numbands(rast);
bandstart int;
bandend int;
newextent text;
newnodata float8;
newpixtype text;
bandi int;
BEGIN
IF rast IS NULL THEN
RETURN null;
END IF;
IF geom IS NULL THEN
RETURN rast;
END IF;
IF band IS NULL THEN
bandstart := 1;
bandend := numband;
ELSEIF ST_HasNoBand(rast, band) THEN
RAISE NOTICE 'Raster do not have band %. Returning null', band;
RETURN null;
ELSE
bandstart := band;
bandend := band;
END IF;
newpixtype := ST_BandPixelType(rast, bandstart);
newnodata := coalesce(nodata, ST_BandNodataValue(rast, bandstart), ST_MinPossibleValue(newpixtype));
newextent := CASE WHEN trimraster THEN 'INTERSECTION' ELSE 'FIRST' END;
--RAISE NOTICE 'newextent=%', newextent;
-- Convert the geometry to a raster
geomrast := ST_AsRaster(geom, rast, ST_BandPixelType(rast, band), 1, newnodata);
-- Set the newnodata
sourceraster := ST_SetBandNodataValue(sourceraster, bandstart, newnodata);
-- Compute the first raster band
newrast := ST_MapAlgebraExpr(sourceraster, bandstart, geomrast, 1, '[rast1.val]', newpixtype, newextent);
FOR bandi IN bandstart+1..bandend LOOP
--RAISE NOTICE 'bandi=%', bandi;
-- for each band we must determine the nodata value
newpixtype := ST_BandPixelType(rast, bandi);
newnodata := coalesce(nodata, ST_BandNodataValue(sourceraster, bandi), ST_MinPossibleValue(newpixtype));
sourceraster := ST_SetBandNodataValue(sourceraster, bandi, newnodata);
newrast := ST_AddBand(newrast, ST_MapAlgebraExpr(sourceraster, bandi, geomrast, 1, '[rast1.val]', newpixtype, newextent));
END LOOP;
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';
-- Variant defaulting to band 1
CREATE OR REPLACE FUNCTION ST_Clip(rast raster, geom geometry, nodata float8 DEFAULT null, trimraster boolean DEFAULT false)
RETURNS raster
AS $$
SELECT ST_Clip($1, 1, $2, $3, $4);
$$ LANGUAGE 'SQL';
-- Variant defaulting nodata to the one of the raster or the min possible value
CREATE OR REPLACE FUNCTION ST_Clip(rast raster, band int, geom geometry, trimraster boolean)
RETURNS raster
AS $$
SELECT ST_Clip($1, $2, $3, null, $4);
$$ LANGUAGE 'SQL';
-- Variant defaulting nodata to the one of the raster or the min possible value
CREATE OR REPLACE FUNCTION ST_Clip(rast raster, geom geometry, trimraster boolean)
RETURNS raster
AS $$
SELECT ST_Clip($1, 1, $2, null, $3);
$$ LANGUAGE 'SQL';
-- 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, 0), '32BF', val, 0);
END;
$$
LANGUAGE 'plpgsql';
SELECT ST_Clip(ST_TestRaster(10, 10, 2), 1, ST_Buffer(ST_MakePoint(8, 5), 4)) rast
-- Test one band raster
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_Clip(ST_TestRaster(10, 10, 2), 1, ST_Buffer(ST_MakePoint(8, 5), 4))) gv
-- Test two bands raster
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_Clip(ST_AddBand(ST_TestRaster(10, 10, 2), '16BUI'::text, 4, 0), null, ST_Buffer(ST_MakePoint(8, 5), 4)), 2) gv
-- Test one band raster with trimraster set to true
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_Clip(ST_TestRaster(10, 10, 2), 1, ST_Buffer(ST_MakePoint(8, 5), 4), null, true)) gv
-- Test two bands raster with trimraster set to true
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_Clip(ST_AddBand(ST_TestRaster(10, 10, 2), '16BUI'::text, 4, 0), null, ST_Buffer(ST_MakePoint(8, 5), 4), null, true), 2) gv
-- Test nodatavalue set by the first raster
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_SetBandNoDataValue(ST_Clip(ST_TestRaster(10, 10, 2), 1, ST_Buffer(ST_MakePoint(8, 5), 4)), null)) gv
-- Test nodatavalue set by the parameter
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_SetBandNoDataValue(ST_Clip(ST_TestRaster(10, 10, 2), 1, ST_Buffer(ST_MakePoint(8, 5), 4), -10), null)) gv
-- Test nodatavalue set by the min possible value
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_SetBandNoDataValue(ST_Clip(ST_SetBandNoDataValue(ST_TestRaster(10, 10, 2), null), 1, ST_Buffer(ST_MakePoint(8, 5), 4)), null)) gv
-- Test the variant defaulting to band 1
SELECT ST_Numbands(ST_Clip(ST_AddBand(ST_TestRaster(10, 10, 2), '16BUI'::text, 4, 0), ST_Buffer(ST_MakePoint(8, 5), 4))) gv
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_Clip(ST_AddBand(ST_TestRaster(10, 10, 2), '16BUI'::text, 4, 0), ST_Buffer(ST_MakePoint(8, 5), 4)), 1) gv
-- Test defaulting to min possible value and band 1
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_SetBandNoDataValue(ST_Clip(ST_SetBandNoDataValue(ST_TestRaster(10, 10, 2), null), ST_Buffer(ST_MakePoint(8, 5), 4)), null)) gv
-- Test defaulting to nodatavalue set by the first raster and band 1
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_SetBandNoDataValue(ST_Clip(ST_TestRaster(10, 10, 2), ST_Buffer(ST_MakePoint(8, 5), 4)), null)) gv
-- Test when band number does not exist
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_Clip(ST_TestRaster(10, 10, 2), 2, ST_Buffer(ST_MakePoint(8, 5), 4))) gv
-- Test point -- bug. The produced raster does not have the same alignment
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_Clip(ST_TestRaster(10, 10, 2), ST_MakePoint(8.5, 5.5))) gv
|
