1 /*
2 *
3 * WKTRaster - Raster Types for PostGIS
4 * http://trac.osgeo.org/postgis/wiki/WKTRaster
5 *
6 * Copyright (C) 2011-2013 Regents of the University of California
7 * <bkpark@ucdavis.edu>
8 * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
9 * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
10 * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
11 * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
12 * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
13 * Copyright (C) 2013 Nathaniel Hunter Clay <clay.nathaniel@gmail.com>
14 *
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * This program is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with this program; if not, write to the Free Software Foundation,
27 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28 *
29 */
30
31 #include <assert.h>
32
33 #include <postgres.h> /* for palloc */
34 #include <fmgr.h>
35 #include <funcapi.h>
36 #include <executor/spi.h>
37 #include <utils/lsyscache.h> /* for get_typlenbyvalalign */
38 #include <utils/array.h> /* for ArrayType */
39 #include <utils/builtins.h> /* for cstring_to_text */
40 #include <catalog/pg_type.h> /* for INT2OID, INT4OID, FLOAT4OID, FLOAT8OID and TEXTOID */
41 #include <executor/executor.h> /* for GetAttributeByName */
42
43 #include "../../postgis_config.h"
44 #include "lwgeom_pg.h"
45
46 #include "rtpostgis.h"
47 #include "rtpg_internal.h"
48
49 /* n-raster MapAlgebra */
50 Datum RASTER_nMapAlgebra(PG_FUNCTION_ARGS);
51 Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS);
52
53 /* raster union aggregate */
54 Datum RASTER_union_transfn(PG_FUNCTION_ARGS);
55 Datum RASTER_union_finalfn(PG_FUNCTION_ARGS);
56
57 /* raster clip */
58 Datum RASTER_clip(PG_FUNCTION_ARGS);
59
60 /* reclassify specified bands of a raster */
61 Datum RASTER_reclass(PG_FUNCTION_ARGS);
62
63 /* apply colormap to specified band of a raster */
64 Datum RASTER_colorMap(PG_FUNCTION_ARGS);
65
66 /* one-raster MapAlgebra */
67 Datum RASTER_mapAlgebraExpr(PG_FUNCTION_ARGS);
68 Datum RASTER_mapAlgebraFct(PG_FUNCTION_ARGS);
69
70 /* one-raster neighborhood MapAlgebra */
71 Datum RASTER_mapAlgebraFctNgb(PG_FUNCTION_ARGS);
72
73 /* two-raster MapAlgebra */
74 Datum RASTER_mapAlgebra2(PG_FUNCTION_ARGS);
75
76 /* ---------------------------------------------------------------- */
77 /* n-raster MapAlgebra */
78 /* ---------------------------------------------------------------- */
79
80 #if defined(__clang__)
81 # pragma clang diagnostic push
82 # pragma clang diagnostic ignored "-Wgnu-variable-sized-type-not-at-end"
83 #endif
84
85 typedef struct {
86 Oid ufc_noid;
87 Oid ufc_rettype;
88 FmgrInfo ufl_info;
89 #if POSTGIS_PGSQL_VERSION < 120
90 FunctionCallInfoData ufc_info;
91 #else
92 /* copied from LOCAL_FCINFO in fmgr.h */
93 union {
94 FunctionCallInfoBaseData fcinfo;
95 char fcinfo_data[SizeForFunctionCallInfo(FUNC_MAX_ARGS)]; /* Could be optimized */
96 } ufc_info_data;
97 FunctionCallInfo ufc_info;
98 #endif
99 } rtpg_nmapalgebra_callback_arg;
100
101 #if defined(__clang__)
102 # pragma clang diagnostic pop
103 #endif
104
105 typedef struct rtpg_nmapalgebra_arg_t *rtpg_nmapalgebra_arg;
106 struct rtpg_nmapalgebra_arg_t {
107 int numraster;
108 rt_pgraster **pgraster;
109 rt_raster *raster;
110 uint8_t *isempty; /* flag indicating if raster is empty */
111 uint8_t *ownsdata; /* is the raster self owned or just a pointer to another raster */
112 int *nband; /* source raster's band index, 0-based */
113 uint8_t *hasband; /* does source raster have band at index nband? */
114
115 rt_pixtype pixtype; /* output raster band's pixel type */
116 int hasnodata; /* NODATA flag */
117 double nodataval; /* NODATA value */
118
119 int distance[2]; /* distance in X and Y axis */
120
121 rt_extenttype extenttype; /* ouput raster's extent type */
122 rt_pgraster *pgcextent; /* custom extent of type rt_pgraster */
123 rt_raster cextent; /* custom extent of type rt_raster */
124 rt_mask mask; /* mask for the nmapalgebra operation */
125
126 rtpg_nmapalgebra_callback_arg callback;
127 };
128
rtpg_nmapalgebra_arg_init()129 static rtpg_nmapalgebra_arg rtpg_nmapalgebra_arg_init() {
130 rtpg_nmapalgebra_arg arg = NULL;
131
132 arg = palloc(sizeof(struct rtpg_nmapalgebra_arg_t));
133 if (arg == NULL) {
134 elog(ERROR, "rtpg_nmapalgebra_arg_init: Could not allocate memory for arguments");
135 return NULL;
136 }
137
138 arg->numraster = 0;
139 arg->pgraster = NULL;
140 arg->raster = NULL;
141 arg->isempty = NULL;
142 arg->ownsdata = NULL;
143 arg->nband = NULL;
144 arg->hasband = NULL;
145
146 arg->pixtype = PT_END;
147 arg->hasnodata = 1;
148 arg->nodataval = 0;
149
150 arg->distance[0] = 0;
151 arg->distance[1] = 0;
152
153 arg->extenttype = ET_INTERSECTION;
154 arg->pgcextent = NULL;
155 arg->cextent = NULL;
156 arg->mask = NULL;
157
158 #if POSTGIS_PGSQL_VERSION >= 120
159 arg->callback.ufc_info = &(arg->callback.ufc_info_data.fcinfo);
160 #endif
161 arg->callback.ufc_noid = InvalidOid;
162 arg->callback.ufc_rettype = InvalidOid;
163
164 return arg;
165 }
166
rtpg_nmapalgebra_arg_destroy(rtpg_nmapalgebra_arg arg)167 static void rtpg_nmapalgebra_arg_destroy(rtpg_nmapalgebra_arg arg) {
168 int i = 0;
169
170 if (arg->raster != NULL) {
171 for (i = 0; i < arg->numraster; i++) {
172 if (arg->raster[i] == NULL || !arg->ownsdata[i])
173 continue;
174
175 rt_raster_destroy(arg->raster[i]);
176 }
177
178 pfree(arg->raster);
179 pfree(arg->pgraster);
180 pfree(arg->isempty);
181 pfree(arg->ownsdata);
182 pfree(arg->nband);
183 }
184
185 if (arg->cextent != NULL)
186 rt_raster_destroy(arg->cextent);
187 if( arg->mask != NULL )
188 pfree(arg->mask);
189
190 pfree(arg);
191 }
192
rtpg_nmapalgebra_rastbandarg_process(rtpg_nmapalgebra_arg arg,ArrayType * array,int * allnull,int * allempty,int * noband)193 static int rtpg_nmapalgebra_rastbandarg_process(rtpg_nmapalgebra_arg arg, ArrayType *array, int *allnull, int *allempty, int *noband) {
194 Oid etype;
195 Datum *e;
196 bool *nulls;
197 int16 typlen;
198 bool typbyval;
199 char typalign;
200 int n = 0;
201
202 HeapTupleHeader tup;
203 bool isnull;
204 Datum tupv;
205
206 int i;
207 int j;
208 int nband;
209
210 if (arg == NULL || array == NULL) {
211 elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: NULL values not permitted for parameters");
212 return 0;
213 }
214
215 etype = ARR_ELEMTYPE(array);
216 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
217
218 deconstruct_array(
219 array,
220 etype,
221 typlen, typbyval, typalign,
222 &e, &nulls, &n
223 );
224
225 if (!n) {
226 elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Invalid argument for rastbandarg");
227 return 0;
228 }
229
230 /* prep arg */
231 arg->numraster = n;
232 arg->pgraster = palloc(sizeof(rt_pgraster *) * arg->numraster);
233 arg->raster = palloc(sizeof(rt_raster) * arg->numraster);
234 arg->isempty = palloc(sizeof(uint8_t) * arg->numraster);
235 arg->ownsdata = palloc(sizeof(uint8_t) * arg->numraster);
236 arg->nband = palloc(sizeof(int) * arg->numraster);
237 arg->hasband = palloc(sizeof(uint8_t) * arg->numraster);
238 arg->mask = palloc(sizeof(struct rt_mask_t));
239 if (
240 arg->pgraster == NULL ||
241 arg->raster == NULL ||
242 arg->isempty == NULL ||
243 arg->ownsdata == NULL ||
244 arg->nband == NULL ||
245 arg->hasband == NULL ||
246 arg->mask == NULL
247 ) {
248 elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Could not allocate memory for processing rastbandarg");
249 return 0;
250 }
251
252 *allnull = 0;
253 *allempty = 0;
254 *noband = 0;
255
256 /* process each element */
257 for (i = 0; i < n; i++) {
258 if (nulls[i]) {
259 arg->numraster--;
260 continue;
261 }
262
263 POSTGIS_RT_DEBUGF(4, "Processing rastbandarg at index %d", i);
264
265 arg->raster[i] = NULL;
266 arg->isempty[i] = 0;
267 arg->ownsdata[i] = 1;
268 arg->nband[i] = 0;
269 arg->hasband[i] = 0;
270
271 /* each element is a tuple */
272 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
273 if (NULL == tup) {
274 elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Invalid argument for rastbandarg at index %d", i);
275 return 0;
276 }
277
278 /* first element, raster */
279 POSTGIS_RT_DEBUG(4, "Processing first element (raster)");
280 tupv = GetAttributeByName(tup, "rast", &isnull);
281 if (isnull) {
282 elog(NOTICE, "First argument (nband) of rastbandarg at index %d is NULL. Assuming NULL raster", i);
283 arg->isempty[i] = 1;
284 arg->ownsdata[i] = 0;
285
286 (*allnull)++;
287 (*allempty)++;
288 (*noband)++;
289
290 continue;
291 }
292
293 arg->pgraster[i] = (rt_pgraster *) PG_DETOAST_DATUM(tupv);
294
295 /* see if this is a copy of an existing pgraster */
296 for (j = 0; j < i; j++) {
297 if (!arg->isempty[j] && (arg->pgraster[i] == arg->pgraster[j])) {
298 POSTGIS_RT_DEBUG(4, "raster matching existing same raster found");
299 arg->raster[i] = arg->raster[j];
300 arg->ownsdata[i] = 0;
301 break;
302 }
303 }
304
305 if (arg->ownsdata[i]) {
306 POSTGIS_RT_DEBUG(4, "deserializing raster");
307 arg->raster[i] = rt_raster_deserialize(arg->pgraster[i], FALSE);
308 if (arg->raster[i] == NULL) {
309 elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Could not deserialize raster at index %d", i);
310 return 0;
311 }
312 }
313
314 /* is raster empty? */
315 arg->isempty[i] = rt_raster_is_empty(arg->raster[i]);
316 if (arg->isempty[i]) {
317 (*allempty)++;
318 (*noband)++;
319
320 continue;
321 }
322
323 /* second element, nband */
324 POSTGIS_RT_DEBUG(4, "Processing second element (nband)");
325 tupv = GetAttributeByName(tup, "nband", &isnull);
326 if (isnull) {
327 nband = 1;
328 elog(NOTICE, "First argument (nband) of rastbandarg at index %d is NULL. Assuming nband = %d", i, nband);
329 }
330 else
331 nband = DatumGetInt32(tupv);
332
333 if (nband < 1) {
334 elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Band number provided for rastbandarg at index %d must be greater than zero (1-based)", i);
335 return 0;
336 }
337
338 arg->nband[i] = nband - 1;
339 arg->hasband[i] = rt_raster_has_band(arg->raster[i], arg->nband[i]);
340 if (!arg->hasband[i]) {
341 (*noband)++;
342 POSTGIS_RT_DEBUGF(4, "Band at index %d not found in raster", nband);
343 }
344 }
345
346 if (arg->numraster < n) {
347 arg->pgraster = repalloc(arg->pgraster, sizeof(rt_pgraster *) * arg->numraster);
348 arg->raster = repalloc(arg->raster, sizeof(rt_raster) * arg->numraster);
349 arg->isempty = repalloc(arg->isempty, sizeof(uint8_t) * arg->numraster);
350 arg->ownsdata = repalloc(arg->ownsdata, sizeof(uint8_t) * arg->numraster);
351 arg->nband = repalloc(arg->nband, sizeof(int) * arg->numraster);
352 arg->hasband = repalloc(arg->hasband, sizeof(uint8_t) * arg->numraster);
353 if (
354 arg->pgraster == NULL ||
355 arg->raster == NULL ||
356 arg->isempty == NULL ||
357 arg->ownsdata == NULL ||
358 arg->nband == NULL ||
359 arg->hasband == NULL
360 ) {
361 elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Could not reallocate memory for processed rastbandarg");
362 return 0;
363 }
364 }
365
366 POSTGIS_RT_DEBUGF(4, "arg->numraster = %d", arg->numraster);
367
368 return 1;
369 }
370
371 /*
372 Callback for RASTER_nMapAlgebra
373 */
rtpg_nmapalgebra_callback(rt_iterator_arg arg,void * userarg,double * value,int * nodata)374 static int rtpg_nmapalgebra_callback(
375 rt_iterator_arg arg, void *userarg,
376 double *value, int *nodata
377 ) {
378 rtpg_nmapalgebra_callback_arg *callback = (rtpg_nmapalgebra_callback_arg *) userarg;
379
380 int16 typlen;
381 bool typbyval;
382 char typalign;
383
384 ArrayType *mdValues = NULL;
385 Datum *_values = NULL;
386 bool *_nodata = NULL;
387
388 ArrayType *mdPos = NULL;
389 Datum *_pos = NULL;
390 bool *_null = NULL;
391
392 int i = 0;
393 uint32_t x = 0;
394 uint32_t y = 0;
395 int z = 0;
396 int dim[3] = {0};
397 int lbound[3] = {1, 1, 1};
398 Datum datum = (Datum) NULL;
399
400 if (arg == NULL)
401 return 0;
402
403 *value = 0;
404 *nodata = 0;
405
406 dim[0] = arg->rasters;
407 dim[1] = arg->rows;
408 dim[2] = arg->columns;
409
410 _values = palloc(sizeof(Datum) * arg->rasters * arg->rows * arg->columns);
411 _nodata = palloc(sizeof(bool) * arg->rasters * arg->rows * arg->columns);
412 if (_values == NULL || _nodata == NULL) {
413 elog(ERROR, "rtpg_nmapalgebra_callback: Could not allocate memory for values array");
414 return 0;
415 }
416
417 /* build mdValues */
418 i = 0;
419 /* raster */
420 for (z = 0; z < arg->rasters; z++) {
421 /* Y axis */
422 for (y = 0; y < arg->rows; y++) {
423 /* X axis */
424 for (x = 0; x < arg->columns; x++) {
425 POSTGIS_RT_DEBUGF(4, "(z, y ,x) = (%d, %d, %d)", z, y, x);
426 POSTGIS_RT_DEBUGF(4, "(value, nodata) = (%f, %d)", arg->values[z][y][x], arg->nodata[z][y][x]);
427
428 _nodata[i] = (bool) arg->nodata[z][y][x];
429 if (!_nodata[i])
430 _values[i] = Float8GetDatum(arg->values[z][y][x]);
431 else
432 _values[i] = (Datum) NULL;
433
434 i++;
435 }
436 }
437 }
438
439 /* info about the type of item in the multi-dimensional array (float8). */
440 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
441
442 /* construct mdValues */
443 mdValues = construct_md_array(
444 _values, _nodata,
445 3, dim, lbound,
446 FLOAT8OID,
447 typlen, typbyval, typalign
448 );
449 pfree(_nodata);
450 pfree(_values);
451
452 _pos = palloc(sizeof(Datum) * (arg->rasters + 1) * 2);
453 _null = palloc(sizeof(bool) * (arg->rasters + 1) * 2);
454 if (_pos == NULL || _null == NULL) {
455 pfree(mdValues);
456 elog(ERROR, "rtpg_nmapalgebra_callback: Could not allocate memory for position array");
457 return 0;
458 }
459 memset(_null, 0, sizeof(bool) * (arg->rasters + 1) * 2);
460
461 /* build mdPos */
462 i = 0;
463 _pos[i] = arg->dst_pixel[0] + 1;
464 i++;
465 _pos[i] = arg->dst_pixel[1] + 1;
466 i++;
467
468 for (z = 0; z < arg->rasters; z++) {
469 _pos[i] = (Datum)arg->src_pixel[z][0] + 1;
470 i++;
471
472 _pos[i] = (Datum)arg->src_pixel[z][1] + 1;
473 i++;
474 }
475
476 /* info about the type of item in the multi-dimensional array (int4). */
477 get_typlenbyvalalign(INT4OID, &typlen, &typbyval, &typalign);
478
479 /* reuse dim and lbound, just tweak to what we need */
480 dim[0] = arg->rasters + 1;
481 dim[1] = 2;
482 lbound[0] = 0;
483
484 /* construct mdPos */
485 mdPos = construct_md_array(
486 _pos, _null,
487 2, dim, lbound,
488 INT4OID,
489 typlen, typbyval, typalign
490 );
491 pfree(_pos);
492 pfree(_null);
493
494 #if POSTGIS_PGSQL_VERSION < 120
495 callback->ufc_info.arg[0] = PointerGetDatum(mdValues);
496 callback->ufc_info.arg[1] = PointerGetDatum(mdPos);
497 #else
498 callback->ufc_info->args[0].value = PointerGetDatum(mdValues);
499 callback->ufc_info->args[1].value = PointerGetDatum(mdPos);
500 #endif
501
502 /* call user callback function */
503 #if POSTGIS_PGSQL_VERSION < 120
504 datum = FunctionCallInvoke(&(callback->ufc_info));
505 #else
506 datum = FunctionCallInvoke(callback->ufc_info);
507 #endif
508 pfree(mdValues);
509 pfree(mdPos);
510
511 /* result is not null*/
512 #if POSTGIS_PGSQL_VERSION < 120
513 if (!callback->ufc_info.isnull) {
514 #else
515 if (!callback->ufc_info->isnull)
516 {
517 #endif
518 switch (callback->ufc_rettype) {
519 case FLOAT8OID:
520 *value = DatumGetFloat8(datum);
521 break;
522 case FLOAT4OID:
523 *value = (double) DatumGetFloat4(datum);
524 break;
525 case INT4OID:
526 *value = (double) DatumGetInt32(datum);
527 break;
528 case INT2OID:
529 *value = (double) DatumGetInt16(datum);
530 break;
531 }
532 }
533 else
534 *nodata = 1;
535
536 return 1;
537 }
538
539 /*
540 ST_MapAlgebra for n rasters
541 */
542 PG_FUNCTION_INFO_V1(RASTER_nMapAlgebra);
543 Datum RASTER_nMapAlgebra(PG_FUNCTION_ARGS)
544 {
545 rtpg_nmapalgebra_arg arg = NULL;
546 rt_iterator itrset;
547 ArrayType *maskArray;
548 Oid etype;
549 Datum *maskElements;
550 bool *maskNulls;
551 int16 typlen;
552 bool typbyval;
553 char typalign;
554 int ndims = 0;
555 int num;
556 int *maskDims;
557 int x,y;
558
559
560 int i = 0;
561 int noerr = 0;
562 int allnull = 0;
563 int allempty = 0;
564 int noband = 0;
565
566 rt_raster raster = NULL;
567 rt_band band = NULL;
568 rt_pgraster *pgraster = NULL;
569
570 POSTGIS_RT_DEBUG(3, "Starting...");
571
572 if (PG_ARGISNULL(0))
573 PG_RETURN_NULL();
574
575 /* init argument struct */
576 arg = rtpg_nmapalgebra_arg_init();
577 if (arg == NULL) {
578 elog(ERROR, "RASTER_nMapAlgebra: Could not initialize argument structure");
579 PG_RETURN_NULL();
580 }
581
582 /* let helper function process rastbandarg (0) */
583 if (!rtpg_nmapalgebra_rastbandarg_process(arg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
584 rtpg_nmapalgebra_arg_destroy(arg);
585 elog(ERROR, "RASTER_nMapAlgebra: Could not process rastbandarg");
586 PG_RETURN_NULL();
587 }
588
589 POSTGIS_RT_DEBUGF(4, "allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
590
591 /* all rasters are NULL, return NULL */
592 if (allnull == arg->numraster) {
593 elog(NOTICE, "All input rasters are NULL. Returning NULL");
594 rtpg_nmapalgebra_arg_destroy(arg);
595 PG_RETURN_NULL();
596 }
597
598 /* pixel type (2) */
599 if (!PG_ARGISNULL(2)) {
600 char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
601
602 /* Get the pixel type index */
603 arg->pixtype = rt_pixtype_index_from_name(pixtypename);
604 if (arg->pixtype == PT_END) {
605 rtpg_nmapalgebra_arg_destroy(arg);
606 elog(ERROR, "RASTER_nMapAlgebra: Invalid pixel type: %s", pixtypename);
607 PG_RETURN_NULL();
608 }
609 }
610
611 /* distancex (3) */
612 if (!PG_ARGISNULL(3)){
613 arg->distance[0] = PG_GETARG_INT32(3);
614 }else{
615 arg->distance[0] = 0;
616 }
617 /* distancey (4) */
618 if (!PG_ARGISNULL(4)){
619 arg->distance[1] = PG_GETARG_INT32(4);
620 }else{
621 arg->distance[1] = 0;
622 }
623 if (arg->distance[0] < 0 || arg->distance[1] < 0) {
624 rtpg_nmapalgebra_arg_destroy(arg);
625 elog(ERROR, "RASTER_nMapAlgebra: Distance for X and Y axis must be greater than or equal to zero");
626 PG_RETURN_NULL();
627 }
628
629 /* extent type (5) */
630 if (!PG_ARGISNULL(5)) {
631 char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(5))));
632 arg->extenttype = rt_util_extent_type(extenttypename);
633 }
634 POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->extenttype);
635
636 /* custom extent (6) */
637 if (arg->extenttype == ET_CUSTOM) {
638 if (PG_ARGISNULL(6)) {
639 elog(NOTICE, "Custom extent is NULL. Returning NULL");
640 rtpg_nmapalgebra_arg_destroy(arg);
641 PG_RETURN_NULL();
642 }
643
644 arg->pgcextent = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(6));
645
646 /* only need the raster header */
647 arg->cextent = rt_raster_deserialize(arg->pgcextent, TRUE);
648 if (arg->cextent == NULL) {
649 rtpg_nmapalgebra_arg_destroy(arg);
650 elog(ERROR, "RASTER_nMapAlgebra: Could not deserialize custom extent");
651 PG_RETURN_NULL();
652 }
653 else if (rt_raster_is_empty(arg->cextent)) {
654 elog(NOTICE, "Custom extent is an empty raster. Returning empty raster");
655 rtpg_nmapalgebra_arg_destroy(arg);
656
657 raster = rt_raster_new(0, 0);
658 if (raster == NULL) {
659 elog(ERROR, "RASTER_nMapAlgebra: Could not create empty raster");
660 PG_RETURN_NULL();
661 }
662
663 pgraster = rt_raster_serialize(raster);
664 rt_raster_destroy(raster);
665 if (!pgraster) PG_RETURN_NULL();
666
667 SET_VARSIZE(pgraster, pgraster->size);
668 PG_RETURN_POINTER(pgraster);
669 }
670 }
671
672 /* mask (7) */
673
674 if( PG_ARGISNULL(7) ){
675 pfree(arg->mask);
676 arg->mask = NULL;
677 }
678 else {
679 maskArray = PG_GETARG_ARRAYTYPE_P(7);
680 etype = ARR_ELEMTYPE(maskArray);
681 get_typlenbyvalalign(etype,&typlen,&typbyval,&typalign);
682
683 switch (etype) {
684 case FLOAT4OID:
685 case FLOAT8OID:
686 break;
687 default:
688 rtpg_nmapalgebra_arg_destroy(arg);
689 elog(ERROR,"RASTER_nMapAlgebra: Mask data type must be FLOAT8 or FLOAT4");
690 PG_RETURN_NULL();
691 }
692
693 ndims = ARR_NDIM(maskArray);
694
695 if (ndims != 2) {
696 elog(ERROR, "RASTER_nMapAlgebra: Mask Must be a 2D array");
697 rtpg_nmapalgebra_arg_destroy(arg);
698 PG_RETURN_NULL();
699 }
700
701 maskDims = ARR_DIMS(maskArray);
702
703 if (maskDims[0] % 2 == 0 || maskDims[1] % 2 == 0) {
704 elog(ERROR,"RASTER_nMapAlgebra: Mask dimensions must be odd");
705 rtpg_nmapalgebra_arg_destroy(arg);
706 PG_RETURN_NULL();
707 }
708
709 deconstruct_array(
710 maskArray,
711 etype,
712 typlen, typbyval,typalign,
713 &maskElements,&maskNulls,&num
714 );
715
716 if (num < 1 || num != (maskDims[0] * maskDims[1])) {
717 if (num) {
718 pfree(maskElements);
719 pfree(maskNulls);
720 }
721 elog(ERROR, "RASTER_nMapAlgebra: Could not deconstruct new values array");
722 rtpg_nmapalgebra_arg_destroy(arg);
723 PG_RETURN_NULL();
724 }
725
726 /* allocate mem for mask array */
727 arg->mask->values = palloc(sizeof(double*)* maskDims[0]);
728 arg->mask->nodata = palloc(sizeof(int*)*maskDims[0]);
729 for (i = 0; i < maskDims[0]; i++) {
730 arg->mask->values[i] = (double*) palloc(sizeof(double) * maskDims[1]);
731 arg->mask->nodata[i] = (int*) palloc(sizeof(int) * maskDims[1]);
732 }
733
734 /* place values in to mask */
735 i = 0;
736 for (y = 0; y < maskDims[0]; y++) {
737 for (x = 0; x < maskDims[1]; x++) {
738 if (maskNulls[i]) {
739 arg->mask->values[y][x] = 0;
740 arg->mask->nodata[y][x] = 1;
741 }
742 else {
743 switch (etype) {
744 case FLOAT4OID:
745 arg->mask->values[y][x] = (double) DatumGetFloat4(maskElements[i]);
746 arg->mask->nodata[y][x] = 0;
747 break;
748 case FLOAT8OID:
749 arg->mask->values[y][x] = (double) DatumGetFloat8(maskElements[i]);
750 arg->mask->nodata[y][x] = 0;
751 }
752 }
753 i++;
754 }
755 }
756
757 /*set mask dimensions*/
758 arg->mask->dimx = maskDims[0];
759 arg->mask->dimy = maskDims[1];
760 if (maskDims[0] == 1 && maskDims[1] == 1) {
761 arg->distance[0] = 0;
762 arg->distance[1] = 0;
763 }
764 else {
765 arg->distance[0] = maskDims[0] % 2;
766 arg->distance[1] = maskDims[1] % 2;
767 }
768 }/*end if else argisnull*/
769
770 /* (8) weighted boolean */
771 if (PG_ARGISNULL(8) || !PG_GETARG_BOOL(8)) {
772 if (arg->mask != NULL)
773 arg->mask->weighted = 0;
774 }else{
775 if(arg->mask !=NULL)
776 arg->mask->weighted = 1;
777 }
778
779 noerr = 1;
780
781 /* all rasters are empty, return empty raster */
782 if (allempty == arg->numraster) {
783 elog(NOTICE, "All input rasters are empty. Returning empty raster");
784 noerr = 0;
785 }
786 /* all rasters don't have indicated band, return empty raster */
787 else if (noband == arg->numraster) {
788 elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
789 noerr = 0;
790 }
791 if (!noerr) {
792 rtpg_nmapalgebra_arg_destroy(arg);
793
794 raster = rt_raster_new(0, 0);
795 if (raster == NULL) {
796 elog(ERROR, "RASTER_nMapAlgebra: Could not create empty raster");
797 PG_RETURN_NULL();
798 }
799
800 pgraster = rt_raster_serialize(raster);
801 rt_raster_destroy(raster);
802 if (!pgraster) PG_RETURN_NULL();
803
804 SET_VARSIZE(pgraster, pgraster->size);
805 PG_RETURN_POINTER(pgraster);
806 }
807
808 /* do regprocedure last (1) */
809 if (!PG_ARGISNULL(1) || get_fn_expr_argtype(fcinfo->flinfo, 1) == REGPROCEDUREOID) {
810 POSTGIS_RT_DEBUG(4, "processing callbackfunc");
811 arg->callback.ufc_noid = PG_GETARG_OID(1);
812
813 /* get function info */
814 fmgr_info(arg->callback.ufc_noid, &(arg->callback.ufl_info));
815
816 /* function cannot return set */
817 noerr = 0;
818 if (arg->callback.ufl_info.fn_retset) {
819 noerr = 1;
820 }
821 /* function should have correct # of args */
822 else if (arg->callback.ufl_info.fn_nargs != 3) {
823 noerr = 2;
824 }
825
826 /* check that callback function return type is supported */
827 if (
828 get_func_result_type(
829 arg->callback.ufc_noid,
830 &(arg->callback.ufc_rettype),
831 NULL
832 ) != TYPEFUNC_SCALAR
833 ) {
834 noerr = 3;
835 }
836
837 if (!(
838 arg->callback.ufc_rettype == FLOAT8OID ||
839 arg->callback.ufc_rettype == FLOAT4OID ||
840 arg->callback.ufc_rettype == INT4OID ||
841 arg->callback.ufc_rettype == INT2OID
842 )) {
843 noerr = 4;
844 }
845
846 /*
847 TODO: consider adding checks of the userfunction parameters
848 should be able to use get_fn_expr_argtype() of fmgr.c
849 */
850
851 if (noerr != 0) {
852 rtpg_nmapalgebra_arg_destroy(arg);
853 switch (noerr) {
854 case 4:
855 elog(ERROR, "RASTER_nMapAlgebra: Function provided must return a double precision, float, int or smallint");
856 break;
857 case 3:
858 elog(ERROR, "RASTER_nMapAlgebra: Function provided must return scalar (double precision, float, int, smallint)");
859 break;
860 case 2:
861 elog(ERROR, "RASTER_nMapAlgebra: Function provided must have three input parameters");
862 break;
863 case 1:
864 elog(ERROR, "RASTER_nMapAlgebra: Function provided must return double precision, not resultset");
865 break;
866 }
867 PG_RETURN_NULL();
868 }
869
870 if (func_volatile(arg->callback.ufc_noid) == 'v')
871 elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
872
873 /* prep function call data */
874 #if POSTGIS_PGSQL_VERSION < 120
875 InitFunctionCallInfoData(arg->callback.ufc_info, &(arg->callback.ufl_info), arg->callback.ufl_info.fn_nargs, InvalidOid, NULL, NULL);
876
877 memset(arg->callback.ufc_info.argnull, FALSE, sizeof(bool) * arg->callback.ufl_info.fn_nargs);
878 #else
879 InitFunctionCallInfoData(*(arg->callback.ufc_info),
880 &(arg->callback.ufl_info),
881 arg->callback.ufl_info.fn_nargs,
882 InvalidOid,
883 NULL,
884 NULL);
885
886 arg->callback.ufc_info->args[0].isnull = FALSE;
887 arg->callback.ufc_info->args[1].isnull = FALSE;
888 arg->callback.ufc_info->args[2].isnull = FALSE;
889 #endif
890
891 /* userargs (7) */
892 if (!PG_ARGISNULL(9))
893 #if POSTGIS_PGSQL_VERSION < 120
894 arg->callback.ufc_info.arg[2] = PG_GETARG_DATUM(9);
895 #else
896 arg->callback.ufc_info->args[2].value = PG_GETARG_DATUM(9);
897 #endif
898 else {
899 if (arg->callback.ufl_info.fn_strict) {
900 /* build and assign an empty TEXT array */
901 /* TODO: manually free the empty array? */
902 #if POSTGIS_PGSQL_VERSION < 120
903 arg->callback.ufc_info.arg[2] = PointerGetDatum(
904 construct_empty_array(TEXTOID)
905 );
906 arg->callback.ufc_info.argnull[2] = FALSE;
907 #else
908 arg->callback.ufc_info->args[2].value = PointerGetDatum(construct_empty_array(TEXTOID));
909 arg->callback.ufc_info->args[2].isnull = FALSE;
910 #endif
911 }
912 else {
913 #if POSTGIS_PGSQL_VERSION < 120
914 arg->callback.ufc_info.arg[2] = (Datum) NULL;
915 arg->callback.ufc_info.argnull[2] = TRUE;
916 #else
917 arg->callback.ufc_info->args[2].value = (Datum)NULL;
918 arg->callback.ufc_info->args[2].isnull = TRUE;
919 #endif
920 }
921 }
922 }
923 else {
924 rtpg_nmapalgebra_arg_destroy(arg);
925 elog(ERROR, "RASTER_nMapAlgebra: callbackfunc must be provided");
926 PG_RETURN_NULL();
927 }
928
929 /* determine nodataval and possibly pixtype */
930 /* band to check */
931 switch (arg->extenttype) {
932 case ET_LAST:
933 i = arg->numraster - 1;
934 break;
935 case ET_SECOND:
936 i = (arg->numraster > 1) ? 1 : 0;
937 break;
938 default:
939 i = 0;
940 break;
941 }
942 /* find first viable band */
943 if (!arg->hasband[i]) {
944 for (i = 0; i < arg->numraster; i++) {
945 if (arg->hasband[i])
946 break;
947 }
948 if (i >= arg->numraster)
949 i = arg->numraster - 1;
950 }
951 band = rt_raster_get_band(arg->raster[i], arg->nband[i]);
952
953 /* set pixel type if PT_END */
954 if (arg->pixtype == PT_END)
955 arg->pixtype = rt_band_get_pixtype(band);
956
957 /* set hasnodata and nodataval */
958 arg->hasnodata = 1;
959 if (rt_band_get_hasnodata_flag(band))
960 rt_band_get_nodata(band, &(arg->nodataval));
961 else
962 arg->nodataval = rt_band_get_min_value(band);
963
964 POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->pixtype), arg->hasnodata, arg->nodataval);
965
966 /* init itrset */
967 itrset = palloc(sizeof(struct rt_iterator_t) * arg->numraster);
968 if (itrset == NULL) {
969 rtpg_nmapalgebra_arg_destroy(arg);
970 elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
971 PG_RETURN_NULL();
972 }
973
974 /* set itrset */
975 for (i = 0; i < arg->numraster; i++) {
976 itrset[i].raster = arg->raster[i];
977 itrset[i].nband = arg->nband[i];
978 itrset[i].nbnodata = 1;
979 }
980
981 /* pass everything to iterator */
982 noerr = rt_raster_iterator(
983 itrset, arg->numraster,
984 arg->extenttype, arg->cextent,
985 arg->pixtype,
986 arg->hasnodata, arg->nodataval,
987 arg->distance[0], arg->distance[1],
988 arg->mask,
989 &(arg->callback),
990 rtpg_nmapalgebra_callback,
991 &raster
992 );
993
994 /* cleanup */
995 pfree(itrset);
996 rtpg_nmapalgebra_arg_destroy(arg);
997
998 if (noerr != ES_NONE) {
999 elog(ERROR, "RASTER_nMapAlgebra: Could not run raster iterator function");
1000 PG_RETURN_NULL();
1001 }
1002 else if (raster == NULL)
1003 PG_RETURN_NULL();
1004
1005 pgraster = rt_raster_serialize(raster);
1006 rt_raster_destroy(raster);
1007
1008 POSTGIS_RT_DEBUG(3, "Finished");
1009
1010 if (!pgraster)
1011 PG_RETURN_NULL();
1012
1013 SET_VARSIZE(pgraster, pgraster->size);
1014 PG_RETURN_POINTER(pgraster);
1015 }
1016
1017 /* ---------------------------------------------------------------- */
1018 /* expression ST_MapAlgebra for n rasters */
1019 /* ---------------------------------------------------------------- */
1020
1021 typedef struct {
1022 int exprcount;
1023
1024 struct {
1025 SPIPlanPtr spi_plan;
1026 uint32_t spi_argcount;
1027 uint8_t *spi_argpos;
1028
1029 int hasval;
1030 double val;
1031 } expr[3];
1032
1033 struct {
1034 int hasval;
1035 double val;
1036 } nodatanodata;
1037
1038 struct {
1039 int count;
1040 char **val;
1041 } kw;
1042
1043 } rtpg_nmapalgebraexpr_callback_arg;
1044
1045 typedef struct rtpg_nmapalgebraexpr_arg_t *rtpg_nmapalgebraexpr_arg;
1046 struct rtpg_nmapalgebraexpr_arg_t {
1047 rtpg_nmapalgebra_arg bandarg;
1048
1049 rtpg_nmapalgebraexpr_callback_arg callback;
1050 };
1051
1052 static rtpg_nmapalgebraexpr_arg rtpg_nmapalgebraexpr_arg_init(int cnt, char **kw) {
1053 rtpg_nmapalgebraexpr_arg arg = NULL;
1054 int i = 0;
1055
1056 arg = palloc(sizeof(struct rtpg_nmapalgebraexpr_arg_t));
1057 if (arg == NULL) {
1058 elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arguments");
1059 return NULL;
1060 }
1061
1062 arg->bandarg = rtpg_nmapalgebra_arg_init();
1063 if (arg->bandarg == NULL) {
1064 elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arg->bandarg");
1065 return NULL;
1066 }
1067
1068 arg->callback.kw.count = cnt;
1069 arg->callback.kw.val = kw;
1070
1071 arg->callback.exprcount = 3;
1072 for (i = 0; i < arg->callback.exprcount; i++) {
1073 arg->callback.expr[i].spi_plan = NULL;
1074 arg->callback.expr[i].spi_argcount = 0;
1075 arg->callback.expr[i].spi_argpos = palloc(cnt * sizeof(uint8_t));
1076 if (arg->callback.expr[i].spi_argpos == NULL) {
1077 elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for spi_argpos");
1078 return NULL;
1079 }
1080 memset(arg->callback.expr[i].spi_argpos, 0, sizeof(uint8_t) * cnt);
1081 arg->callback.expr[i].hasval = 0;
1082 arg->callback.expr[i].val = 0;
1083 }
1084
1085 arg->callback.nodatanodata.hasval = 0;
1086 arg->callback.nodatanodata.val = 0;
1087
1088 return arg;
1089 }
1090
1091 static void rtpg_nmapalgebraexpr_arg_destroy(rtpg_nmapalgebraexpr_arg arg) {
1092 int i = 0;
1093
1094 rtpg_nmapalgebra_arg_destroy(arg->bandarg);
1095
1096 for (i = 0; i < arg->callback.exprcount; i++) {
1097 if (arg->callback.expr[i].spi_plan)
1098 SPI_freeplan(arg->callback.expr[i].spi_plan);
1099 if (arg->callback.kw.count)
1100 pfree(arg->callback.expr[i].spi_argpos);
1101 }
1102
1103 pfree(arg);
1104 }
1105
1106 static int rtpg_nmapalgebraexpr_callback(
1107 rt_iterator_arg arg, void *userarg,
1108 double *value, int *nodata
1109 ) {
1110 rtpg_nmapalgebraexpr_callback_arg *callback = (rtpg_nmapalgebraexpr_callback_arg *) userarg;
1111 SPIPlanPtr plan = NULL;
1112 int i = 0;
1113 uint8_t id = 0;
1114
1115 if (arg == NULL)
1116 return 0;
1117
1118 *value = 0;
1119 *nodata = 0;
1120
1121 /* 2 raster */
1122 if (arg->rasters > 1) {
1123 /* nodata1 = 1 AND nodata2 = 1, nodatanodataval */
1124 if (arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1125 if (callback->nodatanodata.hasval)
1126 *value = callback->nodatanodata.val;
1127 else
1128 *nodata = 1;
1129 }
1130 /* nodata1 = 1 AND nodata2 != 1, nodata1expr */
1131 else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0]) {
1132 id = 1;
1133 if (callback->expr[id].hasval)
1134 *value = callback->expr[id].val;
1135 else if (callback->expr[id].spi_plan)
1136 plan = callback->expr[id].spi_plan;
1137 else
1138 *nodata = 1;
1139 }
1140 /* nodata1 != 1 AND nodata2 = 1, nodata2expr */
1141 else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1142 id = 2;
1143 if (callback->expr[id].hasval)
1144 *value = callback->expr[id].val;
1145 else if (callback->expr[id].spi_plan)
1146 plan = callback->expr[id].spi_plan;
1147 else
1148 *nodata = 1;
1149 }
1150 /* expression */
1151 else {
1152 id = 0;
1153 if (callback->expr[id].hasval)
1154 *value = callback->expr[id].val;
1155 else if (callback->expr[id].spi_plan)
1156 plan = callback->expr[id].spi_plan;
1157 else {
1158 if (callback->nodatanodata.hasval)
1159 *value = callback->nodatanodata.val;
1160 else
1161 *nodata = 1;
1162 }
1163 }
1164 }
1165 /* 1 raster */
1166 else {
1167 /* nodata = 1, nodata1expr */
1168 if (arg->nodata[0][0][0]) {
1169 id = 1;
1170 if (callback->expr[id].hasval)
1171 *value = callback->expr[id].val;
1172 else if (callback->expr[id].spi_plan)
1173 plan = callback->expr[id].spi_plan;
1174 else
1175 *nodata = 1;
1176 }
1177 /* expression */
1178 else {
1179 id = 0;
1180 if (callback->expr[id].hasval)
1181 *value = callback->expr[id].val;
1182 else if (callback->expr[id].spi_plan)
1183 plan = callback->expr[id].spi_plan;
1184 else {
1185 /* see if nodata1expr is available */
1186 id = 1;
1187 if (callback->expr[id].hasval)
1188 *value = callback->expr[id].val;
1189 else if (callback->expr[id].spi_plan)
1190 plan = callback->expr[id].spi_plan;
1191 else
1192 *nodata = 1;
1193 }
1194 }
1195 }
1196
1197 /* run prepared plan */
1198 if (plan != NULL) {
1199 Datum values[12];
1200 char nulls[12];
1201 int err = 0;
1202
1203 TupleDesc tupdesc;
1204 SPITupleTable *tuptable = NULL;
1205 HeapTuple tuple;
1206 Datum datum;
1207 bool isnull = FALSE;
1208
1209 POSTGIS_RT_DEBUGF(4, "Running plan %d", id);
1210
1211 /* init values and nulls */
1212 memset(values, (Datum) NULL, sizeof(Datum) * callback->kw.count);
1213 memset(nulls, FALSE, sizeof(char) * callback->kw.count);
1214
1215 if (callback->expr[id].spi_argcount) {
1216 int idx = 0;
1217
1218 for (i = 0; i < callback->kw.count; i++) {
1219 idx = callback->expr[id].spi_argpos[i];
1220 if (idx < 1) continue;
1221 idx--; /* 1-based now 0-based */
1222
1223 switch (i) {
1224 /* [rast.x] */
1225 case 0:
1226 values[idx] = Int32GetDatum(arg->src_pixel[0][0] + 1);
1227 break;
1228 /* [rast.y] */
1229 case 1:
1230 values[idx] = Int32GetDatum(arg->src_pixel[0][1] + 1);
1231 break;
1232 /* [rast.val] */
1233 case 2:
1234 /* [rast] */
1235 case 3:
1236 if (!arg->nodata[0][0][0])
1237 values[idx] = Float8GetDatum(arg->values[0][0][0]);
1238 else
1239 nulls[idx] = TRUE;
1240 break;
1241
1242 /* [rast1.x] */
1243 case 4:
1244 values[idx] = Int32GetDatum(arg->src_pixel[0][0] + 1);
1245 break;
1246 /* [rast1.y] */
1247 case 5:
1248 values[idx] = Int32GetDatum(arg->src_pixel[0][1] + 1);
1249 break;
1250 /* [rast1.val] */
1251 case 6:
1252 /* [rast1] */
1253 case 7:
1254 if (!arg->nodata[0][0][0])
1255 values[idx] = Float8GetDatum(arg->values[0][0][0]);
1256 else
1257 nulls[idx] = TRUE;
1258 break;
1259
1260 /* [rast2.x] */
1261 case 8:
1262 values[idx] = Int32GetDatum(arg->src_pixel[1][0] + 1);
1263 break;
1264 /* [rast2.y] */
1265 case 9:
1266 values[idx] = Int32GetDatum(arg->src_pixel[1][1] + 1);
1267 break;
1268 /* [rast2.val] */
1269 case 10:
1270 /* [rast2] */
1271 case 11:
1272 if (!arg->nodata[1][0][0])
1273 values[idx] = Float8GetDatum(arg->values[1][0][0]);
1274 else
1275 nulls[idx] = TRUE;
1276 break;
1277 }
1278
1279 }
1280 }
1281
1282 /* run prepared plan */
1283 err = SPI_execute_plan(plan, values, nulls, TRUE, 1);
1284 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1285 elog(ERROR, "rtpg_nmapalgebraexpr_callback: Unexpected error when running prepared statement %d", id);
1286 return 0;
1287 }
1288
1289 /* get output of prepared plan */
1290 tupdesc = SPI_tuptable->tupdesc;
1291 tuptable = SPI_tuptable;
1292 tuple = tuptable->vals[0];
1293
1294 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1295 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1296 if (SPI_tuptable) SPI_freetuptable(tuptable);
1297 elog(ERROR, "rtpg_nmapalgebraexpr_callback: Could not get result of prepared statement %d", id);
1298 return 0;
1299 }
1300
1301 if (!isnull) {
1302 *value = DatumGetFloat8(datum);
1303 POSTGIS_RT_DEBUG(4, "Getting value from Datum");
1304 }
1305 else {
1306 /* 2 raster, check nodatanodataval */
1307 if (arg->rasters > 1) {
1308 if (callback->nodatanodata.hasval)
1309 *value = callback->nodatanodata.val;
1310 else
1311 *nodata = 1;
1312 }
1313 /* 1 raster, check nodataval */
1314 else {
1315 if (callback->expr[1].hasval)
1316 *value = callback->expr[1].val;
1317 else
1318 *nodata = 1;
1319 }
1320 }
1321
1322 if (SPI_tuptable) SPI_freetuptable(tuptable);
1323 }
1324
1325 POSTGIS_RT_DEBUGF(4, "(value, nodata) = (%f, %d)", *value, *nodata);
1326 return 1;
1327 }
1328
1329 PG_FUNCTION_INFO_V1(RASTER_nMapAlgebraExpr);
1330 Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS)
1331 {
1332 MemoryContext mainMemCtx = CurrentMemoryContext;
1333 rtpg_nmapalgebraexpr_arg arg = NULL;
1334 rt_iterator itrset;
1335 uint16_t exprpos[3] = {1, 4, 5};
1336
1337 int i = 0;
1338 int j = 0;
1339 int k = 0;
1340
1341 int numraster = 0;
1342 int err = 0;
1343 int allnull = 0;
1344 int allempty = 0;
1345 int noband = 0;
1346 int len = 0;
1347
1348 TupleDesc tupdesc;
1349 SPITupleTable *tuptable = NULL;
1350 HeapTuple tuple;
1351 Datum datum;
1352 bool isnull = FALSE;
1353
1354 rt_raster raster = NULL;
1355 rt_band band = NULL;
1356 rt_pgraster *pgraster = NULL;
1357
1358 const int argkwcount = 12;
1359 char *argkw[] = {
1360 "[rast.x]",
1361 "[rast.y]",
1362 "[rast.val]",
1363 "[rast]",
1364 "[rast1.x]",
1365 "[rast1.y]",
1366 "[rast1.val]",
1367 "[rast1]",
1368 "[rast2.x]",
1369 "[rast2.y]",
1370 "[rast2.val]",
1371 "[rast2]"
1372 };
1373
1374 if (PG_ARGISNULL(0))
1375 PG_RETURN_NULL();
1376
1377 /* init argument struct */
1378 arg = rtpg_nmapalgebraexpr_arg_init(argkwcount, argkw);
1379 if (arg == NULL) {
1380 elog(ERROR, "RASTER_nMapAlgebraExpr: Could not initialize argument structure");
1381 PG_RETURN_NULL();
1382 }
1383
1384 /* let helper function process rastbandarg (0) */
1385 if (!rtpg_nmapalgebra_rastbandarg_process(arg->bandarg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
1386 rtpg_nmapalgebraexpr_arg_destroy(arg);
1387 elog(ERROR, "RASTER_nMapAlgebra: Could not process rastbandarg");
1388 PG_RETURN_NULL();
1389 }
1390
1391 POSTGIS_RT_DEBUGF(4, "allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
1392
1393 /* all rasters are NULL, return NULL */
1394 if (allnull == arg->bandarg->numraster) {
1395 elog(NOTICE, "All input rasters are NULL. Returning NULL");
1396 rtpg_nmapalgebraexpr_arg_destroy(arg);
1397 PG_RETURN_NULL();
1398 }
1399
1400 /* only work on one or two rasters */
1401 if (arg->bandarg->numraster > 1)
1402 numraster = 2;
1403 else
1404 numraster = 1;
1405
1406 /* pixel type (2) */
1407 if (!PG_ARGISNULL(2)) {
1408 char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
1409
1410 /* Get the pixel type index */
1411 arg->bandarg->pixtype = rt_pixtype_index_from_name(pixtypename);
1412 if (arg->bandarg->pixtype == PT_END) {
1413 rtpg_nmapalgebraexpr_arg_destroy(arg);
1414 elog(ERROR, "RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
1415 PG_RETURN_NULL();
1416 }
1417 }
1418 POSTGIS_RT_DEBUGF(4, "pixeltype: %d", arg->bandarg->pixtype);
1419
1420 /* extent type (3) */
1421 if (!PG_ARGISNULL(3)) {
1422 char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(3))));
1423 arg->bandarg->extenttype = rt_util_extent_type(extenttypename);
1424 }
1425
1426 if (arg->bandarg->extenttype == ET_CUSTOM) {
1427 if (numraster < 2) {
1428 elog(NOTICE, "CUSTOM extent type not supported. Defaulting to FIRST");
1429 arg->bandarg->extenttype = ET_FIRST;
1430 }
1431 else {
1432 elog(NOTICE, "CUSTOM extent type not supported. Defaulting to INTERSECTION");
1433 arg->bandarg->extenttype = ET_INTERSECTION;
1434 }
1435 }
1436 else if (numraster < 2)
1437 arg->bandarg->extenttype = ET_FIRST;
1438
1439 POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->bandarg->extenttype);
1440
1441 /* nodatanodataval (6) */
1442 if (!PG_ARGISNULL(6)) {
1443 arg->callback.nodatanodata.hasval = 1;
1444 arg->callback.nodatanodata.val = PG_GETARG_FLOAT8(6);
1445 }
1446
1447 err = 0;
1448 /* all rasters are empty, return empty raster */
1449 if (allempty == arg->bandarg->numraster) {
1450 elog(NOTICE, "All input rasters are empty. Returning empty raster");
1451 err = 1;
1452 }
1453 /* all rasters don't have indicated band, return empty raster */
1454 else if (noband == arg->bandarg->numraster) {
1455 elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
1456 err = 1;
1457 }
1458 if (err) {
1459 rtpg_nmapalgebraexpr_arg_destroy(arg);
1460
1461 raster = rt_raster_new(0, 0);
1462 if (raster == NULL) {
1463 elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create empty raster");
1464 PG_RETURN_NULL();
1465 }
1466
1467 pgraster = rt_raster_serialize(raster);
1468 rt_raster_destroy(raster);
1469 if (!pgraster) PG_RETURN_NULL();
1470
1471 SET_VARSIZE(pgraster, pgraster->size);
1472 PG_RETURN_POINTER(pgraster);
1473 }
1474
1475 /* connect SPI */
1476 if (SPI_connect() != SPI_OK_CONNECT) {
1477 rtpg_nmapalgebraexpr_arg_destroy(arg);
1478 elog(ERROR, "RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
1479 PG_RETURN_NULL();
1480 }
1481
1482 /*
1483 process expressions
1484
1485 exprpos elements are:
1486 1 - expression => spi_plan[0]
1487 4 - nodata1expr => spi_plan[1]
1488 5 - nodata2expr => spi_plan[2]
1489 */
1490 for (i = 0; i < arg->callback.exprcount; i++) {
1491 char *expr = NULL;
1492 char *tmp = NULL;
1493 char *sql = NULL;
1494 char place[12] = "$1";
1495
1496 if (PG_ARGISNULL(exprpos[i]))
1497 continue;
1498
1499 expr = text_to_cstring(PG_GETARG_TEXT_P(exprpos[i]));
1500 POSTGIS_RT_DEBUGF(3, "raw expr of argument #%d: %s", exprpos[i], expr);
1501
1502 for (j = 0, k = 1; j < argkwcount; j++) {
1503 /* attempt to replace keyword with placeholder */
1504 len = 0;
1505 tmp = rtpg_strreplace(expr, argkw[j], place, &len);
1506 pfree(expr);
1507 expr = tmp;
1508
1509 if (len) {
1510 POSTGIS_RT_DEBUGF(4, "kw #%d (%s) at pos $%d", j, argkw[j], k);
1511 arg->callback.expr[i].spi_argcount++;
1512 arg->callback.expr[i].spi_argpos[j] = k++;
1513
1514 sprintf(place, "$%d", k);
1515 }
1516 else
1517 arg->callback.expr[i].spi_argpos[j] = 0;
1518 }
1519
1520 len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision");
1521 sql = (char *) palloc(len + 1);
1522 if (sql == NULL) {
1523 rtpg_nmapalgebraexpr_arg_destroy(arg);
1524 SPI_finish();
1525 elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
1526 PG_RETURN_NULL();
1527 }
1528
1529 memcpy(sql, "SELECT (", strlen("SELECT ("));
1530 memcpy(sql + strlen("SELECT ("), expr, strlen(expr));
1531 memcpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision"));
1532 sql[len] = '\0';
1533
1534 POSTGIS_RT_DEBUGF(3, "sql #%d: %s", exprpos[i], sql);
1535
1536 /* prepared plan */
1537 if (arg->callback.expr[i].spi_argcount) {
1538 Oid *argtype = (Oid *) palloc(arg->callback.expr[i].spi_argcount * sizeof(Oid));
1539 POSTGIS_RT_DEBUGF(3, "expression parameter %d is a prepared plan", exprpos[i]);
1540 if (argtype == NULL) {
1541 pfree(sql);
1542 rtpg_nmapalgebraexpr_arg_destroy(arg);
1543 SPI_finish();
1544 elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
1545 PG_RETURN_NULL();
1546 }
1547
1548 /* specify datatypes of parameters */
1549 for (j = 0, k = 0; j < argkwcount; j++) {
1550 if (arg->callback.expr[i].spi_argpos[j] < 1) continue;
1551
1552 /* positions are INT4 */
1553 if (
1554 (strstr(argkw[j], "[rast.x]") != NULL) ||
1555 (strstr(argkw[j], "[rast.y]") != NULL) ||
1556 (strstr(argkw[j], "[rast1.x]") != NULL) ||
1557 (strstr(argkw[j], "[rast1.y]") != NULL) ||
1558 (strstr(argkw[j], "[rast2.x]") != NULL) ||
1559 (strstr(argkw[j], "[rast2.y]") != NULL)
1560 )
1561 argtype[k] = INT4OID;
1562 /* everything else is FLOAT8 */
1563 else
1564 argtype[k] = FLOAT8OID;
1565
1566 k++;
1567 }
1568
1569 arg->callback.expr[i].spi_plan = SPI_prepare(sql, arg->callback.expr[i].spi_argcount, argtype);
1570 pfree(argtype);
1571 pfree(sql);
1572
1573 if (arg->callback.expr[i].spi_plan == NULL) {
1574 rtpg_nmapalgebraexpr_arg_destroy(arg);
1575 SPI_finish();
1576 elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
1577 PG_RETURN_NULL();
1578 }
1579 }
1580 /* no args, just execute query */
1581 else {
1582 POSTGIS_RT_DEBUGF(3, "expression parameter %d has no args, simply executing", exprpos[i]);
1583 err = SPI_execute(sql, TRUE, 0);
1584 pfree(sql);
1585
1586 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1587 rtpg_nmapalgebraexpr_arg_destroy(arg);
1588 SPI_finish();
1589 elog(ERROR, "RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
1590 PG_RETURN_NULL();
1591 }
1592
1593 /* get output of prepared plan */
1594 tupdesc = SPI_tuptable->tupdesc;
1595 tuptable = SPI_tuptable;
1596 tuple = tuptable->vals[0];
1597
1598 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1599 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1600 if (SPI_tuptable) SPI_freetuptable(tuptable);
1601 rtpg_nmapalgebraexpr_arg_destroy(arg);
1602 SPI_finish();
1603 elog(ERROR, "RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
1604 PG_RETURN_NULL();
1605 }
1606
1607 if (!isnull) {
1608 arg->callback.expr[i].hasval = 1;
1609 arg->callback.expr[i].val = DatumGetFloat8(datum);
1610 }
1611
1612 if (SPI_tuptable) SPI_freetuptable(tuptable);
1613 }
1614 }
1615
1616 /* determine nodataval and possibly pixtype */
1617 /* band to check */
1618 switch (arg->bandarg->extenttype) {
1619 case ET_LAST:
1620 case ET_SECOND:
1621 if (numraster > 1)
1622 i = 1;
1623 else
1624 i = 0;
1625 break;
1626 default:
1627 i = 0;
1628 break;
1629 }
1630 /* find first viable band */
1631 if (!arg->bandarg->hasband[i]) {
1632 for (i = 0; i < numraster; i++) {
1633 if (arg->bandarg->hasband[i])
1634 break;
1635 }
1636 if (i >= numraster)
1637 i = numraster - 1;
1638 }
1639 band = rt_raster_get_band(arg->bandarg->raster[i], arg->bandarg->nband[i]);
1640
1641 /* set pixel type if PT_END */
1642 if (arg->bandarg->pixtype == PT_END)
1643 arg->bandarg->pixtype = rt_band_get_pixtype(band);
1644
1645 /* set hasnodata and nodataval */
1646 arg->bandarg->hasnodata = 1;
1647 if (rt_band_get_hasnodata_flag(band))
1648 rt_band_get_nodata(band, &(arg->bandarg->nodataval));
1649 else
1650 arg->bandarg->nodataval = rt_band_get_min_value(band);
1651
1652 POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->bandarg->pixtype), arg->bandarg->hasnodata, arg->bandarg->nodataval);
1653
1654 /* init itrset */
1655 itrset = palloc(sizeof(struct rt_iterator_t) * numraster);
1656 if (itrset == NULL) {
1657 rtpg_nmapalgebraexpr_arg_destroy(arg);
1658 SPI_finish();
1659 elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
1660 PG_RETURN_NULL();
1661 }
1662
1663 /* set itrset */
1664 for (i = 0; i < numraster; i++) {
1665 itrset[i].raster = arg->bandarg->raster[i];
1666 itrset[i].nband = arg->bandarg->nband[i];
1667 itrset[i].nbnodata = 1;
1668 }
1669
1670 /* pass everything to iterator */
1671 err = rt_raster_iterator(
1672 itrset, numraster,
1673 arg->bandarg->extenttype, arg->bandarg->cextent,
1674 arg->bandarg->pixtype,
1675 arg->bandarg->hasnodata, arg->bandarg->nodataval,
1676 0, 0,
1677 NULL,
1678 &(arg->callback),
1679 rtpg_nmapalgebraexpr_callback,
1680 &raster
1681 );
1682
1683 pfree(itrset);
1684 rtpg_nmapalgebraexpr_arg_destroy(arg);
1685
1686 if (err != ES_NONE) {
1687 SPI_finish();
1688 elog(ERROR, "RASTER_nMapAlgebraExpr: Could not run raster iterator function");
1689 PG_RETURN_NULL();
1690 }
1691 else if (raster == NULL) {
1692 SPI_finish();
1693 PG_RETURN_NULL();
1694 }
1695
1696 /* switch to prior memory context to ensure memory allocated in correct context */
1697 MemoryContextSwitchTo(mainMemCtx);
1698
1699 pgraster = rt_raster_serialize(raster);
1700 rt_raster_destroy(raster);
1701
1702 /* finish SPI */
1703 SPI_finish();
1704
1705 if (!pgraster)
1706 PG_RETURN_NULL();
1707
1708 SET_VARSIZE(pgraster, pgraster->size);
1709 PG_RETURN_POINTER(pgraster);
1710 }
1711
1712 /* ---------------------------------------------------------------- */
1713 /* ST_Union aggregate functions */
1714 /* ---------------------------------------------------------------- */
1715
1716 typedef enum {
1717 UT_LAST = 0,
1718 UT_FIRST,
1719 UT_MIN,
1720 UT_MAX,
1721 UT_COUNT,
1722 UT_SUM,
1723 UT_MEAN,
1724 UT_RANGE
1725 } rtpg_union_type;
1726
1727 /* internal function translating text of UNION type to enum */
1728 static rtpg_union_type rtpg_uniontype_index_from_name(const char *cutype) {
1729 assert(cutype && strlen(cutype) > 0);
1730
1731 if (strcmp(cutype, "LAST") == 0)
1732 return UT_LAST;
1733 else if (strcmp(cutype, "FIRST") == 0)
1734 return UT_FIRST;
1735 else if (strcmp(cutype, "MIN") == 0)
1736 return UT_MIN;
1737 else if (strcmp(cutype, "MAX") == 0)
1738 return UT_MAX;
1739 else if (strcmp(cutype, "COUNT") == 0)
1740 return UT_COUNT;
1741 else if (strcmp(cutype, "SUM") == 0)
1742 return UT_SUM;
1743 else if (strcmp(cutype, "MEAN") == 0)
1744 return UT_MEAN;
1745 else if (strcmp(cutype, "RANGE") == 0)
1746 return UT_RANGE;
1747
1748 return UT_LAST;
1749 }
1750
1751 typedef struct rtpg_union_band_arg_t *rtpg_union_band_arg;
1752 struct rtpg_union_band_arg_t {
1753 int nband; /* source raster's band index, 0-based */
1754 rtpg_union_type uniontype;
1755
1756 int numraster;
1757 rt_raster *raster;
1758 };
1759
1760 typedef struct rtpg_union_arg_t *rtpg_union_arg;
1761 struct rtpg_union_arg_t {
1762 int numband; /* number of bandargs */
1763 rtpg_union_band_arg bandarg;
1764 };
1765
1766 static void rtpg_union_arg_destroy(rtpg_union_arg arg) {
1767 int i = 0;
1768 int j = 0;
1769 int k = 0;
1770
1771 if (arg->bandarg != NULL) {
1772 for (i = 0; i < arg->numband; i++) {
1773 if (!arg->bandarg[i].numraster)
1774 continue;
1775
1776 for (j = 0; j < arg->bandarg[i].numraster; j++) {
1777 if (arg->bandarg[i].raster[j] == NULL)
1778 continue;
1779
1780 for (k = rt_raster_get_num_bands(arg->bandarg[i].raster[j]) - 1; k >= 0; k--)
1781 rt_band_destroy(rt_raster_get_band(arg->bandarg[i].raster[j], k));
1782 rt_raster_destroy(arg->bandarg[i].raster[j]);
1783 }
1784
1785 pfree(arg->bandarg[i].raster);
1786 }
1787
1788 pfree(arg->bandarg);
1789 }
1790
1791 pfree(arg);
1792 }
1793
1794 static int rtpg_union_callback(
1795 rt_iterator_arg arg, void *userarg,
1796 double *value, int *nodata
1797 ) {
1798 rtpg_union_type *utype = (rtpg_union_type *) userarg;
1799
1800 if (arg == NULL)
1801 return 0;
1802
1803 if (
1804 arg->rasters != 2 ||
1805 arg->rows != 1 ||
1806 arg->columns != 1
1807 ) {
1808 elog(ERROR, "rtpg_union_callback: Invalid arguments passed to callback");
1809 return 0;
1810 }
1811
1812 *value = 0;
1813 *nodata = 0;
1814
1815 /* handle NODATA situations except for COUNT, which is a special case */
1816 if (*utype != UT_COUNT) {
1817 /* both NODATA */
1818 if (arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1819 *nodata = 1;
1820 POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1821 return 1;
1822 }
1823 /* second NODATA */
1824 else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1825 *value = arg->values[0][0][0];
1826 POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1827 return 1;
1828 }
1829 /* first NODATA */
1830 else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0]) {
1831 *value = arg->values[1][0][0];
1832 POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1833 return 1;
1834 }
1835 }
1836
1837 switch (*utype) {
1838 case UT_FIRST:
1839 *value = arg->values[0][0][0];
1840 break;
1841 case UT_MIN:
1842 if (arg->values[0][0][0] < arg->values[1][0][0])
1843 *value = arg->values[0][0][0];
1844 else
1845 *value = arg->values[1][0][0];
1846 break;
1847 case UT_MAX:
1848 if (arg->values[0][0][0] > arg->values[1][0][0])
1849 *value = arg->values[0][0][0];
1850 else
1851 *value = arg->values[1][0][0];
1852 break;
1853 case UT_COUNT:
1854 /* both NODATA */
1855 if (arg->nodata[0][0][0] && arg->nodata[1][0][0])
1856 *value = 0;
1857 /* second NODATA */
1858 else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0])
1859 *value = arg->values[0][0][0];
1860 /* first NODATA */
1861 else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0])
1862 *value = 1;
1863 /* has value, increment */
1864 else
1865 *value = arg->values[0][0][0] + 1;
1866 break;
1867 case UT_SUM:
1868 *value = arg->values[0][0][0] + arg->values[1][0][0];
1869 break;
1870 case UT_MEAN:
1871 case UT_RANGE:
1872 break;
1873 case UT_LAST:
1874 default:
1875 *value = arg->values[1][0][0];
1876 break;
1877 }
1878
1879 POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1880
1881
1882 return 1;
1883 }
1884
1885 static int rtpg_union_mean_callback(
1886 rt_iterator_arg arg, void *userarg,
1887 double *value, int *nodata
1888 ) {
1889 if (arg == NULL)
1890 return 0;
1891
1892 if (
1893 arg->rasters != 2 ||
1894 arg->rows != 1 ||
1895 arg->columns != 1
1896 ) {
1897 elog(ERROR, "rtpg_union_mean_callback: Invalid arguments passed to callback");
1898 return 0;
1899 }
1900
1901 *value = 0;
1902 *nodata = 1;
1903
1904 POSTGIS_RT_DEBUGF(4, "rast0: %f %d", arg->values[0][0][0], arg->nodata[0][0][0]);
1905 POSTGIS_RT_DEBUGF(4, "rast1: %f %d", arg->values[1][0][0], arg->nodata[1][0][0]);
1906
1907 if (!arg->nodata[0][0][0] && FLT_NEQ(arg->values[0][0][0], 0.0) && !arg->nodata[1][0][0])
1908 {
1909 *value = arg->values[1][0][0] / arg->values[0][0][0];
1910 *nodata = 0;
1911 }
1912
1913 POSTGIS_RT_DEBUGF(4, "value, nodata = (%f, %d)", *value, *nodata);
1914
1915 return 1;
1916 }
1917
1918 static int rtpg_union_range_callback(
1919 rt_iterator_arg arg, void *userarg,
1920 double *value, int *nodata
1921 ) {
1922 if (arg == NULL)
1923 return 0;
1924
1925 if (
1926 arg->rasters != 2 ||
1927 arg->rows != 1 ||
1928 arg->columns != 1
1929 ) {
1930 elog(ERROR, "rtpg_union_range_callback: Invalid arguments passed to callback");
1931 return 0;
1932 }
1933
1934 *value = 0;
1935 *nodata = 1;
1936
1937 POSTGIS_RT_DEBUGF(4, "rast0: %f %d", arg->values[0][0][0], arg->nodata[0][0][0]);
1938 POSTGIS_RT_DEBUGF(4, "rast1: %f %d", arg->values[1][0][0], arg->nodata[1][0][0]);
1939
1940 if (
1941 !arg->nodata[0][0][0] &&
1942 !arg->nodata[1][0][0]
1943 ) {
1944 *value = arg->values[1][0][0] - arg->values[0][0][0];
1945 *nodata = 0;
1946 }
1947
1948 POSTGIS_RT_DEBUGF(4, "value, nodata = (%f, %d)", *value, *nodata);
1949
1950 return 1;
1951 }
1952
1953 /* called for ST_Union(raster, unionarg[]) */
1954 static int rtpg_union_unionarg_process(rtpg_union_arg arg, ArrayType *array) {
1955 Oid etype;
1956 Datum *e;
1957 bool *nulls;
1958 int16 typlen;
1959 bool typbyval;
1960 char typalign;
1961 int n = 0;
1962
1963 HeapTupleHeader tup;
1964 bool isnull;
1965 Datum tupv;
1966
1967 int i;
1968 int nband = 1;
1969 char *utypename = NULL;
1970 rtpg_union_type utype = UT_LAST;
1971
1972 etype = ARR_ELEMTYPE(array);
1973 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1974
1975 deconstruct_array(
1976 array,
1977 etype,
1978 typlen, typbyval, typalign,
1979 &e, &nulls, &n
1980 );
1981
1982 if (!n) {
1983 elog(ERROR, "rtpg_union_unionarg_process: Invalid argument for unionarg");
1984 return 0;
1985 }
1986
1987 /* prep arg */
1988 arg->numband = n;
1989 arg->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * arg->numband);
1990 if (arg->bandarg == NULL) {
1991 elog(ERROR, "rtpg_union_unionarg_process: Could not allocate memory for band information");
1992 return 0;
1993 }
1994
1995 /* process each element */
1996 for (i = 0; i < n; i++) {
1997 if (nulls[i]) {
1998 arg->numband--;
1999 continue;
2000 }
2001
2002 POSTGIS_RT_DEBUGF(4, "Processing unionarg at index %d", i);
2003
2004 /* each element is a tuple */
2005 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
2006 if (NULL == tup) {
2007 elog(ERROR, "rtpg_union_unionarg_process: Invalid argument for unionarg");
2008 return 0;
2009 }
2010
2011 /* first element, bandnum */
2012 tupv = GetAttributeByName(tup, "nband", &isnull);
2013 if (isnull) {
2014 nband = i + 1;
2015 elog(NOTICE, "First argument (nband) of unionarg is NULL. Assuming nband = %d", nband);
2016 }
2017 else
2018 nband = DatumGetInt32(tupv);
2019
2020 if (nband < 1) {
2021 elog(ERROR, "rtpg_union_unionarg_process: Band number must be greater than zero (1-based)");
2022 return 0;
2023 }
2024
2025 /* second element, uniontype */
2026 tupv = GetAttributeByName(tup, "uniontype", &isnull);
2027 if (isnull) {
2028 elog(NOTICE, "Second argument (uniontype) of unionarg is NULL. Assuming uniontype = LAST");
2029 utype = UT_LAST;
2030 }
2031 else {
2032 utypename = text_to_cstring((text *) DatumGetPointer(tupv));
2033 utype = rtpg_uniontype_index_from_name(rtpg_strtoupper(utypename));
2034 }
2035
2036 arg->bandarg[i].uniontype = utype;
2037 arg->bandarg[i].nband = nband - 1;
2038 arg->bandarg[i].raster = NULL;
2039
2040 if (
2041 utype != UT_MEAN &&
2042 utype != UT_RANGE
2043 ) {
2044 arg->bandarg[i].numraster = 1;
2045 }
2046 else
2047 arg->bandarg[i].numraster = 2;
2048 }
2049
2050 if (arg->numband < n) {
2051 arg->bandarg = repalloc(arg->bandarg, sizeof(struct rtpg_union_band_arg_t) * arg->numband);
2052 if (arg->bandarg == NULL) {
2053 elog(ERROR, "rtpg_union_unionarg_process: Could not reallocate memory for band information");
2054 return 0;
2055 }
2056 }
2057
2058 return 1;
2059 }
2060
2061 /* called for ST_Union(raster) */
2062 static int rtpg_union_noarg(rtpg_union_arg arg, rt_raster raster) {
2063 int numbands;
2064 int i;
2065
2066 if (rt_raster_is_empty(raster))
2067 return 1;
2068
2069 numbands = rt_raster_get_num_bands(raster);
2070 if (numbands <= arg->numband)
2071 return 1;
2072
2073 /* more bands to process */
2074 POSTGIS_RT_DEBUG(4, "input raster has more bands, adding more bandargs");
2075 if (arg->numband)
2076 arg->bandarg = repalloc(arg->bandarg, sizeof(struct rtpg_union_band_arg_t) * numbands);
2077 else
2078 arg->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * numbands);
2079 if (arg->bandarg == NULL) {
2080 elog(ERROR, "rtpg_union_noarg: Could not reallocate memory for band information");
2081 return 0;
2082 }
2083
2084 i = arg->numband;
2085 arg->numband = numbands;
2086 for (; i < arg->numband; i++) {
2087 POSTGIS_RT_DEBUGF(4, "Adding bandarg for band at index %d", i);
2088 arg->bandarg[i].uniontype = UT_LAST;
2089 arg->bandarg[i].nband = i;
2090 arg->bandarg[i].numraster = 1;
2091
2092 arg->bandarg[i].raster = (rt_raster *) palloc(sizeof(rt_raster) * arg->bandarg[i].numraster);
2093 if (arg->bandarg[i].raster == NULL) {
2094 elog(ERROR, "rtpg_union_noarg: Could not allocate memory for working rasters");
2095 return 0;
2096 }
2097 memset(arg->bandarg[i].raster, 0, sizeof(rt_raster) * arg->bandarg[i].numraster);
2098
2099 /* add new working rt_raster but only if working raster already exists */
2100 if (!rt_raster_is_empty(arg->bandarg[0].raster[0])) {
2101 arg->bandarg[i].raster[0] = rt_raster_clone(arg->bandarg[0].raster[0], 0); /* shallow clone */
2102 if (arg->bandarg[i].raster[0] == NULL) {
2103 elog(ERROR, "rtpg_union_noarg: Could not create working raster");
2104 return 0;
2105 }
2106 }
2107 }
2108
2109 return 1;
2110 }
2111
2112 /* UNION aggregate transition function */
2113 PG_FUNCTION_INFO_V1(RASTER_union_transfn);
2114 Datum RASTER_union_transfn(PG_FUNCTION_ARGS)
2115 {
2116 MemoryContext aggcontext;
2117 MemoryContext oldcontext;
2118 rtpg_union_arg iwr = NULL;
2119 int skiparg = 0;
2120
2121 rt_pgraster *pgraster = NULL;
2122 rt_raster raster = NULL;
2123 rt_raster _raster = NULL;
2124 rt_band _band = NULL;
2125 int nband = 1;
2126 int noerr = 1;
2127 int isempty[2] = {0};
2128 int hasband[2] = {0};
2129 int nargs = 0;
2130 double _offset[4] = {0.};
2131 int nbnodata = 0; /* 1 if adding bands */
2132
2133 int i = 0;
2134 int j = 0;
2135 int k = 0;
2136
2137 rt_iterator itrset;
2138 char *utypename = NULL;
2139 rtpg_union_type utype = UT_LAST;
2140 rt_pixtype pixtype = PT_END;
2141 int hasnodata = 1;
2142 double nodataval = 0;
2143
2144 rt_raster iraster = NULL;
2145 rt_band iband = NULL;
2146 int reuserast = 0;
2147 int y = 0;
2148 uint16_t _dim[2] = {0};
2149 void *vals = NULL;
2150 uint16_t nvals = 0;
2151
2152 POSTGIS_RT_DEBUG(3, "Starting...");
2153
2154 /* cannot be called directly as this is exclusive aggregate function */
2155 if (!AggCheckCallContext(fcinfo, &aggcontext)) {
2156 elog(ERROR, "RASTER_union_transfn: Cannot be called in a non-aggregate context");
2157 PG_RETURN_NULL();
2158 }
2159
2160 /* switch to aggcontext */
2161 oldcontext = MemoryContextSwitchTo(aggcontext);
2162
2163 if (PG_ARGISNULL(0)) {
2164 POSTGIS_RT_DEBUG(3, "Creating state variable");
2165 /* allocate container in aggcontext */
2166 iwr = MemoryContextAlloc(aggcontext, sizeof(struct rtpg_union_arg_t));
2167 if (iwr == NULL) {
2168 MemoryContextSwitchTo(oldcontext);
2169 elog(ERROR, "RASTER_union_transfn: Could not allocate memory for state variable");
2170 PG_RETURN_NULL();
2171 }
2172
2173 iwr->numband = 0;
2174 iwr->bandarg = NULL;
2175
2176 skiparg = 0;
2177 }
2178 else {
2179 POSTGIS_RT_DEBUG(3, "State variable already exists");
2180 iwr = (rtpg_union_arg) PG_GETARG_POINTER(0);
2181 skiparg = 1;
2182 }
2183
2184 /* raster arg is NOT NULL */
2185 if (!PG_ARGISNULL(1)) {
2186 /* deserialize raster */
2187 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2188
2189 /* Get raster object */
2190 raster = rt_raster_deserialize(pgraster, FALSE);
2191 if (raster == NULL) {
2192
2193 rtpg_union_arg_destroy(iwr);
2194 PG_FREE_IF_COPY(pgraster, 1);
2195
2196 MemoryContextSwitchTo(oldcontext);
2197 elog(ERROR, "RASTER_union_transfn: Could not deserialize raster");
2198 PG_RETURN_NULL();
2199 }
2200 }
2201
2202 /* process additional args if needed */
2203 nargs = PG_NARGS();
2204 POSTGIS_RT_DEBUGF(4, "nargs = %d", nargs);
2205 if (nargs > 2) {
2206 POSTGIS_RT_DEBUG(4, "processing additional arguments");
2207
2208 /* if more than 2 arguments, determine the type of argument 3 */
2209 /* band number, UNION type or unionarg */
2210 if (!PG_ARGISNULL(2)) {
2211 Oid calltype = get_fn_expr_argtype(fcinfo->flinfo, 2);
2212
2213 switch (calltype) {
2214 /* UNION type */
2215 case TEXTOID: {
2216 int idx = 0;
2217 int numband = 0;
2218
2219 POSTGIS_RT_DEBUG(4, "Processing arg 3 as UNION type");
2220 nbnodata = 1;
2221
2222 utypename = text_to_cstring(PG_GETARG_TEXT_P(2));
2223 utype = rtpg_uniontype_index_from_name(rtpg_strtoupper(utypename));
2224 POSTGIS_RT_DEBUGF(4, "Union type: %s", utypename);
2225
2226 POSTGIS_RT_DEBUGF(4, "iwr->numband: %d", iwr->numband);
2227 /* see if we need to append new bands */
2228 if (raster) {
2229 idx = iwr->numband;
2230 numband = rt_raster_get_num_bands(raster);
2231 POSTGIS_RT_DEBUGF(4, "numband: %d", numband);
2232
2233 /* only worry about appended bands */
2234 if (numband > iwr->numband)
2235 iwr->numband = numband;
2236 }
2237
2238 if (!iwr->numband)
2239 iwr->numband = 1;
2240 POSTGIS_RT_DEBUGF(4, "iwr->numband: %d", iwr->numband);
2241 POSTGIS_RT_DEBUGF(4, "numband, idx: %d, %d", numband, idx);
2242
2243 /* bandarg set. only possible after the first call to function */
2244 if (iwr->bandarg) {
2245 /* only reallocate if new bands need to be added */
2246 if (numband > idx) {
2247 POSTGIS_RT_DEBUG(4, "Reallocating iwr->bandarg");
2248 iwr->bandarg = repalloc(iwr->bandarg, sizeof(struct rtpg_union_band_arg_t) * iwr->numband);
2249 }
2250 /* prevent initial values step happening below */
2251 else
2252 idx = iwr->numband;
2253 }
2254 /* bandarg not set, first call to function */
2255 else {
2256 POSTGIS_RT_DEBUG(4, "Allocating iwr->bandarg");
2257 iwr->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * iwr->numband);
2258 }
2259 if (iwr->bandarg == NULL) {
2260
2261 rtpg_union_arg_destroy(iwr);
2262 if (raster != NULL) {
2263 rt_raster_destroy(raster);
2264 PG_FREE_IF_COPY(pgraster, 1);
2265 }
2266
2267 MemoryContextSwitchTo(oldcontext);
2268 elog(ERROR, "RASTER_union_transfn: Could not allocate memory for band information");
2269 PG_RETURN_NULL();
2270 }
2271
2272 /* set initial values for bands that are "new" */
2273 for (i = idx; i < iwr->numband; i++) {
2274 iwr->bandarg[i].uniontype = utype;
2275 iwr->bandarg[i].nband = i;
2276
2277 if (
2278 utype == UT_MEAN ||
2279 utype == UT_RANGE
2280 ) {
2281 iwr->bandarg[i].numraster = 2;
2282 }
2283 else
2284 iwr->bandarg[i].numraster = 1;
2285 iwr->bandarg[i].raster = NULL;
2286 }
2287
2288 break;
2289 }
2290 /* band number */
2291 case INT2OID:
2292 case INT4OID:
2293 if (skiparg)
2294 break;
2295
2296 POSTGIS_RT_DEBUG(4, "Processing arg 3 as band number");
2297 nband = PG_GETARG_INT32(2);
2298 if (nband < 1) {
2299
2300 rtpg_union_arg_destroy(iwr);
2301 if (raster != NULL) {
2302 rt_raster_destroy(raster);
2303 PG_FREE_IF_COPY(pgraster, 1);
2304 }
2305
2306 MemoryContextSwitchTo(oldcontext);
2307 elog(ERROR, "RASTER_union_transfn: Band number must be greater than zero (1-based)");
2308 PG_RETURN_NULL();
2309 }
2310
2311 iwr->numband = 1;
2312 iwr->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * iwr->numband);
2313 if (iwr->bandarg == NULL) {
2314
2315 rtpg_union_arg_destroy(iwr);
2316 if (raster != NULL) {
2317 rt_raster_destroy(raster);
2318 PG_FREE_IF_COPY(pgraster, 1);
2319 }
2320
2321 MemoryContextSwitchTo(oldcontext);
2322 elog(ERROR, "RASTER_union_transfn: Could not allocate memory for band information");
2323 PG_RETURN_NULL();
2324 }
2325
2326 iwr->bandarg[0].uniontype = UT_LAST;
2327 iwr->bandarg[0].nband = nband - 1;
2328
2329 iwr->bandarg[0].numraster = 1;
2330 iwr->bandarg[0].raster = NULL;
2331 break;
2332 /* only other type allowed is unionarg */
2333 default:
2334 if (skiparg)
2335 break;
2336
2337 POSTGIS_RT_DEBUG(4, "Processing arg 3 as unionarg");
2338 if (!rtpg_union_unionarg_process(iwr, PG_GETARG_ARRAYTYPE_P(2))) {
2339
2340 rtpg_union_arg_destroy(iwr);
2341 if (raster != NULL) {
2342 rt_raster_destroy(raster);
2343 PG_FREE_IF_COPY(pgraster, 1);
2344 }
2345
2346 MemoryContextSwitchTo(oldcontext);
2347 elog(ERROR, "RASTER_union_transfn: Could not process unionarg");
2348 PG_RETURN_NULL();
2349 }
2350
2351 break;
2352 }
2353 }
2354
2355 /* UNION type */
2356 if (nargs > 3 && !PG_ARGISNULL(3)) {
2357 utypename = text_to_cstring(PG_GETARG_TEXT_P(3));
2358 utype = rtpg_uniontype_index_from_name(rtpg_strtoupper(utypename));
2359 iwr->bandarg[0].uniontype = utype;
2360 POSTGIS_RT_DEBUGF(4, "Union type: %s", utypename);
2361
2362 if (
2363 utype == UT_MEAN ||
2364 utype == UT_RANGE
2365 ) {
2366 iwr->bandarg[0].numraster = 2;
2367 }
2368 }
2369
2370 /* allocate space for pointers to rt_raster */
2371 for (i = 0; i < iwr->numband; i++) {
2372 POSTGIS_RT_DEBUGF(4, "iwr->bandarg[%d].raster @ %p", i, iwr->bandarg[i].raster);
2373
2374 /* no need to allocate */
2375 if (iwr->bandarg[i].raster != NULL)
2376 continue;
2377
2378 POSTGIS_RT_DEBUGF(4, "Allocating space for working rasters of band %d", i);
2379
2380 iwr->bandarg[i].raster = (rt_raster *) palloc(sizeof(rt_raster) * iwr->bandarg[i].numraster);
2381 if (iwr->bandarg[i].raster == NULL) {
2382
2383 rtpg_union_arg_destroy(iwr);
2384 if (raster != NULL) {
2385 rt_raster_destroy(raster);
2386 PG_FREE_IF_COPY(pgraster, 1);
2387 }
2388
2389 MemoryContextSwitchTo(oldcontext);
2390 elog(ERROR, "RASTER_union_transfn: Could not allocate memory for working raster(s)");
2391 PG_RETURN_NULL();
2392 }
2393
2394 memset(iwr->bandarg[i].raster, 0, sizeof(rt_raster) * iwr->bandarg[i].numraster);
2395
2396 /* add new working rt_raster but only if working raster already exists */
2397 if (i > 0 && !rt_raster_is_empty(iwr->bandarg[0].raster[0])) {
2398 for (j = 0; j < iwr->bandarg[i].numraster; j++) {
2399 iwr->bandarg[i].raster[j] = rt_raster_clone(iwr->bandarg[0].raster[0], 0); /* shallow clone */
2400 if (iwr->bandarg[i].raster[j] == NULL) {
2401
2402 rtpg_union_arg_destroy(iwr);
2403 if (raster != NULL) {
2404 rt_raster_destroy(raster);
2405 PG_FREE_IF_COPY(pgraster, 1);
2406 }
2407
2408 MemoryContextSwitchTo(oldcontext);
2409 elog(ERROR, "RASTER_union_transfn: Could not create working raster");
2410 PG_RETURN_NULL();
2411 }
2412 }
2413 }
2414 }
2415 }
2416 /* only raster, no additional args */
2417 /* only do this if raster isn't empty */
2418 else {
2419 POSTGIS_RT_DEBUG(4, "no additional args, checking input raster");
2420 nbnodata = 1;
2421 if (!rtpg_union_noarg(iwr, raster)) {
2422
2423 rtpg_union_arg_destroy(iwr);
2424 if (raster != NULL) {
2425 rt_raster_destroy(raster);
2426 PG_FREE_IF_COPY(pgraster, 1);
2427 }
2428
2429 MemoryContextSwitchTo(oldcontext);
2430 elog(ERROR, "RASTER_union_transfn: Could not check and balance number of bands");
2431 PG_RETURN_NULL();
2432 }
2433 }
2434
2435 /* init itrset */
2436 itrset = palloc(sizeof(struct rt_iterator_t) * 2);
2437 if (itrset == NULL) {
2438
2439 rtpg_union_arg_destroy(iwr);
2440 if (raster != NULL) {
2441 rt_raster_destroy(raster);
2442 PG_FREE_IF_COPY(pgraster, 1);
2443 }
2444
2445 MemoryContextSwitchTo(oldcontext);
2446 elog(ERROR, "RASTER_union_transfn: Could not allocate memory for iterator arguments");
2447 PG_RETURN_NULL();
2448 }
2449
2450 /* by band to UNION */
2451 for (i = 0; i < iwr->numband; i++) {
2452
2453 /* by raster */
2454 for (j = 0; j < iwr->bandarg[i].numraster; j++) {
2455 reuserast = 0;
2456
2457 /* type of union */
2458 utype = iwr->bandarg[i].uniontype;
2459
2460 /* raster flags */
2461 isempty[0] = rt_raster_is_empty(iwr->bandarg[i].raster[j]);
2462 isempty[1] = rt_raster_is_empty(raster);
2463
2464 if (!isempty[0])
2465 hasband[0] = rt_raster_has_band(iwr->bandarg[i].raster[j], 0);
2466 if (!isempty[1])
2467 hasband[1] = rt_raster_has_band(raster, iwr->bandarg[i].nband);
2468
2469 /* determine pixtype, hasnodata and nodataval */
2470 _band = NULL;
2471 if (!isempty[0] && hasband[0])
2472 _band = rt_raster_get_band(iwr->bandarg[i].raster[j], 0);
2473 else if (!isempty[1] && hasband[1])
2474 _band = rt_raster_get_band(raster, iwr->bandarg[i].nband);
2475 else {
2476 pixtype = PT_64BF;
2477 hasnodata = 1;
2478 nodataval = rt_pixtype_get_min_value(pixtype);
2479 }
2480 if (_band != NULL) {
2481 pixtype = rt_band_get_pixtype(_band);
2482 hasnodata = 1;
2483 if (rt_band_get_hasnodata_flag(_band))
2484 rt_band_get_nodata(_band, &nodataval);
2485 else
2486 nodataval = rt_band_get_min_value(_band);
2487 }
2488
2489 /* UT_MEAN and UT_RANGE require two passes */
2490 /* UT_MEAN: first for UT_COUNT and second for UT_SUM */
2491 if (iwr->bandarg[i].uniontype == UT_MEAN) {
2492 /* first pass, UT_COUNT */
2493 if (j < 1)
2494 utype = UT_COUNT;
2495 else
2496 utype = UT_SUM;
2497 }
2498 /* UT_RANGE: first for UT_MIN and second for UT_MAX */
2499 else if (iwr->bandarg[i].uniontype == UT_RANGE) {
2500 /* first pass, UT_MIN */
2501 if (j < 1)
2502 utype = UT_MIN;
2503 else
2504 utype = UT_MAX;
2505 }
2506
2507 /* force band settings for UT_COUNT */
2508 if (utype == UT_COUNT) {
2509 pixtype = PT_32BUI;
2510 hasnodata = 0;
2511 nodataval = 0;
2512 }
2513
2514 POSTGIS_RT_DEBUGF(4, "(pixtype, hasnodata, nodataval) = (%s, %d, %f)", rt_pixtype_name(pixtype), hasnodata, nodataval);
2515
2516 /* set itrset */
2517 itrset[0].raster = iwr->bandarg[i].raster[j];
2518 itrset[0].nband = 0;
2519 itrset[1].raster = raster;
2520 itrset[1].nband = iwr->bandarg[i].nband;
2521
2522 /* allow use NODATA to replace missing bands */
2523 if (nbnodata) {
2524 itrset[0].nbnodata = 1;
2525 itrset[1].nbnodata = 1;
2526 }
2527 /* missing bands are not ignored */
2528 else {
2529 itrset[0].nbnodata = 0;
2530 itrset[1].nbnodata = 0;
2531 }
2532
2533 /* if rasters AND bands are present, use copy approach */
2534 if (!isempty[0] && !isempty[1] && hasband[0] && hasband[1]) {
2535 POSTGIS_RT_DEBUG(3, "using line method");
2536
2537 /* generate empty out raster */
2538 if (rt_raster_from_two_rasters(
2539 iwr->bandarg[i].raster[j], raster,
2540 ET_UNION,
2541 &iraster, _offset
2542 ) != ES_NONE) {
2543
2544 pfree(itrset);
2545 rtpg_union_arg_destroy(iwr);
2546 if (raster != NULL) {
2547 rt_raster_destroy(raster);
2548 PG_FREE_IF_COPY(pgraster, 1);
2549 }
2550
2551 MemoryContextSwitchTo(oldcontext);
2552 elog(ERROR, "RASTER_union_transfn: Could not create internal raster");
2553 PG_RETURN_NULL();
2554 }
2555 POSTGIS_RT_DEBUGF(4, "_offset = %f, %f, %f, %f",
2556 _offset[0], _offset[1], _offset[2], _offset[3]);
2557
2558 /* rasters are spatially the same? */
2559 if (
2560 rt_raster_get_width(iwr->bandarg[i].raster[j]) == rt_raster_get_width(iraster) &&
2561 rt_raster_get_height(iwr->bandarg[i].raster[j]) == rt_raster_get_height(iraster)
2562 ) {
2563 double igt[6] = {0};
2564 double gt[6] = {0};
2565
2566 rt_raster_get_geotransform_matrix(iwr->bandarg[i].raster[j], gt);
2567 rt_raster_get_geotransform_matrix(iraster, igt);
2568
2569 reuserast = rt_util_same_geotransform_matrix(gt, igt);
2570 }
2571
2572 /* use internal raster */
2573 if (!reuserast) {
2574 /* create band of same type */
2575 if (rt_raster_generate_new_band(
2576 iraster,
2577 pixtype,
2578 nodataval,
2579 hasnodata, nodataval,
2580 0
2581 ) == -1) {
2582
2583 pfree(itrset);
2584 rtpg_union_arg_destroy(iwr);
2585 rt_raster_destroy(iraster);
2586 if (raster != NULL) {
2587 rt_raster_destroy(raster);
2588 PG_FREE_IF_COPY(pgraster, 1);
2589 }
2590
2591 MemoryContextSwitchTo(oldcontext);
2592 elog(ERROR, "RASTER_union_transfn: Could not add new band to internal raster");
2593 PG_RETURN_NULL();
2594 }
2595 iband = rt_raster_get_band(iraster, 0);
2596
2597 /* copy working raster to output raster */
2598 _dim[0] = rt_raster_get_width(iwr->bandarg[i].raster[j]);
2599 _dim[1] = rt_raster_get_height(iwr->bandarg[i].raster[j]);
2600 for (y = 0; y < _dim[1]; y++) {
2601 POSTGIS_RT_DEBUGF(4, "Getting pixel line of working raster at (x, y, length) = (0, %d, %d)", y, _dim[0]);
2602 if (rt_band_get_pixel_line(
2603 _band,
2604 0, y,
2605 _dim[0],
2606 &vals, &nvals
2607 ) != ES_NONE) {
2608
2609 pfree(itrset);
2610 rtpg_union_arg_destroy(iwr);
2611 rt_band_destroy(iband);
2612 rt_raster_destroy(iraster);
2613 if (raster != NULL) {
2614 rt_raster_destroy(raster);
2615 PG_FREE_IF_COPY(pgraster, 1);
2616 }
2617
2618 MemoryContextSwitchTo(oldcontext);
2619 elog(ERROR, "RASTER_union_transfn: Could not get pixel line from band of working raster");
2620 PG_RETURN_NULL();
2621 }
2622
2623 POSTGIS_RT_DEBUGF(4, "Setting pixel line at (x, y, length) = (%d, %d, %d)", (int) _offset[0], (int) _offset[1] + y, nvals);
2624 if (rt_band_set_pixel_line(
2625 iband,
2626 (int) _offset[0], (int) _offset[1] + y,
2627 vals, nvals
2628 ) != ES_NONE) {
2629
2630 pfree(itrset);
2631 rtpg_union_arg_destroy(iwr);
2632 rt_band_destroy(iband);
2633 rt_raster_destroy(iraster);
2634 if (raster != NULL) {
2635 rt_raster_destroy(raster);
2636 PG_FREE_IF_COPY(pgraster, 1);
2637 }
2638
2639 MemoryContextSwitchTo(oldcontext);
2640 elog(ERROR, "RASTER_union_transfn: Could not set pixel line to band of internal raster");
2641 PG_RETURN_NULL();
2642 }
2643 }
2644 }
2645 else {
2646 rt_raster_destroy(iraster);
2647 iraster = iwr->bandarg[i].raster[j];
2648 iband = rt_raster_get_band(iraster, 0);
2649 }
2650
2651 /* run iterator for extent of input raster */
2652 noerr = rt_raster_iterator(
2653 itrset, 2,
2654 ET_LAST, NULL,
2655 pixtype,
2656 hasnodata, nodataval,
2657 0, 0,
2658 NULL,
2659 &utype,
2660 rtpg_union_callback,
2661 &_raster
2662 );
2663 if (noerr != ES_NONE) {
2664
2665 pfree(itrset);
2666 rtpg_union_arg_destroy(iwr);
2667 if (!reuserast) {
2668 rt_band_destroy(iband);
2669 rt_raster_destroy(iraster);
2670 }
2671 if (raster != NULL) {
2672 rt_raster_destroy(raster);
2673 PG_FREE_IF_COPY(pgraster, 1);
2674 }
2675
2676 MemoryContextSwitchTo(oldcontext);
2677 elog(ERROR, "RASTER_union_transfn: Could not run raster iterator function");
2678 PG_RETURN_NULL();
2679 }
2680
2681 /* with iterator raster, copy data to output raster */
2682 _band = rt_raster_get_band(_raster, 0);
2683 _dim[0] = rt_raster_get_width(_raster);
2684 _dim[1] = rt_raster_get_height(_raster);
2685 for (y = 0; y < _dim[1]; y++) {
2686 POSTGIS_RT_DEBUGF(4, "Getting pixel line of iterator raster at (x, y, length) = (0, %d, %d)", y, _dim[0]);
2687 if (rt_band_get_pixel_line(
2688 _band,
2689 0, y,
2690 _dim[0],
2691 &vals, &nvals
2692 ) != ES_NONE) {
2693
2694 pfree(itrset);
2695 rtpg_union_arg_destroy(iwr);
2696 if (!reuserast) {
2697 rt_band_destroy(iband);
2698 rt_raster_destroy(iraster);
2699 }
2700 if (raster != NULL) {
2701 rt_raster_destroy(raster);
2702 PG_FREE_IF_COPY(pgraster, 1);
2703 }
2704
2705 MemoryContextSwitchTo(oldcontext);
2706 elog(ERROR, "RASTER_union_transfn: Could not get pixel line from band of working raster");
2707 PG_RETURN_NULL();
2708 }
2709
2710 POSTGIS_RT_DEBUGF(4, "Setting pixel line at (x, y, length) = (%d, %d, %d)", (int) _offset[2], (int) _offset[3] + y, nvals);
2711 if (rt_band_set_pixel_line(
2712 iband,
2713 (int) _offset[2], (int) _offset[3] + y,
2714 vals, nvals
2715 ) != ES_NONE) {
2716
2717 pfree(itrset);
2718 rtpg_union_arg_destroy(iwr);
2719 if (!reuserast) {
2720 rt_band_destroy(iband);
2721 rt_raster_destroy(iraster);
2722 }
2723 if (raster != NULL) {
2724 rt_raster_destroy(raster);
2725 PG_FREE_IF_COPY(pgraster, 1);
2726 }
2727
2728 MemoryContextSwitchTo(oldcontext);
2729 elog(ERROR, "RASTER_union_transfn: Could not set pixel line to band of internal raster");
2730 PG_RETURN_NULL();
2731 }
2732 }
2733
2734 /* free _raster */
2735 rt_band_destroy(_band);
2736 rt_raster_destroy(_raster);
2737
2738 /* replace working raster with output raster */
2739 _raster = iraster;
2740 }
2741 else {
2742 POSTGIS_RT_DEBUG(3, "using pixel method");
2743
2744 /* pass everything to iterator */
2745 noerr = rt_raster_iterator(
2746 itrset, 2,
2747 ET_UNION, NULL,
2748 pixtype,
2749 hasnodata, nodataval,
2750 0, 0,
2751 NULL,
2752 &utype,
2753 rtpg_union_callback,
2754 &_raster
2755 );
2756
2757 if (noerr != ES_NONE) {
2758
2759 pfree(itrset);
2760 rtpg_union_arg_destroy(iwr);
2761 if (raster != NULL) {
2762 rt_raster_destroy(raster);
2763 PG_FREE_IF_COPY(pgraster, 1);
2764 }
2765
2766 MemoryContextSwitchTo(oldcontext);
2767 elog(ERROR, "RASTER_union_transfn: Could not run raster iterator function");
2768 PG_RETURN_NULL();
2769 }
2770 }
2771
2772 /* replace working raster */
2773 if (iwr->bandarg[i].raster[j] != NULL && !reuserast) {
2774 for (k = rt_raster_get_num_bands(iwr->bandarg[i].raster[j]) - 1; k >= 0; k--)
2775 rt_band_destroy(rt_raster_get_band(iwr->bandarg[i].raster[j], k));
2776 rt_raster_destroy(iwr->bandarg[i].raster[j]);
2777 }
2778 iwr->bandarg[i].raster[j] = _raster;
2779 }
2780
2781 }
2782
2783 pfree(itrset);
2784 if (raster != NULL) {
2785 rt_raster_destroy(raster);
2786 PG_FREE_IF_COPY(pgraster, 1);
2787 }
2788
2789 /* switch back to local context */
2790 MemoryContextSwitchTo(oldcontext);
2791
2792 POSTGIS_RT_DEBUG(3, "Finished");
2793
2794 PG_RETURN_POINTER(iwr);
2795 }
2796
2797 /* UNION aggregate final function */
2798 PG_FUNCTION_INFO_V1(RASTER_union_finalfn);
2799 Datum RASTER_union_finalfn(PG_FUNCTION_ARGS)
2800 {
2801 rtpg_union_arg iwr;
2802 rt_raster _rtn = NULL;
2803 rt_raster _raster = NULL;
2804 rt_pgraster *pgraster = NULL;
2805
2806 int i = 0;
2807 int j = 0;
2808 rt_iterator itrset = NULL;
2809 rt_band _band = NULL;
2810 int noerr = 1;
2811 int status = 0;
2812 rt_pixtype pixtype = PT_END;
2813 int hasnodata = 0;
2814 double nodataval = 0;
2815
2816 POSTGIS_RT_DEBUG(3, "Starting...");
2817
2818 /* cannot be called directly as this is exclusive aggregate function */
2819 if (!AggCheckCallContext(fcinfo, NULL)) {
2820 elog(ERROR, "RASTER_union_finalfn: Cannot be called in a non-aggregate context");
2821 PG_RETURN_NULL();
2822 }
2823
2824 /* NULL, return null */
2825 if (PG_ARGISNULL(0))
2826 PG_RETURN_NULL();
2827
2828 iwr = (rtpg_union_arg) PG_GETARG_POINTER(0);
2829
2830 /* init itrset */
2831 itrset = palloc(sizeof(struct rt_iterator_t) * 2);
2832 if (itrset == NULL) {
2833 rtpg_union_arg_destroy(iwr);
2834 elog(ERROR, "RASTER_union_finalfn: Could not allocate memory for iterator arguments");
2835 PG_RETURN_NULL();
2836 }
2837
2838 for (i = 0; i < iwr->numband; i++) {
2839 if (
2840 iwr->bandarg[i].uniontype == UT_MEAN ||
2841 iwr->bandarg[i].uniontype == UT_RANGE
2842 ) {
2843 /* raster containing the SUM or MAX is at index 1 */
2844 _band = rt_raster_get_band(iwr->bandarg[i].raster[1], 0);
2845
2846 pixtype = rt_band_get_pixtype(_band);
2847 hasnodata = rt_band_get_hasnodata_flag(_band);
2848 if (hasnodata)
2849 rt_band_get_nodata(_band, &nodataval);
2850 POSTGIS_RT_DEBUGF(4, "(pixtype, hasnodata, nodataval) = (%s, %d, %f)", rt_pixtype_name(pixtype), hasnodata, nodataval);
2851
2852 itrset[0].raster = iwr->bandarg[i].raster[0];
2853 itrset[0].nband = 0;
2854 itrset[1].raster = iwr->bandarg[i].raster[1];
2855 itrset[1].nband = 0;
2856
2857 /* pass everything to iterator */
2858 if (iwr->bandarg[i].uniontype == UT_MEAN) {
2859 noerr = rt_raster_iterator(
2860 itrset, 2,
2861 ET_UNION, NULL,
2862 pixtype,
2863 hasnodata, nodataval,
2864 0, 0,
2865 NULL,
2866 NULL,
2867 rtpg_union_mean_callback,
2868 &_raster
2869 );
2870 }
2871 else if (iwr->bandarg[i].uniontype == UT_RANGE) {
2872 noerr = rt_raster_iterator(
2873 itrset, 2,
2874 ET_UNION, NULL,
2875 pixtype,
2876 hasnodata, nodataval,
2877 0, 0,
2878 NULL,
2879 NULL,
2880 rtpg_union_range_callback,
2881 &_raster
2882 );
2883 }
2884
2885 if (noerr != ES_NONE) {
2886 pfree(itrset);
2887 rtpg_union_arg_destroy(iwr);
2888 if (_rtn != NULL)
2889 rt_raster_destroy(_rtn);
2890 elog(ERROR, "RASTER_union_finalfn: Could not run raster iterator function");
2891 PG_RETURN_NULL();
2892 }
2893 }
2894 else {
2895 _raster = iwr->bandarg[i].raster[0];
2896 if (_raster == NULL)
2897 continue;
2898 }
2899
2900 /* first band, _rtn doesn't exist */
2901 if (i < 1) {
2902 uint32_t bandNums[1] = {0};
2903 _rtn = rt_raster_from_band(_raster, bandNums, 1);
2904 status = (_rtn == NULL) ? -1 : 0;
2905 }
2906 else
2907 status = rt_raster_copy_band(_rtn, _raster, 0, i);
2908
2909 POSTGIS_RT_DEBUG(4, "destroying source rasters");
2910
2911 /* destroy source rasters */
2912 if (
2913 iwr->bandarg[i].uniontype == UT_MEAN ||
2914 iwr->bandarg[i].uniontype == UT_RANGE
2915 ) {
2916 rt_raster_destroy(_raster);
2917 }
2918
2919 for (j = 0; j < iwr->bandarg[i].numraster; j++) {
2920 if (iwr->bandarg[i].raster[j] == NULL)
2921 continue;
2922 rt_raster_destroy(iwr->bandarg[i].raster[j]);
2923 iwr->bandarg[i].raster[j] = NULL;
2924 }
2925
2926 if (status < 0) {
2927 rtpg_union_arg_destroy(iwr);
2928 rt_raster_destroy(_rtn);
2929 elog(ERROR, "RASTER_union_finalfn: Could not add band to final raster");
2930 PG_RETURN_NULL();
2931 }
2932 }
2933
2934 /* cleanup */
2935 /* For Windowing functions, it is important to leave */
2936 /* the state intact, knowing that the aggcontext will be */
2937 /* freed by PgSQL when the statement is complete. */
2938 /* https://trac.osgeo.org/postgis/ticket/4770 */
2939 // pfree(itrset);
2940 // rtpg_union_arg_destroy(iwr);
2941
2942 if (!_rtn) PG_RETURN_NULL();
2943
2944 pgraster = rt_raster_serialize(_rtn);
2945 rt_raster_destroy(_rtn);
2946
2947 POSTGIS_RT_DEBUG(3, "Finished");
2948
2949 if (!pgraster)
2950 PG_RETURN_NULL();
2951
2952 SET_VARSIZE(pgraster, pgraster->size);
2953 PG_RETURN_POINTER(pgraster);
2954 }
2955
2956 /* ---------------------------------------------------------------- */
2957 /* Clip raster with geometry */
2958 /* ---------------------------------------------------------------- */
2959
2960 typedef struct rtpg_clip_band_t *rtpg_clip_band;
2961 struct rtpg_clip_band_t {
2962 int nband; /* band index */
2963 int hasnodata; /* is there a user-specified NODATA? */
2964 double nodataval; /* user-specified NODATA */
2965 };
2966
2967 typedef struct rtpg_clip_arg_t *rtpg_clip_arg;
2968 struct rtpg_clip_arg_t {
2969 rt_extenttype extenttype;
2970 rt_raster raster;
2971 rt_raster mask;
2972
2973 int numbands; /* number of bandargs */
2974 rtpg_clip_band band;
2975 };
2976
2977 static rtpg_clip_arg rtpg_clip_arg_init() {
2978 rtpg_clip_arg arg = NULL;
2979
2980 arg = palloc(sizeof(struct rtpg_clip_arg_t));
2981 if (arg == NULL) {
2982 elog(ERROR, "rtpg_clip_arg_init: Could not allocate memory for function arguments");
2983 return NULL;
2984 }
2985
2986 arg->extenttype = ET_INTERSECTION;
2987 arg->raster = NULL;
2988 arg->mask = NULL;
2989 arg->numbands = 0;
2990 arg->band = NULL;
2991
2992 return arg;
2993 }
2994
2995 static void rtpg_clip_arg_destroy(rtpg_clip_arg arg) {
2996 if (arg->band != NULL)
2997 pfree(arg->band);
2998
2999 if (arg->raster != NULL)
3000 rt_raster_destroy(arg->raster);
3001 if (arg->mask != NULL)
3002 rt_raster_destroy(arg->mask);
3003
3004 pfree(arg);
3005 }
3006
3007 static int rtpg_clip_callback(
3008 rt_iterator_arg arg, void *userarg,
3009 double *value, int *nodata
3010 ) {
3011 *value = 0;
3012 *nodata = 0;
3013
3014 /* either is NODATA, output is NODATA */
3015 if (arg->nodata[0][0][0] || arg->nodata[1][0][0])
3016 *nodata = 1;
3017 /* set to value */
3018 else
3019 *value = arg->values[0][0][0];
3020
3021 return 1;
3022 }
3023
3024 PG_FUNCTION_INFO_V1(RASTER_clip);
3025 Datum RASTER_clip(PG_FUNCTION_ARGS)
3026 {
3027 rt_pgraster *pgraster = NULL;
3028 LWGEOM *rastgeom = NULL;
3029 double gt[6] = {0};
3030 int32_t srid = SRID_UNKNOWN;
3031
3032 rt_pgraster *pgrtn = NULL;
3033 rt_raster rtn = NULL;
3034
3035 GSERIALIZED *gser = NULL;
3036 LWGEOM *geom = NULL;
3037 lwvarlena_t *wkb = NULL;
3038
3039 ArrayType *array;
3040 Oid etype;
3041 Datum *e;
3042 bool *nulls;
3043
3044 int16 typlen;
3045 bool typbyval;
3046 char typalign;
3047
3048 int i = 0;
3049 int j = 0;
3050 int k = 0;
3051 rtpg_clip_arg arg = NULL;
3052 LWGEOM *tmpgeom = NULL;
3053 rt_iterator itrset;
3054
3055 rt_raster _raster = NULL;
3056 rt_band band = NULL;
3057 rt_pixtype pixtype;
3058 int hasnodata;
3059 double nodataval;
3060 int noerr = 0;
3061
3062 POSTGIS_RT_DEBUG(3, "Starting...");
3063
3064 /* raster or geometry is NULL, return NULL */
3065 if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3066 PG_RETURN_NULL();
3067
3068 /* init arg */
3069 arg = rtpg_clip_arg_init();
3070 if (arg == NULL) {
3071 elog(ERROR, "RASTER_clip: Could not initialize argument structure");
3072 PG_RETURN_NULL();
3073 }
3074
3075 /* raster (0) */
3076 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3077
3078 /* Get raster object */
3079 arg->raster = rt_raster_deserialize(pgraster, FALSE);
3080 if (arg->raster == NULL) {
3081 rtpg_clip_arg_destroy(arg);
3082 PG_FREE_IF_COPY(pgraster, 0);
3083 elog(ERROR, "RASTER_clip: Could not deserialize raster");
3084 PG_RETURN_NULL();
3085 }
3086
3087 /* raster is empty, return empty raster */
3088 if (rt_raster_is_empty(arg->raster) || rt_raster_get_num_bands(arg->raster) == 0) {
3089 elog(NOTICE, "Input raster is empty or has no bands. Returning empty raster");
3090
3091 rtpg_clip_arg_destroy(arg);
3092 PG_FREE_IF_COPY(pgraster, 0);
3093
3094 rtn = rt_raster_new(0, 0);
3095 if (rtn == NULL) {
3096 elog(ERROR, "RASTER_clip: Could not create empty raster");
3097 PG_RETURN_NULL();
3098 }
3099
3100 pgrtn = rt_raster_serialize(rtn);
3101 rt_raster_destroy(rtn);
3102 if (NULL == pgrtn)
3103 PG_RETURN_NULL();
3104
3105 SET_VARSIZE(pgrtn, pgrtn->size);
3106 PG_RETURN_POINTER(pgrtn);
3107 }
3108
3109 /* metadata */
3110 rt_raster_get_geotransform_matrix(arg->raster, gt);
3111 srid = clamp_srid(rt_raster_get_srid(arg->raster));
3112
3113 /* geometry (2) */
3114 gser = PG_GETARG_GSERIALIZED_P(2);
3115 geom = lwgeom_from_gserialized(gser);
3116
3117 /* Get a 2D version of the geometry if necessary */
3118 if (lwgeom_ndims(geom) > 2) {
3119 LWGEOM *geom2d = lwgeom_force_2d(geom);
3120 lwgeom_free(geom);
3121 geom = geom2d;
3122 }
3123
3124 /* check that SRIDs match */
3125 if (srid != clamp_srid(gserialized_get_srid(gser))) {
3126 elog(NOTICE, "Geometry provided does not have the same SRID as the raster. Returning NULL");
3127
3128 rtpg_clip_arg_destroy(arg);
3129 PG_FREE_IF_COPY(pgraster, 0);
3130 lwgeom_free(geom);
3131 PG_FREE_IF_COPY(gser, 2);
3132
3133 PG_RETURN_NULL();
3134 }
3135
3136 /* crop (4) */
3137 if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3138 arg->extenttype = ET_FIRST;
3139
3140 /* get intersection geometry of input raster and input geometry */
3141 if (rt_raster_get_convex_hull(arg->raster, &rastgeom) != ES_NONE) {
3142
3143 rtpg_clip_arg_destroy(arg);
3144 PG_FREE_IF_COPY(pgraster, 0);
3145 lwgeom_free(geom);
3146 PG_FREE_IF_COPY(gser, 2);
3147
3148 elog(ERROR, "RASTER_clip: Could not get convex hull of raster");
3149 PG_RETURN_NULL();
3150 }
3151
3152 tmpgeom = lwgeom_intersection(rastgeom, geom);
3153 lwgeom_free(rastgeom);
3154 lwgeom_free(geom);
3155 PG_FREE_IF_COPY(gser, 2);
3156 geom = tmpgeom;
3157
3158 /* intersection is empty AND extent type is INTERSECTION, return empty */
3159 if (lwgeom_is_empty(geom) && arg->extenttype == ET_INTERSECTION) {
3160 elog(NOTICE, "The input raster and input geometry do not intersect. Returning empty raster");
3161
3162 rtpg_clip_arg_destroy(arg);
3163 PG_FREE_IF_COPY(pgraster, 0);
3164 lwgeom_free(geom);
3165
3166 rtn = rt_raster_new(0, 0);
3167 if (rtn == NULL) {
3168 elog(ERROR, "RASTER_clip: Could not create empty raster");
3169 PG_RETURN_NULL();
3170 }
3171
3172 pgrtn = rt_raster_serialize(rtn);
3173 rt_raster_destroy(rtn);
3174 if (NULL == pgrtn)
3175 PG_RETURN_NULL();
3176
3177 SET_VARSIZE(pgrtn, pgrtn->size);
3178 PG_RETURN_POINTER(pgrtn);
3179 }
3180
3181 /* nband (1) */
3182 if (!PG_ARGISNULL(1)) {
3183 array = PG_GETARG_ARRAYTYPE_P(1);
3184 etype = ARR_ELEMTYPE(array);
3185 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3186
3187 switch (etype) {
3188 case INT2OID:
3189 case INT4OID:
3190 break;
3191 default:
3192 rtpg_clip_arg_destroy(arg);
3193 PG_FREE_IF_COPY(pgraster, 0);
3194 lwgeom_free(geom);
3195 elog(ERROR, "RASTER_clip: Invalid data type for band indexes");
3196 PG_RETURN_NULL();
3197 break;
3198 }
3199
3200 deconstruct_array(
3201 array, etype,
3202 typlen, typbyval, typalign,
3203 &e, &nulls, &(arg->numbands)
3204 );
3205
3206 arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3207 if (arg->band == NULL) {
3208 rtpg_clip_arg_destroy(arg);
3209 PG_FREE_IF_COPY(pgraster, 0);
3210 lwgeom_free(geom);
3211 elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3212 PG_RETURN_NULL();
3213 }
3214
3215 for (i = 0, j = 0; i < arg->numbands; i++) {
3216 if (nulls[i]) continue;
3217
3218 switch (etype) {
3219 case INT2OID:
3220 arg->band[j].nband = DatumGetInt16(e[i]) - 1;
3221 break;
3222 case INT4OID:
3223 arg->band[j].nband = DatumGetInt32(e[i]) - 1;
3224 break;
3225 }
3226
3227 j++;
3228 }
3229
3230 if (j < arg->numbands) {
3231 arg->band = repalloc(arg->band, sizeof(struct rtpg_clip_band_t) * j);
3232 if (arg->band == NULL) {
3233 rtpg_clip_arg_destroy(arg);
3234 PG_FREE_IF_COPY(pgraster, 0);
3235 lwgeom_free(geom);
3236 elog(ERROR, "RASTER_clip: Could not reallocate memory for band arguments");
3237 PG_RETURN_NULL();
3238 }
3239
3240 arg->numbands = j;
3241 }
3242
3243 /* validate band */
3244 for (i = 0; i < arg->numbands; i++) {
3245 if (!rt_raster_has_band(arg->raster, arg->band[i].nband)) {
3246 elog(NOTICE, "Band at index %d not found in raster", arg->band[i].nband + 1);
3247 rtpg_clip_arg_destroy(arg);
3248 PG_FREE_IF_COPY(pgraster, 0);
3249 lwgeom_free(geom);
3250 PG_RETURN_NULL();
3251 }
3252
3253 arg->band[i].hasnodata = 0;
3254 arg->band[i].nodataval = 0;
3255 }
3256 }
3257 else {
3258 arg->numbands = rt_raster_get_num_bands(arg->raster);
3259
3260 /* raster may have no bands */
3261 if (arg->numbands) {
3262 arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3263 if (arg->band == NULL) {
3264
3265 rtpg_clip_arg_destroy(arg);
3266 PG_FREE_IF_COPY(pgraster, 0);
3267 lwgeom_free(geom);
3268
3269 elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3270 PG_RETURN_NULL();
3271 }
3272
3273 for (i = 0; i < arg->numbands; i++) {
3274 arg->band[i].nband = i;
3275 arg->band[i].hasnodata = 0;
3276 arg->band[i].nodataval = 0;
3277 }
3278 }
3279 }
3280
3281 /* nodataval (3) */
3282 if (!PG_ARGISNULL(3)) {
3283 array = PG_GETARG_ARRAYTYPE_P(3);
3284 etype = ARR_ELEMTYPE(array);
3285 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3286
3287 switch (etype) {
3288 case FLOAT4OID:
3289 case FLOAT8OID:
3290 break;
3291 default:
3292 rtpg_clip_arg_destroy(arg);
3293 PG_FREE_IF_COPY(pgraster, 0);
3294 lwgeom_free(geom);
3295 elog(ERROR, "RASTER_clip: Invalid data type for NODATA values");
3296 PG_RETURN_NULL();
3297 break;
3298 }
3299
3300 deconstruct_array(
3301 array, etype,
3302 typlen, typbyval, typalign,
3303 &e, &nulls, &k
3304 );
3305
3306 /* it doesn't matter if there are more nodataval */
3307 for (i = 0, j = 0; i < arg->numbands; i++, j++) {
3308 /* cap j to the last nodataval element */
3309 if (j >= k)
3310 j = k - 1;
3311
3312 if (nulls[j])
3313 continue;
3314
3315 arg->band[i].hasnodata = 1;
3316 switch (etype) {
3317 case FLOAT4OID:
3318 arg->band[i].nodataval = DatumGetFloat4(e[j]);
3319 break;
3320 case FLOAT8OID:
3321 arg->band[i].nodataval = DatumGetFloat8(e[j]);
3322 break;
3323 }
3324 }
3325 }
3326
3327 /* get wkb of geometry */
3328 POSTGIS_RT_DEBUG(3, "getting wkb of geometry");
3329 wkb = lwgeom_to_wkb_varlena(geom, WKB_SFSQL);
3330 lwgeom_free(geom);
3331
3332 /* rasterize geometry */
3333 arg->mask = rt_raster_gdal_rasterize((unsigned char *)wkb->data,
3334 LWSIZE_GET(wkb->size) - LWVARHDRSZ,
3335 NULL,
3336 0,
3337 NULL,
3338 NULL,
3339 NULL,
3340 NULL,
3341 NULL,
3342 NULL,
3343 NULL,
3344 &(gt[1]),
3345 &(gt[5]),
3346 NULL,
3347 NULL,
3348 &(gt[0]),
3349 &(gt[3]),
3350 &(gt[2]),
3351 &(gt[4]),
3352 NULL);
3353
3354 pfree(wkb);
3355 if (arg->mask == NULL) {
3356 rtpg_clip_arg_destroy(arg);
3357 PG_FREE_IF_COPY(pgraster, 0);
3358 elog(ERROR, "RASTER_clip: Could not rasterize intersection geometry");
3359 PG_RETURN_NULL();
3360 }
3361
3362 /* set SRID */
3363 rt_raster_set_srid(arg->mask, srid);
3364
3365 /* run iterator */
3366
3367 /* init itrset */
3368 itrset = palloc(sizeof(struct rt_iterator_t) * 2);
3369 if (itrset == NULL) {
3370 rtpg_clip_arg_destroy(arg);
3371 PG_FREE_IF_COPY(pgraster, 0);
3372 elog(ERROR, "RASTER_clip: Could not allocate memory for iterator arguments");
3373 PG_RETURN_NULL();
3374 }
3375
3376 /* one band at a time */
3377 for (i = 0; i < arg->numbands; i++) {
3378 POSTGIS_RT_DEBUGF(4, "band arg %d (nband, hasnodata, nodataval) = (%d, %d, %f)",
3379 i, arg->band[i].nband, arg->band[i].hasnodata, arg->band[i].nodataval);
3380
3381 band = rt_raster_get_band(arg->raster, arg->band[i].nband);
3382
3383 /* band metadata */
3384 pixtype = rt_band_get_pixtype(band);
3385
3386 if (arg->band[i].hasnodata) {
3387 hasnodata = 1;
3388 nodataval = arg->band[i].nodataval;
3389 }
3390 else if (rt_band_get_hasnodata_flag(band)) {
3391 hasnodata = 1;
3392 rt_band_get_nodata(band, &nodataval);
3393 }
3394 else {
3395 hasnodata = 0;
3396 nodataval = rt_band_get_min_value(band);
3397 }
3398
3399 /* band is NODATA, create NODATA band and continue */
3400 if (rt_band_get_isnodata_flag(band)) {
3401 /* create raster */
3402 if (rtn == NULL) {
3403 noerr = rt_raster_from_two_rasters(arg->raster, arg->mask, arg->extenttype, &rtn, NULL);
3404 if (noerr != ES_NONE) {
3405 rtpg_clip_arg_destroy(arg);
3406 PG_FREE_IF_COPY(pgraster, 0);
3407 elog(ERROR, "RASTER_clip: Could not create output raster");
3408 PG_RETURN_NULL();
3409 }
3410 }
3411
3412 /* create NODATA band */
3413 if (rt_raster_generate_new_band(rtn, pixtype, nodataval, hasnodata, nodataval, i) < 0) {
3414 rt_raster_destroy(rtn);
3415 rtpg_clip_arg_destroy(arg);
3416 PG_FREE_IF_COPY(pgraster, 0);
3417 elog(ERROR, "RASTER_clip: Could not add NODATA band to output raster");
3418 PG_RETURN_NULL();
3419 }
3420
3421 continue;
3422 }
3423
3424 /* raster */
3425 itrset[0].raster = arg->raster;
3426 itrset[0].nband = arg->band[i].nband;
3427 itrset[0].nbnodata = 1;
3428
3429 /* mask */
3430 itrset[1].raster = arg->mask;
3431 itrset[1].nband = 0;
3432 itrset[1].nbnodata = 1;
3433
3434 /* pass to iterator */
3435 noerr = rt_raster_iterator(
3436 itrset, 2,
3437 arg->extenttype, NULL,
3438 pixtype,
3439 hasnodata, nodataval,
3440 0, 0,
3441 NULL,
3442 NULL,
3443 rtpg_clip_callback,
3444 &_raster
3445 );
3446
3447 if (noerr != ES_NONE) {
3448 pfree(itrset);
3449 rtpg_clip_arg_destroy(arg);
3450 if (rtn != NULL) rt_raster_destroy(rtn);
3451 PG_FREE_IF_COPY(pgraster, 0);
3452 elog(ERROR, "RASTER_clip: Could not run raster iterator function");
3453 PG_RETURN_NULL();
3454 }
3455
3456 /* new raster */
3457 if (rtn == NULL)
3458 rtn = _raster;
3459 /* copy band */
3460 else {
3461 band = rt_raster_get_band(_raster, 0);
3462 if (band == NULL) {
3463 pfree(itrset);
3464 rtpg_clip_arg_destroy(arg);
3465 rt_raster_destroy(_raster);
3466 rt_raster_destroy(rtn);
3467 PG_FREE_IF_COPY(pgraster, 0);
3468 elog(NOTICE, "RASTER_clip: Could not get band from working raster");
3469 PG_RETURN_NULL();
3470 }
3471
3472 if (rt_raster_add_band(rtn, band, i) < 0) {
3473 pfree(itrset);
3474 rtpg_clip_arg_destroy(arg);
3475 rt_raster_destroy(_raster);
3476 rt_raster_destroy(rtn);
3477 PG_FREE_IF_COPY(pgraster, 0);
3478 elog(ERROR, "RASTER_clip: Could not add new band to output raster");
3479 PG_RETURN_NULL();
3480 }
3481
3482 rt_raster_destroy(_raster);
3483 }
3484 }
3485
3486 pfree(itrset);
3487 rtpg_clip_arg_destroy(arg);
3488 PG_FREE_IF_COPY(pgraster, 0);
3489
3490 pgrtn = rt_raster_serialize(rtn);
3491 rt_raster_destroy(rtn);
3492
3493 POSTGIS_RT_DEBUG(3, "Finished");
3494
3495 if (!pgrtn)
3496 PG_RETURN_NULL();
3497
3498 SET_VARSIZE(pgrtn, pgrtn->size);
3499 PG_RETURN_POINTER(pgrtn);
3500 }
3501
3502 /**
3503 * Reclassify the specified bands of the raster
3504 */
3505 PG_FUNCTION_INFO_V1(RASTER_reclass);
3506 Datum RASTER_reclass(PG_FUNCTION_ARGS) {
3507 rt_pgraster *pgraster = NULL;
3508 rt_pgraster *pgrtn = NULL;
3509 rt_raster raster = NULL;
3510 rt_band band = NULL;
3511 rt_band newband = NULL;
3512 uint32_t numBands = 0;
3513
3514 ArrayType *array;
3515 Oid etype;
3516 Datum *e;
3517 bool *nulls;
3518 int16 typlen;
3519 bool typbyval;
3520 char typalign;
3521 int n = 0;
3522
3523 int i = 0;
3524 int j = 0;
3525 int k = 0;
3526
3527 uint32_t a = 0;
3528 uint32_t b = 0;
3529 uint32_t c = 0;
3530
3531 rt_reclassexpr *exprset = NULL;
3532 HeapTupleHeader tup;
3533 bool isnull;
3534 Datum tupv;
3535 uint32_t nband = 0;
3536 char *expr = NULL;
3537 text *exprtext = NULL;
3538 double val = 0;
3539 char *junk = NULL;
3540 int inc_val = 0;
3541 int exc_val = 0;
3542 char *pixeltype = NULL;
3543 text *pixeltypetext = NULL;
3544 rt_pixtype pixtype = PT_END;
3545 double nodataval = 0;
3546 bool hasnodata = FALSE;
3547
3548 char **comma_set = NULL;
3549 uint32_t comma_n = 0;
3550 char **colon_set = NULL;
3551 uint32_t colon_n = 0;
3552 char **dash_set = NULL;
3553 uint32_t dash_n = 0;
3554
3555 POSTGIS_RT_DEBUG(3, "RASTER_reclass: Starting");
3556
3557 /* pgraster is null, return null */
3558 if (PG_ARGISNULL(0))
3559 PG_RETURN_NULL();
3560 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3561
3562 /* raster */
3563 raster = rt_raster_deserialize(pgraster, FALSE);
3564 if (!raster) {
3565 PG_FREE_IF_COPY(pgraster, 0);
3566 elog(ERROR, "RASTER_reclass: Could not deserialize raster");
3567 PG_RETURN_NULL();
3568 }
3569 numBands = rt_raster_get_num_bands(raster);
3570 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: %d possible bands to be reclassified", numBands);
3571
3572 /* process set of reclassarg */
3573 POSTGIS_RT_DEBUG(3, "RASTER_reclass: Processing Arg 1 (reclassargset)");
3574 array = PG_GETARG_ARRAYTYPE_P(1);
3575 etype = ARR_ELEMTYPE(array);
3576 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3577
3578 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
3579 &nulls, &n);
3580
3581 if (!n) {
3582 elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
3583
3584 pgrtn = rt_raster_serialize(raster);
3585 rt_raster_destroy(raster);
3586 PG_FREE_IF_COPY(pgraster, 0);
3587 if (!pgrtn)
3588 PG_RETURN_NULL();
3589
3590 SET_VARSIZE(pgrtn, pgrtn->size);
3591 PG_RETURN_POINTER(pgrtn);
3592 }
3593
3594 /*
3595 process each element of reclassarg
3596 each element is the index of the band to process, the set
3597 of reclass ranges and the output band's pixeltype
3598 */
3599 for (i = 0; i < n; i++) {
3600 if (nulls[i]) continue;
3601
3602 /* each element is a tuple */
3603 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
3604 if (NULL == tup) {
3605 elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
3606
3607 pgrtn = rt_raster_serialize(raster);
3608 rt_raster_destroy(raster);
3609 PG_FREE_IF_COPY(pgraster, 0);
3610 if (!pgrtn)
3611 PG_RETURN_NULL();
3612
3613 SET_VARSIZE(pgrtn, pgrtn->size);
3614 PG_RETURN_POINTER(pgrtn);
3615 }
3616
3617 /* band index (1-based) */
3618 tupv = GetAttributeByName(tup, "nband", &isnull);
3619 if (isnull) {
3620 elog(NOTICE, "Invalid argument for reclassargset. Missing value of nband for reclassarg of index %d . Returning original raster", i);
3621
3622 pgrtn = rt_raster_serialize(raster);
3623 rt_raster_destroy(raster);
3624 PG_FREE_IF_COPY(pgraster, 0);
3625 if (!pgrtn)
3626 PG_RETURN_NULL();
3627
3628 SET_VARSIZE(pgrtn, pgrtn->size);
3629 PG_RETURN_POINTER(pgrtn);
3630 }
3631 nband = DatumGetInt32(tupv);
3632 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: expression for band %d", nband);
3633
3634 /* valid band index? */
3635 if (nband < 1 || nband > numBands) {
3636 elog(NOTICE, "Invalid argument for reclassargset. Invalid band index (must use 1-based) for reclassarg of index %d . Returning original raster", i);
3637
3638 pgrtn = rt_raster_serialize(raster);
3639 rt_raster_destroy(raster);
3640 PG_FREE_IF_COPY(pgraster, 0);
3641 if (!pgrtn)
3642 PG_RETURN_NULL();
3643
3644 SET_VARSIZE(pgrtn, pgrtn->size);
3645 PG_RETURN_POINTER(pgrtn);
3646 }
3647
3648 /* reclass expr */
3649 tupv = GetAttributeByName(tup, "reclassexpr", &isnull);
3650 if (isnull) {
3651 elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3652
3653 pgrtn = rt_raster_serialize(raster);
3654 rt_raster_destroy(raster);
3655 PG_FREE_IF_COPY(pgraster, 0);
3656 if (!pgrtn)
3657 PG_RETURN_NULL();
3658
3659 SET_VARSIZE(pgrtn, pgrtn->size);
3660 PG_RETURN_POINTER(pgrtn);
3661 }
3662 exprtext = (text *) DatumGetPointer(tupv);
3663 if (NULL == exprtext) {
3664 elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3665
3666 pgrtn = rt_raster_serialize(raster);
3667 rt_raster_destroy(raster);
3668 PG_FREE_IF_COPY(pgraster, 0);
3669 if (!pgrtn)
3670 PG_RETURN_NULL();
3671
3672 SET_VARSIZE(pgrtn, pgrtn->size);
3673 PG_RETURN_POINTER(pgrtn);
3674 }
3675 expr = text_to_cstring(exprtext);
3676 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: expr (raw) %s", expr);
3677 expr = rtpg_removespaces(expr);
3678 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: expr (clean) %s", expr);
3679
3680 /* split string to its components */
3681 /* comma (,) separating rangesets */
3682 comma_set = rtpg_strsplit(expr, ",", &comma_n);
3683 if (comma_n < 1) {
3684 elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3685
3686 pgrtn = rt_raster_serialize(raster);
3687 rt_raster_destroy(raster);
3688 PG_FREE_IF_COPY(pgraster, 0);
3689 if (!pgrtn)
3690 PG_RETURN_NULL();
3691
3692 SET_VARSIZE(pgrtn, pgrtn->size);
3693 PG_RETURN_POINTER(pgrtn);
3694 }
3695
3696 /* set of reclass expressions */
3697 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: %d possible expressions", comma_n);
3698 exprset = palloc(comma_n * sizeof(rt_reclassexpr));
3699
3700 for (a = 0, j = 0; a < comma_n; a++) {
3701 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: map %s", comma_set[a]);
3702
3703 /* colon (:) separating range "src" and "dst" */
3704 colon_set = rtpg_strsplit(comma_set[a], ":", &colon_n);
3705 if (colon_n != 2) {
3706 elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3707 for (k = 0; k < j; k++) pfree(exprset[k]);
3708 pfree(exprset);
3709
3710 pgrtn = rt_raster_serialize(raster);
3711 rt_raster_destroy(raster);
3712 PG_FREE_IF_COPY(pgraster, 0);
3713 if (!pgrtn)
3714 PG_RETURN_NULL();
3715
3716 SET_VARSIZE(pgrtn, pgrtn->size);
3717 PG_RETURN_POINTER(pgrtn);
3718 }
3719
3720 /* allocate mem for reclass expression */
3721 exprset[j] = palloc(sizeof(struct rt_reclassexpr_t));
3722
3723 for (b = 0; b < colon_n; b++) {
3724 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: range %s", colon_set[b]);
3725
3726 /* dash (-) separating "min" and "max" */
3727 dash_set = rtpg_strsplit(colon_set[b], "-", &dash_n);
3728 if (dash_n < 1 || dash_n > 3) {
3729 elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3730 for (k = 0; k < j; k++) pfree(exprset[k]);
3731 pfree(exprset);
3732
3733 pgrtn = rt_raster_serialize(raster);
3734 rt_raster_destroy(raster);
3735 PG_FREE_IF_COPY(pgraster, 0);
3736 if (!pgrtn)
3737 PG_RETURN_NULL();
3738
3739 SET_VARSIZE(pgrtn, pgrtn->size);
3740 PG_RETURN_POINTER(pgrtn);
3741 }
3742
3743 for (c = 0; c < dash_n; c++) {
3744 /* need to handle: (-9999-100 -> "(", "9999", "100" */
3745 if (
3746 c < 1 &&
3747 strlen(dash_set[c]) == 1 && (
3748 strchr(dash_set[c], '(') != NULL ||
3749 strchr(dash_set[c], '[') != NULL ||
3750 strchr(dash_set[c], ')') != NULL ||
3751 strchr(dash_set[c], ']') != NULL
3752 )
3753 ) {
3754 uint32_t dash_it;
3755 junk = palloc(sizeof(char) * (strlen(dash_set[c + 1]) + 2));
3756 if (NULL == junk) {
3757 for (k = 0; k <= j; k++) pfree(exprset[k]);
3758 pfree(exprset);
3759 rt_raster_destroy(raster);
3760 PG_FREE_IF_COPY(pgraster, 0);
3761
3762 elog(ERROR, "RASTER_reclass: Could not allocate memory");
3763 PG_RETURN_NULL();
3764 }
3765
3766 sprintf(junk, "%s%s", dash_set[c], dash_set[c + 1]);
3767 c++;
3768 dash_set[c] = repalloc(dash_set[c], sizeof(char) * (strlen(junk) + 1));
3769 strcpy(dash_set[c], junk);
3770 pfree(junk);
3771
3772 /* rebuild dash_set */
3773 for (dash_it = 1; dash_it < dash_n; dash_it++) {
3774 dash_set[dash_it - 1] = repalloc(dash_set[dash_it - 1], (strlen(dash_set[dash_it]) + 1) * sizeof(char));
3775 strcpy(dash_set[dash_it - 1], dash_set[dash_it]);
3776 }
3777 dash_n--;
3778 c--;
3779 pfree(dash_set[dash_n]);
3780 dash_set = repalloc(dash_set, sizeof(char *) * dash_n);
3781 }
3782
3783 /* there shouldn't be more than two in dash_n */
3784 if (c < 1 && dash_n > 2) {
3785 elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3786 for (k = 0; k < j; k++) pfree(exprset[k]);
3787 pfree(exprset);
3788
3789 pgrtn = rt_raster_serialize(raster);
3790 rt_raster_destroy(raster);
3791 PG_FREE_IF_COPY(pgraster, 0);
3792 if (!pgrtn)
3793 PG_RETURN_NULL();
3794
3795 SET_VARSIZE(pgrtn, pgrtn->size);
3796 PG_RETURN_POINTER(pgrtn);
3797 }
3798
3799 /* check interval flags */
3800 exc_val = 0;
3801 inc_val = 1;
3802 /* range */
3803 if (dash_n != 1) {
3804 /* min */
3805 if (c < 1) {
3806 if (
3807 strchr(dash_set[c], ')') != NULL ||
3808 strchr(dash_set[c], ']') != NULL
3809 ) {
3810 exc_val = 1;
3811 inc_val = 1;
3812 }
3813 else if (strchr(dash_set[c], '(') != NULL){
3814 inc_val = 0;
3815 }
3816 else {
3817 inc_val = 1;
3818 }
3819 }
3820 /* max */
3821 else {
3822 if (
3823 strrchr(dash_set[c], '(') != NULL ||
3824 strrchr(dash_set[c], '[') != NULL
3825 ) {
3826 exc_val = 1;
3827 inc_val = 0;
3828 }
3829 else if (strrchr(dash_set[c], ']') != NULL) {
3830 inc_val = 1;
3831 }
3832 else {
3833 inc_val = 0;
3834 }
3835 }
3836 }
3837 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: exc_val %d inc_val %d", exc_val, inc_val);
3838
3839 /* remove interval flags */
3840 dash_set[c] = rtpg_chartrim(dash_set[c], "()[]");
3841 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (char) %s", dash_set[c]);
3842
3843 /* value from string to double */
3844 errno = 0;
3845 val = strtod(dash_set[c], &junk);
3846 if (errno != 0 || dash_set[c] == junk) {
3847 elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3848 for (k = 0; k < j; k++) pfree(exprset[k]);
3849 pfree(exprset);
3850
3851 pgrtn = rt_raster_serialize(raster);
3852 rt_raster_destroy(raster);
3853 PG_FREE_IF_COPY(pgraster, 0);
3854 if (!pgrtn)
3855 PG_RETURN_NULL();
3856
3857 SET_VARSIZE(pgrtn, pgrtn->size);
3858 PG_RETURN_POINTER(pgrtn);
3859 }
3860 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (double) %f", val);
3861
3862 /* strsplit removes dash (a.k.a. negative sign), compare now to restore */
3863 if (c < 1)
3864 junk = strstr(colon_set[b], dash_set[c]);
3865 else
3866 junk = rtpg_strrstr(colon_set[b], dash_set[c]);
3867 POSTGIS_RT_DEBUGF(
3868 4,
3869 "(colon_set[%d], dash_set[%d], junk) = (%s, %s, %s)",
3870 b, c, colon_set[b], dash_set[c], junk
3871 );
3872 /* not beginning of string */
3873 if (junk != colon_set[b]) {
3874 /* prior is a dash */
3875 if (*(junk - 1) == '-') {
3876 /* prior is beginning of string or prior - 1 char is dash, negative number */
3877 if (
3878 ((junk - 1) == colon_set[b]) ||
3879 (*(junk - 2) == '-') ||
3880 (*(junk - 2) == '[') ||
3881 (*(junk - 2) == '(')
3882 ) {
3883 val *= -1.;
3884 }
3885 }
3886 }
3887 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (double) %f", val);
3888
3889 /* src */
3890 if (b < 1) {
3891 /* singular value */
3892 if (dash_n == 1) {
3893 exprset[j]->src.exc_min = exprset[j]->src.exc_max = exc_val;
3894 exprset[j]->src.inc_min = exprset[j]->src.inc_max = inc_val;
3895 exprset[j]->src.min = exprset[j]->src.max = val;
3896 }
3897 /* min */
3898 else if (c < 1) {
3899 exprset[j]->src.exc_min = exc_val;
3900 exprset[j]->src.inc_min = inc_val;
3901 exprset[j]->src.min = val;
3902 }
3903 /* max */
3904 else {
3905 exprset[j]->src.exc_max = exc_val;
3906 exprset[j]->src.inc_max = inc_val;
3907 exprset[j]->src.max = val;
3908 }
3909 }
3910 /* dst */
3911 else {
3912 /* singular value */
3913 if (dash_n == 1)
3914 exprset[j]->dst.min = exprset[j]->dst.max = val;
3915 /* min */
3916 else if (c < 1)
3917 exprset[j]->dst.min = val;
3918 /* max */
3919 else
3920 exprset[j]->dst.max = val;
3921 }
3922 }
3923 pfree(dash_set);
3924 }
3925 pfree(colon_set);
3926
3927 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: or: %f - %f nr: %f - %f"
3928 , exprset[j]->src.min
3929 , exprset[j]->src.max
3930 , exprset[j]->dst.min
3931 , exprset[j]->dst.max
3932 );
3933 j++;
3934 }
3935 pfree(comma_set);
3936
3937 /* pixel type */
3938 tupv = GetAttributeByName(tup, "pixeltype", &isnull);
3939 if (isnull) {
3940 elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3941
3942 pgrtn = rt_raster_serialize(raster);
3943 rt_raster_destroy(raster);
3944 PG_FREE_IF_COPY(pgraster, 0);
3945 if (!pgrtn)
3946 PG_RETURN_NULL();
3947
3948 SET_VARSIZE(pgrtn, pgrtn->size);
3949 PG_RETURN_POINTER(pgrtn);
3950 }
3951 pixeltypetext = (text *) DatumGetPointer(tupv);
3952 if (NULL == pixeltypetext) {
3953 elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3954
3955 pgrtn = rt_raster_serialize(raster);
3956 rt_raster_destroy(raster);
3957 PG_FREE_IF_COPY(pgraster, 0);
3958 if (!pgrtn)
3959 PG_RETURN_NULL();
3960
3961 SET_VARSIZE(pgrtn, pgrtn->size);
3962 PG_RETURN_POINTER(pgrtn);
3963 }
3964 pixeltype = text_to_cstring(pixeltypetext);
3965 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: pixeltype %s", pixeltype);
3966 pixtype = rt_pixtype_index_from_name(pixeltype);
3967
3968 /* nodata */
3969 tupv = GetAttributeByName(tup, "nodataval", &isnull);
3970 if (isnull) {
3971 nodataval = 0;
3972 hasnodata = FALSE;
3973 }
3974 else {
3975 nodataval = DatumGetFloat8(tupv);
3976 hasnodata = TRUE;
3977 }
3978 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: nodataval %f", nodataval);
3979 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: hasnodata %d", hasnodata);
3980
3981 /* do reclass */
3982 band = rt_raster_get_band(raster, nband - 1);
3983 if (!band) {
3984 elog(NOTICE, "Could not find raster band of index %d. Returning original raster", nband);
3985 for (k = 0; k < j; k++) pfree(exprset[k]);
3986 pfree(exprset);
3987
3988 pgrtn = rt_raster_serialize(raster);
3989 rt_raster_destroy(raster);
3990 PG_FREE_IF_COPY(pgraster, 0);
3991 if (!pgrtn)
3992 PG_RETURN_NULL();
3993
3994 SET_VARSIZE(pgrtn, pgrtn->size);
3995 PG_RETURN_POINTER(pgrtn);
3996 }
3997 newband = rt_band_reclass(band, pixtype, hasnodata, nodataval, exprset, j);
3998 if (!newband) {
3999 for (k = 0; k < j; k++) pfree(exprset[k]);
4000 pfree(exprset);
4001
4002 rt_raster_destroy(raster);
4003 PG_FREE_IF_COPY(pgraster, 0);
4004
4005 elog(ERROR, "RASTER_reclass: Could not reclassify raster band of index %d", nband);
4006 PG_RETURN_NULL();
4007 }
4008
4009 /* replace old band with new band */
4010 if (rt_raster_replace_band(raster, newband, nband - 1) == NULL) {
4011 for (k = 0; k < j; k++) pfree(exprset[k]);
4012 pfree(exprset);
4013
4014 rt_band_destroy(newband);
4015 rt_raster_destroy(raster);
4016 PG_FREE_IF_COPY(pgraster, 0);
4017
4018 elog(ERROR, "RASTER_reclass: Could not replace raster band of index %d with reclassified band", nband);
4019 PG_RETURN_NULL();
4020 }
4021
4022 /* old band is in the variable band */
4023 rt_band_destroy(band);
4024
4025 /* free exprset */
4026 for (k = 0; k < j; k++) pfree(exprset[k]);
4027 pfree(exprset);
4028 }
4029
4030 pgrtn = rt_raster_serialize(raster);
4031 rt_raster_destroy(raster);
4032 PG_FREE_IF_COPY(pgraster, 0);
4033 if (!pgrtn)
4034 PG_RETURN_NULL();
4035
4036 POSTGIS_RT_DEBUG(3, "RASTER_reclass: Finished");
4037
4038 SET_VARSIZE(pgrtn, pgrtn->size);
4039 PG_RETURN_POINTER(pgrtn);
4040 }
4041
4042 /* ---------------------------------------------------------------- */
4043 /* apply colormap to specified band of a raster */
4044 /* ---------------------------------------------------------------- */
4045
4046 typedef struct rtpg_colormap_arg_t *rtpg_colormap_arg;
4047 struct rtpg_colormap_arg_t {
4048 rt_raster raster;
4049 int nband; /* 1-based */
4050 rt_band band;
4051 rt_bandstats bandstats;
4052
4053 rt_colormap colormap;
4054 int nodataentry;
4055
4056 char **entry;
4057 uint32_t nentry;
4058 char **element;
4059 uint32_t nelement;
4060 };
4061
4062 static rtpg_colormap_arg
4063 rtpg_colormap_arg_init() {
4064 rtpg_colormap_arg arg = NULL;
4065
4066 arg = palloc(sizeof(struct rtpg_colormap_arg_t));
4067 if (arg == NULL) {
4068 elog(ERROR, "rtpg_colormap_arg: Could not allocate memory for function arguments");
4069 return NULL;
4070 }
4071
4072 arg->raster = NULL;
4073 arg->nband = 1;
4074 arg->band = NULL;
4075 arg->bandstats = NULL;
4076
4077 arg->colormap = palloc(sizeof(struct rt_colormap_t));
4078 if (arg->colormap == NULL) {
4079 elog(ERROR, "rtpg_colormap_arg: Could not allocate memory for function arguments");
4080 return NULL;
4081 }
4082 arg->colormap->nentry = 0;
4083 arg->colormap->entry = NULL;
4084 arg->colormap->ncolor = 4; /* assume RGBA */
4085 arg->colormap->method = CM_INTERPOLATE;
4086 arg->nodataentry = -1;
4087
4088 arg->entry = NULL;
4089 arg->nentry = 0;
4090 arg->element = NULL;
4091 arg->nelement = 0;
4092
4093 return arg;
4094 }
4095
4096 static void
4097 rtpg_colormap_arg_destroy(rtpg_colormap_arg arg) {
4098 uint32_t i = 0;
4099 if (arg->raster != NULL)
4100 rt_raster_destroy(arg->raster);
4101
4102 if (arg->bandstats != NULL)
4103 pfree(arg->bandstats);
4104
4105 if (arg->colormap != NULL) {
4106 if (arg->colormap->entry != NULL)
4107 pfree(arg->colormap->entry);
4108 pfree(arg->colormap);
4109 }
4110
4111 if (arg->nentry) {
4112 for (i = 0; i < arg->nentry; i++) {
4113 if (arg->entry[i] != NULL)
4114 pfree(arg->entry[i]);
4115 }
4116 pfree(arg->entry);
4117 }
4118
4119 if (arg->nelement) {
4120 for (i = 0; i < arg->nelement; i++)
4121 pfree(arg->element[i]);
4122 pfree(arg->element);
4123 }
4124
4125 pfree(arg);
4126 arg = NULL;
4127 }
4128
4129 PG_FUNCTION_INFO_V1(RASTER_colorMap);
4130 Datum RASTER_colorMap(PG_FUNCTION_ARGS)
4131 {
4132 rt_pgraster *pgraster = NULL;
4133 rtpg_colormap_arg arg = NULL;
4134 char *junk = NULL;
4135 rt_raster raster = NULL;
4136
4137 POSTGIS_RT_DEBUG(3, "RASTER_colorMap: Starting");
4138
4139 /* pgraster is NULL, return NULL */
4140 if (PG_ARGISNULL(0))
4141 PG_RETURN_NULL();
4142
4143 /* init arg */
4144 arg = rtpg_colormap_arg_init();
4145 if (arg == NULL) {
4146 elog(ERROR, "RASTER_colorMap: Could not initialize argument structure");
4147 PG_RETURN_NULL();
4148 }
4149
4150 /* raster (0) */
4151 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4152
4153 /* raster object */
4154 arg->raster = rt_raster_deserialize(pgraster, FALSE);
4155 if (!arg->raster) {
4156 rtpg_colormap_arg_destroy(arg);
4157 PG_FREE_IF_COPY(pgraster, 0);
4158 elog(ERROR, "RASTER_colorMap: Could not deserialize raster");
4159 PG_RETURN_NULL();
4160 }
4161
4162 /* nband (1) */
4163 if (!PG_ARGISNULL(1))
4164 arg->nband = PG_GETARG_INT32(1);
4165 POSTGIS_RT_DEBUGF(4, "nband = %d", arg->nband);
4166
4167 /* check that band works */
4168 if (!rt_raster_has_band(arg->raster, arg->nband - 1)) {
4169 elog(NOTICE, "Raster does not have band at index %d. Returning empty raster", arg->nband);
4170
4171 raster = rt_raster_clone(arg->raster, 0);
4172 if (raster == NULL) {
4173 rtpg_colormap_arg_destroy(arg);
4174 PG_FREE_IF_COPY(pgraster, 0);
4175 elog(ERROR, "RASTER_colorMap: Could not create empty raster");
4176 PG_RETURN_NULL();
4177 }
4178
4179 rtpg_colormap_arg_destroy(arg);
4180 PG_FREE_IF_COPY(pgraster, 0);
4181
4182 pgraster = rt_raster_serialize(raster);
4183 rt_raster_destroy(raster);
4184 if (pgraster == NULL)
4185 PG_RETURN_NULL();
4186
4187 SET_VARSIZE(pgraster, ((rt_pgraster*) pgraster)->size);
4188 PG_RETURN_POINTER(pgraster);
4189 }
4190
4191 /* get band */
4192 arg->band = rt_raster_get_band(arg->raster, arg->nband - 1);
4193 if (arg->band == NULL) {
4194 int nband = arg->nband;
4195 rtpg_colormap_arg_destroy(arg);
4196 PG_FREE_IF_COPY(pgraster, 0);
4197 elog(ERROR, "RASTER_colorMap: Could not get band at index %d", nband);
4198 PG_RETURN_NULL();
4199 }
4200
4201 /* method (3) */
4202 if (!PG_ARGISNULL(3)) {
4203 char *method = NULL;
4204 char *tmp = text_to_cstring(PG_GETARG_TEXT_P(3));
4205 POSTGIS_RT_DEBUGF(4, "raw method = %s", tmp);
4206
4207 method = rtpg_trim(tmp);
4208 pfree(tmp);
4209 method = rtpg_strtoupper(method);
4210
4211 if (strcmp(method, "INTERPOLATE") == 0)
4212 arg->colormap->method = CM_INTERPOLATE;
4213 else if (strcmp(method, "EXACT") == 0)
4214 arg->colormap->method = CM_EXACT;
4215 else if (strcmp(method, "NEAREST") == 0)
4216 arg->colormap->method = CM_NEAREST;
4217 else {
4218 elog(NOTICE, "Unknown value provided for method. Defaulting to INTERPOLATE");
4219 arg->colormap->method = CM_INTERPOLATE;
4220 }
4221 }
4222 /* default to INTERPOLATE */
4223 else
4224 arg->colormap->method = CM_INTERPOLATE;
4225 POSTGIS_RT_DEBUGF(4, "method = %d", arg->colormap->method);
4226
4227 /* colormap (2) */
4228 if (PG_ARGISNULL(2)) {
4229 rtpg_colormap_arg_destroy(arg);
4230 PG_FREE_IF_COPY(pgraster, 0);
4231 elog(ERROR, "RASTER_colorMap: Value must be provided for colormap");
4232 PG_RETURN_NULL();
4233 }
4234 else {
4235 char *tmp = NULL;
4236 char *colormap = text_to_cstring(PG_GETARG_TEXT_P(2));
4237 char *_entry;
4238 char *_element;
4239 uint32_t i = 0;
4240 uint32_t j = 0;
4241
4242 POSTGIS_RT_DEBUGF(4, "colormap = %s", colormap);
4243
4244 /* empty string */
4245 if (!strlen(colormap)) {
4246 rtpg_colormap_arg_destroy(arg);
4247 PG_FREE_IF_COPY(pgraster, 0);
4248 elog(ERROR, "RASTER_colorMap: Value must be provided for colormap");
4249 PG_RETURN_NULL();
4250 }
4251
4252 arg->entry = rtpg_strsplit(colormap, "\n", &(arg->nentry));
4253 pfree(colormap);
4254 if (arg->nentry < 1) {
4255 rtpg_colormap_arg_destroy(arg);
4256 PG_FREE_IF_COPY(pgraster, 0);
4257 elog(ERROR, "RASTER_colorMap: Could not process the value provided for colormap");
4258 PG_RETURN_NULL();
4259 }
4260
4261 /* allocate the max # of colormap entries */
4262 arg->colormap->entry = palloc(sizeof(struct rt_colormap_entry_t) * arg->nentry);
4263 if (arg->colormap->entry == NULL) {
4264 rtpg_colormap_arg_destroy(arg);
4265 PG_FREE_IF_COPY(pgraster, 0);
4266 elog(ERROR, "RASTER_colorMap: Could not allocate memory for colormap entries");
4267 PG_RETURN_NULL();
4268 }
4269 memset(arg->colormap->entry, 0, sizeof(struct rt_colormap_entry_t) * arg->nentry);
4270
4271 /* each entry */
4272 for (i = 0; i < arg->nentry; i++) {
4273 /* substitute space for other delimiters */
4274 tmp = rtpg_strreplace(arg->entry[i], ":", " ", NULL);
4275 _entry = rtpg_strreplace(tmp, ",", " ", NULL);
4276 pfree(tmp);
4277 tmp = rtpg_strreplace(_entry, "\t", " ", NULL);
4278 pfree(_entry);
4279 _entry = rtpg_trim(tmp);
4280 pfree(tmp);
4281
4282 POSTGIS_RT_DEBUGF(4, "Processing entry[%d] = %s", i, arg->entry[i]);
4283 POSTGIS_RT_DEBUGF(4, "Cleaned entry[%d] = %s", i, _entry);
4284
4285 /* empty entry, continue */
4286 if (!strlen(_entry)) {
4287 POSTGIS_RT_DEBUGF(3, "Skipping empty entry[%d]", i);
4288 pfree(_entry);
4289 continue;
4290 }
4291
4292 arg->element = rtpg_strsplit(_entry, " ", &(arg->nelement));
4293 pfree(_entry);
4294 if (arg->nelement < 2) {
4295 rtpg_colormap_arg_destroy(arg);
4296 PG_FREE_IF_COPY(pgraster, 0);
4297 elog(ERROR, "RASTER_colorMap: Could not process colormap entry %d", i + 1);
4298 PG_RETURN_NULL();
4299 }
4300 else if (arg->nelement > 5) {
4301 elog(NOTICE, "More than five elements in colormap entry %d. Using at most five elements", i + 1);
4302 arg->nelement = 5;
4303 }
4304
4305 /* smallest # of colors */
4306 if (((int)arg->nelement - 1) < arg->colormap->ncolor)
4307 arg->colormap->ncolor = arg->nelement - 1;
4308
4309 /* each element of entry */
4310 for (j = 0; j < arg->nelement; j++) {
4311
4312 _element = rtpg_trim(arg->element[j]);
4313 _element = rtpg_strtoupper(_element);
4314 POSTGIS_RT_DEBUGF(4, "Processing entry[%d][%d] = %s", i, j, arg->element[j]);
4315 POSTGIS_RT_DEBUGF(4, "Cleaned entry[%d][%d] = %s", i, j, _element);
4316
4317 /* first element is ALWAYS a band value, percentage OR "nv" string */
4318 if (j == 0) {
4319 char *percent = NULL;
4320
4321 /* NODATA */
4322 if (
4323 strcmp(_element, "NV") == 0 ||
4324 strcmp(_element, "NULL") == 0 ||
4325 strcmp(_element, "NODATA") == 0
4326 ) {
4327 POSTGIS_RT_DEBUG(4, "Processing NODATA string");
4328
4329 if (arg->nodataentry > -1) {
4330 elog(NOTICE, "More than one NODATA entry found. Using only the first one");
4331 }
4332 else {
4333 arg->colormap->entry[arg->colormap->nentry].isnodata = 1;
4334 /* no need to set value as value comes from band's NODATA */
4335 arg->colormap->entry[arg->colormap->nentry].value = 0;
4336 }
4337 }
4338 /* percent value */
4339 else if ((percent = strchr(_element, '%')) != NULL) {
4340 double value;
4341 POSTGIS_RT_DEBUG(4, "Processing percent string");
4342
4343 /* get the band stats */
4344 if (arg->bandstats == NULL) {
4345 POSTGIS_RT_DEBUG(4, "Getting band stats");
4346
4347 arg->bandstats = rt_band_get_summary_stats(arg->band, 1, 1, 0, NULL, NULL, NULL);
4348 if (arg->bandstats == NULL) {
4349 pfree(_element);
4350 rtpg_colormap_arg_destroy(arg);
4351 PG_FREE_IF_COPY(pgraster, 0);
4352 elog(ERROR, "RASTER_colorMap: Could not get band's summary stats to process percentages");
4353 PG_RETURN_NULL();
4354 }
4355 }
4356
4357 /* get the string before the percent char */
4358 tmp = palloc(sizeof(char) * (percent - _element + 1));
4359 if (tmp == NULL) {
4360 pfree(_element);
4361 rtpg_colormap_arg_destroy(arg);
4362 PG_FREE_IF_COPY(pgraster, 0);
4363 elog(ERROR, "RASTER_colorMap: Could not allocate memory for value of percentage");
4364 PG_RETURN_NULL();
4365 }
4366
4367 memcpy(tmp, _element, percent - _element);
4368 tmp[percent - _element] = '\0';
4369 POSTGIS_RT_DEBUGF(4, "Percent value = %s", tmp);
4370
4371 /* get percentage value */
4372 errno = 0;
4373 value = strtod(tmp, NULL);
4374 pfree(tmp);
4375 if (errno != 0 || _element == junk) {
4376 pfree(_element);
4377 rtpg_colormap_arg_destroy(arg);
4378 PG_FREE_IF_COPY(pgraster, 0);
4379 elog(ERROR, "RASTER_colorMap: Could not process percent string to value");
4380 PG_RETURN_NULL();
4381 }
4382
4383 /* check percentage */
4384 if (value < 0.) {
4385 elog(NOTICE, "Percentage values cannot be less than zero. Defaulting to zero");
4386 value = 0.;
4387 }
4388 else if (value > 100.) {
4389 elog(NOTICE, "Percentage values cannot be greater than 100. Defaulting to 100");
4390 value = 100.;
4391 }
4392
4393 /* get the true pixel value */
4394 /* TODO: should the percentage be quantile based? */
4395 arg->colormap->entry[arg->colormap->nentry].value = ((value / 100.) * (arg->bandstats->max - arg->bandstats->min)) + arg->bandstats->min;
4396 }
4397 /* straight value */
4398 else {
4399 errno = 0;
4400 arg->colormap->entry[arg->colormap->nentry].value = strtod(_element, &junk);
4401 if (errno != 0 || _element == junk) {
4402 pfree(_element);
4403 rtpg_colormap_arg_destroy(arg);
4404 PG_FREE_IF_COPY(pgraster, 0);
4405 elog(ERROR, "RASTER_colorMap: Could not process string to value");
4406 PG_RETURN_NULL();
4407 }
4408 }
4409
4410 }
4411 /* RGB values (0 - 255) */
4412 else {
4413 int value = 0;
4414
4415 errno = 0;
4416 value = (int) strtod(_element, &junk);
4417 if (errno != 0 || _element == junk) {
4418 pfree(_element);
4419 rtpg_colormap_arg_destroy(arg);
4420 PG_FREE_IF_COPY(pgraster, 0);
4421 elog(ERROR, "RASTER_colorMap: Could not process string to value");
4422 PG_RETURN_NULL();
4423 }
4424
4425 if (value > 255) {
4426 elog(NOTICE, "RGBA value cannot be greater than 255. Defaulting to 255");
4427 value = 255;
4428 }
4429 else if (value < 0) {
4430 elog(NOTICE, "RGBA value cannot be less than zero. Defaulting to zero");
4431 value = 0;
4432 }
4433 arg->colormap->entry[arg->colormap->nentry].color[j - 1] = value;
4434 }
4435
4436 pfree(_element);
4437 }
4438
4439 POSTGIS_RT_DEBUGF(4, "colormap->entry[%d] (isnodata, value, R, G, B, A) = (%d, %f, %d, %d, %d, %d)",
4440 arg->colormap->nentry,
4441 arg->colormap->entry[arg->colormap->nentry].isnodata,
4442 arg->colormap->entry[arg->colormap->nentry].value,
4443 arg->colormap->entry[arg->colormap->nentry].color[0],
4444 arg->colormap->entry[arg->colormap->nentry].color[1],
4445 arg->colormap->entry[arg->colormap->nentry].color[2],
4446 arg->colormap->entry[arg->colormap->nentry].color[3]
4447 );
4448
4449 arg->colormap->nentry++;
4450 }
4451
4452 POSTGIS_RT_DEBUGF(4, "colormap->nentry = %d", arg->colormap->nentry);
4453 POSTGIS_RT_DEBUGF(4, "colormap->ncolor = %d", arg->colormap->ncolor);
4454 }
4455
4456 /* call colormap */
4457 raster = rt_raster_colormap(arg->raster, arg->nband - 1, arg->colormap);
4458 if (raster == NULL) {
4459 rtpg_colormap_arg_destroy(arg);
4460 PG_FREE_IF_COPY(pgraster, 0);
4461 elog(ERROR, "RASTER_colorMap: Could not create new raster with applied colormap");
4462 PG_RETURN_NULL();
4463 }
4464
4465 rtpg_colormap_arg_destroy(arg);
4466 PG_FREE_IF_COPY(pgraster, 0);
4467 pgraster = rt_raster_serialize(raster);
4468 rt_raster_destroy(raster);
4469
4470 POSTGIS_RT_DEBUG(3, "RASTER_colorMap: Done");
4471
4472 if (pgraster == NULL)
4473 PG_RETURN_NULL();
4474
4475 SET_VARSIZE(pgraster, ((rt_pgraster*) pgraster)->size);
4476 PG_RETURN_POINTER(pgraster);
4477 }
4478
4479 PG_FUNCTION_INFO_V1(RASTER_mapAlgebraExpr);
4480 Datum RASTER_mapAlgebraExpr(PG_FUNCTION_ARGS)
4481 {
4482 rt_pgraster *pgraster = NULL;
4483 rt_pgraster *pgrtn = NULL;
4484 rt_raster raster = NULL;
4485 rt_raster newrast = NULL;
4486 rt_band band = NULL;
4487 rt_band newband = NULL;
4488 int x, y, nband, width, height;
4489 double r;
4490 double newnodatavalue = 0.0;
4491 double newinitialvalue = 0.0;
4492 double newval = 0.0;
4493 char *newexpr = NULL;
4494 char *initexpr = NULL;
4495 char *expression = NULL;
4496 int hasnodataval = 0;
4497 double nodataval = 0.;
4498 rt_pixtype newpixeltype;
4499 int skipcomputation = 0;
4500 int len = 0;
4501 const int argkwcount = 3;
4502 enum KEYWORDS { kVAL=0, kX=1, kY=2 };
4503 char *argkw[] = {"[rast]", "[rast.x]", "[rast.y]"};
4504 Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
4505 int argcount = 0;
4506 Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
4507 uint8_t argpos[3] = {0};
4508 char place[12];
4509 int idx = 0;
4510 int ret = -1;
4511 TupleDesc tupdesc;
4512 SPIPlanPtr spi_plan = NULL;
4513 SPITupleTable * tuptable = NULL;
4514 HeapTuple tuple;
4515 char * strFromText = NULL;
4516 Datum *values = NULL;
4517 Datum datum = (Datum)NULL;
4518 char *nulls = NULL;
4519 bool isnull = FALSE;
4520 int i = 0;
4521 int j = 0;
4522
4523 POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraExpr: Starting...");
4524
4525 /* Check raster */
4526 if (PG_ARGISNULL(0)) {
4527 elog(NOTICE, "Raster is NULL. Returning NULL");
4528 PG_RETURN_NULL();
4529 }
4530
4531
4532 /* Deserialize raster */
4533 pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4534 raster = rt_raster_deserialize(pgraster, FALSE);
4535 if (NULL == raster) {
4536 PG_FREE_IF_COPY(pgraster, 0);
4537 elog(ERROR, "RASTER_mapAlgebraExpr: Could not deserialize raster");
4538 PG_RETURN_NULL();
4539 }
4540
4541 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting arguments...");
4542
4543 if (PG_ARGISNULL(1))
4544 nband = 1;
4545 else
4546 nband = PG_GETARG_INT32(1);
4547
4548 if (nband < 1)
4549 nband = 1;
4550
4551
4552 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new empty raster...");
4553
4554 /**
4555 * Create a new empty raster with having the same georeference as the
4556 * provided raster
4557 **/
4558 width = rt_raster_get_width(raster);
4559 height = rt_raster_get_height(raster);
4560
4561 newrast = rt_raster_new(width, height);
4562
4563 if ( NULL == newrast ) {
4564 PG_FREE_IF_COPY(pgraster, 0);
4565 elog(ERROR, "RASTER_mapAlgebraExpr: Could not create a new raster");
4566 PG_RETURN_NULL();
4567 }
4568
4569 rt_raster_set_scale(newrast,
4570 rt_raster_get_x_scale(raster),
4571 rt_raster_get_y_scale(raster));
4572
4573 rt_raster_set_offsets(newrast,
4574 rt_raster_get_x_offset(raster),
4575 rt_raster_get_y_offset(raster));
4576
4577 rt_raster_set_skews(newrast,
4578 rt_raster_get_x_skew(raster),
4579 rt_raster_get_y_skew(raster));
4580
4581 rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
4582
4583
4584 /**
4585 * If this new raster is empty (width = 0 OR height = 0) then there is
4586 * nothing to compute and we return it right now
4587 **/
4588 if (rt_raster_is_empty(newrast))
4589 {
4590 elog(NOTICE, "Raster is empty. Returning an empty raster");
4591 rt_raster_destroy(raster);
4592 PG_FREE_IF_COPY(pgraster, 0);
4593
4594 pgrtn = rt_raster_serialize(newrast);
4595 rt_raster_destroy(newrast);
4596 if (NULL == pgrtn) {
4597
4598 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4599 PG_RETURN_NULL();
4600 }
4601
4602 SET_VARSIZE(pgrtn, pgrtn->size);
4603 PG_RETURN_POINTER(pgrtn);
4604 }
4605
4606
4607 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Getting raster band %d...", nband);
4608
4609 /**
4610 * Check if the raster has the required band. Otherwise, return a raster
4611 * without band
4612 **/
4613 if (!rt_raster_has_band(raster, nband - 1)) {
4614 elog(NOTICE, "Raster does not have the required band. Returning a raster "
4615 "without a band");
4616 rt_raster_destroy(raster);
4617 PG_FREE_IF_COPY(pgraster, 0);
4618
4619 pgrtn = rt_raster_serialize(newrast);
4620 rt_raster_destroy(newrast);
4621 if (NULL == pgrtn) {
4622 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4623 PG_RETURN_NULL();
4624 }
4625
4626 SET_VARSIZE(pgrtn, pgrtn->size);
4627 PG_RETURN_POINTER(pgrtn);
4628 }
4629
4630 /* Get the raster band */
4631 band = rt_raster_get_band(raster, nband - 1);
4632 if ( NULL == band ) {
4633 elog(NOTICE, "Could not get the required band. Returning a raster "
4634 "without a band");
4635 rt_raster_destroy(raster);
4636 PG_FREE_IF_COPY(pgraster, 0);
4637
4638 pgrtn = rt_raster_serialize(newrast);
4639 rt_raster_destroy(newrast);
4640 if (NULL == pgrtn) {
4641 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4642 PG_RETURN_NULL();
4643 }
4644
4645 SET_VARSIZE(pgrtn, pgrtn->size);
4646 PG_RETURN_POINTER(pgrtn);
4647 }
4648
4649 /*
4650 * Get NODATA value
4651 */
4652 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting NODATA value for band...");
4653
4654 if (rt_band_get_hasnodata_flag(band)) {
4655 rt_band_get_nodata(band, &newnodatavalue);
4656 }
4657
4658 else {
4659 newnodatavalue = rt_band_get_min_value(band);
4660 }
4661
4662 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: NODATA value for band: = %f",
4663 newnodatavalue);
4664
4665 /**
4666 * We set the initial value of the future band to nodata value. If nodata
4667 * value is null, then the raster will be initialized to
4668 * rt_band_get_min_value but all the values should be recomputed anyway
4669 **/
4670 newinitialvalue = newnodatavalue;
4671
4672 /**
4673 * Set the new pixeltype
4674 **/
4675 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Setting pixeltype...");
4676
4677 if (PG_ARGISNULL(2)) {
4678 newpixeltype = rt_band_get_pixtype(band);
4679 }
4680
4681 else {
4682 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
4683 newpixeltype = rt_pixtype_index_from_name(strFromText);
4684 pfree(strFromText);
4685 if (newpixeltype == PT_END)
4686 newpixeltype = rt_band_get_pixtype(band);
4687 }
4688
4689 if (newpixeltype == PT_END) {
4690 PG_FREE_IF_COPY(pgraster, 0);
4691 elog(ERROR, "RASTER_mapAlgebraExpr: Invalid pixeltype");
4692 PG_RETURN_NULL();
4693 }
4694
4695 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Pixeltype set to %s",
4696 rt_pixtype_name(newpixeltype));
4697
4698
4699 /* Construct expression for raster values */
4700 if (!PG_ARGISNULL(3)) {
4701 expression = text_to_cstring(PG_GETARG_TEXT_P(3));
4702 len = strlen("SELECT (") + strlen(expression) + strlen(")::double precision");
4703 initexpr = (char *)palloc(len + 1);
4704
4705 memcpy(initexpr, "SELECT (", strlen("SELECT ("));
4706 memcpy(initexpr + strlen("SELECT ("), expression, strlen(expression));
4707 memcpy(initexpr + strlen("SELECT (") + strlen(expression), ")::double precision", strlen(")::double precision"));
4708 initexpr[len] = '\0';
4709
4710 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression is %s", initexpr);
4711
4712 /* We don't need this memory */
4713 /*
4714 pfree(expression);
4715 expression = NULL;
4716 */
4717 }
4718
4719
4720
4721 /**
4722 * Optimization: If a nodataval is provided, use it for newinitialvalue.
4723 * Then, we can initialize the raster with this value and skip the
4724 * computation of nodata values one by one in the main computing loop
4725 **/
4726 if (!PG_ARGISNULL(4)) {
4727 hasnodataval = 1;
4728 nodataval = PG_GETARG_FLOAT8(4);
4729 newinitialvalue = nodataval;
4730
4731 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new initial value = %f",
4732 newinitialvalue);
4733 }
4734 else
4735 hasnodataval = 0;
4736
4737
4738
4739 /**
4740 * Optimization: If the raster is only filled with nodata values return
4741 * right now a raster filled with the newinitialvalue
4742 * TODO: Call rt_band_check_isnodata instead?
4743 **/
4744 if (rt_band_get_isnodata_flag(band)) {
4745
4746 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Band is a nodata band, returning "
4747 "a raster filled with nodata");
4748
4749 ret = rt_raster_generate_new_band(newrast, newpixeltype,
4750 newinitialvalue, TRUE, newnodatavalue, 0);
4751
4752 /* Free memory */
4753 if (initexpr)
4754 pfree(initexpr);
4755 rt_raster_destroy(raster);
4756 PG_FREE_IF_COPY(pgraster, 0);
4757
4758 /* Serialize created raster */
4759 pgrtn = rt_raster_serialize(newrast);
4760 rt_raster_destroy(newrast);
4761 if (NULL == pgrtn) {
4762 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4763 PG_RETURN_NULL();
4764 }
4765
4766 SET_VARSIZE(pgrtn, pgrtn->size);
4767 PG_RETURN_POINTER(pgrtn);
4768 }
4769
4770
4771 /**
4772 * Optimization: If expression resume to 'RAST' and hasnodataval is zero,
4773 * we can just return the band from the original raster
4774 **/
4775 if (initexpr != NULL && ( !strcmp(initexpr, "SELECT [rast]") || !strcmp(initexpr, "SELECT [rast.val]") ) && !hasnodataval) {
4776
4777 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression resumes to RAST. "
4778 "Returning raster with band %d from original raster", nband);
4779
4780 POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster has %d bands",
4781 rt_raster_get_num_bands(newrast));
4782
4783 rt_raster_copy_band(newrast, raster, nband - 1, 0);
4784
4785 POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster now has %d bands",
4786 rt_raster_get_num_bands(newrast));
4787
4788 pfree(initexpr);
4789 rt_raster_destroy(raster);
4790 PG_FREE_IF_COPY(pgraster, 0);
4791
4792 /* Serialize created raster */
4793 pgrtn = rt_raster_serialize(newrast);
4794 rt_raster_destroy(newrast);
4795 if (NULL == pgrtn) {
4796 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4797 PG_RETURN_NULL();
4798 }
4799
4800 SET_VARSIZE(pgrtn, pgrtn->size);
4801 PG_RETURN_POINTER(pgrtn);
4802 }
4803
4804 /**
4805 * Optimization: If expression resume to a constant (it does not contain
4806 * [rast)
4807 **/
4808 if (initexpr != NULL && strstr(initexpr, "[rast") == NULL) {
4809 ret = SPI_connect();
4810 if (ret != SPI_OK_CONNECT) {
4811 PG_FREE_IF_COPY(pgraster, 0);
4812 elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4813 PG_RETURN_NULL();
4814 };
4815
4816 /* Execute the expresion into newval */
4817 ret = SPI_execute(initexpr, FALSE, 0);
4818
4819 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
4820
4821 /* Free memory allocated out of the current context */
4822 if (SPI_tuptable)
4823 SPI_freetuptable(tuptable);
4824 PG_FREE_IF_COPY(pgraster, 0);
4825
4826 SPI_finish();
4827 elog(ERROR, "RASTER_mapAlgebraExpr: Invalid construction for expression");
4828 PG_RETURN_NULL();
4829 }
4830
4831 tupdesc = SPI_tuptable->tupdesc;
4832 tuptable = SPI_tuptable;
4833
4834 tuple = tuptable->vals[0];
4835 newexpr = SPI_getvalue(tuple, tupdesc, 1);
4836 if ( ! newexpr ) {
4837 POSTGIS_RT_DEBUG(3, "Constant expression evaluated to NULL, keeping initvalue");
4838 newval = newinitialvalue;
4839 } else {
4840 newval = atof(newexpr);
4841 }
4842
4843 SPI_freetuptable(tuptable);
4844
4845 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: New raster value = %f",
4846 newval);
4847
4848 SPI_finish();
4849
4850 skipcomputation = 1;
4851
4852 /**
4853 * Compute the new value, set it and we will return after creating the
4854 * new raster
4855 **/
4856 if (!hasnodataval) {
4857 newinitialvalue = newval;
4858 skipcomputation = 2;
4859 }
4860
4861 /* Return the new raster as it will be before computing pixel by pixel */
4862 else if (FLT_NEQ(newval, newinitialvalue)) {
4863 skipcomputation = 2;
4864 }
4865 }
4866
4867 /**
4868 * Create the raster receiving all the computed values. Initialize it to the
4869 * new initial value
4870 **/
4871 ret = rt_raster_generate_new_band(newrast, newpixeltype,
4872 newinitialvalue, TRUE, newnodatavalue, 0);
4873
4874 /**
4875 * Optimization: If expression is NULL, or all the pixels could be set in
4876 * one step, return the initialized raster now
4877 **/
4878 /*if (initexpr == NULL || skipcomputation == 2) {*/
4879 if (expression == NULL || skipcomputation == 2) {
4880
4881 /* Free memory */
4882 if (initexpr)
4883 pfree(initexpr);
4884 rt_raster_destroy(raster);
4885 PG_FREE_IF_COPY(pgraster, 0);
4886
4887 /* Serialize created raster */
4888 pgrtn = rt_raster_serialize(newrast);
4889 rt_raster_destroy(newrast);
4890 if (NULL == pgrtn) {
4891 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4892 PG_RETURN_NULL();
4893 }
4894
4895 SET_VARSIZE(pgrtn, pgrtn->size);
4896 PG_RETURN_POINTER(pgrtn);
4897 }
4898
4899 RASTER_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new raster band...");
4900
4901 /* Get the new raster band */
4902 newband = rt_raster_get_band(newrast, 0);
4903 if ( NULL == newband ) {
4904 elog(NOTICE, "Could not modify band for new raster. Returning new "
4905 "raster with the original band");
4906
4907 if (initexpr)
4908 pfree(initexpr);
4909 rt_raster_destroy(raster);
4910 PG_FREE_IF_COPY(pgraster, 0);
4911
4912 /* Serialize created raster */
4913 pgrtn = rt_raster_serialize(newrast);
4914 rt_raster_destroy(newrast);
4915 if (NULL == pgrtn) {
4916 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4917 PG_RETURN_NULL();
4918 }
4919
4920 SET_VARSIZE(pgrtn, pgrtn->size);
4921 PG_RETURN_POINTER(pgrtn);
4922 }
4923
4924 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Main computing loop (%d x %d)",
4925 width, height);
4926
4927 if (initexpr != NULL) {
4928 /* Convert [rast.val] to [rast] */
4929 newexpr = rtpg_strreplace(initexpr, "[rast.val]", "[rast]", NULL);
4930 pfree(initexpr); initexpr=newexpr;
4931
4932 sprintf(place,"$1");
4933 for (i = 0, j = 1; i < argkwcount; i++) {
4934 len = 0;
4935 newexpr = rtpg_strreplace(initexpr, argkw[i], place, &len);
4936 pfree(initexpr); initexpr=newexpr;
4937 if (len > 0) {
4938 argtype[argcount] = argkwtypes[i];
4939 argcount++;
4940 argpos[i] = j++;
4941
4942 sprintf(place, "$%d", j);
4943 }
4944 else {
4945 argpos[i] = 0;
4946 }
4947 }
4948
4949 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: initexpr = %s", initexpr);
4950
4951 /* define values */
4952 values = (Datum *) palloc(sizeof(Datum) * argcount);
4953 if (values == NULL) {
4954
4955 SPI_finish();
4956
4957 rt_raster_destroy(raster);
4958 PG_FREE_IF_COPY(pgraster, 0);
4959 rt_raster_destroy(newrast);
4960
4961 elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
4962 PG_RETURN_NULL();
4963 }
4964
4965 /* define nulls */
4966 nulls = (char *)palloc(argcount);
4967 if (nulls == NULL) {
4968
4969 SPI_finish();
4970
4971 rt_raster_destroy(raster);
4972 PG_FREE_IF_COPY(pgraster, 0);
4973 rt_raster_destroy(newrast);
4974
4975 elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
4976 PG_RETURN_NULL();
4977 }
4978
4979 /* Connect to SPI and prepare the expression */
4980 ret = SPI_connect();
4981 if (ret != SPI_OK_CONNECT) {
4982
4983 if (initexpr)
4984 pfree(initexpr);
4985 rt_raster_destroy(raster);
4986 PG_FREE_IF_COPY(pgraster, 0);
4987 rt_raster_destroy(newrast);
4988
4989 elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4990 PG_RETURN_NULL();
4991 };
4992
4993 /* Type of all arguments is FLOAT8OID */
4994 spi_plan = SPI_prepare(initexpr, argcount, argtype);
4995
4996 if (spi_plan == NULL) {
4997
4998 rt_raster_destroy(raster);
4999 PG_FREE_IF_COPY(pgraster, 0);
5000 rt_raster_destroy(newrast);
5001
5002 SPI_finish();
5003
5004 pfree(initexpr);
5005
5006 elog(ERROR, "RASTER_mapAlgebraExpr: Could not prepare expression");
5007 PG_RETURN_NULL();
5008 }
5009 }
5010
5011 for (x = 0; x < width; x++) {
5012 for(y = 0; y < height; y++) {
5013 ret = rt_band_get_pixel(band, x, y, &r, NULL);
5014
5015 /**
5016 * We compute a value only for the withdata value pixel since the
5017 * nodata value has already been set by the first optimization
5018 **/
5019 if (ret == ES_NONE && FLT_NEQ(r, newnodatavalue)) {
5020 if (skipcomputation == 0) {
5021 if (initexpr != NULL) {
5022 /* Reset the null arg flags. */
5023 memset(nulls, 'n', argcount);
5024
5025 for (i = 0; i < argkwcount; i++) {
5026 idx = argpos[i];
5027 if (idx < 1) continue;
5028 idx--;
5029
5030 if (i == kX ) {
5031 /* x is 0 based index, but SQL expects 1 based index */
5032 values[idx] = Int32GetDatum(x+1);
5033 nulls[idx] = ' ';
5034 }
5035 else if (i == kY) {
5036 /* y is 0 based index, but SQL expects 1 based index */
5037 values[idx] = Int32GetDatum(y+1);
5038 nulls[idx] = ' ';
5039 }
5040 else if (i == kVAL ) {
5041 values[idx] = Float8GetDatum(r);
5042 nulls[idx] = ' ';
5043 }
5044
5045 }
5046
5047 ret = SPI_execute_plan(spi_plan, values, nulls, FALSE, 0);
5048 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
5049 SPI_processed != 1) {
5050 if (SPI_tuptable)
5051 SPI_freetuptable(tuptable);
5052
5053 SPI_freeplan(spi_plan);
5054 SPI_finish();
5055
5056 pfree(values);
5057 pfree(nulls);
5058 pfree(initexpr);
5059
5060 rt_raster_destroy(raster);
5061 PG_FREE_IF_COPY(pgraster, 0);
5062 rt_raster_destroy(newrast);
5063
5064 elog(ERROR, "RASTER_mapAlgebraExpr: Error executing prepared plan");
5065
5066 PG_RETURN_NULL();
5067 }
5068
5069 tupdesc = SPI_tuptable->tupdesc;
5070 tuptable = SPI_tuptable;
5071
5072 tuple = tuptable->vals[0];
5073 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
5074 if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
5075 POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) errored, skip setting", x+1,y+1,r);
5076 newval = newinitialvalue;
5077 }
5078 else if ( isnull ) {
5079 POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting", x+1,y+1,r);
5080 newval = newinitialvalue;
5081 } else {
5082 newval = DatumGetFloat8(datum);
5083 }
5084
5085 SPI_freetuptable(tuptable);
5086 }
5087
5088 else
5089 newval = newinitialvalue;
5090
5091 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new value = %f",
5092 newval);
5093 }
5094
5095
5096 rt_band_set_pixel(newband, x, y, newval, NULL);
5097 }
5098
5099 }
5100 }
5101
5102 if (initexpr != NULL) {
5103 SPI_freeplan(spi_plan);
5104 SPI_finish();
5105
5106 pfree(values);
5107 pfree(nulls);
5108 pfree(initexpr);
5109 }
5110 else {
5111 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: no SPI cleanup");
5112 }
5113
5114
5115 /* The newrast band has been modified */
5116
5117 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster modified, serializing it.");
5118 /* Serialize created raster */
5119
5120 rt_raster_destroy(raster);
5121 PG_FREE_IF_COPY(pgraster, 0);
5122
5123 pgrtn = rt_raster_serialize(newrast);
5124 rt_raster_destroy(newrast);
5125 if (NULL == pgrtn)
5126 PG_RETURN_NULL();
5127
5128 SET_VARSIZE(pgrtn, pgrtn->size);
5129
5130 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster serialized");
5131
5132
5133 POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraExpr: returning raster");
5134
5135
5136 PG_RETURN_POINTER(pgrtn);
5137 }
5138
5139 /*
5140 * One raster user function MapAlgebra.
5141 */
5142 PG_FUNCTION_INFO_V1(RASTER_mapAlgebraFct);
5143 Datum RASTER_mapAlgebraFct(PG_FUNCTION_ARGS)
5144 {
5145 rt_pgraster *pgraster = NULL;
5146 rt_pgraster *pgrtn = NULL;
5147 rt_raster raster = NULL;
5148 rt_raster newrast = NULL;
5149 rt_band band = NULL;
5150 rt_band newband = NULL;
5151 int x, y, nband, width, height;
5152 double r;
5153 double newnodatavalue = 0.0;
5154 double newinitialvalue = 0.0;
5155 double newval = 0.0;
5156 rt_pixtype newpixeltype;
5157 int ret = -1;
5158 Oid oid;
5159 FmgrInfo cbinfo;
5160 #if POSTGIS_PGSQL_VERSION < 120
5161 FunctionCallInfoData cbdata;
5162 #else
5163 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS); /* Could be optimized */
5164 #endif
5165 Datum tmpnewval;
5166 char * strFromText = NULL;
5167 int k = 0;
5168
5169 POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFct: STARTING...");
5170
5171 /* Check raster */
5172 if (PG_ARGISNULL(0)) {
5173 elog(WARNING, "Raster is NULL. Returning NULL");
5174 PG_RETURN_NULL();
5175 }
5176
5177
5178 /* Deserialize raster */
5179 pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5180 raster = rt_raster_deserialize(pgraster, FALSE);
5181 if (NULL == raster) {
5182 PG_FREE_IF_COPY(pgraster, 0);
5183 elog(ERROR, "RASTER_mapAlgebraFct: Could not deserialize raster");
5184 PG_RETURN_NULL();
5185 }
5186
5187 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting arguments...");
5188
5189 /* Get the rest of the arguments */
5190
5191 if (PG_ARGISNULL(1))
5192 nband = 1;
5193 else
5194 nband = PG_GETARG_INT32(1);
5195
5196 if (nband < 1)
5197 nband = 1;
5198
5199 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Creating new empty raster...");
5200
5201 /**
5202 * Create a new empty raster with having the same georeference as the
5203 * provided raster
5204 **/
5205 width = rt_raster_get_width(raster);
5206 height = rt_raster_get_height(raster);
5207
5208 newrast = rt_raster_new(width, height);
5209
5210 if ( NULL == newrast ) {
5211
5212 rt_raster_destroy(raster);
5213 PG_FREE_IF_COPY(pgraster, 0);
5214
5215 elog(ERROR, "RASTER_mapAlgebraFct: Could not create a new raster");
5216 PG_RETURN_NULL();
5217 }
5218
5219 rt_raster_set_scale(newrast,
5220 rt_raster_get_x_scale(raster),
5221 rt_raster_get_y_scale(raster));
5222
5223 rt_raster_set_offsets(newrast,
5224 rt_raster_get_x_offset(raster),
5225 rt_raster_get_y_offset(raster));
5226
5227 rt_raster_set_skews(newrast,
5228 rt_raster_get_x_skew(raster),
5229 rt_raster_get_y_skew(raster));
5230
5231 rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
5232
5233
5234 /**
5235 * If this new raster is empty (width = 0 OR height = 0) then there is
5236 * nothing to compute and we return it right now
5237 **/
5238 if (rt_raster_is_empty(newrast))
5239 {
5240 elog(NOTICE, "Raster is empty. Returning an empty raster");
5241 rt_raster_destroy(raster);
5242 PG_FREE_IF_COPY(pgraster, 0);
5243
5244 pgrtn = rt_raster_serialize(newrast);
5245 rt_raster_destroy(newrast);
5246 if (NULL == pgrtn) {
5247 elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5248 PG_RETURN_NULL();
5249 }
5250
5251 SET_VARSIZE(pgrtn, pgrtn->size);
5252 PG_RETURN_POINTER(pgrtn);
5253 }
5254
5255 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Getting raster band %d...", nband);
5256
5257 /**
5258 * Check if the raster has the required band. Otherwise, return a raster
5259 * without band
5260 **/
5261 if (!rt_raster_has_band(raster, nband - 1)) {
5262 elog(NOTICE, "Raster does not have the required band. Returning a raster "
5263 "without a band");
5264 rt_raster_destroy(raster);
5265 PG_FREE_IF_COPY(pgraster, 0);
5266
5267 pgrtn = rt_raster_serialize(newrast);
5268 rt_raster_destroy(newrast);
5269 if (NULL == pgrtn) {
5270 elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5271 PG_RETURN_NULL();
5272 }
5273
5274 SET_VARSIZE(pgrtn, pgrtn->size);
5275 PG_RETURN_POINTER(pgrtn);
5276 }
5277
5278 /* Get the raster band */
5279 band = rt_raster_get_band(raster, nband - 1);
5280 if ( NULL == band ) {
5281 elog(NOTICE, "Could not get the required band. Returning a raster "
5282 "without a band");
5283 rt_raster_destroy(raster);
5284 PG_FREE_IF_COPY(pgraster, 0);
5285
5286 pgrtn = rt_raster_serialize(newrast);
5287 rt_raster_destroy(newrast);
5288 if (NULL == pgrtn) {
5289 elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5290 PG_RETURN_NULL();
5291 }
5292
5293 SET_VARSIZE(pgrtn, pgrtn->size);
5294 PG_RETURN_POINTER(pgrtn);
5295 }
5296
5297 /*
5298 * Get NODATA value
5299 */
5300 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting NODATA value for band...");
5301
5302 if (rt_band_get_hasnodata_flag(band)) {
5303 rt_band_get_nodata(band, &newnodatavalue);
5304 }
5305
5306 else {
5307 newnodatavalue = rt_band_get_min_value(band);
5308 }
5309
5310 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: NODATA value for band: %f",
5311 newnodatavalue);
5312 /**
5313 * We set the initial value of the future band to nodata value. If nodata
5314 * value is null, then the raster will be initialized to
5315 * rt_band_get_min_value but all the values should be recomputed anyway
5316 **/
5317 newinitialvalue = newnodatavalue;
5318
5319 /**
5320 * Set the new pixeltype
5321 **/
5322 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Setting pixeltype...");
5323
5324 if (PG_ARGISNULL(2)) {
5325 newpixeltype = rt_band_get_pixtype(band);
5326 }
5327
5328 else {
5329 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5330 newpixeltype = rt_pixtype_index_from_name(strFromText);
5331 pfree(strFromText);
5332 if (newpixeltype == PT_END)
5333 newpixeltype = rt_band_get_pixtype(band);
5334 }
5335
5336 if (newpixeltype == PT_END) {
5337
5338 rt_raster_destroy(raster);
5339 PG_FREE_IF_COPY(pgraster, 0);
5340 rt_raster_destroy(newrast);
5341
5342 elog(ERROR, "RASTER_mapAlgebraFct: Invalid pixeltype");
5343 PG_RETURN_NULL();
5344 }
5345
5346 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Pixeltype set to %s",
5347 rt_pixtype_name(newpixeltype));
5348
5349 /* Get the name of the callback user function for raster values */
5350 if (PG_ARGISNULL(3)) {
5351
5352 rt_raster_destroy(raster);
5353 PG_FREE_IF_COPY(pgraster, 0);
5354 rt_raster_destroy(newrast);
5355
5356 elog(ERROR, "RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
5357 PG_RETURN_NULL();
5358 }
5359
5360 oid = PG_GETARG_OID(3);
5361 if (oid == InvalidOid) {
5362
5363 rt_raster_destroy(raster);
5364 PG_FREE_IF_COPY(pgraster, 0);
5365 rt_raster_destroy(newrast);
5366
5367 elog(ERROR, "RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
5368 PG_RETURN_NULL();
5369 }
5370
5371 fmgr_info(oid, &cbinfo);
5372
5373 /* function cannot return set */
5374 if (cbinfo.fn_retset) {
5375
5376 rt_raster_destroy(raster);
5377 PG_FREE_IF_COPY(pgraster, 0);
5378 rt_raster_destroy(newrast);
5379
5380 elog(ERROR, "RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
5381 PG_RETURN_NULL();
5382 }
5383 /* function should have correct # of args */
5384 else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {
5385
5386 rt_raster_destroy(raster);
5387 PG_FREE_IF_COPY(pgraster, 0);
5388 rt_raster_destroy(newrast);
5389
5390 elog(ERROR, "RASTER_mapAlgebraFct: Function does not have two or three input parameters");
5391 PG_RETURN_NULL();
5392 }
5393
5394 if (cbinfo.fn_nargs == 2)
5395 k = 1;
5396 else
5397 k = 2;
5398
5399 if (func_volatile(oid) == 'v') {
5400 elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5401 }
5402
5403 /* prep function call data */
5404 #if POSTGIS_PGSQL_VERSION < 120
5405 InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5406
5407 cbdata.argnull[0] = FALSE;
5408 cbdata.argnull[1] = FALSE;
5409 cbdata.argnull[2] = FALSE;
5410 #else
5411 InitFunctionCallInfoData(*cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5412
5413 cbdata->args[0].isnull = FALSE;
5414 cbdata->args[1].isnull = FALSE;
5415 cbdata->args[2].isnull = FALSE;
5416 #endif
5417
5418 /* check that the function isn't strict if the args are null. */
5419 if (PG_ARGISNULL(4)) {
5420 if (cbinfo.fn_strict) {
5421
5422 rt_raster_destroy(raster);
5423 PG_FREE_IF_COPY(pgraster, 0);
5424 rt_raster_destroy(newrast);
5425
5426 elog(ERROR, "RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
5427 PG_RETURN_NULL();
5428 }
5429
5430 #if POSTGIS_PGSQL_VERSION < 120
5431 cbdata.arg[k] = (Datum)NULL;
5432 cbdata.argnull[k] = TRUE;
5433 #else
5434 cbdata->args[k].value = (Datum)NULL;
5435 cbdata->args[k].isnull = TRUE;
5436 #endif
5437 }
5438 else {
5439 #if POSTGIS_PGSQL_VERSION < 120
5440 cbdata.arg[k] = PG_GETARG_DATUM(4);
5441 #else
5442 cbdata->args[k].value = PG_GETARG_DATUM(4);
5443 #endif
5444 }
5445
5446 /**
5447 * Optimization: If the raster is only filled with nodata values return
5448 * right now a raster filled with the nodatavalueexpr
5449 * TODO: Call rt_band_check_isnodata instead?
5450 **/
5451 if (rt_band_get_isnodata_flag(band)) {
5452
5453 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Band is a nodata band, returning "
5454 "a raster filled with nodata");
5455
5456 ret = rt_raster_generate_new_band(newrast, newpixeltype,
5457 newinitialvalue, TRUE, newnodatavalue, 0);
5458
5459 rt_raster_destroy(raster);
5460 PG_FREE_IF_COPY(pgraster, 0);
5461
5462 /* Serialize created raster */
5463 pgrtn = rt_raster_serialize(newrast);
5464 rt_raster_destroy(newrast);
5465 if (NULL == pgrtn) {
5466 elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5467 PG_RETURN_NULL();
5468 }
5469
5470 SET_VARSIZE(pgrtn, pgrtn->size);
5471 PG_RETURN_POINTER(pgrtn);
5472 }
5473
5474
5475 /**
5476 * Create the raster receiving all the computed values. Initialize it to the
5477 * new initial value
5478 **/
5479 ret = rt_raster_generate_new_band(newrast, newpixeltype,
5480 newinitialvalue, TRUE, newnodatavalue, 0);
5481
5482 /* Get the new raster band */
5483 newband = rt_raster_get_band(newrast, 0);
5484 if ( NULL == newband ) {
5485 elog(NOTICE, "Could not modify band for new raster. Returning new "
5486 "raster with the original band");
5487
5488 rt_raster_destroy(raster);
5489 PG_FREE_IF_COPY(pgraster, 0);
5490
5491 /* Serialize created raster */
5492 pgrtn = rt_raster_serialize(newrast);
5493 rt_raster_destroy(newrast);
5494 if (NULL == pgrtn) {
5495 elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5496 PG_RETURN_NULL();
5497 }
5498
5499 SET_VARSIZE(pgrtn, pgrtn->size);
5500 PG_RETURN_POINTER(pgrtn);
5501 }
5502
5503 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Main computing loop (%d x %d)",
5504 width, height);
5505
5506 for (x = 0; x < width; x++) {
5507 for(y = 0; y < height; y++) {
5508 ret = rt_band_get_pixel(band, x, y, &r, NULL);
5509
5510 /**
5511 * We compute a value only for the withdata value pixel since the
5512 * nodata value has already been set by the first optimization
5513 **/
5514 if (ret == ES_NONE) {
5515 if (FLT_EQ(r, newnodatavalue)) {
5516 if (cbinfo.fn_strict) {
5517 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
5518 continue;
5519 }
5520 #if POSTGIS_PGSQL_VERSION < 120
5521 cbdata.argnull[0] = TRUE;
5522 cbdata.arg[0] = (Datum)NULL;
5523 #else
5524 cbdata->args[0].isnull = TRUE;
5525 cbdata->args[0].value = (Datum)NULL;
5526 #endif
5527 }
5528 else {
5529 #if POSTGIS_PGSQL_VERSION < 120
5530 cbdata.argnull[0] = FALSE;
5531 cbdata.arg[0] = Float8GetDatum(r);
5532 #else
5533 cbdata->args[0].isnull = FALSE;
5534 cbdata->args[0].value = Float8GetDatum(r);
5535 #endif
5536 }
5537
5538 /* Add pixel positions if callback has proper # of args */
5539 if (cbinfo.fn_nargs == 3) {
5540 Datum d[2];
5541 ArrayType *a;
5542
5543 d[0] = Int32GetDatum(x+1);
5544 d[1] = Int32GetDatum(y+1);
5545
5546 a = construct_array(d, 2, INT4OID, sizeof(int32), true, 'i');
5547
5548 #if POSTGIS_PGSQL_VERSION < 120
5549 cbdata.argnull[1] = FALSE;
5550 cbdata.arg[1] = PointerGetDatum(a);
5551 #else
5552 cbdata->args[1].isnull = FALSE;
5553 cbdata->args[1].value = PointerGetDatum(a);
5554 #endif
5555 }
5556
5557 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: (%dx%d), r = %f",
5558 x, y, r);
5559
5560 #if POSTGIS_PGSQL_VERSION < 120
5561 tmpnewval = FunctionCallInvoke(&cbdata);
5562
5563 if (cbdata.isnull) {
5564 newval = newnodatavalue;
5565 }
5566 #else
5567 tmpnewval = FunctionCallInvoke(cbdata);
5568
5569 if (cbdata->isnull)
5570 {
5571 newval = newnodatavalue;
5572 }
5573 #endif
5574 else {
5575 newval = DatumGetFloat8(tmpnewval);
5576 }
5577
5578 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: new value = %f",
5579 newval);
5580
5581 rt_band_set_pixel(newband, x, y, newval, NULL);
5582 }
5583
5584 }
5585 }
5586
5587 /* The newrast band has been modified */
5588
5589 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster modified, serializing it.");
5590 /* Serialize created raster */
5591
5592 rt_raster_destroy(raster);
5593 PG_FREE_IF_COPY(pgraster, 0);
5594
5595 pgrtn = rt_raster_serialize(newrast);
5596 rt_raster_destroy(newrast);
5597 if (NULL == pgrtn)
5598 PG_RETURN_NULL();
5599
5600 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster serialized");
5601
5602 POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFct: returning raster");
5603
5604 SET_VARSIZE(pgrtn, pgrtn->size);
5605 PG_RETURN_POINTER(pgrtn);
5606 }
5607
5608 /**
5609 * One raster neighborhood MapAlgebra
5610 */
5611 PG_FUNCTION_INFO_V1(RASTER_mapAlgebraFctNgb);
5612 Datum RASTER_mapAlgebraFctNgb(PG_FUNCTION_ARGS)
5613 {
5614 rt_pgraster *pgraster = NULL;
5615 rt_pgraster *pgrtn = NULL;
5616 rt_raster raster = NULL;
5617 rt_raster newrast = NULL;
5618 rt_band band = NULL;
5619 rt_band newband = NULL;
5620 int x, y, nband, width, height, ngbwidth, ngbheight, winwidth, winheight, u, v, nIndex, nNullItems;
5621 double r, rpix;
5622 double newnodatavalue = 0.0;
5623 double newinitialvalue = 0.0;
5624 double newval = 0.0;
5625 rt_pixtype newpixeltype;
5626 int ret = -1;
5627 Oid oid;
5628 FmgrInfo cbinfo;
5629 #if POSTGIS_PGSQL_VERSION < 120
5630 FunctionCallInfoData cbdata;
5631 #else
5632 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS); /* Could be optimized */
5633 #endif
5634 Datum tmpnewval;
5635 ArrayType * neighborDatum;
5636 char * strFromText = NULL;
5637 text * txtNodataMode = NULL;
5638 text * txtCallbackParam = NULL;
5639 int intReplace = 0;
5640 float fltReplace = 0;
5641 bool valuereplace = false, pixelreplace, nNodataOnly = true, nNullSkip = false;
5642 Datum * neighborData = NULL;
5643 bool * neighborNulls = NULL;
5644 int neighborDims[2];
5645 int neighborLbs[2];
5646 int16 typlen;
5647 bool typbyval;
5648 char typalign;
5649
5650 POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFctNgb: STARTING...");
5651
5652 /* Check raster */
5653 if (PG_ARGISNULL(0)) {
5654 elog(WARNING, "Raster is NULL. Returning NULL");
5655 PG_RETURN_NULL();
5656 }
5657
5658
5659 /* Deserialize raster */
5660 pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5661 raster = rt_raster_deserialize(pgraster, FALSE);
5662 if (NULL == raster)
5663 {
5664 PG_FREE_IF_COPY(pgraster, 0);
5665 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not deserialize raster");
5666 PG_RETURN_NULL();
5667 }
5668
5669 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Getting arguments...");
5670
5671 /* Get the rest of the arguments */
5672
5673 if (PG_ARGISNULL(1))
5674 nband = 1;
5675 else
5676 nband = PG_GETARG_INT32(1);
5677
5678 if (nband < 1)
5679 nband = 1;
5680
5681 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Creating new empty raster...");
5682
5683 /**
5684 * Create a new empty raster with having the same georeference as the
5685 * provided raster
5686 **/
5687 width = rt_raster_get_width(raster);
5688 height = rt_raster_get_height(raster);
5689
5690 newrast = rt_raster_new(width, height);
5691
5692 if ( NULL == newrast ) {
5693 rt_raster_destroy(raster);
5694 PG_FREE_IF_COPY(pgraster, 0);
5695 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not create a new raster");
5696 PG_RETURN_NULL();
5697 }
5698
5699 rt_raster_set_scale(newrast,
5700 rt_raster_get_x_scale(raster),
5701 rt_raster_get_y_scale(raster));
5702
5703 rt_raster_set_offsets(newrast,
5704 rt_raster_get_x_offset(raster),
5705 rt_raster_get_y_offset(raster));
5706
5707 rt_raster_set_skews(newrast,
5708 rt_raster_get_x_skew(raster),
5709 rt_raster_get_y_skew(raster));
5710
5711 rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
5712
5713
5714 /**
5715 * If this new raster is empty (width = 0 OR height = 0) then there is
5716 * nothing to compute and we return it right now
5717 **/
5718 if (rt_raster_is_empty(newrast))
5719 {
5720 elog(NOTICE, "Raster is empty. Returning an empty raster");
5721 rt_raster_destroy(raster);
5722 PG_FREE_IF_COPY(pgraster, 0);
5723
5724 pgrtn = rt_raster_serialize(newrast);
5725 rt_raster_destroy(newrast);
5726 if (NULL == pgrtn) {
5727 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5728 PG_RETURN_NULL();
5729 }
5730
5731 SET_VARSIZE(pgrtn, pgrtn->size);
5732 PG_RETURN_POINTER(pgrtn);
5733 }
5734
5735 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Getting raster band %d...", nband);
5736
5737 /**
5738 * Check if the raster has the required band. Otherwise, return a raster
5739 * without band
5740 **/
5741 if (!rt_raster_has_band(raster, nband - 1)) {
5742 elog(NOTICE, "Raster does not have the required band. Returning a raster "
5743 "without a band");
5744 rt_raster_destroy(raster);
5745 PG_FREE_IF_COPY(pgraster, 0);
5746
5747 pgrtn = rt_raster_serialize(newrast);
5748 rt_raster_destroy(newrast);
5749 if (NULL == pgrtn) {
5750 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5751 PG_RETURN_NULL();
5752 }
5753
5754 SET_VARSIZE(pgrtn, pgrtn->size);
5755 PG_RETURN_POINTER(pgrtn);
5756 }
5757
5758 /* Get the raster band */
5759 band = rt_raster_get_band(raster, nband - 1);
5760 if ( NULL == band ) {
5761 elog(NOTICE, "Could not get the required band. Returning a raster "
5762 "without a band");
5763 rt_raster_destroy(raster);
5764 PG_FREE_IF_COPY(pgraster, 0);
5765
5766 pgrtn = rt_raster_serialize(newrast);
5767 rt_raster_destroy(newrast);
5768 if (NULL == pgrtn) {
5769 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5770 PG_RETURN_NULL();
5771 }
5772
5773 SET_VARSIZE(pgrtn, pgrtn->size);
5774 PG_RETURN_POINTER(pgrtn);
5775 }
5776
5777 /*
5778 * Get NODATA value
5779 */
5780 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Getting NODATA value for band...");
5781
5782 if (rt_band_get_hasnodata_flag(band)) {
5783 rt_band_get_nodata(band, &newnodatavalue);
5784 }
5785
5786 else {
5787 newnodatavalue = rt_band_get_min_value(band);
5788 }
5789
5790 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: NODATA value for band: %f",
5791 newnodatavalue);
5792 /**
5793 * We set the initial value of the future band to nodata value. If nodata
5794 * value is null, then the raster will be initialized to
5795 * rt_band_get_min_value but all the values should be recomputed anyway
5796 **/
5797 newinitialvalue = newnodatavalue;
5798
5799 /**
5800 * Set the new pixeltype
5801 **/
5802 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Setting pixeltype...");
5803
5804 if (PG_ARGISNULL(2)) {
5805 newpixeltype = rt_band_get_pixtype(band);
5806 }
5807
5808 else {
5809 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5810 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Pixeltype parameter: %s", strFromText);
5811 newpixeltype = rt_pixtype_index_from_name(strFromText);
5812 pfree(strFromText);
5813 if (newpixeltype == PT_END)
5814 newpixeltype = rt_band_get_pixtype(band);
5815 }
5816
5817 if (newpixeltype == PT_END) {
5818
5819 rt_raster_destroy(raster);
5820 PG_FREE_IF_COPY(pgraster, 0);
5821 rt_raster_destroy(newrast);
5822
5823 elog(ERROR, "RASTER_mapAlgebraFctNgb: Invalid pixeltype");
5824 PG_RETURN_NULL();
5825 }
5826
5827 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Pixeltype set to %s (%d)",
5828 rt_pixtype_name(newpixeltype), newpixeltype);
5829
5830 /* Get the name of the callback userfunction */
5831 if (PG_ARGISNULL(5)) {
5832
5833 rt_raster_destroy(raster);
5834 PG_FREE_IF_COPY(pgraster, 0);
5835 rt_raster_destroy(newrast);
5836
5837 elog(ERROR, "RASTER_mapAlgebraFctNgb: Required function is missing");
5838 PG_RETURN_NULL();
5839 }
5840
5841 oid = PG_GETARG_OID(5);
5842 if (oid == InvalidOid) {
5843
5844 rt_raster_destroy(raster);
5845 PG_FREE_IF_COPY(pgraster, 0);
5846 rt_raster_destroy(newrast);
5847
5848 elog(ERROR, "RASTER_mapAlgebraFctNgb: Got invalid function object id");
5849 PG_RETURN_NULL();
5850 }
5851
5852 fmgr_info(oid, &cbinfo);
5853
5854 /* function cannot return set */
5855 if (cbinfo.fn_retset) {
5856
5857 rt_raster_destroy(raster);
5858 PG_FREE_IF_COPY(pgraster, 0);
5859 rt_raster_destroy(newrast);
5860
5861 elog(ERROR, "RASTER_mapAlgebraFctNgb: Function provided must return double precision not resultset");
5862 PG_RETURN_NULL();
5863 }
5864 /* function should have correct # of args */
5865 else if (cbinfo.fn_nargs != 3) {
5866
5867 rt_raster_destroy(raster);
5868 PG_FREE_IF_COPY(pgraster, 0);
5869 rt_raster_destroy(newrast);
5870
5871 elog(ERROR, "RASTER_mapAlgebraFctNgb: Function does not have three input parameters");
5872 PG_RETURN_NULL();
5873 }
5874
5875 if (func_volatile(oid) == 'v') {
5876 elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5877 }
5878
5879 /* prep function call data */
5880 #if POSTGIS_PGSQL_VERSION < 120
5881 InitFunctionCallInfoData(cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5882 memset(cbdata.argnull, FALSE, sizeof(bool) * 3);
5883 #else
5884 InitFunctionCallInfoData(*cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5885 cbdata->args[0].isnull = FALSE;
5886 cbdata->args[1].isnull = FALSE;
5887 cbdata->args[2].isnull = FALSE;
5888 #endif
5889
5890 /* check that the function isn't strict if the args are null. */
5891 if (PG_ARGISNULL(7)) {
5892 if (cbinfo.fn_strict) {
5893
5894 rt_raster_destroy(raster);
5895 PG_FREE_IF_COPY(pgraster, 0);
5896 rt_raster_destroy(newrast);
5897
5898 elog(ERROR, "RASTER_mapAlgebraFctNgb: Strict callback functions cannot have NULL parameters");
5899 PG_RETURN_NULL();
5900 }
5901
5902 #if POSTGIS_PGSQL_VERSION < 120
5903 cbdata.arg[2] = (Datum)NULL;
5904 cbdata.argnull[2] = TRUE;
5905 #else
5906 cbdata->args[2].value = (Datum)NULL;
5907 cbdata->args[2].isnull = TRUE;
5908 #endif
5909 }
5910 else {
5911 #if POSTGIS_PGSQL_VERSION < 120
5912 cbdata.arg[2] = PG_GETARG_DATUM(7);
5913 #else
5914 cbdata->args[2].value = PG_GETARG_DATUM(7);
5915 #endif
5916 }
5917
5918 /**
5919 * Optimization: If the raster is only filled with nodata values return
5920 * right now a raster filled with the nodatavalueexpr
5921 * TODO: Call rt_band_check_isnodata instead?
5922 **/
5923 if (rt_band_get_isnodata_flag(band)) {
5924
5925 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Band is a nodata band, returning "
5926 "a raster filled with nodata");
5927
5928 rt_raster_generate_new_band(newrast, newpixeltype,
5929 newinitialvalue, TRUE, newnodatavalue, 0);
5930
5931 rt_raster_destroy(raster);
5932 PG_FREE_IF_COPY(pgraster, 0);
5933
5934 /* Serialize created raster */
5935 pgrtn = rt_raster_serialize(newrast);
5936 rt_raster_destroy(newrast);
5937 if (NULL == pgrtn) {
5938 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5939 PG_RETURN_NULL();
5940 }
5941
5942 SET_VARSIZE(pgrtn, pgrtn->size);
5943 PG_RETURN_POINTER(pgrtn);
5944 }
5945
5946
5947 /**
5948 * Create the raster receiving all the computed values. Initialize it to the
5949 * new initial value
5950 **/
5951 rt_raster_generate_new_band(newrast, newpixeltype,
5952 newinitialvalue, TRUE, newnodatavalue, 0);
5953
5954 /* Get the new raster band */
5955 newband = rt_raster_get_band(newrast, 0);
5956 if ( NULL == newband ) {
5957 elog(NOTICE, "Could not modify band for new raster. Returning new "
5958 "raster with the original band");
5959
5960 rt_raster_destroy(raster);
5961 PG_FREE_IF_COPY(pgraster, 0);
5962
5963 /* Serialize created raster */
5964 pgrtn = rt_raster_serialize(newrast);
5965 rt_raster_destroy(newrast);
5966 if (NULL == pgrtn) {
5967 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5968 PG_RETURN_NULL();
5969 }
5970
5971 SET_VARSIZE(pgrtn, pgrtn->size);
5972 PG_RETURN_POINTER(pgrtn);
5973 }
5974
5975 /* Get the width of the neighborhood */
5976 if (PG_ARGISNULL(3) || PG_GETARG_INT32(3) <= 0) {
5977 elog(NOTICE, "Neighborhood width is NULL or <= 0. Returning new "
5978 "raster with the original band");
5979
5980 rt_raster_destroy(raster);
5981 PG_FREE_IF_COPY(pgraster, 0);
5982
5983 /* Serialize created raster */
5984 pgrtn = rt_raster_serialize(newrast);
5985 rt_raster_destroy(newrast);
5986 if (NULL == pgrtn) {
5987 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5988 PG_RETURN_NULL();
5989 }
5990
5991 SET_VARSIZE(pgrtn, pgrtn->size);
5992 PG_RETURN_POINTER(pgrtn);
5993 }
5994
5995 ngbwidth = PG_GETARG_INT32(3);
5996 winwidth = ngbwidth * 2 + 1;
5997
5998 /* Get the height of the neighborhood */
5999 if (PG_ARGISNULL(4) || PG_GETARG_INT32(4) <= 0) {
6000 elog(NOTICE, "Neighborhood height is NULL or <= 0. Returning new "
6001 "raster with the original band");
6002
6003 rt_raster_destroy(raster);
6004 PG_FREE_IF_COPY(pgraster, 0);
6005
6006 /* Serialize created raster */
6007 pgrtn = rt_raster_serialize(newrast);
6008 rt_raster_destroy(newrast);
6009 if (NULL == pgrtn) {
6010 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
6011 PG_RETURN_NULL();
6012 }
6013
6014 SET_VARSIZE(pgrtn, pgrtn->size);
6015 PG_RETURN_POINTER(pgrtn);
6016 }
6017
6018 ngbheight = PG_GETARG_INT32(4);
6019 winheight = ngbheight * 2 + 1;
6020
6021 /* Get the type of NODATA behavior for the neighborhoods. */
6022 if (PG_ARGISNULL(6)) {
6023 elog(NOTICE, "Neighborhood NODATA behavior defaulting to 'ignore'");
6024 txtNodataMode = cstring_to_text("ignore");
6025 }
6026 else {
6027 txtNodataMode = PG_GETARG_TEXT_P(6);
6028 }
6029
6030 txtCallbackParam = (text*)palloc(VARSIZE(txtNodataMode));
6031 SET_VARSIZE(txtCallbackParam, VARSIZE(txtNodataMode));
6032 memcpy((void *)VARDATA(txtCallbackParam), (void *)VARDATA(txtNodataMode), VARSIZE(txtNodataMode) - VARHDRSZ);
6033
6034 /* pass the nodata mode into the user function */
6035 #if POSTGIS_PGSQL_VERSION < 120
6036 cbdata.arg[1] = CStringGetDatum(txtCallbackParam);
6037 #else
6038 cbdata->args[1].value = CStringGetDatum(txtCallbackParam);
6039 #endif
6040
6041 strFromText = text_to_cstring(txtNodataMode);
6042 strFromText = rtpg_strtoupper(strFromText);
6043
6044 if (strcmp(strFromText, "VALUE") == 0)
6045 valuereplace = true;
6046 else if (strcmp(strFromText, "IGNORE") != 0 && strcmp(strFromText, "NULL") != 0) {
6047 /* if the text is not "IGNORE" or "NULL", it may be a numerical value */
6048 if (sscanf(strFromText, "%d", &intReplace) <= 0 && sscanf(strFromText, "%f", &fltReplace) <= 0) {
6049 /* the value is NOT an integer NOR a floating point */
6050 elog(NOTICE, "Neighborhood NODATA mode is not recognized. Must be one of 'value', 'ignore', "
6051 "'NULL', or a numeric value. Returning new raster with the original band");
6052
6053 /* clean up the nodatamode string */
6054 pfree(txtCallbackParam);
6055 pfree(strFromText);
6056
6057 rt_raster_destroy(raster);
6058 PG_FREE_IF_COPY(pgraster, 0);
6059
6060 /* Serialize created raster */
6061 pgrtn = rt_raster_serialize(newrast);
6062 rt_raster_destroy(newrast);
6063 if (NULL == pgrtn) {
6064 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
6065 PG_RETURN_NULL();
6066 }
6067
6068 SET_VARSIZE(pgrtn, pgrtn->size);
6069 PG_RETURN_POINTER(pgrtn);
6070 }
6071 }
6072 else if (strcmp(strFromText, "NULL") == 0) {
6073 /* this setting means that the neighborhood should be skipped if any of the values are null */
6074 nNullSkip = true;
6075 }
6076
6077 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Main computing loop (%d x %d)",
6078 width, height);
6079
6080 /* Allocate room for the neighborhood. */
6081 neighborData = (Datum *)palloc(winwidth * winheight * sizeof(Datum));
6082 neighborNulls = (bool *)palloc(winwidth * winheight * sizeof(bool));
6083
6084 /* The dimensions of the neighborhood array, for creating a multi-dimensional array. */
6085 neighborDims[0] = winwidth;
6086 neighborDims[1] = winheight;
6087
6088 /* The lower bounds for the new multi-dimensional array. */
6089 neighborLbs[0] = 1;
6090 neighborLbs[1] = 1;
6091
6092 /* Get information about the type of item in the multi-dimensional array (float8). */
6093 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
6094
6095 for (x = 0 + ngbwidth; x < width - ngbwidth; x++) {
6096 for(y = 0 + ngbheight; y < height - ngbheight; y++) {
6097 /* populate an array with the pixel values in the neighborhood */
6098 nIndex = 0;
6099 nNullItems = 0;
6100 nNodataOnly = true;
6101 pixelreplace = false;
6102 if (valuereplace) {
6103 ret = rt_band_get_pixel(band, x, y, &rpix, NULL);
6104 if (ret == ES_NONE && FLT_NEQ(rpix, newnodatavalue)) {
6105 pixelreplace = true;
6106 }
6107 }
6108 for (u = x - ngbwidth; u <= x + ngbwidth; u++) {
6109 for (v = y - ngbheight; v <= y + ngbheight; v++) {
6110 ret = rt_band_get_pixel(band, u, v, &r, NULL);
6111 if (ret == ES_NONE) {
6112 if (FLT_NEQ(r, newnodatavalue)) {
6113 /* If the pixel value for this neighbor cell is not NODATA */
6114 neighborData[nIndex] = Float8GetDatum((double)r);
6115 neighborNulls[nIndex] = false;
6116 nNodataOnly = false;
6117 }
6118 else {
6119 /* If the pixel value for this neighbor cell is NODATA */
6120 if (valuereplace && pixelreplace) {
6121 /* Replace the NODATA value with the currently processing pixel. */
6122 neighborData[nIndex] = Float8GetDatum((double)rpix);
6123 neighborNulls[nIndex] = false;
6124 /* do not increment nNullItems, since the user requested that the */
6125 /* neighborhood replace NODATA values with the central pixel value */
6126 }
6127 else {
6128 neighborData[nIndex] = PointerGetDatum(NULL);
6129 neighborNulls[nIndex] = true;
6130 nNullItems++;
6131 }
6132 }
6133 }
6134 else {
6135 /* Fill this will NULL if we can't read the raster pixel. */
6136 neighborData[nIndex] = PointerGetDatum(NULL);
6137 neighborNulls[nIndex] = true;
6138 nNullItems++;
6139 }
6140 /* Next neighbor position */
6141 nIndex++;
6142 }
6143 }
6144
6145 /**
6146 * We compute a value only for the withdata value neighborhood since the
6147 * nodata value has already been set by the first optimization
6148 **/
6149 if (!(nNodataOnly || /* neighborhood only contains NODATA -- OR -- */
6150 (nNullSkip && nNullItems > 0) || /* neighborhood should skip any NODATA cells, and a NODATA cell was detected -- OR -- */
6151 (valuereplace && nNullItems > 0))) { /* neighborhood should replace NODATA cells with the central pixel value, and a NODATA cell was detected */
6152 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: (%dx%d), %dx%d neighborhood",
6153 x, y, winwidth, winheight);
6154
6155 neighborDatum = construct_md_array((void *)neighborData, neighborNulls, 2, neighborDims, neighborLbs,
6156 FLOAT8OID, typlen, typbyval, typalign);
6157
6158 #if POSTGIS_PGSQL_VERSION < 120
6159 /* Assign the neighbor matrix as the first argument to the user function */
6160 cbdata.arg[0] = PointerGetDatum(neighborDatum);
6161
6162 /* Invoke the user function */
6163 tmpnewval = FunctionCallInvoke(&cbdata);
6164
6165 /* Get the return value of the user function */
6166 if (cbdata.isnull) {
6167 newval = newnodatavalue;
6168 }
6169 #else
6170 /* Assign the neighbor matrix as the first argument to the user function */
6171 cbdata->args[0].value = PointerGetDatum(neighborDatum);
6172
6173 /* Invoke the user function */
6174 tmpnewval = FunctionCallInvoke(cbdata);
6175
6176 /* Get the return value of the user function */
6177 if (cbdata->isnull)
6178 {
6179 newval = newnodatavalue;
6180 }
6181 #endif
6182 else {
6183 newval = DatumGetFloat8(tmpnewval);
6184 }
6185
6186 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: new value = %f",
6187 newval);
6188
6189 rt_band_set_pixel(newband, x, y, newval, NULL);
6190 }
6191
6192 /* reset the number of null items in the neighborhood */
6193 nNullItems = 0;
6194 }
6195 }
6196
6197
6198 /* clean up */
6199 pfree(neighborNulls);
6200 pfree(neighborData);
6201 pfree(strFromText);
6202 pfree(txtCallbackParam);
6203
6204 rt_raster_destroy(raster);
6205 PG_FREE_IF_COPY(pgraster, 0);
6206
6207 /* The newrast band has been modified */
6208
6209 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: raster modified, serializing it.");
6210 /* Serialize created raster */
6211
6212 pgrtn = rt_raster_serialize(newrast);
6213 rt_raster_destroy(newrast);
6214 if (NULL == pgrtn)
6215 PG_RETURN_NULL();
6216
6217 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: raster serialized");
6218 POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFctNgb: returning raster");
6219
6220 SET_VARSIZE(pgrtn, pgrtn->size);
6221 PG_RETURN_POINTER(pgrtn);
6222 }
6223
6224 #define ARGKWCOUNT 8
6225
6226 /**
6227 * Two raster MapAlgebra
6228 */
6229 PG_FUNCTION_INFO_V1(RASTER_mapAlgebra2);
6230 Datum RASTER_mapAlgebra2(PG_FUNCTION_ARGS)
6231 {
6232 const uint32_t set_count = 2;
6233 rt_pgraster *pgrast[2] = { NULL, NULL };
6234 int pgrastpos[2] = {-1, -1};
6235 rt_pgraster *pgrtn;
6236 rt_raster rast[2] = {NULL};
6237 int _isempty[2] = {0};
6238 uint32_t bandindex[2] = {0};
6239 rt_raster _rast[2] = {NULL};
6240 rt_band _band[2] = {NULL};
6241 int _hasnodata[2] = {0};
6242 double _nodataval[2] = {0};
6243 double _offset[4] = {0.};
6244 double _rastoffset[2][4] = {{0.}};
6245 int _haspixel[2] = {0};
6246 double _pixel[2] = {0};
6247 int _pos[2][2] = {{0}};
6248 uint16_t _dim[2][2] = {{0}};
6249
6250 char *pixtypename = NULL;
6251 rt_pixtype pixtype = PT_END;
6252 char *extenttypename = NULL;
6253 rt_extenttype extenttype = ET_INTERSECTION;
6254
6255 rt_raster raster = NULL;
6256 rt_band band = NULL;
6257 uint16_t dim[2] = {0};
6258 int haspixel = 0;
6259 double pixel = 0.;
6260 double nodataval = 0;
6261 double gt[6] = {0.};
6262
6263 Oid calltype = InvalidOid;
6264
6265 const uint32_t spi_count = 3;
6266 uint16_t spi_exprpos[3] = {4, 7, 8};
6267 uint32_t spi_argcount[3] = {0};
6268 char *expr = NULL;
6269 char *sql = NULL;
6270 SPIPlanPtr spi_plan[3] = {NULL};
6271 uint16_t spi_empty = 0;
6272 Oid *argtype = NULL;
6273 uint8_t argpos[3][8] = {{0}};
6274 char *argkw[] = {"[rast1.x]", "[rast1.y]", "[rast1.val]", "[rast1]", "[rast2.x]", "[rast2.y]", "[rast2.val]", "[rast2]"};
6275 Datum values[ARGKWCOUNT];
6276 char nulls[ARGKWCOUNT];
6277 TupleDesc tupdesc;
6278 SPITupleTable *tuptable = NULL;
6279 HeapTuple tuple;
6280 Datum datum;
6281 bool isnull = FALSE;
6282 int hasargval[3] = {0};
6283 double argval[3] = {0.};
6284 int hasnodatanodataval = 0;
6285 double nodatanodataval = 0;
6286 int isnodata = 0;
6287
6288 Oid ufc_noid = InvalidOid;
6289 FmgrInfo ufl_info;
6290 #if POSTGIS_PGSQL_VERSION < 120
6291 FunctionCallInfoData ufc_info;
6292 #else
6293 LOCAL_FCINFO(ufc_info, FUNC_MAX_ARGS); /* Could be optimized */
6294 #endif
6295 int ufc_nullcount = 0;
6296
6297 int idx = 0;
6298 uint32_t i = 0;
6299 uint32_t j = 0;
6300 uint32_t k = 0;
6301 uint32_t x = 0;
6302 uint32_t y = 0;
6303 int _x = 0;
6304 int _y = 0;
6305 int err;
6306 int aligned = 0;
6307 int len = 0;
6308
6309 POSTGIS_RT_DEBUG(3, "Starting RASTER_mapAlgebra2");
6310
6311 for (i = 0, j = 0; i < set_count; i++) {
6312 if (!PG_ARGISNULL(j)) {
6313 pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
6314 pgrastpos[i] = j;
6315 j++;
6316
6317 /* raster */
6318 rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
6319 if (!rast[i]) {
6320 for (k = 0; k <= i; k++) {
6321 if (k < i && rast[k] != NULL)
6322 rt_raster_destroy(rast[k]);
6323 if (pgrastpos[k] != -1)
6324 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6325 }
6326 elog(ERROR, "RASTER_mapAlgebra2: Could not deserialize the %s raster", i < 1 ? "first" : "second");
6327 PG_RETURN_NULL();
6328 }
6329
6330 /* empty */
6331 _isempty[i] = rt_raster_is_empty(rast[i]);
6332
6333 /* band index */
6334 if (!PG_ARGISNULL(j)) {
6335 bandindex[i] = PG_GETARG_INT32(j);
6336 }
6337 j++;
6338 }
6339 else {
6340 _isempty[i] = 1;
6341 j += 2;
6342 }
6343
6344 POSTGIS_RT_DEBUGF(3, "_isempty[%d] = %d", i, _isempty[i]);
6345 }
6346
6347 /* both rasters are NULL */
6348 if (rast[0] == NULL && rast[1] == NULL) {
6349 elog(NOTICE, "The two rasters provided are NULL. Returning NULL");
6350 for (k = 0; k < set_count; k++) {
6351 if (pgrastpos[k] != -1)
6352 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6353 }
6354 PG_RETURN_NULL();
6355 }
6356
6357 /* both rasters are empty */
6358 if (_isempty[0] && _isempty[1]) {
6359 elog(NOTICE, "The two rasters provided are empty. Returning empty raster");
6360
6361 raster = rt_raster_new(0, 0);
6362 if (raster == NULL) {
6363 for (k = 0; k < set_count; k++) {
6364 if (rast[k] != NULL)
6365 rt_raster_destroy(rast[k]);
6366 if (pgrastpos[k] != -1)
6367 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6368 }
6369 elog(ERROR, "RASTER_mapAlgebra2: Could not create empty raster");
6370 PG_RETURN_NULL();
6371 }
6372 rt_raster_set_scale(raster, 0, 0);
6373
6374 pgrtn = rt_raster_serialize(raster);
6375 rt_raster_destroy(raster);
6376 if (!pgrtn)
6377 PG_RETURN_NULL();
6378
6379 SET_VARSIZE(pgrtn, pgrtn->size);
6380 PG_RETURN_POINTER(pgrtn);
6381 }
6382
6383 /* replace the empty or NULL raster with one matching the other */
6384 if (
6385 (rast[0] == NULL || _isempty[0]) ||
6386 (rast[1] == NULL || _isempty[1])
6387 ) {
6388 /* first raster is empty */
6389 if (rast[0] == NULL || _isempty[0]) {
6390 i = 0;
6391 j = 1;
6392 }
6393 /* second raster is empty */
6394 else {
6395 i = 1;
6396 j = 0;
6397 }
6398
6399 _rast[j] = rast[j];
6400
6401 /* raster is empty, destroy it */
6402 if (_rast[i] != NULL)
6403 rt_raster_destroy(_rast[i]);
6404
6405 _dim[i][0] = rt_raster_get_width(_rast[j]);
6406 _dim[i][1] = rt_raster_get_height(_rast[j]);
6407 _dim[j][0] = rt_raster_get_width(_rast[j]);
6408 _dim[j][1] = rt_raster_get_height(_rast[j]);
6409
6410 _rast[i] = rt_raster_new(
6411 _dim[j][0],
6412 _dim[j][1]
6413 );
6414 if (_rast[i] == NULL) {
6415 rt_raster_destroy(_rast[j]);
6416 for (k = 0; k < set_count; k++) {
6417 if (pgrastpos[k] != -1)
6418 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6419 }
6420 elog(ERROR, "RASTER_mapAlgebra2: Could not create NODATA raster");
6421 PG_RETURN_NULL();
6422 }
6423 rt_raster_set_srid(_rast[i], rt_raster_get_srid(_rast[j]));
6424
6425 rt_raster_get_geotransform_matrix(_rast[j], gt);
6426 rt_raster_set_geotransform_matrix(_rast[i], gt);
6427 }
6428 else {
6429 _rast[0] = rast[0];
6430 _dim[0][0] = rt_raster_get_width(_rast[0]);
6431 _dim[0][1] = rt_raster_get_height(_rast[0]);
6432
6433 _rast[1] = rast[1];
6434 _dim[1][0] = rt_raster_get_width(_rast[1]);
6435 _dim[1][1] = rt_raster_get_height(_rast[1]);
6436 }
6437
6438 /* SRID must match */
6439 /*
6440 if (rt_raster_get_srid(_rast[0]) != rt_raster_get_srid(_rast[1])) {
6441 elog(NOTICE, "The two rasters provided have different SRIDs. Returning NULL");
6442 for (k = 0; k < set_count; k++) {
6443 if (_rast[k] != NULL)
6444 rt_raster_destroy(_rast[k]);
6445 if (pgrastpos[k] != -1)
6446 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6447 }
6448 PG_RETURN_NULL();
6449 }
6450 */
6451
6452 /* same alignment */
6453 if (rt_raster_same_alignment(_rast[0], _rast[1], &aligned, NULL) != ES_NONE) {
6454 for (k = 0; k < set_count; k++) {
6455 if (_rast[k] != NULL)
6456 rt_raster_destroy(_rast[k]);
6457 if (pgrastpos[k] != -1)
6458 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6459 }
6460 elog(ERROR, "RASTER_mapAlgebra2: Could not test for alignment on the two rasters");
6461 PG_RETURN_NULL();
6462 }
6463 if (!aligned) {
6464 elog(NOTICE, "The two rasters provided do not have the same alignment. Returning NULL");
6465 for (k = 0; k < set_count; k++) {
6466 if (_rast[k] != NULL)
6467 rt_raster_destroy(_rast[k]);
6468 if (pgrastpos[k] != -1)
6469 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6470 }
6471 PG_RETURN_NULL();
6472 }
6473
6474 /* pixel type */
6475 if (!PG_ARGISNULL(5)) {
6476 pixtypename = text_to_cstring(PG_GETARG_TEXT_P(5));
6477 /* Get the pixel type index */
6478 pixtype = rt_pixtype_index_from_name(pixtypename);
6479 if (pixtype == PT_END ) {
6480 for (k = 0; k < set_count; k++) {
6481 if (_rast[k] != NULL)
6482 rt_raster_destroy(_rast[k]);
6483 if (pgrastpos[k] != -1)
6484 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6485 }
6486 elog(ERROR, "RASTER_mapAlgebra2: Invalid pixel type: %s", pixtypename);
6487 PG_RETURN_NULL();
6488 }
6489 }
6490
6491 /* extent type */
6492 if (!PG_ARGISNULL(6)) {
6493 extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(6))));
6494 extenttype = rt_util_extent_type(extenttypename);
6495 }
6496 POSTGIS_RT_DEBUGF(3, "extenttype: %d %s", extenttype, extenttypename);
6497
6498 /* computed raster from extent type */
6499 err = rt_raster_from_two_rasters(
6500 _rast[0], _rast[1],
6501 extenttype,
6502 &raster, _offset
6503 );
6504 if (err != ES_NONE) {
6505 for (k = 0; k < set_count; k++) {
6506 if (_rast[k] != NULL)
6507 rt_raster_destroy(_rast[k]);
6508 if (pgrastpos[k] != -1)
6509 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6510 }
6511 elog(ERROR, "RASTER_mapAlgebra2: Could not get output raster of correct extent");
6512 PG_RETURN_NULL();
6513 }
6514
6515 /* copy offsets */
6516 _rastoffset[0][0] = _offset[0];
6517 _rastoffset[0][1] = _offset[1];
6518 _rastoffset[1][0] = _offset[2];
6519 _rastoffset[1][1] = _offset[3];
6520
6521 /* get output raster dimensions */
6522 dim[0] = rt_raster_get_width(raster);
6523 dim[1] = rt_raster_get_height(raster);
6524
6525 i = 2;
6526 /* handle special cases for extent */
6527 switch (extenttype) {
6528 case ET_FIRST:
6529 i = 0;
6530 /* fall through */
6531 case ET_SECOND:
6532 if (i > 1)
6533 i = 1;
6534
6535 if (
6536 _isempty[i] && (
6537 (extenttype == ET_FIRST && i == 0) ||
6538 (extenttype == ET_SECOND && i == 1)
6539 )
6540 ) {
6541 elog(NOTICE, "The %s raster is NULL. Returning NULL", (i != 1 ? "FIRST" : "SECOND"));
6542 for (k = 0; k < set_count; k++) {
6543 if (_rast[k] != NULL)
6544 rt_raster_destroy(_rast[k]);
6545 if (pgrastpos[k] != -1)
6546 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6547 }
6548 rt_raster_destroy(raster);
6549 PG_RETURN_NULL();
6550 }
6551
6552 /* specified band not found */
6553 if (!rt_raster_has_band(_rast[i], bandindex[i] - 1)) {
6554 elog(NOTICE, "The %s raster does not have the band at index %d. Returning no band raster of correct extent",
6555 (i != 1 ? "FIRST" : "SECOND"), bandindex[i]
6556 );
6557
6558 for (k = 0; k < set_count; k++) {
6559 if (_rast[k] != NULL)
6560 rt_raster_destroy(_rast[k]);
6561 if (pgrastpos[k] != -1)
6562 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6563 }
6564
6565 pgrtn = rt_raster_serialize(raster);
6566 rt_raster_destroy(raster);
6567 if (!pgrtn) PG_RETURN_NULL();
6568
6569 SET_VARSIZE(pgrtn, pgrtn->size);
6570 PG_RETURN_POINTER(pgrtn);
6571 }
6572 break;
6573 case ET_UNION:
6574 break;
6575 case ET_INTERSECTION:
6576 /* no intersection */
6577 if (
6578 _isempty[0] || _isempty[1] ||
6579 !dim[0] || !dim[1]
6580 ) {
6581 elog(NOTICE, "The two rasters provided have no intersection. Returning no band raster");
6582
6583 /* raster has dimension, replace with no band raster */
6584 if (dim[0] || dim[1]) {
6585 rt_raster_destroy(raster);
6586
6587 raster = rt_raster_new(0, 0);
6588 if (raster == NULL) {
6589 for (k = 0; k < set_count; k++) {
6590 if (_rast[k] != NULL)
6591 rt_raster_destroy(_rast[k]);
6592 if (pgrastpos[k] != -1)
6593 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6594 }
6595 elog(ERROR, "RASTER_mapAlgebra2: Could not create no band raster");
6596 PG_RETURN_NULL();
6597 }
6598
6599 rt_raster_set_scale(raster, 0, 0);
6600 rt_raster_set_srid(raster, rt_raster_get_srid(_rast[0]));
6601 }
6602
6603 for (k = 0; k < set_count; k++) {
6604 if (_rast[k] != NULL)
6605 rt_raster_destroy(_rast[k]);
6606 if (pgrastpos[k] != -1)
6607 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6608 }
6609
6610 pgrtn = rt_raster_serialize(raster);
6611 rt_raster_destroy(raster);
6612 if (!pgrtn) PG_RETURN_NULL();
6613
6614 SET_VARSIZE(pgrtn, pgrtn->size);
6615 PG_RETURN_POINTER(pgrtn);
6616 }
6617 break;
6618 case ET_LAST:
6619 case ET_CUSTOM:
6620 for (k = 0; k < set_count; k++) {
6621 if (_rast[k] != NULL)
6622 rt_raster_destroy(_rast[k]);
6623 if (pgrastpos[k] != -1)
6624 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6625 }
6626 elog(ERROR, "RASTER_mapAlgebra2: ET_LAST and ET_CUSTOM are not implemented");
6627 PG_RETURN_NULL();
6628 break;
6629 }
6630
6631 /* both rasters do not have specified bands */
6632 if (
6633 (!_isempty[0] && !rt_raster_has_band(_rast[0], bandindex[0] - 1)) &&
6634 (!_isempty[1] && !rt_raster_has_band(_rast[1], bandindex[1] - 1))
6635 ) {
6636 elog(NOTICE, "The two rasters provided do not have the respectively specified band indices. Returning no band raster of correct extent");
6637
6638 for (k = 0; k < set_count; k++) {
6639 if (_rast[k] != NULL)
6640 rt_raster_destroy(_rast[k]);
6641 if (pgrastpos[k] != -1)
6642 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6643 }
6644
6645 pgrtn = rt_raster_serialize(raster);
6646 rt_raster_destroy(raster);
6647 if (!pgrtn) PG_RETURN_NULL();
6648
6649 SET_VARSIZE(pgrtn, pgrtn->size);
6650 PG_RETURN_POINTER(pgrtn);
6651 }
6652
6653 /* get bands */
6654 for (i = 0; i < set_count; i++) {
6655 if (_isempty[i] || !rt_raster_has_band(_rast[i], bandindex[i] - 1)) {
6656 _hasnodata[i] = 1;
6657 _nodataval[i] = 0;
6658
6659 continue;
6660 }
6661
6662 _band[i] = rt_raster_get_band(_rast[i], bandindex[i] - 1);
6663 if (_band[i] == NULL) {
6664 for (k = 0; k < set_count; k++) {
6665 if (_rast[k] != NULL)
6666 rt_raster_destroy(_rast[k]);
6667 if (pgrastpos[k] != -1)
6668 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6669 }
6670 rt_raster_destroy(raster);
6671 elog(ERROR, "RASTER_mapAlgebra2: Could not get band %d of the %s raster",
6672 bandindex[i],
6673 (i < 1 ? "FIRST" : "SECOND")
6674 );
6675 PG_RETURN_NULL();
6676 }
6677
6678 _hasnodata[i] = rt_band_get_hasnodata_flag(_band[i]);
6679 if (_hasnodata[i])
6680 rt_band_get_nodata(_band[i], &(_nodataval[i]));
6681 }
6682
6683 /* pixtype is PT_END, get pixtype based upon extent */
6684 if (pixtype == PT_END) {
6685 if ((extenttype == ET_SECOND && !_isempty[1]) || _isempty[0])
6686 pixtype = rt_band_get_pixtype(_band[1]);
6687 else
6688 pixtype = rt_band_get_pixtype(_band[0]);
6689 }
6690
6691 /* nodata value for new band */
6692 if (extenttype == ET_SECOND && !_isempty[1] && _hasnodata[1]) {
6693 nodataval = _nodataval[1];
6694 }
6695 else if (!_isempty[0] && _hasnodata[0]) {
6696 nodataval = _nodataval[0];
6697 }
6698 else if (!_isempty[1] && _hasnodata[1]) {
6699 nodataval = _nodataval[1];
6700 }
6701 else {
6702 elog(NOTICE, "Neither raster provided has a NODATA value for the specified band indices. NODATA value set to minimum possible for %s", rt_pixtype_name(pixtype));
6703 nodataval = rt_pixtype_get_min_value(pixtype);
6704 }
6705
6706 /* add band to output raster */
6707 if (rt_raster_generate_new_band(
6708 raster,
6709 pixtype,
6710 nodataval,
6711 1, nodataval,
6712 0
6713 ) < 0) {
6714 for (k = 0; k < set_count; k++) {
6715 if (_rast[k] != NULL)
6716 rt_raster_destroy(_rast[k]);
6717 if (pgrastpos[k] != -1)
6718 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6719 }
6720 rt_raster_destroy(raster);
6721 elog(ERROR, "RASTER_mapAlgebra2: Could not add new band to output raster");
6722 PG_RETURN_NULL();
6723 }
6724
6725 /* get output band */
6726 band = rt_raster_get_band(raster, 0);
6727 if (band == NULL) {
6728 for (k = 0; k < set_count; k++) {
6729 if (_rast[k] != NULL)
6730 rt_raster_destroy(_rast[k]);
6731 if (pgrastpos[k] != -1)
6732 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6733 }
6734 rt_raster_destroy(raster);
6735 elog(ERROR, "RASTER_mapAlgebra2: Could not get newly added band of output raster");
6736 PG_RETURN_NULL();
6737 }
6738
6739 POSTGIS_RT_DEBUGF(4, "offsets = (%d, %d, %d, %d)",
6740 (int) _rastoffset[0][0],
6741 (int) _rastoffset[0][1],
6742 (int) _rastoffset[1][0],
6743 (int) _rastoffset[1][1]
6744 );
6745
6746 POSTGIS_RT_DEBUGF(4, "metadata = (%f, %f, %d, %d, %f, %f, %f, %f, %d)",
6747 rt_raster_get_x_offset(raster),
6748 rt_raster_get_y_offset(raster),
6749 rt_raster_get_width(raster),
6750 rt_raster_get_height(raster),
6751 rt_raster_get_x_scale(raster),
6752 rt_raster_get_y_scale(raster),
6753 rt_raster_get_x_skew(raster),
6754 rt_raster_get_y_skew(raster),
6755 rt_raster_get_srid(raster)
6756 );
6757
6758 /*
6759 determine who called this function
6760 Arg 4 will either be text or regprocedure
6761 */
6762 POSTGIS_RT_DEBUG(3, "checking parameter type for arg 4");
6763 calltype = get_fn_expr_argtype(fcinfo->flinfo, 4);
6764
6765 switch(calltype) {
6766 case TEXTOID: {
6767 POSTGIS_RT_DEBUG(3, "arg 4 is \"expression\"!");
6768
6769 /* connect SPI */
6770 if (SPI_connect() != SPI_OK_CONNECT) {
6771 for (k = 0; k < set_count; k++) {
6772 if (_rast[k] != NULL)
6773 rt_raster_destroy(_rast[k]);
6774 if (pgrastpos[k] != -1)
6775 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6776 }
6777 rt_raster_destroy(raster);
6778 elog(ERROR, "RASTER_mapAlgebra2: Could not connect to the SPI manager");
6779 PG_RETURN_NULL();
6780 }
6781
6782 /* reset hasargval */
6783 memset(hasargval, 0, sizeof(int) * spi_count);
6784
6785 /*
6786 process expressions
6787
6788 spi_exprpos elements are:
6789 4 - expression => spi_plan[0]
6790 7 - nodata1expr => spi_plan[1]
6791 8 - nodata2expr => spi_plan[2]
6792 */
6793 for (i = 0; i < spi_count; i++) {
6794 if (!PG_ARGISNULL(spi_exprpos[i])) {
6795 char *tmp = NULL;
6796 char place[5] = "$1";
6797 expr = text_to_cstring(PG_GETARG_TEXT_P(spi_exprpos[i]));
6798 POSTGIS_RT_DEBUGF(3, "raw expr #%d: %s", i, expr);
6799
6800 for (j = 0, k = 1; j < ARGKWCOUNT; j++) {
6801 /* attempt to replace keyword with placeholder */
6802 len = 0;
6803 tmp = rtpg_strreplace(expr, argkw[j], place, &len);
6804 pfree(expr);
6805 expr = tmp;
6806
6807 if (len) {
6808 spi_argcount[i]++;
6809 argpos[i][j] = k++;
6810
6811 sprintf(place, "$%d", k);
6812 }
6813 else
6814 argpos[i][j] = 0;
6815 }
6816
6817 len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision");
6818 sql = (char *) palloc(len + 1);
6819 if (sql == NULL) {
6820
6821 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6822 SPI_finish();
6823
6824 for (k = 0; k < set_count; k++) {
6825 if (_rast[k] != NULL)
6826 rt_raster_destroy(_rast[k]);
6827 if (pgrastpos[k] != -1)
6828 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6829 }
6830 rt_raster_destroy(raster);
6831
6832 elog(ERROR, "RASTER_mapAlgebra2: Could not allocate memory for expression parameter %d", spi_exprpos[i]);
6833 PG_RETURN_NULL();
6834 }
6835
6836 memcpy(sql, "SELECT (", strlen("SELECT ("));
6837 memcpy(sql + strlen("SELECT ("), expr, strlen(expr));
6838 memcpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision"));
6839 sql[len] = '\0';
6840
6841 POSTGIS_RT_DEBUGF(3, "sql #%d: %s", i, sql);
6842
6843 /* create prepared plan */
6844 if (spi_argcount[i]) {
6845 argtype = (Oid *) palloc(spi_argcount[i] * sizeof(Oid));
6846 if (argtype == NULL) {
6847
6848 pfree(sql);
6849 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6850 SPI_finish();
6851
6852 for (k = 0; k < set_count; k++) {
6853 if (_rast[k] != NULL)
6854 rt_raster_destroy(_rast[k]);
6855 if (pgrastpos[k] != -1)
6856 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6857 }
6858 rt_raster_destroy(raster);
6859
6860 elog(ERROR, "RASTER_mapAlgebra2: Could not allocate memory for prepared plan argtypes of expression parameter %d", spi_exprpos[i]);
6861 PG_RETURN_NULL();
6862 }
6863
6864 /* specify datatypes of parameters */
6865 for (j = 0, k = 0; j < ARGKWCOUNT; j++) {
6866 if (argpos[i][j] < 1) continue;
6867
6868 /* positions are INT4 */
6869 if (
6870 (strstr(argkw[j], "[rast1.x]") != NULL) ||
6871 (strstr(argkw[j], "[rast1.y]") != NULL) ||
6872 (strstr(argkw[j], "[rast2.x]") != NULL) ||
6873 (strstr(argkw[j], "[rast2.y]") != NULL)
6874 ) {
6875 argtype[k] = INT4OID;
6876 }
6877 /* everything else is FLOAT8 */
6878 else {
6879 argtype[k] = FLOAT8OID;
6880 }
6881
6882 k++;
6883 }
6884
6885 spi_plan[i] = SPI_prepare(sql, spi_argcount[i], argtype);
6886 pfree(argtype);
6887
6888 if (spi_plan[i] == NULL) {
6889
6890 pfree(sql);
6891 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6892 SPI_finish();
6893
6894 for (k = 0; k < set_count; k++) {
6895 if (_rast[k] != NULL)
6896 rt_raster_destroy(_rast[k]);
6897 if (pgrastpos[k] != -1)
6898 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6899 }
6900 rt_raster_destroy(raster);
6901
6902 elog(ERROR, "RASTER_mapAlgebra2: Could not create prepared plan of expression parameter %d", spi_exprpos[i]);
6903 PG_RETURN_NULL();
6904 }
6905 }
6906 /* no args, just execute query */
6907 else {
6908 err = SPI_execute(sql, TRUE, 0);
6909 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
6910
6911 pfree(sql);
6912 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6913 SPI_finish();
6914
6915 for (k = 0; k < set_count; k++) {
6916 if (_rast[k] != NULL)
6917 rt_raster_destroy(_rast[k]);
6918 if (pgrastpos[k] != -1)
6919 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6920 }
6921 rt_raster_destroy(raster);
6922
6923 elog(ERROR, "RASTER_mapAlgebra2: Could not evaluate expression parameter %d", spi_exprpos[i]);
6924 PG_RETURN_NULL();
6925 }
6926
6927 /* get output of prepared plan */
6928 tupdesc = SPI_tuptable->tupdesc;
6929 tuptable = SPI_tuptable;
6930 tuple = tuptable->vals[0];
6931
6932 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
6933 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
6934
6935 pfree(sql);
6936 if (SPI_tuptable) SPI_freetuptable(tuptable);
6937 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6938 SPI_finish();
6939
6940 for (k = 0; k < set_count; k++) {
6941 if (_rast[k] != NULL)
6942 rt_raster_destroy(_rast[k]);
6943 if (pgrastpos[k] != -1)
6944 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6945 }
6946 rt_raster_destroy(raster);
6947
6948 elog(ERROR, "RASTER_mapAlgebra2: Could not get result of expression parameter %d", spi_exprpos[i]);
6949 PG_RETURN_NULL();
6950 }
6951
6952 if (!isnull) {
6953 hasargval[i] = 1;
6954 argval[i] = DatumGetFloat8(datum);
6955 }
6956
6957 if (SPI_tuptable) SPI_freetuptable(tuptable);
6958 }
6959
6960 pfree(sql);
6961 }
6962 else
6963 spi_empty++;
6964 }
6965
6966 /* nodatanodataval */
6967 if (!PG_ARGISNULL(9)) {
6968 hasnodatanodataval = 1;
6969 nodatanodataval = PG_GETARG_FLOAT8(9);
6970 }
6971 else
6972 hasnodatanodataval = 0;
6973 break;
6974 }
6975 case REGPROCEDUREOID: {
6976 POSTGIS_RT_DEBUG(3, "arg 4 is \"userfunction\"!");
6977 if (!PG_ARGISNULL(4)) {
6978
6979 ufc_nullcount = 0;
6980 ufc_noid = PG_GETARG_OID(4);
6981
6982 /* get function info */
6983 fmgr_info(ufc_noid, &ufl_info);
6984
6985 /* function cannot return set */
6986 err = 0;
6987 if (ufl_info.fn_retset) {
6988 err = 1;
6989 }
6990 /* function should have correct # of args */
6991 else if (ufl_info.fn_nargs < 3 || ufl_info.fn_nargs > 4) {
6992 err = 2;
6993 }
6994
6995 /*
6996 TODO: consider adding checks of the userfunction parameters
6997 should be able to use get_fn_expr_argtype() of fmgr.c
6998 */
6999
7000 if (err > 0) {
7001 for (k = 0; k < set_count; k++) {
7002 if (_rast[k] != NULL)
7003 rt_raster_destroy(_rast[k]);
7004 if (pgrastpos[k] != -1)
7005 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7006 }
7007 rt_raster_destroy(raster);
7008
7009 if (err > 1)
7010 elog(ERROR, "RASTER_mapAlgebra2: Function provided must have three or four input parameters");
7011 else
7012 elog(ERROR, "RASTER_mapAlgebra2: Function provided must return double precision not resultset");
7013 PG_RETURN_NULL();
7014 }
7015
7016 if (func_volatile(ufc_noid) == 'v') {
7017 elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
7018 }
7019
7020 /* prep function call data */
7021 #if POSTGIS_PGSQL_VERSION < 120
7022 InitFunctionCallInfoData(ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7023 memset(ufc_info.argnull, FALSE, sizeof(bool) * ufl_info.fn_nargs);
7024 #else
7025 InitFunctionCallInfoData(
7026 *ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7027 ufc_info->args[0].isnull = FALSE;
7028 ufc_info->args[1].isnull = FALSE;
7029 ufc_info->args[2].isnull = FALSE;
7030 if (ufl_info.fn_nargs == 4)
7031 ufc_info->args[3].isnull = FALSE;
7032 #endif
7033
7034 if (ufl_info.fn_nargs != 4)
7035 k = 2;
7036 else
7037 k = 3;
7038 #if POSTGIS_PGSQL_VERSION < 120
7039 if (!PG_ARGISNULL(7)) {
7040 ufc_info.arg[k] = PG_GETARG_DATUM(7);
7041 }
7042 else {
7043 ufc_info.arg[k] = (Datum) NULL;
7044 ufc_info.argnull[k] = TRUE;
7045 ufc_nullcount++;
7046 }
7047 #else
7048 if (!PG_ARGISNULL(7))
7049 {
7050 ufc_info->args[k].value = PG_GETARG_DATUM(7);
7051 }
7052 else
7053 {
7054 ufc_info->args[k].value = (Datum)NULL;
7055 ufc_info->args[k].isnull = TRUE;
7056 ufc_nullcount++;
7057 }
7058 #endif
7059 }
7060 break;
7061 }
7062 default:
7063 for (k = 0; k < set_count; k++) {
7064 if (_rast[k] != NULL)
7065 rt_raster_destroy(_rast[k]);
7066 if (pgrastpos[k] != -1)
7067 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7068 }
7069 rt_raster_destroy(raster);
7070 elog(ERROR, "RASTER_mapAlgebra2: Invalid data type for expression or userfunction");
7071 PG_RETURN_NULL();
7072 break;
7073 }
7074
7075 /* loop over pixels */
7076 /* if any expression present, run */
7077 if ((
7078 (calltype == TEXTOID) && (
7079 (spi_empty != spi_count) || hasnodatanodataval
7080 )
7081 ) || (
7082 (calltype == REGPROCEDUREOID) && (ufc_noid != InvalidOid)
7083 )) {
7084 for (x = 0; x < dim[0]; x++) {
7085 for (y = 0; y < dim[1]; y++) {
7086
7087 /* get pixel from each raster */
7088 for (i = 0; i < set_count; i++) {
7089 _haspixel[i] = 0;
7090 _pixel[i] = 0;
7091
7092 /* row/column */
7093 _x = (int)x - (int)_rastoffset[i][0];
7094 _y = (int)y - (int)_rastoffset[i][1];
7095
7096 /* store _x and _y in 1-based */
7097 _pos[i][0] = _x + 1;
7098 _pos[i][1] = _y + 1;
7099
7100 /* get pixel value */
7101 if (_band[i] == NULL) {
7102 if (!_hasnodata[i]) {
7103 _haspixel[i] = 1;
7104 _pixel[i] = _nodataval[i];
7105 }
7106 }
7107 else if (
7108 !_isempty[i] &&
7109 (_x >= 0 && _x < _dim[i][0]) &&
7110 (_y >= 0 && _y < _dim[i][1])
7111 ) {
7112 err = rt_band_get_pixel(_band[i], _x, _y, &(_pixel[i]), &isnodata);
7113 if (err != ES_NONE) {
7114
7115 if (calltype == TEXTOID) {
7116 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7117 SPI_finish();
7118 }
7119
7120 for (k = 0; k < set_count; k++) {
7121 if (_rast[k] != NULL)
7122 rt_raster_destroy(_rast[k]);
7123 if (pgrastpos[k] != -1)
7124 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7125 }
7126 rt_raster_destroy(raster);
7127
7128 elog(ERROR, "RASTER_mapAlgebra2: Could not get pixel of %s raster", (i < 1 ? "FIRST" : "SECOND"));
7129 PG_RETURN_NULL();
7130 }
7131
7132 if (!_hasnodata[i] || !isnodata)
7133 _haspixel[i] = 1;
7134 }
7135
7136 POSTGIS_RT_DEBUGF(5, "pixel r%d(%d, %d) = %d, %f",
7137 i,
7138 _x, _y,
7139 _haspixel[i],
7140 _pixel[i]
7141 );
7142 }
7143
7144 haspixel = 0;
7145
7146 switch (calltype) {
7147 case TEXTOID: {
7148 /* which prepared plan to use? */
7149 /* !pixel0 && !pixel1 */
7150 /* use nodatanodataval */
7151 if (!_haspixel[0] && !_haspixel[1])
7152 i = 3;
7153 /* pixel0 && !pixel1 */
7154 /* run spi_plan[2] (nodata2expr) */
7155 else if (_haspixel[0] && !_haspixel[1])
7156 i = 2;
7157 /* !pixel0 && pixel1 */
7158 /* run spi_plan[1] (nodata1expr) */
7159 else if (!_haspixel[0] && _haspixel[1])
7160 i = 1;
7161 /* pixel0 && pixel1 */
7162 /* run spi_plan[0] (expression) */
7163 else
7164 i = 0;
7165
7166 /* process values */
7167 if (i == 3) {
7168 if (hasnodatanodataval) {
7169 haspixel = 1;
7170 pixel = nodatanodataval;
7171 }
7172 }
7173 /* has an evaluated value */
7174 else if (hasargval[i]) {
7175 haspixel = 1;
7176 pixel = argval[i];
7177 }
7178 /* prepared plan exists */
7179 else if (spi_plan[i] != NULL) {
7180 POSTGIS_RT_DEBUGF(4, "Using prepared plan: %d", i);
7181
7182 /* reset values to (Datum) NULL */
7183 memset(values, (Datum) NULL, sizeof(Datum) * ARGKWCOUNT);
7184 /* reset nulls to FALSE */
7185 memset(nulls, FALSE, sizeof(char) * ARGKWCOUNT);
7186
7187 /* expression has argument(s) */
7188 if (spi_argcount[i]) {
7189 /* set values and nulls */
7190 for (j = 0; j < ARGKWCOUNT; j++) {
7191 idx = argpos[i][j];
7192 if (idx < 1) continue;
7193 idx--; /* 1-based becomes 0-based */
7194
7195 if (strstr(argkw[j], "[rast1.x]") != NULL) {
7196 values[idx] = _pos[0][0];
7197 }
7198 else if (strstr(argkw[j], "[rast1.y]") != NULL) {
7199 values[idx] = _pos[0][1];
7200 }
7201 else if (
7202 (strstr(argkw[j], "[rast1.val]") != NULL) ||
7203 (strstr(argkw[j], "[rast1]") != NULL)
7204 ) {
7205 if (_isempty[0] || !_haspixel[0])
7206 nulls[idx] = TRUE;
7207 else
7208 values[idx] = Float8GetDatum(_pixel[0]);
7209 }
7210 else if (strstr(argkw[j], "[rast2.x]") != NULL) {
7211 values[idx] = _pos[1][0];
7212 }
7213 else if (strstr(argkw[j], "[rast2.y]") != NULL) {
7214 values[idx] = _pos[1][1];
7215 }
7216 else if (
7217 (strstr(argkw[j], "[rast2.val]") != NULL) ||
7218 (strstr(argkw[j], "[rast2]") != NULL)
7219 ) {
7220 if (_isempty[1] || !_haspixel[1])
7221 nulls[idx] = TRUE;
7222 else
7223 values[idx] = Float8GetDatum(_pixel[1]);
7224 }
7225 }
7226 }
7227
7228 /* run prepared plan */
7229 err = SPI_execute_plan(spi_plan[i], values, nulls, TRUE, 1);
7230 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
7231
7232 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7233 SPI_finish();
7234
7235 for (k = 0; k < set_count; k++) {
7236 if (_rast[k] != NULL)
7237 rt_raster_destroy(_rast[k]);
7238 if (pgrastpos[k] != -1)
7239 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7240 }
7241 rt_raster_destroy(raster);
7242
7243 elog(ERROR, "RASTER_mapAlgebra2: Unexpected error when running prepared statement %d", i);
7244 PG_RETURN_NULL();
7245 }
7246
7247 /* get output of prepared plan */
7248 tupdesc = SPI_tuptable->tupdesc;
7249 tuptable = SPI_tuptable;
7250 tuple = tuptable->vals[0];
7251
7252 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
7253 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
7254
7255 if (SPI_tuptable) SPI_freetuptable(tuptable);
7256 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7257 SPI_finish();
7258
7259 for (k = 0; k < set_count; k++) {
7260 if (_rast[k] != NULL)
7261 rt_raster_destroy(_rast[k]);
7262 if (pgrastpos[k] != -1)
7263 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7264 }
7265 rt_raster_destroy(raster);
7266
7267 elog(ERROR, "RASTER_mapAlgebra2: Could not get result of prepared statement %d", i);
7268 PG_RETURN_NULL();
7269 }
7270
7271 if (!isnull) {
7272 haspixel = 1;
7273 pixel = DatumGetFloat8(datum);
7274 }
7275
7276 if (SPI_tuptable) SPI_freetuptable(tuptable);
7277 }
7278 } break;
7279 case REGPROCEDUREOID: {
7280 Datum d[4];
7281 ArrayType *a;
7282
7283 /* build fcnarg */
7284 for (i = 0; i < set_count; i++) {
7285 #if POSTGIS_PGSQL_VERSION < 120
7286 ufc_info.arg[i] = Float8GetDatum(_pixel[i]);
7287 #else
7288 ufc_info->args[i].value = Float8GetDatum(_pixel[i]);
7289 #endif
7290
7291 if (_haspixel[i]) {
7292 #if POSTGIS_PGSQL_VERSION < 120
7293 ufc_info.argnull[i] = FALSE;
7294 #else
7295 ufc_info->args[i].isnull = FALSE;
7296 #endif
7297 ufc_nullcount--;
7298 }
7299 else {
7300 #if POSTGIS_PGSQL_VERSION < 120
7301 ufc_info.argnull[i] = TRUE;
7302 #else
7303 ufc_info->args[i].isnull = TRUE;
7304 #endif
7305 ufc_nullcount++;
7306 }
7307 }
7308
7309 /* function is strict and null parameter is passed */
7310 /* http://archives.postgresql.org/pgsql-general/2011-11/msg00424.php */
7311 if (ufl_info.fn_strict && ufc_nullcount)
7312 break;
7313
7314 /* 4 parameters, add position */
7315 if (ufl_info.fn_nargs == 4) {
7316 /* Datum of 4 element array */
7317 /* array is (x1, y1, x2, y2) */
7318 for (i = 0; i < set_count; i++) {
7319 if (i < 1) {
7320 d[0] = Int32GetDatum(_pos[i][0]);
7321 d[1] = Int32GetDatum(_pos[i][1]);
7322 }
7323 else {
7324 d[2] = Int32GetDatum(_pos[i][0]);
7325 d[3] = Int32GetDatum(_pos[i][1]);
7326 }
7327 }
7328
7329 a = construct_array(d, 4, INT4OID, sizeof(int32), true, 'i');
7330 #if POSTGIS_PGSQL_VERSION < 120
7331 ufc_info.arg[2] = PointerGetDatum(a);
7332 ufc_info.argnull[2] = FALSE;
7333 #else
7334 ufc_info->args[2].value = PointerGetDatum(a);
7335 ufc_info->args[2].isnull = FALSE;
7336 #endif
7337 }
7338
7339 #if POSTGIS_PGSQL_VERSION < 120
7340 datum = FunctionCallInvoke(&ufc_info);
7341
7342 /* result is not null*/
7343 if (!ufc_info.isnull) {
7344 haspixel = 1;
7345 pixel = DatumGetFloat8(datum);
7346 }
7347 #else
7348 datum = FunctionCallInvoke(ufc_info);
7349
7350 /* result is not null*/
7351 if (!ufc_info->isnull)
7352 {
7353 haspixel = 1;
7354 pixel = DatumGetFloat8(datum);
7355 }
7356 #endif
7357 } break;
7358 }
7359
7360 /* burn pixel if haspixel != 0 */
7361 if (haspixel) {
7362 if (rt_band_set_pixel(band, x, y, pixel, NULL) != ES_NONE) {
7363
7364 if (calltype == TEXTOID) {
7365 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7366 SPI_finish();
7367 }
7368
7369 for (k = 0; k < set_count; k++) {
7370 if (_rast[k] != NULL)
7371 rt_raster_destroy(_rast[k]);
7372 if (pgrastpos[k] != -1)
7373 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7374 }
7375 rt_raster_destroy(raster);
7376
7377 elog(ERROR, "RASTER_mapAlgebra2: Could not set pixel value of output raster");
7378 PG_RETURN_NULL();
7379 }
7380 }
7381
7382 POSTGIS_RT_DEBUGF(5, "(x, y, val) = (%d, %d, %f)", x, y, haspixel ? pixel : nodataval);
7383
7384 } /* y: height */
7385 } /* x: width */
7386 }
7387
7388 /* CLEANUP */
7389 if (calltype == TEXTOID) {
7390 for (i = 0; i < spi_count; i++) {
7391 if (spi_plan[i] != NULL) SPI_freeplan(spi_plan[i]);
7392 }
7393 SPI_finish();
7394 }
7395
7396 for (k = 0; k < set_count; k++) {
7397 if (_rast[k] != NULL)
7398 rt_raster_destroy(_rast[k]);
7399 if (pgrastpos[k] != -1)
7400 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7401 }
7402
7403 pgrtn = rt_raster_serialize(raster);
7404 rt_raster_destroy(raster);
7405 if (!pgrtn) PG_RETURN_NULL();
7406
7407 POSTGIS_RT_DEBUG(3, "Finished RASTER_mapAlgebra2");
7408
7409 SET_VARSIZE(pgrtn, pgrtn->size);
7410 PG_RETURN_POINTER(pgrtn);
7411 }
7412