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
|
SET client_min_messages TO warning;
DROP TABLE IF EXISTS raster_clip;
CREATE TABLE raster_clip (
rid integer,
rast raster
);
DROP TABLE IF EXISTS geom_clip;
CREATE TABLE geom_clip (
gid integer,
geom geometry
);
DROP TABLE IF EXISTS raster_clip_out;
CREATE TABLE raster_clip_out (
tid integer,
rid integer,
gid integer,
rast raster
);
CREATE OR REPLACE FUNCTION make_test_raster(
rid integer,
width integer DEFAULT 2,
height integer DEFAULT 2,
ul_x double precision DEFAULT 0,
ul_y double precision DEFAULT 0,
skew_x double precision DEFAULT 0,
skew_y double precision DEFAULT 0,
initvalue double precision DEFAULT 1,
nodataval double precision DEFAULT 0
)
RETURNS void
AS $$
DECLARE
x int;
y int;
rast raster;
BEGIN
rast := ST_MakeEmptyRaster(width, height, ul_x, ul_y, 1, -1, skew_x, skew_y, 0);
rast := ST_AddBand(rast, 1, '8BUI', initvalue, nodataval);
INSERT INTO raster_clip VALUES (rid, rast);
RETURN;
END;
$$ LANGUAGE 'plpgsql';
-- Define three rasters
-- The first one 1 band with a novalue defined and one pixel set to nodata value
SELECT make_test_raster(1, 4, 4, 0, 0, 0, 0, 1, 0);
UPDATE raster_clip SET rast = ST_SetValue(rast, 3, 2, NULL) WHERE rid = 1;
-- The second one 3 bands with a novalue defined for the first two band but not set for the third band and one pixel set to nodata value in every band
SELECT make_test_raster(2, 4, 4, 0, 0, 0, 0, 10, 0);
UPDATE raster_clip SET rast = ST_SetValue(rast, 3, 2, NULL) WHERE rid = 2;
UPDATE raster_clip SET rast = ST_AddBand(rast, '8BUI'::text, 2, 0) WHERE rid = 2;
UPDATE raster_clip SET rast = ST_SetValue(rast, 2, 3, 2, NULL) WHERE rid = 2;
UPDATE raster_clip SET rast = ST_AddBand(rast, '8BUI'::text, 3, NULL) WHERE rid = 2;
-- The third one 1 band skewed 40 degree, (Can't test this yet as ST_AsRaster() still produces badly aligned raster. See ticket #1574)
--SELECT make_test_raster(3, 4, 4, 0, 0, 0, 0, 1, 0);
--UPDATE raster_clip SET rast = ST_SetSkew(rast, -0.15, -0.15) WHERE rid = 3;
--UPDATE raster_clip SET rast = ST_SetValue(rast, 3, 2, NULL) WHERE rid = 3;
-- Add a first polygon small and outside the extent of the raster
INSERT INTO geom_clip VALUES (1, ST_Buffer(ST_SetSRID(ST_MakePoint(-1, 1), 0), 0.2));
-- Add a second polygon small, inside the extent of the raster but in the nodata value pixel
INSERT INTO geom_clip VALUES (2, ST_Buffer(ST_SetSRID(ST_MakePoint(2.5, -1.5), 0), 0.2));
-- Add a second polygon small but inside the extent of the raster
INSERT INTO geom_clip VALUES (3, ST_Buffer(ST_SetSRID(ST_MakePoint(1.5, -1.5), 0), 0.2));
-- Add a third polygon big cutting the raster
INSERT INTO geom_clip VALUES (4, ST_Buffer(ST_SetSRID(ST_MakePoint(4, -2.5), 0), 2.8));
-- Add a fourth polygon englobing the two first rasters
INSERT INTO geom_clip VALUES (5, ST_Buffer(ST_SetSRID(ST_MakePoint(2, -2), 0), 3));
DROP FUNCTION make_test_raster(integer, integer, integer, double precision, double precision, double precision, double precision, double precision, double precision);
-- Test 1 without trimming, without defining a nodata value
INSERT INTO raster_clip_out
SELECT 1, rid, gid, ST_Clip(rast, geom, false)
FROM raster_clip, geom_clip;
-- Test 2 with trimming, without defining a nodata value
INSERT INTO raster_clip_out
SELECT 2, rid, gid, ST_Clip(rast, geom, true)
FROM raster_clip, geom_clip;
-- Test 3 without trimming, defining a nodata value
INSERT INTO raster_clip_out
SELECT 3, rid, gid, ST_Clip(rast, geom, ARRAY[255, 254, 253], false)
FROM raster_clip, geom_clip;
-- Test 4 with trimming, defining a nodata value
INSERT INTO raster_clip_out
SELECT 4, rid, gid, ST_Clip(rast, geom, ARRAY[255, 254, 253], true)
FROM raster_clip, geom_clip;
-- Display the metadata of the resulting rasters
SELECT
tid,
rid,
gid,
round(upperleftx::numeric, 3) AS upperleftx,
round(upperlefty::numeric, 3) AS upperlefty,
width,
height,
round(scalex::numeric, 3) AS scalex,
round(scaley::numeric, 3) AS scaley,
round(skewx::numeric, 3) AS skewx,
round(skewy::numeric, 3) AS skewy,
srid,
numbands,
pixeltype,
round(nodatavalue::numeric, 3) AS nodatavalue
FROM (
SELECT tid,
rid,
gid,
(ST_Metadata(rast)).*,
(ST_BandMetadata(rast, 1)).*
FROM raster_clip_out
) AS r;
-- Display the pixels and the values of the resulting rasters (raster 1)
SELECT
tid,
rid,
gid,
(gvxy).x,
(gvxy).y,
(gvxy).val,
ST_AsText((gvxy).geom) geom
FROM (SELECT tid, rid, gid, ST_PixelAsPolygons(rast) gvxy
FROM raster_clip_out
WHERE rid = 1
) foo
ORDER BY 1, 2, 3, 4, 5, 7;
-- Display the pixels and the values of the resulting rasters (raster 2, 3 bands)
SELECT
tid,
rid,
gid,
band,
(gvxy).x,
(gvxy).y,
(gvxy).val,
ST_AsText((gvxy).geom) geom
FROM (SELECT tid, rid, gid, band, ST_PixelAsPolygons(rast, band) gvxy
FROM raster_clip_out, generate_series(1, 3) band
WHERE rid = 2
) foo
ORDER BY 1, 2, 3, 4, 5, 6, 8;
DROP TABLE IF EXISTS geom_clip;
DROP TABLE IF EXISTS raster_clip;
DROP TABLE IF EXISTS raster_clip_out;
|
