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