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