PostGIS  2.2.7dev-r@@SVN_REVISION@@
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@keybit.net>
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>
40 #include <catalog/pg_type.h> /* for INT2OID, INT4OID, FLOAT4OID, FLOAT8OID and TEXTOID */
41 #include <executor/executor.h> /* for GetAttributeByName */
42 
43 #include "../../postgis_config.h"
44 #include "lwgeom_pg.h"
45 
46 #include "rtpostgis.h"
47 #include "rtpg_internal.h"
48 
49 /* n-raster MapAlgebra */
50 Datum RASTER_nMapAlgebra(PG_FUNCTION_ARGS);
51 Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS);
52 
53 /* raster union aggregate */
54 Datum RASTER_union_transfn(PG_FUNCTION_ARGS);
55 Datum RASTER_union_finalfn(PG_FUNCTION_ARGS);
56 
57 /* raster clip */
58 Datum RASTER_clip(PG_FUNCTION_ARGS);
59 
60 /* reclassify specified bands of a raster */
61 Datum RASTER_reclass(PG_FUNCTION_ARGS);
62 
63 /* apply colormap to specified band of a raster */
64 Datum RASTER_colorMap(PG_FUNCTION_ARGS);
65 
66 /* one-raster MapAlgebra */
67 Datum RASTER_mapAlgebraExpr(PG_FUNCTION_ARGS);
68 Datum RASTER_mapAlgebraFct(PG_FUNCTION_ARGS);
69 
70 /* one-raster neighborhood MapAlgebra */
71 Datum RASTER_mapAlgebraFctNgb(PG_FUNCTION_ARGS);
72 
73 /* two-raster MapAlgebra */
74 Datum RASTER_mapAlgebra2(PG_FUNCTION_ARGS);
75 
76 /* ---------------------------------------------------------------- */
77 /* n-raster MapAlgebra */
78 /* ---------------------------------------------------------------- */
79 
80 typedef struct {
81  Oid ufc_noid;
83  FmgrInfo ufl_info;
84  FunctionCallInfoData ufc_info;
86 
89  int numraster;
92  uint8_t *isempty; /* flag indicating if raster is empty */
93  uint8_t *ownsdata; /* is the raster self owned or just a pointer to another raster */
94  int *nband; /* source raster's band index, 0-based */
95  uint8_t *hasband; /* does source raster have band at index nband? */
96 
97  rt_pixtype pixtype; /* output raster band's pixel type */
98  int hasnodata; /* NODATA flag */
99  double nodataval; /* NODATA value */
100 
101  int distance[2]; /* distance in X and Y axis */
102 
103  rt_extenttype extenttype; /* ouput raster's extent type */
104  rt_pgraster *pgcextent; /* custom extent of type rt_pgraster */
105  rt_raster cextent; /* custom extent of type rt_raster */
106  rt_mask mask; /* mask for the nmapalgebra operation */
107 
109 };
110 
111 static rtpg_nmapalgebra_arg rtpg_nmapalgebra_arg_init() {
112  rtpg_nmapalgebra_arg arg = NULL;
113 
114  arg = palloc(sizeof(struct rtpg_nmapalgebra_arg_t));
115  if (arg == NULL) {
116  elog(ERROR, "rtpg_nmapalgebra_arg_init: Could not allocate memory for arguments");
117  return NULL;
118  }
119 
120  arg->numraster = 0;
121  arg->pgraster = NULL;
122  arg->raster = NULL;
123  arg->isempty = NULL;
124  arg->ownsdata = NULL;
125  arg->nband = NULL;
126  arg->hasband = NULL;
127 
128  arg->pixtype = PT_END;
129  arg->hasnodata = 1;
130  arg->nodataval = 0;
131 
132  arg->distance[0] = 0;
133  arg->distance[1] = 0;
134 
136  arg->pgcextent = NULL;
137  arg->cextent = NULL;
138  arg->mask = NULL;
139 
140  arg->callback.ufc_noid = InvalidOid;
141  arg->callback.ufc_rettype = InvalidOid;
142 
143  return arg;
144 }
145 
146 static void rtpg_nmapalgebra_arg_destroy(rtpg_nmapalgebra_arg arg) {
147  int i = 0;
148 
149  if (arg->raster != NULL) {
150  for (i = 0; i < arg->numraster; i++) {
151  if (arg->raster[i] == NULL || !arg->ownsdata[i])
152  continue;
153 
154  rt_raster_destroy(arg->raster[i]);
155  }
156 
157  pfree(arg->raster);
158  pfree(arg->pgraster);
159  pfree(arg->isempty);
160  pfree(arg->ownsdata);
161  pfree(arg->nband);
162  }
163 
164  if (arg->cextent != NULL)
166  if( arg->mask != NULL )
167  pfree(arg->mask);
168 
169  pfree(arg);
170 }
171 
172 static int rtpg_nmapalgebra_rastbandarg_process(rtpg_nmapalgebra_arg arg, ArrayType *array, int *allnull, int *allempty, int *noband) {
173  Oid etype;
174  Datum *e;
175  bool *nulls;
176  int16 typlen;
177  bool typbyval;
178  char typalign;
179  int n = 0;
180 
181  HeapTupleHeader tup;
182  bool isnull;
183  Datum tupv;
184 
185  int i;
186  int j;
187  int nband;
188 
189  if (arg == NULL || array == NULL) {
190  elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: NULL values not permitted for parameters");
191  return 0;
192  }
193 
194  etype = ARR_ELEMTYPE(array);
195  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
196 
197  deconstruct_array(
198  array,
199  etype,
200  typlen, typbyval, typalign,
201  &e, &nulls, &n
202  );
203 
204  if (!n) {
205  elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Invalid argument for rastbandarg");
206  return 0;
207  }
208 
209  /* prep arg */
210  arg->numraster = n;
211  arg->pgraster = palloc(sizeof(rt_pgraster *) * arg->numraster);
212  arg->raster = palloc(sizeof(rt_raster) * arg->numraster);
213  arg->isempty = palloc(sizeof(uint8_t) * arg->numraster);
214  arg->ownsdata = palloc(sizeof(uint8_t) * arg->numraster);
215  arg->nband = palloc(sizeof(int) * arg->numraster);
216  arg->hasband = palloc(sizeof(uint8_t) * arg->numraster);
217  arg->mask = palloc(sizeof(struct rt_mask_t));
218  if (
219  arg->pgraster == NULL ||
220  arg->raster == NULL ||
221  arg->isempty == NULL ||
222  arg->ownsdata == NULL ||
223  arg->nband == NULL ||
224  arg->hasband == NULL ||
225  arg->mask == NULL
226  ) {
227  elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Could not allocate memory for processing rastbandarg");
228  return 0;
229  }
230 
231  *allnull = 0;
232  *allempty = 0;
233  *noband = 0;
234 
235  /* process each element */
236  for (i = 0; i < n; i++) {
237  if (nulls[i]) {
238  arg->numraster--;
239  continue;
240  }
241 
242  POSTGIS_RT_DEBUGF(4, "Processing rastbandarg at index %d", i);
243 
244  arg->raster[i] = NULL;
245  arg->isempty[i] = 0;
246  arg->ownsdata[i] = 1;
247  arg->nband[i] = 0;
248  arg->hasband[i] = 0;
249 
250  /* each element is a tuple */
251  tup = (HeapTupleHeader) DatumGetPointer(e[i]);
252  if (NULL == tup) {
253  elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Invalid argument for rastbandarg at index %d", i);
254  return 0;
255  }
256 
257  /* first element, raster */
258  POSTGIS_RT_DEBUG(4, "Processing first element (raster)");
259  tupv = GetAttributeByName(tup, "rast", &isnull);
260  if (isnull) {
261  elog(NOTICE, "First argument (nband) of rastbandarg at index %d is NULL. Assuming NULL raster", i);
262  arg->isempty[i] = 1;
263  arg->ownsdata[i] = 0;
264 
265  (*allnull)++;
266  (*allempty)++;
267  (*noband)++;
268 
269  continue;
270  }
271 
272  arg->pgraster[i] = (rt_pgraster *) PG_DETOAST_DATUM(tupv);
273 
274  /* see if this is a copy of an existing pgraster */
275  for (j = 0; j < i; j++) {
276  if (!arg->isempty[j] && (arg->pgraster[i] == arg->pgraster[j])) {
277  POSTGIS_RT_DEBUG(4, "raster matching existing same raster found");
278  arg->raster[i] = arg->raster[j];
279  arg->ownsdata[i] = 0;
280  break;
281  }
282  }
283 
284  if (arg->ownsdata[i]) {
285  POSTGIS_RT_DEBUG(4, "deserializing raster");
286  arg->raster[i] = rt_raster_deserialize(arg->pgraster[i], FALSE);
287  if (arg->raster[i] == NULL) {
288  elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Could not deserialize raster at index %d", i);
289  return 0;
290  }
291  }
292 
293  /* is raster empty? */
294  arg->isempty[i] = rt_raster_is_empty(arg->raster[i]);
295  if (arg->isempty[i]) {
296  (*allempty)++;
297  (*noband)++;
298 
299  continue;
300  }
301 
302  /* second element, nband */
303  POSTGIS_RT_DEBUG(4, "Processing second element (nband)");
304  tupv = GetAttributeByName(tup, "nband", &isnull);
305  if (isnull) {
306  nband = 1;
307  elog(NOTICE, "First argument (nband) of rastbandarg at index %d is NULL. Assuming nband = %d", i, nband);
308  }
309  else
310  nband = DatumGetInt32(tupv);
311 
312  if (nband < 1) {
313  elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Band number provided for rastbandarg at index %d must be greater than zero (1-based)", i);
314  return 0;
315  }
316 
317  arg->nband[i] = nband - 1;
318  arg->hasband[i] = rt_raster_has_band(arg->raster[i], arg->nband[i]);
319  if (!arg->hasband[i]) {
320  (*noband)++;
321  POSTGIS_RT_DEBUGF(4, "Band at index %d not found in raster", nband);
322  }
323  }
324 
325  if (arg->numraster < n) {
326  arg->pgraster = repalloc(arg->pgraster, sizeof(rt_pgraster *) * arg->numraster);
327  arg->raster = repalloc(arg->raster, sizeof(rt_raster) * arg->numraster);
328  arg->isempty = repalloc(arg->isempty, sizeof(uint8_t) * arg->numraster);
329  arg->ownsdata = repalloc(arg->ownsdata, sizeof(uint8_t) * arg->numraster);
330  arg->nband = repalloc(arg->nband, sizeof(int) * arg->numraster);
331  arg->hasband = repalloc(arg->hasband, sizeof(uint8_t) * arg->numraster);
332  if (
333  arg->pgraster == NULL ||
334  arg->raster == NULL ||
335  arg->isempty == NULL ||
336  arg->ownsdata == NULL ||
337  arg->nband == NULL ||
338  arg->hasband == NULL
339  ) {
340  elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Could not reallocate memory for processed rastbandarg");
341  return 0;
342  }
343  }
344 
345  POSTGIS_RT_DEBUGF(4, "arg->numraster = %d", arg->numraster);
346 
347  return 1;
348 }
349 
350 /*
351  Callback for RASTER_nMapAlgebra
352 */
354  rt_iterator_arg arg, void *userarg,
355  double *value, int *nodata
356 ) {
358 
359  int16 typlen;
360  bool typbyval;
361  char typalign;
362 
363  ArrayType *mdValues = NULL;
364  Datum *_values = NULL;
365  bool *_nodata = NULL;
366 
367  ArrayType *mdPos = NULL;
368  Datum *_pos = NULL;
369  bool *_null = NULL;
370 
371  int i = 0;
372  int x = 0;
373  int y = 0;
374  int z = 0;
375  int dim[3] = {0};
376  int lbound[3] = {1, 1, 1};
377  Datum datum = (Datum) NULL;
378 
379  if (arg == NULL)
380  return 0;
381 
382  *value = 0;
383  *nodata = 0;
384 
385  dim[0] = arg->rasters;
386  dim[1] = arg->rows;
387  dim[2] = arg->columns;
388 
389  _values = palloc(sizeof(Datum) * arg->rasters * arg->rows * arg->columns);
390  _nodata = palloc(sizeof(bool) * arg->rasters * arg->rows * arg->columns);
391  if (_values == NULL || _nodata == NULL) {
392  elog(ERROR, "rtpg_nmapalgebra_callback: Could not allocate memory for values array");
393  return 0;
394  }
395 
396  /* build mdValues */
397  i = 0;
398  /* raster */
399  for (z = 0; z < arg->rasters; z++) {
400  /* Y axis */
401  for (y = 0; y < arg->rows; y++) {
402  /* X axis */
403  for (x = 0; x < arg->columns; x++) {
404  POSTGIS_RT_DEBUGF(4, "(z, y ,x) = (%d, %d, %d)", z, y, x);
405  POSTGIS_RT_DEBUGF(4, "(value, nodata) = (%f, %d)", arg->values[z][y][x], arg->nodata[z][y][x]);
406 
407  _nodata[i] = (bool) arg->nodata[z][y][x];
408  if (!_nodata[i])
409  _values[i] = Float8GetDatum(arg->values[z][y][x]);
410  else
411  _values[i] = (Datum) NULL;
412 
413  i++;
414  }
415  }
416  }
417 
418  /* info about the type of item in the multi-dimensional array (float8). */
419  get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
420 
421  /* construct mdValues */
422  mdValues = construct_md_array(
423  _values, _nodata,
424  3, dim, lbound,
425  FLOAT8OID,
426  typlen, typbyval, typalign
427  );
428  pfree(_nodata);
429  pfree(_values);
430 
431  _pos = palloc(sizeof(Datum) * (arg->rasters + 1) * 2);
432  _null = palloc(sizeof(bool) * (arg->rasters + 1) * 2);
433  if (_pos == NULL || _null == NULL) {
434  pfree(mdValues);
435  elog(ERROR, "rtpg_nmapalgebra_callback: Could not allocate memory for position array");
436  return 0;
437  }
438  memset(_null, 0, sizeof(bool) * (arg->rasters + 1) * 2);
439 
440  /* build mdPos */
441  i = 0;
442  _pos[i] = arg->dst_pixel[0] + 1;
443  i++;
444  _pos[i] = arg->dst_pixel[1] + 1;
445  i++;
446 
447  for (z = 0; z < arg->rasters; z++) {
448  _pos[i] = arg->src_pixel[z][0] + 1;
449  i++;
450 
451  _pos[i] = arg->src_pixel[z][1] + 1;
452  i++;
453  }
454 
455  /* info about the type of item in the multi-dimensional array (int4). */
456  get_typlenbyvalalign(INT4OID, &typlen, &typbyval, &typalign);
457 
458  /* reuse dim and lbound, just tweak to what we need */
459  dim[0] = arg->rasters + 1;
460  dim[1] = 2;
461  lbound[0] = 0;
462 
463  /* construct mdPos */
464  mdPos = construct_md_array(
465  _pos, _null,
466  2, dim, lbound,
467  INT4OID,
468  typlen, typbyval, typalign
469  );
470  pfree(_pos);
471  pfree(_null);
472 
473  callback->ufc_info.arg[0] = PointerGetDatum(mdValues);
474  callback->ufc_info.arg[1] = PointerGetDatum(mdPos);
475 
476  /* call user callback function */
477  datum = FunctionCallInvoke(&(callback->ufc_info));
478  pfree(mdValues);
479  pfree(mdPos);
480 
481  /* result is not null*/
482  if (!callback->ufc_info.isnull) {
483  switch (callback->ufc_rettype) {
484  case FLOAT8OID:
485  *value = DatumGetFloat8(datum);
486  break;
487  case FLOAT4OID:
488  *value = (double) DatumGetFloat4(datum);
489  break;
490  case INT4OID:
491  *value = (double) DatumGetInt32(datum);
492  break;
493  case INT2OID:
494  *value = (double) DatumGetInt16(datum);
495  break;
496  }
497  }
498  else
499  *nodata = 1;
500 
501  return 1;
502 }
503 
504 /*
505  ST_MapAlgebra for n rasters
506 */
508 Datum RASTER_nMapAlgebra(PG_FUNCTION_ARGS)
509 {
510  rtpg_nmapalgebra_arg arg = NULL;
511  rt_iterator itrset;
512  ArrayType *maskArray;
513  Oid etype;
514  Datum *maskElements;
515  bool *maskNulls;
516  int16 typlen;
517  bool typbyval;
518  char typalign;
519  int ndims = 0;
520  int num;
521  int *maskDims;
522  int x,y;
523 
524 
525  int i = 0;
526  int noerr = 0;
527  int allnull = 0;
528  int allempty = 0;
529  int noband = 0;
530 
531  rt_raster raster = NULL;
532  rt_band band = NULL;
533  rt_pgraster *pgraster = NULL;
534 
535  POSTGIS_RT_DEBUG(3, "Starting...");
536 
537  if (PG_ARGISNULL(0))
538  PG_RETURN_NULL();
539 
540  /* init argument struct */
542  if (arg == NULL) {
543  elog(ERROR, "RASTER_nMapAlgebra: Could not initialize argument structure");
544  PG_RETURN_NULL();
545  }
546 
547  /* let helper function process rastbandarg (0) */
548  if (!rtpg_nmapalgebra_rastbandarg_process(arg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
550  elog(ERROR, "RASTER_nMapAlgebra: Could not process rastbandarg");
551  PG_RETURN_NULL();
552  }
553 
554  POSTGIS_RT_DEBUGF(4, "allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
555 
556  /* all rasters are NULL, return NULL */
557  if (allnull == arg->numraster) {
558  elog(NOTICE, "All input rasters are NULL. Returning NULL");
560  PG_RETURN_NULL();
561  }
562 
563  /* pixel type (2) */
564  if (!PG_ARGISNULL(2)) {
565  char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
566 
567  /* Get the pixel type index */
568  arg->pixtype = rt_pixtype_index_from_name(pixtypename);
569  if (arg->pixtype == PT_END) {
571  elog(ERROR, "RASTER_nMapAlgebra: Invalid pixel type: %s", pixtypename);
572  PG_RETURN_NULL();
573  }
574  }
575 
576  /* distancex (3) */
577  if (!PG_ARGISNULL(3)){
578  arg->distance[0] = PG_GETARG_INT32(3);
579  }else{
580  arg->distance[0] = 0;
581  }
582  /* distancey (4) */
583  if (!PG_ARGISNULL(4)){
584  arg->distance[1] = PG_GETARG_INT32(4);
585  }else{
586  arg->distance[1] = 0;
587  }
588  if (arg->distance[0] < 0 || arg->distance[1] < 0) {
590  elog(ERROR, "RASTER_nMapAlgebra: Distance for X and Y axis must be greater than or equal to zero");
591  PG_RETURN_NULL();
592  }
593 
594  /* extent type (5) */
595  if (!PG_ARGISNULL(5)) {
596  char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(5))));
597  arg->extenttype = rt_util_extent_type(extenttypename);
598  }
599  POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->extenttype);
600 
601  /* custom extent (6) */
602  if (arg->extenttype == ET_CUSTOM) {
603  if (PG_ARGISNULL(6)) {
604  elog(NOTICE, "Custom extent is NULL. Returning NULL");
606  PG_RETURN_NULL();
607  }
608 
609  arg->pgcextent = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(6));
610 
611  /* only need the raster header */
613  if (arg->cextent == NULL) {
615  elog(ERROR, "RASTER_nMapAlgebra: Could not deserialize custom extent");
616  PG_RETURN_NULL();
617  }
618  else if (rt_raster_is_empty(arg->cextent)) {
619  elog(NOTICE, "Custom extent is an empty raster. Returning empty raster");
621 
622  raster = rt_raster_new(0, 0);
623  if (raster == NULL) {
624  elog(ERROR, "RASTER_nMapAlgebra: Could not create empty raster");
625  PG_RETURN_NULL();
626  }
627 
628  pgraster = rt_raster_serialize(raster);
629  rt_raster_destroy(raster);
630  if (!pgraster) PG_RETURN_NULL();
631 
632  SET_VARSIZE(pgraster, pgraster->size);
633  PG_RETURN_POINTER(pgraster);
634  }
635  }
636 
637  /* mask (7) */
638 
639  if( PG_ARGISNULL(7) ){
640  pfree(arg->mask);
641  arg->mask = NULL;
642  }
643  else {
644  maskArray = PG_GETARG_ARRAYTYPE_P(7);
645  etype = ARR_ELEMTYPE(maskArray);
646  get_typlenbyvalalign(etype,&typlen,&typbyval,&typalign);
647 
648  switch (etype) {
649  case FLOAT4OID:
650  case FLOAT8OID:
651  break;
652  default:
654  elog(ERROR,"RASTER_nMapAlgebra: Mask data type must be FLOAT8 or FLOAT4");
655  PG_RETURN_NULL();
656  }
657 
658  ndims = ARR_NDIM(maskArray);
659 
660  if (ndims != 2) {
661  elog(ERROR, "RASTER_nMapAlgebra: Mask Must be a 2D array");
663  PG_RETURN_NULL();
664  }
665 
666  maskDims = ARR_DIMS(maskArray);
667 
668  if (maskDims[0] % 2 == 0 || maskDims[1] % 2 == 0) {
669  elog(ERROR,"RASTER_nMapAlgebra: Mask dimensions must be odd");
671  PG_RETURN_NULL();
672  }
673 
674  deconstruct_array(
675  maskArray,
676  etype,
677  typlen, typbyval,typalign,
678  &maskElements,&maskNulls,&num
679  );
680 
681  if (num < 1 || num != (maskDims[0] * maskDims[1])) {
682  if (num) {
683  pfree(maskElements);
684  pfree(maskNulls);
685  }
686  elog(ERROR, "RASTER_nMapAlgebra: Could not deconstruct new values array");
688  PG_RETURN_NULL();
689  }
690 
691  /* allocate mem for mask array */
692  arg->mask->values = palloc(sizeof(double*)* maskDims[0]);
693  arg->mask->nodata = palloc(sizeof(int*)*maskDims[0]);
694  for (i = 0; i < maskDims[0]; i++) {
695  arg->mask->values[i] = (double*) palloc(sizeof(double) * maskDims[1]);
696  arg->mask->nodata[i] = (int*) palloc(sizeof(int) * maskDims[1]);
697  }
698 
699  /* place values in to mask */
700  i = 0;
701  for (y = 0; y < maskDims[0]; y++) {
702  for (x = 0; x < maskDims[1]; x++) {
703  if (maskNulls[i]) {
704  arg->mask->values[y][x] = 0;
705  arg->mask->nodata[y][x] = 1;
706  }
707  else {
708  switch (etype) {
709  case FLOAT4OID:
710  arg->mask->values[y][x] = (double) DatumGetFloat4(maskElements[i]);
711  arg->mask->nodata[y][x] = 0;
712  break;
713  case FLOAT8OID:
714  arg->mask->values[y][x] = (double) DatumGetFloat8(maskElements[i]);
715  arg->mask->nodata[y][x] = 0;
716  }
717  }
718  i++;
719  }
720  }
721 
722  /*set mask dimensions*/
723  arg->mask->dimx = maskDims[0];
724  arg->mask->dimy = maskDims[1];
725  if (maskDims[0] == 1 && maskDims[1] == 1) {
726  arg->distance[0] = 0;
727  arg->distance[1] = 0;
728  }
729  else {
730  arg->distance[0] = maskDims[0] % 2;
731  arg->distance[1] = maskDims[1] % 2;
732  }
733  }/*end if else argisnull*/
734 
735  /* (8) weighted boolean */
736  if (PG_ARGISNULL(8) || !PG_GETARG_BOOL(8)) {
737  if (arg->mask != NULL)
738  arg->mask->weighted = 0;
739  }else{
740  if(arg->mask !=NULL)
741  arg->mask->weighted = 1;
742  }
743 
744  noerr = 1;
745 
746  /* all rasters are empty, return empty raster */
747  if (allempty == arg->numraster) {
748  elog(NOTICE, "All input rasters are empty. Returning empty raster");
749  noerr = 0;
750  }
751  /* all rasters don't have indicated band, return empty raster */
752  else if (noband == arg->numraster) {
753  elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
754  noerr = 0;
755  }
756  if (!noerr) {
758 
759  raster = rt_raster_new(0, 0);
760  if (raster == NULL) {
761  elog(ERROR, "RASTER_nMapAlgebra: Could not create empty raster");
762  PG_RETURN_NULL();
763  }
764 
765  pgraster = rt_raster_serialize(raster);
766  rt_raster_destroy(raster);
767  if (!pgraster) PG_RETURN_NULL();
768 
769  SET_VARSIZE(pgraster, pgraster->size);
770  PG_RETURN_POINTER(pgraster);
771  }
772 
773  /* do regprocedure last (1) */
774  if (!PG_ARGISNULL(1) || get_fn_expr_argtype(fcinfo->flinfo, 1) == REGPROCEDUREOID) {
775  POSTGIS_RT_DEBUG(4, "processing callbackfunc");
776  arg->callback.ufc_noid = PG_GETARG_OID(1);
777 
778  /* get function info */
779  fmgr_info(arg->callback.ufc_noid, &(arg->callback.ufl_info));
780 
781  /* function cannot return set */
782  noerr = 0;
783  if (arg->callback.ufl_info.fn_retset) {
784  noerr = 1;
785  }
786  /* function should have correct # of args */
787  else if (arg->callback.ufl_info.fn_nargs != 3) {
788  noerr = 2;
789  }
790 
791  /* check that callback function return type is supported */
792  if (
793  get_func_result_type(
794  arg->callback.ufc_noid,
795  &(arg->callback.ufc_rettype),
796  NULL
797  ) != TYPEFUNC_SCALAR
798  ) {
799  noerr = 3;
800  }
801 
802  if (!(
803  arg->callback.ufc_rettype == FLOAT8OID ||
804  arg->callback.ufc_rettype == FLOAT4OID ||
805  arg->callback.ufc_rettype == INT4OID ||
806  arg->callback.ufc_rettype == INT2OID
807  )) {
808  noerr = 4;
809  }
810 
811  /*
812  TODO: consider adding checks of the userfunction parameters
813  should be able to use get_fn_expr_argtype() of fmgr.c
814  */
815 
816  if (noerr != 0) {
818  switch (noerr) {
819  case 4:
820  elog(ERROR, "RASTER_nMapAlgebra: Function provided must return a double precision, float, int or smallint");
821  break;
822  case 3:
823  elog(ERROR, "RASTER_nMapAlgebra: Function provided must return scalar (double precision, float, int, smallint)");
824  break;
825  case 2:
826  elog(ERROR, "RASTER_nMapAlgebra: Function provided must have three input parameters");
827  break;
828  case 1:
829  elog(ERROR, "RASTER_nMapAlgebra: Function provided must return double precision, not resultset");
830  break;
831  }
832  PG_RETURN_NULL();
833  }
834 
835  if (func_volatile(arg->callback.ufc_noid) == 'v')
836  elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
837 
838  /* prep function call data */
839 #if POSTGIS_PGSQL_VERSION > 90
840  InitFunctionCallInfoData(arg->callback.ufc_info, &(arg->callback.ufl_info), arg->callback.ufl_info.fn_nargs, InvalidOid, NULL, NULL);
841 #else
842  InitFunctionCallInfoData(arg->callback.ufc_info, &(arg->callback.ufl_info), arg->callback.ufl_info.fn_nargs, NULL, NULL);
843 #endif
844  memset(arg->callback.ufc_info.argnull, FALSE, sizeof(bool) * arg->callback.ufl_info.fn_nargs);
845 
846  /* userargs (7) */
847  if (!PG_ARGISNULL(9))
848  arg->callback.ufc_info.arg[2] = PG_GETARG_DATUM(9);
849  else {
850  if (arg->callback.ufl_info.fn_strict) {
851  /* build and assign an empty TEXT array */
852  /* TODO: manually free the empty array? */
853  arg->callback.ufc_info.arg[2] = PointerGetDatum(
854  construct_empty_array(TEXTOID)
855  );
856  arg->callback.ufc_info.argnull[2] = FALSE;
857  }
858  else {
859  arg->callback.ufc_info.arg[2] = (Datum) NULL;
860  arg->callback.ufc_info.argnull[2] = TRUE;
861  }
862  }
863  }
864  else {
866  elog(ERROR, "RASTER_nMapAlgebra: callbackfunc must be provided");
867  PG_RETURN_NULL();
868  }
869 
870  /* determine nodataval and possibly pixtype */
871  /* band to check */
872  switch (arg->extenttype) {
873  case ET_LAST:
874  i = arg->numraster - 1;
875  break;
876  case ET_SECOND:
877  if (arg->numraster > 1) {
878  i = 1;
879  break;
880  }
881  default:
882  i = 0;
883  break;
884  }
885  /* find first viable band */
886  if (!arg->hasband[i]) {
887  for (i = 0; i < arg->numraster; i++) {
888  if (arg->hasband[i])
889  break;
890  }
891  if (i >= arg->numraster)
892  i = arg->numraster - 1;
893  }
894  band = rt_raster_get_band(arg->raster[i], arg->nband[i]);
895 
896  /* set pixel type if PT_END */
897  if (arg->pixtype == PT_END)
898  arg->pixtype = rt_band_get_pixtype(band);
899 
900  /* set hasnodata and nodataval */
901  arg->hasnodata = 1;
902  if (rt_band_get_hasnodata_flag(band))
903  rt_band_get_nodata(band, &(arg->nodataval));
904  else
905  arg->nodataval = rt_band_get_min_value(band);
906 
907  POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->pixtype), arg->hasnodata, arg->nodataval);
908 
909  /* init itrset */
910  itrset = palloc(sizeof(struct rt_iterator_t) * arg->numraster);
911  if (itrset == NULL) {
913  elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
914  PG_RETURN_NULL();
915  }
916 
917  /* set itrset */
918  for (i = 0; i < arg->numraster; i++) {
919  itrset[i].raster = arg->raster[i];
920  itrset[i].nband = arg->nband[i];
921  itrset[i].nbnodata = 1;
922  }
923 
924  /* pass everything to iterator */
925  noerr = rt_raster_iterator(
926  itrset, arg->numraster,
927  arg->extenttype, arg->cextent,
928  arg->pixtype,
929  arg->hasnodata, arg->nodataval,
930  arg->distance[0], arg->distance[1],
931  arg->mask,
932  &(arg->callback),
934  &raster
935  );
936 
937  /* cleanup */
938  pfree(itrset);
940 
941  if (noerr != ES_NONE) {
942  elog(ERROR, "RASTER_nMapAlgebra: Could not run raster iterator function");
943  PG_RETURN_NULL();
944  }
945  else if (raster == NULL)
946  PG_RETURN_NULL();
947 
948  pgraster = rt_raster_serialize(raster);
949  rt_raster_destroy(raster);
950 
951  POSTGIS_RT_DEBUG(3, "Finished");
952 
953  if (!pgraster)
954  PG_RETURN_NULL();
955 
956  SET_VARSIZE(pgraster, pgraster->size);
957  PG_RETURN_POINTER(pgraster);
958 }
959 
960 /* ---------------------------------------------------------------- */
961 /* expression ST_MapAlgebra for n rasters */
962 /* ---------------------------------------------------------------- */
963 
964 typedef struct {
966 
967  struct {
968  SPIPlanPtr spi_plan;
969  uint32_t spi_argcount;
970  uint8_t *spi_argpos;
971 
972  int hasval;
973  double val;
974  } expr[3];
975 
976  struct {
977  int hasval;
978  double val;
979  } nodatanodata;
980 
981  struct {
982  int count;
983  char **val;
984  } kw;
985 
987 
990  rtpg_nmapalgebra_arg bandarg;
991 
993 };
994 
995 static rtpg_nmapalgebraexpr_arg rtpg_nmapalgebraexpr_arg_init(int cnt, char **kw) {
996  rtpg_nmapalgebraexpr_arg arg = NULL;
997  int i = 0;
998 
999  arg = palloc(sizeof(struct rtpg_nmapalgebraexpr_arg_t));
1000  if (arg == NULL) {
1001  elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arguments");
1002  return NULL;
1003  }
1004 
1006  if (arg->bandarg == NULL) {
1007  elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arg->bandarg");
1008  return NULL;
1009  }
1010 
1011  arg->callback.kw.count = cnt;
1012  arg->callback.kw.val = kw;
1013 
1014  arg->callback.exprcount = 3;
1015  for (i = 0; i < arg->callback.exprcount; i++) {
1016  arg->callback.expr[i].spi_plan = NULL;
1017  arg->callback.expr[i].spi_argcount = 0;
1018  arg->callback.expr[i].spi_argpos = palloc(cnt * sizeof(uint8_t));
1019  if (arg->callback.expr[i].spi_argpos == NULL) {
1020  elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for spi_argpos");
1021  return NULL;
1022  }
1023  memset(arg->callback.expr[i].spi_argpos, 0, sizeof(uint8_t) * cnt);
1024  arg->callback.expr[i].hasval = 0;
1025  arg->callback.expr[i].val = 0;
1026  }
1027 
1028  arg->callback.nodatanodata.hasval = 0;
1029  arg->callback.nodatanodata.val = 0;
1030 
1031  return arg;
1032 }
1033 
1034 static void rtpg_nmapalgebraexpr_arg_destroy(rtpg_nmapalgebraexpr_arg arg) {
1035  int i = 0;
1036 
1038 
1039  for (i = 0; i < arg->callback.exprcount; i++) {
1040  if (arg->callback.expr[i].spi_plan)
1041  SPI_freeplan(arg->callback.expr[i].spi_plan);
1042  if (arg->callback.kw.count)
1043  pfree(arg->callback.expr[i].spi_argpos);
1044  }
1045 
1046  pfree(arg);
1047 }
1048 
1050  rt_iterator_arg arg, void *userarg,
1051  double *value, int *nodata
1052 ) {
1054  SPIPlanPtr plan = NULL;
1055  int i = 0;
1056  int id = -1;
1057 
1058  if (arg == NULL)
1059  return 0;
1060 
1061  *value = 0;
1062  *nodata = 0;
1063 
1064  /* 2 raster */
1065  if (arg->rasters > 1) {
1066  /* nodata1 = 1 AND nodata2 = 1, nodatanodataval */
1067  if (arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1068  if (callback->nodatanodata.hasval)
1069  *value = callback->nodatanodata.val;
1070  else
1071  *nodata = 1;
1072  }
1073  /* nodata1 = 1 AND nodata2 != 1, nodata1expr */
1074  else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0]) {
1075  id = 1;
1076  if (callback->expr[id].hasval)
1077  *value = callback->expr[id].val;
1078  else if (callback->expr[id].spi_plan)
1079  plan = callback->expr[id].spi_plan;
1080  else
1081  *nodata = 1;
1082  }
1083  /* nodata1 != 1 AND nodata2 = 1, nodata2expr */
1084  else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1085  id = 2;
1086  if (callback->expr[id].hasval)
1087  *value = callback->expr[id].val;
1088  else if (callback->expr[id].spi_plan)
1089  plan = callback->expr[id].spi_plan;
1090  else
1091  *nodata = 1;
1092  }
1093  /* expression */
1094  else {
1095  id = 0;
1096  if (callback->expr[id].hasval)
1097  *value = callback->expr[id].val;
1098  else if (callback->expr[id].spi_plan)
1099  plan = callback->expr[id].spi_plan;
1100  else {
1101  if (callback->nodatanodata.hasval)
1102  *value = callback->nodatanodata.val;
1103  else
1104  *nodata = 1;
1105  }
1106  }
1107  }
1108  /* 1 raster */
1109  else {
1110  /* nodata = 1, nodata1expr */
1111  if (arg->nodata[0][0][0]) {
1112  id = 1;
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  *nodata = 1;
1119  }
1120  /* expression */
1121  else {
1122  id = 0;
1123  if (callback->expr[id].hasval)
1124  *value = callback->expr[id].val;
1125  else if (callback->expr[id].spi_plan)
1126  plan = callback->expr[id].spi_plan;
1127  else {
1128  /* see if nodata1expr is available */
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  }
1138  }
1139 
1140  /* run prepared plan */
1141  if (plan != NULL) {
1142  Datum values[12];
1143  bool nulls[12];
1144  int err = 0;
1145 
1146  TupleDesc tupdesc;
1147  SPITupleTable *tuptable = NULL;
1148  HeapTuple tuple;
1149  Datum datum;
1150  bool isnull = FALSE;
1151 
1152  POSTGIS_RT_DEBUGF(4, "Running plan %d", id);
1153 
1154  /* init values and nulls */
1155  memset(values, (Datum) NULL, sizeof(Datum) * callback->kw.count);
1156  memset(nulls, FALSE, sizeof(bool) * callback->kw.count);
1157 
1158  if (callback->expr[id].spi_argcount) {
1159  int idx = 0;
1160 
1161  for (i = 0; i < callback->kw.count; i++) {
1162  idx = callback->expr[id].spi_argpos[i];
1163  if (idx < 1) continue;
1164  idx--; /* 1-based now 0-based */
1165 
1166  switch (i) {
1167  /* [rast.x] */
1168  case 0:
1169  values[idx] = Int32GetDatum(arg->src_pixel[0][0] + 1);
1170  break;
1171  /* [rast.y] */
1172  case 1:
1173  values[idx] = Int32GetDatum(arg->src_pixel[0][1] + 1);
1174  break;
1175  /* [rast.val] */
1176  case 2:
1177  /* [rast] */
1178  case 3:
1179  if (!arg->nodata[0][0][0])
1180  values[idx] = Float8GetDatum(arg->values[0][0][0]);
1181  else
1182  nulls[idx] = TRUE;
1183  break;
1184 
1185  /* [rast1.x] */
1186  case 4:
1187  values[idx] = Int32GetDatum(arg->src_pixel[0][0] + 1);
1188  break;
1189  /* [rast1.y] */
1190  case 5:
1191  values[idx] = Int32GetDatum(arg->src_pixel[0][1] + 1);
1192  break;
1193  /* [rast1.val] */
1194  case 6:
1195  /* [rast1] */
1196  case 7:
1197  if (!arg->nodata[0][0][0])
1198  values[idx] = Float8GetDatum(arg->values[0][0][0]);
1199  else
1200  nulls[idx] = TRUE;
1201  break;
1202 
1203  /* [rast2.x] */
1204  case 8:
1205  values[idx] = Int32GetDatum(arg->src_pixel[1][0] + 1);
1206  break;
1207  /* [rast2.y] */
1208  case 9:
1209  values[idx] = Int32GetDatum(arg->src_pixel[1][1] + 1);
1210  break;
1211  /* [rast2.val] */
1212  case 10:
1213  /* [rast2] */
1214  case 11:
1215  if (!arg->nodata[1][0][0])
1216  values[idx] = Float8GetDatum(arg->values[1][0][0]);
1217  else
1218  nulls[idx] = TRUE;
1219  break;
1220  }
1221 
1222  }
1223  }
1224 
1225  /* run prepared plan */
1226  err = SPI_execute_plan(plan, values, nulls, TRUE, 1);
1227  if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1228  elog(ERROR, "rtpg_nmapalgebraexpr_callback: Unexpected error when running prepared statement %d", id);
1229  return 0;
1230  }
1231 
1232  /* get output of prepared plan */
1233  tupdesc = SPI_tuptable->tupdesc;
1234  tuptable = SPI_tuptable;
1235  tuple = tuptable->vals[0];
1236 
1237  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1238  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1239  if (SPI_tuptable) SPI_freetuptable(tuptable);
1240  elog(ERROR, "rtpg_nmapalgebraexpr_callback: Could not get result of prepared statement %d", id);
1241  return 0;
1242  }
1243 
1244  if (!isnull) {
1245  *value = DatumGetFloat8(datum);
1246  POSTGIS_RT_DEBUG(4, "Getting value from Datum");
1247  }
1248  else {
1249  /* 2 raster, check nodatanodataval */
1250  if (arg->rasters > 1) {
1251  if (callback->nodatanodata.hasval)
1252  *value = callback->nodatanodata.val;
1253  else
1254  *nodata = 1;
1255  }
1256  /* 1 raster, check nodataval */
1257  else {
1258  if (callback->expr[1].hasval)
1259  *value = callback->expr[1].val;
1260  else
1261  *nodata = 1;
1262  }
1263  }
1264 
1265  if (SPI_tuptable) SPI_freetuptable(tuptable);
1266  }
1267 
1268  POSTGIS_RT_DEBUGF(4, "(value, nodata) = (%f, %d)", *value, *nodata);
1269  return 1;
1270 }
1271 
1273 Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS)
1274 {
1275  MemoryContext mainMemCtx = CurrentMemoryContext;
1276  rtpg_nmapalgebraexpr_arg arg = NULL;
1277  rt_iterator itrset;
1278  uint16_t exprpos[3] = {1, 4, 5};
1279 
1280  int i = 0;
1281  int j = 0;
1282  int k = 0;
1283 
1284  int numraster = 0;
1285  int err = 0;
1286  int allnull = 0;
1287  int allempty = 0;
1288  int noband = 0;
1289  int len = 0;
1290 
1291  TupleDesc tupdesc;
1292  SPITupleTable *tuptable = NULL;
1293  HeapTuple tuple;
1294  Datum datum;
1295  bool isnull = FALSE;
1296 
1297  rt_raster raster = NULL;
1298  rt_band band = NULL;
1299  rt_pgraster *pgraster = NULL;
1300 
1301  const int argkwcount = 12;
1302  char *argkw[] = {
1303  "[rast.x]",
1304  "[rast.y]",
1305  "[rast.val]",
1306  "[rast]",
1307  "[rast1.x]",
1308  "[rast1.y]",
1309  "[rast1.val]",
1310  "[rast1]",
1311  "[rast2.x]",
1312  "[rast2.y]",
1313  "[rast2.val]",
1314  "[rast2]"
1315  };
1316 
1317  if (PG_ARGISNULL(0))
1318  PG_RETURN_NULL();
1319 
1320  /* init argument struct */
1321  arg = rtpg_nmapalgebraexpr_arg_init(argkwcount, argkw);
1322  if (arg == NULL) {
1323  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not initialize argument structure");
1324  PG_RETURN_NULL();
1325  }
1326 
1327  /* let helper function process rastbandarg (0) */
1328  if (!rtpg_nmapalgebra_rastbandarg_process(arg->bandarg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
1330  elog(ERROR, "RASTER_nMapAlgebra: Could not process rastbandarg");
1331  PG_RETURN_NULL();
1332  }
1333 
1334  POSTGIS_RT_DEBUGF(4, "allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
1335 
1336  /* all rasters are NULL, return NULL */
1337  if (allnull == arg->bandarg->numraster) {
1338  elog(NOTICE, "All input rasters are NULL. Returning NULL");
1340  PG_RETURN_NULL();
1341  }
1342 
1343  /* only work on one or two rasters */
1344  if (arg->bandarg->numraster > 1)
1345  numraster = 2;
1346  else
1347  numraster = 1;
1348 
1349  /* pixel type (2) */
1350  if (!PG_ARGISNULL(2)) {
1351  char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
1352 
1353  /* Get the pixel type index */
1354  arg->bandarg->pixtype = rt_pixtype_index_from_name(pixtypename);
1355  if (arg->bandarg->pixtype == PT_END) {
1357  elog(ERROR, "RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
1358  PG_RETURN_NULL();
1359  }
1360  }
1361  POSTGIS_RT_DEBUGF(4, "pixeltype: %d", arg->bandarg->pixtype);
1362 
1363  /* extent type (3) */
1364  if (!PG_ARGISNULL(3)) {
1365  char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(3))));
1366  arg->bandarg->extenttype = rt_util_extent_type(extenttypename);
1367  }
1368 
1369  if (arg->bandarg->extenttype == ET_CUSTOM) {
1370  if (numraster < 2) {
1371  elog(NOTICE, "CUSTOM extent type not supported. Defaulting to FIRST");
1372  arg->bandarg->extenttype = ET_FIRST;
1373  }
1374  else {
1375  elog(NOTICE, "CUSTOM extent type not supported. Defaulting to INTERSECTION");
1377  }
1378  }
1379  else if (numraster < 2)
1380  arg->bandarg->extenttype = ET_FIRST;
1381 
1382  POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->bandarg->extenttype);
1383 
1384  /* nodatanodataval (6) */
1385  if (!PG_ARGISNULL(6)) {
1386  arg->callback.nodatanodata.hasval = 1;
1387  arg->callback.nodatanodata.val = PG_GETARG_FLOAT8(6);
1388  }
1389 
1390  err = 0;
1391  /* all rasters are empty, return empty raster */
1392  if (allempty == arg->bandarg->numraster) {
1393  elog(NOTICE, "All input rasters are empty. Returning empty raster");
1394  err = 1;
1395  }
1396  /* all rasters don't have indicated band, return empty raster */
1397  else if (noband == arg->bandarg->numraster) {
1398  elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
1399  err = 1;
1400  }
1401  if (err) {
1403 
1404  raster = rt_raster_new(0, 0);
1405  if (raster == NULL) {
1406  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create empty raster");
1407  PG_RETURN_NULL();
1408  }
1409 
1410  pgraster = rt_raster_serialize(raster);
1411  rt_raster_destroy(raster);
1412  if (!pgraster) PG_RETURN_NULL();
1413 
1414  SET_VARSIZE(pgraster, pgraster->size);
1415  PG_RETURN_POINTER(pgraster);
1416  }
1417 
1418  /* connect SPI */
1419  if (SPI_connect() != SPI_OK_CONNECT) {
1421  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
1422  PG_RETURN_NULL();
1423  }
1424 
1425  /*
1426  process expressions
1427 
1428  exprpos elements are:
1429  1 - expression => spi_plan[0]
1430  4 - nodata1expr => spi_plan[1]
1431  5 - nodata2expr => spi_plan[2]
1432  */
1433  for (i = 0; i < arg->callback.exprcount; i++) {
1434  char *expr = NULL;
1435  char *tmp = NULL;
1436  char *sql = NULL;
1437  char place[5] = "$1";
1438 
1439  if (PG_ARGISNULL(exprpos[i]))
1440  continue;
1441 
1442  expr = text_to_cstring(PG_GETARG_TEXT_P(exprpos[i]));
1443  POSTGIS_RT_DEBUGF(3, "raw expr of argument #%d: %s", exprpos[i], expr);
1444 
1445  for (j = 0, k = 1; j < argkwcount; j++) {
1446  /* attempt to replace keyword with placeholder */
1447  len = 0;
1448  tmp = rtpg_strreplace(expr, argkw[j], place, &len);
1449  pfree(expr);
1450  expr = tmp;
1451 
1452  if (len) {
1453  POSTGIS_RT_DEBUGF(4, "kw #%d (%s) at pos $%d", j, argkw[j], k);
1454  arg->callback.expr[i].spi_argcount++;
1455  arg->callback.expr[i].spi_argpos[j] = k++;
1456 
1457  sprintf(place, "$%d", k);
1458  }
1459  else
1460  arg->callback.expr[i].spi_argpos[j] = 0;
1461  }
1462 
1463  len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision");
1464  sql = (char *) palloc(len + 1);
1465  if (sql == NULL) {
1467  SPI_finish();
1468  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
1469  PG_RETURN_NULL();
1470  }
1471 
1472  strncpy(sql, "SELECT (", strlen("SELECT ("));
1473  strncpy(sql + strlen("SELECT ("), expr, strlen(expr));
1474  strncpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision"));
1475  sql[len] = '\0';
1476 
1477  POSTGIS_RT_DEBUGF(3, "sql #%d: %s", exprpos[i], sql);
1478 
1479  /* prepared plan */
1480  if (arg->callback.expr[i].spi_argcount) {
1481  Oid *argtype = (Oid *) palloc(arg->callback.expr[i].spi_argcount * sizeof(Oid));
1482  POSTGIS_RT_DEBUGF(3, "expression parameter %d is a prepared plan", exprpos[i]);
1483  if (argtype == NULL) {
1484  pfree(sql);
1486  SPI_finish();
1487  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
1488  PG_RETURN_NULL();
1489  }
1490 
1491  /* specify datatypes of parameters */
1492  for (j = 0, k = 0; j < argkwcount; j++) {
1493  if (arg->callback.expr[i].spi_argpos[j] < 1) continue;
1494 
1495  /* positions are INT4 */
1496  if (
1497  (strstr(argkw[j], "[rast.x]") != NULL) ||
1498  (strstr(argkw[j], "[rast.y]") != NULL) ||
1499  (strstr(argkw[j], "[rast1.x]") != NULL) ||
1500  (strstr(argkw[j], "[rast1.y]") != NULL) ||
1501  (strstr(argkw[j], "[rast2.x]") != NULL) ||
1502  (strstr(argkw[j], "[rast2.y]") != NULL)
1503  )
1504  argtype[k] = INT4OID;
1505  /* everything else is FLOAT8 */
1506  else
1507  argtype[k] = FLOAT8OID;
1508 
1509  k++;
1510  }
1511 
1512  arg->callback.expr[i].spi_plan = SPI_prepare(sql, arg->callback.expr[i].spi_argcount, argtype);
1513  pfree(argtype);
1514  pfree(sql);
1515 
1516  if (arg->callback.expr[i].spi_plan == NULL) {
1518  SPI_finish();
1519  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
1520  PG_RETURN_NULL();
1521  }
1522  }
1523  /* no args, just execute query */
1524  else {
1525  POSTGIS_RT_DEBUGF(3, "expression parameter %d has no args, simply executing", exprpos[i]);
1526  err = SPI_execute(sql, TRUE, 0);
1527  pfree(sql);
1528 
1529  if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1531  SPI_finish();
1532  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
1533  PG_RETURN_NULL();
1534  }
1535 
1536  /* get output of prepared plan */
1537  tupdesc = SPI_tuptable->tupdesc;
1538  tuptable = SPI_tuptable;
1539  tuple = tuptable->vals[0];
1540 
1541  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1542  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1543  if (SPI_tuptable) SPI_freetuptable(tuptable);
1545  SPI_finish();
1546  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
1547  PG_RETURN_NULL();
1548  }
1549 
1550  if (!isnull) {
1551  arg->callback.expr[i].hasval = 1;
1552  arg->callback.expr[i].val = DatumGetFloat8(datum);
1553  }
1554 
1555  if (SPI_tuptable) SPI_freetuptable(tuptable);
1556  }
1557  }
1558 
1559  /* determine nodataval and possibly pixtype */
1560  /* band to check */
1561  switch (arg->bandarg->extenttype) {
1562  case ET_LAST:
1563  case ET_SECOND:
1564  if (numraster > 1)
1565  i = 1;
1566  else
1567  i = 0;
1568  break;
1569  default:
1570  i = 0;
1571  break;
1572  }
1573  /* find first viable band */
1574  if (!arg->bandarg->hasband[i]) {
1575  for (i = 0; i < numraster; i++) {
1576  if (arg->bandarg->hasband[i])
1577  break;
1578  }
1579  if (i >= numraster)
1580  i = numraster - 1;
1581  }
1582  band = rt_raster_get_band(arg->bandarg->raster[i], arg->bandarg->nband[i]);
1583 
1584  /* set pixel type if PT_END */
1585  if (arg->bandarg->pixtype == PT_END)
1586  arg->bandarg->pixtype = rt_band_get_pixtype(band);
1587 
1588  /* set hasnodata and nodataval */
1589  arg->bandarg->hasnodata = 1;
1590  if (rt_band_get_hasnodata_flag(band))
1591  rt_band_get_nodata(band, &(arg->bandarg->nodataval));
1592  else
1593  arg->bandarg->nodataval = rt_band_get_min_value(band);
1594 
1595  POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->bandarg->pixtype), arg->bandarg->hasnodata, arg->bandarg->nodataval);
1596 
1597  /* init itrset */
1598  itrset = palloc(sizeof(struct rt_iterator_t) * numraster);
1599  if (itrset == NULL) {
1601  SPI_finish();
1602  elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
1603  PG_RETURN_NULL();
1604  }
1605 
1606  /* set itrset */
1607  for (i = 0; i < numraster; i++) {
1608  itrset[i].raster = arg->bandarg->raster[i];
1609  itrset[i].nband = arg->bandarg->nband[i];
1610  itrset[i].nbnodata = 1;
1611  }
1612 
1613  /* pass everything to iterator */
1614  err = rt_raster_iterator(
1615  itrset, numraster,
1616  arg->bandarg->extenttype, arg->bandarg->cextent,
1617  arg->bandarg->pixtype,
1618  arg->bandarg->hasnodata, arg->bandarg->nodataval,
1619  0, 0,
1620  NULL,
1621  &(arg->callback),
1623  &raster
1624  );
1625 
1626  pfree(itrset);
1628 
1629  if (err != ES_NONE) {
1630  SPI_finish();
1631  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not run raster iterator function");
1632  PG_RETURN_NULL();
1633  }
1634  else if (raster == NULL) {
1635  SPI_finish();
1636  PG_RETURN_NULL();
1637  }
1638 
1639  /* switch to prior memory context to ensure memory allocated in correct context */
1640  MemoryContextSwitchTo(mainMemCtx);
1641 
1642  pgraster = rt_raster_serialize(raster);
1643  rt_raster_destroy(raster);
1644 
1645  /* finish SPI */
1646  SPI_finish();
1647 
1648  if (!pgraster)
1649  PG_RETURN_NULL();
1650 
1651  SET_VARSIZE(pgraster, pgraster->size);
1652  PG_RETURN_POINTER(pgraster);
1653 }
1654 
1655 /* ---------------------------------------------------------------- */
1656 /* ST_Union aggregate functions */
1657 /* ---------------------------------------------------------------- */
1658 
1659 typedef enum {
1660  UT_LAST = 0,
1668 } rtpg_union_type;
1669 
1670 /* internal function translating text of UNION type to enum */
1672  assert(cutype && strlen(cutype) > 0);
1673 
1674  if (strcmp(cutype, "LAST") == 0)
1675  return UT_LAST;
1676  else if (strcmp(cutype, "FIRST") == 0)
1677  return UT_FIRST;
1678  else if (strcmp(cutype, "MIN") == 0)
1679  return UT_MIN;
1680  else if (strcmp(cutype, "MAX") == 0)
1681  return UT_MAX;
1682  else if (strcmp(cutype, "COUNT") == 0)
1683  return UT_COUNT;
1684  else if (strcmp(cutype, "SUM") == 0)
1685  return UT_SUM;
1686  else if (strcmp(cutype, "MEAN") == 0)
1687  return UT_MEAN;
1688  else if (strcmp(cutype, "RANGE") == 0)
1689  return UT_RANGE;
1690 
1691  return UT_LAST;
1692 }
1693 
1696  int nband; /* source raster's band index, 0-based */
1698 
1701 };
1702 
1705  int numband; /* number of bandargs */
1706  rtpg_union_band_arg bandarg;
1707 };
1708 
1709 static void rtpg_union_arg_destroy(rtpg_union_arg arg) {
1710  int i = 0;
1711  int j = 0;
1712  int k = 0;
1713 
1714  if (arg->bandarg != NULL) {
1715  for (i = 0; i < arg->numband; i++) {
1716  if (!arg->bandarg[i].numraster)
1717  continue;
1718 
1719  for (j = 0; j < arg->bandarg[i].numraster; j++) {
1720  if (arg->bandarg[i].raster[j] == NULL)
1721  continue;
1722 
1723  for (k = rt_raster_get_num_bands(arg->bandarg[i].raster[j]) - 1; k >= 0; k--)
1725  rt_raster_destroy(arg->bandarg[i].raster[j]);
1726  }
1727 
1728  pfree(arg->bandarg[i].raster);
1729  }
1730 
1731  pfree(arg->bandarg);
1732  }
1733 
1734  pfree(arg);
1735 }
1736 
1738  rt_iterator_arg arg, void *userarg,
1739  double *value, int *nodata
1740 ) {
1741  rtpg_union_type *utype = (rtpg_union_type *) userarg;
1742 
1743  if (arg == NULL)
1744  return 0;
1745 
1746  if (
1747  arg->rasters != 2 ||
1748  arg->rows != 1 ||
1749  arg->columns != 1
1750  ) {
1751  elog(ERROR, "rtpg_union_callback: Invalid arguments passed to callback");
1752  return 0;
1753  }
1754 
1755  *value = 0;
1756  *nodata = 0;
1757 
1758  /* handle NODATA situations except for COUNT, which is a special case */
1759  if (*utype != UT_COUNT) {
1760  /* both NODATA */
1761  if (arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1762  *nodata = 1;
1763  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1764  return 1;
1765  }
1766  /* second NODATA */
1767  else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1768  *value = arg->values[0][0][0];
1769  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1770  return 1;
1771  }
1772  /* first NODATA */
1773  else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0]) {
1774  *value = arg->values[1][0][0];
1775  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1776  return 1;
1777  }
1778  }
1779 
1780  switch (*utype) {
1781  case UT_FIRST:
1782  *value = arg->values[0][0][0];
1783  break;
1784  case UT_MIN:
1785  if (arg->values[0][0][0] < arg->values[1][0][0])
1786  *value = arg->values[0][0][0];
1787  else
1788  *value = arg->values[1][0][0];
1789  break;
1790  case UT_MAX:
1791  if (arg->values[0][0][0] > arg->values[1][0][0])
1792  *value = arg->values[0][0][0];
1793  else
1794  *value = arg->values[1][0][0];
1795  break;
1796  case UT_COUNT:
1797  /* both NODATA */
1798  if (arg->nodata[0][0][0] && arg->nodata[1][0][0])
1799  *value = 0;
1800  /* second NODATA */
1801  else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0])
1802  *value = arg->values[0][0][0];
1803  /* first NODATA */
1804  else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0])
1805  *value = 1;
1806  /* has value, increment */
1807  else
1808  *value = arg->values[0][0][0] + 1;
1809  break;
1810  case UT_SUM:
1811  *value = arg->values[0][0][0] + arg->values[1][0][0];
1812  break;
1813  case UT_MEAN:
1814  case UT_RANGE:
1815  break;
1816  case UT_LAST:
1817  default:
1818  *value = arg->values[1][0][0];
1819  break;
1820  }
1821 
1822  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1823 
1824 
1825  return 1;
1826 }
1827 
1829  rt_iterator_arg arg, void *userarg,
1830  double *value, int *nodata
1831 ) {
1832  if (arg == NULL)
1833  return 0;
1834 
1835  if (
1836  arg->rasters != 2 ||
1837  arg->rows != 1 ||
1838  arg->columns != 1
1839  ) {
1840  elog(ERROR, "rtpg_union_mean_callback: Invalid arguments passed to callback");
1841  return 0;
1842  }
1843 
1844  *value = 0;
1845  *nodata = 1;
1846 
1847  POSTGIS_RT_DEBUGF(4, "rast0: %f %d", arg->values[0][0][0], arg->nodata[0][0][0]);
1848  POSTGIS_RT_DEBUGF(4, "rast1: %f %d", arg->values[1][0][0], arg->nodata[1][0][0]);
1849 
1850  if (
1851  !arg->nodata[0][0][0] &&
1852  FLT_NEQ(arg->values[0][0][0], 0) &&
1853  !arg->nodata[1][0][0]
1854  ) {
1855  *value = arg->values[1][0][0] / arg->values[0][0][0];
1856  *nodata = 0;
1857  }
1858 
1859  POSTGIS_RT_DEBUGF(4, "value, nodata = (%f, %d)", *value, *nodata);
1860 
1861  return 1;
1862 }
1863 
1865  rt_iterator_arg arg, void *userarg,
1866  double *value, int *nodata
1867 ) {
1868  if (arg == NULL)
1869  return 0;
1870 
1871  if (
1872  arg->rasters != 2 ||
1873  arg->rows != 1 ||
1874  arg->columns != 1
1875  ) {
1876  elog(ERROR, "rtpg_union_range_callback: Invalid arguments passed to callback");
1877  return 0;
1878  }
1879 
1880  *value = 0;
1881  *nodata = 1;
1882 
1883  POSTGIS_RT_DEBUGF(4, "rast0: %f %d", arg->values[0][0][0], arg->nodata[0][0][0]);
1884  POSTGIS_RT_DEBUGF(4, "rast1: %f %d", arg->values[1][0][0], arg->nodata[1][0][0]);
1885 
1886  if (
1887  !arg->nodata[0][0][0] &&
1888  !arg->nodata[1][0][0]
1889  ) {
1890  *value = arg->values[1][0][0] - arg->values[0][0][0];
1891  *nodata = 0;
1892  }
1893 
1894  POSTGIS_RT_DEBUGF(4, "value, nodata = (%f, %d)", *value, *nodata);
1895 
1896  return 1;
1897 }
1898 
1899 /* called for ST_Union(raster, unionarg[]) */
1900 static int rtpg_union_unionarg_process(rtpg_union_arg arg, ArrayType *array) {
1901  Oid etype;
1902  Datum *e;
1903  bool *nulls;
1904  int16 typlen;
1905  bool typbyval;
1906  char typalign;
1907  int n = 0;
1908 
1909  HeapTupleHeader tup;
1910  bool isnull;
1911  Datum tupv;
1912 
1913  int i;
1914  int nband = 1;
1915  char *utypename = NULL;
1916  rtpg_union_type utype = UT_LAST;
1917 
1918  etype = ARR_ELEMTYPE(array);
1919  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1920 
1921  deconstruct_array(
1922  array,
1923  etype,
1924  typlen, typbyval, typalign,
1925  &e, &nulls, &n
1926  );
1927 
1928  if (!n) {
1929  elog(ERROR, "rtpg_union_unionarg_process: Invalid argument for unionarg");
1930  return 0;
1931  }
1932 
1933  /* prep arg */
1934  arg->numband = n;
1935  arg->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * arg->numband);
1936  if (arg->bandarg == NULL) {
1937  elog(ERROR, "rtpg_union_unionarg_process: Could not allocate memory for band information");
1938  return 0;
1939  }
1940 
1941  /* process each element */
1942  for (i = 0; i < n; i++) {
1943  if (nulls[i]) {
1944  arg->numband--;
1945  continue;
1946  }
1947 
1948  POSTGIS_RT_DEBUGF(4, "Processing unionarg at index %d", i);
1949 
1950  /* each element is a tuple */
1951  tup = (HeapTupleHeader) DatumGetPointer(e[i]);
1952  if (NULL == tup) {
1953  elog(ERROR, "rtpg_union_unionarg_process: Invalid argument for unionarg");
1954  return 0;
1955  }
1956 
1957  /* first element, bandnum */
1958  tupv = GetAttributeByName(tup, "nband", &isnull);
1959  if (isnull) {
1960  nband = i + 1;
1961  elog(NOTICE, "First argument (nband) of unionarg is NULL. Assuming nband = %d", nband);
1962  }
1963  else
1964  nband = DatumGetInt32(tupv);
1965 
1966  if (nband < 1) {
1967  elog(ERROR, "rtpg_union_unionarg_process: Band number must be greater than zero (1-based)");
1968  return 0;
1969  }
1970 
1971  /* second element, uniontype */
1972  tupv = GetAttributeByName(tup, "uniontype", &isnull);
1973  if (isnull) {
1974  elog(NOTICE, "Second argument (uniontype) of unionarg is NULL. Assuming uniontype = LAST");
1975  utype = UT_LAST;
1976  }
1977  else {
1978  utypename = text_to_cstring((text *) DatumGetPointer(tupv));
1979  utype = rtpg_uniontype_index_from_name(rtpg_strtoupper(utypename));
1980  }
1981 
1982  arg->bandarg[i].uniontype = utype;
1983  arg->bandarg[i].nband = nband - 1;
1984  arg->bandarg[i].raster = NULL;
1985 
1986  if (
1987  utype != UT_MEAN &&
1988  utype != UT_RANGE
1989  ) {
1990  arg->bandarg[i].numraster = 1;
1991  }
1992  else
1993  arg->bandarg[i].numraster = 2;
1994  }
1995 
1996  if (arg->numband < n) {
1997  arg->bandarg = repalloc(arg->bandarg, sizeof(struct rtpg_union_band_arg_t) * arg->numband);
1998  if (arg->bandarg == NULL) {
1999  elog(ERROR, "rtpg_union_unionarg_process: Could not reallocate memory for band information");
2000  return 0;
2001  }
2002  }
2003 
2004  return 1;
2005 }
2006 
2007 /* called for ST_Union(raster) */
2008 static int rtpg_union_noarg(rtpg_union_arg arg, rt_raster raster) {
2009  int numbands;
2010  int i;
2011 
2012  if (rt_raster_is_empty(raster))
2013  return 1;
2014 
2015  numbands = rt_raster_get_num_bands(raster);
2016  if (numbands <= arg->numband)
2017  return 1;
2018 
2019  /* more bands to process */
2020  POSTGIS_RT_DEBUG(4, "input raster has more bands, adding more bandargs");
2021  if (arg->numband)
2022  arg->bandarg = repalloc(arg->bandarg, sizeof(struct rtpg_union_band_arg_t) * numbands);
2023  else
2024  arg->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * numbands);
2025  if (arg->bandarg == NULL) {
2026  elog(ERROR, "rtpg_union_noarg: Could not reallocate memory for band information");
2027  return 0;
2028  }
2029 
2030  i = arg->numband;
2031  arg->numband = numbands;
2032  for (; i < arg->numband; i++) {
2033  POSTGIS_RT_DEBUGF(4, "Adding bandarg for band at index %d", i);
2034  arg->bandarg[i].uniontype = UT_LAST;
2035  arg->bandarg[i].nband = i;
2036  arg->bandarg[i].numraster = 1;
2037 
2038  arg->bandarg[i].raster = (rt_raster *) palloc(sizeof(rt_raster) * arg->bandarg[i].numraster);
2039  if (arg->bandarg[i].raster == NULL) {
2040  elog(ERROR, "rtpg_union_noarg: Could not allocate memory for working rasters");
2041  return 0;
2042  }
2043  memset(arg->bandarg[i].raster, 0, sizeof(rt_raster) * arg->bandarg[i].numraster);
2044 
2045  /* add new working rt_raster but only if working raster already exists */
2046  if (!rt_raster_is_empty(arg->bandarg[0].raster[0])) {
2047  arg->bandarg[i].raster[0] = rt_raster_clone(arg->bandarg[0].raster[0], 0); /* shallow clone */
2048  if (arg->bandarg[i].raster[0] == NULL) {
2049  elog(ERROR, "rtpg_union_noarg: Could not create working raster");
2050  return 0;
2051  }
2052  }
2053  }
2054 
2055  return 1;
2056 }
2057 
2058 /* UNION aggregate transition function */
2060 Datum RASTER_union_transfn(PG_FUNCTION_ARGS)
2061 {
2062  MemoryContext aggcontext;
2063  MemoryContext oldcontext;
2064  rtpg_union_arg iwr = NULL;
2065  int skiparg = 0;
2066 
2067  rt_pgraster *pgraster = NULL;
2068  rt_raster raster = NULL;
2069  rt_raster _raster = NULL;
2070  rt_band _band = NULL;
2071  int nband = 1;
2072  int noerr = 1;
2073  int isempty[2] = {0};
2074  int hasband[2] = {0};
2075  int nargs = 0;
2076  double _offset[4] = {0.};
2077  int nbnodata = 0; /* 1 if adding bands */
2078 
2079  int i = 0;
2080  int j = 0;
2081  int k = 0;
2082 
2083  rt_iterator itrset;
2084  char *utypename = NULL;
2085  rtpg_union_type utype = UT_LAST;
2086  rt_pixtype pixtype = PT_END;
2087  int hasnodata = 1;
2088  double nodataval = 0;
2089 
2090  rt_raster iraster = NULL;
2091  rt_band iband = NULL;
2092  int reuserast = 0;
2093  int y = 0;
2094  uint16_t _dim[2] = {0};
2095  void *vals = NULL;
2096  uint16_t nvals = 0;
2097 
2098  POSTGIS_RT_DEBUG(3, "Starting...");
2099 
2100  /* cannot be called directly as this is exclusive aggregate function */
2101  if (!AggCheckCallContext(fcinfo, &aggcontext)) {
2102  elog(ERROR, "RASTER_union_transfn: Cannot be called in a non-aggregate context");
2103  PG_RETURN_NULL();
2104  }
2105 
2106  /* switch to aggcontext */
2107  oldcontext = MemoryContextSwitchTo(aggcontext);
2108 
2109  if (PG_ARGISNULL(0)) {
2110  POSTGIS_RT_DEBUG(3, "Creating state variable");
2111  /* allocate container in aggcontext */
2112  iwr = palloc(sizeof(struct rtpg_union_arg_t));
2113  if (iwr == NULL) {
2114  MemoryContextSwitchTo(oldcontext);
2115  elog(ERROR, "RASTER_union_transfn: Could not allocate memory for state variable");
2116  PG_RETURN_NULL();
2117  }
2118 
2119  iwr->numband = 0;
2120  iwr->bandarg = NULL;
2121 
2122  skiparg = 0;
2123  }
2124  else {
2125  POSTGIS_RT_DEBUG(3, "State variable already exists");
2126  iwr = (rtpg_union_arg) PG_GETARG_POINTER(0);
2127  skiparg = 1;
2128  }
2129 
2130  /* raster arg is NOT NULL */
2131  if (!PG_ARGISNULL(1)) {
2132  /* deserialize raster */
2133  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2134 
2135  /* Get raster object */
2136  raster = rt_raster_deserialize(pgraster, FALSE);
2137  if (raster == NULL) {
2138 
2140  PG_FREE_IF_COPY(pgraster, 1);
2141 
2142  MemoryContextSwitchTo(oldcontext);
2143  elog(ERROR, "RASTER_union_transfn: Could not deserialize raster");
2144  PG_RETURN_NULL();
2145  }
2146  }
2147 
2148  /* process additional args if needed */
2149  nargs = PG_NARGS();
2150  POSTGIS_RT_DEBUGF(4, "nargs = %d", nargs);
2151  if (nargs > 2) {
2152  POSTGIS_RT_DEBUG(4, "processing additional arguments");
2153 
2154  /* if more than 2 arguments, determine the type of argument 3 */
2155  /* band number, UNION type or unionarg */
2156  if (!PG_ARGISNULL(2)) {
2157  Oid calltype = get_fn_expr_argtype(fcinfo->flinfo, 2);
2158 
2159  switch (calltype) {
2160  /* UNION type */
2161  case TEXTOID: {
2162  int idx = 0;
2163  int numband = 0;
2164 
2165  POSTGIS_RT_DEBUG(4, "Processing arg 3 as UNION type");
2166  nbnodata = 1;
2167 
2168  utypename = text_to_cstring(PG_GETARG_TEXT_P(2));
2169  utype = rtpg_uniontype_index_from_name(rtpg_strtoupper(utypename));
2170  POSTGIS_RT_DEBUGF(4, "Union type: %s", utypename);
2171 
2172  POSTGIS_RT_DEBUGF(4, "iwr->numband: %d", iwr->numband);
2173  /* see if we need to append new bands */
2174  if (raster) {
2175  idx = iwr->numband;
2176  numband = rt_raster_get_num_bands(raster);
2177  POSTGIS_RT_DEBUGF(4, "numband: %d", numband);
2178 
2179  /* only worry about appended bands */
2180  if (numband > iwr->numband)
2181  iwr->numband = numband;
2182  }
2183 
2184  if (!iwr->numband)
2185  iwr->numband = 1;
2186  POSTGIS_RT_DEBUGF(4, "iwr->numband: %d", iwr->numband);
2187  POSTGIS_RT_DEBUGF(4, "numband, idx: %d, %d", numband, idx);
2188 
2189  /* bandarg set. only possible after the first call to function */
2190  if (iwr->bandarg) {
2191  /* only reallocate if new bands need to be added */
2192  if (numband > idx) {
2193  POSTGIS_RT_DEBUG(4, "Reallocating iwr->bandarg");
2194  iwr->bandarg = repalloc(iwr->bandarg, sizeof(struct rtpg_union_band_arg_t) * iwr->numband);
2195  }
2196  /* prevent initial values step happening below */
2197  else
2198  idx = iwr->numband;
2199  }
2200  /* bandarg not set, first call to function */
2201  else {
2202  POSTGIS_RT_DEBUG(4, "Allocating iwr->bandarg");
2203  iwr->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * iwr->numband);
2204  }
2205  if (iwr->bandarg == NULL) {
2206 
2208  if (raster != NULL) {
2209  rt_raster_destroy(raster);
2210  PG_FREE_IF_COPY(pgraster, 1);
2211  }
2212 
2213  MemoryContextSwitchTo(oldcontext);
2214  elog(ERROR, "RASTER_union_transfn: Could not allocate memory for band information");
2215  PG_RETURN_NULL();
2216  }
2217 
2218  /* set initial values for bands that are "new" */
2219  for (i = idx; i < iwr->numband; i++) {
2220  iwr->bandarg[i].uniontype = utype;
2221  iwr->bandarg[i].nband = i;
2222 
2223  if (
2224  utype == UT_MEAN ||
2225  utype == UT_RANGE
2226  ) {
2227  iwr->bandarg[i].numraster = 2;
2228  }
2229  else
2230  iwr->bandarg[i].numraster = 1;
2231  iwr->bandarg[i].raster = NULL;
2232  }
2233 
2234  break;
2235  }
2236  /* band number */
2237  case INT2OID:
2238  case INT4OID:
2239  if (skiparg)
2240  break;
2241 
2242  POSTGIS_RT_DEBUG(4, "Processing arg 3 as band number");
2243  nband = PG_GETARG_INT32(2);
2244  if (nband < 1) {
2245 
2247  if (raster != NULL) {
2248  rt_raster_destroy(raster);
2249  PG_FREE_IF_COPY(pgraster, 1);
2250  }
2251 
2252  MemoryContextSwitchTo(oldcontext);
2253  elog(ERROR, "RASTER_union_transfn: Band number must be greater than zero (1-based)");
2254  PG_RETURN_NULL();
2255  }
2256 
2257  iwr->numband = 1;
2258  iwr->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * iwr->numband);
2259  if (iwr->bandarg == NULL) {
2260 
2262  if (raster != NULL) {
2263  rt_raster_destroy(raster);
2264  PG_FREE_IF_COPY(pgraster, 1);
2265  }
2266 
2267  MemoryContextSwitchTo(oldcontext);
2268  elog(ERROR, "RASTER_union_transfn: Could not allocate memory for band information");
2269  PG_RETURN_NULL();
2270  }
2271 
2272  iwr->bandarg[0].uniontype = UT_LAST;
2273  iwr->bandarg[0].nband = nband - 1;
2274 
2275  iwr->bandarg[0].numraster = 1;
2276  iwr->bandarg[0].raster = NULL;
2277  break;
2278  /* only other type allowed is unionarg */
2279  default:
2280  if (skiparg)
2281  break;
2282 
2283  POSTGIS_RT_DEBUG(4, "Processing arg 3 as unionarg");
2284  if (!rtpg_union_unionarg_process(iwr, PG_GETARG_ARRAYTYPE_P(2))) {
2285 
2287  if (raster != NULL) {
2288  rt_raster_destroy(raster);
2289  PG_FREE_IF_COPY(pgraster, 1);
2290  }
2291 
2292  MemoryContextSwitchTo(oldcontext);
2293  elog(ERROR, "RASTER_union_transfn: Could not process unionarg");
2294  PG_RETURN_NULL();
2295  }
2296 
2297  break;
2298  }
2299  }
2300 
2301  /* UNION type */
2302  if (nargs > 3 && !PG_ARGISNULL(3)) {
2303  utypename = text_to_cstring(PG_GETARG_TEXT_P(3));
2304  utype = rtpg_uniontype_index_from_name(rtpg_strtoupper(utypename));
2305  iwr->bandarg[0].uniontype = utype;
2306  POSTGIS_RT_DEBUGF(4, "Union type: %s", utypename);
2307 
2308  if (
2309  utype == UT_MEAN ||
2310  utype == UT_RANGE
2311  ) {
2312  iwr->bandarg[0].numraster = 2;
2313  }
2314  }
2315 
2316  /* allocate space for pointers to rt_raster */
2317  for (i = 0; i < iwr->numband; i++) {
2318  POSTGIS_RT_DEBUGF(4, "iwr->bandarg[%d].raster @ %p", i, iwr->bandarg[i].raster);
2319 
2320  /* no need to allocate */
2321  if (iwr->bandarg[i].raster != NULL)
2322  continue;
2323 
2324  POSTGIS_RT_DEBUGF(4, "Allocating space for working rasters of band %d", i);
2325 
2326  iwr->bandarg[i].raster = (rt_raster *) palloc(sizeof(rt_raster) * iwr->bandarg[i].numraster);
2327  if (iwr->bandarg[i].raster == NULL) {
2328 
2330  if (raster != NULL) {
2331  rt_raster_destroy(raster);
2332  PG_FREE_IF_COPY(pgraster, 1);
2333  }
2334 
2335  MemoryContextSwitchTo(oldcontext);
2336  elog(ERROR, "RASTER_union_transfn: Could not allocate memory for working raster(s)");
2337  PG_RETURN_NULL();
2338  }
2339 
2340  memset(iwr->bandarg[i].raster, 0, sizeof(rt_raster) * iwr->bandarg[i].numraster);
2341 
2342  /* add new working rt_raster but only if working raster already exists */
2343  if (i > 0 && !rt_raster_is_empty(iwr->bandarg[0].raster[0])) {
2344  for (j = 0; j < iwr->bandarg[i].numraster; j++) {
2345  iwr->bandarg[i].raster[j] = rt_raster_clone(iwr->bandarg[0].raster[0], 0); /* shallow clone */
2346  if (iwr->bandarg[i].raster[j] == 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 create working raster");
2356  PG_RETURN_NULL();
2357  }
2358  }
2359  }
2360  }
2361  }
2362  /* only raster, no additional args */
2363  /* only do this if raster isn't empty */
2364  else {
2365  POSTGIS_RT_DEBUG(4, "no additional args, checking input raster");
2366  nbnodata = 1;
2367  if (!rtpg_union_noarg(iwr, raster)) {
2368 
2370  if (raster != NULL) {
2371  rt_raster_destroy(raster);
2372  PG_FREE_IF_COPY(pgraster, 1);
2373  }
2374 
2375  MemoryContextSwitchTo(oldcontext);
2376  elog(ERROR, "RASTER_union_transfn: Could not check and balance number of bands");
2377  PG_RETURN_NULL();
2378  }
2379  }
2380 
2381  /* init itrset */
2382  itrset = palloc(sizeof(struct rt_iterator_t) * 2);
2383  if (itrset == NULL) {
2384 
2386  if (raster != NULL) {
2387  rt_raster_destroy(raster);
2388  PG_FREE_IF_COPY(pgraster, 1);
2389  }
2390 
2391  MemoryContextSwitchTo(oldcontext);
2392  elog(ERROR, "RASTER_union_transfn: Could not allocate memory for iterator arguments");
2393  PG_RETURN_NULL();
2394  }
2395 
2396  /* by band to UNION */
2397  for (i = 0; i < iwr->numband; i++) {
2398 
2399  /* by raster */
2400  for (j = 0; j < iwr->bandarg[i].numraster; j++) {
2401  reuserast = 0;
2402 
2403  /* type of union */
2404  utype = iwr->bandarg[i].uniontype;
2405 
2406  /* raster flags */
2407  isempty[0] = rt_raster_is_empty(iwr->bandarg[i].raster[j]);
2408  isempty[1] = rt_raster_is_empty(raster);
2409 
2410  if (!isempty[0])
2411  hasband[0] = rt_raster_has_band(iwr->bandarg[i].raster[j], 0);
2412  if (!isempty[1])
2413  hasband[1] = rt_raster_has_band(raster, iwr->bandarg[i].nband);
2414 
2415  /* determine pixtype, hasnodata and nodataval */
2416  _band = NULL;
2417  if (!isempty[0] && hasband[0])
2418  _band = rt_raster_get_band(iwr->bandarg[i].raster[j], 0);
2419  else if (!isempty[1] && hasband[1])
2420  _band = rt_raster_get_band(raster, iwr->bandarg[i].nband);
2421  else {
2422  pixtype = PT_64BF;
2423  hasnodata = 1;
2424  nodataval = rt_pixtype_get_min_value(pixtype);
2425  }
2426  if (_band != NULL) {
2427  pixtype = rt_band_get_pixtype(_band);
2428  hasnodata = 1;
2429  if (rt_band_get_hasnodata_flag(_band))
2430  rt_band_get_nodata(_band, &nodataval);
2431  else
2432  nodataval = rt_band_get_min_value(_band);
2433  }
2434 
2435  /* UT_MEAN and UT_RANGE require two passes */
2436  /* UT_MEAN: first for UT_COUNT and second for UT_SUM */
2437  if (iwr->bandarg[i].uniontype == UT_MEAN) {
2438  /* first pass, UT_COUNT */
2439  if (j < 1)
2440  utype = UT_COUNT;
2441  else
2442  utype = UT_SUM;
2443  }
2444  /* UT_RANGE: first for UT_MIN and second for UT_MAX */
2445  else if (iwr->bandarg[i].uniontype == UT_RANGE) {
2446  /* first pass, UT_MIN */
2447  if (j < 1)
2448  utype = UT_MIN;
2449  else
2450  utype = UT_MAX;
2451  }
2452 
2453  /* force band settings for UT_COUNT */
2454  if (utype == UT_COUNT) {
2455  pixtype = PT_32BUI;
2456  hasnodata = 0;
2457  nodataval = 0;
2458  }
2459 
2460  POSTGIS_RT_DEBUGF(4, "(pixtype, hasnodata, nodataval) = (%s, %d, %f)", rt_pixtype_name(pixtype), hasnodata, nodataval);
2461 
2462  /* set itrset */
2463  itrset[0].raster = iwr->bandarg[i].raster[j];
2464  itrset[0].nband = 0;
2465  itrset[1].raster = raster;
2466  itrset[1].nband = iwr->bandarg[i].nband;
2467 
2468  /* allow use NODATA to replace missing bands */
2469  if (nbnodata) {
2470  itrset[0].nbnodata = 1;
2471  itrset[1].nbnodata = 1;
2472  }
2473  /* missing bands are not ignored */
2474  else {
2475  itrset[0].nbnodata = 0;
2476  itrset[1].nbnodata = 0;
2477  }
2478 
2479  /* if rasters AND bands are present, use copy approach */
2480  if (!isempty[0] && !isempty[1] && hasband[0] && hasband[1]) {
2481  POSTGIS_RT_DEBUG(3, "using line method");
2482 
2483  /* generate empty out raster */
2485  iwr->bandarg[i].raster[j], raster,
2486  ET_UNION,
2487  &iraster, _offset
2488  ) != ES_NONE) {
2489 
2490  pfree(itrset);
2492  if (raster != NULL) {
2493  rt_raster_destroy(raster);
2494  PG_FREE_IF_COPY(pgraster, 1);
2495  }
2496 
2497  MemoryContextSwitchTo(oldcontext);
2498  elog(ERROR, "RASTER_union_transfn: Could not create internal raster");
2499  PG_RETURN_NULL();
2500  }
2501  POSTGIS_RT_DEBUGF(4, "_offset = %f, %f, %f, %f",
2502  _offset[0], _offset[1], _offset[2], _offset[3]);
2503 
2504  /* rasters are spatially the same? */
2505  if (
2506  rt_raster_get_width(iwr->bandarg[i].raster[j]) == rt_raster_get_width(iraster) &&
2508  ) {
2509  double igt[6] = {0};
2510  double gt[6] = {0};
2511 
2513  rt_raster_get_geotransform_matrix(iraster, igt);
2514 
2515  reuserast = rt_util_same_geotransform_matrix(gt, igt);
2516  }
2517 
2518  /* use internal raster */
2519  if (!reuserast) {
2520  /* create band of same type */
2522  iraster,
2523  pixtype,
2524  nodataval,
2525  hasnodata, nodataval,
2526  0
2527  ) == -1) {
2528 
2529  pfree(itrset);
2531  rt_raster_destroy(iraster);
2532  if (raster != NULL) {
2533  rt_raster_destroy(raster);
2534  PG_FREE_IF_COPY(pgraster, 1);
2535  }
2536 
2537  MemoryContextSwitchTo(oldcontext);
2538  elog(ERROR, "RASTER_union_transfn: Could not add new band to internal raster");
2539  PG_RETURN_NULL();
2540  }
2541  iband = rt_raster_get_band(iraster, 0);
2542 
2543  /* copy working raster to output raster */
2544  _dim[0] = rt_raster_get_width(iwr->bandarg[i].raster[j]);
2545  _dim[1] = rt_raster_get_height(iwr->bandarg[i].raster[j]);
2546  for (y = 0; y < _dim[1]; y++) {
2547  POSTGIS_RT_DEBUGF(4, "Getting pixel line of working raster at (x, y, length) = (0, %d, %d)", y, _dim[0]);
2549  _band,
2550  0, y,
2551  _dim[0],
2552  &vals, &nvals
2553  ) != ES_NONE) {
2554 
2555  pfree(itrset);
2557  rt_band_destroy(iband);
2558  rt_raster_destroy(iraster);
2559  if (raster != NULL) {
2560  rt_raster_destroy(raster);
2561  PG_FREE_IF_COPY(pgraster, 1);
2562  }
2563 
2564  MemoryContextSwitchTo(oldcontext);
2565  elog(ERROR, "RASTER_union_transfn: Could not get pixel line from band of working raster");
2566  PG_RETURN_NULL();
2567  }
2568 
2569  POSTGIS_RT_DEBUGF(4, "Setting pixel line at (x, y, length) = (%d, %d, %d)", (int) _offset[0], (int) _offset[1] + y, nvals);
2571  iband,
2572  (int) _offset[0], (int) _offset[1] + y,
2573  vals, nvals
2574  ) != ES_NONE) {
2575 
2576  pfree(itrset);
2578  rt_band_destroy(iband);
2579  rt_raster_destroy(iraster);
2580  if (raster != NULL) {
2581  rt_raster_destroy(raster);
2582  PG_FREE_IF_COPY(pgraster, 1);
2583  }
2584 
2585  MemoryContextSwitchTo(oldcontext);
2586  elog(ERROR, "RASTER_union_transfn: Could not set pixel line to band of internal raster");
2587  PG_RETURN_NULL();
2588  }
2589  }
2590  }
2591  else {
2592  rt_raster_destroy(iraster);
2593  iraster = iwr->bandarg[i].raster[j];
2594  iband = rt_raster_get_band(iraster, 0);
2595  }
2596 
2597  /* run iterator for extent of input raster */
2598  noerr = rt_raster_iterator(
2599  itrset, 2,
2600  ET_LAST, NULL,
2601  pixtype,
2602  hasnodata, nodataval,
2603  0, 0,
2604  NULL,
2605  &utype,
2607  &_raster
2608  );
2609  if (noerr != ES_NONE) {
2610 
2611  pfree(itrset);
2613  if (!reuserast) {
2614  rt_band_destroy(iband);
2615  rt_raster_destroy(iraster);
2616  }
2617  if (raster != NULL) {
2618  rt_raster_destroy(raster);
2619  PG_FREE_IF_COPY(pgraster, 1);
2620  }
2621 
2622  MemoryContextSwitchTo(oldcontext);
2623  elog(ERROR, "RASTER_union_transfn: Could not run raster iterator function");
2624  PG_RETURN_NULL();
2625  }
2626 
2627  /* with iterator raster, copy data to output raster */
2628  _band = rt_raster_get_band(_raster, 0);
2629  _dim[0] = rt_raster_get_width(_raster);
2630  _dim[1] = rt_raster_get_height(_raster);
2631  for (y = 0; y < _dim[1]; y++) {
2632  POSTGIS_RT_DEBUGF(4, "Getting pixel line of iterator raster at (x, y, length) = (0, %d, %d)", y, _dim[0]);
2634  _band,
2635  0, y,
2636  _dim[0],
2637  &vals, &nvals
2638  ) != ES_NONE) {
2639 
2640  pfree(itrset);
2642  if (!reuserast) {
2643  rt_band_destroy(iband);
2644  rt_raster_destroy(iraster);
2645  }
2646  if (raster != NULL) {
2647  rt_raster_destroy(raster);
2648  PG_FREE_IF_COPY(pgraster, 1);
2649  }
2650 
2651  MemoryContextSwitchTo(oldcontext);
2652  elog(ERROR, "RASTER_union_transfn: Could not get pixel line from band of working raster");
2653  PG_RETURN_NULL();
2654  }
2655 
2656  POSTGIS_RT_DEBUGF(4, "Setting pixel line at (x, y, length) = (%d, %d, %d)", (int) _offset[2], (int) _offset[3] + y, nvals);
2658  iband,
2659  (int) _offset[2], (int) _offset[3] + y,
2660  vals, nvals
2661  ) != ES_NONE) {
2662 
2663  pfree(itrset);
2665  if (!reuserast) {
2666  rt_band_destroy(iband);
2667  rt_raster_destroy(iraster);
2668  }
2669  if (raster != NULL) {
2670  rt_raster_destroy(raster);
2671  PG_FREE_IF_COPY(pgraster, 1);
2672  }
2673 
2674  MemoryContextSwitchTo(oldcontext);
2675  elog(ERROR, "RASTER_union_transfn: Could not set pixel line to band of internal raster");
2676  PG_RETURN_NULL();
2677  }
2678  }
2679 
2680  /* free _raster */
2681  rt_band_destroy(_band);
2682  rt_raster_destroy(_raster);
2683 
2684  /* replace working raster with output raster */
2685  _raster = iraster;
2686  }
2687  else {
2688  POSTGIS_RT_DEBUG(3, "using pixel method");
2689 
2690  /* pass everything to iterator */
2691  noerr = rt_raster_iterator(
2692  itrset, 2,
2693  ET_UNION, NULL,
2694  pixtype,
2695  hasnodata, nodataval,
2696  0, 0,
2697  NULL,
2698  &utype,
2700  &_raster
2701  );
2702 
2703  if (noerr != ES_NONE) {
2704 
2705  pfree(itrset);
2707  if (raster != NULL) {
2708  rt_raster_destroy(raster);
2709  PG_FREE_IF_COPY(pgraster, 1);
2710  }
2711 
2712  MemoryContextSwitchTo(oldcontext);
2713  elog(ERROR, "RASTER_union_transfn: Could not run raster iterator function");
2714  PG_RETURN_NULL();
2715  }
2716  }
2717 
2718  /* replace working raster */
2719  if (iwr->bandarg[i].raster[j] != NULL && !reuserast) {
2720  for (k = rt_raster_get_num_bands(iwr->bandarg[i].raster[j]) - 1; k >= 0; k--)
2722  rt_raster_destroy(iwr->bandarg[i].raster[j]);
2723  }
2724  iwr->bandarg[i].raster[j] = _raster;
2725  }
2726 
2727  }
2728 
2729  pfree(itrset);
2730  if (raster != NULL) {
2731  rt_raster_destroy(raster);
2732  PG_FREE_IF_COPY(pgraster, 1);
2733  }
2734 
2735  /* switch back to local context */
2736  MemoryContextSwitchTo(oldcontext);
2737 
2738  POSTGIS_RT_DEBUG(3, "Finished");
2739 
2740  PG_RETURN_POINTER(iwr);
2741 }
2742 
2743 /* UNION aggregate final function */
2745 Datum RASTER_union_finalfn(PG_FUNCTION_ARGS)
2746 {
2747  rtpg_union_arg iwr;
2748  rt_raster _rtn = NULL;
2749  rt_raster _raster = NULL;
2750  rt_pgraster *pgraster = NULL;
2751 
2752  int i = 0;
2753  int j = 0;
2754  rt_iterator itrset = NULL;
2755  rt_band _band = NULL;
2756  int noerr = 1;
2757  int status = 0;
2758  rt_pixtype pixtype = PT_END;
2759  int hasnodata = 0;
2760  double nodataval = 0;
2761 
2762  POSTGIS_RT_DEBUG(3, "Starting...");
2763 
2764  /* cannot be called directly as this is exclusive aggregate function */
2765  if (!AggCheckCallContext(fcinfo, NULL)) {
2766  elog(ERROR, "RASTER_union_finalfn: Cannot be called in a non-aggregate context");
2767  PG_RETURN_NULL();
2768  }
2769 
2770  /* NULL, return null */
2771  if (PG_ARGISNULL(0))
2772  PG_RETURN_NULL();
2773 
2774  iwr = (rtpg_union_arg) PG_GETARG_POINTER(0);
2775 
2776  /* init itrset */
2777  itrset = palloc(sizeof(struct rt_iterator_t) * 2);
2778  if (itrset == NULL) {
2780  elog(ERROR, "RASTER_union_finalfn: Could not allocate memory for iterator arguments");
2781  PG_RETURN_NULL();
2782  }
2783 
2784  for (i = 0; i < iwr->numband; i++) {
2785  if (
2786  iwr->bandarg[i].uniontype == UT_MEAN ||
2787  iwr->bandarg[i].uniontype == UT_RANGE
2788  ) {
2789  /* raster containing the SUM or MAX is at index 1 */
2790  _band = rt_raster_get_band(iwr->bandarg[i].raster[1], 0);
2791 
2792  pixtype = rt_band_get_pixtype(_band);
2793  hasnodata = rt_band_get_hasnodata_flag(_band);
2794  if (hasnodata)
2795  rt_band_get_nodata(_band, &nodataval);
2796  POSTGIS_RT_DEBUGF(4, "(pixtype, hasnodata, nodataval) = (%s, %d, %f)", rt_pixtype_name(pixtype), hasnodata, nodataval);
2797 
2798  itrset[0].raster = iwr->bandarg[i].raster[0];
2799  itrset[0].nband = 0;
2800  itrset[1].raster = iwr->bandarg[i].raster[1];
2801  itrset[1].nband = 0;
2802 
2803  /* pass everything to iterator */
2804  if (iwr->bandarg[i].uniontype == UT_MEAN) {
2805  noerr = rt_raster_iterator(
2806  itrset, 2,
2807  ET_UNION, NULL,
2808  pixtype,
2809  hasnodata, nodataval,
2810  0, 0,
2811  NULL,
2812  NULL,
2814  &_raster
2815  );
2816  }
2817  else if (iwr->bandarg[i].uniontype == UT_RANGE) {
2818  noerr = rt_raster_iterator(
2819  itrset, 2,
2820  ET_UNION, NULL,
2821  pixtype,
2822  hasnodata, nodataval,
2823  0, 0,
2824  NULL,
2825  NULL,
2827  &_raster
2828  );
2829  }
2830 
2831  if (noerr != ES_NONE) {
2832  pfree(itrset);
2834  if (_rtn != NULL)
2835  rt_raster_destroy(_rtn);
2836  elog(ERROR, "RASTER_union_finalfn: Could not run raster iterator function");
2837  PG_RETURN_NULL();
2838  }
2839  }
2840  else
2841  _raster = iwr->bandarg[i].raster[0];
2842 
2843  /* first band, _rtn doesn't exist */
2844  if (i < 1) {
2845  uint32_t bandNums[1] = {0};
2846  _rtn = rt_raster_from_band(_raster, bandNums, 1);
2847  status = (_rtn == NULL) ? -1 : 0;
2848  }
2849  else
2850  status = rt_raster_copy_band(_rtn, _raster, 0, i);
2851 
2852  POSTGIS_RT_DEBUG(4, "destroying source rasters");
2853 
2854  /* destroy source rasters */
2855  if (
2856  iwr->bandarg[i].uniontype == UT_MEAN ||
2857  iwr->bandarg[i].uniontype == UT_RANGE
2858  ) {
2859  rt_raster_destroy(_raster);
2860  }
2861 
2862  for (j = 0; j < iwr->bandarg[i].numraster; j++) {
2863  if (iwr->bandarg[i].raster[j] == NULL)
2864  continue;
2865  rt_raster_destroy(iwr->bandarg[i].raster[j]);
2866  iwr->bandarg[i].raster[j] = NULL;
2867  }
2868 
2869  if (status < 0) {
2871  rt_raster_destroy(_rtn);
2872  elog(ERROR, "RASTER_union_finalfn: Could not add band to final raster");
2873  PG_RETURN_NULL();
2874  }
2875  }
2876 
2877  /* cleanup */
2878  pfree(itrset);
2880 
2881  if (!_rtn) PG_RETURN_NULL();
2882 
2883  pgraster = rt_raster_serialize(_rtn);
2884  rt_raster_destroy(_rtn);
2885 
2886  POSTGIS_RT_DEBUG(3, "Finished");
2887 
2888  if (!pgraster)
2889  PG_RETURN_NULL();
2890 
2891  SET_VARSIZE(pgraster, pgraster->size);
2892  PG_RETURN_POINTER(pgraster);
2893 }
2894 
2895 /* ---------------------------------------------------------------- */
2896 /* Clip raster with geometry */
2897 /* ---------------------------------------------------------------- */
2898 
2901  int nband; /* band index */
2902  int hasnodata; /* is there a user-specified NODATA? */
2903  double nodataval; /* user-specified NODATA */
2904 };
2905 
2911 
2912  int numbands; /* number of bandargs */
2913  rtpg_clip_band band;
2914 };
2915 
2916 static rtpg_clip_arg rtpg_clip_arg_init() {
2917  rtpg_clip_arg arg = NULL;
2918 
2919  arg = palloc(sizeof(struct rtpg_clip_arg_t));
2920  if (arg == NULL) {
2921  elog(ERROR, "rtpg_clip_arg_init: Could not allocate memory for function arguments");
2922  return NULL;
2923  }
2924 
2925  arg->extenttype = ET_INTERSECTION;
2926  arg->raster = NULL;
2927  arg->mask = NULL;
2928  arg->numbands = 0;
2929  arg->band = NULL;
2930 
2931  return arg;
2932 }
2933 
2934 static void rtpg_clip_arg_destroy(rtpg_clip_arg arg) {
2935  if (arg->band != NULL)
2936  pfree(arg->band);
2937 
2938  if (arg->raster != NULL)
2939  rt_raster_destroy(arg->raster);
2940  if (arg->mask != NULL)
2941  rt_raster_destroy(arg->mask);
2942 
2943  pfree(arg);
2944 }
2945 
2947  rt_iterator_arg arg, void *userarg,
2948  double *value, int *nodata
2949 ) {
2950  *value = 0;
2951  *nodata = 0;
2952 
2953  /* either is NODATA, output is NODATA */
2954  if (arg->nodata[0][0][0] || arg->nodata[1][0][0])
2955  *nodata = 1;
2956  /* set to value */
2957  else
2958  *value = arg->values[0][0][0];
2959 
2960  return 1;
2961 }
2962 
2964 Datum RASTER_clip(PG_FUNCTION_ARGS)
2965 {
2966  rt_pgraster *pgraster = NULL;
2967  LWGEOM *rastgeom = NULL;
2968  double gt[6] = {0};
2969  int srid = SRID_UNKNOWN;
2970 
2971  rt_pgraster *pgrtn = NULL;
2972  rt_raster rtn = NULL;
2973 
2974  GSERIALIZED *gser = NULL;
2975  LWGEOM *geom = NULL;
2976  unsigned char *wkb = NULL;
2977  size_t wkb_len;
2978 
2979  ArrayType *array;
2980  Oid etype;
2981  Datum *e;
2982  bool *nulls;
2983 
2984  int16 typlen;
2985  bool typbyval;
2986  char typalign;
2987 
2988  int i = 0;
2989  int j = 0;
2990  int k = 0;
2991  rtpg_clip_arg arg = NULL;
2992  LWGEOM *tmpgeom = NULL;
2993  rt_iterator itrset;
2994 
2995  rt_raster _raster = NULL;
2996  rt_band band = NULL;
2997  rt_pixtype pixtype;
2998  int hasnodata;
2999  double nodataval;
3000  int noerr = 0;
3001 
3002  POSTGIS_RT_DEBUG(3, "Starting...");
3003 
3004  /* raster or geometry is NULL, return NULL */
3005  if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3006  PG_RETURN_NULL();
3007 
3008  /* init arg */
3009  arg = rtpg_clip_arg_init();
3010  if (arg == NULL) {
3011  elog(ERROR, "RASTER_clip: Could not initialize argument structure");
3012  PG_RETURN_NULL();
3013  }
3014 
3015  /* raster (0) */
3016  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3017 
3018  /* Get raster object */
3019  arg->raster = rt_raster_deserialize(pgraster, FALSE);
3020  if (arg->raster == NULL) {
3021  rtpg_clip_arg_destroy(arg);
3022  PG_FREE_IF_COPY(pgraster, 0);
3023  elog(ERROR, "RASTER_clip: Could not deserialize raster");
3024  PG_RETURN_NULL();
3025  }
3026 
3027  /* raster is empty, return empty raster */
3028  if (rt_raster_is_empty(arg->raster)) {
3029  elog(NOTICE, "Input raster is empty. Returning empty raster");
3030 
3031  rtpg_clip_arg_destroy(arg);
3032  PG_FREE_IF_COPY(pgraster, 0);
3033 
3034  rtn = rt_raster_new(0, 0);
3035  if (rtn == NULL) {
3036  elog(ERROR, "RASTER_clip: Could not create empty raster");
3037  PG_RETURN_NULL();
3038  }
3039 
3040  pgrtn = rt_raster_serialize(rtn);
3041  rt_raster_destroy(rtn);
3042  if (NULL == pgrtn)
3043  PG_RETURN_NULL();
3044 
3045  SET_VARSIZE(pgrtn, pgrtn->size);
3046  PG_RETURN_POINTER(pgrtn);
3047  }
3048 
3049  /* metadata */
3051  srid = clamp_srid(rt_raster_get_srid(arg->raster));
3052 
3053  /* geometry (2) */
3054  gser = PG_GETARG_GSERIALIZED_P(2);
3055  geom = lwgeom_from_gserialized(gser);
3056 
3057  /* Get a 2D version of the geometry if necessary */
3058  if (lwgeom_ndims(geom) > 2) {
3059  LWGEOM *geom2d = lwgeom_force_2d(geom);
3060  lwgeom_free(geom);
3061  geom = geom2d;
3062  }
3063 
3064  /* check that SRIDs match */
3065  if (srid != clamp_srid(gserialized_get_srid(gser))) {
3066  elog(NOTICE, "Geometry provided does not have the same SRID as the raster. Returning NULL");
3067 
3068  rtpg_clip_arg_destroy(arg);
3069  PG_FREE_IF_COPY(pgraster, 0);
3070  lwgeom_free(geom);
3071  PG_FREE_IF_COPY(gser, 2);
3072 
3073  PG_RETURN_NULL();
3074  }
3075 
3076  /* crop (4) */
3077  if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3078  arg->extenttype = ET_FIRST;
3079 
3080  /* get intersection geometry of input raster and input geometry */
3081  if (rt_raster_get_convex_hull(arg->raster, &rastgeom) != ES_NONE) {
3082 
3083  rtpg_clip_arg_destroy(arg);
3084  PG_FREE_IF_COPY(pgraster, 0);
3085  lwgeom_free(geom);
3086  PG_FREE_IF_COPY(gser, 2);
3087 
3088  elog(ERROR, "RASTER_clip: Could not get convex hull of raster");
3089  PG_RETURN_NULL();
3090  }
3091 
3092  tmpgeom = lwgeom_intersection(rastgeom, geom);
3093  lwgeom_free(rastgeom);
3094  lwgeom_free(geom);
3095  PG_FREE_IF_COPY(gser, 2);
3096  geom = tmpgeom;
3097 
3098  /* intersection is empty AND extent type is INTERSECTION, return empty */
3099  if (lwgeom_is_empty(geom) && arg->extenttype == ET_INTERSECTION) {
3100  elog(NOTICE, "The input raster and input geometry do not intersect. Returning empty raster");
3101 
3102  rtpg_clip_arg_destroy(arg);
3103  PG_FREE_IF_COPY(pgraster, 0);
3104  lwgeom_free(geom);
3105 
3106  rtn = rt_raster_new(0, 0);
3107  if (rtn == NULL) {
3108  elog(ERROR, "RASTER_clip: Could not create empty raster");
3109  PG_RETURN_NULL();
3110  }
3111 
3112  pgrtn = rt_raster_serialize(rtn);
3113  rt_raster_destroy(rtn);
3114  if (NULL == pgrtn)
3115  PG_RETURN_NULL();
3116 
3117  SET_VARSIZE(pgrtn, pgrtn->size);
3118  PG_RETURN_POINTER(pgrtn);
3119  }
3120 
3121  /* nband (1) */
3122  if (!PG_ARGISNULL(1)) {
3123  array = PG_GETARG_ARRAYTYPE_P(1);
3124  etype = ARR_ELEMTYPE(array);
3125  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3126 
3127  switch (etype) {
3128  case INT2OID:
3129  case INT4OID:
3130  break;
3131  default:
3132  rtpg_clip_arg_destroy(arg);
3133  PG_FREE_IF_COPY(pgraster, 0);
3134  lwgeom_free(geom);
3135  elog(ERROR, "RASTER_clip: Invalid data type for band indexes");
3136  PG_RETURN_NULL();
3137  break;
3138  }
3139 
3140  deconstruct_array(
3141  array, etype,
3142  typlen, typbyval, typalign,
3143  &e, &nulls, &(arg->numbands)
3144  );
3145 
3146  arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3147  if (arg->band == NULL) {
3148  rtpg_clip_arg_destroy(arg);
3149  PG_FREE_IF_COPY(pgraster, 0);
3150  lwgeom_free(geom);
3151  elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3152  PG_RETURN_NULL();
3153  }
3154 
3155  for (i = 0, j = 0; i < arg->numbands; i++) {
3156  if (nulls[i]) continue;
3157 
3158  switch (etype) {
3159  case INT2OID:
3160  arg->band[j].nband = DatumGetInt16(e[i]) - 1;
3161  break;
3162  case INT4OID:
3163  arg->band[j].nband = DatumGetInt32(e[i]) - 1;
3164  break;
3165  }
3166 
3167  j++;
3168  }
3169 
3170  if (j < arg->numbands) {
3171  arg->band = repalloc(arg->band, sizeof(struct rtpg_clip_band_t) * j);
3172  if (arg->band == NULL) {
3173  rtpg_clip_arg_destroy(arg);
3174  PG_FREE_IF_COPY(pgraster, 0);
3175  lwgeom_free(geom);
3176  elog(ERROR, "RASTER_clip: Could not reallocate memory for band arguments");
3177  PG_RETURN_NULL();
3178  }
3179 
3180  arg->numbands = j;
3181  }
3182 
3183  /* validate band */
3184  for (i = 0; i < arg->numbands; i++) {
3185  if (!rt_raster_has_band(arg->raster, arg->band[i].nband)) {
3186  elog(NOTICE, "Band at index %d not found in raster", arg->band[i].nband + 1);
3187  rtpg_clip_arg_destroy(arg);
3188  PG_FREE_IF_COPY(pgraster, 0);
3189  lwgeom_free(geom);
3190  PG_RETURN_NULL();
3191  }
3192 
3193  arg->band[i].hasnodata = 0;
3194  arg->band[i].nodataval = 0;
3195  }
3196  }
3197  else {
3199 
3200  /* raster may have no bands */
3201  if (arg->numbands) {
3202  arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3203  if (arg->band == NULL) {
3204 
3205  rtpg_clip_arg_destroy(arg);
3206  PG_FREE_IF_COPY(pgraster, 0);
3207  lwgeom_free(geom);
3208 
3209  elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3210  PG_RETURN_NULL();
3211  }
3212 
3213  for (i = 0; i < arg->numbands; i++) {
3214  arg->band[i].nband = i;
3215  arg->band[i].hasnodata = 0;
3216  arg->band[i].nodataval = 0;
3217  }
3218  }
3219  }
3220 
3221  /* nodataval (3) */
3222  if (!PG_ARGISNULL(3)) {
3223  array = PG_GETARG_ARRAYTYPE_P(3);
3224  etype = ARR_ELEMTYPE(array);
3225  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3226 
3227  switch (etype) {
3228  case FLOAT4OID:
3229  case FLOAT8OID:
3230  break;
3231  default:
3232  rtpg_clip_arg_destroy(arg);
3233  PG_FREE_IF_COPY(pgraster, 0);
3234  lwgeom_free(geom);
3235  elog(ERROR, "RASTER_clip: Invalid data type for NODATA values");
3236  PG_RETURN_NULL();
3237  break;
3238  }
3239 
3240  deconstruct_array(
3241  array, etype,
3242  typlen, typbyval, typalign,
3243  &e, &nulls, &k
3244  );
3245 
3246  /* it doesn't matter if there are more nodataval */
3247  for (i = 0, j = 0; i < arg->numbands; i++, j++) {
3248  /* cap j to the last nodataval element */
3249  if (j >= k)
3250  j = k - 1;
3251 
3252  if (nulls[j])
3253  continue;
3254 
3255  arg->band[i].hasnodata = 1;
3256  switch (etype) {
3257  case FLOAT4OID:
3258  arg->band[i].nodataval = DatumGetFloat4(e[j]);
3259  break;
3260  case FLOAT8OID:
3261  arg->band[i].nodataval = DatumGetFloat8(e[j]);
3262  break;
3263  }
3264  }
3265  }
3266 
3267  /* get wkb of geometry */
3268  POSTGIS_RT_DEBUG(3, "getting wkb of geometry");
3269  wkb = lwgeom_to_wkb(geom, WKB_SFSQL, &wkb_len);
3270  lwgeom_free(geom);
3271 
3272  /* rasterize geometry */
3274  wkb, wkb_len,
3275  NULL,
3276  0, NULL,
3277  NULL, NULL,
3278  NULL, NULL,
3279  NULL, NULL,
3280  &(gt[1]), &(gt[5]),
3281  NULL, NULL,
3282  &(gt[0]), &(gt[3]),
3283  &(gt[2]), &(gt[4]),
3284  NULL
3285  );
3286 
3287  pfree(wkb);
3288  if (arg->mask == NULL) {
3289  rtpg_clip_arg_destroy(arg);
3290  PG_FREE_IF_COPY(pgraster, 0);
3291  elog(ERROR, "RASTER_clip: Could not rasterize intersection geometry");
3292  PG_RETURN_NULL();
3293  }
3294 
3295  /* set SRID */
3296  rt_raster_set_srid(arg->mask, srid);
3297 
3298  /* run iterator */
3299 
3300  /* init itrset */
3301  itrset = palloc(sizeof(struct rt_iterator_t) * 2);
3302  if (itrset == NULL) {
3303  rtpg_clip_arg_destroy(arg);
3304  PG_FREE_IF_COPY(pgraster, 0);
3305  elog(ERROR, "RASTER_clip: Could not allocate memory for iterator arguments");
3306  PG_RETURN_NULL();
3307  }
3308 
3309  /* one band at a time */
3310  for (i = 0; i < arg->numbands; i++) {
3311  POSTGIS_RT_DEBUGF(4, "band arg %d (nband, hasnodata, nodataval) = (%d, %d, %f)",
3312  i, arg->band[i].nband, arg->band[i].hasnodata, arg->band[i].nodataval);
3313 
3314  band = rt_raster_get_band(arg->raster, arg->band[i].nband);
3315 
3316  /* band metadata */
3317  pixtype = rt_band_get_pixtype(band);
3318 
3319  if (arg->band[i].hasnodata) {
3320  hasnodata = 1;
3321  nodataval = arg->band[i].nodataval;
3322  }
3323  else if (rt_band_get_hasnodata_flag(band)) {
3324  hasnodata = 1;
3325  rt_band_get_nodata(band, &nodataval);
3326  }
3327  else {
3328  hasnodata = 0;
3329  nodataval = rt_band_get_min_value(band);
3330  }
3331 
3332  /* band is NODATA, create NODATA band and continue */
3333  if (rt_band_get_isnodata_flag(band)) {
3334  /* create raster */
3335  if (rtn == NULL) {
3336  noerr = rt_raster_from_two_rasters(arg->raster, arg->mask, arg->extenttype, &rtn, NULL);
3337  if (noerr != ES_NONE) {
3338  rtpg_clip_arg_destroy(arg);
3339  PG_FREE_IF_COPY(pgraster, 0);
3340  elog(ERROR, "RASTER_clip: Could not create output raster");
3341  PG_RETURN_NULL();
3342  }
3343  }
3344 
3345  /* create NODATA band */
3346  if (rt_raster_generate_new_band(rtn, pixtype, nodataval, hasnodata, nodataval, i) < 0) {
3347  rt_raster_destroy(rtn);
3348  rtpg_clip_arg_destroy(arg);
3349  PG_FREE_IF_COPY(pgraster, 0);
3350  elog(ERROR, "RASTER_clip: Could not add NODATA band to output raster");
3351  PG_RETURN_NULL();
3352  }
3353 
3354  continue;
3355  }
3356 
3357  /* raster */
3358  itrset[0].raster = arg->raster;
3359  itrset[0].nband = arg->band[i].nband;
3360  itrset[0].nbnodata = 1;
3361 
3362  /* mask */
3363  itrset[1].raster = arg->mask;
3364  itrset[1].nband = 0;
3365  itrset[1].nbnodata = 1;
3366 
3367  /* pass to iterator */
3368  noerr = rt_raster_iterator(
3369  itrset, 2,
3370  arg->extenttype, NULL,
3371  pixtype,
3372  hasnodata, nodataval,
3373  0, 0,
3374  NULL,
3375  NULL,
3377  &_raster
3378  );
3379 
3380  if (noerr != ES_NONE) {
3381  pfree(itrset);
3382  rtpg_clip_arg_destroy(arg);
3383  if (rtn != NULL) rt_raster_destroy(rtn);
3384  PG_FREE_IF_COPY(pgraster, 0);
3385  elog(ERROR, "RASTER_clip: Could not run raster iterator function");
3386  PG_RETURN_NULL();
3387  }
3388 
3389  /* new raster */
3390  if (rtn == NULL)
3391  rtn = _raster;
3392  /* copy band */
3393  else {
3394  band = rt_raster_get_band(_raster, 0);
3395  if (band == NULL) {
3396  pfree(itrset);
3397  rtpg_clip_arg_destroy(arg);
3398  rt_raster_destroy(_raster);
3399  rt_raster_destroy(rtn);
3400  PG_FREE_IF_COPY(pgraster, 0);
3401  elog(ERROR, "RASTER_clip: Could not get band from working raster");
3402  PG_RETURN_NULL();
3403  }
3404 
3405  if (rt_raster_add_band(rtn, band, i) < 0) {
3406  pfree(itrset);
3407  rtpg_clip_arg_destroy(arg);
3408  rt_raster_destroy(_raster);
3409  rt_raster_destroy(rtn);
3410  PG_FREE_IF_COPY(pgraster, 0);
3411  elog(ERROR, "RASTER_clip: Could not add new band to output raster");
3412  PG_RETURN_NULL();
3413  }
3414 
3415  rt_raster_destroy(_raster);
3416  }
3417  }
3418 
3419  pfree(itrset);
3420  rtpg_clip_arg_destroy(arg);
3421  PG_FREE_IF_COPY(pgraster, 0);
3422 
3423  pgrtn = rt_raster_serialize(rtn);
3424  rt_raster_destroy(rtn);
3425 
3426  POSTGIS_RT_DEBUG(3, "Finished");
3427 
3428  if (!pgrtn)
3429  PG_RETURN_NULL();
3430 
3431  SET_VARSIZE(pgrtn, pgrtn->size);
3432  PG_RETURN_POINTER(pgrtn);
3433 }
3434 
3439 Datum RASTER_reclass(PG_FUNCTION_ARGS) {
3440  rt_pgraster *pgraster = NULL;
3441  rt_pgraster *pgrtn = NULL;
3442  rt_raster raster = NULL;
3443  rt_band band = NULL;
3444  rt_band newband = NULL;
3445  uint32_t numBands = 0;
3446 
3447  ArrayType *array;
3448  Oid etype;
3449  Datum *e;
3450  bool *nulls;
3451  int16 typlen;
3452  bool typbyval;
3453  char typalign;
3454  int n = 0;
3455 
3456  int i = 0;
3457  int j = 0;
3458  int k = 0;
3459 
3460  int a = 0;
3461  int b = 0;
3462  int c = 0;
3463 
3464  rt_reclassexpr *exprset = NULL;
3465  HeapTupleHeader tup;
3466  bool isnull;
3467  Datum tupv;
3468  uint32_t nband = 0;
3469  char *expr = NULL;
3470  text *exprtext = NULL;
3471  double val = 0;
3472  char *junk = NULL;
3473  int inc_val = 0;
3474  int exc_val = 0;
3475  char *pixeltype = NULL;
3476  text *pixeltypetext = NULL;
3477  rt_pixtype pixtype = PT_END;
3478  double nodataval = 0;
3479  bool hasnodata = FALSE;
3480 
3481  char **comma_set = NULL;
3482  int comma_n = 0;
3483  char **colon_set = NULL;
3484  int colon_n = 0;
3485  char **dash_set = NULL;
3486  int dash_n = 0;
3487 
3488  POSTGIS_RT_DEBUG(3, "RASTER_reclass: Starting");
3489 
3490  /* pgraster is null, return null */
3491  if (PG_ARGISNULL(0))
3492  PG_RETURN_NULL();
3493  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3494 
3495  /* raster */
3496  raster = rt_raster_deserialize(pgraster, FALSE);
3497  if (!raster) {
3498  PG_FREE_IF_COPY(pgraster, 0);
3499  elog(ERROR, "RASTER_reclass: Could not deserialize raster");
3500  PG_RETURN_NULL();
3501  }
3502  numBands = rt_raster_get_num_bands(raster);
3503  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: %d possible bands to be reclassified", numBands);
3504 
3505  /* process set of reclassarg */
3506  POSTGIS_RT_DEBUG(3, "RASTER_reclass: Processing Arg 1 (reclassargset)");
3507  array = PG_GETARG_ARRAYTYPE_P(1);
3508  etype = ARR_ELEMTYPE(array);
3509  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3510 
3511  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
3512  &nulls, &n);
3513 
3514  if (!n) {
3515  elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
3516 
3517  pgrtn = rt_raster_serialize(raster);
3518  rt_raster_destroy(raster);
3519  PG_FREE_IF_COPY(pgraster, 0);
3520  if (!pgrtn)
3521  PG_RETURN_NULL();
3522 
3523  SET_VARSIZE(pgrtn, pgrtn->size);
3524  PG_RETURN_POINTER(pgrtn);
3525  }
3526 
3527  /*
3528  process each element of reclassarg
3529  each element is the index of the band to process, the set
3530  of reclass ranges and the output band's pixeltype
3531  */
3532  for (i = 0; i < n; i++) {
3533  if (nulls[i]) continue;
3534 
3535  /* each element is a tuple */
3536  tup = (HeapTupleHeader) DatumGetPointer(e[i]);
3537  if (NULL == tup) {
3538  elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
3539 
3540  pgrtn = rt_raster_serialize(raster);
3541  rt_raster_destroy(raster);
3542  PG_FREE_IF_COPY(pgraster, 0);
3543  if (!pgrtn)
3544  PG_RETURN_NULL();
3545 
3546  SET_VARSIZE(pgrtn, pgrtn->size);
3547  PG_RETURN_POINTER(pgrtn);
3548  }
3549 
3550  /* band index (1-based) */
3551  tupv = GetAttributeByName(tup, "nband", &isnull);
3552  if (isnull) {
3553  elog(NOTICE, "Invalid argument for reclassargset. Missing value of nband for reclassarg of index %d . Returning original raster", i);
3554 
3555  pgrtn = rt_raster_serialize(raster);
3556  rt_raster_destroy(raster);
3557  PG_FREE_IF_COPY(pgraster, 0);
3558  if (!pgrtn)
3559  PG_RETURN_NULL();
3560 
3561  SET_VARSIZE(pgrtn, pgrtn->size);
3562  PG_RETURN_POINTER(pgrtn);
3563  }
3564  nband = DatumGetInt32(tupv);
3565  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: expression for band %d", nband);
3566 
3567  /* valid band index? */
3568  if (nband < 1 || nband > numBands) {
3569  elog(NOTICE, "Invalid argument for reclassargset. Invalid band index (must use 1-based) for reclassarg of index %d . Returning original raster", i);
3570 
3571  pgrtn = rt_raster_serialize(raster);
3572  rt_raster_destroy(raster);
3573  PG_FREE_IF_COPY(pgraster, 0);
3574  if (!pgrtn)
3575  PG_RETURN_NULL();
3576 
3577  SET_VARSIZE(pgrtn, pgrtn->size);
3578  PG_RETURN_POINTER(pgrtn);
3579  }
3580 
3581  /* reclass expr */
3582  tupv = GetAttributeByName(tup, "reclassexpr", &isnull);
3583  if (isnull) {
3584  elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3585 
3586  pgrtn = rt_raster_serialize(raster);
3587  rt_raster_destroy(raster);
3588  PG_FREE_IF_COPY(pgraster, 0);
3589  if (!pgrtn)
3590  PG_RETURN_NULL();
3591 
3592  SET_VARSIZE(pgrtn, pgrtn->size);
3593  PG_RETURN_POINTER(pgrtn);
3594  }
3595  exprtext = (text *) DatumGetPointer(tupv);
3596  if (NULL == exprtext) {
3597  elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3598 
3599  pgrtn = rt_raster_serialize(raster);
3600  rt_raster_destroy(raster);
3601  PG_FREE_IF_COPY(pgraster, 0);
3602  if (!pgrtn)
3603  PG_RETURN_NULL();
3604 
3605  SET_VARSIZE(pgrtn, pgrtn->size);
3606  PG_RETURN_POINTER(pgrtn);
3607  }
3608  expr = text_to_cstring(exprtext);
3609  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: expr (raw) %s", expr);
3610  expr = rtpg_removespaces(expr);
3611  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: expr (clean) %s", expr);
3612 
3613  /* split string to its components */
3614  /* comma (,) separating rangesets */
3615  comma_set = rtpg_strsplit(expr, ",", &comma_n);
3616  if (comma_n < 1) {
3617  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3618 
3619  pgrtn = rt_raster_serialize(raster);
3620  rt_raster_destroy(raster);
3621  PG_FREE_IF_COPY(pgraster, 0);
3622  if (!pgrtn)
3623  PG_RETURN_NULL();
3624 
3625  SET_VARSIZE(pgrtn, pgrtn->size);
3626  PG_RETURN_POINTER(pgrtn);
3627  }
3628 
3629  /* set of reclass expressions */
3630  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: %d possible expressions", comma_n);
3631  exprset = palloc(comma_n * sizeof(rt_reclassexpr));
3632 
3633  for (a = 0, j = 0; a < comma_n; a++) {
3634  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: map %s", comma_set[a]);
3635 
3636  /* colon (:) separating range "src" and "dst" */
3637  colon_set = rtpg_strsplit(comma_set[a], ":", &colon_n);
3638  if (colon_n != 2) {
3639  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3640  for (k = 0; k < j; k++) pfree(exprset[k]);
3641  pfree(exprset);
3642 
3643  pgrtn = rt_raster_serialize(raster);
3644  rt_raster_destroy(raster);
3645  PG_FREE_IF_COPY(pgraster, 0);
3646  if (!pgrtn)
3647  PG_RETURN_NULL();
3648 
3649  SET_VARSIZE(pgrtn, pgrtn->size);
3650  PG_RETURN_POINTER(pgrtn);
3651  }
3652 
3653  /* allocate mem for reclass expression */
3654  exprset[j] = palloc(sizeof(struct rt_reclassexpr_t));
3655 
3656  for (b = 0; b < colon_n; b++) {
3657  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: range %s", colon_set[b]);
3658 
3659  /* dash (-) separating "min" and "max" */
3660  dash_set = rtpg_strsplit(colon_set[b], "-", &dash_n);
3661  if (dash_n < 1 || dash_n > 3) {
3662  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3663  for (k = 0; k < j; k++) pfree(exprset[k]);
3664  pfree(exprset);
3665 
3666  pgrtn = rt_raster_serialize(raster);
3667  rt_raster_destroy(raster);
3668  PG_FREE_IF_COPY(pgraster, 0);
3669  if (!pgrtn)
3670  PG_RETURN_NULL();
3671 
3672  SET_VARSIZE(pgrtn, pgrtn->size);
3673  PG_RETURN_POINTER(pgrtn);
3674  }
3675 
3676  for (c = 0; c < dash_n; c++) {
3677  /* need to handle: (-9999-100 -> "(", "9999", "100" */
3678  if (
3679  c < 1 &&
3680  strlen(dash_set[c]) == 1 && (
3681  strchr(dash_set[c], '(') != NULL ||
3682  strchr(dash_set[c], '[') != NULL ||
3683  strchr(dash_set[c], ')') != NULL ||
3684  strchr(dash_set[c], ']') != NULL
3685  )
3686  ) {
3687  junk = palloc(sizeof(char) * (strlen(dash_set[c + 1]) + 2));
3688  if (NULL == junk) {
3689  for (k = 0; k <= j; k++) pfree(exprset[k]);
3690  pfree(exprset);
3691  rt_raster_destroy(raster);
3692  PG_FREE_IF_COPY(pgraster, 0);
3693 
3694  elog(ERROR, "RASTER_reclass: Could not allocate memory");
3695  PG_RETURN_NULL();
3696  }
3697 
3698  sprintf(junk, "%s%s", dash_set[c], dash_set[c + 1]);
3699  c++;
3700  dash_set[c] = repalloc(dash_set[c], sizeof(char) * (strlen(junk) + 1));
3701  strcpy(dash_set[c], junk);
3702  pfree(junk);
3703 
3704  /* rebuild dash_set */
3705  for (k = 1; k < dash_n; k++) {
3706  dash_set[k - 1] = repalloc(dash_set[k - 1], (strlen(dash_set[k]) + 1) * sizeof(char));
3707  strcpy(dash_set[k - 1], dash_set[k]);
3708  }
3709  dash_n--;
3710  c--;
3711  pfree(dash_set[dash_n]);
3712  dash_set = repalloc(dash_set, sizeof(char *) * dash_n);
3713  }
3714 
3715  /* there shouldn't be more than two in dash_n */
3716  if (c < 1 && dash_n > 2) {
3717  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3718  for (k = 0; k < j; k++) pfree(exprset[k]);
3719  pfree(exprset);
3720 
3721  pgrtn = rt_raster_serialize(raster);
3722  rt_raster_destroy(raster);
3723  PG_FREE_IF_COPY(pgraster, 0);
3724  if (!pgrtn)
3725  PG_RETURN_NULL();
3726 
3727  SET_VARSIZE(pgrtn, pgrtn->size);
3728  PG_RETURN_POINTER(pgrtn);
3729  }
3730 
3731  /* check interval flags */
3732  exc_val = 0;
3733  inc_val = 1;
3734  /* range */
3735  if (dash_n != 1) {
3736  /* min */
3737  if (c < 1) {
3738  if (
3739  strchr(dash_set[c], ')') != NULL ||
3740  strchr(dash_set[c], ']') != NULL
3741  ) {
3742  exc_val = 1;
3743  inc_val = 1;
3744  }
3745  else if (strchr(dash_set[c], '(') != NULL){
3746  inc_val = 0;
3747  }
3748  else {
3749  inc_val = 1;
3750  }
3751  }
3752  /* max */
3753  else {
3754  if (
3755  strrchr(dash_set[c], '(') != NULL ||
3756  strrchr(dash_set[c], '[') != NULL
3757  ) {
3758  exc_val = 1;
3759  inc_val = 0;
3760  }
3761  else if (strrchr(dash_set[c], ']') != NULL) {
3762  inc_val = 1;
3763  }
3764  else {
3765  inc_val = 0;
3766  }
3767  }
3768  }
3769  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: exc_val %d inc_val %d", exc_val, inc_val);
3770 
3771  /* remove interval flags */
3772  dash_set[c] = rtpg_chartrim(dash_set[c], "()[]");
3773  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (char) %s", dash_set[c]);
3774 
3775  /* value from string to double */
3776  errno = 0;
3777  val = strtod(dash_set[c], &junk);
3778  if (errno != 0 || dash_set[c] == junk) {
3779  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3780  for (k = 0; k < j; k++) pfree(exprset[k]);
3781  pfree(exprset);
3782 
3783  pgrtn = rt_raster_serialize(raster);
3784  rt_raster_destroy(raster);
3785  PG_FREE_IF_COPY(pgraster, 0);
3786  if (!pgrtn)
3787  PG_RETURN_NULL();
3788 
3789  SET_VARSIZE(pgrtn, pgrtn->size);
3790  PG_RETURN_POINTER(pgrtn);
3791  }
3792  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (double) %f", val);
3793 
3794  /* strsplit removes dash (a.k.a. negative sign), compare now to restore */
3795  if (c < 1)
3796  junk = strstr(colon_set[b], dash_set[c]);
3797  else
3798  junk = rtpg_strrstr(colon_set[b], dash_set[c]);
3800  4,
3801  "(colon_set[%d], dash_set[%d], junk) = (%s, %s, %s)",
3802  b, c, colon_set[b], dash_set[c], junk
3803  );
3804  /* not beginning of string */
3805  if (junk != colon_set[b]) {
3806  /* prior is a dash */
3807  if (*(junk - 1) == '-') {
3808  /* prior is beginning of string or prior - 1 char is dash, negative number */
3809  if (
3810  ((junk - 1) == colon_set[b]) ||
3811  (*(junk - 2) == '-') ||
3812  (*(junk - 2) == '[') ||
3813  (*(junk - 2) == '(')
3814  ) {
3815  val *= -1.;
3816  }
3817  }
3818  }
3819  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (double) %f", val);
3820 
3821  /* src */
3822  if (b < 1) {
3823  /* singular value */
3824  if (dash_n == 1) {
3825  exprset[j]->src.exc_min = exprset[j]->src.exc_max = exc_val;
3826  exprset[j]->src.inc_min = exprset[j]->src.inc_max = inc_val;
3827  exprset[j]->src.min = exprset[j]->src.max = val;
3828  }
3829  /* min */
3830  else if (c < 1) {
3831  exprset[j]->src.exc_min = exc_val;
3832  exprset[j]->src.inc_min = inc_val;
3833  exprset[j]->src.min = val;
3834  }
3835  /* max */
3836  else {
3837  exprset[j]->src.exc_max = exc_val;
3838  exprset[j]->src.inc_max = inc_val;
3839  exprset[j]->src.max = val;
3840  }
3841  }
3842  /* dst */
3843  else {
3844  /* singular value */
3845  if (dash_n == 1)
3846  exprset[j]->dst.min = exprset[j]->dst.max = val;
3847  /* min */
3848  else if (c < 1)
3849  exprset[j]->dst.min = val;
3850  /* max */
3851  else
3852  exprset[j]->dst.max = val;
3853  }
3854  }
3855  pfree(dash_set);
3856  }
3857  pfree(colon_set);
3858 
3859  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: or: %f - %f nr: %f - %f"
3860  , exprset[j]->src.min
3861  , exprset[j]->src.max
3862  , exprset[j]->dst.min
3863  , exprset[j]->dst.max
3864  );
3865  j++;
3866  }
3867  pfree(comma_set);
3868 
3869  /* pixel type */
3870  tupv = GetAttributeByName(tup, "pixeltype", &isnull);
3871  if (isnull) {
3872  elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3873 
3874  pgrtn = rt_raster_serialize(raster);
3875  rt_raster_destroy(raster);
3876  PG_FREE_IF_COPY(pgraster, 0);
3877  if (!pgrtn)
3878  PG_RETURN_NULL();
3879 
3880  SET_VARSIZE(pgrtn, pgrtn->size);
3881  PG_RETURN_POINTER(pgrtn);
3882  }
3883  pixeltypetext = (text *) DatumGetPointer(tupv);
3884  if (NULL == pixeltypetext) {
3885  elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3886 
3887  pgrtn = rt_raster_serialize(raster);
3888  rt_raster_destroy(raster);
3889  PG_FREE_IF_COPY(pgraster, 0);
3890  if (!pgrtn)
3891  PG_RETURN_NULL();
3892 
3893  SET_VARSIZE(pgrtn, pgrtn->size);
3894  PG_RETURN_POINTER(pgrtn);
3895  }
3896  pixeltype = text_to_cstring(pixeltypetext);
3897  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: pixeltype %s", pixeltype);
3898  pixtype = rt_pixtype_index_from_name(pixeltype);
3899 
3900  /* nodata */
3901  tupv = GetAttributeByName(tup, "nodataval", &isnull);
3902  if (isnull) {
3903  nodataval = 0;
3904  hasnodata = FALSE;
3905  }
3906  else {
3907  nodataval = DatumGetFloat8(tupv);
3908  hasnodata = TRUE;
3909  }
3910  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: nodataval %f", nodataval);
3911  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: hasnodata %d", hasnodata);
3912 
3913  /* do reclass */
3914  band = rt_raster_get_band(raster, nband - 1);
3915  if (!band) {
3916  elog(NOTICE, "Could not find raster band of index %d. Returning original raster", nband);
3917  for (k = 0; k < j; k++) pfree(exprset[k]);
3918  pfree(exprset);
3919 
3920  pgrtn = rt_raster_serialize(raster);
3921  rt_raster_destroy(raster);
3922  PG_FREE_IF_COPY(pgraster, 0);
3923  if (!pgrtn)
3924  PG_RETURN_NULL();
3925 
3926  SET_VARSIZE(pgrtn, pgrtn->size);
3927  PG_RETURN_POINTER(pgrtn);
3928  }
3929  newband = rt_band_reclass(band, pixtype, hasnodata, nodataval, exprset, j);
3930  if (!newband) {
3931  for (k = 0; k < j; k++) pfree(exprset[k]);
3932  pfree(exprset);
3933 
3934  rt_raster_destroy(raster);
3935  PG_FREE_IF_COPY(pgraster, 0);
3936 
3937  elog(ERROR, "RASTER_reclass: Could not reclassify raster band of index %d", nband);
3938  PG_RETURN_NULL();
3939  }
3940 
3941  /* replace old band with new band */
3942  if (rt_raster_replace_band(raster, newband, nband - 1) == NULL) {
3943  for (k = 0; k < j; k++) pfree(exprset[k]);
3944  pfree(exprset);
3945 
3946  rt_band_destroy(newband);
3947  rt_raster_destroy(raster);
3948  PG_FREE_IF_COPY(pgraster, 0);
3949 
3950  elog(ERROR, "RASTER_reclass: Could not replace raster band of index %d with reclassified band", nband);
3951  PG_RETURN_NULL();
3952  }
3953 
3954  /* old band is in the variable band */
3955  rt_band_destroy(band);
3956 
3957  /* free exprset */
3958  for (k = 0; k < j; k++) pfree(exprset[k]);
3959  pfree(exprset);
3960  }
3961 
3962  pgrtn = rt_raster_serialize(raster);
3963  rt_raster_destroy(raster);
3964  PG_FREE_IF_COPY(pgraster, 0);
3965  if (!pgrtn)
3966  PG_RETURN_NULL();
3967 
3968  POSTGIS_RT_DEBUG(3, "RASTER_reclass: Finished");
3969 
3970  SET_VARSIZE(pgrtn, pgrtn->size);
3971  PG_RETURN_POINTER(pgrtn);
3972 }
3973 
3974 /* ---------------------------------------------------------------- */
3975 /* apply colormap to specified band of a raster */
3976 /* ---------------------------------------------------------------- */
3977 
3981  int nband; /* 1-based */
3984 
3987 
3988  char **entry;
3989  int nentry;
3990  char **element;
3992 };
3993 
3994 static rtpg_colormap_arg
3996  rtpg_colormap_arg arg = NULL;
3997 
3998  arg = palloc(sizeof(struct rtpg_colormap_arg_t));
3999  if (arg == NULL) {
4000  elog(ERROR, "rtpg_colormap_arg: Could not allocate memory for function arguments");
4001  return NULL;
4002  }
4003 
4004  arg->raster = NULL;
4005  arg->nband = 1;
4006  arg->band = NULL;
4007  arg->bandstats = NULL;
4008 
4009  arg->colormap = palloc(sizeof(struct rt_colormap_t));
4010  if (arg->colormap == NULL) {
4011  elog(ERROR, "rtpg_colormap_arg: Could not allocate memory for function arguments");
4012  return NULL;
4013  }
4014  arg->colormap->nentry = 0;
4015  arg->colormap->entry = NULL;
4016  arg->colormap->ncolor = 4; /* assume RGBA */
4017  arg->colormap->method = CM_INTERPOLATE;
4018  arg->nodataentry = -1;
4019 
4020  arg->entry = NULL;
4021  arg->nentry = 0;
4022  arg->element = NULL;
4023  arg->nelement = 0;
4024 
4025  return arg;
4026 }
4027 
4028 static void
4029 rtpg_colormap_arg_destroy(rtpg_colormap_arg arg) {
4030  int i = 0;
4031  if (arg->raster != NULL)
4032  rt_raster_destroy(arg->raster);
4033 
4034  if (arg->bandstats != NULL)
4035  pfree(arg->bandstats);
4036 
4037  if (arg->colormap != NULL) {
4038  if (arg->colormap->entry != NULL)
4039  pfree(arg->colormap->entry);
4040  pfree(arg->colormap);
4041  }
4042 
4043  if (arg->nentry) {
4044  for (i = 0; i < arg->nentry; i++) {
4045  if (arg->entry[i] != NULL)
4046  pfree(arg->entry[i]);
4047  }
4048  pfree(arg->entry);
4049  }
4050 
4051  if (arg->nelement) {
4052  for (i = 0; i < arg->nelement; i++)
4053  pfree(arg->element[i]);
4054  pfree(arg->element);
4055  }
4056 
4057  pfree(arg);
4058  arg = NULL;
4059 }
4060 
4062 Datum RASTER_colorMap(PG_FUNCTION_ARGS)
4063 {
4064  rt_pgraster *pgraster = NULL;
4065  rtpg_colormap_arg arg = NULL;
4066  char *junk = NULL;
4067  rt_raster raster = NULL;
4068 
4069  POSTGIS_RT_DEBUG(3, "RASTER_colorMap: Starting");
4070 
4071  /* pgraster is NULL, return NULL */
4072  if (PG_ARGISNULL(0))
4073  PG_RETURN_NULL();
4074 
4075  /* init arg */
4076  arg = rtpg_colormap_arg_init();
4077  if (arg == NULL) {
4078  elog(ERROR, "RASTER_colorMap: Could not initialize argument structure");
4079  PG_RETURN_NULL();
4080  }
4081 
4082  /* raster (0) */
4083  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4084 
4085  /* raster object */
4086  arg->raster = rt_raster_deserialize(pgraster, FALSE);
4087  if (!arg->raster) {
4089  PG_FREE_IF_COPY(pgraster, 0);
4090  elog(ERROR, "RASTER_colorMap: Could not deserialize raster");
4091  PG_RETURN_NULL();
4092  }
4093 
4094  /* nband (1) */
4095  if (!PG_ARGISNULL(1))
4096  arg->nband = PG_GETARG_INT32(1);
4097  POSTGIS_RT_DEBUGF(4, "nband = %d", arg->nband);
4098 
4099  /* check that band works */
4100  if (!rt_raster_has_band(arg->raster, arg->nband - 1)) {
4101  elog(NOTICE, "Raster does not have band at index %d. Returning empty raster", arg->nband);
4102 
4103  raster = rt_raster_clone(arg->raster, 0);
4104  if (raster == NULL) {
4106  PG_FREE_IF_COPY(pgraster, 0);
4107  elog(ERROR, "RASTER_colorMap: Could not create empty raster");
4108  PG_RETURN_NULL();
4109  }
4110 
4112  PG_FREE_IF_COPY(pgraster, 0);
4113 
4114  pgraster = rt_raster_serialize(raster);
4115  rt_raster_destroy(raster);
4116  if (pgraster == NULL)
4117  PG_RETURN_NULL();
4118 
4119  SET_VARSIZE(pgraster, ((rt_pgraster*) pgraster)->size);
4120  PG_RETURN_POINTER(pgraster);
4121  }
4122 
4123  /* get band */
4124  arg->band = rt_raster_get_band(arg->raster, arg->nband - 1);
4125  if (arg->band == NULL) {
4126  int nband = arg->nband;
4128  PG_FREE_IF_COPY(pgraster, 0);
4129  elog(ERROR, "RASTER_colorMap: Could not get band at index %d", nband);
4130  PG_RETURN_NULL();
4131  }
4132 
4133  /* method (3) */
4134  if (!PG_ARGISNULL(3)) {
4135  char *method = NULL;
4136  char *tmp = text_to_cstring(PG_GETARG_TEXT_P(3));
4137  POSTGIS_RT_DEBUGF(4, "raw method = %s", tmp);
4138 
4139  method = rtpg_trim(tmp);
4140  pfree(tmp);
4141  method = rtpg_strtoupper(method);
4142 
4143  if (strcmp(method, "INTERPOLATE") == 0)
4144  arg->colormap->method = CM_INTERPOLATE;
4145  else if (strcmp(method, "EXACT") == 0)
4146  arg->colormap->method = CM_EXACT;
4147  else if (strcmp(method, "NEAREST") == 0)
4148  arg->colormap->method = CM_NEAREST;
4149  else {
4150  elog(NOTICE, "Unknown value provided for method. Defaulting to INTERPOLATE");
4151  arg->colormap->method = CM_INTERPOLATE;
4152  }
4153  }
4154  /* default to INTERPOLATE */
4155  else
4156  arg->colormap->method = CM_INTERPOLATE;
4157  POSTGIS_RT_DEBUGF(4, "method = %d", arg->colormap->method);
4158 
4159  /* colormap (2) */
4160  if (PG_ARGISNULL(2)) {
4162  PG_FREE_IF_COPY(pgraster, 0);
4163  elog(ERROR, "RASTER_colorMap: Value must be provided for colormap");
4164  PG_RETURN_NULL();
4165  }
4166  else {
4167  char *tmp = NULL;
4168  char *colormap = text_to_cstring(PG_GETARG_TEXT_P(2));
4169  char *_entry;
4170  char *_element;
4171  int i = 0;
4172  int j = 0;
4173 
4174  POSTGIS_RT_DEBUGF(4, "colormap = %s", colormap);
4175 
4176  /* empty string */
4177  if (!strlen(colormap)) {
4179  PG_FREE_IF_COPY(pgraster, 0);
4180  elog(ERROR, "RASTER_colorMap: Value must be provided for colormap");
4181  PG_RETURN_NULL();
4182  }
4183 
4184  arg->entry = rtpg_strsplit(colormap, "\n", &(arg->nentry));
4185  pfree(colormap);
4186  if (arg->nentry < 1) {
4188  PG_FREE_IF_COPY(pgraster, 0);
4189  elog(ERROR, "RASTER_colorMap: Could not process the value provided for colormap");
4190  PG_RETURN_NULL();
4191  }
4192 
4193  /* allocate the max # of colormap entries */
4194  arg->colormap->entry = palloc(sizeof(struct rt_colormap_entry_t) * arg->nentry);
4195  if (arg->colormap->entry == NULL) {
4197  PG_FREE_IF_COPY(pgraster, 0);
4198  elog(ERROR, "RASTER_colorMap: Could not allocate memory for colormap entries");
4199  PG_RETURN_NULL();
4200  }
4201  memset(arg->colormap->entry, 0, sizeof(struct rt_colormap_entry_t) * arg->nentry);
4202 
4203  /* each entry */
4204  for (i = 0; i < arg->nentry; i++) {
4205  /* substitute space for other delimiters */
4206  tmp = rtpg_strreplace(arg->entry[i], ":", " ", NULL);
4207  _entry = rtpg_strreplace(tmp, ",", " ", NULL);
4208  pfree(tmp);
4209  tmp = rtpg_strreplace(_entry, "\t", " ", NULL);
4210  pfree(_entry);
4211  _entry = rtpg_trim(tmp);
4212  pfree(tmp);
4213 
4214  POSTGIS_RT_DEBUGF(4, "Processing entry[%d] = %s", i, arg->entry[i]);
4215  POSTGIS_RT_DEBUGF(4, "Cleaned entry[%d] = %s", i, _entry);
4216 
4217  /* empty entry, continue */
4218  if (!strlen(_entry)) {
4219  POSTGIS_RT_DEBUGF(3, "Skipping empty entry[%d]", i);
4220  pfree(_entry);
4221  continue;
4222  }
4223 
4224  arg->element = rtpg_strsplit(_entry, " ", &(arg->nelement));
4225  pfree(_entry);
4226  if (arg->nelement < 2) {
4228  PG_FREE_IF_COPY(pgraster, 0);
4229  elog(ERROR, "RASTER_colorMap: Could not process colormap entry %d", i + 1);
4230  PG_RETURN_NULL();
4231  }
4232  else if (arg->nelement > 5) {
4233  elog(NOTICE, "More than five elements in colormap entry %d. Using at most five elements", i + 1);
4234  arg->nelement = 5;
4235  }
4236 
4237  /* smallest # of colors */
4238  if ((arg->nelement - 1) < arg->colormap->ncolor)
4239  arg->colormap->ncolor = arg->nelement - 1;
4240 
4241  /* each element of entry */
4242  for (j = 0; j < arg->nelement; j++) {
4243 
4244  _element = rtpg_trim(arg->element[j]);
4245  _element = rtpg_strtoupper(_element);
4246  POSTGIS_RT_DEBUGF(4, "Processing entry[%d][%d] = %s", i, j, arg->element[j]);
4247  POSTGIS_RT_DEBUGF(4, "Cleaned entry[%d][%d] = %s", i, j, _element);
4248 
4249  /* first element is ALWAYS a band value, percentage OR "nv" string */
4250  if (j == 0) {
4251  char *percent = NULL;
4252 
4253  /* NODATA */
4254  if (
4255  strcmp(_element, "NV") == 0 ||
4256  strcmp(_element, "NULL") == 0 ||
4257  strcmp(_element, "NODATA") == 0
4258  ) {
4259  POSTGIS_RT_DEBUG(4, "Processing NODATA string");
4260 
4261  if (arg->nodataentry > -1) {
4262  elog(NOTICE, "More than one NODATA entry found. Using only the first one");
4263  }
4264  else {
4265  arg->colormap->entry[arg->colormap->nentry].isnodata = 1;
4266  /* no need to set value as value comes from band's NODATA */
4267  arg->colormap->entry[arg->colormap->nentry].value = 0;
4268  }
4269  }
4270  /* percent value */
4271  else if ((percent = strchr(_element, '%')) != NULL) {
4272  double value;
4273  POSTGIS_RT_DEBUG(4, "Processing percent string");
4274 
4275  /* get the band stats */
4276  if (arg->bandstats == NULL) {
4277  POSTGIS_RT_DEBUG(4, "Getting band stats");
4278 
4279  arg->bandstats = rt_band_get_summary_stats(arg->band, 1, 1, 0, NULL, NULL, NULL);
4280  if (arg->bandstats == NULL) {
4281  pfree(_element);
4283  PG_FREE_IF_COPY(pgraster, 0);
4284  elog(ERROR, "RASTER_colorMap: Could not get band's summary stats to process percentages");
4285  PG_RETURN_NULL();
4286  }
4287  }
4288 
4289  /* get the string before the percent char */
4290  tmp = palloc(sizeof(char) * (percent - _element + 1));
4291  if (tmp == NULL) {
4292  pfree(_element);
4294  PG_FREE_IF_COPY(pgraster, 0);
4295  elog(ERROR, "RASTER_colorMap: Could not allocate memory for value of percentage");
4296  PG_RETURN_NULL();
4297  }
4298 
4299  memcpy(tmp, _element, percent - _element);
4300  tmp[percent - _element] = '\0';
4301  POSTGIS_RT_DEBUGF(4, "Percent value = %s", tmp);
4302 
4303  /* get percentage value */
4304  errno = 0;
4305  value = strtod(tmp, NULL);
4306  pfree(tmp);
4307  if (errno != 0 || _element == junk) {
4308  pfree(_element);
4310  PG_FREE_IF_COPY(pgraster, 0);
4311  elog(ERROR, "RASTER_colorMap: Could not process percent string to value");
4312  PG_RETURN_NULL();
4313  }
4314 
4315  /* check percentage */
4316  if (value < 0.) {
4317  elog(NOTICE, "Percentage values cannot be less than zero. Defaulting to zero");
4318  value = 0.;
4319  }
4320  else if (value > 100.) {
4321  elog(NOTICE, "Percentage values cannot be greater than 100. Defaulting to 100");
4322  value = 100.;
4323  }
4324 
4325  /* get the true pixel value */
4326  /* TODO: should the percentage be quantile based? */
4327  arg->colormap->entry[arg->colormap->nentry].value = ((value / 100.) * (arg->bandstats->max - arg->bandstats->min)) + arg->bandstats->min;
4328  }
4329  /* straight value */
4330  else {
4331  errno = 0;
4332  arg->colormap->entry[arg->colormap->nentry].value = strtod(_element, &junk);
4333  if (errno != 0 || _element == junk) {
4334  pfree(_element);
4336  PG_FREE_IF_COPY(pgraster, 0);
4337  elog(ERROR, "RASTER_colorMap: Could not process string to value");
4338  PG_RETURN_NULL();
4339  }
4340  }
4341 
4342  }
4343  /* RGB values (0 - 255) */
4344  else {
4345  int value = 0;
4346 
4347  errno = 0;
4348  value = (int) strtod(_element, &junk);
4349  if (errno != 0 || _element == junk) {
4350  pfree(_element);
4352  PG_FREE_IF_COPY(pgraster, 0);
4353  elog(ERROR, "RASTER_colorMap: Could not process string to value");
4354  PG_RETURN_NULL();
4355  }
4356 
4357  if (value > 255) {
4358  elog(NOTICE, "RGBA value cannot be greater than 255. Defaulting to 255");
4359  value = 255;
4360  }
4361  else if (value < 0) {
4362  elog(NOTICE, "RGBA value cannot be less than zero. Defaulting to zero");
4363  value = 0;
4364  }
4365  arg->colormap->entry[arg->colormap->nentry].color[j - 1] = value;
4366  }
4367 
4368  pfree(_element);
4369  }
4370 
4371  POSTGIS_RT_DEBUGF(4, "colormap->entry[%d] (isnodata, value, R, G, B, A) = (%d, %f, %d, %d, %d, %d)",
4372  arg->colormap->nentry,
4373  arg->colormap->entry[arg->colormap->nentry].isnodata,
4374  arg->colormap->entry[arg->colormap->nentry].value,
4375  arg->colormap->entry[arg->colormap->nentry].color[0],
4376  arg->colormap->entry[arg->colormap->nentry].color[1],
4377  arg->colormap->entry[arg->colormap->nentry].color[2],
4378  arg->colormap->entry[arg->colormap->nentry].color[3]
4379  );
4380 
4381  arg->colormap->nentry++;
4382  }
4383 
4384  POSTGIS_RT_DEBUGF(4, "colormap->nentry = %d", arg->colormap->nentry);
4385  POSTGIS_RT_DEBUGF(4, "colormap->ncolor = %d", arg->colormap->ncolor);
4386  }
4387 
4388  /* call colormap */
4389  raster = rt_raster_colormap(arg->raster, arg->nband - 1, arg->colormap);
4390  if (raster == NULL) {
4392  PG_FREE_IF_COPY(pgraster, 0);
4393  elog(ERROR, "RASTER_colorMap: Could not create new raster with applied colormap");
4394  PG_RETURN_NULL();
4395  }
4396 
4398  PG_FREE_IF_COPY(pgraster, 0);
4399  pgraster = rt_raster_serialize(raster);
4400  rt_raster_destroy(raster);
4401 
4402  POSTGIS_RT_DEBUG(3, "RASTER_colorMap: Done");
4403 
4404  if (pgraster == NULL)
4405  PG_RETURN_NULL();
4406 
4407  SET_VARSIZE(pgraster, ((rt_pgraster*) pgraster)->size);
4408  PG_RETURN_POINTER(pgraster);
4409 }
4410 
4412 Datum RASTER_mapAlgebraExpr(PG_FUNCTION_ARGS)
4413 {
4414  rt_pgraster *pgraster = NULL;
4415  rt_pgraster *pgrtn = NULL;
4416  rt_raster raster = NULL;
4417  rt_raster newrast = NULL;
4418  rt_band band = NULL;
4419  rt_band newband = NULL;
4420  int x, y, nband, width, height;
4421  double r;
4422  double newnodatavalue = 0.0;
4423  double newinitialvalue = 0.0;
4424  double newval = 0.0;
4425  char *newexpr = NULL;
4426  char *initexpr = NULL;
4427  char *expression = NULL;
4428  int hasnodataval = 0;
4429  double nodataval = 0.;
4430  rt_pixtype newpixeltype;
4431  int skipcomputation = 0;
4432  int len = 0;
4433  const int argkwcount = 3;
4434  enum KEYWORDS { kVAL=0, kX=1, kY=2 };
4435  char *argkw[] = {"[rast]", "[rast.x]", "[rast.y]"};
4436  Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
4437  int argcount = 0;
4438  Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
4439  uint8_t argpos[3] = {0};
4440  char place[5];
4441  int idx = 0;
4442  int ret = -1;
4443  TupleDesc tupdesc;
4444  SPIPlanPtr spi_plan = NULL;
4445  SPITupleTable * tuptable = NULL;
4446  HeapTuple tuple;
4447  char * strFromText = NULL;
4448  Datum *values = NULL;
4449  Datum datum = (Datum)NULL;
4450  char *nulls = NULL;
4451  bool isnull = FALSE;
4452  int i = 0;
4453  int j = 0;
4454 
4455  POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraExpr: Starting...");
4456 
4457  /* Check raster */
4458  if (PG_ARGISNULL(0)) {
4459  elog(NOTICE, "Raster is NULL. Returning NULL");
4460  PG_RETURN_NULL();
4461  }
4462 
4463 
4464  /* Deserialize raster */
4465  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4466  raster = rt_raster_deserialize(pgraster, FALSE);
4467  if (NULL == raster) {
4468  PG_FREE_IF_COPY(pgraster, 0);
4469  elog(ERROR, "RASTER_mapAlgebraExpr: Could not deserialize raster");
4470  PG_RETURN_NULL();
4471  }
4472 
4473  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting arguments...");
4474 
4475  if (PG_ARGISNULL(1))
4476  nband = 1;
4477  else
4478  nband = PG_GETARG_INT32(1);
4479 
4480  if (nband < 1)
4481  nband = 1;
4482 
4483 
4484  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new empty raster...");
4485 
4490  width = rt_raster_get_width(raster);
4491  height = rt_raster_get_height(raster);
4492 
4493  newrast = rt_raster_new(width, height);
4494 
4495  if ( NULL == newrast ) {
4496  PG_FREE_IF_COPY(pgraster, 0);
4497  elog(ERROR, "RASTER_mapAlgebraExpr: Could not create a new raster");
4498  PG_RETURN_NULL();
4499  }
4500 
4501  rt_raster_set_scale(newrast,
4502  rt_raster_get_x_scale(raster),
4503  rt_raster_get_y_scale(raster));
4504 
4505  rt_raster_set_offsets(newrast,
4506  rt_raster_get_x_offset(raster),
4507  rt_raster_get_y_offset(raster));
4508 
4509  rt_raster_set_skews(newrast,
4510  rt_raster_get_x_skew(raster),
4511  rt_raster_get_y_skew(raster));
4512 
4513  rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
4514 
4515 
4520  if (rt_raster_is_empty(newrast))
4521  {
4522  elog(NOTICE, "Raster is empty. Returning an empty raster");
4523  rt_raster_destroy(raster);
4524  PG_FREE_IF_COPY(pgraster, 0);
4525 
4526  pgrtn = rt_raster_serialize(newrast);
4527  rt_raster_destroy(newrast);
4528  if (NULL == pgrtn) {
4529 
4530  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4531  PG_RETURN_NULL();
4532  }
4533 
4534  SET_VARSIZE(pgrtn, pgrtn->size);
4535  PG_RETURN_POINTER(pgrtn);
4536  }
4537 
4538 
4539  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Getting raster band %d...", nband);
4540 
4545  if (!rt_raster_has_band(raster, nband - 1)) {
4546  elog(NOTICE, "Raster does not have the required band. Returning a raster "
4547  "without a band");
4548  rt_raster_destroy(raster);
4549  PG_FREE_IF_COPY(pgraster, 0);
4550 
4551  pgrtn = rt_raster_serialize(newrast);
4552  rt_raster_destroy(newrast);
4553  if (NULL == pgrtn) {
4554  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4555  PG_RETURN_NULL();
4556  }
4557 
4558  SET_VARSIZE(pgrtn, pgrtn->size);
4559  PG_RETURN_POINTER(pgrtn);
4560  }
4561 
4562  /* Get the raster band */
4563  band = rt_raster_get_band(raster, nband - 1);
4564  if ( NULL == band ) {
4565  elog(NOTICE, "Could not get the required band. Returning a raster "
4566  "without a band");
4567  rt_raster_destroy(raster);
4568  PG_FREE_IF_COPY(pgraster, 0);
4569 
4570  pgrtn = rt_raster_serialize(newrast);
4571  rt_raster_destroy(newrast);
4572  if (NULL == pgrtn) {
4573  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4574  PG_RETURN_NULL();
4575  }
4576 
4577  SET_VARSIZE(pgrtn, pgrtn->size);
4578  PG_RETURN_POINTER(pgrtn);
4579  }
4580 
4581  /*
4582  * Get NODATA value
4583  */
4584  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting NODATA value for band...");
4585 
4586  if (rt_band_get_hasnodata_flag(band)) {
4587  rt_band_get_nodata(band, &newnodatavalue);
4588  }
4589 
4590  else {
4591  newnodatavalue = rt_band_get_min_value(band);
4592  }
4593 
4594  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: NODATA value for band: = %f",
4595  newnodatavalue);
4596 
4602  newinitialvalue = newnodatavalue;
4603 
4607  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Setting pixeltype...");
4608 
4609  if (PG_ARGISNULL(2)) {
4610  newpixeltype = rt_band_get_pixtype(band);
4611  }
4612 
4613  else {
4614  strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
4615  newpixeltype = rt_pixtype_index_from_name(strFromText);
4616  pfree(strFromText);
4617  if (newpixeltype == PT_END)
4618  newpixeltype = rt_band_get_pixtype(band);
4619  }
4620 
4621  if (newpixeltype == PT_END) {
4622  PG_FREE_IF_COPY(pgraster, 0);
4623  elog(ERROR, "RASTER_mapAlgebraExpr: Invalid pixeltype");
4624  PG_RETURN_NULL();
4625  }
4626 
4627  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Pixeltype set to %s",
4628  rt_pixtype_name(newpixeltype));
4629 
4630 
4631  /* Construct expression for raster values */
4632  if (!PG_ARGISNULL(3)) {
4633  expression = text_to_cstring(PG_GETARG_TEXT_P(3));
4634  len = strlen("SELECT (") + strlen(expression) + strlen(")::double precision");
4635  initexpr = (char *)palloc(len + 1);
4636 
4637  strncpy(initexpr, "SELECT (", strlen("SELECT ("));
4638  strncpy(initexpr + strlen("SELECT ("), expression, strlen(expression));
4639  strncpy(initexpr + strlen("SELECT (") + strlen(expression), ")::double precision", strlen(")::double precision"));
4640  initexpr[len] = '\0';
4641 
4642  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression is %s", initexpr);
4643 
4644  /* We don't need this memory */
4645  /*
4646  pfree(expression);
4647  expression = NULL;
4648  */
4649  }
4650 
4651 
4652 
4658  if (!PG_ARGISNULL(4)) {
4659  hasnodataval = 1;
4660  nodataval = PG_GETARG_FLOAT8(4);
4661  newinitialvalue = nodataval;
4662 
4663  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new initial value = %f",
4664  newinitialvalue);
4665  }
4666  else
4667  hasnodataval = 0;
4668 
4669 
4670 
4676  if (rt_band_get_isnodata_flag(band)) {
4677 
4678  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Band is a nodata band, returning "
4679  "a raster filled with nodata");
4680 
4681  ret = rt_raster_generate_new_band(newrast, newpixeltype,
4682  newinitialvalue, TRUE, newnodatavalue, 0);
4683 
4684  /* Free memory */
4685  if (initexpr)
4686  pfree(initexpr);
4687  rt_raster_destroy(raster);
4688  PG_FREE_IF_COPY(pgraster, 0);
4689 
4690  /* Serialize created raster */
4691  pgrtn = rt_raster_serialize(newrast);
4692  rt_raster_destroy(newrast);
4693  if (NULL == pgrtn) {
4694  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4695  PG_RETURN_NULL();
4696  }
4697 
4698  SET_VARSIZE(pgrtn, pgrtn->size);
4699  PG_RETURN_POINTER(pgrtn);
4700  }
4701 
4702 
4707  if (initexpr != NULL && ( !strcmp(initexpr, "SELECT [rast]") || !strcmp(initexpr, "SELECT [rast.val]") ) && !hasnodataval) {
4708 
4709  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression resumes to RAST. "
4710  "Returning raster with band %d from original raster", nband);
4711 
4712  POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster has %d bands",
4713  rt_raster_get_num_bands(newrast));
4714 
4715  rt_raster_copy_band(newrast, raster, nband - 1, 0);
4716 
4717  POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster now has %d bands",
4718  rt_raster_get_num_bands(newrast));
4719 
4720  if (initexpr)
4721  pfree(initexpr);
4722  rt_raster_destroy(raster);
4723  PG_FREE_IF_COPY(pgraster, 0);
4724 
4725  /* Serialize created raster */
4726  pgrtn = rt_raster_serialize(newrast);
4727  rt_raster_destroy(newrast);
4728  if (NULL == pgrtn) {
4729  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4730  PG_RETURN_NULL();
4731  }
4732 
4733  SET_VARSIZE(pgrtn, pgrtn->size);
4734  PG_RETURN_POINTER(pgrtn);
4735  }
4736 
4741  if (initexpr != NULL && strstr(initexpr, "[rast") == NULL) {
4742  ret = SPI_connect();
4743  if (ret != SPI_OK_CONNECT) {
4744  PG_FREE_IF_COPY(pgraster, 0);
4745  elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4746  PG_RETURN_NULL();
4747  };
4748 
4749  /* Execute the expresion into newval */
4750  ret = SPI_execute(initexpr, FALSE, 0);
4751 
4752  if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
4753 
4754  /* Free memory allocated out of the current context */
4755  if (SPI_tuptable)
4756  SPI_freetuptable(tuptable);
4757  PG_FREE_IF_COPY(pgraster, 0);
4758 
4759  SPI_finish();
4760  elog(ERROR, "RASTER_mapAlgebraExpr: Invalid construction for expression");
4761  PG_RETURN_NULL();
4762  }
4763 
4764  tupdesc = SPI_tuptable->tupdesc;
4765  tuptable = SPI_tuptable;
4766 
4767  tuple = tuptable->vals[0];
4768  newexpr = SPI_getvalue(tuple, tupdesc, 1);
4769  if ( ! newexpr ) {
4770  POSTGIS_RT_DEBUG(3, "Constant expression evaluated to NULL, keeping initvalue");
4771  newval = newinitialvalue;
4772  } else {
4773  newval = atof(newexpr);
4774  }
4775 
4776  SPI_freetuptable(tuptable);
4777 
4778  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: New raster value = %f",
4779  newval);
4780 
4781  SPI_finish();
4782 
4783  skipcomputation = 1;
4784 
4789  if (!hasnodataval) {
4790  newinitialvalue = newval;
4791  skipcomputation = 2;
4792  }
4793 
4794  /* Return the new raster as it will be before computing pixel by pixel */
4795  else if (FLT_NEQ(newval, newinitialvalue)) {
4796  skipcomputation = 2;
4797  }
4798  }
4799 
4804  ret = rt_raster_generate_new_band(newrast, newpixeltype,
4805  newinitialvalue, TRUE, newnodatavalue, 0);
4806 
4811  /*if (initexpr == NULL || skipcomputation == 2) {*/
4812  if (expression == NULL || skipcomputation == 2) {
4813 
4814  /* Free memory */
4815  if (initexpr)
4816  pfree(initexpr);
4817  rt_raster_destroy(raster);
4818  PG_FREE_IF_COPY(pgraster, 0);
4819 
4820  /* Serialize created raster */
4821  pgrtn = rt_raster_serialize(newrast);
4822  rt_raster_destroy(newrast);
4823  if (NULL == pgrtn) {
4824  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4825  PG_RETURN_NULL();
4826  }
4827 
4828  SET_VARSIZE(pgrtn, pgrtn->size);
4829  PG_RETURN_POINTER(pgrtn);
4830  }
4831 
4832  RASTER_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new raster band...");
4833 
4834  /* Get the new raster band */
4835  newband = rt_raster_get_band(newrast, 0);
4836  if ( NULL == newband ) {
4837  elog(NOTICE, "Could not modify band for new raster. Returning new "
4838  "raster with the original band");
4839 
4840  if (initexpr)
4841  pfree(initexpr);
4842  rt_raster_destroy(raster);
4843  PG_FREE_IF_COPY(pgraster, 0);
4844 
4845  /* Serialize created raster */
4846  pgrtn = rt_raster_serialize(newrast);
4847  rt_raster_destroy(newrast);
4848  if (NULL == pgrtn) {
4849  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4850  PG_RETURN_NULL();
4851  }
4852 
4853  SET_VARSIZE(pgrtn, pgrtn->size);
4854  PG_RETURN_POINTER(pgrtn);
4855  }
4856 
4857  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Main computing loop (%d x %d)",
4858  width, height);
4859 
4860  if (initexpr != NULL) {
4861  /* Convert [rast.val] to [rast] */
4862  newexpr = rtpg_strreplace(initexpr, "[rast.val]", "[rast]", NULL);
4863  pfree(initexpr); initexpr=newexpr;
4864 
4865  sprintf(place,"$1");
4866  for (i = 0, j = 1; i < argkwcount; i++) {
4867  len = 0;
4868  newexpr = rtpg_strreplace(initexpr, argkw[i], place, &len);
4869  pfree(initexpr); initexpr=newexpr;
4870  if (len > 0) {
4871  argtype[argcount] = argkwtypes[i];
4872  argcount++;
4873  argpos[i] = j++;
4874 
4875  sprintf(place, "$%d", j);
4876  }
4877  else {
4878  argpos[i] = 0;
4879  }
4880  }
4881 
4882  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: initexpr = %s", initexpr);
4883 
4884  /* define values */
4885  values = (Datum *) palloc(sizeof(Datum) * argcount);
4886  if (values == NULL) {
4887 
4888  SPI_finish();
4889 
4890  rt_raster_destroy(raster);
4891  PG_FREE_IF_COPY(pgraster, 0);
4892  rt_raster_destroy(newrast);
4893 
4894  elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
4895  PG_RETURN_NULL();
4896  }
4897 
4898  /* define nulls */
4899  nulls = (char *)palloc(argcount);
4900  if (nulls == NULL) {
4901 
4902  SPI_finish();
4903 
4904  rt_raster_destroy(raster);
4905  PG_FREE_IF_COPY(pgraster, 0);
4906  rt_raster_destroy(newrast);
4907 
4908  elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
4909  PG_RETURN_NULL();
4910  }
4911 
4912  /* Connect to SPI and prepare the expression */
4913  ret = SPI_connect();
4914  if (ret != SPI_OK_CONNECT) {
4915 
4916  if (initexpr)
4917  pfree(initexpr);
4918  rt_raster_destroy(raster);
4919  PG_FREE_IF_COPY(pgraster, 0);
4920  rt_raster_destroy(newrast);
4921 
4922  elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4923  PG_RETURN_NULL();
4924  };
4925 
4926  /* Type of all arguments is FLOAT8OID */
4927  spi_plan = SPI_prepare(initexpr, argcount, argtype);
4928 
4929  if (spi_plan == NULL) {
4930 
4931  rt_raster_destroy(raster);
4932  PG_FREE_IF_COPY(pgraster, 0);
4933  rt_raster_destroy(newrast);
4934 
4935  SPI_finish();
4936 
4937  pfree(initexpr);
4938 
4939  elog(ERROR, "RASTER_mapAlgebraExpr: Could not prepare expression");
4940  PG_RETURN_NULL();
4941  }
4942  }
4943 
4944  for (x = 0; x < width; x++) {
4945  for(y = 0; y < height; y++) {
4946  ret = rt_band_get_pixel(band, x, y, &r, NULL);
4947 
4952  if (ret == ES_NONE && FLT_NEQ(r, newnodatavalue)) {
4953  if (skipcomputation == 0) {
4954  if (initexpr != NULL) {
4955  /* Reset the null arg flags. */
4956  memset(nulls, 'n', argcount);
4957 
4958  for (i = 0; i < argkwcount; i++) {
4959  idx = argpos[i];
4960  if (idx < 1) continue;
4961  idx--;
4962 
4963  if (i == kX ) {
4964  /* x is 0 based index, but SQL expects 1 based index */
4965  values[idx] = Int32GetDatum(x+1);
4966  nulls[idx] = ' ';
4967  }
4968  else if (i == kY) {
4969  /* y is 0 based index, but SQL expects 1 based index */
4970  values[idx] = Int32GetDatum(y+1);
4971  nulls[idx] = ' ';
4972  }
4973  else if (i == kVAL ) {
4974  values[idx] = Float8GetDatum(r);
4975  nulls[idx] = ' ';
4976  }
4977 
4978  }
4979 
4980  ret = SPI_execute_plan(spi_plan, values, nulls, FALSE, 0);
4981  if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
4982  SPI_processed != 1) {
4983  if (SPI_tuptable)
4984  SPI_freetuptable(tuptable);
4985 
4986  SPI_freeplan(spi_plan);
4987  SPI_finish();
4988 
4989  pfree(values);
4990  pfree(nulls);
4991  pfree(initexpr);
4992 
4993  rt_raster_destroy(raster);
4994  PG_FREE_IF_COPY(pgraster, 0);
4995  rt_raster_destroy(newrast);
4996 
4997  elog(ERROR, "RASTER_mapAlgebraExpr: Error executing prepared plan");
4998 
4999  PG_RETURN_NULL();
5000  }
5001 
5002  tupdesc = SPI_tuptable->tupdesc;
5003  tuptable = SPI_tuptable;
5004 
5005  tuple = tuptable->vals[0];
5006  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
5007  if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
5008  POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) errored, skip setting", x+1,y+1,r);
5009  newval = newinitialvalue;
5010  }
5011  else if ( isnull ) {
5012  POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting", x+1,y+1,r);
5013  newval = newinitialvalue;
5014  } else {
5015  newval = DatumGetFloat8(datum);
5016  }
5017 
5018  SPI_freetuptable(tuptable);
5019  }
5020 
5021  else
5022  newval = newinitialvalue;
5023 
5024  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new value = %f",
5025  newval);
5026  }
5027 
5028 
5029  rt_band_set_pixel(newband, x, y, newval, NULL);
5030  }
5031 
5032  }
5033  }
5034 
5035  if (initexpr != NULL) {
5036  SPI_freeplan(spi_plan);
5037  SPI_finish();
5038 
5039  pfree(values);
5040  pfree(nulls);
5041  pfree(initexpr);
5042  }
5043  else {
5044  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: no SPI cleanup");
5045  }
5046 
5047 
5048  /* The newrast band has been modified */
5049 
5050  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster modified, serializing it.");
5051  /* Serialize created raster */
5052 
5053  rt_raster_destroy(raster);
5054  PG_FREE_IF_COPY(pgraster, 0);
5055 
5056  pgrtn = rt_raster_serialize(newrast);
5057  rt_raster_destroy(newrast);
5058  if (NULL == pgrtn)
5059  PG_RETURN_NULL();
5060 
5061  SET_VARSIZE(pgrtn, pgrtn->size);
5062 
5063  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster serialized");
5064 
5065 
5066  POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraExpr: returning raster");
5067 
5068 
5069  PG_RETURN_POINTER(pgrtn);
5070 }
5071 
5072 /*
5073  * One raster user function MapAlgebra.
5074  */
5076 Datum RASTER_mapAlgebraFct(PG_FUNCTION_ARGS)
5077 {
5078  rt_pgraster *pgraster = NULL;
5079  rt_pgraster *pgrtn = NULL;
5080  rt_raster raster = NULL;
5081  rt_raster newrast = NULL;
5082  rt_band band = NULL;
5083  rt_band newband = NULL;
5084  int x, y, nband, width, height;
5085  double r;
5086  double newnodatavalue = 0.0;
5087  double newinitialvalue = 0.0;
5088  double newval = 0.0;
5089  rt_pixtype newpixeltype;
5090  int ret = -1;
5091  Oid oid;
5092  FmgrInfo cbinfo;
5093  FunctionCallInfoData cbdata;
5094  Datum tmpnewval;
5095  char * strFromText = NULL;
5096  int k = 0;
5097 
5098  POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFct: STARTING...");
5099 
5100  /* Check raster */
5101  if (PG_ARGISNULL(0)) {
5102  elog(WARNING, "Raster is NULL. Returning NULL");
5103  PG_RETURN_NULL();
5104  }
5105 
5106 
5107  /* Deserialize raster */
5108  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5109  raster = rt_raster_deserialize(pgraster, FALSE);
5110  if (NULL == raster) {
5111  PG_FREE_IF_COPY(pgraster, 0);
5112  elog(ERROR, "RASTER_mapAlgebraFct: Could not deserialize raster");
5113  PG_RETURN_NULL();
5114  }
5115 
5116  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting arguments...");
5117 
5118  /* Get the rest of the arguments */
5119 
5120  if (PG_ARGISNULL(1))
5121  nband = 1;
5122  else
5123  nband = PG_GETARG_INT32(1);
5124 
5125  if (nband < 1)
5126  nband = 1;
5127 
5128  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Creating new empty raster...");
5129 
5134  width = rt_raster_get_width(raster);
5135  height = rt_raster_get_height(raster);
5136 
5137  newrast = rt_raster_new(width, height);
5138 
5139  if ( NULL == newrast ) {
5140 
5141  rt_raster_destroy(raster);
5142  PG_FREE_IF_COPY(pgraster, 0);
5143 
5144  elog(ERROR, "RASTER_mapAlgebraFct: Could not create a new raster");
5145  PG_RETURN_NULL();
5146  }
5147 
5148  rt_raster_set_scale(newrast,
5149  rt_raster_get_x_scale(raster),
5150  rt_raster_get_y_scale(raster));
5151 
5152  rt_raster_set_offsets(newrast,
5153  rt_raster_get_x_offset(raster),
5154  rt_raster_get_y_offset(raster));
5155 
5156  rt_raster_set_skews(newrast,
5157  rt_raster_get_x_skew(raster),
5158  rt_raster_get_y_skew(raster));
5159 
5160  rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
5161 
5162 
5167  if (rt_raster_is_empty(newrast))
5168  {
5169  elog(NOTICE, "Raster is empty. Returning an empty raster");
5170  rt_raster_destroy(raster);
5171  PG_FREE_IF_COPY(pgraster, 0);
5172 
5173  pgrtn = rt_raster_serialize(newrast);
5174  rt_raster_destroy(newrast);
5175  if (NULL == pgrtn) {
5176  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5177  PG_RETURN_NULL();
5178  }
5179 
5180  SET_VARSIZE(pgrtn, pgrtn->size);
5181  PG_RETURN_POINTER(pgrtn);
5182  }
5183 
5184  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Getting raster band %d...", nband);
5185 
5190  if (!rt_raster_has_band(raster, nband - 1)) {
5191  elog(NOTICE, "Raster does not have the required band. Returning a raster "
5192  "without a band");
5193  rt_raster_destroy(raster);
5194  PG_FREE_IF_COPY(pgraster, 0);
5195 
5196  pgrtn = rt_raster_serialize(newrast);
5197  rt_raster_destroy(newrast);
5198  if (NULL == pgrtn) {
5199  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5200  PG_RETURN_NULL();
5201  }
5202 
5203  SET_VARSIZE(pgrtn, pgrtn->size);
5204  PG_RETURN_POINTER(pgrtn);
5205  }
5206 
5207  /* Get the raster band */
5208  band = rt_raster_get_band(raster, nband - 1);
5209  if ( NULL == band ) {
5210  elog(NOTICE, "Could not get the required band. Returning a raster "
5211  "without a band");
5212  rt_raster_destroy(raster);
5213  PG_FREE_IF_COPY(pgraster, 0);
5214 
5215  pgrtn = rt_raster_serialize(newrast);
5216  rt_raster_destroy(newrast);
5217  if (NULL == pgrtn) {
5218  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5219  PG_RETURN_NULL();
5220  }
5221 
5222  SET_VARSIZE(pgrtn, pgrtn->size);
5223  PG_RETURN_POINTER(pgrtn);
5224  }
5225 
5226  /*
5227  * Get NODATA value
5228  */
5229  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting NODATA value for band...");
5230 
5231  if (rt_band_get_hasnodata_flag(band)) {
5232  rt_band_get_nodata(band, &newnodatavalue);
5233  }
5234 
5235  else {
5236  newnodatavalue = rt_band_get_min_value(band);
5237  }
5238 
5239  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: NODATA value for band: %f",
5240  newnodatavalue);
5246  newinitialvalue = newnodatavalue;
5247 
5251  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Setting pixeltype...");
5252 
5253  if (PG_ARGISNULL(2)) {
5254  newpixeltype = rt_band_get_pixtype(band);
5255  }
5256 
5257  else {
5258  strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5259  newpixeltype = rt_pixtype_index_from_name(strFromText);
5260  pfree(strFromText);
5261  if (newpixeltype == PT_END)
5262  newpixeltype = rt_band_get_pixtype(band);
5263  }
5264