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