1SET client_min_messages TO warning;
2
3DROP TABLE IF EXISTS raster_clip;
4CREATE TABLE raster_clip (
5	rid integer,
6	rast raster
7);
8DROP TABLE IF EXISTS geom_clip;
9CREATE TABLE geom_clip (
10	gid integer,
11	geom geometry
12);
13DROP TABLE IF EXISTS raster_clip_out;
14CREATE TABLE raster_clip_out (
15	tid integer,
16	rid integer,
17	gid integer,
18	rast raster
19);
20
21CREATE OR REPLACE FUNCTION make_test_raster(
22	rid integer,
23	width integer DEFAULT 2,
24	height integer DEFAULT 2,
25	ul_x double precision DEFAULT 0,
26	ul_y double precision DEFAULT 0,
27	skew_x double precision DEFAULT 0,
28	skew_y double precision DEFAULT 0,
29	initvalue double precision DEFAULT 1,
30	nodataval double precision DEFAULT 0
31)
32	RETURNS void
33	AS $$
34	DECLARE
35		x int;
36		y int;
37		rast raster;
38	BEGIN
39		rast := ST_MakeEmptyRaster(width, height, ul_x, ul_y, 1, -1, skew_x, skew_y, 0);
40		rast := ST_AddBand(rast, 1, '8BUI', initvalue, nodataval);
41
42		INSERT INTO raster_clip VALUES (rid, rast);
43
44		RETURN;
45	END;
46	$$ LANGUAGE 'plpgsql';
47
48-- Define three rasters
49-- The first one 1 band with a novalue defined and one pixel set to nodata value
50SELECT make_test_raster(1, 4, 4, 0, 0, 0, 0, 1, 0);
51UPDATE raster_clip SET rast = ST_SetValue(rast, 3, 2, NULL) WHERE rid = 1;
52
53-- 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
54SELECT make_test_raster(2, 4, 4, 0, 0, 0, 0, 10, 0);
55UPDATE raster_clip SET rast = ST_SetValue(rast, 3, 2, NULL) WHERE rid = 2;
56UPDATE raster_clip SET rast = ST_AddBand(rast, '8BUI'::text, 2, 0) WHERE rid = 2;
57UPDATE raster_clip SET rast = ST_SetValue(rast, 2, 3, 2, NULL) WHERE rid = 2;
58UPDATE raster_clip SET rast = ST_AddBand(rast, '8BUI'::text, 3, NULL) WHERE rid = 2;
59
60-- The third one 1 band skewed 40 degree, (Can't test this yet as ST_AsRaster() still produces badly aligned raster. See ticket #1574)
61--SELECT make_test_raster(3, 4, 4, 0, 0, 0, 0, 1, 0);
62--UPDATE raster_clip SET rast = ST_SetSkew(rast, -0.15, -0.15) WHERE rid = 3;
63--UPDATE raster_clip SET rast = ST_SetValue(rast, 3, 2, NULL) WHERE rid = 3;
64
65-- Add a first polygon small and outside the extent of the raster
66INSERT INTO geom_clip VALUES (1, ST_Buffer(ST_SetSRID(ST_MakePoint(-1, 1), 0), 0.2));
67-- Add a second polygon small, inside the extent of the raster but in the nodata value pixel
68INSERT INTO geom_clip VALUES (2, ST_Buffer(ST_SetSRID(ST_MakePoint(2.5, -1.5), 0), 0.2));
69-- Add a second polygon small but inside the extent of the raster
70INSERT INTO geom_clip VALUES (3, ST_Buffer(ST_SetSRID(ST_MakePoint(1.5, -1.5), 0), 0.2));
71-- Add a third polygon big cutting the raster
72INSERT INTO geom_clip VALUES (4, ST_Buffer(ST_SetSRID(ST_MakePoint(4, -2.5), 0), 2.8));
73-- Add a fourth polygon englobing the two first rasters
74INSERT INTO geom_clip VALUES (5, ST_Buffer(ST_SetSRID(ST_MakePoint(2, -2), 0), 3));
75
76DROP FUNCTION make_test_raster(integer, integer, integer, double precision, double precision, double precision, double precision, double precision, double precision);
77
78-- Test 1 without trimming, without defining a nodata value
79INSERT INTO raster_clip_out
80SELECT 1, rid, gid, ST_Clip(rast, geom, false)
81FROM raster_clip, geom_clip;
82
83-- Test 2 with trimming, without defining a nodata value
84INSERT INTO raster_clip_out
85SELECT 2, rid, gid, ST_Clip(rast, geom, true)
86FROM raster_clip, geom_clip;
87
88-- Test 3 without trimming, defining a nodata value
89INSERT INTO raster_clip_out
90SELECT 3, rid, gid, ST_Clip(rast, geom, ARRAY[255, 254, 253], false)
91FROM raster_clip, geom_clip;
92
93-- Test 4 with trimming, defining a nodata value
94INSERT INTO raster_clip_out
95SELECT 4, rid, gid, ST_Clip(rast, geom, ARRAY[255, 254, 253], true)
96FROM raster_clip, geom_clip;
97
98-- Display the metadata of the resulting rasters
99SELECT
100	tid,
101	rid,
102	gid,
103	round(upperleftx::numeric, 3) AS upperleftx,
104	round(upperlefty::numeric, 3) AS upperlefty,
105	width,
106	height,
107	round(scalex::numeric, 3) AS scalex,
108	round(scaley::numeric, 3) AS scaley,
109	round(skewx::numeric, 3) AS skewx,
110	round(skewy::numeric, 3) AS skewy,
111	srid,
112	numbands,
113	pixeltype,
114	round(nodatavalue::numeric, 3) AS nodatavalue
115FROM (
116	SELECT  tid,
117	        rid,
118		gid,
119		(ST_Metadata(rast)).*,
120		(ST_BandMetadata(rast, 1)).*
121	FROM raster_clip_out
122) AS r;
123
124-- Display the pixels and the values of the resulting rasters (raster 1)
125SELECT
126	tid,
127	rid,
128	gid,
129	(gvxy).x,
130	(gvxy).y,
131	(gvxy).val,
132	ST_AsText((gvxy).geom) geom
133FROM (SELECT tid, rid, gid, ST_PixelAsPolygons(rast) gvxy
134      FROM raster_clip_out
135      WHERE rid = 1
136) foo
137ORDER BY 1, 2, 3, 4, 5, 7;
138
139-- Display the pixels and the values of the resulting rasters (raster 2, 3 bands)
140SELECT
141	tid,
142	rid,
143	gid,
144	band,
145	(gvxy).x,
146	(gvxy).y,
147	(gvxy).val,
148	ST_AsText((gvxy).geom) geom
149FROM (SELECT tid, rid, gid, band, ST_PixelAsPolygons(rast, band) gvxy
150      FROM raster_clip_out, generate_series(1, 3) band
151      WHERE rid = 2
152) foo
153ORDER BY 1, 2, 3, 4, 5, 6, 8;
154
155DROP TABLE IF EXISTS geom_clip;
156DROP TABLE IF EXISTS raster_clip;
157DROP TABLE IF EXISTS raster_clip_out;
158