PostGIS  2.5.0dev-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@kbt.io>
13  * Copyright (C) 2013 Nathaniel Hunter Clay <clay.nathaniel@gmail.com>
14  *
15  * This program is free software; you can redistribute it and/or
16  * modify it under the terms of the GNU General Public License
17  * as published by the Free Software Foundation; either version 2
18  * of the License, or (at your option) any later version.
19  *
20  * This program is distributed in the hope that it will be useful,
21  * but WITHOUT ANY WARRANTY; without even the implied warranty of
22  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  * GNU General Public License for more details.
24  *
25  * You should have received a copy of the GNU General Public License
26  * along with this program; if not, write to the Free Software Foundation,
27  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28  *
29  */
30 
31 #include <assert.h>
32 
33 #include <postgres.h> /* for palloc */
34 #include <fmgr.h>
35 #include <funcapi.h>
36 #include <executor/spi.h>
37 #include <utils/lsyscache.h> /* for get_typlenbyvalalign */
38 #include <utils/array.h> /* for ArrayType */
39 #include <utils/builtins.h> /* for cstring_to_text */
40 #include <catalog/pg_type.h> /* for INT2OID, INT4OID, FLOAT4OID, FLOAT8OID and TEXTOID */
41 #include <executor/executor.h> /* for GetAttributeByName */
42 
43 #include "../../postgis_config.h"
44 #include "lwgeom_pg.h"
45 
46 #include "rtpostgis.h"
47 #include "rtpg_internal.h"
48 
49 /* n-raster MapAlgebra */
50 Datum RASTER_nMapAlgebra(PG_FUNCTION_ARGS);
51 Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS);
52 
53 /* raster union aggregate */
54 Datum RASTER_union_transfn(PG_FUNCTION_ARGS);
55 Datum RASTER_union_finalfn(PG_FUNCTION_ARGS);
56 
57 /* raster clip */
58 Datum RASTER_clip(PG_FUNCTION_ARGS);
59 
60 /* reclassify specified bands of a raster */
61 Datum RASTER_reclass(PG_FUNCTION_ARGS);
62 
63 /* apply colormap to specified band of a raster */
64 Datum RASTER_colorMap(PG_FUNCTION_ARGS);
65 
66 /* one-raster MapAlgebra */
67 Datum RASTER_mapAlgebraExpr(PG_FUNCTION_ARGS);
68 Datum RASTER_mapAlgebraFct(PG_FUNCTION_ARGS);
69 
70 /* one-raster neighborhood MapAlgebra */
71 Datum RASTER_mapAlgebraFctNgb(PG_FUNCTION_ARGS);
72 
73 /* two-raster MapAlgebra */
74 Datum RASTER_mapAlgebra2(PG_FUNCTION_ARGS);
75 
76 /* ---------------------------------------------------------------- */
77 /* n-raster MapAlgebra */
78 /* ---------------------------------------------------------------- */
79 
80 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  uint32_t x = 0;
373  uint32_t 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  InitFunctionCallInfoData(arg->callback.ufc_info, &(arg->callback.ufl_info), arg->callback.ufl_info.fn_nargs, InvalidOid, NULL, NULL);
840 
841  memset(arg->callback.ufc_info.argnull, FALSE, sizeof(bool) * arg->callback.ufl_info.fn_nargs);
842 
843  /* userargs (7) */
844  if (!PG_ARGISNULL(9))
845  arg->callback.ufc_info.arg[2] = PG_GETARG_DATUM(9);
846  else {
847  if (arg->callback.ufl_info.fn_strict) {
848  /* build and assign an empty TEXT array */
849  /* TODO: manually free the empty array? */
850  arg->callback.ufc_info.arg[2] = PointerGetDatum(
851  construct_empty_array(TEXTOID)
852  );
853  arg->callback.ufc_info.argnull[2] = FALSE;
854  }
855  else {
856  arg->callback.ufc_info.arg[2] = (Datum) NULL;
857  arg->callback.ufc_info.argnull[2] = TRUE;
858  }
859  }
860  }
861  else {
863  elog(ERROR, "RASTER_nMapAlgebra: callbackfunc must be provided");
864  PG_RETURN_NULL();
865  }
866 
867  /* determine nodataval and possibly pixtype */
868  /* band to check */
869  switch (arg->extenttype) {
870  case ET_LAST:
871  i = arg->numraster - 1;
872  break;
873  case ET_SECOND:
874  if (arg->numraster > 1) {
875  i = 1;
876  break;
877  }
878  default:
879  i = 0;
880  break;
881  }
882  /* find first viable band */
883  if (!arg->hasband[i]) {
884  for (i = 0; i < arg->numraster; i++) {
885  if (arg->hasband[i])
886  break;
887  }
888  if (i >= arg->numraster)
889  i = arg->numraster - 1;
890  }
891  band = rt_raster_get_band(arg->raster[i], arg->nband[i]);
892 
893  /* set pixel type if PT_END */
894  if (arg->pixtype == PT_END)
895  arg->pixtype = rt_band_get_pixtype(band);
896 
897  /* set hasnodata and nodataval */
898  arg->hasnodata = 1;
899  if (rt_band_get_hasnodata_flag(band))
900  rt_band_get_nodata(band, &(arg->nodataval));
901  else
902  arg->nodataval = rt_band_get_min_value(band);
903 
904  POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->pixtype), arg->hasnodata, arg->nodataval);
905 
906  /* init itrset */
907  itrset = palloc(sizeof(struct rt_iterator_t) * arg->numraster);
908  if (itrset == NULL) {
910  elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
911  PG_RETURN_NULL();
912  }
913 
914  /* set itrset */
915  for (i = 0; i < arg->numraster; i++) {
916  itrset[i].raster = arg->raster[i];
917  itrset[i].nband = arg->nband[i];
918  itrset[i].nbnodata = 1;
919  }
920 
921  /* pass everything to iterator */
922  noerr = rt_raster_iterator(
923  itrset, arg->numraster,
924  arg->extenttype, arg->cextent,
925  arg->pixtype,
926  arg->hasnodata, arg->nodataval,
927  arg->distance[0], arg->distance[1],
928  arg->mask,
929  &(arg->callback),
931  &raster
932  );
933 
934  /* cleanup */
935  pfree(itrset);
937 
938  if (noerr != ES_NONE) {
939  elog(ERROR, "RASTER_nMapAlgebra: Could not run raster iterator function");
940  PG_RETURN_NULL();
941  }
942  else if (raster == NULL)
943  PG_RETURN_NULL();
944 
945  pgraster = rt_raster_serialize(raster);
946  rt_raster_destroy(raster);
947 
948  POSTGIS_RT_DEBUG(3, "Finished");
949 
950  if (!pgraster)
951  PG_RETURN_NULL();
952 
953  SET_VARSIZE(pgraster, pgraster->size);
954  PG_RETURN_POINTER(pgraster);
955 }
956 
957 /* ---------------------------------------------------------------- */
958 /* expression ST_MapAlgebra for n rasters */
959 /* ---------------------------------------------------------------- */
960 
961 typedef struct {
963 
964  struct {
965  SPIPlanPtr spi_plan;
968 
969  int hasval;
970  double val;
971  } expr[3];
972 
973  struct {
974  int hasval;
975  double val;
976  } nodatanodata;
977 
978  struct {
979  int count;
980  char **val;
981  } kw;
982 
984 
987  rtpg_nmapalgebra_arg bandarg;
988 
990 };
991 
992 static rtpg_nmapalgebraexpr_arg rtpg_nmapalgebraexpr_arg_init(int cnt, char **kw) {
993  rtpg_nmapalgebraexpr_arg arg = NULL;
994  int i = 0;
995 
996  arg = palloc(sizeof(struct rtpg_nmapalgebraexpr_arg_t));
997  if (arg == NULL) {
998  elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arguments");
999  return NULL;
1000  }
1001 
1003  if (arg->bandarg == NULL) {
1004  elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arg->bandarg");
1005  return NULL;
1006  }
1007 
1008  arg->callback.kw.count = cnt;
1009  arg->callback.kw.val = kw;
1010 
1011  arg->callback.exprcount = 3;
1012  for (i = 0; i < arg->callback.exprcount; i++) {
1013  arg->callback.expr[i].spi_plan = NULL;
1014  arg->callback.expr[i].spi_argcount = 0;
1015  arg->callback.expr[i].spi_argpos = palloc(cnt * sizeof(uint8_t));
1016  if (arg->callback.expr[i].spi_argpos == NULL) {
1017  elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for spi_argpos");
1018  return NULL;
1019  }
1020  memset(arg->callback.expr[i].spi_argpos, 0, sizeof(uint8_t) * cnt);
1021  arg->callback.expr[i].hasval = 0;
1022  arg->callback.expr[i].val = 0;
1023  }
1024 
1025  arg->callback.nodatanodata.hasval = 0;
1026  arg->callback.nodatanodata.val = 0;
1027 
1028  return arg;
1029 }
1030 
1031 static void rtpg_nmapalgebraexpr_arg_destroy(rtpg_nmapalgebraexpr_arg arg) {
1032  int i = 0;
1033 
1035 
1036  for (i = 0; i < arg->callback.exprcount; i++) {
1037  if (arg->callback.expr[i].spi_plan)
1038  SPI_freeplan(arg->callback.expr[i].spi_plan);
1039  if (arg->callback.kw.count)
1040  pfree(arg->callback.expr[i].spi_argpos);
1041  }
1042 
1043  pfree(arg);
1044 }
1045 
1047  rt_iterator_arg arg, void *userarg,
1048  double *value, int *nodata
1049 ) {
1051  SPIPlanPtr plan = NULL;
1052  int i = 0;
1053  int id = -1;
1054 
1055  if (arg == NULL)
1056  return 0;
1057 
1058  *value = 0;
1059  *nodata = 0;
1060 
1061  /* 2 raster */
1062  if (arg->rasters > 1) {
1063  /* nodata1 = 1 AND nodata2 = 1, nodatanodataval */
1064  if (arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1065  if (callback->nodatanodata.hasval)
1066  *value = callback->nodatanodata.val;
1067  else
1068  *nodata = 1;
1069  }
1070  /* nodata1 = 1 AND nodata2 != 1, nodata1expr */
1071  else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0]) {
1072  id = 1;
1073  if (callback->expr[id].hasval)
1074  *value = callback->expr[id].val;
1075  else if (callback->expr[id].spi_plan)
1076  plan = callback->expr[id].spi_plan;
1077  else
1078  *nodata = 1;
1079  }
1080  /* nodata1 != 1 AND nodata2 = 1, nodata2expr */
1081  else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1082  id = 2;
1083  if (callback->expr[id].hasval)
1084  *value = callback->expr[id].val;
1085  else if (callback->expr[id].spi_plan)
1086  plan = callback->expr[id].spi_plan;
1087  else
1088  *nodata = 1;
1089  }
1090  /* expression */
1091  else {
1092  id = 0;
1093  if (callback->expr[id].hasval)
1094  *value = callback->expr[id].val;
1095  else if (callback->expr[id].spi_plan)
1096  plan = callback->expr[id].spi_plan;
1097  else {
1098  if (callback->nodatanodata.hasval)
1099  *value = callback->nodatanodata.val;
1100  else
1101  *nodata = 1;
1102  }
1103  }
1104  }
1105  /* 1 raster */
1106  else {
1107  /* nodata = 1, nodata1expr */
1108  if (arg->nodata[0][0][0]) {
1109  id = 1;
1110  if (callback->expr[id].hasval)
1111  *value = callback->expr[id].val;
1112  else if (callback->expr[id].spi_plan)
1113  plan = callback->expr[id].spi_plan;
1114  else
1115  *nodata = 1;
1116  }
1117  /* expression */
1118  else {
1119  id = 0;
1120  if (callback->expr[id].hasval)
1121  *value = callback->expr[id].val;
1122  else if (callback->expr[id].spi_plan)
1123  plan = callback->expr[id].spi_plan;
1124  else {
1125  /* see if nodata1expr is available */
1126  id = 1;
1127  if (callback->expr[id].hasval)
1128  *value = callback->expr[id].val;
1129  else if (callback->expr[id].spi_plan)
1130  plan = callback->expr[id].spi_plan;
1131  else
1132  *nodata = 1;
1133  }
1134  }
1135  }
1136 
1137  /* run prepared plan */
1138  if (plan != NULL) {
1139  Datum values[12];
1140  bool nulls[12];
1141  int err = 0;
1142 
1143  TupleDesc tupdesc;
1144  SPITupleTable *tuptable = NULL;
1145  HeapTuple tuple;
1146  Datum datum;
1147  bool isnull = FALSE;
1148 
1149  POSTGIS_RT_DEBUGF(4, "Running plan %d", id);
1150 
1151  /* init values and nulls */
1152  memset(values, (Datum) NULL, sizeof(Datum) * callback->kw.count);
1153  memset(nulls, FALSE, sizeof(bool) * callback->kw.count);
1154 
1155  if (callback->expr[id].spi_argcount) {
1156  int idx = 0;
1157 
1158  for (i = 0; i < callback->kw.count; i++) {
1159  idx = callback->expr[id].spi_argpos[i];
1160  if (idx < 1) continue;
1161  idx--; /* 1-based now 0-based */
1162 
1163  switch (i) {
1164  /* [rast.x] */
1165  case 0:
1166  values[idx] = Int32GetDatum(arg->src_pixel[0][0] + 1);
1167  break;
1168  /* [rast.y] */
1169  case 1:
1170  values[idx] = Int32GetDatum(arg->src_pixel[0][1] + 1);
1171  break;
1172  /* [rast.val] */
1173  case 2:
1174  /* [rast] */
1175  case 3:
1176  if (!arg->nodata[0][0][0])
1177  values[idx] = Float8GetDatum(arg->values[0][0][0]);
1178  else
1179  nulls[idx] = TRUE;
1180  break;
1181 
1182  /* [rast1.x] */
1183  case 4:
1184  values[idx] = Int32GetDatum(arg->src_pixel[0][0] + 1);
1185  break;
1186  /* [rast1.y] */
1187  case 5:
1188  values[idx] = Int32GetDatum(arg->src_pixel[0][1] + 1);
1189  break;
1190  /* [rast1.val] */
1191  case 6:
1192  /* [rast1] */
1193  case 7:
1194  if (!arg->nodata[0][0][0])
1195  values[idx] = Float8GetDatum(arg->values[0][0][0]);
1196  else
1197  nulls[idx] = TRUE;
1198  break;
1199 
1200  /* [rast2.x] */
1201  case 8:
1202  values[idx] = Int32GetDatum(arg->src_pixel[1][0] + 1);
1203  break;
1204  /* [rast2.y] */
1205  case 9:
1206  values[idx] = Int32GetDatum(arg->src_pixel[1][1] + 1);
1207  break;
1208  /* [rast2.val] */
1209  case 10:
1210  /* [rast2] */
1211  case 11:
1212  if (!arg->nodata[1][0][0])
1213  values[idx] = Float8GetDatum(arg->values[1][0][0]);
1214  else
1215  nulls[idx] = TRUE;
1216  break;
1217  }
1218 
1219  }
1220  }
1221 
1222  /* run prepared plan */
1223  err = SPI_execute_plan(plan, values, nulls, TRUE, 1);
1224  if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1225  elog(ERROR, "rtpg_nmapalgebraexpr_callback: Unexpected error when running prepared statement %d", id);
1226  return 0;
1227  }
1228 
1229  /* get output of prepared plan */
1230  tupdesc = SPI_tuptable->tupdesc;
1231  tuptable = SPI_tuptable;
1232  tuple = tuptable->vals[0];
1233 
1234  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1235  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1236  if (SPI_tuptable) SPI_freetuptable(tuptable);
1237  elog(ERROR, "rtpg_nmapalgebraexpr_callback: Could not get result of prepared statement %d", id);
1238  return 0;
1239  }
1240 
1241  if (!isnull) {
1242  *value = DatumGetFloat8(datum);
1243  POSTGIS_RT_DEBUG(4, "Getting value from Datum");
1244  }
1245  else {
1246  /* 2 raster, check nodatanodataval */
1247  if (arg->rasters > 1) {
1248  if (callback->nodatanodata.hasval)
1249  *value = callback->nodatanodata.val;
1250  else
1251  *nodata = 1;
1252  }
1253  /* 1 raster, check nodataval */
1254  else {
1255  if (callback->expr[1].hasval)
1256  *value = callback->expr[1].val;
1257  else
1258  *nodata = 1;
1259  }
1260  }
1261 
1262  if (SPI_tuptable) SPI_freetuptable(tuptable);
1263  }
1264 
1265  POSTGIS_RT_DEBUGF(4, "(value, nodata) = (%f, %d)", *value, *nodata);
1266  return 1;
1267 }
1268 
1270 Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS)
1271 {
1272  MemoryContext mainMemCtx = CurrentMemoryContext;
1273  rtpg_nmapalgebraexpr_arg arg = NULL;
1274  rt_iterator itrset;
1275  uint16_t exprpos[3] = {1, 4, 5};
1276 
1277  int i = 0;
1278  int j = 0;
1279  int k = 0;
1280 
1281  int numraster = 0;
1282  int err = 0;
1283  int allnull = 0;
1284  int allempty = 0;
1285  int noband = 0;
1286  int len = 0;
1287 
1288  TupleDesc tupdesc;
1289  SPITupleTable *tuptable = NULL;
1290  HeapTuple tuple;
1291  Datum datum;
1292  bool isnull = FALSE;
1293 
1294  rt_raster raster = NULL;
1295  rt_band band = NULL;
1296  rt_pgraster *pgraster = NULL;
1297 
1298  const int argkwcount = 12;
1299  char *argkw[] = {
1300  "[rast.x]",
1301  "[rast.y]",
1302  "[rast.val]",
1303  "[rast]",
1304  "[rast1.x]",
1305  "[rast1.y]",
1306  "[rast1.val]",
1307  "[rast1]",
1308  "[rast2.x]",
1309  "[rast2.y]",
1310  "[rast2.val]",
1311  "[rast2]"
1312  };
1313 
1314  if (PG_ARGISNULL(0))
1315  PG_RETURN_NULL();
1316 
1317  /* init argument struct */
1318  arg = rtpg_nmapalgebraexpr_arg_init(argkwcount, argkw);
1319  if (arg == NULL) {
1320  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not initialize argument structure");
1321  PG_RETURN_NULL();
1322  }
1323 
1324  /* let helper function process rastbandarg (0) */
1325  if (!rtpg_nmapalgebra_rastbandarg_process(arg->bandarg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
1327  elog(ERROR, "RASTER_nMapAlgebra: Could not process rastbandarg");
1328  PG_RETURN_NULL();
1329  }
1330 
1331  POSTGIS_RT_DEBUGF(4, "allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
1332 
1333  /* all rasters are NULL, return NULL */
1334  if (allnull == arg->bandarg->numraster) {
1335  elog(NOTICE, "All input rasters are NULL. Returning NULL");
1337  PG_RETURN_NULL();
1338  }
1339 
1340  /* only work on one or two rasters */
1341  if (arg->bandarg->numraster > 1)
1342  numraster = 2;
1343  else
1344  numraster = 1;
1345 
1346  /* pixel type (2) */
1347  if (!PG_ARGISNULL(2)) {
1348  char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
1349 
1350  /* Get the pixel type index */
1351  arg->bandarg->pixtype = rt_pixtype_index_from_name(pixtypename);
1352  if (arg->bandarg->pixtype == PT_END) {
1354  elog(ERROR, "RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
1355  PG_RETURN_NULL();
1356  }
1357  }
1358  POSTGIS_RT_DEBUGF(4, "pixeltype: %d", arg->bandarg->pixtype);
1359 
1360  /* extent type (3) */
1361  if (!PG_ARGISNULL(3)) {
1362  char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(3))));
1363  arg->bandarg->extenttype = rt_util_extent_type(extenttypename);
1364  }
1365 
1366  if (arg->bandarg->extenttype == ET_CUSTOM) {
1367  if (numraster < 2) {
1368  elog(NOTICE, "CUSTOM extent type not supported. Defaulting to FIRST");
1369  arg->bandarg->extenttype = ET_FIRST;
1370  }
1371  else {
1372  elog(NOTICE, "CUSTOM extent type not supported. Defaulting to INTERSECTION");
1374  }
1375  }
1376  else if (numraster < 2)
1377  arg->bandarg->extenttype = ET_FIRST;
1378 
1379  POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->bandarg->extenttype);
1380 
1381  /* nodatanodataval (6) */
1382  if (!PG_ARGISNULL(6)) {
1383  arg->callback.nodatanodata.hasval = 1;
1384  arg->callback.nodatanodata.val = PG_GETARG_FLOAT8(6);
1385  }
1386 
1387  err = 0;
1388  /* all rasters are empty, return empty raster */
1389  if (allempty == arg->bandarg->numraster) {
1390  elog(NOTICE, "All input rasters are empty. Returning empty raster");
1391  err = 1;
1392  }
1393  /* all rasters don't have indicated band, return empty raster */
1394  else if (noband == arg->bandarg->numraster) {
1395  elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
1396  err = 1;
1397  }
1398  if (err) {
1400 
1401  raster = rt_raster_new(0, 0);
1402  if (raster == NULL) {
1403  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create empty raster");
1404  PG_RETURN_NULL();
1405  }
1406 
1407  pgraster = rt_raster_serialize(raster);
1408  rt_raster_destroy(raster);
1409  if (!pgraster) PG_RETURN_NULL();
1410 
1411  SET_VARSIZE(pgraster, pgraster->size);
1412  PG_RETURN_POINTER(pgraster);
1413  }
1414 
1415  /* connect SPI */
1416  if (SPI_connect() != SPI_OK_CONNECT) {
1418  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
1419  PG_RETURN_NULL();
1420  }
1421 
1422  /*
1423  process expressions
1424 
1425  exprpos elements are:
1426  1 - expression => spi_plan[0]
1427  4 - nodata1expr => spi_plan[1]
1428  5 - nodata2expr => spi_plan[2]
1429  */
1430  for (i = 0; i < arg->callback.exprcount; i++) {
1431  char *expr = NULL;
1432  char *tmp = NULL;
1433  char *sql = NULL;
1434  char place[12] = "$1";
1435 
1436  if (PG_ARGISNULL(exprpos[i]))
1437  continue;
1438 
1439  expr = text_to_cstring(PG_GETARG_TEXT_P(exprpos[i]));
1440  POSTGIS_RT_DEBUGF(3, "raw expr of argument #%d: %s", exprpos[i], expr);
1441 
1442  for (j = 0, k = 1; j < argkwcount; j++) {
1443  /* attempt to replace keyword with placeholder */
1444  len = 0;
1445  tmp = rtpg_strreplace(expr, argkw[j], place, &len);
1446  pfree(expr);
1447  expr = tmp;
1448 
1449  if (len) {
1450  POSTGIS_RT_DEBUGF(4, "kw #%d (%s) at pos $%d", j, argkw[j], k);
1451  arg->callback.expr[i].spi_argcount++;
1452  arg->callback.expr[i].spi_argpos[j] = k++;
1453 
1454  sprintf(place, "$%d", k);
1455  }
1456  else
1457  arg->callback.expr[i].spi_argpos[j] = 0;
1458  }
1459 
1460  len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision");
1461  sql = (char *) palloc(len + 1);
1462  if (sql == NULL) {
1464  SPI_finish();
1465  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
1466  PG_RETURN_NULL();
1467  }
1468 
1469  strncpy(sql, "SELECT (", strlen("SELECT ("));
1470  strncpy(sql + strlen("SELECT ("), expr, strlen(expr));
1471  strncpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision"));
1472  sql[len] = '\0';
1473 
1474  POSTGIS_RT_DEBUGF(3, "sql #%d: %s", exprpos[i], sql);
1475 
1476  /* prepared plan */
1477  if (arg->callback.expr[i].spi_argcount) {
1478  Oid *argtype = (Oid *) palloc(arg->callback.expr[i].spi_argcount * sizeof(Oid));
1479  POSTGIS_RT_DEBUGF(3, "expression parameter %d is a prepared plan", exprpos[i]);
1480  if (argtype == NULL) {
1481  pfree(sql);
1483  SPI_finish();
1484  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
1485  PG_RETURN_NULL();
1486  }
1487 
1488  /* specify datatypes of parameters */
1489  for (j = 0, k = 0; j < argkwcount; j++) {
1490  if (arg->callback.expr[i].spi_argpos[j] < 1) continue;
1491 
1492  /* positions are INT4 */
1493  if (
1494  (strstr(argkw[j], "[rast.x]") != NULL) ||
1495  (strstr(argkw[j], "[rast.y]") != NULL) ||
1496  (strstr(argkw[j], "[rast1.x]") != NULL) ||
1497  (strstr(argkw[j], "[rast1.y]") != NULL) ||
1498  (strstr(argkw[j], "[rast2.x]") != NULL) ||
1499  (strstr(argkw[j], "[rast2.y]") != NULL)
1500  )
1501  argtype[k] = INT4OID;
1502  /* everything else is FLOAT8 */
1503  else
1504  argtype[k] = FLOAT8OID;
1505 
1506  k++;
1507  }
1508 
1509  arg->callback.expr[i].spi_plan = SPI_prepare(sql, arg->callback.expr[i].spi_argcount, argtype);
1510  pfree(argtype);
1511  pfree(sql);
1512 
1513  if (arg->callback.expr[i].spi_plan == NULL) {
1515  SPI_finish();
1516  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
1517  PG_RETURN_NULL();
1518  }
1519  }
1520  /* no args, just execute query */
1521  else {
1522  POSTGIS_RT_DEBUGF(3, "expression parameter %d has no args, simply executing", exprpos[i]);
1523  err = SPI_execute(sql, TRUE, 0);
1524  pfree(sql);
1525 
1526  if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1528  SPI_finish();
1529  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
1530  PG_RETURN_NULL();
1531  }
1532 
1533  /* get output of prepared plan */
1534  tupdesc = SPI_tuptable->tupdesc;
1535  tuptable = SPI_tuptable;
1536  tuple = tuptable->vals[0];
1537 
1538  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1539  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1540  if (SPI_tuptable) SPI_freetuptable(tuptable);
1542  SPI_finish();
1543  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
1544  PG_RETURN_NULL();
1545  }
1546 
1547  if (!isnull) {
1548  arg->callback.expr[i].hasval = 1;
1549  arg->callback.expr[i].val = DatumGetFloat8(datum);
1550  }
1551 
1552  if (SPI_tuptable) SPI_freetuptable(tuptable);
1553  }
1554  }
1555 
1556  /* determine nodataval and possibly pixtype */
1557  /* band to check */
1558  switch (arg->bandarg->extenttype) {
1559  case ET_LAST:
1560  case ET_SECOND:
1561  if (numraster > 1)
1562  i = 1;
1563  else
1564  i = 0;
1565  break;
1566  default:
1567  i = 0;
1568  break;
1569  }
1570  /* find first viable band */
1571  if (!arg->bandarg->hasband[i]) {
1572  for (i = 0; i < numraster; i++) {
1573  if (arg->bandarg->hasband[i])
1574  break;
1575  }
1576  if (i >= numraster)
1577  i = numraster - 1;
1578  }
1579  band = rt_raster_get_band(arg->bandarg->raster[i], arg->bandarg->nband[i]);
1580 
1581  /* set pixel type if PT_END */
1582  if (arg->bandarg->pixtype == PT_END)
1583  arg->bandarg->pixtype = rt_band_get_pixtype(band);
1584 
1585  /* set hasnodata and nodataval */
1586  arg->bandarg->hasnodata = 1;
1587  if (rt_band_get_hasnodata_flag(band))
1588  rt_band_get_nodata(band, &(arg->bandarg->nodataval));
1589  else
1590  arg->bandarg->nodataval = rt_band_get_min_value(band);
1591 
1592  POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->bandarg->pixtype), arg->bandarg->hasnodata, arg->bandarg->nodataval);
1593 
1594  /* init itrset */
1595  itrset = palloc(sizeof(struct rt_iterator_t) * numraster);
1596  if (itrset == NULL) {
1598  SPI_finish();
1599  elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
1600  PG_RETURN_NULL();
1601  }
1602 
1603  /* set itrset */
1604  for (i = 0; i < numraster; i++) {
1605  itrset[i].raster = arg->bandarg->raster[i];
1606  itrset[i].nband = arg->bandarg->nband[i];
1607  itrset[i].nbnodata = 1;
1608  }
1609 
1610  /* pass everything to iterator */
1611  err = rt_raster_iterator(
1612  itrset, numraster,
1613  arg->bandarg->extenttype, arg->bandarg->cextent,
1614  arg->bandarg->pixtype,
1615  arg->bandarg->hasnodata, arg->bandarg->nodataval,
1616  0, 0,
1617  NULL,
1618  &(arg->callback),
1620  &raster
1621  );
1622 
1623  pfree(itrset);
1625 
1626  if (err != ES_NONE) {
1627  SPI_finish();
1628  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not run raster iterator function");
1629  PG_RETURN_NULL();
1630  }
1631  else if (raster == NULL) {
1632  SPI_finish();
1633  PG_RETURN_NULL();
1634  }
1635 
1636  /* switch to prior memory context to ensure memory allocated in correct context */
1637  MemoryContextSwitchTo(mainMemCtx);
1638 
1639  pgraster = rt_raster_serialize(raster);
1640  rt_raster_destroy(raster);
1641 
1642  /* finish SPI */
1643  SPI_finish();
1644 
1645  if (!pgraster)
1646  PG_RETURN_NULL();
1647 
1648  SET_VARSIZE(pgraster, pgraster->size);
1649  PG_RETURN_POINTER(pgraster);
1650 }
1651 
1652 /* ---------------------------------------------------------------- */
1653 /* ST_Union aggregate functions */
1654 /* ---------------------------------------------------------------- */
1655 
1656 typedef enum {
1657  UT_LAST = 0,
1665 } rtpg_union_type;
1666 
1667 /* internal function translating text of UNION type to enum */
1669  assert(cutype && strlen(cutype) > 0);
1670 
1671  if (strcmp(cutype, "LAST") == 0)
1672  return UT_LAST;
1673  else if (strcmp(cutype, "FIRST") == 0)
1674  return UT_FIRST;
1675  else if (strcmp(cutype, "MIN") == 0)
1676  return UT_MIN;
1677  else if (strcmp(cutype, "MAX") == 0)
1678  return UT_MAX;
1679  else if (strcmp(cutype, "COUNT") == 0)
1680  return UT_COUNT;
1681  else if (strcmp(cutype, "SUM") == 0)
1682  return UT_SUM;
1683  else if (strcmp(cutype, "MEAN") == 0)
1684  return UT_MEAN;
1685  else if (strcmp(cutype, "RANGE") == 0)
1686  return UT_RANGE;
1687 
1688  return UT_LAST;
1689 }
1690 
1693  int nband; /* source raster's band index, 0-based */
1695 
1698 };
1699 
1702  int numband; /* number of bandargs */
1703  rtpg_union_band_arg bandarg;
1704 };
1705 
1706 static void rtpg_union_arg_destroy(rtpg_union_arg arg) {
1707  int i = 0;
1708  int j = 0;
1709  int k = 0;
1710 
1711  if (arg->bandarg != NULL) {
1712  for (i = 0; i < arg->numband; i++) {
1713  if (!arg->bandarg[i].numraster)
1714  continue;
1715 
1716  for (j = 0; j < arg->bandarg[i].numraster; j++) {
1717  if (arg->bandarg[i].raster[j] == NULL)
1718  continue;
1719 
1720  for (k = rt_raster_get_num_bands(arg->bandarg[i].raster[j]) - 1; k >= 0; k--)
1722  rt_raster_destroy(arg->bandarg[i].raster[j]);
1723  }
1724 
1725  pfree(arg->bandarg[i].raster);
1726  }
1727 
1728  pfree(arg->bandarg);
1729  }
1730 
1731  pfree(arg);
1732 }
1733 
1735  rt_iterator_arg arg, void *userarg,
1736  double *value, int *nodata
1737 ) {
1738  rtpg_union_type *utype = (rtpg_union_type *) userarg;
1739 
1740  if (arg == NULL)
1741  return 0;
1742 
1743  if (
1744  arg->rasters != 2 ||
1745  arg->rows != 1 ||
1746  arg->columns != 1
1747  ) {
1748  elog(ERROR, "rtpg_union_callback: Invalid arguments passed to callback");
1749  return 0;
1750  }
1751 
1752  *value = 0;
1753  *nodata = 0;
1754 
1755  /* handle NODATA situations except for COUNT, which is a special case */
1756  if (*utype != UT_COUNT) {
1757  /* both NODATA */
1758  if (arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1759  *nodata = 1;
1760  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1761  return 1;
1762  }
1763  /* second NODATA */
1764  else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1765  *value = arg->values[0][0][0];
1766  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1767  return 1;
1768  }
1769  /* first NODATA */
1770  else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0]) {
1771  *value = arg->values[1][0][0];
1772  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1773  return 1;
1774  }
1775  }
1776 
1777  switch (*utype) {
1778  case UT_FIRST:
1779  *value = arg->values[0][0][0];
1780  break;
1781  case UT_MIN:
1782  if (arg->values[0][0][0] < arg->values[1][0][0])
1783  *value = arg->values[0][0][0];
1784  else
1785  *value = arg->values[1][0][0];
1786  break;
1787  case UT_MAX:
1788  if (arg->values[0][0][0] > arg->values[1][0][0])
1789  *value = arg->values[0][0][0];
1790  else
1791  *value = arg->values[1][0][0];
1792  break;
1793  case UT_COUNT:
1794  /* both NODATA */
1795  if (arg->nodata[0][0][0] && arg->nodata[1][0][0])
1796  *value = 0;
1797  /* second NODATA */
1798  else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0])
1799  *value = arg->values[0][0][0];
1800  /* first NODATA */
1801  else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0])
1802  *value = 1;
1803  /* has value, increment */
1804  else
1805  *value = arg->values[0][0][0] + 1;
1806  break;
1807  case UT_SUM:
1808  *value = arg->values[0][0][0] + arg->values[1][0][0];
1809  break;
1810  case UT_MEAN:
1811  case UT_RANGE:
1812  break;
1813  case UT_LAST:
1814  default:
1815  *value = arg->values[1][0][0];
1816  break;
1817  }
1818 
1819  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1820 
1821 
1822  return 1;
1823 }
1824 
1826  rt_iterator_arg arg, void *userarg,
1827  double *value, int *nodata
1828 ) {
1829  if (arg == NULL)
1830  return 0;
1831 
1832  if (
1833  arg->rasters != 2 ||
1834  arg->rows != 1 ||
1835  arg->columns != 1
1836  ) {
1837  elog(ERROR, "rtpg_union_mean_callback: Invalid arguments passed to callback");
1838  return 0;
1839  }
1840 
1841  *value = 0;
1842  *nodata = 1;
1843 
1844  POSTGIS_RT_DEBUGF(4, "rast0: %f %d", arg->values[0][0][0], arg->nodata[0][0][0]);
1845  POSTGIS_RT_DEBUGF(4, "rast1: %f %d", arg->values[1][0][0], arg->nodata[1][0][0]);
1846 
1847  if (
1848  !arg->nodata[0][0][0] &&
1849  FLT_NEQ(arg->values[0][0][0], 0) &&
1850  !arg->nodata[1][0][0]
1851  ) {
1852  *value = arg->values[1][0][0] / arg->values[0][0][0];
1853  *nodata = 0;
1854  }
1855 
1856  POSTGIS_RT_DEBUGF(4, "value, nodata = (%f, %d)", *value, *nodata);
1857 
1858  return 1;
1859 }
1860 
1862  rt_iterator_arg arg, void *userarg,
1863  double *value, int *nodata
1864 ) {
1865  if (arg == NULL)
1866  return 0;
1867 
1868  if (
1869  arg->rasters != 2 ||
1870  arg->rows != 1 ||
1871  arg->columns != 1
1872  ) {
1873  elog(ERROR, "rtpg_union_range_callback: Invalid arguments passed to callback");
1874  return 0;
1875  }
1876 
1877  *value = 0;
1878  *nodata = 1;
1879 
1880  POSTGIS_RT_DEBUGF(4, "rast0: %f %d", arg->values[0][0][0], arg->nodata[0][0][0]);
1881  POSTGIS_RT_DEBUGF(4, "rast1: %f %d", arg->values[1][0][0], arg->nodata[1][0][0]);
1882 
1883  if (
1884  !arg->nodata[0][0][0] &&
1885  !arg->nodata[1][0][0]
1886  ) {
1887  *value = arg->values[1][0][0] - arg->values[0][0][0];
1888  *nodata = 0;
1889  }
1890 
1891  POSTGIS_RT_DEBUGF(4, "value, nodata = (%f, %d)", *value, *nodata);
1892 
1893  return 1;
1894 }
1895 
1896 /* called for ST_Union(raster, unionarg[]) */
1897 static int rtpg_union_unionarg_process(rtpg_union_arg arg, ArrayType *array) {
1898  Oid etype;
1899  Datum *e;
1900  bool *nulls;
1901  int16 typlen;
1902  bool typbyval;
1903  char typalign;
1904  int n = 0;
1905 
1906  HeapTupleHeader tup;
1907  bool isnull;
1908  Datum tupv;
1909 
1910  int i;
1911  int nband = 1;
1912  char *utypename = NULL;
1913  rtpg_union_type utype = UT_LAST;
1914 
1915  etype = ARR_ELEMTYPE(array);
1916  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1917 
1918  deconstruct_array(
1919  array,
1920  etype,
1921  typlen, typbyval, typalign,
1922  &e, &nulls, &n
1923  );
1924 
1925  if (!n) {
1926  elog(ERROR, "rtpg_union_unionarg_process: Invalid argument for unionarg");
1927  return 0;
1928  }
1929 
1930  /* prep arg */
1931  arg->numband = n;
1932  arg->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * arg->numband);
1933  if (arg->bandarg == NULL) {
1934  elog(ERROR, "rtpg_union_unionarg_process: Could not allocate memory for band information");
1935  return 0;
1936  }
1937 
1938  /* process each element */
1939  for (i = 0; i < n; i++) {
1940  if (nulls[i]) {
1941  arg->numband--;
1942  continue;
1943  }
1944 
1945  POSTGIS_RT_DEBUGF(4, "Processing unionarg at index %d", i);
1946 
1947  /* each element is a tuple */
1948  tup = (HeapTupleHeader) DatumGetPointer(e[i]);
1949  if (NULL == tup) {
1950  elog(ERROR, "rtpg_union_unionarg_process: Invalid argument for unionarg");
1951  return 0;
1952  }
1953 
1954  /* first element, bandnum */
1955  tupv = GetAttributeByName(tup, "nband", &isnull);
1956  if (isnull) {
1957  nband = i + 1;
1958  elog(NOTICE, "First argument (nband) of unionarg is NULL. Assuming nband = %d", nband);
1959  }
1960  else
1961  nband = DatumGetInt32(tupv);
1962 
1963  if (nband < 1) {
1964  elog(ERROR, "rtpg_union_unionarg_process: Band number must be greater than zero (1-based)");
1965  return 0;
1966  }
1967 
1968  /* second element, uniontype */
1969  tupv = GetAttributeByName(tup, "uniontype", &isnull);
1970  if (isnull) {
1971  elog(NOTICE, "Second argument (uniontype) of unionarg is NULL. Assuming uniontype = LAST");
1972  utype = UT_LAST;
1973  }
1974  else {
1975  utypename = text_to_cstring((text *) DatumGetPointer(tupv));
1976  utype = rtpg_uniontype_index_from_name(rtpg_strtoupper(utypename));
1977  }
1978 
1979  arg->bandarg[i].uniontype = utype;
1980  arg->bandarg[i].nband = nband - 1;
1981  arg->bandarg[i].raster = NULL;
1982 
1983  if (
1984  utype != UT_MEAN &&
1985  utype != UT_RANGE
1986  ) {
1987  arg->bandarg[i].numraster = 1;
1988  }
1989  else
1990  arg->bandarg[i].numraster = 2;
1991  }
1992 
1993  if (arg->numband < n) {
1994  arg->bandarg = repalloc(arg->bandarg, sizeof(struct rtpg_union_band_arg_t) * arg->numband);
1995  if (arg->bandarg == NULL) {
1996  elog(ERROR, "rtpg_union_unionarg_process: Could not reallocate memory for band information");
1997  return 0;
1998  }
1999  }
2000 
2001  return 1;
2002 }
2003 
2004 /* called for ST_Union(raster) */
2005 static int rtpg_union_noarg(rtpg_union_arg arg, rt_raster raster) {
2006  int numbands;
2007  int i;
2008 
2009  if (rt_raster_is_empty(raster))
2010  return 1;
2011 
2012  numbands = rt_raster_get_num_bands(raster);
2013  if (numbands <= arg->numband)
2014  return 1;
2015 
2016  /* more bands to process */
2017  POSTGIS_RT_DEBUG(4, "input raster has more bands, adding more bandargs");
2018  if (arg->numband)
2019  arg->bandarg = repalloc(arg->bandarg, sizeof(struct rtpg_union_band_arg_t) * numbands);
2020  else
2021  arg->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * numbands);
2022  if (arg->bandarg == NULL) {
2023  elog(ERROR, "rtpg_union_noarg: Could not reallocate memory for band information");
2024  return 0;
2025  }
2026 
2027  i = arg->numband;
2028  arg->numband = numbands;
2029  for (; i < arg->numband; i++) {
2030  POSTGIS_RT_DEBUGF(4, "Adding bandarg for band at index %d", i);
2031  arg->bandarg[i].uniontype = UT_LAST;
2032  arg->bandarg[i].nband = i;
2033  arg->bandarg[i].numraster = 1;
2034 
2035  arg->bandarg[i].raster = (rt_raster *) palloc(sizeof(rt_raster) * arg->bandarg[i].numraster);
2036  if (arg->bandarg[i].raster == NULL) {
2037  elog(ERROR, "rtpg_union_noarg: Could not allocate memory for working rasters");
2038  return 0;
2039  }
2040  memset(arg->bandarg[i].raster, 0, sizeof(rt_raster) * arg->bandarg[i].numraster);
2041 
2042  /* add new working rt_raster but only if working raster already exists */
2043  if (!rt_raster_is_empty(arg->bandarg[0].raster[0])) {
2044  arg->bandarg[i].raster[0] = rt_raster_clone(arg->bandarg[0].raster[0], 0); /* shallow clone */
2045  if (arg->bandarg[i].raster[0] == NULL) {
2046  elog(ERROR, "rtpg_union_noarg: Could not create working raster");
2047  return 0;
2048  }
2049  }
2050  }
2051 
2052  return 1;
2053 }
2054 
2055 /* UNION aggregate transition function */
2057 Datum RASTER_union_transfn(PG_FUNCTION_ARGS)
2058 {
2059  MemoryContext aggcontext;
2060  MemoryContext oldcontext;
2061  rtpg_union_arg iwr = NULL;
2062  int skiparg = 0;
2063 
2064  rt_pgraster *pgraster = NULL;
2065  rt_raster raster = NULL;
2066  rt_raster _raster = NULL;
2067  rt_band _band = NULL;
2068  int nband = 1;
2069  int noerr = 1;
2070  int isempty[2] = {0};
2071  int hasband[2] = {0};
2072  int nargs = 0;
2073  double _offset[4] = {0.};
2074  int nbnodata = 0; /* 1 if adding bands */
2075 
2076  int i = 0;
2077  int j = 0;
2078  int k = 0;
2079 
2080  rt_iterator itrset;
2081  char *utypename = NULL;
2082  rtpg_union_type utype = UT_LAST;
2083  rt_pixtype pixtype = PT_END;
2084  int hasnodata = 1;
2085  double nodataval = 0;
2086 
2087  rt_raster iraster = NULL;
2088  rt_band iband = NULL;
2089  int reuserast = 0;
2090  int y = 0;
2091  uint16_t _dim[2] = {0};
2092  void *vals = NULL;
2093  uint16_t nvals = 0;
2094 
2095  POSTGIS_RT_DEBUG(3, "Starting...");
2096 
2097  /* cannot be called directly as this is exclusive aggregate function */
2098  if (!AggCheckCallContext(fcinfo, &aggcontext)) {
2099  elog(ERROR, "RASTER_union_transfn: Cannot be called in a non-aggregate context");
2100  PG_RETURN_NULL();
2101  }
2102 
2103  /* switch to aggcontext */
2104  oldcontext = MemoryContextSwitchTo(aggcontext);
2105 
2106  if (PG_ARGISNULL(0)) {
2107  POSTGIS_RT_DEBUG(3, "Creating state variable");
2108  /* allocate container in aggcontext */
2109  iwr = palloc(sizeof(struct rtpg_union_arg_t));
2110  if (iwr == NULL) {
2111  MemoryContextSwitchTo(oldcontext);
2112  elog(ERROR, "RASTER_union_transfn: Could not allocate memory for state variable");
2113  PG_RETURN_NULL();
2114  }
2115 
2116  iwr->numband = 0;
2117  iwr->bandarg = NULL;
2118 
2119  skiparg = 0;
2120  }
2121  else {
2122  POSTGIS_RT_DEBUG(3, "State variable already exists");
2123  iwr = (rtpg_union_arg) PG_GETARG_POINTER(0);
2124  skiparg = 1;
2125  }
2126 
2127  /* raster arg is NOT NULL */
2128  if (!PG_ARGISNULL(1)) {
2129  /* deserialize raster */
2130  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2131 
2132  /* Get raster object */
2133  raster = rt_raster_deserialize(pgraster, FALSE);
2134  if (raster == NULL) {
2135 
2137  PG_FREE_IF_COPY(pgraster, 1);
2138 
2139  MemoryContextSwitchTo(oldcontext);
2140  elog(ERROR, "RASTER_union_transfn: Could not deserialize raster");
2141  PG_RETURN_NULL();
2142  }
2143  }
2144 
2145  /* process additional args if needed */
2146  nargs = PG_NARGS();
2147  POSTGIS_RT_DEBUGF(4, "nargs = %d", nargs);
2148  if (nargs > 2) {
2149  POSTGIS_RT_DEBUG(4, "processing additional arguments");
2150 
2151  /* if more than 2 arguments, determine the type of argument 3 */
2152  /* band number, UNION type or unionarg */
2153  if (!PG_ARGISNULL(2)) {
2154  Oid calltype = get_fn_expr_argtype(fcinfo->flinfo, 2);
2155 
2156  switch (calltype) {
2157  /* UNION type */
2158  case TEXTOID: {
2159  int idx = 0;
2160  int numband = 0;
2161 
2162  POSTGIS_RT_DEBUG(4, "Processing arg 3 as UNION type");
2163  nbnodata = 1;
2164 
2165  utypename = text_to_cstring(PG_GETARG_TEXT_P(2));
2166  utype = rtpg_uniontype_index_from_name(rtpg_strtoupper(utypename));
2167  POSTGIS_RT_DEBUGF(4, "Union type: %s", utypename);
2168 
2169  POSTGIS_RT_DEBUGF(4, "iwr->numband: %d", iwr->numband);
2170  /* see if we need to append new bands */
2171  if (raster) {
2172  idx = iwr->numband;
2173  numband = rt_raster_get_num_bands(raster);
2174  POSTGIS_RT_DEBUGF(4, "numband: %d", numband);
2175 
2176  /* only worry about appended bands */
2177  if (numband > iwr->numband)
2178  iwr->numband = numband;
2179  }
2180 
2181  if (!iwr->numband)
2182  iwr->numband = 1;
2183  POSTGIS_RT_DEBUGF(4, "iwr->numband: %d", iwr->numband);
2184  POSTGIS_RT_DEBUGF(4, "numband, idx: %d, %d", numband, idx);
2185 
2186  /* bandarg set. only possible after the first call to function */
2187  if (iwr->bandarg) {
2188  /* only reallocate if new bands need to be added */
2189  if (numband > idx) {
2190  POSTGIS_RT_DEBUG(4, "Reallocating iwr->bandarg");
2191  iwr->bandarg = repalloc(iwr->bandarg, sizeof(struct rtpg_union_band_arg_t) * iwr->numband);
2192  }
2193  /* prevent initial values step happening below */
2194  else
2195  idx = iwr->numband;
2196  }
2197  /* bandarg not set, first call to function */
2198  else {
2199  POSTGIS_RT_DEBUG(4, "Allocating iwr->bandarg");
2200  iwr->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * iwr->numband);
2201  }
2202  if (iwr->bandarg == NULL) {
2203 
2205  if (raster != NULL) {
2206  rt_raster_destroy(raster);
2207  PG_FREE_IF_COPY(pgraster, 1);
2208  }
2209 
2210  MemoryContextSwitchTo(oldcontext);
2211  elog(ERROR, "RASTER_union_transfn: Could not allocate memory for band information");
2212  PG_RETURN_NULL();
2213  }
2214 
2215  /* set initial values for bands that are "new" */
2216  for (i = idx; i < iwr->numband; i++) {
2217  iwr->bandarg[i].uniontype = utype;
2218  iwr->bandarg[i].nband = i;
2219 
2220  if (
2221  utype == UT_MEAN ||
2222  utype == UT_RANGE
2223  ) {
2224  iwr->bandarg[i].numraster = 2;
2225  }
2226  else
2227  iwr->bandarg[i].numraster = 1;
2228  iwr->bandarg[i].raster = NULL;
2229  }
2230 
2231  break;
2232  }
2233  /* band number */
2234  case INT2OID:
2235  case INT4OID:
2236  if (skiparg)
2237  break;
2238 
2239  POSTGIS_RT_DEBUG(4, "Processing arg 3 as band number");
2240  nband = PG_GETARG_INT32(2);
2241  if (nband < 1) {
2242 
2244  if (raster != NULL) {
2245  rt_raster_destroy(raster);
2246  PG_FREE_IF_COPY(pgraster, 1);
2247  }
2248 
2249  MemoryContextSwitchTo(oldcontext);
2250  elog(ERROR, "RASTER_union_transfn: Band number must be greater than zero (1-based)");
2251  PG_RETURN_NULL();
2252  }
2253 
2254  iwr->numband = 1;
2255  iwr->bandarg = palloc(sizeof(struct rtpg_union_band_arg_t) * iwr->numband);
2256  if (iwr->bandarg == NULL) {
2257 
2259  if (raster != NULL) {
2260  rt_raster_destroy(raster);
2261  PG_FREE_IF_COPY(pgraster, 1);
2262  }
2263 
2264  MemoryContextSwitchTo(oldcontext);
2265  elog(ERROR, "RASTER_union_transfn: Could not allocate memory for band information");
2266  PG_RETURN_NULL();
2267  }
2268 
2269  iwr->bandarg[0].uniontype = UT_LAST;
2270  iwr->bandarg[0].nband = nband - 1;
2271 
2272  iwr->bandarg[0].numraster = 1;
2273  iwr->bandarg[0].raster = NULL;
2274  break;
2275  /* only other type allowed is unionarg */
2276  default:
2277  if (skiparg)
2278  break;
2279 
2280  POSTGIS_RT_DEBUG(4, "Processing arg 3 as unionarg");
2281  if (!rtpg_union_unionarg_process(iwr, PG_GETARG_ARRAYTYPE_P(2))) {
2282 
2284  if (raster != NULL) {
2285  rt_raster_destroy(raster);
2286  PG_FREE_IF_COPY(pgraster, 1);
2287  }
2288 
2289  MemoryContextSwitchTo(oldcontext);
2290  elog(ERROR, "RASTER_union_transfn: Could not process unionarg");
2291  PG_RETURN_NULL();
2292  }
2293 
2294  break;
2295  }
2296  }
2297 
2298  /* UNION type */
2299  if (nargs > 3 && !PG_ARGISNULL(3)) {
2300  utypename = text_to_cstring(PG_GETARG_TEXT_P(3));
2301  utype = rtpg_uniontype_index_from_name(rtpg_strtoupper(utypename));
2302  iwr->bandarg[0].uniontype = utype;
2303  POSTGIS_RT_DEBUGF(4, "Union type: %s", utypename);
2304 
2305  if (
2306  utype == UT_MEAN ||
2307  utype == UT_RANGE
2308  ) {
2309  iwr->bandarg[0].numraster = 2;
2310  }
2311  }
2312 
2313  /* allocate space for pointers to rt_raster */
2314  for (i = 0; i < iwr->numband; i++) {
2315  POSTGIS_RT_DEBUGF(4, "iwr->bandarg[%d].raster @ %p", i, iwr->bandarg[i].raster);
2316 
2317  /* no need to allocate */
2318  if (iwr->bandarg[i].raster != NULL)
2319  continue;
2320 
2321  POSTGIS_RT_DEBUGF(4, "Allocating space for working rasters of band %d", i);
2322 
2323  iwr->bandarg[i].raster = (rt_raster *) palloc(sizeof(rt_raster) * iwr->bandarg[i].numraster);
2324  if (iwr->bandarg[i].raster == NULL) {
2325 
2327  if (raster != NULL) {
2328  rt_raster_destroy(raster);
2329  PG_FREE_IF_COPY(pgraster, 1);
2330  }
2331 
2332  MemoryContextSwitchTo(oldcontext);
2333  elog(ERROR, "RASTER_union_transfn: Could not allocate memory for working raster(s)");
2334  PG_RETURN_NULL();
2335  }
2336 
2337  memset(iwr->bandarg[i].raster, 0, sizeof(rt_raster) * iwr->bandarg[i].numraster);
2338 
2339  /* add new working rt_raster but only if working raster already exists */
2340  if (i > 0 && !rt_raster_is_empty(iwr->bandarg[0].raster[0])) {
2341  for (j = 0; j < iwr->bandarg[i].numraster; j++) {
2342  iwr->bandarg[i].raster[j] = rt_raster_clone(iwr->bandarg[0].raster[0], 0); /* shallow clone */
2343  if (iwr->bandarg[i].raster[j] == NULL) {
2344 
2346  if (raster != NULL) {
2347  rt_raster_destroy(raster);
2348  PG_FREE_IF_COPY(pgraster, 1);
2349  }
2350 
2351  MemoryContextSwitchTo(oldcontext);
2352  elog(ERROR, "RASTER_union_transfn: Could not create working raster");
2353  PG_RETURN_NULL();
2354  }
2355  }
2356  }
2357  }
2358  }
2359  /* only raster, no additional args */
2360  /* only do this if raster isn't empty */
2361  else {
2362  POSTGIS_RT_DEBUG(4, "no additional args, checking input raster");
2363  nbnodata = 1;
2364  if (!rtpg_union_noarg(iwr, raster)) {
2365 
2367  if (raster != NULL) {
2368  rt_raster_destroy(raster);
2369  PG_FREE_IF_COPY(pgraster, 1);
2370  }
2371 
2372  MemoryContextSwitchTo(oldcontext);
2373  elog(ERROR, "RASTER_union_transfn: Could not check and balance number of bands");
2374  PG_RETURN_NULL();
2375  }
2376  }
2377 
2378  /* init itrset */
2379  itrset = palloc(sizeof(struct rt_iterator_t) * 2);
2380  if (itrset == NULL) {
2381 
2383  if (raster != NULL) {
2384  rt_raster_destroy(raster);
2385  PG_FREE_IF_COPY(pgraster, 1);
2386  }
2387 
2388  MemoryContextSwitchTo(oldcontext);
2389  elog(ERROR, "RASTER_union_transfn: Could not allocate memory for iterator arguments");
2390  PG_RETURN_NULL();
2391  }
2392 
2393  /* by band to UNION */
2394  for (i = 0; i < iwr->numband; i++) {
2395 
2396  /* by raster */
2397  for (j = 0; j < iwr->bandarg[i].numraster; j++) {
2398  reuserast = 0;
2399 
2400  /* type of union */
2401  utype = iwr->bandarg[i].uniontype;
2402 
2403  /* raster flags */
2404  isempty[0] = rt_raster_is_empty(iwr->bandarg[i].raster[j]);
2405  isempty[1] = rt_raster_is_empty(raster);
2406 
2407  if (!isempty[0])
2408  hasband[0] = rt_raster_has_band(iwr->bandarg[i].raster[j], 0);
2409  if (!isempty[1])
2410  hasband[1] = rt_raster_has_band(raster, iwr->bandarg[i].nband);
2411 
2412  /* determine pixtype, hasnodata and nodataval */
2413  _band = NULL;
2414  if (!isempty[0] && hasband[0])
2415  _band = rt_raster_get_band(iwr->bandarg[i].raster[j], 0);
2416  else if (!isempty[1] && hasband[1])
2417  _band = rt_raster_get_band(raster, iwr->bandarg[i].nband);
2418  else {
2419  pixtype = PT_64BF;
2420  hasnodata = 1;
2421  nodataval = rt_pixtype_get_min_value(pixtype);
2422  }
2423  if (_band != NULL) {
2424  pixtype = rt_band_get_pixtype(_band);
2425  hasnodata = 1;
2426  if (rt_band_get_hasnodata_flag(_band))
2427  rt_band_get_nodata(_band, &nodataval);
2428  else
2429  nodataval = rt_band_get_min_value(_band);
2430  }
2431 
2432  /* UT_MEAN and UT_RANGE require two passes */
2433  /* UT_MEAN: first for UT_COUNT and second for UT_SUM */
2434  if (iwr->bandarg[i].uniontype == UT_MEAN) {
2435  /* first pass, UT_COUNT */
2436  if (j < 1)
2437  utype = UT_COUNT;
2438  else
2439  utype = UT_SUM;
2440  }
2441  /* UT_RANGE: first for UT_MIN and second for UT_MAX */
2442  else if (iwr->bandarg[i].uniontype == UT_RANGE) {
2443  /* first pass, UT_MIN */
2444  if (j < 1)
2445  utype = UT_MIN;
2446  else
2447  utype = UT_MAX;
2448  }
2449 
2450  /* force band settings for UT_COUNT */
2451  if (utype == UT_COUNT) {
2452  pixtype = PT_32BUI;
2453  hasnodata = 0;
2454  nodataval = 0;
2455  }
2456 
2457  POSTGIS_RT_DEBUGF(4, "(pixtype, hasnodata, nodataval) = (%s, %d, %f)", rt_pixtype_name(pixtype), hasnodata, nodataval);
2458 
2459  /* set itrset */
2460  itrset[0].raster = iwr->bandarg[i].raster[j];
2461  itrset[0].nband = 0;
2462  itrset[1].raster = raster;
2463  itrset[1].nband = iwr->bandarg[i].nband;
2464 
2465  /* allow use NODATA to replace missing bands */
2466  if (nbnodata) {
2467  itrset[0].nbnodata = 1;
2468  itrset[1].nbnodata = 1;
2469  }
2470  /* missing bands are not ignored */
2471  else {
2472  itrset[0].nbnodata = 0;
2473  itrset[1].nbnodata = 0;
2474  }
2475 
2476  /* if rasters AND bands are present, use copy approach */
2477  if (!isempty[0] && !isempty[1] && hasband[0] && hasband[1]) {
2478  POSTGIS_RT_DEBUG(3, "using line method");
2479 
2480  /* generate empty out raster */
2482  iwr->bandarg[i].raster[j], raster,
2483  ET_UNION,
2484  &iraster, _offset
2485  ) != ES_NONE) {
2486 
2487  pfree(itrset);
2489  if (raster != NULL) {
2490  rt_raster_destroy(raster);
2491  PG_FREE_IF_COPY(pgraster, 1);
2492  }
2493 
2494  MemoryContextSwitchTo(oldcontext);
2495  elog(ERROR, "RASTER_union_transfn: Could not create internal raster");
2496  PG_RETURN_NULL();
2497  }
2498  POSTGIS_RT_DEBUGF(4, "_offset = %f, %f, %f, %f",
2499  _offset[0], _offset[1], _offset[2], _offset[3]);
2500 
2501  /* rasters are spatially the same? */
2502  if (
2503  rt_raster_get_width(iwr->bandarg[i].raster[j]) == rt_raster_get_width(iraster) &&
2505  ) {
2506  double igt[6] = {0};
2507  double gt[6] = {0};
2508 
2510  rt_raster_get_geotransform_matrix(iraster, igt);
2511 
2512  reuserast = rt_util_same_geotransform_matrix(gt, igt);
2513  }
2514 
2515  /* use internal raster */
2516  if (!reuserast) {
2517  /* create band of same type */
2519  iraster,
2520  pixtype,
2521  nodataval,
2522  hasnodata, nodataval,
2523  0
2524  ) == -1) {
2525 
2526  pfree(itrset);
2528  rt_raster_destroy(iraster);
2529  if (raster != NULL) {
2530  rt_raster_destroy(raster);
2531  PG_FREE_IF_COPY(pgraster, 1);
2532  }
2533 
2534  MemoryContextSwitchTo(oldcontext);
2535  elog(ERROR, "RASTER_union_transfn: Could not add new band to internal raster");
2536  PG_RETURN_NULL();
2537  }
2538  iband = rt_raster_get_band(iraster, 0);
2539 
2540  /* copy working raster to output raster */
2541  _dim[0] = rt_raster_get_width(iwr->bandarg[i].raster[j]);
2542  _dim[1] = rt_raster_get_height(iwr->bandarg[i].raster[j]);
2543  for (y = 0; y < _dim[1]; y++) {
2544  POSTGIS_RT_DEBUGF(4, "Getting pixel line of working raster at (x, y, length) = (0, %d, %d)", y, _dim[0]);
2546  _band,
2547  0, y,
2548  _dim[0],
2549  &vals, &nvals
2550  ) != ES_NONE) {
2551 
2552  pfree(itrset);
2554  rt_band_destroy(iband);
2555  rt_raster_destroy(iraster);
2556  if (raster != NULL) {
2557  rt_raster_destroy(raster);
2558  PG_FREE_IF_COPY(pgraster, 1);
2559  }
2560 
2561  MemoryContextSwitchTo(oldcontext);
2562  elog(ERROR, "RASTER_union_transfn: Could not get pixel line from band of working raster");
2563  PG_RETURN_NULL();
2564  }
2565 
2566  POSTGIS_RT_DEBUGF(4, "Setting pixel line at (x, y, length) = (%d, %d, %d)", (int) _offset[0], (int) _offset[1] + y, nvals);
2568  iband,
2569  (int) _offset[0], (int) _offset[1] + y,
2570  vals, nvals
2571  ) != ES_NONE) {
2572 
2573  pfree(itrset);
2575  rt_band_destroy(iband);
2576  rt_raster_destroy(iraster);
2577  if (raster != NULL) {
2578  rt_raster_destroy(raster);
2579  PG_FREE_IF_COPY(pgraster, 1);
2580  }
2581 
2582  MemoryContextSwitchTo(oldcontext);
2583  elog(ERROR, "RASTER_union_transfn: Could not set pixel line to band of internal raster");
2584  PG_RETURN_NULL();
2585  }
2586  }
2587  }
2588  else {
2589  rt_raster_destroy(iraster);
2590  iraster = iwr->bandarg[i].raster[j];
2591  iband = rt_raster_get_band(iraster, 0);
2592  }
2593 
2594  /* run iterator for extent of input raster */
2595  noerr = rt_raster_iterator(
2596  itrset, 2,
2597  ET_LAST, NULL,
2598  pixtype,
2599  hasnodata, nodataval,
2600  0, 0,
2601  NULL,
2602  &utype,
2604  &_raster
2605  );
2606  if (noerr != ES_NONE) {
2607 
2608  pfree(itrset);
2610  if (!reuserast) {
2611  rt_band_destroy(iband);
2612  rt_raster_destroy(iraster);
2613  }
2614  if (raster != NULL) {
2615  rt_raster_destroy(raster);
2616  PG_FREE_IF_COPY(pgraster, 1);
2617  }
2618 
2619  MemoryContextSwitchTo(oldcontext);
2620  elog(ERROR, "RASTER_union_transfn: Could not run raster iterator function");
2621  PG_RETURN_NULL();
2622  }
2623 
2624  /* with iterator raster, copy data to output raster */
2625  _band = rt_raster_get_band(_raster, 0);
2626  _dim[0] = rt_raster_get_width(_raster);
2627  _dim[1] = rt_raster_get_height(_raster);
2628  for (y = 0; y < _dim[1]; y++) {
2629  POSTGIS_RT_DEBUGF(4, "Getting pixel line of iterator raster at (x, y, length) = (0, %d, %d)", y, _dim[0]);
2631  _band,
2632  0, y,
2633  _dim[0],
2634  &vals, &nvals
2635  ) != ES_NONE) {
2636 
2637  pfree(itrset);
2639  if (!reuserast) {
2640  rt_band_destroy(iband);
2641  rt_raster_destroy(iraster);
2642  }
2643  if (raster != NULL) {
2644  rt_raster_destroy(raster);
2645  PG_FREE_IF_COPY(pgraster, 1);
2646  }
2647 
2648  MemoryContextSwitchTo(oldcontext);
2649  elog(ERROR, "RASTER_union_transfn: Could not get pixel line from band of working raster");
2650  PG_RETURN_NULL();
2651  }
2652 
2653  POSTGIS_RT_DEBUGF(4, "Setting pixel line at (x, y, length) = (%d, %d, %d)", (int) _offset[2], (int) _offset[3] + y, nvals);
2655  iband,
2656  (int) _offset[2], (int) _offset[3] + y,
2657  vals, nvals
2658  ) != ES_NONE) {
2659 
2660  pfree(itrset);
2662  if (!reuserast) {
2663  rt_band_destroy(iband);
2664  rt_raster_destroy(iraster);
2665  }
2666  if (raster != NULL) {
2667  rt_raster_destroy(raster);
2668  PG_FREE_IF_COPY(pgraster, 1);
2669  }
2670 
2671  MemoryContextSwitchTo(oldcontext);
2672  elog(ERROR, "RASTER_union_transfn: Could not set pixel line to band of internal raster");
2673  PG_RETURN_NULL();
2674  }
2675  }
2676 
2677  /* free _raster */
2678  rt_band_destroy(_band);
2679  rt_raster_destroy(_raster);
2680 
2681  /* replace working raster with output raster */
2682  _raster = iraster;
2683  }
2684  else {
2685  POSTGIS_RT_DEBUG(3, "using pixel method");
2686 
2687  /* pass everything to iterator */
2688  noerr = rt_raster_iterator(
2689  itrset, 2,
2690  ET_UNION, NULL,
2691  pixtype,
2692  hasnodata, nodataval,
2693  0, 0,
2694  NULL,
2695  &utype,
2697  &_raster
2698  );
2699 
2700  if (noerr != ES_NONE) {
2701 
2702  pfree(itrset);
2704  if (raster != NULL) {
2705  rt_raster_destroy(raster);
2706  PG_FREE_IF_COPY(pgraster, 1);
2707  }
2708 
2709  MemoryContextSwitchTo(oldcontext);
2710  elog(ERROR, "RASTER_union_transfn: Could not run raster iterator function");
2711  PG_RETURN_NULL();
2712  }
2713  }
2714 
2715  /* replace working raster */
2716  if (iwr->bandarg[i].raster[j] != NULL && !reuserast) {
2717  for (k = rt_raster_get_num_bands(iwr->bandarg[i].raster[j]) - 1; k >= 0; k--)
2719  rt_raster_destroy(iwr->bandarg[i].raster[j]);
2720  }
2721  iwr->bandarg[i].raster[j] = _raster;
2722  }
2723 
2724  }
2725 
2726  pfree(itrset);
2727  if (raster != NULL) {
2728  rt_raster_destroy(raster);
2729  PG_FREE_IF_COPY(pgraster, 1);
2730  }
2731 
2732  /* switch back to local context */
2733  MemoryContextSwitchTo(oldcontext);
2734 
2735  POSTGIS_RT_DEBUG(3, "Finished");
2736 
2737  PG_RETURN_POINTER(iwr);
2738 }
2739 
2740 /* UNION aggregate final function */
2742 Datum RASTER_union_finalfn(PG_FUNCTION_ARGS)
2743 {
2744  rtpg_union_arg iwr;
2745  rt_raster _rtn = NULL;
2746  rt_raster _raster = NULL;
2747  rt_pgraster *pgraster = NULL;
2748 
2749  int i = 0;
2750  int j = 0;
2751  rt_iterator itrset = NULL;
2752  rt_band _band = NULL;
2753  int noerr = 1;
2754  int status = 0;
2755  rt_pixtype pixtype = PT_END;
2756  int hasnodata = 0;
2757  double nodataval = 0;
2758 
2759  POSTGIS_RT_DEBUG(3, "Starting...");
2760 
2761  /* cannot be called directly as this is exclusive aggregate function */
2762  if (!AggCheckCallContext(fcinfo, NULL)) {
2763  elog(ERROR, "RASTER_union_finalfn: Cannot be called in a non-aggregate context");
2764  PG_RETURN_NULL();
2765  }
2766 
2767  /* NULL, return null */
2768  if (PG_ARGISNULL(0))
2769  PG_RETURN_NULL();
2770 
2771  iwr = (rtpg_union_arg) PG_GETARG_POINTER(0);
2772 
2773  /* init itrset */
2774  itrset = palloc(sizeof(struct rt_iterator_t) * 2);
2775  if (itrset == NULL) {
2777  elog(ERROR, "RASTER_union_finalfn: Could not allocate memory for iterator arguments");
2778  PG_RETURN_NULL();
2779  }
2780 
2781  for (i = 0; i < iwr->numband; i++) {
2782  if (
2783  iwr->bandarg[i].uniontype == UT_MEAN ||
2784  iwr->bandarg[i].uniontype == UT_RANGE
2785  ) {
2786  /* raster containing the SUM or MAX is at index 1 */
2787  _band = rt_raster_get_band(iwr->bandarg[i].raster[1], 0);
2788 
2789  pixtype = rt_band_get_pixtype(_band);
2790  hasnodata = rt_band_get_hasnodata_flag(_band);
2791  if (hasnodata)
2792  rt_band_get_nodata(_band, &nodataval);
2793  POSTGIS_RT_DEBUGF(4, "(pixtype, hasnodata, nodataval) = (%s, %d, %f)", rt_pixtype_name(pixtype), hasnodata, nodataval);
2794 
2795  itrset[0].raster = iwr->bandarg[i].raster[0];
2796  itrset[0].nband = 0;
2797  itrset[1].raster = iwr->bandarg[i].raster[1];
2798  itrset[1].nband = 0;
2799 
2800  /* pass everything to iterator */
2801  if (iwr->bandarg[i].uniontype == UT_MEAN) {
2802  noerr = rt_raster_iterator(
2803  itrset, 2,
2804  ET_UNION, NULL,
2805  pixtype,
2806  hasnodata, nodataval,
2807  0, 0,
2808  NULL,
2809  NULL,
2811  &_raster
2812  );
2813  }
2814  else if (iwr->bandarg[i].uniontype == UT_RANGE) {
2815  noerr = rt_raster_iterator(
2816  itrset, 2,
2817  ET_UNION, NULL,
2818  pixtype,
2819  hasnodata, nodataval,
2820  0, 0,
2821  NULL,
2822  NULL,
2824  &_raster
2825  );
2826  }
2827 
2828  if (noerr != ES_NONE) {
2829  pfree(itrset);
2831  if (_rtn != NULL)
2832  rt_raster_destroy(_rtn);
2833  elog(ERROR, "RASTER_union_finalfn: Could not run raster iterator function");
2834  PG_RETURN_NULL();
2835  }
2836  }
2837  else
2838  _raster = iwr->bandarg[i].raster[0];
2839 
2840  /* first band, _rtn doesn't exist */
2841  if (i < 1) {
2842  uint32_t bandNums[1] = {0};
2843  _rtn = rt_raster_from_band(_raster, bandNums, 1);
2844  status = (_rtn == NULL) ? -1 : 0;
2845  }
2846  else
2847  status = rt_raster_copy_band(_rtn, _raster, 0, i);
2848 
2849  POSTGIS_RT_DEBUG(4, "destroying source rasters");
2850 
2851  /* destroy source rasters */
2852  if (
2853  iwr->bandarg[i].uniontype == UT_MEAN ||
2854  iwr->bandarg[i].uniontype == UT_RANGE
2855  ) {
2856  rt_raster_destroy(_raster);
2857  }
2858 
2859  for (j = 0; j < iwr->bandarg[i].numraster; j++) {
2860  if (iwr->bandarg[i].raster[j] == NULL)
2861  continue;
2862  rt_raster_destroy(iwr->bandarg[i].raster[j]);
2863  iwr->bandarg[i].raster[j] = NULL;
2864  }
2865 
2866  if (status < 0) {
2868  rt_raster_destroy(_rtn);
2869  elog(ERROR, "RASTER_union_finalfn: Could not add band to final raster");
2870  PG_RETURN_NULL();
2871  }
2872  }
2873 
2874  /* cleanup */
2875  pfree(itrset);
2877 
2878  if (!_rtn) PG_RETURN_NULL();
2879 
2880  pgraster = rt_raster_serialize(_rtn);
2881  rt_raster_destroy(_rtn);
2882 
2883  POSTGIS_RT_DEBUG(3, "Finished");
2884 
2885  if (!pgraster)
2886  PG_RETURN_NULL();
2887 
2888  SET_VARSIZE(pgraster, pgraster->size);
2889  PG_RETURN_POINTER(pgraster);
2890 }
2891 
2892 /* ---------------------------------------------------------------- */
2893 /* Clip raster with geometry */
2894 /* ---------------------------------------------------------------- */
2895 
2898  int nband; /* band index */
2899  int hasnodata; /* is there a user-specified NODATA? */
2900  double nodataval; /* user-specified NODATA */
2901 };
2902 
2908 
2909  int numbands; /* number of bandargs */
2910  rtpg_clip_band band;
2911 };
2912 
2913 static rtpg_clip_arg rtpg_clip_arg_init() {
2914  rtpg_clip_arg arg = NULL;
2915 
2916  arg = palloc(sizeof(struct rtpg_clip_arg_t));
2917  if (arg == NULL) {
2918  elog(ERROR, "rtpg_clip_arg_init: Could not allocate memory for function arguments");
2919  return NULL;
2920  }
2921 
2922  arg->extenttype = ET_INTERSECTION;
2923  arg->raster = NULL;
2924  arg->mask = NULL;
2925  arg->numbands = 0;
2926  arg->band = NULL;
2927 
2928  return arg;
2929 }
2930 
2931 static void rtpg_clip_arg_destroy(rtpg_clip_arg arg) {
2932  if (arg->band != NULL)
2933  pfree(arg->band);
2934 
2935  if (arg->raster != NULL)
2936  rt_raster_destroy(arg->raster);
2937  if (arg->mask != NULL)
2938  rt_raster_destroy(arg->mask);
2939 
2940  pfree(arg);
2941 }
2942 
2944  rt_iterator_arg arg, void *userarg,
2945  double *value, int *nodata
2946 ) {
2947  *value = 0;
2948  *nodata = 0;
2949 
2950  /* either is NODATA, output is NODATA */
2951  if (arg->nodata[0][0][0] || arg->nodata[1][0][0])
2952  *nodata = 1;
2953  /* set to value */
2954  else
2955  *value = arg->values[0][0][0];
2956 
2957  return 1;
2958 }
2959 
2961 Datum RASTER_clip(PG_FUNCTION_ARGS)
2962 {
2963  rt_pgraster *pgraster = NULL;
2964  LWGEOM *rastgeom = NULL;
2965  double gt[6] = {0};
2966  int srid = SRID_UNKNOWN;
2967 
2968  rt_pgraster *pgrtn = NULL;
2969  rt_raster rtn = NULL;
2970 
2971  GSERIALIZED *gser = NULL;
2972  LWGEOM *geom = NULL;
2973  unsigned char *wkb = NULL;
2974  size_t wkb_len;
2975 
2976  ArrayType *array;
2977  Oid etype;
2978  Datum *e;
2979  bool *nulls;
2980 
2981  int16 typlen;
2982  bool typbyval;
2983  char typalign;
2984 
2985  int i = 0;
2986  int j = 0;
2987  int k = 0;
2988  rtpg_clip_arg arg = NULL;
2989  LWGEOM *tmpgeom = NULL;
2990  rt_iterator itrset;
2991 
2992  rt_raster _raster = NULL;
2993  rt_band band = NULL;
2994  rt_pixtype pixtype;
2995  int hasnodata;
2996  double nodataval;
2997  int noerr = 0;
2998 
2999  POSTGIS_RT_DEBUG(3, "Starting...");
3000 
3001  /* raster or geometry is NULL, return NULL */
3002  if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3003  PG_RETURN_NULL();
3004 
3005  /* init arg */
3006  arg = rtpg_clip_arg_init();
3007  if (arg == NULL) {
3008  elog(ERROR, "RASTER_clip: Could not initialize argument structure");
3009  PG_RETURN_NULL();
3010  }
3011 
3012  /* raster (0) */
3013  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3014 
3015  /* Get raster object */
3016  arg->raster = rt_raster_deserialize(pgraster, FALSE);
3017  if (arg->raster == NULL) {
3018  rtpg_clip_arg_destroy(arg);
3019  PG_FREE_IF_COPY(pgraster, 0);
3020  elog(ERROR, "RASTER_clip: Could not deserialize raster");
3021  PG_RETURN_NULL();
3022  }
3023 
3024  /* raster is empty, return empty raster */
3025  if (rt_raster_is_empty(arg->raster) || rt_raster_get_num_bands(arg->raster) == 0) {
3026  elog(NOTICE, "Input raster is empty or has no bands. Returning empty raster");
3027 
3028  rtpg_clip_arg_destroy(arg);
3029  PG_FREE_IF_COPY(pgraster, 0);
3030 
3031  rtn = rt_raster_new(0, 0);
3032  if (rtn == NULL) {
3033  elog(ERROR, "RASTER_clip: Could not create empty raster");
3034  PG_RETURN_NULL();
3035  }
3036 
3037  pgrtn = rt_raster_serialize(rtn);
3038  rt_raster_destroy(rtn);
3039  if (NULL == pgrtn)
3040  PG_RETURN_NULL();
3041 
3042  SET_VARSIZE(pgrtn, pgrtn->size);
3043  PG_RETURN_POINTER(pgrtn);
3044  }
3045 
3046  /* metadata */
3048  srid = clamp_srid(rt_raster_get_srid(arg->raster));
3049 
3050  /* geometry (2) */
3051  gser = PG_GETARG_GSERIALIZED_P(2);
3052  geom = lwgeom_from_gserialized(gser);
3053 
3054  /* Get a 2D version of the geometry if necessary */
3055  if (lwgeom_ndims(geom) > 2) {
3056  LWGEOM *geom2d = lwgeom_force_2d(geom);
3057  lwgeom_free(geom);
3058  geom = geom2d;
3059  }
3060 
3061  /* check that SRIDs match */
3062  if (srid != clamp_srid(gserialized_get_srid(gser))) {
3063  elog(NOTICE, "Geometry provided does not have the same SRID as the raster. Returning NULL");
3064 
3065  rtpg_clip_arg_destroy(arg);
3066  PG_FREE_IF_COPY(pgraster, 0);
3067  lwgeom_free(geom);
3068  PG_FREE_IF_COPY(gser, 2);
3069 
3070  PG_RETURN_NULL();
3071  }
3072 
3073  /* crop (4) */
3074  if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3075  arg->extenttype = ET_FIRST;
3076 
3077  /* get intersection geometry of input raster and input geometry */
3078  if (rt_raster_get_convex_hull(arg->raster, &rastgeom) != ES_NONE) {
3079 
3080  rtpg_clip_arg_destroy(arg);
3081  PG_FREE_IF_COPY(pgraster, 0);
3082  lwgeom_free(geom);
3083  PG_FREE_IF_COPY(gser, 2);
3084 
3085  elog(ERROR, "RASTER_clip: Could not get convex hull of raster");
3086  PG_RETURN_NULL();
3087  }
3088 
3089  tmpgeom = lwgeom_intersection(rastgeom, geom);
3090  lwgeom_free(rastgeom);
3091  lwgeom_free(geom);
3092  PG_FREE_IF_COPY(gser, 2);
3093  geom = tmpgeom;
3094 
3095  /* intersection is empty AND extent type is INTERSECTION, return empty */
3096  if (lwgeom_is_empty(geom) && arg->extenttype == ET_INTERSECTION) {
3097  elog(NOTICE, "The input raster and input geometry do not intersect. Returning empty raster");
3098 
3099  rtpg_clip_arg_destroy(arg);
3100  PG_FREE_IF_COPY(pgraster, 0);
3101  lwgeom_free(geom);
3102 
3103  rtn = rt_raster_new(0, 0);
3104  if (rtn == NULL) {
3105  elog(ERROR, "RASTER_clip: Could not create empty raster");
3106  PG_RETURN_NULL();
3107  }
3108 
3109  pgrtn = rt_raster_serialize(rtn);
3110  rt_raster_destroy(rtn);
3111  if (NULL == pgrtn)
3112  PG_RETURN_NULL();
3113 
3114  SET_VARSIZE(pgrtn, pgrtn->size);
3115  PG_RETURN_POINTER(pgrtn);
3116  }
3117 
3118  /* nband (1) */
3119  if (!PG_ARGISNULL(1)) {
3120  array = PG_GETARG_ARRAYTYPE_P(1);
3121  etype = ARR_ELEMTYPE(array);
3122  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3123 
3124  switch (etype) {
3125  case INT2OID:
3126  case INT4OID:
3127  break;
3128  default:
3129  rtpg_clip_arg_destroy(arg);
3130  PG_FREE_IF_COPY(pgraster, 0);
3131  lwgeom_free(geom);
3132  elog(ERROR, "RASTER_clip: Invalid data type for band indexes");
3133  PG_RETURN_NULL();
3134  break;
3135  }
3136 
3137  deconstruct_array(
3138  array, etype,
3139  typlen, typbyval, typalign,
3140  &e, &nulls, &(arg->numbands)
3141  );
3142 
3143  arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3144  if (arg->band == NULL) {
3145  rtpg_clip_arg_destroy(arg);
3146  PG_FREE_IF_COPY(pgraster, 0);
3147  lwgeom_free(geom);
3148  elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3149  PG_RETURN_NULL();
3150  }
3151 
3152  for (i = 0, j = 0; i < arg->numbands; i++) {
3153  if (nulls[i]) continue;
3154 
3155  switch (etype) {
3156  case INT2OID:
3157  arg->band[j].nband = DatumGetInt16(e[i]) - 1;
3158  break;
3159  case INT4OID:
3160  arg->band[j].nband = DatumGetInt32(e[i]) - 1;
3161  break;
3162  }
3163 
3164  j++;
3165  }
3166 
3167  if (j < arg->numbands) {
3168  arg->band = repalloc(arg->band, sizeof(struct rtpg_clip_band_t) * j);
3169  if (arg->band == NULL) {
3170  rtpg_clip_arg_destroy(arg);
3171  PG_FREE_IF_COPY(pgraster, 0);
3172  lwgeom_free(geom);
3173  elog(ERROR, "RASTER_clip: Could not reallocate memory for band arguments");
3174  PG_RETURN_NULL();
3175  }
3176 
3177  arg->numbands = j;
3178  }
3179 
3180  /* validate band */
3181  for (i = 0; i < arg->numbands; i++) {
3182  if (!rt_raster_has_band(arg->raster, arg->band[i].nband)) {
3183  elog(NOTICE, "Band at index %d not found in raster", arg->band[i].nband + 1);
3184  rtpg_clip_arg_destroy(arg);
3185  PG_FREE_IF_COPY(pgraster, 0);
3186  lwgeom_free(geom);
3187  PG_RETURN_NULL();
3188  }
3189 
3190  arg->band[i].hasnodata = 0;
3191  arg->band[i].nodataval = 0;
3192  }
3193  }
3194  else {
3196 
3197  /* raster may have no bands */
3198  if (arg->numbands) {
3199  arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3200  if (arg->band == NULL) {
3201 
3202  rtpg_clip_arg_destroy(arg);
3203  PG_FREE_IF_COPY(pgraster, 0);
3204  lwgeom_free(geom);
3205 
3206  elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3207  PG_RETURN_NULL();
3208  }
3209 
3210  for (i = 0; i < arg->numbands; i++) {
3211  arg->band[i].nband = i;
3212  arg->band[i].hasnodata = 0;
3213  arg->band[i].nodataval = 0;
3214  }
3215  }
3216  }
3217 
3218  /* nodataval (3) */
3219  if (!PG_ARGISNULL(3)) {
3220  array = PG_GETARG_ARRAYTYPE_P(3);
3221  etype = ARR_ELEMTYPE(array);
3222  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3223 
3224  switch (etype) {
3225  case FLOAT4OID:
3226  case FLOAT8OID:
3227  break;
3228  default:
3229  rtpg_clip_arg_destroy(arg);
3230  PG_FREE_IF_COPY(pgraster, 0);
3231  lwgeom_free(geom);
3232  elog(ERROR, "RASTER_clip: Invalid data type for NODATA values");
3233  PG_RETURN_NULL();
3234  break;
3235  }
3236 
3237  deconstruct_array(
3238  array, etype,
3239  typlen, typbyval, typalign,
3240  &e, &nulls, &k
3241  );
3242 
3243  /* it doesn't matter if there are more nodataval */
3244  for (i = 0, j = 0; i < arg->numbands; i++, j++) {
3245  /* cap j to the last nodataval element */
3246  if (j >= k)
3247  j = k - 1;
3248 
3249  if (nulls[j])
3250  continue;
3251 
3252  arg->band[i].hasnodata = 1;
3253  switch (etype) {
3254  case FLOAT4OID:
3255  arg->band[i].nodataval = DatumGetFloat4(e[j]);
3256  break;
3257  case FLOAT8OID:
3258  arg->band[i].nodataval = DatumGetFloat8(e[j]);
3259  break;
3260  }
3261  }
3262  }
3263 
3264  /* get wkb of geometry */
3265  POSTGIS_RT_DEBUG(3, "getting wkb of geometry");
3266  wkb = lwgeom_to_wkb(geom, WKB_SFSQL, &wkb_len);
3267  lwgeom_free(geom);
3268 
3269  /* rasterize geometry */
3271  wkb, wkb_len,
3272  NULL,
3273  0, NULL,
3274  NULL, NULL,
3275  NULL, NULL,
3276  NULL, NULL,
3277  &(gt[1]), &(gt[5]),
3278  NULL, NULL,
3279  &(gt[0]), &(gt[3]),
3280  &(gt[2]), &(gt[4]),
3281  NULL
3282  );
3283 
3284  pfree(wkb);
3285  if (arg->mask == NULL) {
3286  rtpg_clip_arg_destroy(arg);
3287  PG_FREE_IF_COPY(pgraster, 0);
3288  elog(ERROR, "RASTER_clip: Could not rasterize intersection geometry");
3289  PG_RETURN_NULL();
3290  }
3291 
3292  /* set SRID */
3293  rt_raster_set_srid(arg->mask, srid);
3294 
3295  /* run iterator */
3296 
3297  /* init itrset */
3298  itrset = palloc(sizeof(struct rt_iterator_t) * 2);
3299  if (itrset == NULL) {
3300  rtpg_clip_arg_destroy(arg);
3301  PG_FREE_IF_COPY(pgraster, 0);
3302  elog(ERROR, "RASTER_clip: Could not allocate memory for iterator arguments");
3303  PG_RETURN_NULL();
3304  }
3305 
3306  /* one band at a time */
3307  for (i = 0; i < arg->numbands; i++) {
3308  POSTGIS_RT_DEBUGF(4, "band arg %d (nband, hasnodata, nodataval) = (%d, %d, %f)",
3309  i, arg->band[i].nband, arg->band[i].hasnodata, arg->band[i].nodataval);
3310 
3311  band = rt_raster_get_band(arg->raster, arg->band[i].nband);
3312 
3313  /* band metadata */
3314  pixtype = rt_band_get_pixtype(band);
3315 
3316  if (arg->band[i].hasnodata) {
3317  hasnodata = 1;
3318  nodataval = arg->band[i].nodataval;
3319  }
3320  else if (rt_band_get_hasnodata_flag(band)) {
3321  hasnodata = 1;
3322  rt_band_get_nodata(band, &nodataval);
3323  }
3324  else {
3325  hasnodata = 0;
3326  nodataval = rt_band_get_min_value(band);
3327  }
3328 
3329  /* band is NODATA, create NODATA band and continue */
3330  if (rt_band_get_isnodata_flag(band)) {
3331  /* create raster */
3332  if (rtn == NULL) {
3333  noerr = rt_raster_from_two_rasters(arg->raster, arg->mask, arg->extenttype, &rtn, NULL);
3334  if (noerr != ES_NONE) {
3335  rtpg_clip_arg_destroy(arg);
3336  PG_FREE_IF_COPY(pgraster, 0);
3337  elog(ERROR, "RASTER_clip: Could not create output raster");
3338  PG_RETURN_NULL();
3339  }
3340  }
3341 
3342  /* create NODATA band */
3343  if (rt_raster_generate_new_band(rtn, pixtype, nodataval, hasnodata, nodataval, i) < 0) {
3344  rt_raster_destroy(rtn);
3345  rtpg_clip_arg_destroy(arg);
3346  PG_FREE_IF_COPY(pgraster, 0);
3347  elog(ERROR, "RASTER_clip: Could not add NODATA band to output raster");
3348  PG_RETURN_NULL();
3349  }
3350 
3351  continue;
3352  }
3353 
3354  /* raster */
3355  itrset[0].raster = arg->raster;
3356  itrset[0].nband = arg->band[i].nband;
3357  itrset[0].nbnodata = 1;
3358 
3359  /* mask */
3360  itrset[1].raster = arg->mask;
3361  itrset[1].nband = 0;
3362  itrset[1].nbnodata = 1;
3363 
3364  /* pass to iterator */
3365  noerr = rt_raster_iterator(
3366  itrset, 2,
3367  arg->extenttype, NULL,
3368  pixtype,
3369  hasnodata, nodataval,
3370  0, 0,
3371  NULL,
3372  NULL,
3374  &_raster
3375  );
3376 
3377  if (noerr != ES_NONE) {
3378  pfree(itrset);
3379  rtpg_clip_arg_destroy(arg);
3380  if (rtn != NULL) rt_raster_destroy(rtn);
3381  PG_FREE_IF_COPY(pgraster, 0);
3382  elog(ERROR, "RASTER_clip: Could not run raster iterator function");
3383  PG_RETURN_NULL();
3384  }
3385 
3386  /* new raster */
3387  if (rtn == NULL)
3388  rtn = _raster;
3389  /* copy band */
3390  else {
3391  band = rt_raster_get_band(_raster, 0);
3392  if (band == NULL) {
3393  pfree(itrset);
3394  rtpg_clip_arg_destroy(arg);
3395  rt_raster_destroy(_raster);
3396  rt_raster_destroy(rtn);
3397  PG_FREE_IF_COPY(pgraster, 0);
3398  elog(NOTICE, "RASTER_clip: Could not get band from working raster");
3399  PG_RETURN_NULL();
3400  }
3401 
3402  if (rt_raster_add_band(rtn, band, i) < 0) {
3403  pfree(itrset);
3404  rtpg_clip_arg_destroy(arg);
3405  rt_raster_destroy(_raster);
3406  rt_raster_destroy(rtn);
3407  PG_FREE_IF_COPY(pgraster, 0);
3408  elog(ERROR, "RASTER_clip: Could not add new band to output raster");
3409  PG_RETURN_NULL();
3410  }
3411 
3412  rt_raster_destroy(_raster);
3413  }
3414  }
3415 
3416  pfree(itrset);
3417  rtpg_clip_arg_destroy(arg);
3418  PG_FREE_IF_COPY(pgraster, 0);
3419 
3420  pgrtn = rt_raster_serialize(rtn);
3421  rt_raster_destroy(rtn);
3422 
3423  POSTGIS_RT_DEBUG(3, "Finished");
3424 
3425  if (!pgrtn)
3426  PG_RETURN_NULL();
3427 
3428  SET_VARSIZE(pgrtn, pgrtn->size);
3429  PG_RETURN_POINTER(pgrtn);
3430 }
3431 
3436 Datum RASTER_reclass(PG_FUNCTION_ARGS) {
3437  rt_pgraster *pgraster = NULL;
3438  rt_pgraster *pgrtn = NULL;
3439  rt_raster raster = NULL;
3440  rt_band band = NULL;
3441  rt_band newband = NULL;
3442  uint32_t numBands = 0;
3443 
3444  ArrayType *array;
3445  Oid etype;
3446  Datum *e;
3447  bool *nulls;
3448  int16 typlen;
3449  bool typbyval;
3450  char typalign;
3451  int n = 0;
3452 
3453  int i = 0;
3454  int j = 0;
3455  int k = 0;
3456 
3457  uint32_t a = 0;
3458  uint32_t b = 0;
3459  uint32_t c = 0;
3460 
3461  rt_reclassexpr *exprset = NULL;
3462  HeapTupleHeader tup;
3463  bool isnull;
3464  Datum tupv;
3465  uint32_t nband = 0;
3466  char *expr = NULL;
3467  text *exprtext = NULL;
3468  double val = 0;
3469  char *junk = NULL;
3470  int inc_val = 0;
3471  int exc_val = 0;
3472  char *pixeltype = NULL;
3473  text *pixeltypetext = NULL;
3474  rt_pixtype pixtype = PT_END;
3475  double nodataval = 0;
3476  bool hasnodata = FALSE;
3477 
3478  char **comma_set = NULL;
3479  uint32_t comma_n = 0;
3480  char **colon_set = NULL;
3481  uint32_t colon_n = 0;
3482  char **dash_set = NULL;
3483  uint32_t dash_n = 0;
3484 
3485  POSTGIS_RT_DEBUG(3, "RASTER_reclass: Starting");
3486 
3487  /* pgraster is null, return null */
3488  if (PG_ARGISNULL(0))
3489  PG_RETURN_NULL();
3490  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3491 
3492  /* raster */
3493  raster = rt_raster_deserialize(pgraster, FALSE);
3494  if (!raster) {
3495  PG_FREE_IF_COPY(pgraster, 0);
3496  elog(ERROR, "RASTER_reclass: Could not deserialize raster");
3497  PG_RETURN_NULL();
3498  }
3499  numBands = rt_raster_get_num_bands(raster);
3500  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: %d possible bands to be reclassified", numBands);
3501 
3502  /* process set of reclassarg */
3503  POSTGIS_RT_DEBUG(3, "RASTER_reclass: Processing Arg 1 (reclassargset)");
3504  array = PG_GETARG_ARRAYTYPE_P(1);
3505  etype = ARR_ELEMTYPE(array);
3506  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3507 
3508  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
3509  &nulls, &n);
3510 
3511  if (!n) {
3512  elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
3513 
3514  pgrtn = rt_raster_serialize(raster);
3515  rt_raster_destroy(raster);
3516  PG_FREE_IF_COPY(pgraster, 0);
3517  if (!pgrtn)
3518  PG_RETURN_NULL();
3519 
3520  SET_VARSIZE(pgrtn, pgrtn->size);
3521  PG_RETURN_POINTER(pgrtn);
3522  }
3523 
3524  /*
3525  process each element of reclassarg
3526  each element is the index of the band to process, the set
3527  of reclass ranges and the output band's pixeltype
3528  */
3529  for (i = 0; i < n; i++) {
3530  if (nulls[i]) continue;
3531 
3532  /* each element is a tuple */
3533  tup = (HeapTupleHeader) DatumGetPointer(e[i]);
3534  if (NULL == tup) {
3535  elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
3536 
3537  pgrtn = rt_raster_serialize(raster);
3538  rt_raster_destroy(raster);
3539  PG_FREE_IF_COPY(pgraster, 0);
3540  if (!pgrtn)
3541  PG_RETURN_NULL();
3542 
3543  SET_VARSIZE(pgrtn, pgrtn->size);
3544  PG_RETURN_POINTER(pgrtn);
3545  }
3546 
3547  /* band index (1-based) */
3548  tupv = GetAttributeByName(tup, "nband", &isnull);
3549  if (isnull) {
3550  elog(NOTICE, "Invalid argument for reclassargset. Missing value of nband for reclassarg of index %d . Returning original raster", i);
3551 
3552  pgrtn = rt_raster_serialize(raster);
3553  rt_raster_destroy(raster);
3554  PG_FREE_IF_COPY(pgraster, 0);
3555  if (!pgrtn)
3556  PG_RETURN_NULL();
3557 
3558  SET_VARSIZE(pgrtn, pgrtn->size);
3559  PG_RETURN_POINTER(pgrtn);
3560  }
3561  nband = DatumGetInt32(tupv);
3562  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: expression for band %d", nband);
3563 
3564  /* valid band index? */
3565  if (nband < 1 || nband > numBands) {
3566  elog(NOTICE, "Invalid argument for reclassargset. Invalid band index (must use 1-based) for reclassarg of index %d . Returning original raster", i);
3567 
3568  pgrtn = rt_raster_serialize(raster);
3569  rt_raster_destroy(raster);
3570  PG_FREE_IF_COPY(pgraster, 0);
3571  if (!pgrtn)
3572  PG_RETURN_NULL();
3573 
3574  SET_VARSIZE(pgrtn, pgrtn->size);
3575  PG_RETURN_POINTER(pgrtn);
3576  }
3577 
3578  /* reclass expr */
3579  tupv = GetAttributeByName(tup, "reclassexpr", &isnull);
3580  if (isnull) {
3581  elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3582 
3583  pgrtn = rt_raster_serialize(raster);
3584  rt_raster_destroy(raster);
3585  PG_FREE_IF_COPY(pgraster, 0);
3586  if (!pgrtn)
3587  PG_RETURN_NULL();
3588 
3589  SET_VARSIZE(pgrtn, pgrtn->size);
3590  PG_RETURN_POINTER(pgrtn);
3591  }
3592  exprtext = (text *) DatumGetPointer(tupv);
3593  if (NULL == exprtext) {
3594  elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3595 
3596  pgrtn = rt_raster_serialize(raster);
3597  rt_raster_destroy(raster);
3598  PG_FREE_IF_COPY(pgraster, 0);
3599  if (!pgrtn)
3600  PG_RETURN_NULL();
3601 
3602  SET_VARSIZE(pgrtn, pgrtn->size);
3603  PG_RETURN_POINTER(pgrtn);
3604  }
3605  expr = text_to_cstring(exprtext);
3606  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: expr (raw) %s", expr);
3607  expr = rtpg_removespaces(expr);
3608  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: expr (clean) %s", expr);
3609 
3610  /* split string to its components */
3611  /* comma (,) separating rangesets */
3612  comma_set = rtpg_strsplit(expr, ",", &comma_n);
3613  if (comma_n < 1) {
3614  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3615 
3616  pgrtn = rt_raster_serialize(raster);
3617  rt_raster_destroy(raster);
3618  PG_FREE_IF_COPY(pgraster, 0);
3619  if (!pgrtn)
3620  PG_RETURN_NULL();
3621 
3622  SET_VARSIZE(pgrtn, pgrtn->size);
3623  PG_RETURN_POINTER(pgrtn);
3624  }
3625 
3626  /* set of reclass expressions */
3627  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: %d possible expressions", comma_n);
3628  exprset = palloc(comma_n * sizeof(rt_reclassexpr));
3629 
3630  for (a = 0, j = 0; a < comma_n; a++) {
3631  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: map %s", comma_set[a]);
3632 
3633  /* colon (:) separating range "src" and "dst" */
3634  colon_set = rtpg_strsplit(comma_set[a], ":", &colon_n);
3635  if (colon_n != 2) {
3636  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3637  for (k = 0; k < j; k++) pfree(exprset[k]);
3638  pfree(exprset);
3639 
3640  pgrtn = rt_raster_serialize(raster);
3641  rt_raster_destroy(raster);
3642  PG_FREE_IF_COPY(pgraster, 0);
3643  if (!pgrtn)
3644  PG_RETURN_NULL();
3645 
3646  SET_VARSIZE(pgrtn, pgrtn->size);
3647  PG_RETURN_POINTER(pgrtn);
3648  }
3649 
3650  /* allocate mem for reclass expression */
3651  exprset[j] = palloc(sizeof(struct rt_reclassexpr_t));
3652 
3653  for (b = 0; b < colon_n; b++) {
3654  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: range %s", colon_set[b]);
3655 
3656  /* dash (-) separating "min" and "max" */
3657  dash_set = rtpg_strsplit(colon_set[b], "-", &dash_n);
3658  if (dash_n < 1 || dash_n > 3) {
3659  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3660  for (k = 0; k < j; k++) pfree(exprset[k]);
3661  pfree(exprset);
3662 
3663  pgrtn = rt_raster_serialize(raster);
3664  rt_raster_destroy(raster);
3665  PG_FREE_IF_COPY(pgraster, 0);
3666  if (!pgrtn)
3667  PG_RETURN_NULL();
3668 
3669  SET_VARSIZE(pgrtn, pgrtn->size);
3670  PG_RETURN_POINTER(pgrtn);
3671  }
3672 
3673  for (c = 0; c < dash_n; c++) {
3674  /* need to handle: (-9999-100 -> "(", "9999", "100" */
3675  if (
3676  c < 1 &&
3677  strlen(dash_set[c]) == 1 && (
3678  strchr(dash_set[c], '(') != NULL ||
3679  strchr(dash_set[c], '[') != NULL ||
3680  strchr(dash_set[c], ')') != NULL ||
3681  strchr(dash_set[c], ']') != NULL
3682  )
3683  ) {
3684  uint32_t dash_it;
3685  junk = palloc(sizeof(char) * (strlen(dash_set[c + 1]) + 2));
3686  if (NULL == junk) {
3687  for (k = 0; k <= j; k++) pfree(exprset[k]);
3688  pfree(exprset);
3689  rt_raster_destroy(raster);
3690  PG_FREE_IF_COPY(pgraster, 0);
3691 
3692  elog(ERROR, "RASTER_reclass: Could not allocate memory");
3693  PG_RETURN_NULL();
3694  }
3695 
3696  sprintf(junk, "%s%s", dash_set[c], dash_set[c + 1]);
3697  c++;
3698  dash_set[c] = repalloc(dash_set[c], sizeof(char) * (strlen(junk) + 1));
3699  strcpy(dash_set[c], junk);
3700  pfree(junk);
3701 
3702  /* rebuild dash_set */
3703  for (dash_it = 1; dash_it < dash_n; dash_it++) {
3704  dash_set[dash_it - 1] = repalloc(dash_set[dash_it - 1], (strlen(dash_set[dash_it]) + 1) * sizeof(char));
3705  strcpy(dash_set[dash_it - 1], dash_set[dash_it]);
3706  }
3707  dash_n--;
3708  c--;
3709  pfree(dash_set[dash_n]);
3710  dash_set = repalloc(dash_set, sizeof(char *) * dash_n);
3711  }
3712 
3713  /* there shouldn't be more than two in dash_n */
3714  if (c < 1 && dash_n > 2) {
3715  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3716  for (k = 0; k < j; k++) pfree(exprset[k]);
3717  pfree(exprset);
3718 
3719  pgrtn = rt_raster_serialize(raster);
3720  rt_raster_destroy(raster);
3721  PG_FREE_IF_COPY(pgraster, 0);
3722  if (!pgrtn)
3723  PG_RETURN_NULL();
3724 
3725  SET_VARSIZE(pgrtn, pgrtn->size);
3726  PG_RETURN_POINTER(pgrtn);
3727  }
3728 
3729  /* check interval flags */
3730  exc_val = 0;
3731  inc_val = 1;
3732  /* range */
3733  if (dash_n != 1) {
3734  /* min */
3735  if (c < 1) {
3736  if (
3737  strchr(dash_set[c], ')') != NULL ||
3738  strchr(dash_set[c], ']') != NULL
3739  ) {
3740  exc_val = 1;
3741  inc_val = 1;
3742  }
3743  else if (strchr(dash_set[c], '(') != NULL){
3744  inc_val = 0;
3745  }
3746  else {
3747  inc_val = 1;
3748  }
3749  }
3750  /* max */
3751  else {
3752  if (
3753  strrchr(dash_set[c], '(') != NULL ||
3754  strrchr(dash_set[c], '[') != NULL
3755  ) {
3756  exc_val = 1;
3757  inc_val = 0;
3758  }
3759  else if (strrchr(dash_set[c], ']') != NULL) {
3760  inc_val = 1;
3761  }
3762  else {
3763  inc_val = 0;
3764  }
3765  }
3766  }
3767  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: exc_val %d inc_val %d", exc_val, inc_val);
3768 
3769  /* remove interval flags */
3770  dash_set[c] = rtpg_chartrim(dash_set[c], "()[]");
3771  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (char) %s", dash_set[c]);
3772 
3773  /* value from string to double */
3774  errno = 0;
3775  val = strtod(dash_set[c], &junk);
3776  if (errno != 0 || dash_set[c] == junk) {
3777  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3778  for (k = 0; k < j; k++) pfree(exprset[k]);
3779  pfree(exprset);
3780 
3781  pgrtn = rt_raster_serialize(raster);
3782  rt_raster_destroy(raster);
3783  PG_FREE_IF_COPY(pgraster, 0);
3784  if (!pgrtn)
3785  PG_RETURN_NULL();
3786 
3787  SET_VARSIZE(pgrtn, pgrtn->size);
3788  PG_RETURN_POINTER(pgrtn);
3789  }
3790  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (double) %f", val);
3791 
3792  /* strsplit removes dash (a.k.a. negative sign), compare now to restore */
3793  if (c < 1)
3794  junk = strstr(colon_set[b], dash_set[c]);
3795  else
3796  junk = rtpg_strrstr(colon_set[b], dash_set[c]);
3798  4,
3799  "(colon_set[%d], dash_set[%d], junk) = (%s, %s, %s)",
3800  b, c, colon_set[b], dash_set[c], junk
3801  );
3802  /* not beginning of string */
3803  if (junk != colon_set[b]) {
3804  /* prior is a dash */
3805  if (*(junk - 1) == '-') {
3806  /* prior is beginning of string or prior - 1 char is dash, negative number */
3807  if (
3808  ((junk - 1) == colon_set[b]) ||
3809  (*(junk - 2) == '-') ||
3810  (*(junk - 2) == '[') ||
3811  (*(junk - 2) == '(')
3812  ) {
3813  val *= -1.;
3814  }
3815  }
3816  }
3817  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (double) %f", val);
3818 
3819  /* src */
3820  if (b < 1) {
3821  /* singular value */
3822  if (dash_n == 1) {
3823  exprset[j]->src.exc_min = exprset[j]->src.exc_max = exc_val;
3824  exprset[j]->src.inc_min = exprset[j]->src.inc_max = inc_val;
3825  exprset[j]->src.min = exprset[j]->src.max = val;
3826  }
3827  /* min */
3828  else if (c < 1) {
3829  exprset[j]->src.exc_min = exc_val;
3830  exprset[j]->src.inc_min = inc_val;
3831  exprset[j]->src.min = val;
3832  }
3833  /* max */
3834  else {
3835  exprset[j]->src.exc_max = exc_val;
3836  exprset[j]->src.inc_max = inc_val;
3837  exprset[j]->src.max = val;
3838  }
3839  }
3840  /* dst */
3841  else {
3842  /* singular value */
3843  if (dash_n == 1)
3844  exprset[j]->dst.min = exprset[j]->dst.max = val;
3845  /* min */
3846  else if (c < 1)
3847  exprset[j]->dst.min = val;
3848  /* max */
3849  else
3850  exprset[j]->dst.max = val;
3851  }
3852  }
3853  pfree(dash_set);
3854  }
3855  pfree(colon_set);
3856 
3857  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: or: %f - %f nr: %f - %f"
3858  , exprset[j]->src.min
3859  , exprset[j]->src.max
3860  , exprset[j]->dst.min
3861  , exprset[j]->dst.max
3862  );
3863  j++;
3864  }
3865  pfree(comma_set);
3866 
3867  /* pixel type */
3868  tupv = GetAttributeByName(tup, "pixeltype", &isnull);
3869  if (isnull) {
3870  elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3871 
3872  pgrtn = rt_raster_serialize(raster);
3873  rt_raster_destroy(raster);
3874  PG_FREE_IF_COPY(pgraster, 0);
3875  if (!pgrtn)
3876  PG_RETURN_NULL();
3877 
3878  SET_VARSIZE(pgrtn, pgrtn->size);
3879  PG_RETURN_POINTER(pgrtn);
3880  }
3881  pixeltypetext = (text *) DatumGetPointer(tupv);
3882  if (NULL == pixeltypetext) {
3883  elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3884 
3885  pgrtn = rt_raster_serialize(raster);
3886  rt_raster_destroy(raster);
3887  PG_FREE_IF_COPY(pgraster, 0);
3888  if (!pgrtn)
3889  PG_RETURN_NULL();
3890 
3891  SET_VARSIZE(pgrtn, pgrtn->size);
3892  PG_RETURN_POINTER(pgrtn);
3893  }
3894  pixeltype = text_to_cstring(pixeltypetext);
3895  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: pixeltype %s", pixeltype);
3896  pixtype = rt_pixtype_index_from_name(pixeltype);
3897 
3898  /* nodata */
3899  tupv = GetAttributeByName(tup, "nodataval", &isnull);
3900  if (isnull) {
3901  nodataval = 0;
3902  hasnodata = FALSE;
3903  }
3904  else {
3905  nodataval = DatumGetFloat8(tupv);
3906  hasnodata = TRUE;
3907  }
3908  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: nodataval %f", nodataval);
3909  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: hasnodata %d", hasnodata);
3910 
3911  /* do reclass */
3912  band = rt_raster_get_band(raster, nband - 1);
3913  if (!band) {
3914  elog(NOTICE, "Could not find raster band of index %d. Returning original raster", nband);
3915  for (k = 0; k < j; k++) pfree(exprset[k]);
3916  pfree(exprset);
3917 
3918  pgrtn = rt_raster_serialize(raster);
3919  rt_raster_destroy(raster);
3920  PG_FREE_IF_COPY(pgraster, 0);
3921  if (!pgrtn)
3922  PG_RETURN_NULL();
3923 
3924  SET_VARSIZE(pgrtn, pgrtn->size);
3925  PG_RETURN_POINTER(pgrtn);
3926  }
3927  newband = rt_band_reclass(band, pixtype, hasnodata, nodataval, exprset, j);
3928  if (!newband) {
3929  for (k = 0; k < j; k++) pfree(exprset[k]);
3930  pfree(exprset);
3931 
3932  rt_raster_destroy(raster);
3933  PG_FREE_IF_COPY(pgraster, 0);
3934 
3935  elog(ERROR, "RASTER_reclass: Could not reclassify raster band of index %d", nband);
3936  PG_RETURN_NULL();
3937  }
3938 
3939  /* replace old band with new band */
3940  if (rt_raster_replace_band(raster, newband, nband - 1) == NULL) {
3941  for (k = 0; k < j; k++) pfree(exprset[k]);
3942  pfree(exprset);
3943 
3944  rt_band_destroy(newband);
3945  rt_raster_destroy(raster);
3946  PG_FREE_IF_COPY(pgraster, 0);
3947 
3948  elog(ERROR, "RASTER_reclass: Could not replace raster band of index %d with reclassified band", nband);
3949  PG_RETURN_NULL();
3950  }
3951 
3952  /* old band is in the variable band */
3953  rt_band_destroy(band);
3954 
3955  /* free exprset */
3956  for (k = 0; k < j; k++) pfree(exprset[k]);
3957  pfree(exprset);
3958  }
3959 
3960  pgrtn = rt_raster_serialize(raster);
3961  rt_raster_destroy(raster);
3962  PG_FREE_IF_COPY(pgraster, 0);
3963  if (!pgrtn)
3964  PG_RETURN_NULL();
3965 
3966  POSTGIS_RT_DEBUG(3, "RASTER_reclass: Finished");
3967 
3968  SET_VARSIZE(pgrtn, pgrtn->size);
3969  PG_RETURN_POINTER(pgrtn);
3970 }
3971 
3972 /* ---------------------------------------------------------------- */
3973 /* apply colormap to specified band of a raster */
3974 /* ---------------------------------------------------------------- */
3975 
3979  int nband; /* 1-based */
3982 
3985 
3986  char **entry;
3988  char **element;
3990 };
3991 
3992 static rtpg_colormap_arg
3994  rtpg_colormap_arg arg = NULL;
3995 
3996  arg = palloc(sizeof(struct rtpg_colormap_arg_t));
3997  if (arg == NULL) {
3998  elog(ERROR, "rtpg_colormap_arg: Could not allocate memory for function arguments");
3999  return NULL;
4000  }
4001 
4002  arg->raster = NULL;
4003  arg->nband = 1;
4004  arg->band = NULL;
4005  arg->bandstats = NULL;
4006 
4007  arg->colormap = palloc(sizeof(struct rt_colormap_t));
4008  if (arg->colormap == NULL) {
4009  elog(ERROR, "rtpg_colormap_arg: Could not allocate memory for function arguments");
4010  return NULL;
4011  }
4012  arg->colormap->nentry = 0;
4013  arg->colormap->entry = NULL;
4014  arg->colormap->ncolor = 4; /* assume RGBA */
4015  arg->colormap->method = CM_INTERPOLATE;
4016  arg->nodataentry = -1;
4017 
4018  arg->entry = NULL;
4019  arg->nentry = 0;
4020  arg->element = NULL;
4021  arg->nelement = 0;
4022 
4023  return arg;
4024 }
4025 
4026 static void
4027 rtpg_colormap_arg_destroy(rtpg_colormap_arg arg) {
4028  uint32_t i = 0;
4029  if (arg->raster != NULL)
4030  rt_raster_destroy(arg->raster);
4031 
4032  if (arg->bandstats != NULL)
4033  pfree(arg->bandstats);
4034 
4035  if (arg->colormap != NULL) {
4036  if (arg->colormap->entry != NULL)
4037  pfree(arg->colormap->entry);
4038  pfree(arg->colormap);
4039  }
4040 
4041  if (arg->nentry) {
4042  for (i = 0; i < arg->nentry; i++) {
4043  if (arg->entry[i] != NULL)
4044  pfree(arg->entry[i]);
4045  }
4046  pfree(arg->entry);
4047  }
4048 
4049  if (arg->nelement) {
4050  for (i = 0; i < arg->nelement; i++)
4051  pfree(arg->element[i]);
4052  pfree(arg->element);
4053  }
4054 
4055  pfree(arg);
4056  arg = NULL;
4057 }
4058 
4060 Datum RASTER_colorMap(PG_FUNCTION_ARGS)
4061 {
4062  rt_pgraster *pgraster = NULL;
4063  rtpg_colormap_arg arg = NULL;
4064  char *junk = NULL;
4065  rt_raster raster = NULL;
4066 
4067  POSTGIS_RT_DEBUG(3, "RASTER_colorMap: Starting");
4068 
4069  /* pgraster is NULL, return NULL */
4070  if (PG_ARGISNULL(0))
4071  PG_RETURN_NULL();
4072 
4073  /* init arg */
4074  arg = rtpg_colormap_arg_init();
4075  if (arg == NULL) {
4076  elog(ERROR, "RASTER_colorMap: Could not initialize argument structure");
4077  PG_RETURN_NULL();
4078  }
4079 
4080  /* raster (0) */
4081  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4082 
4083  /* raster object */
4084  arg->raster = rt_raster_deserialize(pgraster, FALSE);
4085  if (!arg->raster) {
4087  PG_FREE_IF_COPY(pgraster, 0);
4088  elog(ERROR, "RASTER_colorMap: Could not deserialize raster");
4089  PG_RETURN_NULL();
4090  }
4091 
4092  /* nband (1) */
4093  if (!PG_ARGISNULL(1))
4094  arg->nband = PG_GETARG_INT32(1);
4095  POSTGIS_RT_DEBUGF(4, "nband = %d", arg->nband);
4096 
4097  /* check that band works */
4098  if (!rt_raster_has_band(arg->raster, arg->nband - 1)) {
4099  elog(NOTICE, "Raster does not have band at index %d. Returning empty raster", arg->nband);
4100 
4101  raster = rt_raster_clone(arg->raster, 0);
4102  if (raster == NULL) {
4104  PG_FREE_IF_COPY(pgraster, 0);
4105  elog(ERROR, "RASTER_colorMap: Could not create empty raster");
4106  PG_RETURN_NULL();
4107  }
4108 
4110  PG_FREE_IF_COPY(pgraster, 0);
4111 
4112  pgraster = rt_raster_serialize(raster);
4113  rt_raster_destroy(raster);
4114  if (pgraster == NULL)
4115  PG_RETURN_NULL();
4116 
4117  SET_VARSIZE(pgraster, ((rt_pgraster*) pgraster)->size);
4118  PG_RETURN_POINTER(pgraster);
4119  }
4120 
4121  /* get band */
4122  arg->band = rt_raster_get_band(arg->raster, arg->nband - 1);
4123  if (arg->band == NULL) {
4124  int nband = arg->nband;
4126  PG_FREE_IF_COPY(pgraster, 0);
4127  elog(ERROR, "RASTER_colorMap: Could not get band at index %d", nband);
4128  PG_RETURN_NULL();
4129  }
4130 
4131  /* method (3) */
4132  if (!PG_ARGISNULL(3)) {
4133  char *method = NULL;
4134  char *tmp = text_to_cstring(PG_GETARG_TEXT_P(3));
4135  POSTGIS_RT_DEBUGF(4, "raw method = %s", tmp);
4136 
4137  method = rtpg_trim(tmp);
4138  pfree(tmp);
4139  method = rtpg_strtoupper(method);
4140 
4141  if (strcmp(method, "INTERPOLATE") == 0)
4142  arg->colormap->method = CM_INTERPOLATE;
4143  else if (strcmp(method, "EXACT") == 0)
4144  arg->colormap->method = CM_EXACT;
4145  else if (strcmp(method, "NEAREST") == 0)
4146  arg->colormap->method = CM_NEAREST;
4147  else {
4148  elog(NOTICE, "Unknown value provided for method. Defaulting to INTERPOLATE");
4149  arg->colormap->method = CM_INTERPOLATE;
4150  }
4151  }
4152  /* default to INTERPOLATE */
4153  else
4154  arg->colormap->method = CM_INTERPOLATE;
4155  POSTGIS_RT_DEBUGF(4, "method = %d", arg->colormap->method);
4156 
4157  /* colormap (2) */
4158  if (PG_ARGISNULL(2)) {
4160  PG_FREE_IF_COPY(pgraster, 0);
4161  elog(ERROR, "RASTER_colorMap: Value must be provided for colormap");
4162  PG_RETURN_NULL();
4163  }
4164  else {
4165  char *tmp = NULL;
4166  char *colormap = text_to_cstring(PG_GETARG_TEXT_P(2));
4167  char *_entry;
4168  char *_element;
4169  uint32_t i = 0;
4170  uint32_t j = 0;
4171 
4172  POSTGIS_RT_DEBUGF(4, "colormap = %s", colormap);
4173 
4174  /* empty string */
4175  if (!strlen(colormap)) {
4177  PG_FREE_IF_COPY(pgraster, 0);
4178  elog(ERROR, "RASTER_colorMap: Value must be provided for colormap");
4179  PG_RETURN_NULL();
4180  }
4181 
4182  arg->entry = rtpg_strsplit(colormap, "\n", &(arg->nentry));
4183  pfree(colormap);
4184  if (arg->nentry < 1) {
4186  PG_FREE_IF_COPY(pgraster, 0);
4187  elog(ERROR, "RASTER_colorMap: Could not process the value provided for colormap");
4188  PG_RETURN_NULL();
4189  }
4190 
4191  /* allocate the max # of colormap entries */
4192  arg->colormap->entry = palloc(sizeof(struct rt_colormap_entry_t) * arg->nentry);
4193  if (arg->colormap->entry == NULL) {
4195  PG_FREE_IF_COPY(pgraster, 0);
4196  elog(ERROR, "RASTER_colorMap: Could not allocate memory for colormap entries");
4197  PG_RETURN_NULL();
4198  }
4199  memset(arg->colormap->entry, 0, sizeof(struct rt_colormap_entry_t) * arg->nentry);
4200 
4201  /* each entry */
4202  for (i = 0; i < arg->nentry; i++) {
4203  /* substitute space for other delimiters */
4204  tmp = rtpg_strreplace(arg->entry[i], ":", " ", NULL);
4205  _entry = rtpg_strreplace(tmp, ",", " ", NULL);
4206  pfree(tmp);
4207  tmp = rtpg_strreplace(_entry, "\t", " ", NULL);
4208  pfree(_entry);
4209  _entry = rtpg_trim(tmp);
4210  pfree(tmp);
4211 
4212  POSTGIS_RT_DEBUGF(4, "Processing entry[%d] = %s", i, arg->entry[i]);
4213  POSTGIS_RT_DEBUGF(4, "Cleaned entry[%d] = %s", i, _entry);
4214 
4215  /* empty entry, continue */
4216  if (!strlen(_entry)) {
4217  POSTGIS_RT_DEBUGF(3, "Skipping empty entry[%d]", i);
4218  pfree(_entry);
4219  continue;
4220  }
4221 
4222  arg->element = rtpg_strsplit(_entry, " ", &(arg->nelement));
4223  pfree(_entry);
4224  if (arg->nelement < 2) {
4226  PG_FREE_IF_COPY(pgraster, 0);
4227  elog(ERROR, "RASTER_colorMap: Could not process colormap entry %d", i + 1);
4228  PG_RETURN_NULL();
4229  }
4230  else if (arg->nelement > 5) {
4231  elog(NOTICE, "More than five elements in colormap entry %d. Using at most five elements", i + 1);
4232  arg->nelement = 5;
4233  }
4234 
4235  /* smallest # of colors */
4236  if (((int)arg->nelement - 1) < arg->colormap->ncolor)
4237  arg->colormap->ncolor = arg->nelement - 1;
4238 
4239  /* each element of entry */
4240  for (j = 0; j < arg->nelement; j++) {
4241 
4242  _element = rtpg_trim(arg->element[j]);
4243  _element = rtpg_strtoupper(_element);
4244  POSTGIS_RT_DEBUGF(4, "Processing entry[%d][%d] = %s", i, j, arg->element[j]);
4245  POSTGIS_RT_DEBUGF(4, "Cleaned entry[%d][%d] = %s", i, j, _element);
4246 
4247  /* first element is ALWAYS a band value, percentage OR "nv" string */
4248  if (j == 0) {
4249  char *percent = NULL;
4250 
4251  /* NODATA */
4252  if (
4253  strcmp(_element, "NV") == 0 ||
4254  strcmp(_element, "NULL") == 0 ||
4255  strcmp(_element, "NODATA") == 0
4256  ) {
4257  POSTGIS_RT_DEBUG(4, "Processing NODATA string");
4258 
4259  if (arg->nodataentry > -1) {
4260  elog(NOTICE, "More than one NODATA entry found. Using only the first one");
4261  }
4262  else {
4263  arg->colormap->entry[arg->colormap->nentry].isnodata = 1;
4264  /* no need to set value as value comes from band's NODATA */
4265  arg->colormap->entry[arg->colormap->nentry].value = 0;
4266  }
4267  }
4268  /* percent value */
4269  else if ((percent = strchr(_element, '%')) != NULL) {
4270  double value;
4271  POSTGIS_RT_DEBUG(4, "Processing percent string");
4272 
4273  /* get the band stats */
4274  if (arg->bandstats == NULL) {
4275  POSTGIS_RT_DEBUG(4, "Getting band stats");
4276 
4277  arg->bandstats = rt_band_get_summary_stats(arg->band, 1, 1, 0, NULL, NULL, NULL);
4278  if (arg->bandstats == NULL) {
4279  pfree(_element);
4281  PG_FREE_IF_COPY(pgraster, 0);
4282  elog(ERROR, "RASTER_colorMap: Could not get band's summary stats to process percentages");
4283  PG_RETURN_NULL();
4284  }
4285  }
4286 
4287  /* get the string before the percent char */
4288  tmp = palloc(sizeof(char) * (percent - _element + 1));
4289  if (tmp == NULL) {
4290  pfree(_element);
4292  PG_FREE_IF_COPY(pgraster, 0);
4293  elog(ERROR, "RASTER_colorMap: Could not allocate memory for value of percentage");
4294  PG_RETURN_NULL();
4295  }
4296 
4297  memcpy(tmp, _element, percent - _element);
4298  tmp[percent - _element] = '\0';
4299  POSTGIS_RT_DEBUGF(4, "Percent value = %s", tmp);
4300 
4301  /* get percentage value */
4302  errno = 0;
4303  value = strtod(tmp, NULL);
4304  pfree(tmp);
4305  if (errno != 0 || _element == junk) {
4306  pfree(_element);
4308  PG_FREE_IF_COPY(pgraster, 0);
4309  elog(ERROR, "RASTER_colorMap: Could not process percent string to value");
4310  PG_RETURN_NULL();
4311  }
4312 
4313  /* check percentage */
4314  if (value < 0.) {
4315  elog(NOTICE, "Percentage values cannot be less than zero. Defaulting to zero");
4316  value = 0.;
4317  }
4318  else if (value > 100.) {
4319  elog(NOTICE, "Percentage values cannot be greater than 100. Defaulting to 100");
4320  value = 100.;
4321  }
4322 
4323  /* get the true pixel value */
4324  /* TODO: should the percentage be quantile based? */
4325  arg->colormap->entry[arg->colormap->nentry].value = ((value / 100.) * (arg->bandstats->max - arg->bandstats->min)) + arg->bandstats->min;
4326  }
4327  /* straight value */
4328  else {
4329  errno = 0;
4330  arg->colormap->entry[arg->colormap->nentry].value = strtod(_element, &junk);
4331  if (errno != 0 || _element == junk) {
4332  pfree(_element);
4334  PG_FREE_IF_COPY(pgraster, 0);
4335  elog(ERROR, "RASTER_colorMap: Could not process string to value");
4336  PG_RETURN_NULL();
4337  }
4338  }
4339 
4340  }
4341  /* RGB values (0 - 255) */
4342  else {
4343  int value = 0;
4344 
4345  errno = 0;
4346  value = (int) strtod(_element, &junk);
4347  if (errno != 0 || _element == junk) {
4348  pfree(_element);
4350  PG_FREE_IF_COPY(pgraster, 0);
4351  elog(ERROR, "RASTER_colorMap: Could not process string to value");
4352  PG_RETURN_NULL();
4353  }
4354 
4355  if (value > 255) {
4356  elog(NOTICE, "RGBA value cannot be greater than 255. Defaulting to 255");
4357  value = 255;
4358  }
4359  else if (value < 0) {
4360  elog(NOTICE, "RGBA value cannot be less than zero. Defaulting to zero");
4361  value = 0;
4362  }
4363  arg->colormap->entry[arg->colormap->nentry].color[j - 1] = value;
4364  }
4365 
4366  pfree(_element);
4367  }
4368 
4369  POSTGIS_RT_DEBUGF(4, "colormap->entry[%d] (isnodata, value, R, G, B, A) = (%d, %f, %d, %d, %d, %d)",
4370  arg->colormap->nentry,
4371  arg->colormap->entry[arg->colormap->nentry].isnodata,
4372  arg->colormap->entry[arg->colormap->nentry].value,
4373  arg->colormap->entry[arg->colormap->nentry].color[0],
4374  arg->colormap->entry[arg->colormap->nentry].color[1],
4375  arg->colormap->entry[arg->colormap->nentry].color[2],
4376  arg->colormap->entry[arg->colormap->nentry].color[3]
4377  );
4378 
4379  arg->colormap->nentry++;
4380  }
4381 
4382  POSTGIS_RT_DEBUGF(4, "colormap->nentry = %d", arg->colormap->nentry);
4383  POSTGIS_RT_DEBUGF(4, "colormap->ncolor = %d", arg->colormap->ncolor);
4384  }
4385 
4386  /* call colormap */
4387  raster = rt_raster_colormap(arg->raster, arg->nband - 1, arg->colormap);
4388  if (raster == NULL) {
4390  PG_FREE_IF_COPY(pgraster, 0);
4391  elog(ERROR, "RASTER_colorMap: Could not create new raster with applied colormap");
4392  PG_RETURN_NULL();
4393  }
4394 
4396  PG_FREE_IF_COPY(pgraster, 0);
4397  pgraster = rt_raster_serialize(raster);
4398  rt_raster_destroy(raster);
4399 
4400  POSTGIS_RT_DEBUG(3, "RASTER_colorMap: Done");
4401 
4402  if (pgraster == NULL)
4403  PG_RETURN_NULL();
4404 
4405  SET_VARSIZE(pgraster, ((rt_pgraster*) pgraster)->size);
4406  PG_RETURN_POINTER(pgraster);
4407 }
4408 
4410 Datum RASTER_mapAlgebraExpr(PG_FUNCTION_ARGS)
4411 {
4412  rt_pgraster *pgraster = NULL;
4413  rt_pgraster *pgrtn = NULL;
4414  rt_raster raster = NULL;
4415  rt_raster newrast = NULL;
4416  rt_band band = NULL;
4417  rt_band newband = NULL;
4418  int x, y, nband, width, height;
4419  double r;
4420  double newnodatavalue = 0.0;
4421  double newinitialvalue = 0.0;
4422  double newval = 0.0;
4423  char *newexpr = NULL;
4424  char *initexpr = NULL;
4425  char *expression = NULL;
4426  int hasnodataval = 0;
4427  double nodataval = 0.;
4428  rt_pixtype newpixeltype;
4429  int skipcomputation = 0;
4430  int len = 0;
4431  const int argkwcount = 3;
4432  enum KEYWORDS { kVAL=0, kX=1, kY=2 };
4433  char *argkw[] = {"[rast]", "[rast.x]", "[rast.y]"};
4434  Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
4435  int argcount = 0;
4436  Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
4437  uint8_t argpos[3] = {0};
4438  char place[12];
4439  int idx = 0;
4440  int ret = -1;
4441  TupleDesc tupdesc;
4442  SPIPlanPtr spi_plan = NULL;
4443  SPITupleTable * tuptable = NULL;
4444  HeapTuple tuple;
4445  char * strFromText = NULL;
4446  Datum *values = NULL;
4447  Datum datum = (Datum)NULL;
4448  char *nulls = NULL;
4449  bool isnull = FALSE;
4450  int i = 0;
4451  int j = 0;
4452 
4453  POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraExpr: Starting...");
4454 
4455  /* Check raster */
4456  if (PG_ARGISNULL(0)) {
4457  elog(NOTICE, "Raster is NULL. Returning NULL");
4458  PG_RETURN_NULL();
4459  }
4460 
4461 
4462  /* Deserialize raster */
4463  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4464  raster = rt_raster_deserialize(pgraster, FALSE);
4465  if (NULL == raster) {
4466  PG_FREE_IF_COPY(pgraster, 0);
4467  elog(ERROR, "RASTER_mapAlgebraExpr: Could not deserialize raster");
4468  PG_RETURN_NULL();
4469  }
4470 
4471  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting arguments...");
4472 
4473  if (PG_ARGISNULL(1))
4474  nband = 1;
4475  else
4476  nband = PG_GETARG_INT32(1);
4477 
4478  if (nband < 1)
4479  nband = 1;
4480 
4481 
4482  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new empty raster...");
4483 
4488  width = rt_raster_get_width(raster);
4489  height = rt_raster_get_height(raster);
4490 
4491  newrast = rt_raster_new(width, height);
4492 
4493  if ( NULL == newrast ) {
4494  PG_FREE_IF_COPY(pgraster, 0);
4495  elog(ERROR, "RASTER_mapAlgebraExpr: Could not create a new raster");
4496  PG_RETURN_NULL();
4497  }
4498 
4499  rt_raster_set_scale(newrast,
4500  rt_raster_get_x_scale(raster),
4501  rt_raster_get_y_scale(raster));
4502 
4503  rt_raster_set_offsets(newrast,
4504  rt_raster_get_x_offset(raster),
4505  rt_raster_get_y_offset(raster));
4506 
4507  rt_raster_set_skews(newrast,
4508  rt_raster_get_x_skew(raster),
4509  rt_raster_get_y_skew(raster));
4510 
4511  rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
4512 
4513 
4518  if (rt_raster_is_empty(newrast))
4519  {
4520  elog(NOTICE, "Raster is empty. Returning an empty raster");
4521  rt_raster_destroy(raster);
4522  PG_FREE_IF_COPY(pgraster, 0);
4523 
4524  pgrtn = rt_raster_serialize(newrast);
4525  rt_raster_destroy(newrast);
4526  if (NULL == pgrtn) {
4527 
4528  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4529  PG_RETURN_NULL();
4530  }
4531 
4532  SET_VARSIZE(pgrtn, pgrtn->size);
4533  PG_RETURN_POINTER(pgrtn);
4534  }
4535 
4536 
4537  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Getting raster band %d...", nband);
4538 
4543  if (!rt_raster_has_band(raster, nband - 1)) {
4544  elog(NOTICE, "Raster does not have the required band. Returning a raster "
4545  "without a band");
4546  rt_raster_destroy(raster);
4547  PG_FREE_IF_COPY(pgraster, 0);
4548 
4549  pgrtn = rt_raster_serialize(newrast);
4550  rt_raster_destroy(newrast);
4551  if (NULL == pgrtn) {
4552  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4553  PG_RETURN_NULL();
4554  }
4555 
4556  SET_VARSIZE(pgrtn, pgrtn->size);
4557  PG_RETURN_POINTER(pgrtn);
4558  }
4559 
4560  /* Get the raster band */
4561  band = rt_raster_get_band(raster, nband - 1);
4562  if ( NULL == band ) {
4563  elog(NOTICE, "Could not get the required band. Returning a raster "
4564  "without a band");
4565  rt_raster_destroy(raster);
4566  PG_FREE_IF_COPY(pgraster, 0);
4567 
4568  pgrtn = rt_raster_serialize(newrast);
4569  rt_raster_destroy(newrast);
4570  if (NULL == pgrtn) {
4571  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4572  PG_RETURN_NULL();
4573  }
4574 
4575  SET_VARSIZE(pgrtn, pgrtn->size);
4576  PG_RETURN_POINTER(pgrtn);
4577  }
4578 
4579  /*
4580  * Get NODATA value
4581  */
4582  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting NODATA value for band...");
4583 
4584  if (rt_band_get_hasnodata_flag(band)) {
4585  rt_band_get_nodata(band, &newnodatavalue);
4586  }
4587 
4588  else {
4589  newnodatavalue = rt_band_get_min_value(band);
4590  }
4591 
4592  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: NODATA value for band: = %f",
4593  newnodatavalue);
4594 
4600  newinitialvalue = newnodatavalue;
4601 
4605  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Setting pixeltype...");
4606 
4607  if (PG_ARGISNULL(2)) {
4608  newpixeltype = rt_band_get_pixtype(band);
4609  }
4610 
4611  else {
4612  strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
4613  newpixeltype = rt_pixtype_index_from_name(strFromText);
4614  pfree(strFromText);
4615  if (newpixeltype == PT_END)
4616  newpixeltype = rt_band_get_pixtype(band);
4617  }
4618 
4619  if (newpixeltype == PT_END) {
4620  PG_FREE_IF_COPY(pgraster, 0);
4621  elog(ERROR, "RASTER_mapAlgebraExpr: Invalid pixeltype");
4622  PG_RETURN_NULL();
4623  }
4624 
4625  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Pixeltype set to %s",
4626  rt_pixtype_name(newpixeltype));
4627 
4628 
4629  /* Construct expression for raster values */
4630  if (!PG_ARGISNULL(3)) {
4631  expression = text_to_cstring(PG_GETARG_TEXT_P(3));
4632  len = strlen("SELECT (") + strlen(expression) + strlen(")::double precision");
4633  initexpr = (char *)palloc(len + 1);
4634 
4635  strncpy(initexpr, "SELECT (", strlen("SELECT ("));
4636  strncpy(initexpr + strlen("SELECT ("), expression, strlen(expression));
4637  strncpy(initexpr + strlen("SELECT (") + strlen(expression), ")::double precision", strlen(")::double precision"));
4638  initexpr[len] = '\0';
4639 
4640  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression is %s", initexpr);
4641 
4642  /* We don't need this memory */
4643  /*
4644  pfree(expression);
4645  expression = NULL;
4646  */
4647  }
4648 
4649 
4650 
4656  if (!PG_ARGISNULL(4)) {
4657  hasnodataval = 1;
4658  nodataval = PG_GETARG_FLOAT8(4);
4659  newinitialvalue = nodataval;
4660 
4661  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new initial value = %f",
4662  newinitialvalue);
4663  }
4664  else
4665  hasnodataval = 0;
4666 
4667 
4668 
4674  if (rt_band_get_isnodata_flag(band)) {
4675 
4676  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Band is a nodata band, returning "
4677  "a raster filled with nodata");
4678 
4679  ret = rt_raster_generate_new_band(newrast, newpixeltype,
4680  newinitialvalue, TRUE, newnodatavalue, 0);
4681 
4682  /* Free memory */
4683  if (initexpr)
4684  pfree(initexpr);
4685  rt_raster_destroy(raster);
4686  PG_FREE_IF_COPY(pgraster, 0);
4687 
4688  /* Serialize created raster */
4689  pgrtn = rt_raster_serialize(newrast);
4690  rt_raster_destroy(newrast);
4691  if (NULL == pgrtn) {
4692  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4693  PG_RETURN_NULL();
4694  }
4695 
4696  SET_VARSIZE(pgrtn, pgrtn->size);
4697  PG_RETURN_POINTER(pgrtn);
4698  }
4699 
4700 
4705  if (initexpr != NULL && ( !strcmp(initexpr, "SELECT [rast]") || !strcmp(initexpr, "SELECT [rast.val]") ) && !hasnodataval) {
4706 
4707  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression resumes to RAST. "
4708  "Returning raster with band %d from original raster", nband);
4709 
4710  POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster has %d bands",
4711  rt_raster_get_num_bands(newrast));
4712 
4713  rt_raster_copy_band(newrast, raster, nband - 1, 0);
4714 
4715  POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster now has %d bands",
4716  rt_raster_get_num_bands(newrast));
4717 
4718  if (initexpr)
4719  pfree(initexpr);
4720  rt_raster_destroy(raster);
4721  PG_FREE_IF_COPY(pgraster, 0);
4722 
4723  /* Serialize created raster */
4724  pgrtn = rt_raster_serialize(newrast);
4725  rt_raster_destroy(newrast);
4726  if (NULL == pgrtn) {
4727  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4728  PG_RETURN_NULL();
4729  }
4730 
4731  SET_VARSIZE(pgrtn, pgrtn->size);
4732  PG_RETURN_POINTER(pgrtn);
4733  }
4734 
4739  if (initexpr != NULL && strstr(initexpr, "[rast") == NULL) {
4740  ret = SPI_connect();
4741  if (ret != SPI_OK_CONNECT) {
4742  PG_FREE_IF_COPY(pgraster, 0);
4743  elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4744  PG_RETURN_NULL();
4745  };
4746 
4747  /* Execute the expresion into newval */
4748  ret = SPI_execute(initexpr, FALSE, 0);
4749 
4750  if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
4751 
4752  /* Free memory allocated out of the current context */
4753  if (SPI_tuptable)
4754  SPI_freetuptable(tuptable);
4755  PG_FREE_IF_COPY(pgraster, 0);
4756 
4757  SPI_finish();
4758  elog(ERROR, "RASTER_mapAlgebraExpr: Invalid construction for expression");
4759  PG_RETURN_NULL();
4760  }
4761 
4762  tupdesc = SPI_tuptable->tupdesc;
4763  tuptable = SPI_tuptable;
4764 
4765  tuple = tuptable->vals[0];
4766  newexpr = SPI_getvalue(tuple, tupdesc, 1);
4767  if ( ! newexpr ) {
4768  POSTGIS_RT_DEBUG(3, "Constant expression evaluated to NULL, keeping initvalue");
4769  newval = newinitialvalue;
4770  } else {
4771  newval = atof(newexpr);
4772  }
4773 
4774  SPI_freetuptable(tuptable);
4775 
4776  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: New raster value = %f",
4777  newval);
4778 
4779  SPI_finish();
4780 
4781  skipcomputation = 1;
4782 
4787  if (!hasnodataval) {
4788  newinitialvalue = newval;
4789  skipcomputation = 2;
4790  }
4791 
4792  /* Return the new raster as it will be before computing pixel by pixel */
4793  else if (FLT_NEQ(newval, newinitialvalue)) {
4794  skipcomputation = 2;
4795  }
4796  }
4797 
4802  ret = rt_raster_generate_new_band(newrast, newpixeltype,
4803  newinitialvalue, TRUE, newnodatavalue, 0);
4804 
4809  /*if (initexpr == NULL || skipcomputation == 2) {*/
4810  if (expression == NULL || skipcomputation == 2) {
4811 
4812  /* Free memory */
4813  if (initexpr)
4814  pfree(initexpr);
4815  rt_raster_destroy(raster);
4816  PG_FREE_IF_COPY(pgraster, 0);
4817 
4818  /* Serialize created raster */
4819  pgrtn = rt_raster_serialize(newrast);
4820  rt_raster_destroy(newrast);
4821  if (NULL == pgrtn) {
4822  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4823  PG_RETURN_NULL();
4824  }
4825 
4826  SET_VARSIZE(pgrtn, pgrtn->size);
4827  PG_RETURN_POINTER(pgrtn);
4828  }
4829 
4830  RASTER_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new raster band...");
4831 
4832  /* Get the new raster band */
4833  newband = rt_raster_get_band(newrast, 0);
4834  if ( NULL == newband ) {
4835  elog(NOTICE, "Could not modify band for new raster. Returning new "
4836  "raster with the original band");
4837 
4838  if (initexpr)
4839  pfree(initexpr);
4840  rt_raster_destroy(raster);
4841  PG_FREE_IF_COPY(pgraster, 0);
4842 
4843  /* Serialize created raster */
4844  pgrtn = rt_raster_serialize(newrast);
4845  rt_raster_destroy(newrast);
4846  if (NULL == pgrtn) {
4847  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4848  PG_RETURN_NULL();
4849  }
4850 
4851  SET_VARSIZE(pgrtn, pgrtn->size);
4852  PG_RETURN_POINTER(pgrtn);
4853  }
4854 
4855  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Main computing loop (%d x %d)",
4856  width, height);
4857 
4858  if (initexpr != NULL) {
4859  /* Convert [rast.val] to [rast] */
4860  newexpr = rtpg_strreplace(initexpr, "[rast.val]", "[rast]", NULL);
4861  pfree(initexpr); initexpr=newexpr;
4862 
4863  sprintf(place,"$1");
4864  for (i = 0, j = 1; i < argkwcount; i++) {
4865  len = 0;
4866  newexpr = rtpg_strreplace(initexpr, argkw[i], place, &len);
4867  pfree(initexpr); initexpr=newexpr;
4868  if (len > 0) {
4869  argtype[argcount] = argkwtypes[i];
4870  argcount++;
4871  argpos[i] = j++;
4872 
4873  sprintf(place, "$%d", j);
4874  }
4875  else {
4876  argpos[i] = 0;
4877  }
4878  }
4879 
4880  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: initexpr = %s", initexpr);
4881 
4882  /* define values */
4883  values = (Datum *) palloc(sizeof(Datum) * argcount);
4884  if (values == NULL) {
4885 
4886  SPI_finish();
4887 
4888  rt_raster_destroy(raster);
4889  PG_FREE_IF_COPY(pgraster, 0);
4890  rt_raster_destroy(newrast);
4891 
4892  elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
4893  PG_RETURN_NULL();
4894  }
4895 
4896  /* define nulls */
4897  nulls = (char *)palloc(argcount);
4898  if (nulls == NULL) {
4899 
4900  SPI_finish();
4901 
4902  rt_raster_destroy(raster);
4903  PG_FREE_IF_COPY(pgraster, 0);
4904  rt_raster_destroy(newrast);
4905 
4906  elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
4907  PG_RETURN_NULL();
4908  }
4909 
4910  /* Connect to SPI and prepare the expression */
4911  ret = SPI_connect();
4912  if (ret != SPI_OK_CONNECT) {
4913 
4914  if (initexpr)
4915  pfree(initexpr);
4916  rt_raster_destroy(raster);
4917  PG_FREE_IF_COPY(pgraster, 0);
4918  rt_raster_destroy(newrast);
4919 
4920  elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4921  PG_RETURN_NULL();
4922  };
4923 
4924  /* Type of all arguments is FLOAT8OID */
4925  spi_plan = SPI_prepare(initexpr, argcount, argtype);
4926 
4927  if (spi_plan == NULL) {
4928 
4929  rt_raster_destroy(raster);
4930  PG_FREE_IF_COPY(pgraster, 0);
4931  rt_raster_destroy(newrast);
4932 
4933  SPI_finish();
4934 
4935  pfree(initexpr);
4936 
4937  elog(ERROR, "RASTER_mapAlgebraExpr: Could not prepare expression");
4938  PG_RETURN_NULL();
4939  }
4940  }
4941 
4942  for (x = 0; x < width; x++) {
4943  for(y = 0; y < height; y++) {
4944  ret = rt_band_get_pixel(band, x, y, &r, NULL);
4945 
4950  if (ret == ES_NONE && FLT_NEQ(r, newnodatavalue)) {
4951  if (skipcomputation == 0) {
4952  if (initexpr != NULL) {
4953  /* Reset the null arg flags. */
4954  memset(nulls, 'n', argcount);
4955 
4956  for (i = 0; i < argkwcount; i++) {
4957  idx = argpos[i];
4958  if (idx < 1) continue;
4959  idx--;
4960 
4961  if (i == kX ) {
4962  /* x is 0 based index, but SQL expects 1 based index */
4963  values[idx] = Int32GetDatum(x+1);
4964  nulls[idx] = ' ';
4965  }
4966  else if (i == kY) {
4967  /* y is 0 based index, but SQL expects 1 based index */
4968  values[idx] = Int32GetDatum(y+1);
4969  nulls[idx] = ' ';
4970  }
4971  else if (i == kVAL ) {
4972  values[idx] = Float8GetDatum(r);
4973  nulls[idx] = ' ';
4974  }
4975 
4976  }
4977 
4978  ret = SPI_execute_plan(spi_plan, values, nulls, FALSE, 0);
4979  if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
4980  SPI_processed != 1) {
4981  if (SPI_tuptable)
4982  SPI_freetuptable(tuptable);
4983 
4984  SPI_freeplan(spi_plan);
4985  SPI_finish();
4986 
4987  pfree(values);
4988  pfree(nulls);
4989  pfree(initexpr);
4990 
4991  rt_raster_destroy(raster);
4992  PG_FREE_IF_COPY(pgraster, 0);
4993  rt_raster_destroy(newrast);
4994 
4995  elog(ERROR, "RASTER_mapAlgebraExpr: Error executing prepared plan");
4996 
4997  PG_RETURN_NULL();
4998  }
4999 
5000  tupdesc = SPI_tuptable->tupdesc;
5001  tuptable = SPI_tuptable;
5002 
5003  tuple = tuptable->vals[0];
5004  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
5005  if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
5006  POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) errored, skip setting", x+1,y+1,r);
5007  newval = newinitialvalue;
5008  }
5009  else if ( isnull ) {
5010  POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting", x+1,y+1,r);
5011  newval = newinitialvalue;
5012  } else {
5013  newval = DatumGetFloat8(datum);
5014  }
5015 
5016  SPI_freetuptable(tuptable);
5017  }
5018 
5019  else
5020  newval = newinitialvalue;
5021 
5022  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new value = %f",
5023  newval);
5024  }
5025 
5026 
5027  rt_band_set_pixel(newband, x, y, newval, NULL);
5028  }
5029 
5030  }
5031  }
5032 
5033  if (initexpr != NULL) {
5034  SPI_freeplan(spi_plan);
5035  SPI_finish();
5036 
5037  pfree(values);
5038  pfree(nulls);
5039  pfree(initexpr);
5040  }
5041  else {
5042  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: no SPI cleanup");
5043  }
5044 
5045 
5046  /* The newrast band has been modified */
5047 
5048  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster modified, serializing it.");
5049  /* Serialize created raster */
5050 
5051  rt_raster_destroy(raster);
5052  PG_FREE_IF_COPY(pgraster, 0);
5053 
5054  pgrtn = rt_raster_serialize(newrast);
5055  rt_raster_destroy(newrast);
5056  if (NULL == pgrtn)
5057  PG_RETURN_NULL();
5058 
5059  SET_VARSIZE(pgrtn, pgrtn->size);
5060 
5061  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster serialized");
5062 
5063 
5064  POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraExpr: returning raster");
5065 
5066 
5067  PG_RETURN_POINTER(pgrtn);
5068 }
5069 
5070 /*
5071  * One raster user function MapAlgebra.
5072  */
5074 Datum RASTER_mapAlgebraFct(PG_FUNCTION_ARGS)
5075 {
5076  rt_pgraster *pgraster = NULL;
5077  rt_pgraster *pgrtn = NULL;
5078  rt_raster raster = NULL;
5079  rt_raster newrast = NULL;
5080  rt_band band = NULL;
5081  rt_band newband = NULL;
5082  int x, y, nband, width, height;
5083  double r;
5084  double newnodatavalue = 0.0;
5085  double newinitialvalue = 0.0;
5086  double newval = 0.0;
5087  rt_pixtype newpixeltype;
5088  int ret = -1;
5089  Oid oid;
5090  FmgrInfo cbinfo;
5091  FunctionCallInfoData cbdata;
5092  Datum tmpnewval;
5093  char * strFromText = NULL;
5094  int k = 0;
5095 
5096  POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFct: STARTING...");
5097 
5098  /* Check raster */
5099  if (PG_ARGISNULL(0)) {
5100  elog(WARNING, "Raster is NULL. Returning NULL");
5101  PG_RETURN_NULL();
5102  }
5103 
5104 
5105  /* Deserialize raster */
5106  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5107  raster = rt_raster_deserialize(pgraster, FALSE);
5108  if (NULL == raster) {
5109  PG_FREE_IF_COPY(pgraster, 0);
5110  elog(ERROR, "RASTER_mapAlgebraFct: Could not deserialize raster");
5111  PG_RETURN_NULL();
5112  }
5113 
5114  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting arguments...");
5115 
5116  /* Get the rest of the arguments */
5117 
5118  if (PG_ARGISNULL(1))
5119  nband = 1;
5120  else
5121  nband = PG_GETARG_INT32(1);
5122 
5123  if (nband < 1)
5124  nband = 1;
5125 
5126  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Creating new empty raster...");
5127 
5132  width = rt_raster_get_width(raster);
5133  height = rt_raster_get_height(raster);
5134 
5135  newrast = rt_raster_new(width, height);
5136 
5137  if ( NULL == newrast ) {
5138 
5139  rt_raster_destroy(raster);
5140  PG_FREE_IF_COPY(pgraster, 0);
5141 
5142  elog(ERROR, "RASTER_mapAlgebraFct: Could not create a new raster");
5143  PG_RETURN_NULL();
5144  }
5145 
5146  rt_raster_set_scale(newrast,
5147  rt_raster_get_x_scale(raster),
5148  rt_raster_get_y_scale(raster));
5149 
5150  rt_raster_set_offsets(newrast,
5151  rt_raster_get_x_offset(raster),
5152  rt_raster_get_y_offset(raster));
5153 
5154  rt_raster_set_skews(newrast,
5155  rt_raster_get_x_skew(raster),
5156  rt_raster_get_y_skew(raster));
5157 
5158  rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
5159 
5160 
5165  if (rt_raster_is_empty(newrast))
5166  {
5167  elog(NOTICE, "Raster is empty. Returning an empty raster");
5168  rt_raster_destroy(raster);
5169  PG_FREE_IF_COPY(pgraster, 0);
5170 
5171  pgrtn = rt_raster_serialize(newrast);
5172  rt_raster_destroy(newrast);
5173  if (NULL == pgrtn) {
5174  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5175  PG_RETURN_NULL();
5176  }
5177 
5178  SET_VARSIZE(pgrtn, pgrtn->size);
5179  PG_RETURN_POINTER(pgrtn);
5180  }
5181 
5182  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Getting raster band %d...", nband);
5183 
5188  if (!rt_raster_has_band(raster, nband - 1)) {
5189  elog(NOTICE, "Raster does not have the required band. Returning a raster "
5190  "without a band");
5191  rt_raster_destroy(raster);
5192  PG_FREE_IF_COPY(pgraster, 0);
5193 
5194  pgrtn = rt_raster_serialize(newrast);
5195  rt_raster_destroy(newrast);
5196  if (NULL == pgrtn) {
5197  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5198  PG_RETURN_NULL();
5199  }
5200 
5201  SET_VARSIZE(pgrtn, pgrtn->size);
5202  PG_RETURN_POINTER(pgrtn);
5203  }
5204 
5205  /* Get the raster band */
5206  band = rt_raster_get_band(raster, nband - 1);
5207  if ( NULL == band ) {
5208  elog(NOTICE, "Could not get the required band. Returning a raster "
5209  "without a band");
5210  rt_raster_destroy(raster);
5211  PG_FREE_IF_COPY(pgraster, 0);
5212 
5213  pgrtn = rt_raster_serialize(newrast);
5214  rt_raster_destroy(newrast);
5215  if (NULL == pgrtn) {
5216  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5217  PG_RETURN_NULL();
5218  }
5219 
5220  SET_VARSIZE(pgrtn, pgrtn->size);
5221  PG_RETURN_POINTER(pgrtn);
5222  }
5223 
5224  /*
5225  * Get NODATA value
5226  */
5227  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting NODATA value for band...");
5228 
5229  if (rt_band_get_hasnodata_flag(band)) {
5230  rt_band_get_nodata(band, &newnodatavalue);
5231  }
5232 
5233  else {
5234  newnodatavalue = rt_band_get_min_value(band);
5235  }
5236 
5237  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: NODATA value for band: %f",
5238  newnodatavalue);
5244  newinitialvalue = newnodatavalue;
5245 
5249  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Setting pixeltype...");
5250 
5251  if (PG_ARGISNULL(2)) {
5252  newpixeltype = rt_band_get_pixtype(band);
5253  }
5254 
5255  else {
5256  strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5257  newpixeltype = rt_pixtype_index_from_name(strFromText);
5258  pfree(strFromText);
5259  if (newpixeltype == PT_END)
5260  newpixeltype = rt_band_get_pixtype(band);
5261  }
5262 
5263  if (newpixeltype == PT_END) {
5264 
5265