1---------------------------------------------------------------------- 2-- 3-- 4-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca> 5-- 6---------------------------------------------------------------------- 7 8-- Note: The functions provided in this script are in developement. Do not use. 9 10-- Note: this script is dependent on 11-- _MapAlgebraParts(r1x int, r1y int, r1w int, r1h int, r2x int, r2y int, r2w int, r2h int) 12-- ST_SameAlignment(rast1ulx float8, rast1uly float8, rast1scalex float8, rast1scaley float8, rast1skewx float8, rast1skewy float8, rast2ulx float8, rast2uly float8, rast2scalex float8, rast2scaley float8, rast2skewx float8, rast2skewy float8) 13-- ST_IsEmpty(raster) 14-- ST_HasNoBand(raster, int) 15-- ST_BandIsNoData(raster, int) 16-- to be found in the script/plpgsql folder 17 18-------------------------------------------------------------------- 19-- ST_MapAlgebra - (two rasters version) Return a raster which 20-- values are the result of an SQL expression involving 21-- pixel values from input rasters bands. 22-- Arguments 23-- rast1 raster - First raster referred by rast1 in the expression. 24-- band1 integer - Band number of the first raster. Default to 1. 25-- rast2 raster - Second raster referred by rast2 in the expression. 26-- band2 integer - Band number of the second raster. Default to 1. 27-- expression text - SQL expression. Ex.: "rast1 + 2 * rast2" 28-- pixeltype text - Pixeltype assigned to the resulting raster. Expression 29-- results are truncated to this type. Default to the 30-- pixeltype of the first raster. 31-- extentexpr text - Raster extent of the result. Can be: 32-- -FIRST: Same extent as the first raster. Default. 33-- -SECOND: Same extent as the second) raster. Default. 34-- -INTERSECTION: Intersection of extent of the two rasters. 35-- -UNION: Union oof extent of the two rasters. 36-- nodata1expr text - Expression used when only rast1 pixel value are nodata or absent, i.e. rast2 pixel value is with data. 37-- nodata2expr text - Expression used when only rast2 pixel value are nodata or absent, i.e. rast1 pixel value is with data. 38-- nodatanodataexpr text - Expression used when both pixel values are nodata values or absent. 39 40-- Further enhancements: 41-- -Move the expression parameter in seventh position just before other expression parameters. 42-- -Optimization for UNION & FIRST. We might not have to iterate over all the new raster. See st_mapalgebra_optimized.sql 43-- -Add the possibility to assign the x or y coordinate of the pixel to the pixel (do the same to the one raster version). 44-- -Resample the second raster when necessary (Require ST_Resample) 45-- -More test with rotated images 46-------------------------------------------------------------------- 47CREATE OR REPLACE FUNCTION ST_MapAlgebra2(rast1 raster, 48 band1 integer, 49 rast2 raster, 50 band2 integer, 51 expression text, 52 pixeltype text, 53 extentexpr text, 54 nodata1expr text, 55 nodata2expr text, 56 nodatanodataexpr text) 57 RETURNS raster AS 58 $$ 59 DECLARE 60 x integer; 61 y integer; 62 r1 float; 63 r2 float; 64 rast1ulx float8; 65 rast1uly float8; 66 rast1width int; 67 rast1height int; 68 rast1scalex float8; 69 rast1scaley float8; 70 rast1skewx float8; 71 rast1skewy float8; 72 rast1nodataval float8; 73 rast1srid int; 74 75 rast1offsetx int; 76 rast1offsety int; 77 78 rast2ulx float8; 79 rast2uly float8; 80 rast2width int; 81 rast2height int; 82 rast2scalex float8; 83 rast2scaley float8; 84 rast2skewx float8; 85 rast2skewy float8; 86 rast2nodataval float8; 87 rast2srid int; 88 89 r1x int; 90 r1y int; 91 r1w int; 92 r1h int; 93 r2x int; 94 r2y int; 95 r2w int; 96 r2h int; 97 98 newrx int; 99 newry int; 100 101 newrast raster; 102 tmprast raster; 103 newsrid int; 104 105 newscalex float8; 106 newscaley float8; 107 newskewx float8; 108 newskewy float8; 109 newnodatavalue float8; 110 newpixeltype text; 111 newulx float8; 112 newuly float8; 113 newwidth int; 114 newheight int; 115 newoffsetx1 int; 116 newoffsety1 int; 117 newoffsetx2 int; 118 newoffsety2 int; 119 120 newval float; 121 newexpr text; 122 upnodatanodataexpr text; 123 upnodata1expr text; 124 upnodata2expr text; 125 upexpression text; 126 nodatanodataval float; 127 skipcomputation int; 128 129 zones int[]; 130 z11x int; 131 z11y int; 132 z11w int; 133 z11h int; 134 z12x int; 135 z12y int; 136 z12w int; 137 z12h int; 138 z13x int; 139 z13y int; 140 z13w int; 141 z13h int; 142 z14x int; 143 z14y int; 144 z14w int; 145 z14h int; 146 z21x int; 147 z21y int; 148 z21w int; 149 z21h int; 150 z22x int; 151 z22y int; 152 z22w int; 153 z22h int; 154 z23x int; 155 z23y int; 156 z23w int; 157 z23h int; 158 z24x int; 159 z24y int; 160 z24w int; 161 z24h int; 162 zcx int; 163 zcy int; 164 zcw int; 165 zch int; 166 167 BEGIN 168 -- We have to deal with NULL, empty, hasnoband and hasnodatavalue band rasters... 169 -- These are respectively tested by "IS NULL", "ST_IsEmpty()", "ST_HasNoBand()" and "ST_BandIsNodata()" 170 171 -- If both raster are null, we return NULL. ok 172 -- If both raster do not have extent (are empty), we return an empty raster. ok 173 -- If both raster do not have the specified band, 174 -- we return a no band raster with the correct extent (according to the extent expression). ok 175 -- If both raster bands are nodatavalue and there is no replacement value, we return a nodata value band. ok 176 177 -- If only one raster is null or empty or has no band or hasnodata band we treat it as a nodata band raster. 178 -- If there is a replacement value we replace the missing raster values with this replacement value. ok 179 -- If there is no replacement value, we return a nodata value band. ok 180 181 -- What to do when only one raster is NULL or empty 182 -- If the extent expression is FIRST and the first raster is null we return NULL. ok 183 -- If the extent expression is FIRST and the first raster do not have extent (is empty), we return an empty raster. ok 184 -- If the extent expression is SECOND and the second raster is null we return NULL. ok 185 -- If the extent expression is SECOND and the second raster do not have extent (is empty), we return an empty raster. ok 186 -- If the extent expression is INTERSECTION and one raster is null or do not have extent (is empty), we return an empty raster. ok 187 -- If the extent expression is UNION and one raster is null or do not have extent (is empty), 188 -- we return a raster having the extent and the band characteristics of the other raster. ok 189 190 -- What to do when only one raster do not have the required band. 191 -- If the extent expression is FIRST and the first raster do not have the specified band, 192 -- we return a no band raster with the correct extent (according to the extent expression). ok 193 -- If the extent expression is SECOND and the second raster do not have the specified band, 194 -- we return a no band raster with the correct extent (according to the extent expression). ok 195 -- If the extent expression is INTERSECTION and one raster do not have the specified band, 196 -- we treat it as a nodata raster band. ok 197 -- If the extent expression is UNION and one raster do not have the specified band, 198 -- we treat it as a nodata raster band. ok 199 200 -- In all those cases, we make a warning. 201 202 -- Check if both rasters are NULL 203RAISE NOTICE 'ST_MapAlgebra2 000'; 204 205 IF rast1 IS NULL AND rast2 IS NULL THEN 206 RAISE NOTICE 'ST_MapAlgebra: Both raster are NULL. Returning NULL'; 207 RETURN NULL; 208 END IF; 209 210 -- Check if both rasters are empty (width or height = 0) 211 IF ST_IsEmpty(rast1) AND ST_IsEmpty(rast2) THEN 212 RAISE NOTICE 'ST_MapAlgebra: Both raster are empty. Returning an empty raster'; 213 RETURN ST_MakeEmptyRaster(0, 0, 0, 0, 0, 0, 0, 0, -1); 214 END IF; 215 216 rast1ulx := ST_UpperLeftX(rast1); 217 rast1uly := ST_UpperLeftY(rast1); 218 rast1width := ST_Width(rast1); 219 rast1height := ST_Height(rast1); 220 rast1scalex := ST_ScaleX(rast1); 221 rast1scaley := ST_ScaleY(rast1); 222 rast1skewx := ST_SkewX(rast1); 223 rast1skewy := ST_SkewY(rast1); 224 rast1srid := ST_SRID(rast1); 225 226 rast2ulx := ST_UpperLeftX(rast2); 227 rast2uly := ST_UpperLeftY(rast2); 228 rast2width := ST_Width(rast2); 229 rast2height := ST_Height(rast2); 230 rast2scalex := ST_ScaleX(rast2); 231 rast2scaley := ST_ScaleY(rast2); 232 rast2skewx := ST_SkewX(rast2); 233 rast2skewy := ST_SkewY(rast2); 234 rast2srid := ST_SRID(rast2); 235 236 -- Check if the first raster is NULL or empty 237 IF rast1 IS NULL OR ST_IsEmpty(rast1) THEN 238 rast1ulx := rast2ulx; 239 rast1uly := rast2uly; 240 rast1width := rast2width; 241 rast1height := rast2height; 242 rast1scalex := rast2scalex; 243 rast1scaley := rast2scaley; 244 rast1skewx := rast2skewx; 245 rast1skewy := rast2skewy; 246 rast1srid := rast2srid; 247 END IF; 248 -- Check if the second raster is NULL or empty 249 IF rast2 IS NULL OR ST_IsEmpty(rast2) THEN 250 rast2ulx := rast1ulx; 251 rast2uly := rast1uly; 252 rast2width := rast1width; 253 rast2height := rast1height; 254 rast2scalex := rast1scalex; 255 rast2scaley := rast1scaley; 256 rast2skewx := rast1skewx; 257 rast2skewy := rast1skewy; 258 rast2srid := rast1srid; 259 END IF; 260 261 -- Check for SRID 262 IF rast1srid != rast2srid THEN 263 RAISE EXCEPTION 'ST_MapAlgebra: Provided raster with different SRID. Aborting'; 264 END IF; 265 newsrid := rast1srid; 266 267 -- Check for alignment. (Rotation problem here) 268 IF NOT ST_SameAlignment(rast1ulx, rast1uly, rast1scalex, rast1scaley, rast1skewx, rast1skewy, rast2ulx, rast2uly, rast2scalex, rast2scaley, rast2skewx, rast2skewy) THEN 269 -- For now print an error message, but a more robust implementation should resample the second raster to the alignment of the first raster. 270 RAISE EXCEPTION 'ST_MapAlgebra: Provided raster do not have the same alignment. Aborting'; 271 END IF; 272 273 -- Set new pixel size and skew. We set it to the rast1 scale and skew 274 -- since both rasters are aligned and thus have the same scale and skew 275 newscalex := rast1scalex; 276 newscaley := rast1scaley; 277 newskewx := rast1skewx; 278 newskewy := rast1skewy; 279 280 --r1x & r2x are the offset of each rasters relatively to global extent 281 r1x := 0; 282 r2x := -st_world2rastercoordx(rast2, rast1ulx, rast1uly) + 1; 283 IF r2x < 0 THEN 284 r1x := -r2x; 285 r2x := 0; 286 END IF; 287 r1y := 0; 288 r2y := -st_world2rastercoordy(rast2, rast1ulx, rast1uly) + 1; 289 IF r2y < 0 THEN 290 r1y := -r2y; 291 r2y := 0; 292 END IF; 293 294 r1w := rast1width; 295 r1h := rast1height; 296 r2w := rast2width; 297 r2h := rast2height; 298 299 zones := _MapAlgebraParts(r1x + 1, r1y + 1, r1w, r1h, r2x + 1, r2y + 1, r2w, r2h); 300 z11x := zones[1]; 301 z11y := zones[2]; 302 z11w := zones[3]; 303 z11h := zones[4]; 304 z12x := zones[5]; 305 z12y := zones[6]; 306 z12w := zones[7]; 307 z12h := zones[8]; 308 z13x := zones[9]; 309 z13y := zones[10]; 310 z13w := zones[11]; 311 z13h := zones[12]; 312 z14x := zones[13]; 313 z14y := zones[14]; 314 z14w := zones[15]; 315 z14h := zones[16]; 316 z21x := zones[17]; 317 z21y := zones[18]; 318 z21w := zones[19]; 319 z21h := zones[20]; 320 z22x := zones[21]; 321 z22y := zones[22]; 322 z22w := zones[23]; 323 z22h := zones[24]; 324 z23x := zones[25]; 325 z23y := zones[26]; 326 z23w := zones[27]; 327 z23h := zones[28]; 328 z24x := zones[29]; 329 z24y := zones[30]; 330 z24w := zones[31]; 331 z24h := zones[32]; 332 zcx := zones[33]; 333 zcy := zones[34]; 334 zcw := zones[35]; 335 zch := zones[36]; 336 337 -- Compute x and y relative index of master and slave according to the extent expression (FIRST, SECOND, INTERSECTION or UNION) 338 IF extentexpr IS NULL OR upper(extentexpr) = 'FIRST' THEN 339 340 -- Check if rast1 is NULL 341 IF rast1 IS NULL THEN 342 RAISE NOTICE 'ST_MapAlgebra: FIRST raster is NULL. Returning NULL'; 343 RETURN NULL; 344 END IF; 345 346 -- Check if rast1 is empty 347 IF ST_IsEmpty(rast1) THEN 348 RAISE NOTICE 'ST_MapAlgebra: FIRST raster is empty. Returning an empty raster'; 349 RETURN ST_MakeEmptyRaster(0, 0, 0, 0, 0, 0, 0, 0, newsrid); 350 END IF; 351 352 -- Check if rast1 has the required band 353 IF ST_HasNoBand(rast1, band1) THEN 354 RAISE NOTICE 'ST_MapAlgebra: FIRST raster has no band. Returning a raster without band'; 355 RETURN ST_MakeEmptyRaster(rast1width, rast1height, rast1ulx, rast1uly, rast1scalex, rast1scaley, rast1skewx, rast1skewy, rast1srid); 356 END IF; 357 358 newulx := rast1ulx; 359 newuly := rast1uly; 360 newwidth := rast1width; 361 newheight := rast1height; 362 363 newrx := r1x; 364 newry := r1y; 365 z21w := 0; 366 z22w := 0; 367 z23w := 0; 368 z24w := 0; 369 370 ELSIF upper(extentexpr) = 'SECOND' THEN 371 372 -- Check if rast2 is NULL 373 IF rast2 IS NULL THEN 374 RAISE NOTICE 'ST_MapAlgebra: SECOND raster is NULL. Returning NULL'; 375 RETURN NULL; 376 END IF; 377 378 -- Check if rast2 is empty 379 IF ST_IsEmpty(rast2) THEN 380 RAISE NOTICE 'ST_MapAlgebra: SECOND raster is empty. Returning an empty raster'; 381 RETURN ST_MakeEmptyRaster(0, 0, 0, 0, 0, 0, 0, 0, newsrid); 382 END IF; 383 384 -- Check if rast2 has the required band 385 IF ST_HasNoBand(rast2, band2) THEN 386 RAISE NOTICE 'ST_MapAlgebra: SECOND raster has no band. Returning an empty raster'; 387 RETURN ST_MakeEmptyRaster(rast2width, rast2height, rast2ulx, rast2uly, rast2scalex, rast2scaley, rast2skewx, rast2skewy, rast2srid); 388 END IF; 389 390 newulx := rast2ulx; 391 newuly := rast2uly; 392 newwidth := rast2width; 393 newheight := rast2height; 394 395 newrx := r2x; 396 newry := r2y; 397 z11w := 0; 398 z12w := 0; 399 z13w := 0; 400 z14w := 0; 401 402 ELSIF upper(extentexpr) = 'INTERSECTION' THEN 403 404 -- Check if the intersection is empty. 405 IF zcw = 0 OR zch = 0 OR 406 rast1 IS NULL OR ST_IsEmpty(rast1) OR 407 rast2 IS NULL OR ST_IsEmpty(rast2) THEN 408 RAISE NOTICE 'ST_MapAlgebra: INTERSECTION of provided rasters is empty. Returning an empty raster'; 409 RETURN ST_MakeEmptyRaster(0, 0, 0, 0, 0, 0, 0, 0, newsrid); 410 END IF; 411 412 -- Compute the new ulx and uly 413 newulx := st_raster2worldcoordx(rast1, zcx - r1x + 1, zcy - r1y + 1); 414 newuly := st_raster2worldcoordy(rast1, zcx - r1x + 1, zcy - r1y + 1); 415 newwidth := zcw; 416 newheight := zch; 417 418 newrx := zcx; 419 newry := zcy; 420 z11w := 0; 421 z12w := 0; 422 z13w := 0; 423 z14w := 0; 424 z21w := 0; 425 z22w := 0; 426 z23w := 0; 427 z24w := 0; 428 429 ELSIF upper(extentexpr) = 'UNION' THEN 430 431 IF rast1 IS NULL OR ST_IsEmpty(rast1) THEN 432 newulx := rast2ulx; 433 newuly := rast2uly; 434 newwidth := rast2width; 435 newheight := rast2height; 436 437 newrx := r2x; 438 newry := r2y; 439 z21w := 0; 440 z22w := 0; 441 z23w := 0; 442 z24w := 0; 443 ELSIF rast2 IS NULL OR ST_IsEmpty(rast2) THEN 444 newulx := rast1ulx; 445 newuly := rast1uly; 446 newwidth := rast1width; 447 newheight := rast1height; 448 449 newrx := r1x; 450 newry := r1y; 451 z11w := 0; 452 z12w := 0; 453 z13w := 0; 454 z14w := 0; 455 ELSE 456 newrx := 0; 457 newry := 0; 458 459 newulx := st_raster2worldcoordx(rast1, r1x + 1, r1y + 1); 460 newuly := st_raster2worldcoordy(rast1, r1x + 1, r1y + 1); 461 newwidth := max(r1x + r1w, r2x + r2w); 462 newheight := max(r1y + r1h, r2y + r2h); 463 464 END IF; 465 ELSE 466 RAISE EXCEPTION 'ST_MapAlgebra: Unhandled extent expression "%". Only MASTER, INTERSECTION and UNION are accepted. Aborting.', upper(extentexpr); 467 END IF; 468 469 -- Check if both rasters do not have the specified band. 470 IF ST_HasNoband(rast1, band1) AND ST_HasNoband(rast2, band2) THEN 471 RAISE NOTICE 'ST_MapAlgebra: Both raster do not have the specified band. Returning a no band raster with the correct extent'; 472 RETURN ST_MakeEmptyRaster(newwidth, newheight, newulx, newuly, newscalex, newscaley, newskewx, newskewy, newsrid); 473 END IF; 474 475 -- Check newpixeltype 476 newpixeltype := pixeltype; 477 IF newpixeltype NOTNULL AND newpixeltype != '1BB' AND newpixeltype != '2BUI' AND newpixeltype != '4BUI' AND newpixeltype != '8BSI' AND newpixeltype != '8BUI' AND 478 newpixeltype != '16BSI' AND newpixeltype != '16BUI' AND newpixeltype != '32BSI' AND newpixeltype != '32BUI' AND newpixeltype != '32BF' AND newpixeltype != '64BF' THEN 479 RAISE EXCEPTION 'ST_MapAlgebra: Invalid pixeltype "%". Aborting.', newpixeltype; 480 END IF; 481 482 -- If no newpixeltype was provided, get it from the provided rasters. 483 IF newpixeltype IS NULL THEN 484 IF (upper(extentexpr) = 'SECOND' AND NOT ST_HasNoBand(rast2, band2)) OR ST_HasNoBand(rast1, band1) THEN 485 newpixeltype := ST_BandPixelType(rast2, band2); 486 ELSE 487 newpixeltype := ST_BandPixelType(rast1, band1); 488 END IF; 489 END IF; 490 491 -- Get the nodata value for first raster 492 IF NOT ST_HasNoBand(rast1, band1) AND ST_BandHasNodataValue(rast1, band1) THEN 493 rast1nodataval := ST_BandNodatavalue(rast1, band1); 494 ELSE 495 rast1nodataval := NULL; 496 END IF; 497 -- Get the nodata value for second raster 498 IF NOT ST_HasNoBand(rast2, band2) AND ST_BandHasNodatavalue(rast2, band2) THEN 499 rast2nodataval := ST_BandNodatavalue(rast2, band2); 500 ELSE 501 rast2nodataval := NULL; 502 END IF; 503 504 -- Determine new notadavalue 505 IF (upper(extentexpr) = 'SECOND' AND NOT rast2nodataval IS NULL) THEN 506 newnodatavalue := rast2nodataval; 507 ELSEIF NOT rast1nodataval IS NULL THEN 508 newnodatavalue := rast1nodataval; 509 ELSE 510 RAISE NOTICE 'ST_MapAlgebra: Both source rasters do not have a nodata value, nodata value for new raster set to the minimum value possible'; 511 newnodatavalue := ST_MinPossibleValue(newrast); 512 END IF; 513 514 upnodatanodataexpr := upper(nodatanodataexpr); 515 upnodata1expr := upper(nodata1expr); 516 upnodata2expr := upper(nodata2expr); 517 upexpression := upper(expression); 518 519 -- Initialise newrast with nodata-nodata values. Then we don't have anymore to set values for nodata-nodata pixels. 520 IF upnodatanodataexpr IS NULL THEN 521 nodatanodataval := newnodatavalue; 522 ELSE 523 EXECUTE 'SELECT ' || upnodatanodataexpr INTO nodatanodataval; 524 IF nodatanodataval IS NULL THEN 525 nodatanodataval := newnodatavalue; 526 ELSE 527 newnodatavalue := nodatanodataval; 528 END IF; 529 END IF; 530 531 ------------------------------------------------------------------- 532 --Create the raster receiving all the computed values. Initialize it to the new nodatavalue. 533 newrast := ST_AddBand(ST_MakeEmptyRaster(newwidth, newheight, newulx, newuly, newscalex, newscaley, newskewx, newskewy, newsrid), newpixeltype, newnodatavalue, newnodatavalue); 534 ------------------------------------------------------------------- 535 536RAISE NOTICE 'ST_MapAlgebra2 111 z11x=%, z11y=%, z11w=%, z11h=%', z11x, z11y, z11w, z11h; 537 -- First zone with only rast1 (z11) 538 IF z11w > 0 AND z11h > 0 AND NOT ST_BandIsNodata(rast1, band1) AND NOT nodata2expr IS NULL THEN 539 IF upnodata2expr = 'RAST' THEN 540 541 -- IF rast1nodataval != nodatanodataval THEN 542RAISE NOTICE 'ST_MapAlgebra2 222'; 543 -- newrast := ST_SetValues(newrast, 1, z11x, z11y, z11w, z11h, nodatanodataval); 544 -- END IF; 545RAISE NOTICE 'ST_MapAlgebra2 333 z11x=%, z11y=%, z11w=%, z11h=%', z11x, z11y, z11w, z11h; 546 newrast := ST_SetValues(newrast, 1, z11x, z11y, z11w, z11h, rast1, band1, NULL, TRUE); 547 548 ELSE 549RAISE NOTICE 'ST_MapAlgebra2 444'; 550 551 tmprast := ST_Clip(rast1, z11x - r1x + 1, z11y - r1y + 1, z11w, z11h); 552 newrast := ST_Mapalgebra2(newrast, 1, tmprast, 1, replace(upnodata2expr, 'RAST', 'RAST2'), NULL, 'FIRST', NULL, 'RAST', NULL); 553 554 END IF; 555 END IF; 556 557RAISE NOTICE 'ST_MapAlgebra2 555'; 558 559 -- Common zone (zc) 560 skipcomputation = 0; 561 IF zcw > 0 AND zch > 0 AND (NOT ST_BandIsNodata(rast1, band1) OR NOT ST_BandIsNodata(rast2, band2)) THEN 562 563RAISE NOTICE 'ST_MapAlgebra2 666'; 564 -- Initialize the zone with nodatavalue. We will not further compute nodata nodata pixels 565 -- newrast := ST_SetValues(newrast, 1, zcx + 1, zcy + 1, zcw, zch, newnodatavalue); 566 567 -- If rast1 is only nodata values, apply nodata1expr 568 IF ST_BandIsNodata(rast1, band1) THEN 569 570 IF nodata1expr IS NULL THEN 571 572 -- Do nothing 573 skipcomputation = 0; 574 575 ELSEIF upnodata1expr = 'RAST' THEN 576 577 -- Copy rast2 into newrast 578 newrast := ST_SetValues(newrast, 1, zcx, zcy, zcw, zch, , rast2, band2, NULL, 'KEEP'); 579 580 ELSE 581 -- If nodata1expr resume to a constant 582 IF position('RAST' in upnodata1expr) = 0 THEN 583 584 EXECUTE 'SELECT ' || upnodata1expr INTO newval; 585 IF newval IS NULL OR newval = newnodatavalue THEN 586 -- The constant is equal to nodata. We have nothing to compute since newrast was already initialized to nodata 587 skipcomputation := 2; 588 ELSEIF newnodatavalue IS NULL THEN 589 -- We can globally initialize to the constant only if there was no newnodatavalue. 590 newrast := ST_SetValues(newrast, 1, zcx, zcy, zcw, zch, newval); 591 skipcomputation := 2; 592 ELSE 593 -- We will assign the constant pixel by pixel. 594 skipcomputation := 1; 595 END IF; 596 END IF; 597 IF skipcomputation < 2 THEN 598 FOR x IN 1..zcw LOOP 599 FOR y IN 1..zch LOOP 600 r2 := ST_Value(rast2, band2, x + r2x, y + r2y); 601 IF (r2 IS NULL OR r2 = rast2nodataval) THEN 602 -- Do nothing since the raster have already been all set to this value 603 ELSE 604 IF skipcomputation < 1 THEN 605 newexpr := replace('SELECT ' || upnodata1expr, 'RAST', r2); 606 EXECUTE newexpr INTO newval; 607 IF newval IS NULL THEN 608 newval := newnodatavalue; 609 END IF; 610 END IF; 611 newrast = ST_SetValue(newrast, 1, x + zcx, y + zcy, newval); 612 END IF; 613 END LOOP; 614 END LOOP; 615 END IF; 616 END IF; 617 ELSEIF ST_BandIsNodata(rast2, band2) THEN 618 ELSE 619 END IF; 620 END IF; 621 622 RETURN newrast; 623 END; 624 $$ 625 LANGUAGE 'plpgsql'; 626 627CREATE OR REPLACE FUNCTION ST_TestRaster(ulx float8, uly float8, val float8) 628 RETURNS raster AS 629 $$ 630 DECLARE 631 BEGIN 632 RETURN ST_AddBand(ST_MakeEmptyRaster(5, 5, ulx, uly, 1, -1, 0, 0, -1), '32BF', val, -1); 633 END; 634 $$ 635 LANGUAGE 'plpgsql'; 636 637SELECT asbinary((gv).geom), (gv).val 638FROM st_pixelaspolygons(ST_TestRaster(-10, 2, 1)) gv; 639 640SELECT asbinary(_MapAlgebraAllPartsGeom(0, 0, 2, 1, 1, 0, 2, 1)) 641 642SELECT asbinary(pix.geom) as geom, pix.val 643FROM st_pixelaspolygons(ST_MapAlgebra2(ST_TestRaster(0, 1, 1), 1, ST_TestRaster(1, 0, 1), 1, '(rast1 + rast2) / 2', NULL, 'union', '2*rast', 'rast', NULL), 1) as pix 644 645