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 unsigned char *wkb = NULL;
3038 size_t wkb_len;
3039
3040 ArrayType *array;
3041 Oid etype;
3042 Datum *e;
3043 bool *nulls;
3044
3045 int16 typlen;
3046 bool typbyval;
3047 char typalign;
3048
3049 int i = 0;
3050 int j = 0;
3051 int k = 0;
3052 rtpg_clip_arg arg = NULL;
3053 LWGEOM *tmpgeom = NULL;
3054 rt_iterator itrset;
3055
3056 rt_raster _raster = NULL;
3057 rt_band band = NULL;
3058 rt_pixtype pixtype;
3059 int hasnodata;
3060 double nodataval;
3061 int noerr = 0;
3062
3063 POSTGIS_RT_DEBUG(3, "Starting...");
3064
3065 /* raster or geometry is NULL, return NULL */
3066 if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3067 PG_RETURN_NULL();
3068
3069 /* init arg */
3070 arg = rtpg_clip_arg_init();
3071 if (arg == NULL) {
3072 elog(ERROR, "RASTER_clip: Could not initialize argument structure");
3073 PG_RETURN_NULL();
3074 }
3075
3076 /* raster (0) */
3077 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3078
3079 /* Get raster object */
3080 arg->raster = rt_raster_deserialize(pgraster, FALSE);
3081 if (arg->raster == NULL) {
3082 rtpg_clip_arg_destroy(arg);
3083 PG_FREE_IF_COPY(pgraster, 0);
3084 elog(ERROR, "RASTER_clip: Could not deserialize raster");
3085 PG_RETURN_NULL();
3086 }
3087
3088 /* raster is empty, return empty raster */
3089 if (rt_raster_is_empty(arg->raster) || rt_raster_get_num_bands(arg->raster) == 0) {
3090 elog(NOTICE, "Input raster is empty or has no bands. Returning empty raster");
3091
3092 rtpg_clip_arg_destroy(arg);
3093 PG_FREE_IF_COPY(pgraster, 0);
3094
3095 rtn = rt_raster_new(0, 0);
3096 if (rtn == NULL) {
3097 elog(ERROR, "RASTER_clip: Could not create empty raster");
3098 PG_RETURN_NULL();
3099 }
3100
3101 pgrtn = rt_raster_serialize(rtn);
3102 rt_raster_destroy(rtn);
3103 if (NULL == pgrtn)
3104 PG_RETURN_NULL();
3105
3106 SET_VARSIZE(pgrtn, pgrtn->size);
3107 PG_RETURN_POINTER(pgrtn);
3108 }
3109
3110 /* metadata */
3111 rt_raster_get_geotransform_matrix(arg->raster, gt);
3112 srid = clamp_srid(rt_raster_get_srid(arg->raster));
3113
3114 /* geometry (2) */
3115 gser = PG_GETARG_GSERIALIZED_P(2);
3116 geom = lwgeom_from_gserialized(gser);
3117
3118 /* Get a 2D version of the geometry if necessary */
3119 if (lwgeom_ndims(geom) > 2) {
3120 LWGEOM *geom2d = lwgeom_force_2d(geom);
3121 lwgeom_free(geom);
3122 geom = geom2d;
3123 }
3124
3125 /* check that SRIDs match */
3126 if (srid != clamp_srid(gserialized_get_srid(gser))) {
3127 elog(NOTICE, "Geometry provided does not have the same SRID as the raster. Returning NULL");
3128
3129 rtpg_clip_arg_destroy(arg);
3130 PG_FREE_IF_COPY(pgraster, 0);
3131 lwgeom_free(geom);
3132 PG_FREE_IF_COPY(gser, 2);
3133
3134 PG_RETURN_NULL();
3135 }
3136
3137 /* crop (4) */
3138 if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3139 arg->extenttype = ET_FIRST;
3140
3141 /* get intersection geometry of input raster and input geometry */
3142 if (rt_raster_get_convex_hull(arg->raster, &rastgeom) != ES_NONE) {
3143
3144 rtpg_clip_arg_destroy(arg);
3145 PG_FREE_IF_COPY(pgraster, 0);
3146 lwgeom_free(geom);
3147 PG_FREE_IF_COPY(gser, 2);
3148
3149 elog(ERROR, "RASTER_clip: Could not get convex hull of raster");
3150 PG_RETURN_NULL();
3151 }
3152
3153 tmpgeom = lwgeom_intersection(rastgeom, geom);
3154 lwgeom_free(rastgeom);
3155 lwgeom_free(geom);
3156 PG_FREE_IF_COPY(gser, 2);
3157 geom = tmpgeom;
3158
3159 /* intersection is empty AND extent type is INTERSECTION, return empty */
3160 if (lwgeom_is_empty(geom) && arg->extenttype == ET_INTERSECTION) {
3161 elog(NOTICE, "The input raster and input geometry do not intersect. Returning empty raster");
3162
3163 rtpg_clip_arg_destroy(arg);
3164 PG_FREE_IF_COPY(pgraster, 0);
3165 lwgeom_free(geom);
3166
3167 rtn = rt_raster_new(0, 0);
3168 if (rtn == NULL) {
3169 elog(ERROR, "RASTER_clip: Could not create empty raster");
3170 PG_RETURN_NULL();
3171 }
3172
3173 pgrtn = rt_raster_serialize(rtn);
3174 rt_raster_destroy(rtn);
3175 if (NULL == pgrtn)
3176 PG_RETURN_NULL();
3177
3178 SET_VARSIZE(pgrtn, pgrtn->size);
3179 PG_RETURN_POINTER(pgrtn);
3180 }
3181
3182 /* nband (1) */
3183 if (!PG_ARGISNULL(1)) {
3184 array = PG_GETARG_ARRAYTYPE_P(1);
3185 etype = ARR_ELEMTYPE(array);
3186 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3187
3188 switch (etype) {
3189 case INT2OID:
3190 case INT4OID:
3191 break;
3192 default:
3193 rtpg_clip_arg_destroy(arg);
3194 PG_FREE_IF_COPY(pgraster, 0);
3195 lwgeom_free(geom);
3196 elog(ERROR, "RASTER_clip: Invalid data type for band indexes");
3197 PG_RETURN_NULL();
3198 break;
3199 }
3200
3201 deconstruct_array(
3202 array, etype,
3203 typlen, typbyval, typalign,
3204 &e, &nulls, &(arg->numbands)
3205 );
3206
3207 arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3208 if (arg->band == NULL) {
3209 rtpg_clip_arg_destroy(arg);
3210 PG_FREE_IF_COPY(pgraster, 0);
3211 lwgeom_free(geom);
3212 elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3213 PG_RETURN_NULL();
3214 }
3215
3216 for (i = 0, j = 0; i < arg->numbands; i++) {
3217 if (nulls[i]) continue;
3218
3219 switch (etype) {
3220 case INT2OID:
3221 arg->band[j].nband = DatumGetInt16(e[i]) - 1;
3222 break;
3223 case INT4OID:
3224 arg->band[j].nband = DatumGetInt32(e[i]) - 1;
3225 break;
3226 }
3227
3228 j++;
3229 }
3230
3231 if (j < arg->numbands) {
3232 arg->band = repalloc(arg->band, sizeof(struct rtpg_clip_band_t) * j);
3233 if (arg->band == NULL) {
3234 rtpg_clip_arg_destroy(arg);
3235 PG_FREE_IF_COPY(pgraster, 0);
3236 lwgeom_free(geom);
3237 elog(ERROR, "RASTER_clip: Could not reallocate memory for band arguments");
3238 PG_RETURN_NULL();
3239 }
3240
3241 arg->numbands = j;
3242 }
3243
3244 /* validate band */
3245 for (i = 0; i < arg->numbands; i++) {
3246 if (!rt_raster_has_band(arg->raster, arg->band[i].nband)) {
3247 elog(NOTICE, "Band at index %d not found in raster", arg->band[i].nband + 1);
3248 rtpg_clip_arg_destroy(arg);
3249 PG_FREE_IF_COPY(pgraster, 0);
3250 lwgeom_free(geom);
3251 PG_RETURN_NULL();
3252 }
3253
3254 arg->band[i].hasnodata = 0;
3255 arg->band[i].nodataval = 0;
3256 }
3257 }
3258 else {
3259 arg->numbands = rt_raster_get_num_bands(arg->raster);
3260
3261 /* raster may have no bands */
3262 if (arg->numbands) {
3263 arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3264 if (arg->band == NULL) {
3265
3266 rtpg_clip_arg_destroy(arg);
3267 PG_FREE_IF_COPY(pgraster, 0);
3268 lwgeom_free(geom);
3269
3270 elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3271 PG_RETURN_NULL();
3272 }
3273
3274 for (i = 0; i < arg->numbands; i++) {
3275 arg->band[i].nband = i;
3276 arg->band[i].hasnodata = 0;
3277 arg->band[i].nodataval = 0;
3278 }
3279 }
3280 }
3281
3282 /* nodataval (3) */
3283 if (!PG_ARGISNULL(3)) {
3284 array = PG_GETARG_ARRAYTYPE_P(3);
3285 etype = ARR_ELEMTYPE(array);
3286 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3287
3288 switch (etype) {
3289 case FLOAT4OID:
3290 case FLOAT8OID:
3291 break;
3292 default:
3293 rtpg_clip_arg_destroy(arg);
3294 PG_FREE_IF_COPY(pgraster, 0);
3295 lwgeom_free(geom);
3296 elog(ERROR, "RASTER_clip: Invalid data type for NODATA values");
3297 PG_RETURN_NULL();
3298 break;
3299 }
3300
3301 deconstruct_array(
3302 array, etype,
3303 typlen, typbyval, typalign,
3304 &e, &nulls, &k
3305 );
3306
3307 /* it doesn't matter if there are more nodataval */
3308 for (i = 0, j = 0; i < arg->numbands; i++, j++) {
3309 /* cap j to the last nodataval element */
3310 if (j >= k)
3311 j = k - 1;
3312
3313 if (nulls[j])
3314 continue;
3315
3316 arg->band[i].hasnodata = 1;
3317 switch (etype) {
3318 case FLOAT4OID:
3319 arg->band[i].nodataval = DatumGetFloat4(e[j]);
3320 break;
3321 case FLOAT8OID:
3322 arg->band[i].nodataval = DatumGetFloat8(e[j]);
3323 break;
3324 }
3325 }
3326 }
3327
3328 /* get wkb of geometry */
3329 POSTGIS_RT_DEBUG(3, "getting wkb of geometry");
3330 wkb = lwgeom_to_wkb(geom, WKB_SFSQL, &wkb_len);
3331 lwgeom_free(geom);
3332
3333 /* rasterize geometry */
3334 arg->mask = rt_raster_gdal_rasterize(
3335 wkb, wkb_len,
3336 NULL,
3337 0, NULL,
3338 NULL, NULL,
3339 NULL, NULL,
3340 NULL, NULL,
3341 &(gt[1]), &(gt[5]),
3342 NULL, NULL,
3343 &(gt[0]), &(gt[3]),
3344 &(gt[2]), &(gt[4]),
3345 NULL
3346 );
3347
3348 pfree(wkb);
3349 if (arg->mask == NULL) {
3350 rtpg_clip_arg_destroy(arg);
3351 PG_FREE_IF_COPY(pgraster, 0);
3352 elog(ERROR, "RASTER_clip: Could not rasterize intersection geometry");
3353 PG_RETURN_NULL();
3354 }
3355
3356 /* set SRID */
3357 rt_raster_set_srid(arg->mask, srid);
3358
3359 /* run iterator */
3360
3361 /* init itrset */
3362 itrset = palloc(sizeof(struct rt_iterator_t) * 2);
3363 if (itrset == NULL) {
3364 rtpg_clip_arg_destroy(arg);
3365 PG_FREE_IF_COPY(pgraster, 0);
3366 elog(ERROR, "RASTER_clip: Could not allocate memory for iterator arguments");
3367 PG_RETURN_NULL();
3368 }
3369
3370 /* one band at a time */
3371 for (i = 0; i < arg->numbands; i++) {
3372 POSTGIS_RT_DEBUGF(4, "band arg %d (nband, hasnodata, nodataval) = (%d, %d, %f)",
3373 i, arg->band[i].nband, arg->band[i].hasnodata, arg->band[i].nodataval);
3374
3375 band = rt_raster_get_band(arg->raster, arg->band[i].nband);
3376
3377 /* band metadata */
3378 pixtype = rt_band_get_pixtype(band);
3379
3380 if (arg->band[i].hasnodata) {
3381 hasnodata = 1;
3382 nodataval = arg->band[i].nodataval;
3383 }
3384 else if (rt_band_get_hasnodata_flag(band)) {
3385 hasnodata = 1;
3386 rt_band_get_nodata(band, &nodataval);
3387 }
3388 else {
3389 hasnodata = 0;
3390 nodataval = rt_band_get_min_value(band);
3391 }
3392
3393 /* band is NODATA, create NODATA band and continue */
3394 if (rt_band_get_isnodata_flag(band)) {
3395 /* create raster */
3396 if (rtn == NULL) {
3397 noerr = rt_raster_from_two_rasters(arg->raster, arg->mask, arg->extenttype, &rtn, NULL);
3398 if (noerr != ES_NONE) {
3399 rtpg_clip_arg_destroy(arg);
3400 PG_FREE_IF_COPY(pgraster, 0);
3401 elog(ERROR, "RASTER_clip: Could not create output raster");
3402 PG_RETURN_NULL();
3403 }
3404 }
3405
3406 /* create NODATA band */
3407 if (rt_raster_generate_new_band(rtn, pixtype, nodataval, hasnodata, nodataval, i) < 0) {
3408 rt_raster_destroy(rtn);
3409 rtpg_clip_arg_destroy(arg);
3410 PG_FREE_IF_COPY(pgraster, 0);
3411 elog(ERROR, "RASTER_clip: Could not add NODATA band to output raster");
3412 PG_RETURN_NULL();
3413 }
3414
3415 continue;
3416 }
3417
3418 /* raster */
3419 itrset[0].raster = arg->raster;
3420 itrset[0].nband = arg->band[i].nband;
3421 itrset[0].nbnodata = 1;
3422
3423 /* mask */
3424 itrset[1].raster = arg->mask;
3425 itrset[1].nband = 0;
3426 itrset[1].nbnodata = 1;
3427
3428 /* pass to iterator */
3429 noerr = rt_raster_iterator(
3430 itrset, 2,
3431 arg->extenttype, NULL,
3432 pixtype,
3433 hasnodata, nodataval,
3434 0, 0,
3435 NULL,
3436 NULL,
3437 rtpg_clip_callback,
3438 &_raster
3439 );
3440
3441 if (noerr != ES_NONE) {
3442 pfree(itrset);
3443 rtpg_clip_arg_destroy(arg);
3444 if (rtn != NULL) rt_raster_destroy(rtn);
3445 PG_FREE_IF_COPY(pgraster, 0);
3446 elog(ERROR, "RASTER_clip: Could not run raster iterator function");
3447 PG_RETURN_NULL();
3448 }
3449
3450 /* new raster */
3451 if (rtn == NULL)
3452 rtn = _raster;
3453 /* copy band */
3454 else {
3455 band = rt_raster_get_band(_raster, 0);
3456 if (band == NULL) {
3457 pfree(itrset);
3458 rtpg_clip_arg_destroy(arg);
3459 rt_raster_destroy(_raster);
3460 rt_raster_destroy(rtn);
3461 PG_FREE_IF_COPY(pgraster, 0);
3462 elog(NOTICE, "RASTER_clip: Could not get band from working raster");
3463 PG_RETURN_NULL();
3464 }
3465
3466 if (rt_raster_add_band(rtn, band, i) < 0) {
3467 pfree(itrset);
3468 rtpg_clip_arg_destroy(arg);
3469 rt_raster_destroy(_raster);
3470 rt_raster_destroy(rtn);
3471 PG_FREE_IF_COPY(pgraster, 0);
3472 elog(ERROR, "RASTER_clip: Could not add new band to output raster");
3473 PG_RETURN_NULL();
3474 }
3475
3476 rt_raster_destroy(_raster);
3477 }
3478 }
3479
3480 pfree(itrset);
3481 rtpg_clip_arg_destroy(arg);
3482 PG_FREE_IF_COPY(pgraster, 0);
3483
3484 pgrtn = rt_raster_serialize(rtn);
3485 rt_raster_destroy(rtn);
3486
3487 POSTGIS_RT_DEBUG(3, "Finished");
3488
3489 if (!pgrtn)
3490 PG_RETURN_NULL();
3491
3492 SET_VARSIZE(pgrtn, pgrtn->size);
3493 PG_RETURN_POINTER(pgrtn);
3494 }
3495
3496 /**
3497 * Reclassify the specified bands of the raster
3498 */
3499 PG_FUNCTION_INFO_V1(RASTER_reclass);
3500 Datum RASTER_reclass(PG_FUNCTION_ARGS) {
3501 rt_pgraster *pgraster = NULL;
3502 rt_pgraster *pgrtn = NULL;
3503 rt_raster raster = NULL;
3504 rt_band band = NULL;
3505 rt_band newband = NULL;
3506 uint32_t numBands = 0;
3507
3508 ArrayType *array;
3509 Oid etype;
3510 Datum *e;
3511 bool *nulls;
3512 int16 typlen;
3513 bool typbyval;
3514 char typalign;
3515 int n = 0;
3516
3517 int i = 0;
3518 int j = 0;
3519 int k = 0;
3520
3521 uint32_t a = 0;
3522 uint32_t b = 0;
3523 uint32_t c = 0;
3524
3525 rt_reclassexpr *exprset = NULL;
3526 HeapTupleHeader tup;
3527 bool isnull;
3528 Datum tupv;
3529 uint32_t nband = 0;
3530 char *expr = NULL;
3531 text *exprtext = NULL;
3532 double val = 0;
3533 char *junk = NULL;
3534 int inc_val = 0;
3535 int exc_val = 0;
3536 char *pixeltype = NULL;
3537 text *pixeltypetext = NULL;
3538 rt_pixtype pixtype = PT_END;
3539 double nodataval = 0;
3540 bool hasnodata = FALSE;
3541
3542 char **comma_set = NULL;
3543 uint32_t comma_n = 0;
3544 char **colon_set = NULL;
3545 uint32_t colon_n = 0;
3546 char **dash_set = NULL;
3547 uint32_t dash_n = 0;
3548
3549 POSTGIS_RT_DEBUG(3, "RASTER_reclass: Starting");
3550
3551 /* pgraster is null, return null */
3552 if (PG_ARGISNULL(0))
3553 PG_RETURN_NULL();
3554 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3555
3556 /* raster */
3557 raster = rt_raster_deserialize(pgraster, FALSE);
3558 if (!raster) {
3559 PG_FREE_IF_COPY(pgraster, 0);
3560 elog(ERROR, "RASTER_reclass: Could not deserialize raster");
3561 PG_RETURN_NULL();
3562 }
3563 numBands = rt_raster_get_num_bands(raster);
3564 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: %d possible bands to be reclassified", numBands);
3565
3566 /* process set of reclassarg */
3567 POSTGIS_RT_DEBUG(3, "RASTER_reclass: Processing Arg 1 (reclassargset)");
3568 array = PG_GETARG_ARRAYTYPE_P(1);
3569 etype = ARR_ELEMTYPE(array);
3570 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3571
3572 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
3573 &nulls, &n);
3574
3575 if (!n) {
3576 elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
3577
3578 pgrtn = rt_raster_serialize(raster);
3579 rt_raster_destroy(raster);
3580 PG_FREE_IF_COPY(pgraster, 0);
3581 if (!pgrtn)
3582 PG_RETURN_NULL();
3583
3584 SET_VARSIZE(pgrtn, pgrtn->size);
3585 PG_RETURN_POINTER(pgrtn);
3586 }
3587
3588 /*
3589 process each element of reclassarg
3590 each element is the index of the band to process, the set
3591 of reclass ranges and the output band's pixeltype
3592 */
3593 for (i = 0; i < n; i++) {
3594 if (nulls[i]) continue;
3595
3596 /* each element is a tuple */
3597 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
3598 if (NULL == tup) {
3599 elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
3600
3601 pgrtn = rt_raster_serialize(raster);
3602 rt_raster_destroy(raster);
3603 PG_FREE_IF_COPY(pgraster, 0);
3604 if (!pgrtn)
3605 PG_RETURN_NULL();
3606
3607 SET_VARSIZE(pgrtn, pgrtn->size);
3608 PG_RETURN_POINTER(pgrtn);
3609 }
3610
3611 /* band index (1-based) */
3612 tupv = GetAttributeByName(tup, "nband", &isnull);
3613 if (isnull) {
3614 elog(NOTICE, "Invalid argument for reclassargset. Missing value of nband for reclassarg of index %d . Returning original raster", i);
3615
3616 pgrtn = rt_raster_serialize(raster);
3617 rt_raster_destroy(raster);
3618 PG_FREE_IF_COPY(pgraster, 0);
3619 if (!pgrtn)
3620 PG_RETURN_NULL();
3621
3622 SET_VARSIZE(pgrtn, pgrtn->size);
3623 PG_RETURN_POINTER(pgrtn);
3624 }
3625 nband = DatumGetInt32(tupv);
3626 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: expression for band %d", nband);
3627
3628 /* valid band index? */
3629 if (nband < 1 || nband > numBands) {
3630 elog(NOTICE, "Invalid argument for reclassargset. Invalid band index (must use 1-based) for reclassarg of index %d . Returning original raster", i);
3631
3632 pgrtn = rt_raster_serialize(raster);
3633 rt_raster_destroy(raster);
3634 PG_FREE_IF_COPY(pgraster, 0);
3635 if (!pgrtn)
3636 PG_RETURN_NULL();
3637
3638 SET_VARSIZE(pgrtn, pgrtn->size);
3639 PG_RETURN_POINTER(pgrtn);
3640 }
3641
3642 /* reclass expr */
3643 tupv = GetAttributeByName(tup, "reclassexpr", &isnull);
3644 if (isnull) {
3645 elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3646
3647 pgrtn = rt_raster_serialize(raster);
3648 rt_raster_destroy(raster);
3649 PG_FREE_IF_COPY(pgraster, 0);
3650 if (!pgrtn)
3651 PG_RETURN_NULL();
3652
3653 SET_VARSIZE(pgrtn, pgrtn->size);
3654 PG_RETURN_POINTER(pgrtn);
3655 }
3656 exprtext = (text *) DatumGetPointer(tupv);
3657 if (NULL == exprtext) {
3658 elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3659
3660 pgrtn = rt_raster_serialize(raster);
3661 rt_raster_destroy(raster);
3662 PG_FREE_IF_COPY(pgraster, 0);
3663 if (!pgrtn)
3664 PG_RETURN_NULL();
3665
3666 SET_VARSIZE(pgrtn, pgrtn->size);
3667 PG_RETURN_POINTER(pgrtn);
3668 }
3669 expr = text_to_cstring(exprtext);
3670 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: expr (raw) %s", expr);
3671 expr = rtpg_removespaces(expr);
3672 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: expr (clean) %s", expr);
3673
3674 /* split string to its components */
3675 /* comma (,) separating rangesets */
3676 comma_set = rtpg_strsplit(expr, ",", &comma_n);
3677 if (comma_n < 1) {
3678 elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3679
3680 pgrtn = rt_raster_serialize(raster);
3681 rt_raster_destroy(raster);
3682 PG_FREE_IF_COPY(pgraster, 0);
3683 if (!pgrtn)
3684 PG_RETURN_NULL();
3685
3686 SET_VARSIZE(pgrtn, pgrtn->size);
3687 PG_RETURN_POINTER(pgrtn);
3688 }
3689
3690 /* set of reclass expressions */
3691 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: %d possible expressions", comma_n);
3692 exprset = palloc(comma_n * sizeof(rt_reclassexpr));
3693
3694 for (a = 0, j = 0; a < comma_n; a++) {
3695 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: map %s", comma_set[a]);
3696
3697 /* colon (:) separating range "src" and "dst" */
3698 colon_set = rtpg_strsplit(comma_set[a], ":", &colon_n);
3699 if (colon_n != 2) {
3700 elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3701 for (k = 0; k < j; k++) pfree(exprset[k]);
3702 pfree(exprset);
3703
3704 pgrtn = rt_raster_serialize(raster);
3705 rt_raster_destroy(raster);
3706 PG_FREE_IF_COPY(pgraster, 0);
3707 if (!pgrtn)
3708 PG_RETURN_NULL();
3709
3710 SET_VARSIZE(pgrtn, pgrtn->size);
3711 PG_RETURN_POINTER(pgrtn);
3712 }
3713
3714 /* allocate mem for reclass expression */
3715 exprset[j] = palloc(sizeof(struct rt_reclassexpr_t));
3716
3717 for (b = 0; b < colon_n; b++) {
3718 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: range %s", colon_set[b]);
3719
3720 /* dash (-) separating "min" and "max" */
3721 dash_set = rtpg_strsplit(colon_set[b], "-", &dash_n);
3722 if (dash_n < 1 || dash_n > 3) {
3723 elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3724 for (k = 0; k < j; k++) pfree(exprset[k]);
3725 pfree(exprset);
3726
3727 pgrtn = rt_raster_serialize(raster);
3728 rt_raster_destroy(raster);
3729 PG_FREE_IF_COPY(pgraster, 0);
3730 if (!pgrtn)
3731 PG_RETURN_NULL();
3732
3733 SET_VARSIZE(pgrtn, pgrtn->size);
3734 PG_RETURN_POINTER(pgrtn);
3735 }
3736
3737 for (c = 0; c < dash_n; c++) {
3738 /* need to handle: (-9999-100 -> "(", "9999", "100" */
3739 if (
3740 c < 1 &&
3741 strlen(dash_set[c]) == 1 && (
3742 strchr(dash_set[c], '(') != NULL ||
3743 strchr(dash_set[c], '[') != NULL ||
3744 strchr(dash_set[c], ')') != NULL ||
3745 strchr(dash_set[c], ']') != NULL
3746 )
3747 ) {
3748 uint32_t dash_it;
3749 junk = palloc(sizeof(char) * (strlen(dash_set[c + 1]) + 2));
3750 if (NULL == junk) {
3751 for (k = 0; k <= j; k++) pfree(exprset[k]);
3752 pfree(exprset);
3753 rt_raster_destroy(raster);
3754 PG_FREE_IF_COPY(pgraster, 0);
3755
3756 elog(ERROR, "RASTER_reclass: Could not allocate memory");
3757 PG_RETURN_NULL();
3758 }
3759
3760 sprintf(junk, "%s%s", dash_set[c], dash_set[c + 1]);
3761 c++;
3762 dash_set[c] = repalloc(dash_set[c], sizeof(char) * (strlen(junk) + 1));
3763 strcpy(dash_set[c], junk);
3764 pfree(junk);
3765
3766 /* rebuild dash_set */
3767 for (dash_it = 1; dash_it < dash_n; dash_it++) {
3768 dash_set[dash_it - 1] = repalloc(dash_set[dash_it - 1], (strlen(dash_set[dash_it]) + 1) * sizeof(char));
3769 strcpy(dash_set[dash_it - 1], dash_set[dash_it]);
3770 }
3771 dash_n--;
3772 c--;
3773 pfree(dash_set[dash_n]);
3774 dash_set = repalloc(dash_set, sizeof(char *) * dash_n);
3775 }
3776
3777 /* there shouldn't be more than two in dash_n */
3778 if (c < 1 && dash_n > 2) {
3779 elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3780 for (k = 0; k < j; k++) pfree(exprset[k]);
3781 pfree(exprset);
3782
3783 pgrtn = rt_raster_serialize(raster);
3784 rt_raster_destroy(raster);
3785 PG_FREE_IF_COPY(pgraster, 0);
3786 if (!pgrtn)
3787 PG_RETURN_NULL();
3788
3789 SET_VARSIZE(pgrtn, pgrtn->size);
3790 PG_RETURN_POINTER(pgrtn);
3791 }
3792
3793 /* check interval flags */
3794 exc_val = 0;
3795 inc_val = 1;
3796 /* range */
3797 if (dash_n != 1) {
3798 /* min */
3799 if (c < 1) {
3800 if (
3801 strchr(dash_set[c], ')') != NULL ||
3802 strchr(dash_set[c], ']') != NULL
3803 ) {
3804 exc_val = 1;
3805 inc_val = 1;
3806 }
3807 else if (strchr(dash_set[c], '(') != NULL){
3808 inc_val = 0;
3809 }
3810 else {
3811 inc_val = 1;
3812 }
3813 }
3814 /* max */
3815 else {
3816 if (
3817 strrchr(dash_set[c], '(') != NULL ||
3818 strrchr(dash_set[c], '[') != NULL
3819 ) {
3820 exc_val = 1;
3821 inc_val = 0;
3822 }
3823 else if (strrchr(dash_set[c], ']') != NULL) {
3824 inc_val = 1;
3825 }
3826 else {
3827 inc_val = 0;
3828 }
3829 }
3830 }
3831 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: exc_val %d inc_val %d", exc_val, inc_val);
3832
3833 /* remove interval flags */
3834 dash_set[c] = rtpg_chartrim(dash_set[c], "()[]");
3835 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (char) %s", dash_set[c]);
3836
3837 /* value from string to double */
3838 errno = 0;
3839 val = strtod(dash_set[c], &junk);
3840 if (errno != 0 || dash_set[c] == junk) {
3841 elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3842 for (k = 0; k < j; k++) pfree(exprset[k]);
3843 pfree(exprset);
3844
3845 pgrtn = rt_raster_serialize(raster);
3846 rt_raster_destroy(raster);
3847 PG_FREE_IF_COPY(pgraster, 0);
3848 if (!pgrtn)
3849 PG_RETURN_NULL();
3850
3851 SET_VARSIZE(pgrtn, pgrtn->size);
3852 PG_RETURN_POINTER(pgrtn);
3853 }
3854 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (double) %f", val);
3855
3856 /* strsplit removes dash (a.k.a. negative sign), compare now to restore */
3857 if (c < 1)
3858 junk = strstr(colon_set[b], dash_set[c]);
3859 else
3860 junk = rtpg_strrstr(colon_set[b], dash_set[c]);
3861 POSTGIS_RT_DEBUGF(
3862 4,
3863 "(colon_set[%d], dash_set[%d], junk) = (%s, %s, %s)",
3864 b, c, colon_set[b], dash_set[c], junk
3865 );
3866 /* not beginning of string */
3867 if (junk != colon_set[b]) {
3868 /* prior is a dash */
3869 if (*(junk - 1) == '-') {
3870 /* prior is beginning of string or prior - 1 char is dash, negative number */
3871 if (
3872 ((junk - 1) == colon_set[b]) ||
3873 (*(junk - 2) == '-') ||
3874 (*(junk - 2) == '[') ||
3875 (*(junk - 2) == '(')
3876 ) {
3877 val *= -1.;
3878 }
3879 }
3880 }
3881 POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (double) %f", val);
3882
3883 /* src */
3884 if (b < 1) {
3885 /* singular value */
3886 if (dash_n == 1) {
3887 exprset[j]->src.exc_min = exprset[j]->src.exc_max = exc_val;
3888 exprset[j]->src.inc_min = exprset[j]->src.inc_max = inc_val;
3889 exprset[j]->src.min = exprset[j]->src.max = val;
3890 }
3891 /* min */
3892 else if (c < 1) {
3893 exprset[j]->src.exc_min = exc_val;
3894 exprset[j]->src.inc_min = inc_val;
3895 exprset[j]->src.min = val;
3896 }
3897 /* max */
3898 else {
3899 exprset[j]->src.exc_max = exc_val;
3900 exprset[j]->src.inc_max = inc_val;
3901 exprset[j]->src.max = val;
3902 }
3903 }
3904 /* dst */
3905 else {
3906 /* singular value */
3907 if (dash_n == 1)
3908 exprset[j]->dst.min = exprset[j]->dst.max = val;
3909 /* min */
3910 else if (c < 1)
3911 exprset[j]->dst.min = val;
3912 /* max */
3913 else
3914 exprset[j]->dst.max = val;
3915 }
3916 }
3917 pfree(dash_set);
3918 }
3919 pfree(colon_set);
3920
3921 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: or: %f - %f nr: %f - %f"
3922 , exprset[j]->src.min
3923 , exprset[j]->src.max
3924 , exprset[j]->dst.min
3925 , exprset[j]->dst.max
3926 );
3927 j++;
3928 }
3929 pfree(comma_set);
3930
3931 /* pixel type */
3932 tupv = GetAttributeByName(tup, "pixeltype", &isnull);
3933 if (isnull) {
3934 elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3935
3936 pgrtn = rt_raster_serialize(raster);
3937 rt_raster_destroy(raster);
3938 PG_FREE_IF_COPY(pgraster, 0);
3939 if (!pgrtn)
3940 PG_RETURN_NULL();
3941
3942 SET_VARSIZE(pgrtn, pgrtn->size);
3943 PG_RETURN_POINTER(pgrtn);
3944 }
3945 pixeltypetext = (text *) DatumGetPointer(tupv);
3946 if (NULL == pixeltypetext) {
3947 elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3948
3949 pgrtn = rt_raster_serialize(raster);
3950 rt_raster_destroy(raster);
3951 PG_FREE_IF_COPY(pgraster, 0);
3952 if (!pgrtn)
3953 PG_RETURN_NULL();
3954
3955 SET_VARSIZE(pgrtn, pgrtn->size);
3956 PG_RETURN_POINTER(pgrtn);
3957 }
3958 pixeltype = text_to_cstring(pixeltypetext);
3959 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: pixeltype %s", pixeltype);
3960 pixtype = rt_pixtype_index_from_name(pixeltype);
3961
3962 /* nodata */
3963 tupv = GetAttributeByName(tup, "nodataval", &isnull);
3964 if (isnull) {
3965 nodataval = 0;
3966 hasnodata = FALSE;
3967 }
3968 else {
3969 nodataval = DatumGetFloat8(tupv);
3970 hasnodata = TRUE;
3971 }
3972 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: nodataval %f", nodataval);
3973 POSTGIS_RT_DEBUGF(3, "RASTER_reclass: hasnodata %d", hasnodata);
3974
3975 /* do reclass */
3976 band = rt_raster_get_band(raster, nband - 1);
3977 if (!band) {
3978 elog(NOTICE, "Could not find raster band of index %d. Returning original raster", nband);
3979 for (k = 0; k < j; k++) pfree(exprset[k]);
3980 pfree(exprset);
3981
3982 pgrtn = rt_raster_serialize(raster);
3983 rt_raster_destroy(raster);
3984 PG_FREE_IF_COPY(pgraster, 0);
3985 if (!pgrtn)
3986 PG_RETURN_NULL();
3987
3988 SET_VARSIZE(pgrtn, pgrtn->size);
3989 PG_RETURN_POINTER(pgrtn);
3990 }
3991 newband = rt_band_reclass(band, pixtype, hasnodata, nodataval, exprset, j);
3992 if (!newband) {
3993 for (k = 0; k < j; k++) pfree(exprset[k]);
3994 pfree(exprset);
3995
3996 rt_raster_destroy(raster);
3997 PG_FREE_IF_COPY(pgraster, 0);
3998
3999 elog(ERROR, "RASTER_reclass: Could not reclassify raster band of index %d", nband);
4000 PG_RETURN_NULL();
4001 }
4002
4003 /* replace old band with new band */
4004 if (rt_raster_replace_band(raster, newband, nband - 1) == NULL) {
4005 for (k = 0; k < j; k++) pfree(exprset[k]);
4006 pfree(exprset);
4007
4008 rt_band_destroy(newband);
4009 rt_raster_destroy(raster);
4010 PG_FREE_IF_COPY(pgraster, 0);
4011
4012 elog(ERROR, "RASTER_reclass: Could not replace raster band of index %d with reclassified band", nband);
4013 PG_RETURN_NULL();
4014 }
4015
4016 /* old band is in the variable band */
4017 rt_band_destroy(band);
4018
4019 /* free exprset */
4020 for (k = 0; k < j; k++) pfree(exprset[k]);
4021 pfree(exprset);
4022 }
4023
4024 pgrtn = rt_raster_serialize(raster);
4025 rt_raster_destroy(raster);
4026 PG_FREE_IF_COPY(pgraster, 0);
4027 if (!pgrtn)
4028 PG_RETURN_NULL();
4029
4030 POSTGIS_RT_DEBUG(3, "RASTER_reclass: Finished");
4031
4032 SET_VARSIZE(pgrtn, pgrtn->size);
4033 PG_RETURN_POINTER(pgrtn);
4034 }
4035
4036 /* ---------------------------------------------------------------- */
4037 /* apply colormap to specified band of a raster */
4038 /* ---------------------------------------------------------------- */
4039
4040 typedef struct rtpg_colormap_arg_t *rtpg_colormap_arg;
4041 struct rtpg_colormap_arg_t {
4042 rt_raster raster;
4043 int nband; /* 1-based */
4044 rt_band band;
4045 rt_bandstats bandstats;
4046
4047 rt_colormap colormap;
4048 int nodataentry;
4049
4050 char **entry;
4051 uint32_t nentry;
4052 char **element;
4053 uint32_t nelement;
4054 };
4055
4056 static rtpg_colormap_arg
4057 rtpg_colormap_arg_init() {
4058 rtpg_colormap_arg arg = NULL;
4059
4060 arg = palloc(sizeof(struct rtpg_colormap_arg_t));
4061 if (arg == NULL) {
4062 elog(ERROR, "rtpg_colormap_arg: Could not allocate memory for function arguments");
4063 return NULL;
4064 }
4065
4066 arg->raster = NULL;
4067 arg->nband = 1;
4068 arg->band = NULL;
4069 arg->bandstats = NULL;
4070
4071 arg->colormap = palloc(sizeof(struct rt_colormap_t));
4072 if (arg->colormap == NULL) {
4073 elog(ERROR, "rtpg_colormap_arg: Could not allocate memory for function arguments");
4074 return NULL;
4075 }
4076 arg->colormap->nentry = 0;
4077 arg->colormap->entry = NULL;
4078 arg->colormap->ncolor = 4; /* assume RGBA */
4079 arg->colormap->method = CM_INTERPOLATE;
4080 arg->nodataentry = -1;
4081
4082 arg->entry = NULL;
4083 arg->nentry = 0;
4084 arg->element = NULL;
4085 arg->nelement = 0;
4086
4087 return arg;
4088 }
4089
4090 static void
4091 rtpg_colormap_arg_destroy(rtpg_colormap_arg arg) {
4092 uint32_t i = 0;
4093 if (arg->raster != NULL)
4094 rt_raster_destroy(arg->raster);
4095
4096 if (arg->bandstats != NULL)
4097 pfree(arg->bandstats);
4098
4099 if (arg->colormap != NULL) {
4100 if (arg->colormap->entry != NULL)
4101 pfree(arg->colormap->entry);
4102 pfree(arg->colormap);
4103 }
4104
4105 if (arg->nentry) {
4106 for (i = 0; i < arg->nentry; i++) {
4107 if (arg->entry[i] != NULL)
4108 pfree(arg->entry[i]);
4109 }
4110 pfree(arg->entry);
4111 }
4112
4113 if (arg->nelement) {
4114 for (i = 0; i < arg->nelement; i++)
4115 pfree(arg->element[i]);
4116 pfree(arg->element);
4117 }
4118
4119 pfree(arg);
4120 arg = NULL;
4121 }
4122
4123 PG_FUNCTION_INFO_V1(RASTER_colorMap);
4124 Datum RASTER_colorMap(PG_FUNCTION_ARGS)
4125 {
4126 rt_pgraster *pgraster = NULL;
4127 rtpg_colormap_arg arg = NULL;
4128 char *junk = NULL;
4129 rt_raster raster = NULL;
4130
4131 POSTGIS_RT_DEBUG(3, "RASTER_colorMap: Starting");
4132
4133 /* pgraster is NULL, return NULL */
4134 if (PG_ARGISNULL(0))
4135 PG_RETURN_NULL();
4136
4137 /* init arg */
4138 arg = rtpg_colormap_arg_init();
4139 if (arg == NULL) {
4140 elog(ERROR, "RASTER_colorMap: Could not initialize argument structure");
4141 PG_RETURN_NULL();
4142 }
4143
4144 /* raster (0) */
4145 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4146
4147 /* raster object */
4148 arg->raster = rt_raster_deserialize(pgraster, FALSE);
4149 if (!arg->raster) {
4150 rtpg_colormap_arg_destroy(arg);
4151 PG_FREE_IF_COPY(pgraster, 0);
4152 elog(ERROR, "RASTER_colorMap: Could not deserialize raster");
4153 PG_RETURN_NULL();
4154 }
4155
4156 /* nband (1) */
4157 if (!PG_ARGISNULL(1))
4158 arg->nband = PG_GETARG_INT32(1);
4159 POSTGIS_RT_DEBUGF(4, "nband = %d", arg->nband);
4160
4161 /* check that band works */
4162 if (!rt_raster_has_band(arg->raster, arg->nband - 1)) {
4163 elog(NOTICE, "Raster does not have band at index %d. Returning empty raster", arg->nband);
4164
4165 raster = rt_raster_clone(arg->raster, 0);
4166 if (raster == NULL) {
4167 rtpg_colormap_arg_destroy(arg);
4168 PG_FREE_IF_COPY(pgraster, 0);
4169 elog(ERROR, "RASTER_colorMap: Could not create empty raster");
4170 PG_RETURN_NULL();
4171 }
4172
4173 rtpg_colormap_arg_destroy(arg);
4174 PG_FREE_IF_COPY(pgraster, 0);
4175
4176 pgraster = rt_raster_serialize(raster);
4177 rt_raster_destroy(raster);
4178 if (pgraster == NULL)
4179 PG_RETURN_NULL();
4180
4181 SET_VARSIZE(pgraster, ((rt_pgraster*) pgraster)->size);
4182 PG_RETURN_POINTER(pgraster);
4183 }
4184
4185 /* get band */
4186 arg->band = rt_raster_get_band(arg->raster, arg->nband - 1);
4187 if (arg->band == NULL) {
4188 int nband = arg->nband;
4189 rtpg_colormap_arg_destroy(arg);
4190 PG_FREE_IF_COPY(pgraster, 0);
4191 elog(ERROR, "RASTER_colorMap: Could not get band at index %d", nband);
4192 PG_RETURN_NULL();
4193 }
4194
4195 /* method (3) */
4196 if (!PG_ARGISNULL(3)) {
4197 char *method = NULL;
4198 char *tmp = text_to_cstring(PG_GETARG_TEXT_P(3));
4199 POSTGIS_RT_DEBUGF(4, "raw method = %s", tmp);
4200
4201 method = rtpg_trim(tmp);
4202 pfree(tmp);
4203 method = rtpg_strtoupper(method);
4204
4205 if (strcmp(method, "INTERPOLATE") == 0)
4206 arg->colormap->method = CM_INTERPOLATE;
4207 else if (strcmp(method, "EXACT") == 0)
4208 arg->colormap->method = CM_EXACT;
4209 else if (strcmp(method, "NEAREST") == 0)
4210 arg->colormap->method = CM_NEAREST;
4211 else {
4212 elog(NOTICE, "Unknown value provided for method. Defaulting to INTERPOLATE");
4213 arg->colormap->method = CM_INTERPOLATE;
4214 }
4215 }
4216 /* default to INTERPOLATE */
4217 else
4218 arg->colormap->method = CM_INTERPOLATE;
4219 POSTGIS_RT_DEBUGF(4, "method = %d", arg->colormap->method);
4220
4221 /* colormap (2) */
4222 if (PG_ARGISNULL(2)) {
4223 rtpg_colormap_arg_destroy(arg);
4224 PG_FREE_IF_COPY(pgraster, 0);
4225 elog(ERROR, "RASTER_colorMap: Value must be provided for colormap");
4226 PG_RETURN_NULL();
4227 }
4228 else {
4229 char *tmp = NULL;
4230 char *colormap = text_to_cstring(PG_GETARG_TEXT_P(2));
4231 char *_entry;
4232 char *_element;
4233 uint32_t i = 0;
4234 uint32_t j = 0;
4235
4236 POSTGIS_RT_DEBUGF(4, "colormap = %s", colormap);
4237
4238 /* empty string */
4239 if (!strlen(colormap)) {
4240 rtpg_colormap_arg_destroy(arg);
4241 PG_FREE_IF_COPY(pgraster, 0);
4242 elog(ERROR, "RASTER_colorMap: Value must be provided for colormap");
4243 PG_RETURN_NULL();
4244 }
4245
4246 arg->entry = rtpg_strsplit(colormap, "\n", &(arg->nentry));
4247 pfree(colormap);
4248 if (arg->nentry < 1) {
4249 rtpg_colormap_arg_destroy(arg);
4250 PG_FREE_IF_COPY(pgraster, 0);
4251 elog(ERROR, "RASTER_colorMap: Could not process the value provided for colormap");
4252 PG_RETURN_NULL();
4253 }
4254
4255 /* allocate the max # of colormap entries */
4256 arg->colormap->entry = palloc(sizeof(struct rt_colormap_entry_t) * arg->nentry);
4257 if (arg->colormap->entry == NULL) {
4258 rtpg_colormap_arg_destroy(arg);
4259 PG_FREE_IF_COPY(pgraster, 0);
4260 elog(ERROR, "RASTER_colorMap: Could not allocate memory for colormap entries");
4261 PG_RETURN_NULL();
4262 }
4263 memset(arg->colormap->entry, 0, sizeof(struct rt_colormap_entry_t) * arg->nentry);
4264
4265 /* each entry */
4266 for (i = 0; i < arg->nentry; i++) {
4267 /* substitute space for other delimiters */
4268 tmp = rtpg_strreplace(arg->entry[i], ":", " ", NULL);
4269 _entry = rtpg_strreplace(tmp, ",", " ", NULL);
4270 pfree(tmp);
4271 tmp = rtpg_strreplace(_entry, "\t", " ", NULL);
4272 pfree(_entry);
4273 _entry = rtpg_trim(tmp);
4274 pfree(tmp);
4275
4276 POSTGIS_RT_DEBUGF(4, "Processing entry[%d] = %s", i, arg->entry[i]);
4277 POSTGIS_RT_DEBUGF(4, "Cleaned entry[%d] = %s", i, _entry);
4278
4279 /* empty entry, continue */
4280 if (!strlen(_entry)) {
4281 POSTGIS_RT_DEBUGF(3, "Skipping empty entry[%d]", i);
4282 pfree(_entry);
4283 continue;
4284 }
4285
4286 arg->element = rtpg_strsplit(_entry, " ", &(arg->nelement));
4287 pfree(_entry);
4288 if (arg->nelement < 2) {
4289 rtpg_colormap_arg_destroy(arg);
4290 PG_FREE_IF_COPY(pgraster, 0);
4291 elog(ERROR, "RASTER_colorMap: Could not process colormap entry %d", i + 1);
4292 PG_RETURN_NULL();
4293 }
4294 else if (arg->nelement > 5) {
4295 elog(NOTICE, "More than five elements in colormap entry %d. Using at most five elements", i + 1);
4296 arg->nelement = 5;
4297 }
4298
4299 /* smallest # of colors */
4300 if (((int)arg->nelement - 1) < arg->colormap->ncolor)
4301 arg->colormap->ncolor = arg->nelement - 1;
4302
4303 /* each element of entry */
4304 for (j = 0; j < arg->nelement; j++) {
4305
4306 _element = rtpg_trim(arg->element[j]);
4307 _element = rtpg_strtoupper(_element);
4308 POSTGIS_RT_DEBUGF(4, "Processing entry[%d][%d] = %s", i, j, arg->element[j]);
4309 POSTGIS_RT_DEBUGF(4, "Cleaned entry[%d][%d] = %s", i, j, _element);
4310
4311 /* first element is ALWAYS a band value, percentage OR "nv" string */
4312 if (j == 0) {
4313 char *percent = NULL;
4314
4315 /* NODATA */
4316 if (
4317 strcmp(_element, "NV") == 0 ||
4318 strcmp(_element, "NULL") == 0 ||
4319 strcmp(_element, "NODATA") == 0
4320 ) {
4321 POSTGIS_RT_DEBUG(4, "Processing NODATA string");
4322
4323 if (arg->nodataentry > -1) {
4324 elog(NOTICE, "More than one NODATA entry found. Using only the first one");
4325 }
4326 else {
4327 arg->colormap->entry[arg->colormap->nentry].isnodata = 1;
4328 /* no need to set value as value comes from band's NODATA */
4329 arg->colormap->entry[arg->colormap->nentry].value = 0;
4330 }
4331 }
4332 /* percent value */
4333 else if ((percent = strchr(_element, '%')) != NULL) {
4334 double value;
4335 POSTGIS_RT_DEBUG(4, "Processing percent string");
4336
4337 /* get the band stats */
4338 if (arg->bandstats == NULL) {
4339 POSTGIS_RT_DEBUG(4, "Getting band stats");
4340
4341 arg->bandstats = rt_band_get_summary_stats(arg->band, 1, 1, 0, NULL, NULL, NULL);
4342 if (arg->bandstats == NULL) {
4343 pfree(_element);
4344 rtpg_colormap_arg_destroy(arg);
4345 PG_FREE_IF_COPY(pgraster, 0);
4346 elog(ERROR, "RASTER_colorMap: Could not get band's summary stats to process percentages");
4347 PG_RETURN_NULL();
4348 }
4349 }
4350
4351 /* get the string before the percent char */
4352 tmp = palloc(sizeof(char) * (percent - _element + 1));
4353 if (tmp == NULL) {
4354 pfree(_element);
4355 rtpg_colormap_arg_destroy(arg);
4356 PG_FREE_IF_COPY(pgraster, 0);
4357 elog(ERROR, "RASTER_colorMap: Could not allocate memory for value of percentage");
4358 PG_RETURN_NULL();
4359 }
4360
4361 memcpy(tmp, _element, percent - _element);
4362 tmp[percent - _element] = '\0';
4363 POSTGIS_RT_DEBUGF(4, "Percent value = %s", tmp);
4364
4365 /* get percentage value */
4366 errno = 0;
4367 value = strtod(tmp, NULL);
4368 pfree(tmp);
4369 if (errno != 0 || _element == junk) {
4370 pfree(_element);
4371 rtpg_colormap_arg_destroy(arg);
4372 PG_FREE_IF_COPY(pgraster, 0);
4373 elog(ERROR, "RASTER_colorMap: Could not process percent string to value");
4374 PG_RETURN_NULL();
4375 }
4376
4377 /* check percentage */
4378 if (value < 0.) {
4379 elog(NOTICE, "Percentage values cannot be less than zero. Defaulting to zero");
4380 value = 0.;
4381 }
4382 else if (value > 100.) {
4383 elog(NOTICE, "Percentage values cannot be greater than 100. Defaulting to 100");
4384 value = 100.;
4385 }
4386
4387 /* get the true pixel value */
4388 /* TODO: should the percentage be quantile based? */
4389 arg->colormap->entry[arg->colormap->nentry].value = ((value / 100.) * (arg->bandstats->max - arg->bandstats->min)) + arg->bandstats->min;
4390 }
4391 /* straight value */
4392 else {
4393 errno = 0;
4394 arg->colormap->entry[arg->colormap->nentry].value = strtod(_element, &junk);
4395 if (errno != 0 || _element == junk) {
4396 pfree(_element);
4397 rtpg_colormap_arg_destroy(arg);
4398 PG_FREE_IF_COPY(pgraster, 0);
4399 elog(ERROR, "RASTER_colorMap: Could not process string to value");
4400 PG_RETURN_NULL();
4401 }
4402 }
4403
4404 }
4405 /* RGB values (0 - 255) */
4406 else {
4407 int value = 0;
4408
4409 errno = 0;
4410 value = (int) strtod(_element, &junk);
4411 if (errno != 0 || _element == junk) {
4412 pfree(_element);
4413 rtpg_colormap_arg_destroy(arg);
4414 PG_FREE_IF_COPY(pgraster, 0);
4415 elog(ERROR, "RASTER_colorMap: Could not process string to value");
4416 PG_RETURN_NULL();
4417 }
4418
4419 if (value > 255) {
4420 elog(NOTICE, "RGBA value cannot be greater than 255. Defaulting to 255");
4421 value = 255;
4422 }
4423 else if (value < 0) {
4424 elog(NOTICE, "RGBA value cannot be less than zero. Defaulting to zero");
4425 value = 0;
4426 }
4427 arg->colormap->entry[arg->colormap->nentry].color[j - 1] = value;
4428 }
4429
4430 pfree(_element);
4431 }
4432
4433 POSTGIS_RT_DEBUGF(4, "colormap->entry[%d] (isnodata, value, R, G, B, A) = (%d, %f, %d, %d, %d, %d)",
4434 arg->colormap->nentry,
4435 arg->colormap->entry[arg->colormap->nentry].isnodata,
4436 arg->colormap->entry[arg->colormap->nentry].value,
4437 arg->colormap->entry[arg->colormap->nentry].color[0],
4438 arg->colormap->entry[arg->colormap->nentry].color[1],
4439 arg->colormap->entry[arg->colormap->nentry].color[2],
4440 arg->colormap->entry[arg->colormap->nentry].color[3]
4441 );
4442
4443 arg->colormap->nentry++;
4444 }
4445
4446 POSTGIS_RT_DEBUGF(4, "colormap->nentry = %d", arg->colormap->nentry);
4447 POSTGIS_RT_DEBUGF(4, "colormap->ncolor = %d", arg->colormap->ncolor);
4448 }
4449
4450 /* call colormap */
4451 raster = rt_raster_colormap(arg->raster, arg->nband - 1, arg->colormap);
4452 if (raster == NULL) {
4453 rtpg_colormap_arg_destroy(arg);
4454 PG_FREE_IF_COPY(pgraster, 0);
4455 elog(ERROR, "RASTER_colorMap: Could not create new raster with applied colormap");
4456 PG_RETURN_NULL();
4457 }
4458
4459 rtpg_colormap_arg_destroy(arg);
4460 PG_FREE_IF_COPY(pgraster, 0);
4461 pgraster = rt_raster_serialize(raster);
4462 rt_raster_destroy(raster);
4463
4464 POSTGIS_RT_DEBUG(3, "RASTER_colorMap: Done");
4465
4466 if (pgraster == NULL)
4467 PG_RETURN_NULL();
4468
4469 SET_VARSIZE(pgraster, ((rt_pgraster*) pgraster)->size);
4470 PG_RETURN_POINTER(pgraster);
4471 }
4472
4473 PG_FUNCTION_INFO_V1(RASTER_mapAlgebraExpr);
4474 Datum RASTER_mapAlgebraExpr(PG_FUNCTION_ARGS)
4475 {
4476 rt_pgraster *pgraster = NULL;
4477 rt_pgraster *pgrtn = NULL;
4478 rt_raster raster = NULL;
4479 rt_raster newrast = NULL;
4480 rt_band band = NULL;
4481 rt_band newband = NULL;
4482 int x, y, nband, width, height;
4483 double r;
4484 double newnodatavalue = 0.0;
4485 double newinitialvalue = 0.0;
4486 double newval = 0.0;
4487 char *newexpr = NULL;
4488 char *initexpr = NULL;
4489 char *expression = NULL;
4490 int hasnodataval = 0;
4491 double nodataval = 0.;
4492 rt_pixtype newpixeltype;
4493 int skipcomputation = 0;
4494 int len = 0;
4495 const int argkwcount = 3;
4496 enum KEYWORDS { kVAL=0, kX=1, kY=2 };
4497 char *argkw[] = {"[rast]", "[rast.x]", "[rast.y]"};
4498 Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
4499 int argcount = 0;
4500 Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
4501 uint8_t argpos[3] = {0};
4502 char place[12];
4503 int idx = 0;
4504 int ret = -1;
4505 TupleDesc tupdesc;
4506 SPIPlanPtr spi_plan = NULL;
4507 SPITupleTable * tuptable = NULL;
4508 HeapTuple tuple;
4509 char * strFromText = NULL;
4510 Datum *values = NULL;
4511 Datum datum = (Datum)NULL;
4512 char *nulls = NULL;
4513 bool isnull = FALSE;
4514 int i = 0;
4515 int j = 0;
4516
4517 POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraExpr: Starting...");
4518
4519 /* Check raster */
4520 if (PG_ARGISNULL(0)) {
4521 elog(NOTICE, "Raster is NULL. Returning NULL");
4522 PG_RETURN_NULL();
4523 }
4524
4525
4526 /* Deserialize raster */
4527 pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4528 raster = rt_raster_deserialize(pgraster, FALSE);
4529 if (NULL == raster) {
4530 PG_FREE_IF_COPY(pgraster, 0);
4531 elog(ERROR, "RASTER_mapAlgebraExpr: Could not deserialize raster");
4532 PG_RETURN_NULL();
4533 }
4534
4535 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting arguments...");
4536
4537 if (PG_ARGISNULL(1))
4538 nband = 1;
4539 else
4540 nband = PG_GETARG_INT32(1);
4541
4542 if (nband < 1)
4543 nband = 1;
4544
4545
4546 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new empty raster...");
4547
4548 /**
4549 * Create a new empty raster with having the same georeference as the
4550 * provided raster
4551 **/
4552 width = rt_raster_get_width(raster);
4553 height = rt_raster_get_height(raster);
4554
4555 newrast = rt_raster_new(width, height);
4556
4557 if ( NULL == newrast ) {
4558 PG_FREE_IF_COPY(pgraster, 0);
4559 elog(ERROR, "RASTER_mapAlgebraExpr: Could not create a new raster");
4560 PG_RETURN_NULL();
4561 }
4562
4563 rt_raster_set_scale(newrast,
4564 rt_raster_get_x_scale(raster),
4565 rt_raster_get_y_scale(raster));
4566
4567 rt_raster_set_offsets(newrast,
4568 rt_raster_get_x_offset(raster),
4569 rt_raster_get_y_offset(raster));
4570
4571 rt_raster_set_skews(newrast,
4572 rt_raster_get_x_skew(raster),
4573 rt_raster_get_y_skew(raster));
4574
4575 rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
4576
4577
4578 /**
4579 * If this new raster is empty (width = 0 OR height = 0) then there is
4580 * nothing to compute and we return it right now
4581 **/
4582 if (rt_raster_is_empty(newrast))
4583 {
4584 elog(NOTICE, "Raster is empty. Returning an empty raster");
4585 rt_raster_destroy(raster);
4586 PG_FREE_IF_COPY(pgraster, 0);
4587
4588 pgrtn = rt_raster_serialize(newrast);
4589 rt_raster_destroy(newrast);
4590 if (NULL == pgrtn) {
4591
4592 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4593 PG_RETURN_NULL();
4594 }
4595
4596 SET_VARSIZE(pgrtn, pgrtn->size);
4597 PG_RETURN_POINTER(pgrtn);
4598 }
4599
4600
4601 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Getting raster band %d...", nband);
4602
4603 /**
4604 * Check if the raster has the required band. Otherwise, return a raster
4605 * without band
4606 **/
4607 if (!rt_raster_has_band(raster, nband - 1)) {
4608 elog(NOTICE, "Raster does not have the required band. Returning a raster "
4609 "without a band");
4610 rt_raster_destroy(raster);
4611 PG_FREE_IF_COPY(pgraster, 0);
4612
4613 pgrtn = rt_raster_serialize(newrast);
4614 rt_raster_destroy(newrast);
4615 if (NULL == pgrtn) {
4616 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4617 PG_RETURN_NULL();
4618 }
4619
4620 SET_VARSIZE(pgrtn, pgrtn->size);
4621 PG_RETURN_POINTER(pgrtn);
4622 }
4623
4624 /* Get the raster band */
4625 band = rt_raster_get_band(raster, nband - 1);
4626 if ( NULL == band ) {
4627 elog(NOTICE, "Could not get the required band. Returning a raster "
4628 "without a band");
4629 rt_raster_destroy(raster);
4630 PG_FREE_IF_COPY(pgraster, 0);
4631
4632 pgrtn = rt_raster_serialize(newrast);
4633 rt_raster_destroy(newrast);
4634 if (NULL == pgrtn) {
4635 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4636 PG_RETURN_NULL();
4637 }
4638
4639 SET_VARSIZE(pgrtn, pgrtn->size);
4640 PG_RETURN_POINTER(pgrtn);
4641 }
4642
4643 /*
4644 * Get NODATA value
4645 */
4646 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting NODATA value for band...");
4647
4648 if (rt_band_get_hasnodata_flag(band)) {
4649 rt_band_get_nodata(band, &newnodatavalue);
4650 }
4651
4652 else {
4653 newnodatavalue = rt_band_get_min_value(band);
4654 }
4655
4656 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: NODATA value for band: = %f",
4657 newnodatavalue);
4658
4659 /**
4660 * We set the initial value of the future band to nodata value. If nodata
4661 * value is null, then the raster will be initialized to
4662 * rt_band_get_min_value but all the values should be recomputed anyway
4663 **/
4664 newinitialvalue = newnodatavalue;
4665
4666 /**
4667 * Set the new pixeltype
4668 **/
4669 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Setting pixeltype...");
4670
4671 if (PG_ARGISNULL(2)) {
4672 newpixeltype = rt_band_get_pixtype(band);
4673 }
4674
4675 else {
4676 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
4677 newpixeltype = rt_pixtype_index_from_name(strFromText);
4678 pfree(strFromText);
4679 if (newpixeltype == PT_END)
4680 newpixeltype = rt_band_get_pixtype(band);
4681 }
4682
4683 if (newpixeltype == PT_END) {
4684 PG_FREE_IF_COPY(pgraster, 0);
4685 elog(ERROR, "RASTER_mapAlgebraExpr: Invalid pixeltype");
4686 PG_RETURN_NULL();
4687 }
4688
4689 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Pixeltype set to %s",
4690 rt_pixtype_name(newpixeltype));
4691
4692
4693 /* Construct expression for raster values */
4694 if (!PG_ARGISNULL(3)) {
4695 expression = text_to_cstring(PG_GETARG_TEXT_P(3));
4696 len = strlen("SELECT (") + strlen(expression) + strlen(")::double precision");
4697 initexpr = (char *)palloc(len + 1);
4698
4699 memcpy(initexpr, "SELECT (", strlen("SELECT ("));
4700 memcpy(initexpr + strlen("SELECT ("), expression, strlen(expression));
4701 memcpy(initexpr + strlen("SELECT (") + strlen(expression), ")::double precision", strlen(")::double precision"));
4702 initexpr[len] = '\0';
4703
4704 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression is %s", initexpr);
4705
4706 /* We don't need this memory */
4707 /*
4708 pfree(expression);
4709 expression = NULL;
4710 */
4711 }
4712
4713
4714
4715 /**
4716 * Optimization: If a nodataval is provided, use it for newinitialvalue.
4717 * Then, we can initialize the raster with this value and skip the
4718 * computation of nodata values one by one in the main computing loop
4719 **/
4720 if (!PG_ARGISNULL(4)) {
4721 hasnodataval = 1;
4722 nodataval = PG_GETARG_FLOAT8(4);
4723 newinitialvalue = nodataval;
4724
4725 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new initial value = %f",
4726 newinitialvalue);
4727 }
4728 else
4729 hasnodataval = 0;
4730
4731
4732
4733 /**
4734 * Optimization: If the raster is only filled with nodata values return
4735 * right now a raster filled with the newinitialvalue
4736 * TODO: Call rt_band_check_isnodata instead?
4737 **/
4738 if (rt_band_get_isnodata_flag(band)) {
4739
4740 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Band is a nodata band, returning "
4741 "a raster filled with nodata");
4742
4743 ret = rt_raster_generate_new_band(newrast, newpixeltype,
4744 newinitialvalue, TRUE, newnodatavalue, 0);
4745
4746 /* Free memory */
4747 if (initexpr)
4748 pfree(initexpr);
4749 rt_raster_destroy(raster);
4750 PG_FREE_IF_COPY(pgraster, 0);
4751
4752 /* Serialize created raster */
4753 pgrtn = rt_raster_serialize(newrast);
4754 rt_raster_destroy(newrast);
4755 if (NULL == pgrtn) {
4756 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4757 PG_RETURN_NULL();
4758 }
4759
4760 SET_VARSIZE(pgrtn, pgrtn->size);
4761 PG_RETURN_POINTER(pgrtn);
4762 }
4763
4764
4765 /**
4766 * Optimization: If expression resume to 'RAST' and hasnodataval is zero,
4767 * we can just return the band from the original raster
4768 **/
4769 if (initexpr != NULL && ( !strcmp(initexpr, "SELECT [rast]") || !strcmp(initexpr, "SELECT [rast.val]") ) && !hasnodataval) {
4770
4771 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression resumes to RAST. "
4772 "Returning raster with band %d from original raster", nband);
4773
4774 POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster has %d bands",
4775 rt_raster_get_num_bands(newrast));
4776
4777 rt_raster_copy_band(newrast, raster, nband - 1, 0);
4778
4779 POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster now has %d bands",
4780 rt_raster_get_num_bands(newrast));
4781
4782 pfree(initexpr);
4783 rt_raster_destroy(raster);
4784 PG_FREE_IF_COPY(pgraster, 0);
4785
4786 /* Serialize created raster */
4787 pgrtn = rt_raster_serialize(newrast);
4788 rt_raster_destroy(newrast);
4789 if (NULL == pgrtn) {
4790 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4791 PG_RETURN_NULL();
4792 }
4793
4794 SET_VARSIZE(pgrtn, pgrtn->size);
4795 PG_RETURN_POINTER(pgrtn);
4796 }
4797
4798 /**
4799 * Optimization: If expression resume to a constant (it does not contain
4800 * [rast)
4801 **/
4802 if (initexpr != NULL && strstr(initexpr, "[rast") == NULL) {
4803 ret = SPI_connect();
4804 if (ret != SPI_OK_CONNECT) {
4805 PG_FREE_IF_COPY(pgraster, 0);
4806 elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4807 PG_RETURN_NULL();
4808 };
4809
4810 /* Execute the expresion into newval */
4811 ret = SPI_execute(initexpr, FALSE, 0);
4812
4813 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
4814
4815 /* Free memory allocated out of the current context */
4816 if (SPI_tuptable)
4817 SPI_freetuptable(tuptable);
4818 PG_FREE_IF_COPY(pgraster, 0);
4819
4820 SPI_finish();
4821 elog(ERROR, "RASTER_mapAlgebraExpr: Invalid construction for expression");
4822 PG_RETURN_NULL();
4823 }
4824
4825 tupdesc = SPI_tuptable->tupdesc;
4826 tuptable = SPI_tuptable;
4827
4828 tuple = tuptable->vals[0];
4829 newexpr = SPI_getvalue(tuple, tupdesc, 1);
4830 if ( ! newexpr ) {
4831 POSTGIS_RT_DEBUG(3, "Constant expression evaluated to NULL, keeping initvalue");
4832 newval = newinitialvalue;
4833 } else {
4834 newval = atof(newexpr);
4835 }
4836
4837 SPI_freetuptable(tuptable);
4838
4839 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: New raster value = %f",
4840 newval);
4841
4842 SPI_finish();
4843
4844 skipcomputation = 1;
4845
4846 /**
4847 * Compute the new value, set it and we will return after creating the
4848 * new raster
4849 **/
4850 if (!hasnodataval) {
4851 newinitialvalue = newval;
4852 skipcomputation = 2;
4853 }
4854
4855 /* Return the new raster as it will be before computing pixel by pixel */
4856 else if (FLT_NEQ(newval, newinitialvalue)) {
4857 skipcomputation = 2;
4858 }
4859 }
4860
4861 /**
4862 * Create the raster receiving all the computed values. Initialize it to the
4863 * new initial value
4864 **/
4865 ret = rt_raster_generate_new_band(newrast, newpixeltype,
4866 newinitialvalue, TRUE, newnodatavalue, 0);
4867
4868 /**
4869 * Optimization: If expression is NULL, or all the pixels could be set in
4870 * one step, return the initialized raster now
4871 **/
4872 /*if (initexpr == NULL || skipcomputation == 2) {*/
4873 if (expression == NULL || skipcomputation == 2) {
4874
4875 /* Free memory */
4876 if (initexpr)
4877 pfree(initexpr);
4878 rt_raster_destroy(raster);
4879 PG_FREE_IF_COPY(pgraster, 0);
4880
4881 /* Serialize created raster */
4882 pgrtn = rt_raster_serialize(newrast);
4883 rt_raster_destroy(newrast);
4884 if (NULL == pgrtn) {
4885 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4886 PG_RETURN_NULL();
4887 }
4888
4889 SET_VARSIZE(pgrtn, pgrtn->size);
4890 PG_RETURN_POINTER(pgrtn);
4891 }
4892
4893 RASTER_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new raster band...");
4894
4895 /* Get the new raster band */
4896 newband = rt_raster_get_band(newrast, 0);
4897 if ( NULL == newband ) {
4898 elog(NOTICE, "Could not modify band for new raster. Returning new "
4899 "raster with the original band");
4900
4901 if (initexpr)
4902 pfree(initexpr);
4903 rt_raster_destroy(raster);
4904 PG_FREE_IF_COPY(pgraster, 0);
4905
4906 /* Serialize created raster */
4907 pgrtn = rt_raster_serialize(newrast);
4908 rt_raster_destroy(newrast);
4909 if (NULL == pgrtn) {
4910 elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4911 PG_RETURN_NULL();
4912 }
4913
4914 SET_VARSIZE(pgrtn, pgrtn->size);
4915 PG_RETURN_POINTER(pgrtn);
4916 }
4917
4918 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Main computing loop (%d x %d)",
4919 width, height);
4920
4921 if (initexpr != NULL) {
4922 /* Convert [rast.val] to [rast] */
4923 newexpr = rtpg_strreplace(initexpr, "[rast.val]", "[rast]", NULL);
4924 pfree(initexpr); initexpr=newexpr;
4925
4926 sprintf(place,"$1");
4927 for (i = 0, j = 1; i < argkwcount; i++) {
4928 len = 0;
4929 newexpr = rtpg_strreplace(initexpr, argkw[i], place, &len);
4930 pfree(initexpr); initexpr=newexpr;
4931 if (len > 0) {
4932 argtype[argcount] = argkwtypes[i];
4933 argcount++;
4934 argpos[i] = j++;
4935
4936 sprintf(place, "$%d", j);
4937 }
4938 else {
4939 argpos[i] = 0;
4940 }
4941 }
4942
4943 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: initexpr = %s", initexpr);
4944
4945 /* define values */
4946 values = (Datum *) palloc(sizeof(Datum) * argcount);
4947 if (values == NULL) {
4948
4949 SPI_finish();
4950
4951 rt_raster_destroy(raster);
4952 PG_FREE_IF_COPY(pgraster, 0);
4953 rt_raster_destroy(newrast);
4954
4955 elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
4956 PG_RETURN_NULL();
4957 }
4958
4959 /* define nulls */
4960 nulls = (char *)palloc(argcount);
4961 if (nulls == NULL) {
4962
4963 SPI_finish();
4964
4965 rt_raster_destroy(raster);
4966 PG_FREE_IF_COPY(pgraster, 0);
4967 rt_raster_destroy(newrast);
4968
4969 elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
4970 PG_RETURN_NULL();
4971 }
4972
4973 /* Connect to SPI and prepare the expression */
4974 ret = SPI_connect();
4975 if (ret != SPI_OK_CONNECT) {
4976
4977 if (initexpr)
4978 pfree(initexpr);
4979 rt_raster_destroy(raster);
4980 PG_FREE_IF_COPY(pgraster, 0);
4981 rt_raster_destroy(newrast);
4982
4983 elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4984 PG_RETURN_NULL();
4985 };
4986
4987 /* Type of all arguments is FLOAT8OID */
4988 spi_plan = SPI_prepare(initexpr, argcount, argtype);
4989
4990 if (spi_plan == NULL) {
4991
4992 rt_raster_destroy(raster);
4993 PG_FREE_IF_COPY(pgraster, 0);
4994 rt_raster_destroy(newrast);
4995
4996 SPI_finish();
4997
4998 pfree(initexpr);
4999
5000 elog(ERROR, "RASTER_mapAlgebraExpr: Could not prepare expression");
5001 PG_RETURN_NULL();
5002 }
5003 }
5004
5005 for (x = 0; x < width; x++) {
5006 for(y = 0; y < height; y++) {
5007 ret = rt_band_get_pixel(band, x, y, &r, NULL);
5008
5009 /**
5010 * We compute a value only for the withdata value pixel since the
5011 * nodata value has already been set by the first optimization
5012 **/
5013 if (ret == ES_NONE && FLT_NEQ(r, newnodatavalue)) {
5014 if (skipcomputation == 0) {
5015 if (initexpr != NULL) {
5016 /* Reset the null arg flags. */
5017 memset(nulls, 'n', argcount);
5018
5019 for (i = 0; i < argkwcount; i++) {
5020 idx = argpos[i];
5021 if (idx < 1) continue;
5022 idx--;
5023
5024 if (i == kX ) {
5025 /* x is 0 based index, but SQL expects 1 based index */
5026 values[idx] = Int32GetDatum(x+1);
5027 nulls[idx] = ' ';
5028 }
5029 else if (i == kY) {
5030 /* y is 0 based index, but SQL expects 1 based index */
5031 values[idx] = Int32GetDatum(y+1);
5032 nulls[idx] = ' ';
5033 }
5034 else if (i == kVAL ) {
5035 values[idx] = Float8GetDatum(r);
5036 nulls[idx] = ' ';
5037 }
5038
5039 }
5040
5041 ret = SPI_execute_plan(spi_plan, values, nulls, FALSE, 0);
5042 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
5043 SPI_processed != 1) {
5044 if (SPI_tuptable)
5045 SPI_freetuptable(tuptable);
5046
5047 SPI_freeplan(spi_plan);
5048 SPI_finish();
5049
5050 pfree(values);
5051 pfree(nulls);
5052 pfree(initexpr);
5053
5054 rt_raster_destroy(raster);
5055 PG_FREE_IF_COPY(pgraster, 0);
5056 rt_raster_destroy(newrast);
5057
5058 elog(ERROR, "RASTER_mapAlgebraExpr: Error executing prepared plan");
5059
5060 PG_RETURN_NULL();
5061 }
5062
5063 tupdesc = SPI_tuptable->tupdesc;
5064 tuptable = SPI_tuptable;
5065
5066 tuple = tuptable->vals[0];
5067 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
5068 if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
5069 POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) errored, skip setting", x+1,y+1,r);
5070 newval = newinitialvalue;
5071 }
5072 else if ( isnull ) {
5073 POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting", x+1,y+1,r);
5074 newval = newinitialvalue;
5075 } else {
5076 newval = DatumGetFloat8(datum);
5077 }
5078
5079 SPI_freetuptable(tuptable);
5080 }
5081
5082 else
5083 newval = newinitialvalue;
5084
5085 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new value = %f",
5086 newval);
5087 }
5088
5089
5090 rt_band_set_pixel(newband, x, y, newval, NULL);
5091 }
5092
5093 }
5094 }
5095
5096 if (initexpr != NULL) {
5097 SPI_freeplan(spi_plan);
5098 SPI_finish();
5099
5100 pfree(values);
5101 pfree(nulls);
5102 pfree(initexpr);
5103 }
5104 else {
5105 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: no SPI cleanup");
5106 }
5107
5108
5109 /* The newrast band has been modified */
5110
5111 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster modified, serializing it.");
5112 /* Serialize created raster */
5113
5114 rt_raster_destroy(raster);
5115 PG_FREE_IF_COPY(pgraster, 0);
5116
5117 pgrtn = rt_raster_serialize(newrast);
5118 rt_raster_destroy(newrast);
5119 if (NULL == pgrtn)
5120 PG_RETURN_NULL();
5121
5122 SET_VARSIZE(pgrtn, pgrtn->size);
5123
5124 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster serialized");
5125
5126
5127 POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraExpr: returning raster");
5128
5129
5130 PG_RETURN_POINTER(pgrtn);
5131 }
5132
5133 /*
5134 * One raster user function MapAlgebra.
5135 */
5136 PG_FUNCTION_INFO_V1(RASTER_mapAlgebraFct);
5137 Datum RASTER_mapAlgebraFct(PG_FUNCTION_ARGS)
5138 {
5139 rt_pgraster *pgraster = NULL;
5140 rt_pgraster *pgrtn = NULL;
5141 rt_raster raster = NULL;
5142 rt_raster newrast = NULL;
5143 rt_band band = NULL;
5144 rt_band newband = NULL;
5145 int x, y, nband, width, height;
5146 double r;
5147 double newnodatavalue = 0.0;
5148 double newinitialvalue = 0.0;
5149 double newval = 0.0;
5150 rt_pixtype newpixeltype;
5151 int ret = -1;
5152 Oid oid;
5153 FmgrInfo cbinfo;
5154 #if POSTGIS_PGSQL_VERSION < 120
5155 FunctionCallInfoData cbdata;
5156 #else
5157 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS); /* Could be optimized */
5158 #endif
5159 Datum tmpnewval;
5160 char * strFromText = NULL;
5161 int k = 0;
5162
5163 POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFct: STARTING...");
5164
5165 /* Check raster */
5166 if (PG_ARGISNULL(0)) {
5167 elog(WARNING, "Raster is NULL. Returning NULL");
5168 PG_RETURN_NULL();
5169 }
5170
5171
5172 /* Deserialize raster */
5173 pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5174 raster = rt_raster_deserialize(pgraster, FALSE);
5175 if (NULL == raster) {
5176 PG_FREE_IF_COPY(pgraster, 0);
5177 elog(ERROR, "RASTER_mapAlgebraFct: Could not deserialize raster");
5178 PG_RETURN_NULL();
5179 }
5180
5181 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting arguments...");
5182
5183 /* Get the rest of the arguments */
5184
5185 if (PG_ARGISNULL(1))
5186 nband = 1;
5187 else
5188 nband = PG_GETARG_INT32(1);
5189
5190 if (nband < 1)
5191 nband = 1;
5192
5193 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Creating new empty raster...");
5194
5195 /**
5196 * Create a new empty raster with having the same georeference as the
5197 * provided raster
5198 **/
5199 width = rt_raster_get_width(raster);
5200 height = rt_raster_get_height(raster);
5201
5202 newrast = rt_raster_new(width, height);
5203
5204 if ( NULL == newrast ) {
5205
5206 rt_raster_destroy(raster);
5207 PG_FREE_IF_COPY(pgraster, 0);
5208
5209 elog(ERROR, "RASTER_mapAlgebraFct: Could not create a new raster");
5210 PG_RETURN_NULL();
5211 }
5212
5213 rt_raster_set_scale(newrast,
5214 rt_raster_get_x_scale(raster),
5215 rt_raster_get_y_scale(raster));
5216
5217 rt_raster_set_offsets(newrast,
5218 rt_raster_get_x_offset(raster),
5219 rt_raster_get_y_offset(raster));
5220
5221 rt_raster_set_skews(newrast,
5222 rt_raster_get_x_skew(raster),
5223 rt_raster_get_y_skew(raster));
5224
5225 rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
5226
5227
5228 /**
5229 * If this new raster is empty (width = 0 OR height = 0) then there is
5230 * nothing to compute and we return it right now
5231 **/
5232 if (rt_raster_is_empty(newrast))
5233 {
5234 elog(NOTICE, "Raster is empty. Returning an empty raster");
5235 rt_raster_destroy(raster);
5236 PG_FREE_IF_COPY(pgraster, 0);
5237
5238 pgrtn = rt_raster_serialize(newrast);
5239 rt_raster_destroy(newrast);
5240 if (NULL == pgrtn) {
5241 elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5242 PG_RETURN_NULL();
5243 }
5244
5245 SET_VARSIZE(pgrtn, pgrtn->size);
5246 PG_RETURN_POINTER(pgrtn);
5247 }
5248
5249 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Getting raster band %d...", nband);
5250
5251 /**
5252 * Check if the raster has the required band. Otherwise, return a raster
5253 * without band
5254 **/
5255 if (!rt_raster_has_band(raster, nband - 1)) {
5256 elog(NOTICE, "Raster does not have the required band. Returning a raster "
5257 "without a band");
5258 rt_raster_destroy(raster);
5259 PG_FREE_IF_COPY(pgraster, 0);
5260
5261 pgrtn = rt_raster_serialize(newrast);
5262 rt_raster_destroy(newrast);
5263 if (NULL == pgrtn) {
5264 elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5265 PG_RETURN_NULL();
5266 }
5267
5268 SET_VARSIZE(pgrtn, pgrtn->size);
5269 PG_RETURN_POINTER(pgrtn);
5270 }
5271
5272 /* Get the raster band */
5273 band = rt_raster_get_band(raster, nband - 1);
5274 if ( NULL == band ) {
5275 elog(NOTICE, "Could not get the required band. Returning a raster "
5276 "without a band");
5277 rt_raster_destroy(raster);
5278 PG_FREE_IF_COPY(pgraster, 0);
5279
5280 pgrtn = rt_raster_serialize(newrast);
5281 rt_raster_destroy(newrast);
5282 if (NULL == pgrtn) {
5283 elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5284 PG_RETURN_NULL();
5285 }
5286
5287 SET_VARSIZE(pgrtn, pgrtn->size);
5288 PG_RETURN_POINTER(pgrtn);
5289 }
5290
5291 /*
5292 * Get NODATA value
5293 */
5294 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting NODATA value for band...");
5295
5296 if (rt_band_get_hasnodata_flag(band)) {
5297 rt_band_get_nodata(band, &newnodatavalue);
5298 }
5299
5300 else {
5301 newnodatavalue = rt_band_get_min_value(band);
5302 }
5303
5304 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: NODATA value for band: %f",
5305 newnodatavalue);
5306 /**
5307 * We set the initial value of the future band to nodata value. If nodata
5308 * value is null, then the raster will be initialized to
5309 * rt_band_get_min_value but all the values should be recomputed anyway
5310 **/
5311 newinitialvalue = newnodatavalue;
5312
5313 /**
5314 * Set the new pixeltype
5315 **/
5316 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Setting pixeltype...");
5317
5318 if (PG_ARGISNULL(2)) {
5319 newpixeltype = rt_band_get_pixtype(band);
5320 }
5321
5322 else {
5323 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5324 newpixeltype = rt_pixtype_index_from_name(strFromText);
5325 pfree(strFromText);
5326 if (newpixeltype == PT_END)
5327 newpixeltype = rt_band_get_pixtype(band);
5328 }
5329
5330 if (newpixeltype == PT_END) {
5331
5332 rt_raster_destroy(raster);
5333 PG_FREE_IF_COPY(pgraster, 0);
5334 rt_raster_destroy(newrast);
5335
5336 elog(ERROR, "RASTER_mapAlgebraFct: Invalid pixeltype");
5337 PG_RETURN_NULL();
5338 }
5339
5340 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Pixeltype set to %s",
5341 rt_pixtype_name(newpixeltype));
5342
5343 /* Get the name of the callback user function for raster values */
5344 if (PG_ARGISNULL(3)) {
5345
5346 rt_raster_destroy(raster);
5347 PG_FREE_IF_COPY(pgraster, 0);
5348 rt_raster_destroy(newrast);
5349
5350 elog(ERROR, "RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
5351 PG_RETURN_NULL();
5352 }
5353
5354 oid = PG_GETARG_OID(3);
5355 if (oid == InvalidOid) {
5356
5357 rt_raster_destroy(raster);
5358 PG_FREE_IF_COPY(pgraster, 0);
5359 rt_raster_destroy(newrast);
5360
5361 elog(ERROR, "RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
5362 PG_RETURN_NULL();
5363 }
5364
5365 fmgr_info(oid, &cbinfo);
5366
5367 /* function cannot return set */
5368 if (cbinfo.fn_retset) {
5369
5370 rt_raster_destroy(raster);
5371 PG_FREE_IF_COPY(pgraster, 0);
5372 rt_raster_destroy(newrast);
5373
5374 elog(ERROR, "RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
5375 PG_RETURN_NULL();
5376 }
5377 /* function should have correct # of args */
5378 else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {
5379
5380 rt_raster_destroy(raster);
5381 PG_FREE_IF_COPY(pgraster, 0);
5382 rt_raster_destroy(newrast);
5383
5384 elog(ERROR, "RASTER_mapAlgebraFct: Function does not have two or three input parameters");
5385 PG_RETURN_NULL();
5386 }
5387
5388 if (cbinfo.fn_nargs == 2)
5389 k = 1;
5390 else
5391 k = 2;
5392
5393 if (func_volatile(oid) == 'v') {
5394 elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5395 }
5396
5397 /* prep function call data */
5398 #if POSTGIS_PGSQL_VERSION < 120
5399 InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5400
5401 cbdata.argnull[0] = FALSE;
5402 cbdata.argnull[1] = FALSE;
5403 cbdata.argnull[2] = FALSE;
5404 #else
5405 InitFunctionCallInfoData(*cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5406
5407 cbdata->args[0].isnull = FALSE;
5408 cbdata->args[1].isnull = FALSE;
5409 cbdata->args[2].isnull = FALSE;
5410 #endif
5411
5412 /* check that the function isn't strict if the args are null. */
5413 if (PG_ARGISNULL(4)) {
5414 if (cbinfo.fn_strict) {
5415
5416 rt_raster_destroy(raster);
5417 PG_FREE_IF_COPY(pgraster, 0);
5418 rt_raster_destroy(newrast);
5419
5420 elog(ERROR, "RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
5421 PG_RETURN_NULL();
5422 }
5423
5424 #if POSTGIS_PGSQL_VERSION < 120
5425 cbdata.arg[k] = (Datum)NULL;
5426 cbdata.argnull[k] = TRUE;
5427 #else
5428 cbdata->args[k].value = (Datum)NULL;
5429 cbdata->args[k].isnull = TRUE;
5430 #endif
5431 }
5432 else {
5433 #if POSTGIS_PGSQL_VERSION < 120
5434 cbdata.arg[k] = PG_GETARG_DATUM(4);
5435 #else
5436 cbdata->args[k].value = PG_GETARG_DATUM(4);
5437 #endif
5438 }
5439
5440 /**
5441 * Optimization: If the raster is only filled with nodata values return
5442 * right now a raster filled with the nodatavalueexpr
5443 * TODO: Call rt_band_check_isnodata instead?
5444 **/
5445 if (rt_band_get_isnodata_flag(band)) {
5446
5447 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Band is a nodata band, returning "
5448 "a raster filled with nodata");
5449
5450 ret = rt_raster_generate_new_band(newrast, newpixeltype,
5451 newinitialvalue, TRUE, newnodatavalue, 0);
5452
5453 rt_raster_destroy(raster);
5454 PG_FREE_IF_COPY(pgraster, 0);
5455
5456 /* Serialize created raster */
5457 pgrtn = rt_raster_serialize(newrast);
5458 rt_raster_destroy(newrast);
5459 if (NULL == pgrtn) {
5460 elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5461 PG_RETURN_NULL();
5462 }
5463
5464 SET_VARSIZE(pgrtn, pgrtn->size);
5465 PG_RETURN_POINTER(pgrtn);
5466 }
5467
5468
5469 /**
5470 * Create the raster receiving all the computed values. Initialize it to the
5471 * new initial value
5472 **/
5473 ret = rt_raster_generate_new_band(newrast, newpixeltype,
5474 newinitialvalue, TRUE, newnodatavalue, 0);
5475
5476 /* Get the new raster band */
5477 newband = rt_raster_get_band(newrast, 0);
5478 if ( NULL == newband ) {
5479 elog(NOTICE, "Could not modify band for new raster. Returning new "
5480 "raster with the original band");
5481
5482 rt_raster_destroy(raster);
5483 PG_FREE_IF_COPY(pgraster, 0);
5484
5485 /* Serialize created raster */
5486 pgrtn = rt_raster_serialize(newrast);
5487 rt_raster_destroy(newrast);
5488 if (NULL == pgrtn) {
5489 elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5490 PG_RETURN_NULL();
5491 }
5492
5493 SET_VARSIZE(pgrtn, pgrtn->size);
5494 PG_RETURN_POINTER(pgrtn);
5495 }
5496
5497 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Main computing loop (%d x %d)",
5498 width, height);
5499
5500 for (x = 0; x < width; x++) {
5501 for(y = 0; y < height; y++) {
5502 ret = rt_band_get_pixel(band, x, y, &r, NULL);
5503
5504 /**
5505 * We compute a value only for the withdata value pixel since the
5506 * nodata value has already been set by the first optimization
5507 **/
5508 if (ret == ES_NONE) {
5509 if (FLT_EQ(r, newnodatavalue)) {
5510 if (cbinfo.fn_strict) {
5511 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
5512 continue;
5513 }
5514 #if POSTGIS_PGSQL_VERSION < 120
5515 cbdata.argnull[0] = TRUE;
5516 cbdata.arg[0] = (Datum)NULL;
5517 #else
5518 cbdata->args[0].isnull = TRUE;
5519 cbdata->args[0].value = (Datum)NULL;
5520 #endif
5521 }
5522 else {
5523 #if POSTGIS_PGSQL_VERSION < 120
5524 cbdata.argnull[0] = FALSE;
5525 cbdata.arg[0] = Float8GetDatum(r);
5526 #else
5527 cbdata->args[0].isnull = FALSE;
5528 cbdata->args[0].value = Float8GetDatum(r);
5529 #endif
5530 }
5531
5532 /* Add pixel positions if callback has proper # of args */
5533 if (cbinfo.fn_nargs == 3) {
5534 Datum d[2];
5535 ArrayType *a;
5536
5537 d[0] = Int32GetDatum(x+1);
5538 d[1] = Int32GetDatum(y+1);
5539
5540 a = construct_array(d, 2, INT4OID, sizeof(int32), true, 'i');
5541
5542 #if POSTGIS_PGSQL_VERSION < 120
5543 cbdata.argnull[1] = FALSE;
5544 cbdata.arg[1] = PointerGetDatum(a);
5545 #else
5546 cbdata->args[1].isnull = FALSE;
5547 cbdata->args[1].value = PointerGetDatum(a);
5548 #endif
5549 }
5550
5551 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: (%dx%d), r = %f",
5552 x, y, r);
5553
5554 #if POSTGIS_PGSQL_VERSION < 120
5555 tmpnewval = FunctionCallInvoke(&cbdata);
5556
5557 if (cbdata.isnull) {
5558 newval = newnodatavalue;
5559 }
5560 #else
5561 tmpnewval = FunctionCallInvoke(cbdata);
5562
5563 if (cbdata->isnull)
5564 {
5565 newval = newnodatavalue;
5566 }
5567 #endif
5568 else {
5569 newval = DatumGetFloat8(tmpnewval);
5570 }
5571
5572 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: new value = %f",
5573 newval);
5574
5575 rt_band_set_pixel(newband, x, y, newval, NULL);
5576 }
5577
5578 }
5579 }
5580
5581 /* The newrast band has been modified */
5582
5583 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster modified, serializing it.");
5584 /* Serialize created raster */
5585
5586 rt_raster_destroy(raster);
5587 PG_FREE_IF_COPY(pgraster, 0);
5588
5589 pgrtn = rt_raster_serialize(newrast);
5590 rt_raster_destroy(newrast);
5591 if (NULL == pgrtn)
5592 PG_RETURN_NULL();
5593
5594 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster serialized");
5595
5596 POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFct: returning raster");
5597
5598 SET_VARSIZE(pgrtn, pgrtn->size);
5599 PG_RETURN_POINTER(pgrtn);
5600 }
5601
5602 /**
5603 * One raster neighborhood MapAlgebra
5604 */
5605 PG_FUNCTION_INFO_V1(RASTER_mapAlgebraFctNgb);
5606 Datum RASTER_mapAlgebraFctNgb(PG_FUNCTION_ARGS)
5607 {
5608 rt_pgraster *pgraster = NULL;
5609 rt_pgraster *pgrtn = NULL;
5610 rt_raster raster = NULL;
5611 rt_raster newrast = NULL;
5612 rt_band band = NULL;
5613 rt_band newband = NULL;
5614 int x, y, nband, width, height, ngbwidth, ngbheight, winwidth, winheight, u, v, nIndex, nNullItems;
5615 double r, rpix;
5616 double newnodatavalue = 0.0;
5617 double newinitialvalue = 0.0;
5618 double newval = 0.0;
5619 rt_pixtype newpixeltype;
5620 int ret = -1;
5621 Oid oid;
5622 FmgrInfo cbinfo;
5623 #if POSTGIS_PGSQL_VERSION < 120
5624 FunctionCallInfoData cbdata;
5625 #else
5626 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS); /* Could be optimized */
5627 #endif
5628 Datum tmpnewval;
5629 ArrayType * neighborDatum;
5630 char * strFromText = NULL;
5631 text * txtNodataMode = NULL;
5632 text * txtCallbackParam = NULL;
5633 int intReplace = 0;
5634 float fltReplace = 0;
5635 bool valuereplace = false, pixelreplace, nNodataOnly = true, nNullSkip = false;
5636 Datum * neighborData = NULL;
5637 bool * neighborNulls = NULL;
5638 int neighborDims[2];
5639 int neighborLbs[2];
5640 int16 typlen;
5641 bool typbyval;
5642 char typalign;
5643
5644 POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFctNgb: STARTING...");
5645
5646 /* Check raster */
5647 if (PG_ARGISNULL(0)) {
5648 elog(WARNING, "Raster is NULL. Returning NULL");
5649 PG_RETURN_NULL();
5650 }
5651
5652
5653 /* Deserialize raster */
5654 pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5655 raster = rt_raster_deserialize(pgraster, FALSE);
5656 if (NULL == raster)
5657 {
5658 PG_FREE_IF_COPY(pgraster, 0);
5659 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not deserialize raster");
5660 PG_RETURN_NULL();
5661 }
5662
5663 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Getting arguments...");
5664
5665 /* Get the rest of the arguments */
5666
5667 if (PG_ARGISNULL(1))
5668 nband = 1;
5669 else
5670 nband = PG_GETARG_INT32(1);
5671
5672 if (nband < 1)
5673 nband = 1;
5674
5675 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Creating new empty raster...");
5676
5677 /**
5678 * Create a new empty raster with having the same georeference as the
5679 * provided raster
5680 **/
5681 width = rt_raster_get_width(raster);
5682 height = rt_raster_get_height(raster);
5683
5684 newrast = rt_raster_new(width, height);
5685
5686 if ( NULL == newrast ) {
5687 rt_raster_destroy(raster);
5688 PG_FREE_IF_COPY(pgraster, 0);
5689 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not create a new raster");
5690 PG_RETURN_NULL();
5691 }
5692
5693 rt_raster_set_scale(newrast,
5694 rt_raster_get_x_scale(raster),
5695 rt_raster_get_y_scale(raster));
5696
5697 rt_raster_set_offsets(newrast,
5698 rt_raster_get_x_offset(raster),
5699 rt_raster_get_y_offset(raster));
5700
5701 rt_raster_set_skews(newrast,
5702 rt_raster_get_x_skew(raster),
5703 rt_raster_get_y_skew(raster));
5704
5705 rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
5706
5707
5708 /**
5709 * If this new raster is empty (width = 0 OR height = 0) then there is
5710 * nothing to compute and we return it right now
5711 **/
5712 if (rt_raster_is_empty(newrast))
5713 {
5714 elog(NOTICE, "Raster is empty. Returning an empty raster");
5715 rt_raster_destroy(raster);
5716 PG_FREE_IF_COPY(pgraster, 0);
5717
5718 pgrtn = rt_raster_serialize(newrast);
5719 rt_raster_destroy(newrast);
5720 if (NULL == pgrtn) {
5721 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5722 PG_RETURN_NULL();
5723 }
5724
5725 SET_VARSIZE(pgrtn, pgrtn->size);
5726 PG_RETURN_POINTER(pgrtn);
5727 }
5728
5729 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Getting raster band %d...", nband);
5730
5731 /**
5732 * Check if the raster has the required band. Otherwise, return a raster
5733 * without band
5734 **/
5735 if (!rt_raster_has_band(raster, nband - 1)) {
5736 elog(NOTICE, "Raster does not have the required band. Returning a raster "
5737 "without a band");
5738 rt_raster_destroy(raster);
5739 PG_FREE_IF_COPY(pgraster, 0);
5740
5741 pgrtn = rt_raster_serialize(newrast);
5742 rt_raster_destroy(newrast);
5743 if (NULL == pgrtn) {
5744 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5745 PG_RETURN_NULL();
5746 }
5747
5748 SET_VARSIZE(pgrtn, pgrtn->size);
5749 PG_RETURN_POINTER(pgrtn);
5750 }
5751
5752 /* Get the raster band */
5753 band = rt_raster_get_band(raster, nband - 1);
5754 if ( NULL == band ) {
5755 elog(NOTICE, "Could not get the required band. Returning a raster "
5756 "without a band");
5757 rt_raster_destroy(raster);
5758 PG_FREE_IF_COPY(pgraster, 0);
5759
5760 pgrtn = rt_raster_serialize(newrast);
5761 rt_raster_destroy(newrast);
5762 if (NULL == pgrtn) {
5763 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5764 PG_RETURN_NULL();
5765 }
5766
5767 SET_VARSIZE(pgrtn, pgrtn->size);
5768 PG_RETURN_POINTER(pgrtn);
5769 }
5770
5771 /*
5772 * Get NODATA value
5773 */
5774 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Getting NODATA value for band...");
5775
5776 if (rt_band_get_hasnodata_flag(band)) {
5777 rt_band_get_nodata(band, &newnodatavalue);
5778 }
5779
5780 else {
5781 newnodatavalue = rt_band_get_min_value(band);
5782 }
5783
5784 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: NODATA value for band: %f",
5785 newnodatavalue);
5786 /**
5787 * We set the initial value of the future band to nodata value. If nodata
5788 * value is null, then the raster will be initialized to
5789 * rt_band_get_min_value but all the values should be recomputed anyway
5790 **/
5791 newinitialvalue = newnodatavalue;
5792
5793 /**
5794 * Set the new pixeltype
5795 **/
5796 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Setting pixeltype...");
5797
5798 if (PG_ARGISNULL(2)) {
5799 newpixeltype = rt_band_get_pixtype(band);
5800 }
5801
5802 else {
5803 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5804 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Pixeltype parameter: %s", strFromText);
5805 newpixeltype = rt_pixtype_index_from_name(strFromText);
5806 pfree(strFromText);
5807 if (newpixeltype == PT_END)
5808 newpixeltype = rt_band_get_pixtype(band);
5809 }
5810
5811 if (newpixeltype == PT_END) {
5812
5813 rt_raster_destroy(raster);
5814 PG_FREE_IF_COPY(pgraster, 0);
5815 rt_raster_destroy(newrast);
5816
5817 elog(ERROR, "RASTER_mapAlgebraFctNgb: Invalid pixeltype");
5818 PG_RETURN_NULL();
5819 }
5820
5821 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Pixeltype set to %s (%d)",
5822 rt_pixtype_name(newpixeltype), newpixeltype);
5823
5824 /* Get the name of the callback userfunction */
5825 if (PG_ARGISNULL(5)) {
5826
5827 rt_raster_destroy(raster);
5828 PG_FREE_IF_COPY(pgraster, 0);
5829 rt_raster_destroy(newrast);
5830
5831 elog(ERROR, "RASTER_mapAlgebraFctNgb: Required function is missing");
5832 PG_RETURN_NULL();
5833 }
5834
5835 oid = PG_GETARG_OID(5);
5836 if (oid == InvalidOid) {
5837
5838 rt_raster_destroy(raster);
5839 PG_FREE_IF_COPY(pgraster, 0);
5840 rt_raster_destroy(newrast);
5841
5842 elog(ERROR, "RASTER_mapAlgebraFctNgb: Got invalid function object id");
5843 PG_RETURN_NULL();
5844 }
5845
5846 fmgr_info(oid, &cbinfo);
5847
5848 /* function cannot return set */
5849 if (cbinfo.fn_retset) {
5850
5851 rt_raster_destroy(raster);
5852 PG_FREE_IF_COPY(pgraster, 0);
5853 rt_raster_destroy(newrast);
5854
5855 elog(ERROR, "RASTER_mapAlgebraFctNgb: Function provided must return double precision not resultset");
5856 PG_RETURN_NULL();
5857 }
5858 /* function should have correct # of args */
5859 else if (cbinfo.fn_nargs != 3) {
5860
5861 rt_raster_destroy(raster);
5862 PG_FREE_IF_COPY(pgraster, 0);
5863 rt_raster_destroy(newrast);
5864
5865 elog(ERROR, "RASTER_mapAlgebraFctNgb: Function does not have three input parameters");
5866 PG_RETURN_NULL();
5867 }
5868
5869 if (func_volatile(oid) == 'v') {
5870 elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5871 }
5872
5873 /* prep function call data */
5874 #if POSTGIS_PGSQL_VERSION < 120
5875 InitFunctionCallInfoData(cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5876 memset(cbdata.argnull, FALSE, sizeof(bool) * 3);
5877 #else
5878 InitFunctionCallInfoData(*cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5879 cbdata->args[0].isnull = FALSE;
5880 cbdata->args[1].isnull = FALSE;
5881 cbdata->args[2].isnull = FALSE;
5882 #endif
5883
5884 /* check that the function isn't strict if the args are null. */
5885 if (PG_ARGISNULL(7)) {
5886 if (cbinfo.fn_strict) {
5887
5888 rt_raster_destroy(raster);
5889 PG_FREE_IF_COPY(pgraster, 0);
5890 rt_raster_destroy(newrast);
5891
5892 elog(ERROR, "RASTER_mapAlgebraFctNgb: Strict callback functions cannot have NULL parameters");
5893 PG_RETURN_NULL();
5894 }
5895
5896 #if POSTGIS_PGSQL_VERSION < 120
5897 cbdata.arg[2] = (Datum)NULL;
5898 cbdata.argnull[2] = TRUE;
5899 #else
5900 cbdata->args[2].value = (Datum)NULL;
5901 cbdata->args[2].isnull = TRUE;
5902 #endif
5903 }
5904 else {
5905 #if POSTGIS_PGSQL_VERSION < 120
5906 cbdata.arg[2] = PG_GETARG_DATUM(7);
5907 #else
5908 cbdata->args[2].value = PG_GETARG_DATUM(7);
5909 #endif
5910 }
5911
5912 /**
5913 * Optimization: If the raster is only filled with nodata values return
5914 * right now a raster filled with the nodatavalueexpr
5915 * TODO: Call rt_band_check_isnodata instead?
5916 **/
5917 if (rt_band_get_isnodata_flag(band)) {
5918
5919 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Band is a nodata band, returning "
5920 "a raster filled with nodata");
5921
5922 rt_raster_generate_new_band(newrast, newpixeltype,
5923 newinitialvalue, TRUE, newnodatavalue, 0);
5924
5925 rt_raster_destroy(raster);
5926 PG_FREE_IF_COPY(pgraster, 0);
5927
5928 /* Serialize created raster */
5929 pgrtn = rt_raster_serialize(newrast);
5930 rt_raster_destroy(newrast);
5931 if (NULL == pgrtn) {
5932 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5933 PG_RETURN_NULL();
5934 }
5935
5936 SET_VARSIZE(pgrtn, pgrtn->size);
5937 PG_RETURN_POINTER(pgrtn);
5938 }
5939
5940
5941 /**
5942 * Create the raster receiving all the computed values. Initialize it to the
5943 * new initial value
5944 **/
5945 rt_raster_generate_new_band(newrast, newpixeltype,
5946 newinitialvalue, TRUE, newnodatavalue, 0);
5947
5948 /* Get the new raster band */
5949 newband = rt_raster_get_band(newrast, 0);
5950 if ( NULL == newband ) {
5951 elog(NOTICE, "Could not modify band for new raster. Returning new "
5952 "raster with the original band");
5953
5954 rt_raster_destroy(raster);
5955 PG_FREE_IF_COPY(pgraster, 0);
5956
5957 /* Serialize created raster */
5958 pgrtn = rt_raster_serialize(newrast);
5959 rt_raster_destroy(newrast);
5960 if (NULL == pgrtn) {
5961 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5962 PG_RETURN_NULL();
5963 }
5964
5965 SET_VARSIZE(pgrtn, pgrtn->size);
5966 PG_RETURN_POINTER(pgrtn);
5967 }
5968
5969 /* Get the width of the neighborhood */
5970 if (PG_ARGISNULL(3) || PG_GETARG_INT32(3) <= 0) {
5971 elog(NOTICE, "Neighborhood width is NULL or <= 0. Returning new "
5972 "raster with the original band");
5973
5974 rt_raster_destroy(raster);
5975 PG_FREE_IF_COPY(pgraster, 0);
5976
5977 /* Serialize created raster */
5978 pgrtn = rt_raster_serialize(newrast);
5979 rt_raster_destroy(newrast);
5980 if (NULL == pgrtn) {
5981 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5982 PG_RETURN_NULL();
5983 }
5984
5985 SET_VARSIZE(pgrtn, pgrtn->size);
5986 PG_RETURN_POINTER(pgrtn);
5987 }
5988
5989 ngbwidth = PG_GETARG_INT32(3);
5990 winwidth = ngbwidth * 2 + 1;
5991
5992 /* Get the height of the neighborhood */
5993 if (PG_ARGISNULL(4) || PG_GETARG_INT32(4) <= 0) {
5994 elog(NOTICE, "Neighborhood height is NULL or <= 0. Returning new "
5995 "raster with the original band");
5996
5997 rt_raster_destroy(raster);
5998 PG_FREE_IF_COPY(pgraster, 0);
5999
6000 /* Serialize created raster */
6001 pgrtn = rt_raster_serialize(newrast);
6002 rt_raster_destroy(newrast);
6003 if (NULL == pgrtn) {
6004 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
6005 PG_RETURN_NULL();
6006 }
6007
6008 SET_VARSIZE(pgrtn, pgrtn->size);
6009 PG_RETURN_POINTER(pgrtn);
6010 }
6011
6012 ngbheight = PG_GETARG_INT32(4);
6013 winheight = ngbheight * 2 + 1;
6014
6015 /* Get the type of NODATA behavior for the neighborhoods. */
6016 if (PG_ARGISNULL(6)) {
6017 elog(NOTICE, "Neighborhood NODATA behavior defaulting to 'ignore'");
6018 txtNodataMode = cstring_to_text("ignore");
6019 }
6020 else {
6021 txtNodataMode = PG_GETARG_TEXT_P(6);
6022 }
6023
6024 txtCallbackParam = (text*)palloc(VARSIZE(txtNodataMode));
6025 SET_VARSIZE(txtCallbackParam, VARSIZE(txtNodataMode));
6026 memcpy((void *)VARDATA(txtCallbackParam), (void *)VARDATA(txtNodataMode), VARSIZE(txtNodataMode) - VARHDRSZ);
6027
6028 /* pass the nodata mode into the user function */
6029 #if POSTGIS_PGSQL_VERSION < 120
6030 cbdata.arg[1] = CStringGetDatum(txtCallbackParam);
6031 #else
6032 cbdata->args[1].value = CStringGetDatum(txtCallbackParam);
6033 #endif
6034
6035 strFromText = text_to_cstring(txtNodataMode);
6036 strFromText = rtpg_strtoupper(strFromText);
6037
6038 if (strcmp(strFromText, "VALUE") == 0)
6039 valuereplace = true;
6040 else if (strcmp(strFromText, "IGNORE") != 0 && strcmp(strFromText, "NULL") != 0) {
6041 /* if the text is not "IGNORE" or "NULL", it may be a numerical value */
6042 if (sscanf(strFromText, "%d", &intReplace) <= 0 && sscanf(strFromText, "%f", &fltReplace) <= 0) {
6043 /* the value is NOT an integer NOR a floating point */
6044 elog(NOTICE, "Neighborhood NODATA mode is not recognized. Must be one of 'value', 'ignore', "
6045 "'NULL', or a numeric value. Returning new raster with the original band");
6046
6047 /* clean up the nodatamode string */
6048 pfree(txtCallbackParam);
6049 pfree(strFromText);
6050
6051 rt_raster_destroy(raster);
6052 PG_FREE_IF_COPY(pgraster, 0);
6053
6054 /* Serialize created raster */
6055 pgrtn = rt_raster_serialize(newrast);
6056 rt_raster_destroy(newrast);
6057 if (NULL == pgrtn) {
6058 elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
6059 PG_RETURN_NULL();
6060 }
6061
6062 SET_VARSIZE(pgrtn, pgrtn->size);
6063 PG_RETURN_POINTER(pgrtn);
6064 }
6065 }
6066 else if (strcmp(strFromText, "NULL") == 0) {
6067 /* this setting means that the neighborhood should be skipped if any of the values are null */
6068 nNullSkip = true;
6069 }
6070
6071 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Main computing loop (%d x %d)",
6072 width, height);
6073
6074 /* Allocate room for the neighborhood. */
6075 neighborData = (Datum *)palloc(winwidth * winheight * sizeof(Datum));
6076 neighborNulls = (bool *)palloc(winwidth * winheight * sizeof(bool));
6077
6078 /* The dimensions of the neighborhood array, for creating a multi-dimensional array. */
6079 neighborDims[0] = winwidth;
6080 neighborDims[1] = winheight;
6081
6082 /* The lower bounds for the new multi-dimensional array. */
6083 neighborLbs[0] = 1;
6084 neighborLbs[1] = 1;
6085
6086 /* Get information about the type of item in the multi-dimensional array (float8). */
6087 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
6088
6089 for (x = 0 + ngbwidth; x < width - ngbwidth; x++) {
6090 for(y = 0 + ngbheight; y < height - ngbheight; y++) {
6091 /* populate an array with the pixel values in the neighborhood */
6092 nIndex = 0;
6093 nNullItems = 0;
6094 nNodataOnly = true;
6095 pixelreplace = false;
6096 if (valuereplace) {
6097 ret = rt_band_get_pixel(band, x, y, &rpix, NULL);
6098 if (ret == ES_NONE && FLT_NEQ(rpix, newnodatavalue)) {
6099 pixelreplace = true;
6100 }
6101 }
6102 for (u = x - ngbwidth; u <= x + ngbwidth; u++) {
6103 for (v = y - ngbheight; v <= y + ngbheight; v++) {
6104 ret = rt_band_get_pixel(band, u, v, &r, NULL);
6105 if (ret == ES_NONE) {
6106 if (FLT_NEQ(r, newnodatavalue)) {
6107 /* If the pixel value for this neighbor cell is not NODATA */
6108 neighborData[nIndex] = Float8GetDatum((double)r);
6109 neighborNulls[nIndex] = false;
6110 nNodataOnly = false;
6111 }
6112 else {
6113 /* If the pixel value for this neighbor cell is NODATA */
6114 if (valuereplace && pixelreplace) {
6115 /* Replace the NODATA value with the currently processing pixel. */
6116 neighborData[nIndex] = Float8GetDatum((double)rpix);
6117 neighborNulls[nIndex] = false;
6118 /* do not increment nNullItems, since the user requested that the */
6119 /* neighborhood replace NODATA values with the central pixel value */
6120 }
6121 else {
6122 neighborData[nIndex] = PointerGetDatum(NULL);
6123 neighborNulls[nIndex] = true;
6124 nNullItems++;
6125 }
6126 }
6127 }
6128 else {
6129 /* Fill this will NULL if we can't read the raster pixel. */
6130 neighborData[nIndex] = PointerGetDatum(NULL);
6131 neighborNulls[nIndex] = true;
6132 nNullItems++;
6133 }
6134 /* Next neighbor position */
6135 nIndex++;
6136 }
6137 }
6138
6139 /**
6140 * We compute a value only for the withdata value neighborhood since the
6141 * nodata value has already been set by the first optimization
6142 **/
6143 if (!(nNodataOnly || /* neighborhood only contains NODATA -- OR -- */
6144 (nNullSkip && nNullItems > 0) || /* neighborhood should skip any NODATA cells, and a NODATA cell was detected -- OR -- */
6145 (valuereplace && nNullItems > 0))) { /* neighborhood should replace NODATA cells with the central pixel value, and a NODATA cell was detected */
6146 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: (%dx%d), %dx%d neighborhood",
6147 x, y, winwidth, winheight);
6148
6149 neighborDatum = construct_md_array((void *)neighborData, neighborNulls, 2, neighborDims, neighborLbs,
6150 FLOAT8OID, typlen, typbyval, typalign);
6151
6152 #if POSTGIS_PGSQL_VERSION < 120
6153 /* Assign the neighbor matrix as the first argument to the user function */
6154 cbdata.arg[0] = PointerGetDatum(neighborDatum);
6155
6156 /* Invoke the user function */
6157 tmpnewval = FunctionCallInvoke(&cbdata);
6158
6159 /* Get the return value of the user function */
6160 if (cbdata.isnull) {
6161 newval = newnodatavalue;
6162 }
6163 #else
6164 /* Assign the neighbor matrix as the first argument to the user function */
6165 cbdata->args[0].value = PointerGetDatum(neighborDatum);
6166
6167 /* Invoke the user function */
6168 tmpnewval = FunctionCallInvoke(cbdata);
6169
6170 /* Get the return value of the user function */
6171 if (cbdata->isnull)
6172 {
6173 newval = newnodatavalue;
6174 }
6175 #endif
6176 else {
6177 newval = DatumGetFloat8(tmpnewval);
6178 }
6179
6180 POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: new value = %f",
6181 newval);
6182
6183 rt_band_set_pixel(newband, x, y, newval, NULL);
6184 }
6185
6186 /* reset the number of null items in the neighborhood */
6187 nNullItems = 0;
6188 }
6189 }
6190
6191
6192 /* clean up */
6193 pfree(neighborNulls);
6194 pfree(neighborData);
6195 pfree(strFromText);
6196 pfree(txtCallbackParam);
6197
6198 rt_raster_destroy(raster);
6199 PG_FREE_IF_COPY(pgraster, 0);
6200
6201 /* The newrast band has been modified */
6202
6203 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: raster modified, serializing it.");
6204 /* Serialize created raster */
6205
6206 pgrtn = rt_raster_serialize(newrast);
6207 rt_raster_destroy(newrast);
6208 if (NULL == pgrtn)
6209 PG_RETURN_NULL();
6210
6211 POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: raster serialized");
6212 POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFctNgb: returning raster");
6213
6214 SET_VARSIZE(pgrtn, pgrtn->size);
6215 PG_RETURN_POINTER(pgrtn);
6216 }
6217
6218 #define ARGKWCOUNT 8
6219
6220 /**
6221 * Two raster MapAlgebra
6222 */
6223 PG_FUNCTION_INFO_V1(RASTER_mapAlgebra2);
6224 Datum RASTER_mapAlgebra2(PG_FUNCTION_ARGS)
6225 {
6226 const uint32_t set_count = 2;
6227 rt_pgraster *pgrast[2] = { NULL, NULL };
6228 int pgrastpos[2] = {-1, -1};
6229 rt_pgraster *pgrtn;
6230 rt_raster rast[2] = {NULL};
6231 int _isempty[2] = {0};
6232 uint32_t bandindex[2] = {0};
6233 rt_raster _rast[2] = {NULL};
6234 rt_band _band[2] = {NULL};
6235 int _hasnodata[2] = {0};
6236 double _nodataval[2] = {0};
6237 double _offset[4] = {0.};
6238 double _rastoffset[2][4] = {{0.}};
6239 int _haspixel[2] = {0};
6240 double _pixel[2] = {0};
6241 int _pos[2][2] = {{0}};
6242 uint16_t _dim[2][2] = {{0}};
6243
6244 char *pixtypename = NULL;
6245 rt_pixtype pixtype = PT_END;
6246 char *extenttypename = NULL;
6247 rt_extenttype extenttype = ET_INTERSECTION;
6248
6249 rt_raster raster = NULL;
6250 rt_band band = NULL;
6251 uint16_t dim[2] = {0};
6252 int haspixel = 0;
6253 double pixel = 0.;
6254 double nodataval = 0;
6255 double gt[6] = {0.};
6256
6257 Oid calltype = InvalidOid;
6258
6259 const uint32_t spi_count = 3;
6260 uint16_t spi_exprpos[3] = {4, 7, 8};
6261 uint32_t spi_argcount[3] = {0};
6262 char *expr = NULL;
6263 char *sql = NULL;
6264 SPIPlanPtr spi_plan[3] = {NULL};
6265 uint16_t spi_empty = 0;
6266 Oid *argtype = NULL;
6267 uint8_t argpos[3][8] = {{0}};
6268 char *argkw[] = {"[rast1.x]", "[rast1.y]", "[rast1.val]", "[rast1]", "[rast2.x]", "[rast2.y]", "[rast2.val]", "[rast2]"};
6269 Datum values[ARGKWCOUNT];
6270 char nulls[ARGKWCOUNT];
6271 TupleDesc tupdesc;
6272 SPITupleTable *tuptable = NULL;
6273 HeapTuple tuple;
6274 Datum datum;
6275 bool isnull = FALSE;
6276 int hasargval[3] = {0};
6277 double argval[3] = {0.};
6278 int hasnodatanodataval = 0;
6279 double nodatanodataval = 0;
6280 int isnodata = 0;
6281
6282 Oid ufc_noid = InvalidOid;
6283 FmgrInfo ufl_info;
6284 #if POSTGIS_PGSQL_VERSION < 120
6285 FunctionCallInfoData ufc_info;
6286 #else
6287 LOCAL_FCINFO(ufc_info, FUNC_MAX_ARGS); /* Could be optimized */
6288 #endif
6289 int ufc_nullcount = 0;
6290
6291 int idx = 0;
6292 uint32_t i = 0;
6293 uint32_t j = 0;
6294 uint32_t k = 0;
6295 uint32_t x = 0;
6296 uint32_t y = 0;
6297 int _x = 0;
6298 int _y = 0;
6299 int err;
6300 int aligned = 0;
6301 int len = 0;
6302
6303 POSTGIS_RT_DEBUG(3, "Starting RASTER_mapAlgebra2");
6304
6305 for (i = 0, j = 0; i < set_count; i++) {
6306 if (!PG_ARGISNULL(j)) {
6307 pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
6308 pgrastpos[i] = j;
6309 j++;
6310
6311 /* raster */
6312 rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
6313 if (!rast[i]) {
6314 for (k = 0; k <= i; k++) {
6315 if (k < i && rast[k] != NULL)
6316 rt_raster_destroy(rast[k]);
6317 if (pgrastpos[k] != -1)
6318 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6319 }
6320 elog(ERROR, "RASTER_mapAlgebra2: Could not deserialize the %s raster", i < 1 ? "first" : "second");
6321 PG_RETURN_NULL();
6322 }
6323
6324 /* empty */
6325 _isempty[i] = rt_raster_is_empty(rast[i]);
6326
6327 /* band index */
6328 if (!PG_ARGISNULL(j)) {
6329 bandindex[i] = PG_GETARG_INT32(j);
6330 }
6331 j++;
6332 }
6333 else {
6334 _isempty[i] = 1;
6335 j += 2;
6336 }
6337
6338 POSTGIS_RT_DEBUGF(3, "_isempty[%d] = %d", i, _isempty[i]);
6339 }
6340
6341 /* both rasters are NULL */
6342 if (rast[0] == NULL && rast[1] == NULL) {
6343 elog(NOTICE, "The two rasters provided are NULL. Returning NULL");
6344 for (k = 0; k < set_count; k++) {
6345 if (pgrastpos[k] != -1)
6346 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6347 }
6348 PG_RETURN_NULL();
6349 }
6350
6351 /* both rasters are empty */
6352 if (_isempty[0] && _isempty[1]) {
6353 elog(NOTICE, "The two rasters provided are empty. Returning empty raster");
6354
6355 raster = rt_raster_new(0, 0);
6356 if (raster == NULL) {
6357 for (k = 0; k < set_count; k++) {
6358 if (rast[k] != NULL)
6359 rt_raster_destroy(rast[k]);
6360 if (pgrastpos[k] != -1)
6361 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6362 }
6363 elog(ERROR, "RASTER_mapAlgebra2: Could not create empty raster");
6364 PG_RETURN_NULL();
6365 }
6366 rt_raster_set_scale(raster, 0, 0);
6367
6368 pgrtn = rt_raster_serialize(raster);
6369 rt_raster_destroy(raster);
6370 if (!pgrtn)
6371 PG_RETURN_NULL();
6372
6373 SET_VARSIZE(pgrtn, pgrtn->size);
6374 PG_RETURN_POINTER(pgrtn);
6375 }
6376
6377 /* replace the empty or NULL raster with one matching the other */
6378 if (
6379 (rast[0] == NULL || _isempty[0]) ||
6380 (rast[1] == NULL || _isempty[1])
6381 ) {
6382 /* first raster is empty */
6383 if (rast[0] == NULL || _isempty[0]) {
6384 i = 0;
6385 j = 1;
6386 }
6387 /* second raster is empty */
6388 else {
6389 i = 1;
6390 j = 0;
6391 }
6392
6393 _rast[j] = rast[j];
6394
6395 /* raster is empty, destroy it */
6396 if (_rast[i] != NULL)
6397 rt_raster_destroy(_rast[i]);
6398
6399 _dim[i][0] = rt_raster_get_width(_rast[j]);
6400 _dim[i][1] = rt_raster_get_height(_rast[j]);
6401 _dim[j][0] = rt_raster_get_width(_rast[j]);
6402 _dim[j][1] = rt_raster_get_height(_rast[j]);
6403
6404 _rast[i] = rt_raster_new(
6405 _dim[j][0],
6406 _dim[j][1]
6407 );
6408 if (_rast[i] == NULL) {
6409 rt_raster_destroy(_rast[j]);
6410 for (k = 0; k < set_count; k++) {
6411 if (pgrastpos[k] != -1)
6412 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6413 }
6414 elog(ERROR, "RASTER_mapAlgebra2: Could not create NODATA raster");
6415 PG_RETURN_NULL();
6416 }
6417 rt_raster_set_srid(_rast[i], rt_raster_get_srid(_rast[j]));
6418
6419 rt_raster_get_geotransform_matrix(_rast[j], gt);
6420 rt_raster_set_geotransform_matrix(_rast[i], gt);
6421 }
6422 else {
6423 _rast[0] = rast[0];
6424 _dim[0][0] = rt_raster_get_width(_rast[0]);
6425 _dim[0][1] = rt_raster_get_height(_rast[0]);
6426
6427 _rast[1] = rast[1];
6428 _dim[1][0] = rt_raster_get_width(_rast[1]);
6429 _dim[1][1] = rt_raster_get_height(_rast[1]);
6430 }
6431
6432 /* SRID must match */
6433 /*
6434 if (rt_raster_get_srid(_rast[0]) != rt_raster_get_srid(_rast[1])) {
6435 elog(NOTICE, "The two rasters provided have different SRIDs. Returning NULL");
6436 for (k = 0; k < set_count; k++) {
6437 if (_rast[k] != NULL)
6438 rt_raster_destroy(_rast[k]);
6439 if (pgrastpos[k] != -1)
6440 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6441 }
6442 PG_RETURN_NULL();
6443 }
6444 */
6445
6446 /* same alignment */
6447 if (rt_raster_same_alignment(_rast[0], _rast[1], &aligned, NULL) != ES_NONE) {
6448 for (k = 0; k < set_count; k++) {
6449 if (_rast[k] != NULL)
6450 rt_raster_destroy(_rast[k]);
6451 if (pgrastpos[k] != -1)
6452 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6453 }
6454 elog(ERROR, "RASTER_mapAlgebra2: Could not test for alignment on the two rasters");
6455 PG_RETURN_NULL();
6456 }
6457 if (!aligned) {
6458 elog(NOTICE, "The two rasters provided do not have the same alignment. Returning NULL");
6459 for (k = 0; k < set_count; k++) {
6460 if (_rast[k] != NULL)
6461 rt_raster_destroy(_rast[k]);
6462 if (pgrastpos[k] != -1)
6463 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6464 }
6465 PG_RETURN_NULL();
6466 }
6467
6468 /* pixel type */
6469 if (!PG_ARGISNULL(5)) {
6470 pixtypename = text_to_cstring(PG_GETARG_TEXT_P(5));
6471 /* Get the pixel type index */
6472 pixtype = rt_pixtype_index_from_name(pixtypename);
6473 if (pixtype == PT_END ) {
6474 for (k = 0; k < set_count; k++) {
6475 if (_rast[k] != NULL)
6476 rt_raster_destroy(_rast[k]);
6477 if (pgrastpos[k] != -1)
6478 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6479 }
6480 elog(ERROR, "RASTER_mapAlgebra2: Invalid pixel type: %s", pixtypename);
6481 PG_RETURN_NULL();
6482 }
6483 }
6484
6485 /* extent type */
6486 if (!PG_ARGISNULL(6)) {
6487 extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(6))));
6488 extenttype = rt_util_extent_type(extenttypename);
6489 }
6490 POSTGIS_RT_DEBUGF(3, "extenttype: %d %s", extenttype, extenttypename);
6491
6492 /* computed raster from extent type */
6493 err = rt_raster_from_two_rasters(
6494 _rast[0], _rast[1],
6495 extenttype,
6496 &raster, _offset
6497 );
6498 if (err != ES_NONE) {
6499 for (k = 0; k < set_count; k++) {
6500 if (_rast[k] != NULL)
6501 rt_raster_destroy(_rast[k]);
6502 if (pgrastpos[k] != -1)
6503 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6504 }
6505 elog(ERROR, "RASTER_mapAlgebra2: Could not get output raster of correct extent");
6506 PG_RETURN_NULL();
6507 }
6508
6509 /* copy offsets */
6510 _rastoffset[0][0] = _offset[0];
6511 _rastoffset[0][1] = _offset[1];
6512 _rastoffset[1][0] = _offset[2];
6513 _rastoffset[1][1] = _offset[3];
6514
6515 /* get output raster dimensions */
6516 dim[0] = rt_raster_get_width(raster);
6517 dim[1] = rt_raster_get_height(raster);
6518
6519 i = 2;
6520 /* handle special cases for extent */
6521 switch (extenttype) {
6522 case ET_FIRST:
6523 i = 0;
6524 /* fall through */
6525 case ET_SECOND:
6526 if (i > 1)
6527 i = 1;
6528
6529 if (
6530 _isempty[i] && (
6531 (extenttype == ET_FIRST && i == 0) ||
6532 (extenttype == ET_SECOND && i == 1)
6533 )
6534 ) {
6535 elog(NOTICE, "The %s raster is NULL. Returning NULL", (i != 1 ? "FIRST" : "SECOND"));
6536 for (k = 0; k < set_count; k++) {
6537 if (_rast[k] != NULL)
6538 rt_raster_destroy(_rast[k]);
6539 if (pgrastpos[k] != -1)
6540 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6541 }
6542 rt_raster_destroy(raster);
6543 PG_RETURN_NULL();
6544 }
6545
6546 /* specified band not found */
6547 if (!rt_raster_has_band(_rast[i], bandindex[i] - 1)) {
6548 elog(NOTICE, "The %s raster does not have the band at index %d. Returning no band raster of correct extent",
6549 (i != 1 ? "FIRST" : "SECOND"), bandindex[i]
6550 );
6551
6552 for (k = 0; k < set_count; k++) {
6553 if (_rast[k] != NULL)
6554 rt_raster_destroy(_rast[k]);
6555 if (pgrastpos[k] != -1)
6556 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6557 }
6558
6559 pgrtn = rt_raster_serialize(raster);
6560 rt_raster_destroy(raster);
6561 if (!pgrtn) PG_RETURN_NULL();
6562
6563 SET_VARSIZE(pgrtn, pgrtn->size);
6564 PG_RETURN_POINTER(pgrtn);
6565 }
6566 break;
6567 case ET_UNION:
6568 break;
6569 case ET_INTERSECTION:
6570 /* no intersection */
6571 if (
6572 _isempty[0] || _isempty[1] ||
6573 !dim[0] || !dim[1]
6574 ) {
6575 elog(NOTICE, "The two rasters provided have no intersection. Returning no band raster");
6576
6577 /* raster has dimension, replace with no band raster */
6578 if (dim[0] || dim[1]) {
6579 rt_raster_destroy(raster);
6580
6581 raster = rt_raster_new(0, 0);
6582 if (raster == NULL) {
6583 for (k = 0; k < set_count; k++) {
6584 if (_rast[k] != NULL)
6585 rt_raster_destroy(_rast[k]);
6586 if (pgrastpos[k] != -1)
6587 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6588 }
6589 elog(ERROR, "RASTER_mapAlgebra2: Could not create no band raster");
6590 PG_RETURN_NULL();
6591 }
6592
6593 rt_raster_set_scale(raster, 0, 0);
6594 rt_raster_set_srid(raster, rt_raster_get_srid(_rast[0]));
6595 }
6596
6597 for (k = 0; k < set_count; k++) {
6598 if (_rast[k] != NULL)
6599 rt_raster_destroy(_rast[k]);
6600 if (pgrastpos[k] != -1)
6601 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6602 }
6603
6604 pgrtn = rt_raster_serialize(raster);
6605 rt_raster_destroy(raster);
6606 if (!pgrtn) PG_RETURN_NULL();
6607
6608 SET_VARSIZE(pgrtn, pgrtn->size);
6609 PG_RETURN_POINTER(pgrtn);
6610 }
6611 break;
6612 case ET_LAST:
6613 case ET_CUSTOM:
6614 for (k = 0; k < set_count; k++) {
6615 if (_rast[k] != NULL)
6616 rt_raster_destroy(_rast[k]);
6617 if (pgrastpos[k] != -1)
6618 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6619 }
6620 elog(ERROR, "RASTER_mapAlgebra2: ET_LAST and ET_CUSTOM are not implemented");
6621 PG_RETURN_NULL();
6622 break;
6623 }
6624
6625 /* both rasters do not have specified bands */
6626 if (
6627 (!_isempty[0] && !rt_raster_has_band(_rast[0], bandindex[0] - 1)) &&
6628 (!_isempty[1] && !rt_raster_has_band(_rast[1], bandindex[1] - 1))
6629 ) {
6630 elog(NOTICE, "The two rasters provided do not have the respectively specified band indices. Returning no band raster of correct extent");
6631
6632 for (k = 0; k < set_count; k++) {
6633 if (_rast[k] != NULL)
6634 rt_raster_destroy(_rast[k]);
6635 if (pgrastpos[k] != -1)
6636 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6637 }
6638
6639 pgrtn = rt_raster_serialize(raster);
6640 rt_raster_destroy(raster);
6641 if (!pgrtn) PG_RETURN_NULL();
6642
6643 SET_VARSIZE(pgrtn, pgrtn->size);
6644 PG_RETURN_POINTER(pgrtn);
6645 }
6646
6647 /* get bands */
6648 for (i = 0; i < set_count; i++) {
6649 if (_isempty[i] || !rt_raster_has_band(_rast[i], bandindex[i] - 1)) {
6650 _hasnodata[i] = 1;
6651 _nodataval[i] = 0;
6652
6653 continue;
6654 }
6655
6656 _band[i] = rt_raster_get_band(_rast[i], bandindex[i] - 1);
6657 if (_band[i] == NULL) {
6658 for (k = 0; k < set_count; k++) {
6659 if (_rast[k] != NULL)
6660 rt_raster_destroy(_rast[k]);
6661 if (pgrastpos[k] != -1)
6662 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6663 }
6664 rt_raster_destroy(raster);
6665 elog(ERROR, "RASTER_mapAlgebra2: Could not get band %d of the %s raster",
6666 bandindex[i],
6667 (i < 1 ? "FIRST" : "SECOND")
6668 );
6669 PG_RETURN_NULL();
6670 }
6671
6672 _hasnodata[i] = rt_band_get_hasnodata_flag(_band[i]);
6673 if (_hasnodata[i])
6674 rt_band_get_nodata(_band[i], &(_nodataval[i]));
6675 }
6676
6677 /* pixtype is PT_END, get pixtype based upon extent */
6678 if (pixtype == PT_END) {
6679 if ((extenttype == ET_SECOND && !_isempty[1]) || _isempty[0])
6680 pixtype = rt_band_get_pixtype(_band[1]);
6681 else
6682 pixtype = rt_band_get_pixtype(_band[0]);
6683 }
6684
6685 /* nodata value for new band */
6686 if (extenttype == ET_SECOND && !_isempty[1] && _hasnodata[1]) {
6687 nodataval = _nodataval[1];
6688 }
6689 else if (!_isempty[0] && _hasnodata[0]) {
6690 nodataval = _nodataval[0];
6691 }
6692 else if (!_isempty[1] && _hasnodata[1]) {
6693 nodataval = _nodataval[1];
6694 }
6695 else {
6696 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));
6697 nodataval = rt_pixtype_get_min_value(pixtype);
6698 }
6699
6700 /* add band to output raster */
6701 if (rt_raster_generate_new_band(
6702 raster,
6703 pixtype,
6704 nodataval,
6705 1, nodataval,
6706 0
6707 ) < 0) {
6708 for (k = 0; k < set_count; k++) {
6709 if (_rast[k] != NULL)
6710 rt_raster_destroy(_rast[k]);
6711 if (pgrastpos[k] != -1)
6712 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6713 }
6714 rt_raster_destroy(raster);
6715 elog(ERROR, "RASTER_mapAlgebra2: Could not add new band to output raster");
6716 PG_RETURN_NULL();
6717 }
6718
6719 /* get output band */
6720 band = rt_raster_get_band(raster, 0);
6721 if (band == NULL) {
6722 for (k = 0; k < set_count; k++) {
6723 if (_rast[k] != NULL)
6724 rt_raster_destroy(_rast[k]);
6725 if (pgrastpos[k] != -1)
6726 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6727 }
6728 rt_raster_destroy(raster);
6729 elog(ERROR, "RASTER_mapAlgebra2: Could not get newly added band of output raster");
6730 PG_RETURN_NULL();
6731 }
6732
6733 POSTGIS_RT_DEBUGF(4, "offsets = (%d, %d, %d, %d)",
6734 (int) _rastoffset[0][0],
6735 (int) _rastoffset[0][1],
6736 (int) _rastoffset[1][0],
6737 (int) _rastoffset[1][1]
6738 );
6739
6740 POSTGIS_RT_DEBUGF(4, "metadata = (%f, %f, %d, %d, %f, %f, %f, %f, %d)",
6741 rt_raster_get_x_offset(raster),
6742 rt_raster_get_y_offset(raster),
6743 rt_raster_get_width(raster),
6744 rt_raster_get_height(raster),
6745 rt_raster_get_x_scale(raster),
6746 rt_raster_get_y_scale(raster),
6747 rt_raster_get_x_skew(raster),
6748 rt_raster_get_y_skew(raster),
6749 rt_raster_get_srid(raster)
6750 );
6751
6752 /*
6753 determine who called this function
6754 Arg 4 will either be text or regprocedure
6755 */
6756 POSTGIS_RT_DEBUG(3, "checking parameter type for arg 4");
6757 calltype = get_fn_expr_argtype(fcinfo->flinfo, 4);
6758
6759 switch(calltype) {
6760 case TEXTOID: {
6761 POSTGIS_RT_DEBUG(3, "arg 4 is \"expression\"!");
6762
6763 /* connect SPI */
6764 if (SPI_connect() != SPI_OK_CONNECT) {
6765 for (k = 0; k < set_count; k++) {
6766 if (_rast[k] != NULL)
6767 rt_raster_destroy(_rast[k]);
6768 if (pgrastpos[k] != -1)
6769 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6770 }
6771 rt_raster_destroy(raster);
6772 elog(ERROR, "RASTER_mapAlgebra2: Could not connect to the SPI manager");
6773 PG_RETURN_NULL();
6774 }
6775
6776 /* reset hasargval */
6777 memset(hasargval, 0, sizeof(int) * spi_count);
6778
6779 /*
6780 process expressions
6781
6782 spi_exprpos elements are:
6783 4 - expression => spi_plan[0]
6784 7 - nodata1expr => spi_plan[1]
6785 8 - nodata2expr => spi_plan[2]
6786 */
6787 for (i = 0; i < spi_count; i++) {
6788 if (!PG_ARGISNULL(spi_exprpos[i])) {
6789 char *tmp = NULL;
6790 char place[5] = "$1";
6791 expr = text_to_cstring(PG_GETARG_TEXT_P(spi_exprpos[i]));
6792 POSTGIS_RT_DEBUGF(3, "raw expr #%d: %s", i, expr);
6793
6794 for (j = 0, k = 1; j < ARGKWCOUNT; j++) {
6795 /* attempt to replace keyword with placeholder */
6796 len = 0;
6797 tmp = rtpg_strreplace(expr, argkw[j], place, &len);
6798 pfree(expr);
6799 expr = tmp;
6800
6801 if (len) {
6802 spi_argcount[i]++;
6803 argpos[i][j] = k++;
6804
6805 sprintf(place, "$%d", k);
6806 }
6807 else
6808 argpos[i][j] = 0;
6809 }
6810
6811 len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision");
6812 sql = (char *) palloc(len + 1);
6813 if (sql == NULL) {
6814
6815 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6816 SPI_finish();
6817
6818 for (k = 0; k < set_count; k++) {
6819 if (_rast[k] != NULL)
6820 rt_raster_destroy(_rast[k]);
6821 if (pgrastpos[k] != -1)
6822 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6823 }
6824 rt_raster_destroy(raster);
6825
6826 elog(ERROR, "RASTER_mapAlgebra2: Could not allocate memory for expression parameter %d", spi_exprpos[i]);
6827 PG_RETURN_NULL();
6828 }
6829
6830 memcpy(sql, "SELECT (", strlen("SELECT ("));
6831 memcpy(sql + strlen("SELECT ("), expr, strlen(expr));
6832 memcpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision"));
6833 sql[len] = '\0';
6834
6835 POSTGIS_RT_DEBUGF(3, "sql #%d: %s", i, sql);
6836
6837 /* create prepared plan */
6838 if (spi_argcount[i]) {
6839 argtype = (Oid *) palloc(spi_argcount[i] * sizeof(Oid));
6840 if (argtype == NULL) {
6841
6842 pfree(sql);
6843 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6844 SPI_finish();
6845
6846 for (k = 0; k < set_count; k++) {
6847 if (_rast[k] != NULL)
6848 rt_raster_destroy(_rast[k]);
6849 if (pgrastpos[k] != -1)
6850 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6851 }
6852 rt_raster_destroy(raster);
6853
6854 elog(ERROR, "RASTER_mapAlgebra2: Could not allocate memory for prepared plan argtypes of expression parameter %d", spi_exprpos[i]);
6855 PG_RETURN_NULL();
6856 }
6857
6858 /* specify datatypes of parameters */
6859 for (j = 0, k = 0; j < ARGKWCOUNT; j++) {
6860 if (argpos[i][j] < 1) continue;
6861
6862 /* positions are INT4 */
6863 if (
6864 (strstr(argkw[j], "[rast1.x]") != NULL) ||
6865 (strstr(argkw[j], "[rast1.y]") != NULL) ||
6866 (strstr(argkw[j], "[rast2.x]") != NULL) ||
6867 (strstr(argkw[j], "[rast2.y]") != NULL)
6868 ) {
6869 argtype[k] = INT4OID;
6870 }
6871 /* everything else is FLOAT8 */
6872 else {
6873 argtype[k] = FLOAT8OID;
6874 }
6875
6876 k++;
6877 }
6878
6879 spi_plan[i] = SPI_prepare(sql, spi_argcount[i], argtype);
6880 pfree(argtype);
6881
6882 if (spi_plan[i] == NULL) {
6883
6884 pfree(sql);
6885 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6886 SPI_finish();
6887
6888 for (k = 0; k < set_count; k++) {
6889 if (_rast[k] != NULL)
6890 rt_raster_destroy(_rast[k]);
6891 if (pgrastpos[k] != -1)
6892 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6893 }
6894 rt_raster_destroy(raster);
6895
6896 elog(ERROR, "RASTER_mapAlgebra2: Could not create prepared plan of expression parameter %d", spi_exprpos[i]);
6897 PG_RETURN_NULL();
6898 }
6899 }
6900 /* no args, just execute query */
6901 else {
6902 err = SPI_execute(sql, TRUE, 0);
6903 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
6904
6905 pfree(sql);
6906 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6907 SPI_finish();
6908
6909 for (k = 0; k < set_count; k++) {
6910 if (_rast[k] != NULL)
6911 rt_raster_destroy(_rast[k]);
6912 if (pgrastpos[k] != -1)
6913 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6914 }
6915 rt_raster_destroy(raster);
6916
6917 elog(ERROR, "RASTER_mapAlgebra2: Could not evaluate expression parameter %d", spi_exprpos[i]);
6918 PG_RETURN_NULL();
6919 }
6920
6921 /* get output of prepared plan */
6922 tupdesc = SPI_tuptable->tupdesc;
6923 tuptable = SPI_tuptable;
6924 tuple = tuptable->vals[0];
6925
6926 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
6927 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
6928
6929 pfree(sql);
6930 if (SPI_tuptable) SPI_freetuptable(tuptable);
6931 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6932 SPI_finish();
6933
6934 for (k = 0; k < set_count; k++) {
6935 if (_rast[k] != NULL)
6936 rt_raster_destroy(_rast[k]);
6937 if (pgrastpos[k] != -1)
6938 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6939 }
6940 rt_raster_destroy(raster);
6941
6942 elog(ERROR, "RASTER_mapAlgebra2: Could not get result of expression parameter %d", spi_exprpos[i]);
6943 PG_RETURN_NULL();
6944 }
6945
6946 if (!isnull) {
6947 hasargval[i] = 1;
6948 argval[i] = DatumGetFloat8(datum);
6949 }
6950
6951 if (SPI_tuptable) SPI_freetuptable(tuptable);
6952 }
6953
6954 pfree(sql);
6955 }
6956 else
6957 spi_empty++;
6958 }
6959
6960 /* nodatanodataval */
6961 if (!PG_ARGISNULL(9)) {
6962 hasnodatanodataval = 1;
6963 nodatanodataval = PG_GETARG_FLOAT8(9);
6964 }
6965 else
6966 hasnodatanodataval = 0;
6967 break;
6968 }
6969 case REGPROCEDUREOID: {
6970 POSTGIS_RT_DEBUG(3, "arg 4 is \"userfunction\"!");
6971 if (!PG_ARGISNULL(4)) {
6972
6973 ufc_nullcount = 0;
6974 ufc_noid = PG_GETARG_OID(4);
6975
6976 /* get function info */
6977 fmgr_info(ufc_noid, &ufl_info);
6978
6979 /* function cannot return set */
6980 err = 0;
6981 if (ufl_info.fn_retset) {
6982 err = 1;
6983 }
6984 /* function should have correct # of args */
6985 else if (ufl_info.fn_nargs < 3 || ufl_info.fn_nargs > 4) {
6986 err = 2;
6987 }
6988
6989 /*
6990 TODO: consider adding checks of the userfunction parameters
6991 should be able to use get_fn_expr_argtype() of fmgr.c
6992 */
6993
6994 if (err > 0) {
6995 for (k = 0; k < set_count; k++) {
6996 if (_rast[k] != NULL)
6997 rt_raster_destroy(_rast[k]);
6998 if (pgrastpos[k] != -1)
6999 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7000 }
7001 rt_raster_destroy(raster);
7002
7003 if (err > 1)
7004 elog(ERROR, "RASTER_mapAlgebra2: Function provided must have three or four input parameters");
7005 else
7006 elog(ERROR, "RASTER_mapAlgebra2: Function provided must return double precision not resultset");
7007 PG_RETURN_NULL();
7008 }
7009
7010 if (func_volatile(ufc_noid) == 'v') {
7011 elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
7012 }
7013
7014 /* prep function call data */
7015 #if POSTGIS_PGSQL_VERSION < 120
7016 InitFunctionCallInfoData(ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7017 memset(ufc_info.argnull, FALSE, sizeof(bool) * ufl_info.fn_nargs);
7018 #else
7019 InitFunctionCallInfoData(
7020 *ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7021 ufc_info->args[0].isnull = FALSE;
7022 ufc_info->args[1].isnull = FALSE;
7023 ufc_info->args[2].isnull = FALSE;
7024 if (ufl_info.fn_nargs == 4)
7025 ufc_info->args[3].isnull = FALSE;
7026 #endif
7027
7028 if (ufl_info.fn_nargs != 4)
7029 k = 2;
7030 else
7031 k = 3;
7032 #if POSTGIS_PGSQL_VERSION < 120
7033 if (!PG_ARGISNULL(7)) {
7034 ufc_info.arg[k] = PG_GETARG_DATUM(7);
7035 }
7036 else {
7037 ufc_info.arg[k] = (Datum) NULL;
7038 ufc_info.argnull[k] = TRUE;
7039 ufc_nullcount++;
7040 }
7041 #else
7042 if (!PG_ARGISNULL(7))
7043 {
7044 ufc_info->args[k].value = PG_GETARG_DATUM(7);
7045 }
7046 else
7047 {
7048 ufc_info->args[k].value = (Datum)NULL;
7049 ufc_info->args[k].isnull = TRUE;
7050 ufc_nullcount++;
7051 }
7052 #endif
7053 }
7054 break;
7055 }
7056 default:
7057 for (k = 0; k < set_count; k++) {
7058 if (_rast[k] != NULL)
7059 rt_raster_destroy(_rast[k]);
7060 if (pgrastpos[k] != -1)
7061 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7062 }
7063 rt_raster_destroy(raster);
7064 elog(ERROR, "RASTER_mapAlgebra2: Invalid data type for expression or userfunction");
7065 PG_RETURN_NULL();
7066 break;
7067 }
7068
7069 /* loop over pixels */
7070 /* if any expression present, run */
7071 if ((
7072 (calltype == TEXTOID) && (
7073 (spi_empty != spi_count) || hasnodatanodataval
7074 )
7075 ) || (
7076 (calltype == REGPROCEDUREOID) && (ufc_noid != InvalidOid)
7077 )) {
7078 for (x = 0; x < dim[0]; x++) {
7079 for (y = 0; y < dim[1]; y++) {
7080
7081 /* get pixel from each raster */
7082 for (i = 0; i < set_count; i++) {
7083 _haspixel[i] = 0;
7084 _pixel[i] = 0;
7085
7086 /* row/column */
7087 _x = (int)x - (int)_rastoffset[i][0];
7088 _y = (int)y - (int)_rastoffset[i][1];
7089
7090 /* store _x and _y in 1-based */
7091 _pos[i][0] = _x + 1;
7092 _pos[i][1] = _y + 1;
7093
7094 /* get pixel value */
7095 if (_band[i] == NULL) {
7096 if (!_hasnodata[i]) {
7097 _haspixel[i] = 1;
7098 _pixel[i] = _nodataval[i];
7099 }
7100 }
7101 else if (
7102 !_isempty[i] &&
7103 (_x >= 0 && _x < _dim[i][0]) &&
7104 (_y >= 0 && _y < _dim[i][1])
7105 ) {
7106 err = rt_band_get_pixel(_band[i], _x, _y, &(_pixel[i]), &isnodata);
7107 if (err != ES_NONE) {
7108
7109 if (calltype == TEXTOID) {
7110 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7111 SPI_finish();
7112 }
7113
7114 for (k = 0; k < set_count; k++) {
7115 if (_rast[k] != NULL)
7116 rt_raster_destroy(_rast[k]);
7117 if (pgrastpos[k] != -1)
7118 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7119 }
7120 rt_raster_destroy(raster);
7121
7122 elog(ERROR, "RASTER_mapAlgebra2: Could not get pixel of %s raster", (i < 1 ? "FIRST" : "SECOND"));
7123 PG_RETURN_NULL();
7124 }
7125
7126 if (!_hasnodata[i] || !isnodata)
7127 _haspixel[i] = 1;
7128 }
7129
7130 POSTGIS_RT_DEBUGF(5, "pixel r%d(%d, %d) = %d, %f",
7131 i,
7132 _x, _y,
7133 _haspixel[i],
7134 _pixel[i]
7135 );
7136 }
7137
7138 haspixel = 0;
7139
7140 switch (calltype) {
7141 case TEXTOID: {
7142 /* which prepared plan to use? */
7143 /* !pixel0 && !pixel1 */
7144 /* use nodatanodataval */
7145 if (!_haspixel[0] && !_haspixel[1])
7146 i = 3;
7147 /* pixel0 && !pixel1 */
7148 /* run spi_plan[2] (nodata2expr) */
7149 else if (_haspixel[0] && !_haspixel[1])
7150 i = 2;
7151 /* !pixel0 && pixel1 */
7152 /* run spi_plan[1] (nodata1expr) */
7153 else if (!_haspixel[0] && _haspixel[1])
7154 i = 1;
7155 /* pixel0 && pixel1 */
7156 /* run spi_plan[0] (expression) */
7157 else
7158 i = 0;
7159
7160 /* process values */
7161 if (i == 3) {
7162 if (hasnodatanodataval) {
7163 haspixel = 1;
7164 pixel = nodatanodataval;
7165 }
7166 }
7167 /* has an evaluated value */
7168 else if (hasargval[i]) {
7169 haspixel = 1;
7170 pixel = argval[i];
7171 }
7172 /* prepared plan exists */
7173 else if (spi_plan[i] != NULL) {
7174 POSTGIS_RT_DEBUGF(4, "Using prepared plan: %d", i);
7175
7176 /* reset values to (Datum) NULL */
7177 memset(values, (Datum) NULL, sizeof(Datum) * ARGKWCOUNT);
7178 /* reset nulls to FALSE */
7179 memset(nulls, FALSE, sizeof(char) * ARGKWCOUNT);
7180
7181 /* expression has argument(s) */
7182 if (spi_argcount[i]) {
7183 /* set values and nulls */
7184 for (j = 0; j < ARGKWCOUNT; j++) {
7185 idx = argpos[i][j];
7186 if (idx < 1) continue;
7187 idx--; /* 1-based becomes 0-based */
7188
7189 if (strstr(argkw[j], "[rast1.x]") != NULL) {
7190 values[idx] = _pos[0][0];
7191 }
7192 else if (strstr(argkw[j], "[rast1.y]") != NULL) {
7193 values[idx] = _pos[0][1];
7194 }
7195 else if (
7196 (strstr(argkw[j], "[rast1.val]") != NULL) ||
7197 (strstr(argkw[j], "[rast1]") != NULL)
7198 ) {
7199 if (_isempty[0] || !_haspixel[0])
7200 nulls[idx] = TRUE;
7201 else
7202 values[idx] = Float8GetDatum(_pixel[0]);
7203 }
7204 else if (strstr(argkw[j], "[rast2.x]") != NULL) {
7205 values[idx] = _pos[1][0];
7206 }
7207 else if (strstr(argkw[j], "[rast2.y]") != NULL) {
7208 values[idx] = _pos[1][1];
7209 }
7210 else if (
7211 (strstr(argkw[j], "[rast2.val]") != NULL) ||
7212 (strstr(argkw[j], "[rast2]") != NULL)
7213 ) {
7214 if (_isempty[1] || !_haspixel[1])
7215 nulls[idx] = TRUE;
7216 else
7217 values[idx] = Float8GetDatum(_pixel[1]);
7218 }
7219 }
7220 }
7221
7222 /* run prepared plan */
7223 err = SPI_execute_plan(spi_plan[i], values, nulls, TRUE, 1);
7224 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
7225
7226 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7227 SPI_finish();
7228
7229 for (k = 0; k < set_count; k++) {
7230 if (_rast[k] != NULL)
7231 rt_raster_destroy(_rast[k]);
7232 if (pgrastpos[k] != -1)
7233 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7234 }
7235 rt_raster_destroy(raster);
7236
7237 elog(ERROR, "RASTER_mapAlgebra2: Unexpected error when running prepared statement %d", i);
7238 PG_RETURN_NULL();
7239 }
7240
7241 /* get output of prepared plan */
7242 tupdesc = SPI_tuptable->tupdesc;
7243 tuptable = SPI_tuptable;
7244 tuple = tuptable->vals[0];
7245
7246 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
7247 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
7248
7249 if (SPI_tuptable) SPI_freetuptable(tuptable);
7250 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7251 SPI_finish();
7252
7253 for (k = 0; k < set_count; k++) {
7254 if (_rast[k] != NULL)
7255 rt_raster_destroy(_rast[k]);
7256 if (pgrastpos[k] != -1)
7257 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7258 }
7259 rt_raster_destroy(raster);
7260
7261 elog(ERROR, "RASTER_mapAlgebra2: Could not get result of prepared statement %d", i);
7262 PG_RETURN_NULL();
7263 }
7264
7265 if (!isnull) {
7266 haspixel = 1;
7267 pixel = DatumGetFloat8(datum);
7268 }
7269
7270 if (SPI_tuptable) SPI_freetuptable(tuptable);
7271 }
7272 } break;
7273 case REGPROCEDUREOID: {
7274 Datum d[4];
7275 ArrayType *a;
7276
7277 /* build fcnarg */
7278 for (i = 0; i < set_count; i++) {
7279 #if POSTGIS_PGSQL_VERSION < 120
7280 ufc_info.arg[i] = Float8GetDatum(_pixel[i]);
7281 #else
7282 ufc_info->args[i].value = Float8GetDatum(_pixel[i]);
7283 #endif
7284
7285 if (_haspixel[i]) {
7286 #if POSTGIS_PGSQL_VERSION < 120
7287 ufc_info.argnull[i] = FALSE;
7288 #else
7289 ufc_info->args[i].isnull = FALSE;
7290 #endif
7291 ufc_nullcount--;
7292 }
7293 else {
7294 #if POSTGIS_PGSQL_VERSION < 120
7295 ufc_info.argnull[i] = TRUE;
7296 #else
7297 ufc_info->args[i].isnull = TRUE;
7298 #endif
7299 ufc_nullcount++;
7300 }
7301 }
7302
7303 /* function is strict and null parameter is passed */
7304 /* http://archives.postgresql.org/pgsql-general/2011-11/msg00424.php */
7305 if (ufl_info.fn_strict && ufc_nullcount)
7306 break;
7307
7308 /* 4 parameters, add position */
7309 if (ufl_info.fn_nargs == 4) {
7310 /* Datum of 4 element array */
7311 /* array is (x1, y1, x2, y2) */
7312 for (i = 0; i < set_count; i++) {
7313 if (i < 1) {
7314 d[0] = Int32GetDatum(_pos[i][0]);
7315 d[1] = Int32GetDatum(_pos[i][1]);
7316 }
7317 else {
7318 d[2] = Int32GetDatum(_pos[i][0]);
7319 d[3] = Int32GetDatum(_pos[i][1]);
7320 }
7321 }
7322
7323 a = construct_array(d, 4, INT4OID, sizeof(int32), true, 'i');
7324 #if POSTGIS_PGSQL_VERSION < 120
7325 ufc_info.arg[2] = PointerGetDatum(a);
7326 ufc_info.argnull[2] = FALSE;
7327 #else
7328 ufc_info->args[2].value = PointerGetDatum(a);
7329 ufc_info->args[2].isnull = FALSE;
7330 #endif
7331 }
7332
7333 #if POSTGIS_PGSQL_VERSION < 120
7334 datum = FunctionCallInvoke(&ufc_info);
7335
7336 /* result is not null*/
7337 if (!ufc_info.isnull) {
7338 haspixel = 1;
7339 pixel = DatumGetFloat8(datum);
7340 }
7341 #else
7342 datum = FunctionCallInvoke(ufc_info);
7343
7344 /* result is not null*/
7345 if (!ufc_info->isnull)
7346 {
7347 haspixel = 1;
7348 pixel = DatumGetFloat8(datum);
7349 }
7350 #endif
7351 } break;
7352 }
7353
7354 /* burn pixel if haspixel != 0 */
7355 if (haspixel) {
7356 if (rt_band_set_pixel(band, x, y, pixel, NULL) != ES_NONE) {
7357
7358 if (calltype == TEXTOID) {
7359 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7360 SPI_finish();
7361 }
7362
7363 for (k = 0; k < set_count; k++) {
7364 if (_rast[k] != NULL)
7365 rt_raster_destroy(_rast[k]);
7366 if (pgrastpos[k] != -1)
7367 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7368 }
7369 rt_raster_destroy(raster);
7370
7371 elog(ERROR, "RASTER_mapAlgebra2: Could not set pixel value of output raster");
7372 PG_RETURN_NULL();
7373 }
7374 }
7375
7376 POSTGIS_RT_DEBUGF(5, "(x, y, val) = (%d, %d, %f)", x, y, haspixel ? pixel : nodataval);
7377
7378 } /* y: height */
7379 } /* x: width */
7380 }
7381
7382 /* CLEANUP */
7383 if (calltype == TEXTOID) {
7384 for (i = 0; i < spi_count; i++) {
7385 if (spi_plan[i] != NULL) SPI_freeplan(spi_plan[i]);
7386 }
7387 SPI_finish();
7388 }
7389
7390 for (k = 0; k < set_count; k++) {
7391 if (_rast[k] != NULL)
7392 rt_raster_destroy(_rast[k]);
7393 if (pgrastpos[k] != -1)
7394 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7395 }
7396
7397 pgrtn = rt_raster_serialize(raster);
7398 rt_raster_destroy(raster);
7399 if (!pgrtn) PG_RETURN_NULL();
7400
7401 POSTGIS_RT_DEBUG(3, "Finished RASTER_mapAlgebra2");
7402
7403 SET_VARSIZE(pgrtn, pgrtn->size);
7404 PG_RETURN_POINTER(pgrtn);
7405 }
7406