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