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