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