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