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