1DROP TABLE IF EXISTS raster_mapalgebra; 2CREATE TABLE raster_mapalgebra ( 3 rid integer, 4 rast raster 5); 6DROP TABLE IF EXISTS raster_mapalgebra_out; 7CREATE TABLE raster_mapalgebra_out ( 8 rid1 integer, 9 rid2 integer, 10 extent varchar, 11 rast raster 12); 13CREATE OR REPLACE FUNCTION make_test_raster( 14 rid integer, 15 width integer DEFAULT 2, 16 height integer DEFAULT 2, 17 ul_x double precision DEFAULT 0, 18 ul_y double precision DEFAULT 0, 19 skew_x double precision DEFAULT 0, 20 skew_y double precision DEFAULT 0, 21 initvalue double precision DEFAULT 1, 22 nodataval double precision DEFAULT 0 23) 24 RETURNS void 25 AS $$ 26 DECLARE 27 x int; 28 y int; 29 rast raster; 30 BEGIN 31 rast := ST_MakeEmptyRaster(width, height, ul_x, ul_y, 1, 1, skew_x, skew_y, 0); 32 rast := ST_AddBand(rast, 1, '8BUI', initvalue, nodataval); 33 34 INSERT INTO raster_mapalgebra VALUES (rid, rast); 35 36 RETURN; 37 END; 38 $$ LANGUAGE 'plpgsql'; 39-- no skew 40SELECT make_test_raster(0, 4, 4, -2, -2); 41SELECT make_test_raster(1, 2, 2, 0, 0, 0, 0, 2); 42SELECT make_test_raster(2, 2, 2, 1, -1, 0, 0, 3); 43SELECT make_test_raster(3, 2, 2, 1, 1, 0, 0, 4); 44SELECT make_test_raster(4, 2, 2, 2, 2, 0, 0, 5); 45 46-- skew 47SELECT make_test_raster(10, 4, 4, -2, -2, 1, -1); 48SELECT make_test_raster(11, 2, 2, 0, 0, 1, -1, 2); 49SELECT make_test_raster(12, 2, 2, 1, -1, 1, -1, 3); 50SELECT make_test_raster(13, 2, 2, 1, 1, 1, -1, 4); 51SELECT make_test_raster(14, 2, 2, 2, 2, 1, -1, 5); 52 53DROP FUNCTION IF EXISTS make_test_raster(integer, integer, integer, double precision, double precision, double precision, double precision, double precision, double precision); 54 55CREATE OR REPLACE FUNCTION raster_mapalgebra_intersection( 56 rast1 double precision, 57 rast2 double precision, 58 xy int[], 59 VARIADIC userargs text[] 60) 61 RETURNS double precision 62 AS $$ 63 DECLARE 64 BEGIN 65 IF rast1 IS NOT NULL AND rast2 IS NOT NULL THEN 66 RETURN rast1; 67 ELSE 68 RETURN NULL; 69 END IF; 70 71 RETURN NULL; 72 END; 73 $$ LANGUAGE 'plpgsql' IMMUTABLE; 74 75CREATE OR REPLACE FUNCTION raster_mapalgebra_union( 76 rast1 double precision, 77 rast2 double precision, 78 VARIADIC userargs text[] 79) 80 RETURNS double precision 81 AS $$ 82 DECLARE 83 BEGIN 84 CASE 85 WHEN rast1 IS NOT NULL AND rast2 IS NOT NULL THEN 86 RETURN ((rast1 + rast2)/2.); 87 WHEN rast1 IS NULL AND rast2 IS NULL THEN 88 RETURN NULL; 89 WHEN rast1 IS NULL THEN 90 RETURN rast2; 91 ELSE 92 RETURN rast1; 93 END CASE; 94 95 RETURN NULL; 96 END; 97 $$ LANGUAGE 'plpgsql' IMMUTABLE; 98 99CREATE OR REPLACE FUNCTION raster_mapalgebra_first( 100 rast1 double precision, 101 rast2 double precision, 102 VARIADIC userargs text[] 103) 104 RETURNS double precision 105 AS $$ 106 DECLARE 107 BEGIN 108 CASE 109 WHEN rast1 IS NOT NULL AND rast2 IS NOT NULL THEN 110 RETURN NULL; 111 WHEN rast1 IS NOT NULL THEN 112 RETURN rast1; 113 ELSE 114 RETURN NULL; 115 END CASE; 116 117 RETURN NULL; 118 END; 119 $$ LANGUAGE 'plpgsql' IMMUTABLE; 120 121CREATE OR REPLACE FUNCTION raster_mapalgebra_second( 122 rast1 double precision, 123 rast2 double precision, 124 VARIADIC userargs text[] 125) 126 RETURNS double precision 127 AS $$ 128 DECLARE 129 BEGIN 130 CASE 131 WHEN rast1 IS NOT NULL AND rast2 IS NOT NULL THEN 132 RETURN NULL; 133 WHEN rast2 IS NOT NULL THEN 134 RETURN rast2; 135 ELSE 136 RETURN NULL; 137 END CASE; 138 139 RETURN NULL; 140 END; 141 $$ LANGUAGE 'plpgsql'; 142 143-- INTERSECTION 144INSERT INTO raster_mapalgebra_out 145 (SELECT r1.rid, r2.rid, 'INTERSECTION', st_mapalgebrafct( 146 r1.rast, r2.rast, 'raster_mapalgebra_intersection(double precision, double precision, int[], text[])'::regprocedure, '32BF', 'INTERSECTION' 147 ) 148 FROM raster_mapalgebra r1 149 JOIN raster_mapalgebra r2 150 ON r1.rid != r2.rid 151 WHERE r1.rid = 0 152 AND r2.rid BETWEEN 1 AND 9 153 ) UNION ALL ( 154 SELECT r1.rid, r2.rid, 'INTERSECTION', st_mapalgebrafct( 155 r1.rast, r2.rast, 'raster_mapalgebra_intersection(double precision, double precision, int[], text[])'::regprocedure, '32BF', 'INTERSECTION' 156 ) 157 FROM raster_mapalgebra r1 158 JOIN raster_mapalgebra r2 159 ON r1.rid != r2.rid 160 WHERE r1.rid = 10 161 AND r2.rid BETWEEN 11 AND 19) 162; 163 164INSERT INTO raster_mapalgebra_out 165 SELECT NULL AS rid, rid, 'INTERSECTION', st_mapalgebrafct( 166 NULL::raster, rast, 'raster_mapalgebra_intersection(double precision, double precision, int[], text[])'::regprocedure, '32BF', 'INTERSECTION' 167 ) 168 FROM raster_mapalgebra 169; 170 171INSERT INTO raster_mapalgebra_out 172 SELECT rid, NULL AS rid, 'INTERSECTION', st_mapalgebrafct( 173 rast, NULL::raster, 'raster_mapalgebra_intersection(double precision, double precision, int[], text[])'::regprocedure, '32BF', 'INTERSECTION' 174 ) 175 FROM raster_mapalgebra 176; 177 178INSERT INTO raster_mapalgebra_out 179 SELECT NULL AS rid, NULL AS rid, 'INTERSECTION', st_mapalgebrafct( 180 NULL::raster, NULL::raster, 'raster_mapalgebra_intersection(double precision, double precision, int[], text[])'::regprocedure, '32BF', 'INTERSECTION' 181 ) 182; 183 184-- UNION 185INSERT INTO raster_mapalgebra_out 186 (SELECT r1.rid, r2.rid, 'UNION', st_mapalgebrafct( 187 r1.rast, r2.rast, 'raster_mapalgebra_union(double precision, double precision, text[])'::regprocedure, '32BF', 'UNION' 188 ) 189 FROM raster_mapalgebra r1 190 JOIN raster_mapalgebra r2 191 ON r1.rid != r2.rid 192 WHERE r1.rid = 0 193 AND r2.rid BETWEEN 1 AND 9 194 ) UNION ALL ( 195 SELECT r1.rid, r2.rid, 'UNION', st_mapalgebrafct( 196 r1.rast, r2.rast, 'raster_mapalgebra_union(double precision, double precision, text[])'::regprocedure, '32BF', 'UNION' 197 ) 198 FROM raster_mapalgebra r1 199 JOIN raster_mapalgebra r2 200 ON r1.rid != r2.rid 201 WHERE r1.rid = 10 202 AND r2.rid BETWEEN 11 AND 19) 203; 204 205INSERT INTO raster_mapalgebra_out 206 SELECT NULL AS rid, rid, 'UNION', st_mapalgebrafct( 207 NULL::raster, rast, 'raster_mapalgebra_union(double precision, double precision, text[])'::regprocedure, '32BF', 'UNION' 208 ) 209 FROM raster_mapalgebra 210; 211 212INSERT INTO raster_mapalgebra_out 213 SELECT rid, NULL AS rid, 'UNION', st_mapalgebrafct( 214 rast, NULL::raster, 'raster_mapalgebra_union(double precision, double precision, text[])'::regprocedure, '32BF', 'UNION' 215 ) 216 FROM raster_mapalgebra 217; 218 219INSERT INTO raster_mapalgebra_out 220 SELECT NULL AS rid, NULL AS rid, 'UNION', st_mapalgebrafct( 221 NULL::raster, NULL::raster, 'raster_mapalgebra_union(double precision, double precision, text[])'::regprocedure, '32BF', 'UNION' 222 ) 223; 224 225-- FIRST 226INSERT INTO raster_mapalgebra_out 227 (SELECT r1.rid, r2.rid, 'FIRST', st_mapalgebrafct( 228 r1.rast, r2.rast, 'raster_mapalgebra_first(double precision, double precision, text[])'::regprocedure, '32BF', 'FIRST' 229 ) 230 FROM raster_mapalgebra r1 231 JOIN raster_mapalgebra r2 232 ON r1.rid != r2.rid 233 WHERE r1.rid = 0 234 AND r2.rid BETWEEN 1 AND 9 235 ) UNION ALL ( 236 SELECT r1.rid, r2.rid, 'FIRST', st_mapalgebrafct( 237 r1.rast, r2.rast, 'raster_mapalgebra_first(double precision, double precision, text[])'::regprocedure, '32BF', 'FIRST' 238 ) 239 FROM raster_mapalgebra r1 240 JOIN raster_mapalgebra r2 241 ON r1.rid != r2.rid 242 WHERE r1.rid = 10 243 AND r2.rid BETWEEN 11 AND 19) 244; 245 246INSERT INTO raster_mapalgebra_out 247 SELECT NULL AS rid, rid, 'FIRST', st_mapalgebrafct( 248 NULL::raster, rast, 'raster_mapalgebra_first(double precision, double precision, text[])'::regprocedure, '32BF', 'FIRST' 249 ) 250 FROM raster_mapalgebra 251; 252 253INSERT INTO raster_mapalgebra_out 254 SELECT rid, NULL AS rid, 'FIRST', st_mapalgebrafct( 255 rast, NULL::raster, 'raster_mapalgebra_first(double precision, double precision, text[])'::regprocedure, '32BF', 'FIRST' 256 ) 257 FROM raster_mapalgebra 258; 259 260INSERT INTO raster_mapalgebra_out 261 SELECT NULL AS rid, NULL AS rid, 'FIRST', st_mapalgebrafct( 262 NULL::raster, NULL::raster, 'raster_mapalgebra_first(double precision, double precision, text[])'::regprocedure, '32BF', 'FIRST' 263 ) 264; 265 266-- SECOND 267INSERT INTO raster_mapalgebra_out 268 (SELECT r1.rid, r2.rid, 'SECOND', st_mapalgebrafct( 269 r1.rast, r2.rast, 'raster_mapalgebra_second(double precision, double precision, text[])'::regprocedure, '32BF', 'SECOND' 270 ) 271 FROM raster_mapalgebra r1 272 JOIN raster_mapalgebra r2 273 ON r1.rid != r2.rid 274 WHERE r1.rid = 0 275 AND r2.rid BETWEEN 1 AND 9 276 ) UNION ALL ( 277 SELECT r1.rid, r2.rid, 'SECOND', st_mapalgebrafct( 278 r1.rast, r2.rast, 'raster_mapalgebra_second(double precision, double precision, text[])'::regprocedure, '32BF', 'SECOND' 279 ) 280 FROM raster_mapalgebra r1 281 JOIN raster_mapalgebra r2 282 ON r1.rid != r2.rid 283 WHERE r1.rid = 10 284 AND r2.rid BETWEEN 11 AND 19) 285; 286 287INSERT INTO raster_mapalgebra_out 288 SELECT NULL AS rid, rid, 'SECOND', st_mapalgebrafct( 289 NULL::raster, rast, 'raster_mapalgebra_second(double precision, double precision, text[])'::regprocedure, '32BF', 'SECOND' 290 ) 291 FROM raster_mapalgebra 292; 293 294INSERT INTO raster_mapalgebra_out 295 SELECT rid, NULL AS rid, 'SECOND', st_mapalgebrafct( 296 rast, NULL::raster, 'raster_mapalgebra_second(double precision, double precision, text[])'::regprocedure, '32BF', 'SECOND' 297 ) 298 FROM raster_mapalgebra 299; 300 301INSERT INTO raster_mapalgebra_out 302 SELECT NULL AS rid, NULL AS rid, 'SECOND', st_mapalgebrafct( 303 NULL::raster, NULL::raster, 'raster_mapalgebra_second(double precision, double precision, text[])'::regprocedure, '32BF', 'SECOND' 304 ) 305; 306 307-- output 308SELECT 309 rid1, 310 rid2, 311 extent, 312 round(upperleftx::numeric, 3) AS upperleftx, 313 round(upperlefty::numeric, 3) AS upperlefty, 314 width, 315 height, 316 round(scalex::numeric, 3) AS scalex, 317 round(scaley::numeric, 3) AS scaley, 318 round(skewx::numeric, 3) AS skewx, 319 round(skewy::numeric, 3) AS skewy, 320 srid, 321 numbands, 322 pixeltype, 323 round(nodatavalue::numeric, 3) AS nodatavalue, 324 round(firstvalue::numeric, 3) AS firstvalue, 325 round(lastvalue::numeric, 3) AS lastvalue 326FROM ( 327 SELECT 328 rid1, 329 rid2, 330 extent, 331 mda.*, 332 bmd.*, 333 ST_Value(rast, 1, 1, 1) AS firstvalue, 334 ST_Value(rast, 1, ST_Width(rast), ST_Height(rast)) AS lastvalue 335 FROM raster_mapalgebra_out 336 LEFT JOIN LATERAL ST_Metadata(rast) AS mda ON true 337 LEFT JOIN LATERAL ST_BandMetadata(rast,1) AS bmd ON true 338) AS r; 339 340DROP FUNCTION IF EXISTS raster_mapalgebra_intersection(double precision, double precision, int[], VARIADIC text[]); 341DROP FUNCTION IF EXISTS raster_mapalgebra_union(double precision, double precision, VARIADIC text[]); 342DROP FUNCTION IF EXISTS raster_mapalgebra_first(double precision, double precision, VARIADIC text[]); 343DROP FUNCTION IF EXISTS raster_mapalgebra_second(double precision, double precision, VARIADIC text[]); 344 345DROP TABLE IF EXISTS raster_mapalgebra; 346DROP TABLE IF EXISTS raster_mapalgebra_out; 347