PostGIS  3.3.9dev-r@@SVN_REVISION@@
rtpg_mapalgebra.c
Go to the documentation of this file.
1 /*
2  *
3  * WKTRaster - Raster Types for PostGIS
4  * http://trac.osgeo.org/postgis/wiki/WKTRaster
5  *
6  * Copyright (C) 2011-2013 Regents of the University of California
7  * <bkpark@ucdavis.edu>
8  * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
9  * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
10  * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
11  * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
12  * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
13  * Copyright (C) 2013 Nathaniel Hunter Clay <clay.nathaniel@gmail.com>
14  *
15  * This program is free software; you can redistribute it and/or
16  * modify it under the terms of the GNU General Public License
17  * as published by the Free Software Foundation; either version 2
18  * of the License, or (at your option) any later version.
19  *
20  * This program is distributed in the hope that it will be useful,
21  * but WITHOUT ANY WARRANTY; without even the implied warranty of
22  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  * GNU General Public License for more details.
24  *
25  * You should have received a copy of the GNU General Public License
26  * along with this program; if not, write to the Free Software Foundation,
27  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28  *
29  */
30 
31 #include <assert.h>
32 
33 #include <postgres.h> /* for palloc */
34 #include <fmgr.h>
35 #include <funcapi.h>
36 #include <executor/spi.h>
37 #include <utils/lsyscache.h> /* for get_typlenbyvalalign */
38 #include <utils/array.h> /* for ArrayType */
39 #include <utils/builtins.h> /* 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
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] = (Datum)arg->src_pixel[z][0] + 1;
470  i++;
471 
472  _pos[i] = (Datum)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  i = (arg->numraster > 1) ? 1 : 0;
937  break;
938  default:
939  i = 0;
940  break;
941  }
942  /* find first viable band */
943  if (!arg->hasband[i]) {
944  for (i = 0; i < arg->numraster; i++) {
945  if (arg->hasband[i])
946  break;
947  }
948  if (i >= arg->numraster)
949  i = arg->numraster - 1;
950  }
951  band = rt_raster_get_band(arg->raster[i], arg->nband[i]);
952 
953  /* set pixel type if PT_END */
954  if (arg->pixtype == PT_END)
956 
957  /* set hasnodata and nodataval */
958  arg->hasnodata = 1;
961  else
963 
964  POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->pixtype), arg->hasnodata, arg->nodataval);
965 
966  /* init itrset */
967  itrset = palloc(sizeof(struct rt_iterator_t) * arg->numraster);
968  if (itrset == NULL) {
970  elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
971  PG_RETURN_NULL();
972  }
973 
974  /* set itrset */
975  for (i = 0; i < arg->numraster; i++) {
976  itrset[i].raster = arg->raster[i];
977  itrset[i].nband = arg->nband[i];
978  itrset[i].nbnodata = 1;
979  }
980 
981  /* pass everything to iterator */
982  noerr = rt_raster_iterator(
983  itrset, arg->numraster,
984  arg->extenttype, arg->cextent,
985  arg->pixtype,
986  arg->hasnodata, arg->nodataval,
987  arg->distance[0], arg->distance[1],
988  arg->mask,
989  &(arg->callback),
991  &raster
992  );
993 
994  /* cleanup */
995  pfree(itrset);
997 
998  if (noerr != ES_NONE) {
999  elog(ERROR, "RASTER_nMapAlgebra: Could not run raster iterator function");
1000  PG_RETURN_NULL();
1001  }
1002  else if (raster == NULL)
1003  PG_RETURN_NULL();
1004 
1005  pgraster = rt_raster_serialize(raster);
1007 
1008  POSTGIS_RT_DEBUG(3, "Finished");
1009 
1010  if (!pgraster)
1011  PG_RETURN_NULL();
1012 
1013  SET_VARSIZE(pgraster, pgraster->size);
1014  PG_RETURN_POINTER(pgraster);
1015 }
1016 
1017 /* ---------------------------------------------------------------- */
1018 /* expression ST_MapAlgebra for n rasters */
1019 /* ---------------------------------------------------------------- */
1020 
1021 typedef struct {
1023 
1024  struct {
1025  SPIPlanPtr spi_plan;
1026  uint32_t spi_argcount;
1027  uint8_t *spi_argpos;
1028 
1029  int hasval;
1030  double val;
1031  } expr[3];
1032 
1033  struct {
1034  int hasval;
1035  double val;
1036  } nodatanodata;
1037 
1038  struct {
1039  int count;
1040  char **val;
1041  } kw;
1042 
1044 
1048 
1050 };
1051 
1053  rtpg_nmapalgebraexpr_arg arg = NULL;
1054  int i = 0;
1055 
1056  arg = palloc(sizeof(struct rtpg_nmapalgebraexpr_arg_t));
1057  if (arg == NULL) {
1058  elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arguments");
1059  return NULL;
1060  }
1061 
1063  if (arg->bandarg == NULL) {
1064  elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arg->bandarg");
1065  return NULL;
1066  }
1067 
1068  arg->callback.kw.count = cnt;
1069  arg->callback.kw.val = kw;
1070 
1071  arg->callback.exprcount = 3;
1072  for (i = 0; i < arg->callback.exprcount; i++) {
1073  arg->callback.expr[i].spi_plan = NULL;
1074  arg->callback.expr[i].spi_argcount = 0;
1075  arg->callback.expr[i].spi_argpos = palloc(cnt * sizeof(uint8_t));
1076  if (arg->callback.expr[i].spi_argpos == NULL) {
1077  elog(ERROR, "rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for spi_argpos");
1078  return NULL;
1079  }
1080  memset(arg->callback.expr[i].spi_argpos, 0, sizeof(uint8_t) * cnt);
1081  arg->callback.expr[i].hasval = 0;
1082  arg->callback.expr[i].val = 0;
1083  }
1084 
1085  arg->callback.nodatanodata.hasval = 0;
1086  arg->callback.nodatanodata.val = 0;
1087 
1088  return arg;
1089 }
1090 
1092  int i = 0;
1093 
1095 
1096  for (i = 0; i < arg->callback.exprcount; i++) {
1097  if (arg->callback.expr[i].spi_plan)
1098  SPI_freeplan(arg->callback.expr[i].spi_plan);
1099  if (arg->callback.kw.count)
1100  pfree(arg->callback.expr[i].spi_argpos);
1101  }
1102 
1103  pfree(arg);
1104 }
1105 
1107  rt_iterator_arg arg, void *userarg,
1108  double *value, int *nodata
1109 ) {
1111  SPIPlanPtr plan = NULL;
1112  int i = 0;
1113  uint8_t id = 0;
1114 
1115  if (arg == NULL)
1116  return 0;
1117 
1118  *value = 0;
1119  *nodata = 0;
1120 
1121  /* 2 raster */
1122  if (arg->rasters > 1) {
1123  /* nodata1 = 1 AND nodata2 = 1, nodatanodataval */
1124  if (arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1125  if (callback->nodatanodata.hasval)
1126  *value = callback->nodatanodata.val;
1127  else
1128  *nodata = 1;
1129  }
1130  /* nodata1 = 1 AND nodata2 != 1, nodata1expr */
1131  else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0]) {
1132  id = 1;
1133  if (callback->expr[id].hasval)
1134  *value = callback->expr[id].val;
1135  else if (callback->expr[id].spi_plan)
1136  plan = callback->expr[id].spi_plan;
1137  else
1138  *nodata = 1;
1139  }
1140  /* nodata1 != 1 AND nodata2 = 1, nodata2expr */
1141  else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1142  id = 2;
1143  if (callback->expr[id].hasval)
1144  *value = callback->expr[id].val;
1145  else if (callback->expr[id].spi_plan)
1146  plan = callback->expr[id].spi_plan;
1147  else
1148  *nodata = 1;
1149  }
1150  /* expression */
1151  else {
1152  id = 0;
1153  if (callback->expr[id].hasval)
1154  *value = callback->expr[id].val;
1155  else if (callback->expr[id].spi_plan)
1156  plan = callback->expr[id].spi_plan;
1157  else {
1158  if (callback->nodatanodata.hasval)
1159  *value = callback->nodatanodata.val;
1160  else
1161  *nodata = 1;
1162  }
1163  }
1164  }
1165  /* 1 raster */
1166  else {
1167  /* nodata = 1, nodata1expr */
1168  if (arg->nodata[0][0][0]) {
1169  id = 1;
1170  if (callback->expr[id].hasval)
1171  *value = callback->expr[id].val;
1172  else if (callback->expr[id].spi_plan)
1173  plan = callback->expr[id].spi_plan;
1174  else
1175  *nodata = 1;
1176  }
1177  /* expression */
1178  else {
1179  id = 0;
1180  if (callback->expr[id].hasval)
1181  *value = callback->expr[id].val;
1182  else if (callback->expr[id].spi_plan)
1183  plan = callback->expr[id].spi_plan;
1184  else {
1185  /* see if nodata1expr is available */
1186  id = 1;
1187  if (callback->expr[id].hasval)
1188  *value = callback->expr[id].val;
1189  else if (callback->expr[id].spi_plan)
1190  plan = callback->expr[id].spi_plan;
1191  else
1192  *nodata = 1;
1193  }
1194  }
1195  }
1196 
1197  /* run prepared plan */
1198  if (plan != NULL) {
1199  Datum values[12];
1200  char nulls[12];
1201  int err = 0;
1202 
1203  TupleDesc tupdesc;
1204  SPITupleTable *tuptable = NULL;
1205  HeapTuple tuple;
1206  Datum datum;
1207  bool isnull = FALSE;
1208 
1209  POSTGIS_RT_DEBUGF(4, "Running plan %d", id);
1210 
1211  /* init values and nulls */
1212  memset(values, (Datum) NULL, sizeof(Datum) * callback->kw.count);
1213  memset(nulls, FALSE, sizeof(char) * callback->kw.count);
1214 
1215  if (callback->expr[id].spi_argcount) {
1216  int idx = 0;
1217 
1218  for (i = 0; i < callback->kw.count; i++) {
1219  idx = callback->expr[id].spi_argpos[i];
1220  if (idx < 1) continue;
1221  idx--; /* 1-based now 0-based */
1222 
1223  if (arg->rasters == 1 && i > 7) {
1224  elog(ERROR, "rtpg_nmapalgebraexpr_callback: rast2 argument specified in single-raster invocation");
1225  return 0;
1226  }
1227 
1228  switch (i) {
1229  /* [rast.x] */
1230  case 0:
1231  values[idx] = Int32GetDatum(arg->src_pixel[0][0] + 1);
1232  break;
1233  /* [rast.y] */
1234  case 1:
1235  values[idx] = Int32GetDatum(arg->src_pixel[0][1] + 1);
1236  break;
1237  /* [rast.val] */
1238  case 2:
1239  /* [rast] */
1240  case 3:
1241  if (!arg->nodata[0][0][0])
1242  values[idx] = Float8GetDatum(arg->values[0][0][0]);
1243  else
1244  nulls[idx] = TRUE;
1245  break;
1246 
1247  /* [rast1.x] */
1248  case 4:
1249  values[idx] = Int32GetDatum(arg->src_pixel[0][0] + 1);
1250  break;
1251  /* [rast1.y] */
1252  case 5:
1253  values[idx] = Int32GetDatum(arg->src_pixel[0][1] + 1);
1254  break;
1255  /* [rast1.val] */
1256  case 6:
1257  /* [rast1] */
1258  case 7:
1259  if (!arg->nodata[0][0][0])
1260  values[idx] = Float8GetDatum(arg->values[0][0][0]);
1261  else
1262  nulls[idx] = TRUE;
1263  break;
1264 
1265  /* [rast2.x] */
1266  case 8:
1267  values[idx] = Int32GetDatum(arg->src_pixel[1][0] + 1);
1268  break;
1269  /* [rast2.y] */
1270  case 9:
1271  values[idx] = Int32GetDatum(arg->src_pixel[1][1] + 1);
1272  break;
1273  /* [rast2.val] */
1274  case 10:
1275  /* [rast2] */
1276  case 11:
1277  if (!arg->nodata[1][0][0])
1278  values[idx] = Float8GetDatum(arg->values[1][0][0]);
1279  else
1280  nulls[idx] = TRUE;
1281  break;
1282  }
1283 
1284  }
1285  }
1286 
1287  /* run prepared plan */
1288  err = SPI_execute_plan(plan, values, nulls, TRUE, 1);
1289  if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1290  elog(ERROR, "rtpg_nmapalgebraexpr_callback: Unexpected error when running prepared statement %d", id);
1291  return 0;
1292  }
1293 
1294  /* get output of prepared plan */
1295  tupdesc = SPI_tuptable->tupdesc;
1296  tuptable = SPI_tuptable;
1297  tuple = tuptable->vals[0];
1298 
1299  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1300  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1301  if (SPI_tuptable) SPI_freetuptable(tuptable);
1302  elog(ERROR, "rtpg_nmapalgebraexpr_callback: Could not get result of prepared statement %d", id);
1303  return 0;
1304  }
1305 
1306  if (!isnull) {
1307  *value = DatumGetFloat8(datum);
1308  POSTGIS_RT_DEBUG(4, "Getting value from Datum");
1309  }
1310  else {
1311  /* 2 raster, check nodatanodataval */
1312  if (arg->rasters > 1) {
1313  if (callback->nodatanodata.hasval)
1314  *value = callback->nodatanodata.val;
1315  else
1316  *nodata = 1;
1317  }
1318  /* 1 raster, check nodataval */
1319  else {
1320  if (callback->expr[1].hasval)
1321  *value = callback->expr[1].val;
1322  else
1323  *nodata = 1;
1324  }
1325  }
1326 
1327  if (SPI_tuptable) SPI_freetuptable(tuptable);
1328  }
1329 
1330  POSTGIS_RT_DEBUGF(4, "(value, nodata) = (%f, %d)", *value, *nodata);
1331  return 1;
1332 }
1333 
1335 Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS)
1336 {
1337  MemoryContext mainMemCtx = CurrentMemoryContext;
1338  rtpg_nmapalgebraexpr_arg arg = NULL;
1339  rt_iterator itrset;
1340  uint16_t exprpos[3] = {1, 4, 5};
1341 
1342  int i = 0;
1343  int j = 0;
1344  int k = 0;
1345 
1346  int numraster = 0;
1347  int err = 0;
1348  int allnull = 0;
1349  int allempty = 0;
1350  int noband = 0;
1351  int len = 0;
1352 
1353  TupleDesc tupdesc;
1354  SPITupleTable *tuptable = NULL;
1355  HeapTuple tuple;
1356  Datum datum;
1357  bool isnull = FALSE;
1358 
1359  rt_raster raster = NULL;
1360  rt_band band = NULL;
1361  rt_pgraster *pgraster = NULL;
1362 
1363  const int argkwcount = 12;
1364  char *argkw[] = {
1365  "[rast.x]",
1366  "[rast.y]",
1367  "[rast.val]",
1368  "[rast]",
1369  "[rast1.x]",
1370  "[rast1.y]",
1371  "[rast1.val]",
1372  "[rast1]",
1373  "[rast2.x]",
1374  "[rast2.y]",
1375  "[rast2.val]",
1376  "[rast2]"
1377  };
1378 
1379  if (PG_ARGISNULL(0))
1380  PG_RETURN_NULL();
1381 
1382  /* init argument struct */
1383  arg = rtpg_nmapalgebraexpr_arg_init(argkwcount, argkw);
1384  if (arg == NULL) {
1385  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not initialize argument structure");
1386  PG_RETURN_NULL();
1387  }
1388 
1389  /* let helper function process rastbandarg (0) */
1390  if (!rtpg_nmapalgebra_rastbandarg_process(arg->bandarg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
1392  elog(ERROR, "RASTER_nMapAlgebra: Could not process rastbandarg");
1393  PG_RETURN_NULL();
1394  }
1395 
1396  POSTGIS_RT_DEBUGF(4, "allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
1397 
1398  /* all rasters are NULL, return NULL */
1399  if (allnull == arg->bandarg->numraster) {
1400  elog(NOTICE, "All input rasters are NULL. Returning NULL");
1402  PG_RETURN_NULL();
1403  }
1404 
1405  /* only work on one or two rasters */
1406  if (arg->bandarg->numraster > 1)
1407  numraster = 2;
1408  else
1409  numraster = 1;
1410 
1411  /* pixel type (2) */
1412  if (!PG_ARGISNULL(2)) {
1413  char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
1414 
1415  /* Get the pixel type index */
1416  arg->bandarg->pixtype = rt_pixtype_index_from_name(pixtypename);
1417  if (arg->bandarg->pixtype == PT_END) {
1419  elog(ERROR, "RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
1420  PG_RETURN_NULL();
1421  }
1422  }
1423  POSTGIS_RT_DEBUGF(4, "pixeltype: %d", arg->bandarg->pixtype);
1424 
1425  /* extent type (3) */
1426  if (!PG_ARGISNULL(3)) {
1427  char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(3))));
1428  arg->bandarg->extenttype = rt_util_extent_type(extenttypename);
1429  }
1430 
1431  if (arg->bandarg->extenttype == ET_CUSTOM) {
1432  if (numraster < 2) {
1433  elog(NOTICE, "CUSTOM extent type not supported. Defaulting to FIRST");
1434  arg->bandarg->extenttype = ET_FIRST;
1435  }
1436  else {
1437  elog(NOTICE, "CUSTOM extent type not supported. Defaulting to INTERSECTION");
1439  }
1440  }
1441  else if (numraster < 2)
1442  arg->bandarg->extenttype = ET_FIRST;
1443 
1444  POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->bandarg->extenttype);
1445 
1446  /* nodatanodataval (6) */
1447  if (!PG_ARGISNULL(6)) {
1448  arg->callback.nodatanodata.hasval = 1;
1449  arg->callback.nodatanodata.val = PG_GETARG_FLOAT8(6);
1450  }
1451 
1452  err = 0;
1453  /* all rasters are empty, return empty raster */
1454  if (allempty == arg->bandarg->numraster) {
1455  elog(NOTICE, "All input rasters are empty. Returning empty raster");
1456  err = 1;
1457  }
1458  /* all rasters don't have indicated band, return empty raster */
1459  else if (noband == arg->bandarg->numraster) {
1460  elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
1461  err = 1;
1462  }
1463  if (err) {
1465 
1466  raster = rt_raster_new(0, 0);
1467  if (raster == NULL) {
1468  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create empty raster");
1469  PG_RETURN_NULL();
1470  }
1471 
1472  pgraster = rt_raster_serialize(raster);
1474  if (!pgraster) PG_RETURN_NULL();
1475 
1476  SET_VARSIZE(pgraster, pgraster->size);
1477  PG_RETURN_POINTER(pgraster);
1478  }
1479 
1480  /* connect SPI */
1481  if (SPI_connect() != SPI_OK_CONNECT) {
1483  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
1484  PG_RETURN_NULL();
1485  }
1486 
1487  /*
1488  process expressions
1489 
1490  exprpos elements are:
1491  1 - expression => spi_plan[0]
1492  4 - nodata1expr => spi_plan[1]
1493  5 - nodata2expr => spi_plan[2]
1494  */
1495  for (i = 0; i < arg->callback.exprcount; i++) {
1496  char *expr = NULL;
1497  char *tmp = NULL;
1498  char *sql = NULL;
1499  char place[12] = "$1";
1500 
1501  if (PG_ARGISNULL(exprpos[i]))
1502  continue;
1503 
1504  expr = text_to_cstring(PG_GETARG_TEXT_P(exprpos[i]));
1505  POSTGIS_RT_DEBUGF(3, "raw expr of argument #%d: %s", exprpos[i], expr);
1506 
1507  for (j = 0, k = 1; j < argkwcount; j++) {
1508  /* attempt to replace keyword with placeholder */
1509  len = 0;
1510  tmp = rtpg_strreplace(expr, argkw[j], place, &len);
1511  pfree(expr);
1512  expr = tmp;
1513 
1514  if (len) {
1515  POSTGIS_RT_DEBUGF(4, "kw #%d (%s) at pos $%d", j, argkw[j], k);
1516  arg->callback.expr[i].spi_argcount++;
1517  arg->callback.expr[i].spi_argpos[j] = k++;
1518 
1519  sprintf(place, "$%d", k);
1520  }
1521  else
1522  arg->callback.expr[i].spi_argpos[j] = 0;
1523  }
1524 
1525  len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision");
1526  sql = (char *) palloc(len + 1);
1527  if (sql == NULL) {
1529  SPI_finish();
1530  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
1531  PG_RETURN_NULL();
1532  }
1533 
1534  memcpy(sql, "SELECT (", strlen("SELECT ("));
1535  memcpy(sql + strlen("SELECT ("), expr, strlen(expr));
1536  memcpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision"));
1537  sql[len] = '\0';
1538 
1539  POSTGIS_RT_DEBUGF(3, "sql #%d: %s", exprpos[i], sql);
1540 
1541  /* prepared plan */
1542  if (arg->callback.expr[i].spi_argcount) {
1543  Oid *argtype = (Oid *) palloc(arg->callback.expr[i].spi_argcount * sizeof(Oid));
1544  POSTGIS_RT_DEBUGF(3, "expression parameter %d is a prepared plan", exprpos[i]);
1545  if (argtype == NULL) {
1546  pfree(sql);
1548  SPI_finish();
1549  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
1550  PG_RETURN_NULL();
1551  }
1552 
1553  /* specify datatypes of parameters */
1554  for (j = 0, k = 0; j < argkwcount; j++) {
1555  if (arg->callback.expr[i].spi_argpos[j] < 1) continue;
1556 
1557  /* positions are INT4 */
1558  if (
1559  (strstr(argkw[j], "[rast.x]") != NULL) ||
1560  (strstr(argkw[j], "[rast.y]") != NULL) ||
1561  (strstr(argkw[j], "[rast1.x]") != NULL) ||
1562  (strstr(argkw[j], "[rast1.y]") != NULL) ||
1563  (strstr(argkw[j], "[rast2.x]") != NULL) ||
1564  (strstr(argkw[j], "[rast2.y]") != NULL)
1565  )
1566  argtype[k] = INT4OID;
1567  /* everything else is FLOAT8 */
1568  else
1569  argtype[k] = FLOAT8OID;
1570 
1571  k++;
1572  }
1573 
1574  arg->callback.expr[i].spi_plan = SPI_prepare(sql, arg->callback.expr[i].spi_argcount, argtype);
1575  pfree(argtype);
1576  pfree(sql);
1577 
1578  if (arg->callback.expr[i].spi_plan == NULL) {
1580  SPI_finish();
1581  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
1582  PG_RETURN_NULL();
1583  }
1584  }
1585  /* no args, just execute query */
1586  else {
1587  POSTGIS_RT_DEBUGF(3, "expression parameter %d has no args, simply executing", exprpos[i]);
1588  err = SPI_execute(sql, TRUE, 0);
1589  pfree(sql);
1590 
1591  if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1593  SPI_finish();
1594  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
1595  PG_RETURN_NULL();
1596  }
1597 
1598  /* get output of prepared plan */
1599  tupdesc = SPI_tuptable->tupdesc;
1600  tuptable = SPI_tuptable;
1601  tuple = tuptable->vals[0];
1602 
1603  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1604  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1605  if (SPI_tuptable) SPI_freetuptable(tuptable);
1607  SPI_finish();
1608  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
1609  PG_RETURN_NULL();
1610  }
1611 
1612  if (!isnull) {
1613  arg->callback.expr[i].hasval = 1;
1614  arg->callback.expr[i].val = DatumGetFloat8(datum);
1615  }
1616 
1617  if (SPI_tuptable) SPI_freetuptable(tuptable);
1618  }
1619  }
1620 
1621  /* determine nodataval and possibly pixtype */
1622  /* band to check */
1623  switch (arg->bandarg->extenttype) {
1624  case ET_LAST:
1625  case ET_SECOND:
1626  if (numraster > 1)
1627  i = 1;
1628  else
1629  i = 0;
1630  break;
1631  default:
1632  i = 0;
1633  break;
1634  }
1635  /* find first viable band */
1636  if (!arg->bandarg->hasband[i]) {
1637  for (i = 0; i < numraster; i++) {
1638  if (arg->bandarg->hasband[i])
1639  break;
1640  }
1641  if (i >= numraster)
1642  i = numraster - 1;
1643  }
1644  band = rt_raster_get_band(arg->bandarg->raster[i], arg->bandarg->nband[i]);
1645 
1646  /* set pixel type if PT_END */
1647  if (arg->bandarg->pixtype == PT_END)
1649 
1650  /* set hasnodata and nodataval */
1651  arg->bandarg->hasnodata = 1;
1654  else
1656 
1657  POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->bandarg->pixtype), arg->bandarg->hasnodata, arg->bandarg->nodataval);
1658 
1659  /* init itrset */
1660  itrset = palloc(sizeof(struct rt_iterator_t) * numraster);
1661  if (itrset == NULL) {
1663  SPI_finish();
1664  elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
1665  PG_RETURN_NULL();
1666  }
1667 
1668  /* set itrset */
1669  for (i = 0; i < numraster; i++) {
1670  itrset[i].raster = arg->bandarg->raster[i];
1671  itrset[i].nband = arg->bandarg->nband[i];
1672  itrset[i].nbnodata = 1;
1673  }
1674 
1675  /* pass everything to iterator */
1676  err = rt_raster_iterator(
1677  itrset, numraster,
1678  arg->bandarg->extenttype, arg->bandarg->cextent,
1679  arg->bandarg->pixtype,
1680  arg->bandarg->hasnodata, arg->bandarg->nodataval,
1681  0, 0,
1682  NULL,
1683  &(arg->callback),
1685  &raster
1686  );
1687 
1688  pfree(itrset);
1690 
1691  if (err != ES_NONE) {
1692  SPI_finish();
1693  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not run raster iterator function");
1694  PG_RETURN_NULL();
1695  }
1696  else if (raster == NULL) {
1697  SPI_finish();
1698  PG_RETURN_NULL();
1699  }
1700 
1701  /* switch to prior memory context to ensure memory allocated in correct context */
1702  MemoryContextSwitchTo(mainMemCtx);
1703 
1704  pgraster = rt_raster_serialize(raster);
1706 
1707  /* finish SPI */
1708  SPI_finish();
1709 
1710  if (!pgraster)
1711  PG_RETURN_NULL();
1712 
1713  SET_VARSIZE(pgraster, pgraster->size);
1714  PG_RETURN_POINTER(pgraster);
1715 }
1716 
1717 /* ---------------------------------------------------------------- */
1718 /* ST_Union aggregate functions */
1719 /* ---------------------------------------------------------------- */
1720 
1721 typedef enum {
1722  UT_LAST = 0,
1729  UT_RANGE
1731 
1732 /* internal function translating text of UNION type to enum */
1734  assert(cutype && strlen(cutype) > 0);
1735 
1736  if (strcmp(cutype, "LAST") == 0)
1737  return UT_LAST;
1738  else if (strcmp(cutype, "FIRST") == 0)
1739  return UT_FIRST;
1740  else if (strcmp(cutype, "MIN") == 0)
1741  return UT_MIN;
1742  else if (strcmp(cutype, "MAX") == 0)
1743  return UT_MAX;
1744  else if (strcmp(cutype, "COUNT") == 0)
1745  return UT_COUNT;
1746  else if (strcmp(cutype, "SUM") == 0)
1747  return UT_SUM;
1748  else if (strcmp(cutype, "MEAN") == 0)
1749  return UT_MEAN;
1750  else if (strcmp(cutype, "RANGE") == 0)
1751  return UT_RANGE;
1752 
1753  return UT_LAST;
1754 }
1755 
1758  int nband; /* source raster's band index, 0-based */
1760 
1763 };
1764 
1767  int numband; /* number of bandargs */
1769 };
1770 
1772  int i = 0;
1773  int j = 0;
1774  int k = 0;
1775 
1776  if (arg->bandarg != NULL) {
1777  for (i = 0; i < arg->numband; i++) {
1778  if (!arg->bandarg[i].numraster)
1779  continue;
1780 
1781  for (j = 0; j < arg->bandarg[i].numraster; j++) {
1782  if (arg->bandarg[i].raster[j] == NULL)
1783  continue;
1784 
1785  for (k = rt_raster_get_num_bands(arg->bandarg[i].raster[j]) - 1; k >= 0; k--)
1787  rt_raster_destroy(arg->bandarg[i].raster[j]);
1788  }
1789 
1790  pfree(arg->bandarg[i].raster);
1791  }
1792 
1793  pfree(arg->bandarg);
1794  }
1795 
1796  pfree(arg);
1797 }
1798 
1800  rt_iterator_arg arg, void *userarg,
1801  double *value, int *nodata
1802 ) {
1803  rtpg_union_type *utype = (rtpg_union_type *) userarg;
1804 
1805  if (arg == NULL)
1806  return 0;
1807 
1808  if (
1809  arg->rasters != 2 ||
1810  arg->rows != 1 ||
1811  arg->columns != 1
1812  ) {
1813  elog(ERROR, "rtpg_union_callback: Invalid arguments passed to callback");
1814  return 0;
1815  }
1816 
1817  *value = 0;
1818  *nodata = 0;
1819 
1820  /* handle NODATA situations except for COUNT, which is a special case */
1821  if (*utype != UT_COUNT) {
1822  /* both NODATA */
1823  if (arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1824  *nodata = 1;
1825  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1826  return 1;
1827  }
1828  /* second NODATA */
1829  else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0]) {
1830  *value = arg->values[0][0][0];
1831  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1832  return 1;
1833  }
1834  /* first NODATA */
1835  else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0]) {
1836  *value = arg->values[1][0][0];
1837  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1838  return 1;
1839  }
1840  }
1841 
1842  switch (*utype) {
1843  case UT_FIRST:
1844  *value = arg->values[0][0][0];
1845  break;
1846  case UT_MIN:
1847  if (arg->values[0][0][0] < arg->values[1][0][0])
1848  *value = arg->values[0][0][0];
1849  else
1850  *value = arg->values[1][0][0];
1851  break;
1852  case UT_MAX:
1853  if (arg->values[0][0][0] > arg->values[1][0][0])
1854  *value = arg->values[0][0][0];
1855  else
1856  *value = arg->values[1][0][0];
1857  break;
1858  case UT_COUNT:
1859  /* both NODATA */
1860  if (arg->nodata[0][0][0] && arg->nodata[1][0][0])
1861  *value = 0;
1862  /* second NODATA */
1863  else if (!arg->nodata[0][0][0] && arg->nodata[1][0][0])
1864  *value = arg->values[0][0][0];
1865  /* first NODATA */
1866  else if (arg->nodata[0][0][0] && !arg->nodata[1][0][0])
1867  *value = 1;
1868  /* has value, increment */
1869  else
1870  *value = arg->values[0][0][0] + 1;
1871  break;
1872  case UT_SUM:
1873  *value = arg->values[0][0][0] + arg->values[1][0][0];
1874  break;
1875  case UT_MEAN:
1876  case UT_RANGE:
1877  break;
1878  case UT_LAST:
1879  default:
1880  *value = arg->values[1][0][0];
1881  break;
1882  }
1883 
1884  POSTGIS_RT_DEBUGF(4, "value, nodata = %f, %d", *value, *nodata);
1885 
1886 
1887  return 1;
1888 }
1889 
1891  rt_iterator_arg arg, void *userarg,
1892  double *value, int *nodata
1893 ) {
1894  if (arg == NULL)
1895  return 0;
1896 
1897  if (
1898  arg->rasters != 2 ||
1899  arg->rows != 1 ||
1900  arg->columns != 1
1901  ) {
1902  elog(ERROR, "rtpg_union_mean_callback: Invalid arguments passed to callback");
1903  return 0;
1904  }
1905 
1906  *value = 0;
1907  *nodata = 1;
1908 
1909  POSTGIS_RT_DEBUGF(4, "rast0: %f %d", arg->values[0][0][0], arg->nodata[0][0][0]);
1910  POSTGIS_RT_DEBUGF(4, "rast1: %f %d", arg->values[1][0][0], arg->nodata[1][0][0]);
1911 
1912  if (!arg->nodata[0][0][0] && FLT_NEQ(arg->values[0][0][0], 0.0) && !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 = MemoryContextAlloc(aggcontext, 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  /* For Windowing functions, it is important to leave */
2941  /* the state intact, knowing that the aggcontext will be */
2942  /* freed by PgSQL when the statement is complete. */
2943  /* https://trac.osgeo.org/postgis/ticket/4770 */
2944  // pfree(itrset);
2945  // rtpg_union_arg_destroy(iwr);
2946 
2947  if (!_rtn) PG_RETURN_NULL();
2948 
2949  pgraster = rt_raster_serialize(_rtn);
2950  rt_raster_destroy(_rtn);
2951 
2952  POSTGIS_RT_DEBUG(3, "Finished");
2953 
2954  if (!pgraster)
2955  PG_RETURN_NULL();
2956 
2957  SET_VARSIZE(pgraster, pgraster->size);
2958  PG_RETURN_POINTER(pgraster);
2959 }
2960 
2961 /* ---------------------------------------------------------------- */
2962 /* Clip raster with geometry */
2963 /* ---------------------------------------------------------------- */
2964 
2967  int nband; /* band index */
2968  int hasnodata; /* is there a user-specified NODATA? */
2969  double nodataval; /* user-specified NODATA */
2970 };
2971 
2977 
2978  int numbands; /* number of bandargs */
2980 };
2981 
2983  rtpg_clip_arg arg = NULL;
2984 
2985  arg = palloc(sizeof(struct rtpg_clip_arg_t));
2986  if (arg == NULL) {
2987  elog(ERROR, "rtpg_clip_arg_init: Could not allocate memory for function arguments");
2988  return NULL;
2989  }
2990 
2991  arg->extenttype = ET_INTERSECTION;
2992  arg->raster = NULL;
2993  arg->mask = NULL;
2994  arg->numbands = 0;
2995  arg->band = NULL;
2996 
2997  return arg;
2998 }
2999 
3001  if (arg->band != NULL)
3002  pfree(arg->band);
3003 
3004  if (arg->raster != NULL)
3005  rt_raster_destroy(arg->raster);
3006  if (arg->mask != NULL)
3007  rt_raster_destroy(arg->mask);
3008 
3009  pfree(arg);
3010 }
3011 
3013 Datum RASTER_clip(PG_FUNCTION_ARGS)
3014 {
3015  rt_pgraster *pgraster = NULL;
3016  LWGEOM *rastgeom = NULL;
3017  double gt[6] = {0};
3018  int32_t srid = SRID_UNKNOWN;
3019 
3020  rt_pgraster *pgrtn = NULL;
3021  rt_raster rtn = NULL;
3022 
3023  GSERIALIZED *gser = NULL;
3024  LWGEOM *geom = NULL;
3025  lwvarlena_t *wkb = NULL;
3026 
3027  ArrayType *array;
3028  Oid etype;
3029  Datum *e;
3030  bool *nulls;
3031 
3032  int16 typlen;
3033  bool typbyval;
3034  char typalign;
3035 
3036  int i = 0;
3037  int j = 0;
3038  int k = 0;
3039  rtpg_clip_arg arg = NULL;
3040  LWGEOM *tmpgeom = NULL;
3041 
3042  rt_pixtype pixtype;
3043  int hasnodata;
3044  double nodataval;
3045 
3046  double offset[4] = {0.};
3047  int input_x = 0;
3048  int input_y = 0;
3049  int mask_x = 0;
3050  int mask_y = 0;
3051  int x = 0;
3052  int y = 0;
3053  int width = 0;
3054  int height = 0;
3055  int mask_width = 0;
3056  int mask_height = 0;
3057  rt_band input_band = NULL;
3058  rt_band mask_band = NULL;
3059  rt_band output_band = NULL;
3060  double value;
3061  int isnodata;
3062 
3063  POSTGIS_RT_DEBUG(3, "Starting...");
3064 
3065  /* raster or geometry is NULL, return NULL */
3066  if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3067  PG_RETURN_NULL();
3068 
3069  /* init arg */
3070  arg = rtpg_clip_arg_init();
3071  if (arg == NULL) {
3072  elog(ERROR, "RASTER_clip: Could not initialize argument structure");
3073  PG_RETURN_NULL();
3074  }
3075 
3076  /* raster (0) */
3077  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3078 
3079  /* Get raster object */
3080  arg->raster = rt_raster_deserialize(pgraster, FALSE);
3081  if (arg->raster == NULL) {
3082  rtpg_clip_arg_destroy(arg);
3083  PG_FREE_IF_COPY(pgraster, 0);
3084  elog(ERROR, "RASTER_clip: Could not deserialize raster");
3085  PG_RETURN_NULL();
3086  }
3087 
3088  /* raster is empty, return empty raster */
3089  if (rt_raster_is_empty(arg->raster) || rt_raster_get_num_bands(arg->raster) == 0) {
3090  elog(NOTICE, "Input raster is empty or has no bands. Returning empty raster");
3091 
3092  rtpg_clip_arg_destroy(arg);
3093  PG_FREE_IF_COPY(pgraster, 0);
3094 
3095  rtn = rt_raster_new(0, 0);
3096  if (rtn == NULL) {
3097  elog(ERROR, "RASTER_clip: Could not create empty raster");
3098  PG_RETURN_NULL();
3099  }
3100 
3101  pgrtn = rt_raster_serialize(rtn);
3102  rt_raster_destroy(rtn);
3103  if (NULL == pgrtn)
3104  PG_RETURN_NULL();
3105 
3106  SET_VARSIZE(pgrtn, pgrtn->size);
3107  PG_RETURN_POINTER(pgrtn);
3108  }
3109 
3110  /* metadata */
3112  srid = clamp_srid(rt_raster_get_srid(arg->raster));
3113 
3114  /* geometry (2) */
3115  gser = PG_GETARG_GSERIALIZED_P(2);
3116  geom = lwgeom_from_gserialized(gser);
3117 
3118  /* Get a 2D version of the geometry if necessary */
3119  if (lwgeom_ndims(geom) > 2) {
3120  LWGEOM *geom2d = lwgeom_force_2d(geom);
3121  lwgeom_free(geom);
3122  geom = geom2d;
3123  }
3124 
3125  /* check that SRIDs match */
3126  if (srid != clamp_srid(gserialized_get_srid(gser))) {
3127  elog(NOTICE, "Geometry provided does not have the same SRID as the raster. Returning NULL");
3128 
3129  rtpg_clip_arg_destroy(arg);
3130  PG_FREE_IF_COPY(pgraster, 0);
3131  lwgeom_free(geom);
3132  PG_FREE_IF_COPY(gser, 2);
3133 
3134  PG_RETURN_NULL();
3135  }
3136 
3137  /* crop (4) */
3138  if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3139  arg->extenttype = ET_FIRST;
3140 
3141  /* get intersection geometry of input raster and input geometry */
3142  if (rt_raster_get_convex_hull(arg->raster, &rastgeom) != ES_NONE) {
3143 
3144  rtpg_clip_arg_destroy(arg);
3145  PG_FREE_IF_COPY(pgraster, 0);
3146  lwgeom_free(geom);
3147  PG_FREE_IF_COPY(gser, 2);
3148 
3149  elog(ERROR, "RASTER_clip: Could not get convex hull of raster");
3150  PG_RETURN_NULL();
3151  }
3152 
3153  tmpgeom = lwgeom_intersection(rastgeom, geom);
3154  lwgeom_free(rastgeom);
3155  lwgeom_free(geom);
3156  PG_FREE_IF_COPY(gser, 2);
3157  geom = tmpgeom;
3158 
3159  /* intersection is empty AND extent type is INTERSECTION, return empty */
3160  if (lwgeom_is_empty(geom) && arg->extenttype == ET_INTERSECTION) {
3161  elog(NOTICE, "The input raster and input geometry do not intersect. Returning empty raster");
3162 
3163  rtpg_clip_arg_destroy(arg);
3164  PG_FREE_IF_COPY(pgraster, 0);
3165  lwgeom_free(geom);
3166 
3167  rtn = rt_raster_new(0, 0);
3168  if (rtn == NULL) {
3169  elog(ERROR, "RASTER_clip: Could not create empty raster");
3170  PG_RETURN_NULL();
3171  }
3172 
3173  pgrtn = rt_raster_serialize(rtn);
3174  rt_raster_destroy(rtn);
3175  if (NULL == pgrtn)
3176  PG_RETURN_NULL();
3177 
3178  SET_VARSIZE(pgrtn, pgrtn->size);
3179  PG_RETURN_POINTER(pgrtn);
3180  }
3181 
3182  /* nband (1) */
3183  if (!PG_ARGISNULL(1)) {
3184  array = PG_GETARG_ARRAYTYPE_P(1);
3185  etype = ARR_ELEMTYPE(array);
3186  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3187 
3188  switch (etype) {
3189  case INT2OID:
3190  case INT4OID:
3191  break;
3192  default:
3193  rtpg_clip_arg_destroy(arg);
3194  PG_FREE_IF_COPY(pgraster, 0);
3195  lwgeom_free(geom);
3196  elog(ERROR, "RASTER_clip: Invalid data type for band indexes");
3197  PG_RETURN_NULL();
3198  break;
3199  }
3200 
3201  deconstruct_array(
3202  array, etype,
3203  typlen, typbyval, typalign,
3204  &e, &nulls, &(arg->numbands)
3205  );
3206 
3207  arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3208  if (arg->band == NULL) {
3209  rtpg_clip_arg_destroy(arg);
3210  PG_FREE_IF_COPY(pgraster, 0);
3211  lwgeom_free(geom);
3212  elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3213  PG_RETURN_NULL();
3214  }
3215 
3216  for (i = 0, j = 0; i < arg->numbands; i++) {
3217  if (nulls[i]) continue;
3218 
3219  switch (etype) {
3220  case INT2OID:
3221  arg->band[j].nband = DatumGetInt16(e[i]) - 1;
3222  break;
3223  case INT4OID:
3224  arg->band[j].nband = DatumGetInt32(e[i]) - 1;
3225  break;
3226  }
3227 
3228  j++;
3229  }
3230 
3231  if (j < arg->numbands) {
3232  arg->band = repalloc(arg->band, sizeof(struct rtpg_clip_band_t) * j);
3233  if (arg->band == NULL) {
3234  rtpg_clip_arg_destroy(arg);
3235  PG_FREE_IF_COPY(pgraster, 0);
3236  lwgeom_free(geom);
3237  elog(ERROR, "RASTER_clip: Could not reallocate memory for band arguments");
3238  PG_RETURN_NULL();
3239  }
3240 
3241  arg->numbands = j;
3242  }
3243 
3244  /* validate band */
3245  for (i = 0; i < arg->numbands; i++) {
3246  if (!rt_raster_has_band(arg->raster, arg->band[i].nband)) {
3247  elog(NOTICE, "Band at index %d not found in raster", arg->band[i].nband + 1);
3248  rtpg_clip_arg_destroy(arg);
3249  PG_FREE_IF_COPY(pgraster, 0);
3250  lwgeom_free(geom);
3251  PG_RETURN_NULL();
3252  }
3253 
3254  arg->band[i].hasnodata = 0;
3255  arg->band[i].nodataval = 0;
3256  }
3257  }
3258  else {
3260 
3261  /* raster may have no bands */
3262  if (arg->numbands) {
3263  arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3264  if (arg->band == NULL) {
3265 
3266  rtpg_clip_arg_destroy(arg);
3267  PG_FREE_IF_COPY(pgraster, 0);
3268  lwgeom_free(geom);
3269 
3270  elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3271  PG_RETURN_NULL();
3272  }
3273 
3274  for (i = 0; i < arg->numbands; i++) {
3275  arg->band[i].nband = i;
3276  arg->band[i].hasnodata = 0;
3277  arg->band[i].nodataval = 0;
3278  }
3279  }
3280  }
3281 
3282  /* nodataval (3) */
3283  if (!PG_ARGISNULL(3)) {
3284  array = PG_GETARG_ARRAYTYPE_P(3);
3285  etype = ARR_ELEMTYPE(array);
3286  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3287 
3288  switch (etype) {
3289  case FLOAT4OID:
3290  case FLOAT8OID:
3291  break;
3292  default:
3293  rtpg_clip_arg_destroy(arg);
3294  PG_FREE_IF_COPY(pgraster, 0);
3295  lwgeom_free(geom);
3296  elog(ERROR, "RASTER_clip: Invalid data type for NODATA values");
3297  PG_RETURN_NULL();
3298  break;
3299  }
3300 
3301  deconstruct_array(
3302  array, etype,
3303  typlen, typbyval, typalign,
3304  &e, &nulls, &k
3305  );
3306 
3307  /* it doesn't matter if there are more nodataval */
3308  for (i = 0, j = 0; i < arg->numbands; i++, j++) {
3309  /* cap j to the last nodataval element */
3310  if (j >= k)
3311  j = k - 1;
3312 
3313  if (nulls[j])
3314  continue;
3315 
3316  arg->band[i].hasnodata = 1;
3317  switch (etype) {
3318  case FLOAT4OID:
3319  arg->band[i].nodataval = DatumGetFloat4(e[j]);
3320  break;
3321  case FLOAT8OID:
3322  arg->band[i].nodataval = DatumGetFloat8(e[j]);
3323  break;
3324  }
3325  }
3326  }
3327 
3328  /* get wkb of geometry */
3329  POSTGIS_RT_DEBUG(3, "getting wkb of geometry");
3330  wkb = lwgeom_to_wkb_varlena(geom, WKB_SFSQL);
3331  lwgeom_free(geom);
3332 
3333  /* rasterize geometry */
3334  arg->mask = rt_raster_gdal_rasterize((unsigned char *)wkb->data,
3335  LWSIZE_GET(wkb->size) - LWVARHDRSZ,
3336  NULL,
3337  0,
3338  NULL,
3339  NULL,
3340  NULL,
3341  NULL,
3342  NULL,
3343  NULL,
3344  NULL,
3345  &(gt[1]),
3346  &(gt[5]),
3347  NULL,
3348  NULL,
3349  &(gt[0]),
3350  &(gt[3]),
3351  &(gt[2]),
3352  &(gt[4]),
3353  NULL);
3354 
3355  pfree(wkb);
3356  if (arg->mask == NULL) {
3357  rtpg_clip_arg_destroy(arg);
3358  PG_FREE_IF_COPY(pgraster, 0);
3359  elog(ERROR, "RASTER_clip: Could not rasterize intersection geometry");
3360  PG_RETURN_NULL();
3361  }
3362 
3363  /* set SRID */
3364  rt_raster_set_srid(arg->mask, srid);
3365 
3366  mask_width = rt_raster_get_width(arg->mask);
3367  mask_height = rt_raster_get_height(arg->mask);
3368 
3369  if (rt_raster_from_two_rasters(arg->raster, arg->mask, arg->extenttype, &rtn, offset) != ES_NONE) {
3370  rtpg_clip_arg_destroy(arg);
3371  PG_FREE_IF_COPY(pgraster, 0);
3372  elog(ERROR, "RASTER_clip: Could not compute extent of rasters");
3373  PG_RETURN_NULL();
3374  }
3375 
3376  width = rt_raster_get_width(rtn);
3377  height = rt_raster_get_height(rtn);
3378 
3379  mask_band = rt_raster_get_band(arg->mask, 0);
3380 
3381  for (i = 0; i < arg->numbands; i++) {
3382  input_band = rt_raster_get_band(arg->raster, arg->band[i].nband);
3383 
3384  /* band metadata */
3385  pixtype = rt_band_get_pixtype(input_band);
3386 
3387  if (arg->band[i].hasnodata) {
3388  hasnodata = 1;
3389  nodataval = arg->band[i].nodataval;
3390  }
3391  else if (rt_band_get_hasnodata_flag(input_band)) {
3392  hasnodata = 1;
3393  rt_band_get_nodata(input_band, &nodataval);
3394  }
3395  else {
3396  hasnodata = 0;
3397  nodataval = rt_band_get_min_value(input_band);
3398  }
3399 
3400  if (rt_raster_generate_new_band(rtn, pixtype, nodataval, hasnodata, nodataval, i) < 0) {
3401  rtpg_clip_arg_destroy(arg);
3402  PG_FREE_IF_COPY(pgraster, 0);
3403  elog(ERROR, "RASTER_clip: Could not add new band in output raster");
3404  PG_RETURN_NULL();
3405  }
3406 
3407  if (rt_band_get_isnodata_flag(input_band)) {
3408  continue;
3409  }
3410 
3411  output_band = rt_raster_get_band(rtn, arg->band[i].nband);
3412 
3413  if (!mask_band) {
3414  continue;
3415  }
3416 
3417  for (y = 0; y < height; y++) {
3418  for (x = 0; x < width; x++) {
3419  mask_x = x - (int)offset[2];
3420  mask_y = y - (int)offset[3];
3421 
3422  if (!(
3423  mask_x >= 0 &&
3424  mask_x < mask_width &&
3425  mask_y >= 0 &&
3426  mask_y < mask_height
3427  )) {
3428  continue;
3429  }
3430 
3431  if (rt_band_get_pixel(mask_band, mask_x, mask_y, &value, &isnodata) != ES_NONE) {
3432  rtpg_clip_arg_destroy(arg);
3433  PG_FREE_IF_COPY(pgraster, 0);
3434  elog(ERROR, "RASTER_clip: Could not get pixel value");
3435  PG_RETURN_NULL();
3436  }
3437 
3438  if (isnodata) {
3439  continue;
3440  }
3441 
3442  input_x = x - (int)offset[0];
3443  input_y = y - (int)offset[1];
3444 
3445  if (rt_band_get_pixel(input_band, input_x, input_y, &value, &isnodata) != ES_NONE) {
3446  rtpg_clip_arg_destroy(arg);
3447  PG_FREE_IF_COPY(pgraster, 0);
3448  elog(ERROR, "RASTER_clip: Could not get pixel value");
3449  PG_RETURN_NULL();
3450  }
3451 
3452  if (isnodata) {
3453  continue;
3454  }
3455 
3456  if (rt_band_set_pixel(output_band, x, y, value, NULL)) {
3457  rtpg_clip_arg_destroy(arg);
3458  PG_FREE_IF_COPY(pgraster, 0);
3459  elog(ERROR, "RASTER_clip: Could not set pixel value");
3460  PG_RETURN_NULL();
3461  }
3462  }
3463  }
3464  }
3465 
3466  rtpg_clip_arg_destroy(arg);
3467  PG_FREE_IF_COPY(pgraster, 0);
3468 
3469  pgrtn = rt_raster_serialize(rtn);
3470  rt_raster_destroy(rtn);
3471 
3472  POSTGIS_RT_DEBUG(3, "Finished");
3473 
3474  if (!pgrtn)
3475  PG_RETURN_NULL();
3476 
3477  SET_VARSIZE(pgrtn, pgrtn->size);
3478  PG_RETURN_POINTER(pgrtn);
3479 }
3480 
3485 Datum RASTER_reclass(PG_FUNCTION_ARGS) {
3486  rt_pgraster *pgraster = NULL;
3487  rt_pgraster *pgrtn = NULL;
3488  rt_raster raster = NULL;
3489  rt_band band = NULL;
3490  rt_band newband = NULL;
3491  uint32_t numBands = 0;
3492 
3493  ArrayType *array;
3494  Oid etype;
3495  Datum *e;
3496  bool *nulls;
3497  int16 typlen;
3498  bool typbyval;
3499  char typalign;
3500  int n = 0;
3501 
3502  int i = 0;
3503  int j = 0;
3504  int k = 0;
3505 
3506  uint32_t a = 0;
3507  uint32_t b = 0;
3508  uint32_t c = 0;
3509 
3510  rt_reclassexpr *exprset = NULL;
3511  HeapTupleHeader tup;
3512  bool isnull;
3513  Datum tupv;
3514  uint32_t nband = 0;
3515  char *expr = NULL;
3516  text *exprtext = NULL;
3517  double val = 0;
3518  char *junk = NULL;
3519  int inc_val = 0;
3520  int exc_val = 0;
3521  char *pixeltype = NULL;
3522  text *pixeltypetext = NULL;
3523  rt_pixtype pixtype = PT_END;
3524  double nodataval = 0;
3525  bool hasnodata = FALSE;
3526 
3527  char **comma_set = NULL;
3528  uint32_t comma_n = 0;
3529  char **colon_set = NULL;
3530  uint32_t colon_n = 0;
3531  char **dash_set = NULL;
3532  uint32_t dash_n = 0;
3533 
3534  POSTGIS_RT_DEBUG(3, "RASTER_reclass: Starting");
3535 
3536  /* pgraster is null, return null */
3537  if (PG_ARGISNULL(0))
3538  PG_RETURN_NULL();
3539  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3540 
3541  /* raster */
3542  raster = rt_raster_deserialize(pgraster, FALSE);
3543  if (!raster) {
3544  PG_FREE_IF_COPY(pgraster, 0);
3545  elog(ERROR, "RASTER_reclass: Could not deserialize raster");
3546  PG_RETURN_NULL();
3547  }
3548  numBands = rt_raster_get_num_bands(raster);
3549  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: %d possible bands to be reclassified", numBands);
3550 
3551  /* process set of reclassarg */
3552  POSTGIS_RT_DEBUG(3, "RASTER_reclass: Processing Arg 1 (reclassargset)");
3553  array = PG_GETARG_ARRAYTYPE_P(1);
3554  etype = ARR_ELEMTYPE(array);
3555  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3556 
3557  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
3558  &nulls, &n);
3559 
3560  if (!n) {
3561  elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
3562 
3563  pgrtn = rt_raster_serialize(raster);
3565  PG_FREE_IF_COPY(pgraster, 0);
3566  if (!pgrtn)
3567  PG_RETURN_NULL();
3568 
3569  SET_VARSIZE(pgrtn, pgrtn->size);
3570  PG_RETURN_POINTER(pgrtn);
3571  }
3572 
3573  /*
3574  process each element of reclassarg
3575  each element is the index of the band to process, the set
3576  of reclass ranges and the output band's pixeltype
3577  */
3578  for (i = 0; i < n; i++) {
3579  if (nulls[i]) continue;
3580 
3581  /* each element is a tuple */
3582  tup = (HeapTupleHeader) DatumGetPointer(e[i]);
3583  if (NULL == tup) {
3584  elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
3585 
3586  pgrtn = rt_raster_serialize(raster);
3588  PG_FREE_IF_COPY(pgraster, 0);
3589  if (!pgrtn)
3590  PG_RETURN_NULL();
3591 
3592  SET_VARSIZE(pgrtn, pgrtn->size);
3593  PG_RETURN_POINTER(pgrtn);
3594  }
3595 
3596  /* band index (1-based) */
3597  tupv = GetAttributeByName(tup, "nband", &isnull);
3598  if (isnull) {
3599  elog(NOTICE, "Invalid argument for reclassargset. Missing value of nband for reclassarg of index %d . Returning original raster", i);
3600 
3601  pgrtn = rt_raster_serialize(raster);
3603  PG_FREE_IF_COPY(pgraster, 0);
3604  if (!pgrtn)
3605  PG_RETURN_NULL();
3606 
3607  SET_VARSIZE(pgrtn, pgrtn->size);
3608  PG_RETURN_POINTER(pgrtn);
3609  }
3610  nband = DatumGetInt32(tupv);
3611  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: expression for band %d", nband);
3612 
3613  /* valid band index? */
3614  if (nband < 1 || nband > numBands) {
3615  elog(NOTICE, "Invalid argument for reclassargset. Invalid band index (must use 1-based) 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 
3627  /* reclass expr */
3628  tupv = GetAttributeByName(tup, "reclassexpr", &isnull);
3629  if (isnull) {
3630  elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3631 
3632  pgrtn = rt_raster_serialize(raster);
3634  PG_FREE_IF_COPY(pgraster, 0);
3635  if (!pgrtn)
3636  PG_RETURN_NULL();
3637 
3638  SET_VARSIZE(pgrtn, pgrtn->size);
3639  PG_RETURN_POINTER(pgrtn);
3640  }
3641  exprtext = (text *) DatumGetPointer(tupv);
3642  if (NULL == exprtext) {
3643  elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3644 
3645  pgrtn = rt_raster_serialize(raster);
3647  PG_FREE_IF_COPY(pgraster, 0);
3648  if (!pgrtn)
3649  PG_RETURN_NULL();
3650 
3651  SET_VARSIZE(pgrtn, pgrtn->size);
3652  PG_RETURN_POINTER(pgrtn);
3653  }
3654  expr = text_to_cstring(exprtext);
3655  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: expr (raw) %s", expr);
3656  expr = rtpg_removespaces(expr);
3657  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: expr (clean) %s", expr);
3658 
3659  /* split string to its components */
3660  /* comma (,) separating rangesets */
3661  comma_set = rtpg_strsplit(expr, ",", &comma_n);
3662  if (comma_n < 1) {
3663  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3664 
3665  pgrtn = rt_raster_serialize(raster);
3667  PG_FREE_IF_COPY(pgraster, 0);
3668  if (!pgrtn)
3669  PG_RETURN_NULL();
3670 
3671  SET_VARSIZE(pgrtn, pgrtn->size);
3672  PG_RETURN_POINTER(pgrtn);
3673  }
3674 
3675  /* set of reclass expressions */
3676  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: %d possible expressions", comma_n);
3677  exprset = palloc(comma_n * sizeof(rt_reclassexpr));
3678 
3679  for (a = 0, j = 0; a < comma_n; a++) {
3680  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: map %s", comma_set[a]);
3681 
3682  /* colon (:) separating range "src" and "dst" */
3683  colon_set = rtpg_strsplit(comma_set[a], ":", &colon_n);
3684  if (colon_n != 2) {
3685  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3686  for (k = 0; k < j; k++) pfree(exprset[k]);
3687  pfree(exprset);
3688 
3689  pgrtn = rt_raster_serialize(raster);
3691  PG_FREE_IF_COPY(pgraster, 0);
3692  if (!pgrtn)
3693  PG_RETURN_NULL();
3694 
3695  SET_VARSIZE(pgrtn, pgrtn->size);
3696  PG_RETURN_POINTER(pgrtn);
3697  }
3698 
3699  /* allocate mem for reclass expression */
3700  exprset[j] = palloc(sizeof(struct rt_reclassexpr_t));
3701 
3702  for (b = 0; b < colon_n; b++) {
3703  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: range %s", colon_set[b]);
3704 
3705  /* dash (-) separating "min" and "max" */
3706  dash_set = rtpg_strsplit(colon_set[b], "-", &dash_n);
3707  if (dash_n < 1 || dash_n > 3) {
3708  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3709  for (k = 0; k < j; k++) pfree(exprset[k]);
3710  pfree(exprset);
3711 
3712  pgrtn = rt_raster_serialize(raster);
3714  PG_FREE_IF_COPY(pgraster, 0);
3715  if (!pgrtn)
3716  PG_RETURN_NULL();
3717 
3718  SET_VARSIZE(pgrtn, pgrtn->size);
3719  PG_RETURN_POINTER(pgrtn);
3720  }
3721 
3722  for (c = 0; c < dash_n; c++) {
3723  /* need to handle: (-9999-100 -> "(", "9999", "100" */
3724  if (
3725  c < 1 &&
3726  strlen(dash_set[c]) == 1 && (
3727  strchr(dash_set[c], '(') != NULL ||
3728  strchr(dash_set[c], '[') != NULL ||
3729  strchr(dash_set[c], ')') != NULL ||
3730  strchr(dash_set[c], ']') != NULL
3731  )
3732  ) {
3733  uint32_t dash_it;
3734  junk = palloc(sizeof(char) * (strlen(dash_set[c + 1]) + 2));
3735  if (NULL == junk) {
3736  for (k = 0; k <= j; k++) pfree(exprset[k]);
3737  pfree(exprset);
3739  PG_FREE_IF_COPY(pgraster, 0);
3740 
3741  elog(ERROR, "RASTER_reclass: Could not allocate memory");
3742  PG_RETURN_NULL();
3743  }
3744 
3745  sprintf(junk, "%s%s", dash_set[c], dash_set[c + 1]);
3746  c++;
3747  dash_set[c] = repalloc(dash_set[c], sizeof(char) * (strlen(junk) + 1));
3748  strcpy(dash_set[c], junk);
3749  pfree(junk);
3750 
3751  /* rebuild dash_set */
3752  for (dash_it = 1; dash_it < dash_n; dash_it++) {
3753  dash_set[dash_it - 1] = repalloc(dash_set[dash_it - 1], (strlen(dash_set[dash_it]) + 1) * sizeof(char));
3754  strcpy(dash_set[dash_it - 1], dash_set[dash_it]);
3755  }
3756  dash_n--;
3757  c--;
3758  pfree(dash_set[dash_n]);
3759  dash_set = repalloc(dash_set, sizeof(char *) * dash_n);
3760  }
3761 
3762  /* there shouldn't be more than two in dash_n */
3763  if (c < 1 && dash_n > 2) {
3764  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3765  for (k = 0; k < j; k++) pfree(exprset[k]);
3766  pfree(exprset);
3767 
3768  pgrtn = rt_raster_serialize(raster);
3770  PG_FREE_IF_COPY(pgraster, 0);
3771  if (!pgrtn)
3772  PG_RETURN_NULL();
3773 
3774  SET_VARSIZE(pgrtn, pgrtn->size);
3775  PG_RETURN_POINTER(pgrtn);
3776  }
3777 
3778  /* check interval flags */
3779  exc_val = 0;
3780  inc_val = 1;
3781  /* range */
3782  if (dash_n != 1) {
3783  /* min */
3784  if (c < 1) {
3785  if (
3786  strchr(dash_set[c], ')') != NULL ||
3787  strchr(dash_set[c], ']') != NULL
3788  ) {
3789  exc_val = 1;
3790  inc_val = 1;
3791  }
3792  else if (strchr(dash_set[c], '(') != NULL){
3793  inc_val = 0;
3794  }
3795  else {
3796  inc_val = 1;
3797  }
3798  }
3799  /* max */
3800  else {
3801  if (
3802  strrchr(dash_set[c], '(') != NULL ||
3803  strrchr(dash_set[c], '[') != NULL
3804  ) {
3805  exc_val = 1;
3806  inc_val = 0;
3807  }
3808  else if (strrchr(dash_set[c], ']') != NULL) {
3809  inc_val = 1;
3810  }
3811  else {
3812  inc_val = 0;
3813  }
3814  }
3815  }
3816  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: exc_val %d inc_val %d", exc_val, inc_val);
3817 
3818  /* remove interval flags */
3819  dash_set[c] = rtpg_chartrim(dash_set[c], "()[]");
3820  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (char) %s", dash_set[c]);
3821 
3822  /* value from string to double */
3823  errno = 0;
3824  val = strtod(dash_set[c], &junk);
3825  if (errno != 0 || dash_set[c] == junk) {
3826  elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3827  for (k = 0; k < j; k++) pfree(exprset[k]);
3828  pfree(exprset);
3829 
3830  pgrtn = rt_raster_serialize(raster);
3832  PG_FREE_IF_COPY(pgraster, 0);
3833  if (!pgrtn)
3834  PG_RETURN_NULL();
3835 
3836  SET_VARSIZE(pgrtn, pgrtn->size);
3837  PG_RETURN_POINTER(pgrtn);
3838  }
3839  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (double) %f", val);
3840 
3841  /* strsplit removes dash (a.k.a. negative sign), compare now to restore */
3842  if (c < 1)
3843  junk = strstr(colon_set[b], dash_set[c]);
3844  else
3845  junk = rtpg_strrstr(colon_set[b], dash_set[c]);
3847  4,
3848  "(colon_set[%d], dash_set[%d], junk) = (%s, %s, %s)",
3849  b, c, colon_set[b], dash_set[c], junk
3850  );
3851  /* not beginning of string */
3852  if (junk != colon_set[b]) {
3853  /* prior is a dash */
3854  if (*(junk - 1) == '-') {
3855  /* prior is beginning of string or prior - 1 char is dash, negative number */
3856  if (
3857  ((junk - 1) == colon_set[b]) ||
3858  (*(junk - 2) == '-') ||
3859  (*(junk - 2) == '[') ||
3860  (*(junk - 2) == '(')
3861  ) {
3862  val *= -1.;
3863  }
3864  }
3865  }
3866  POSTGIS_RT_DEBUGF(4, "RASTER_reclass: min/max (double) %f", val);
3867 
3868  /* src */
3869  if (b < 1) {
3870  /* singular value */
3871  if (dash_n == 1) {
3872  exprset[j]->src.exc_min = exprset[j]->src.exc_max = exc_val;
3873  exprset[j]->src.inc_min = exprset[j]->src.inc_max = inc_val;
3874  exprset[j]->src.min = exprset[j]->src.max = val;
3875  }
3876  /* min */
3877  else if (c < 1) {
3878  exprset[j]->src.exc_min = exc_val;
3879  exprset[j]->src.inc_min = inc_val;
3880  exprset[j]->src.min = val;
3881  }
3882  /* max */
3883  else {
3884  exprset[j]->src.exc_max = exc_val;
3885  exprset[j]->src.inc_max = inc_val;
3886  exprset[j]->src.max = val;
3887  }
3888  }
3889  /* dst */
3890  else {
3891  /* singular value */
3892  if (dash_n == 1)
3893  exprset[j]->dst.min = exprset[j]->dst.max = val;
3894  /* min */
3895  else if (c < 1)
3896  exprset[j]->dst.min = val;
3897  /* max */
3898  else
3899  exprset[j]->dst.max = val;
3900  }
3901  }
3902  pfree(dash_set);
3903  }
3904  pfree(colon_set);
3905 
3906  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: or: %f - %f nr: %f - %f"
3907  , exprset[j]->src.min
3908  , exprset[j]->src.max
3909  , exprset[j]->dst.min
3910  , exprset[j]->dst.max
3911  );
3912  j++;
3913  }
3914  pfree(comma_set);
3915 
3916  /* pixel type */
3917  tupv = GetAttributeByName(tup, "pixeltype", &isnull);
3918  if (isnull) {
3919  elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3920 
3921  pgrtn = rt_raster_serialize(raster);
3923  PG_FREE_IF_COPY(pgraster, 0);
3924  if (!pgrtn)
3925  PG_RETURN_NULL();
3926 
3927  SET_VARSIZE(pgrtn, pgrtn->size);
3928  PG_RETURN_POINTER(pgrtn);
3929  }
3930  pixeltypetext = (text *) DatumGetPointer(tupv);
3931  if (NULL == pixeltypetext) {
3932  elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3933 
3934  pgrtn = rt_raster_serialize(raster);
3936  PG_FREE_IF_COPY(pgraster, 0);
3937  if (!pgrtn)
3938  PG_RETURN_NULL();
3939 
3940  SET_VARSIZE(pgrtn, pgrtn->size);
3941  PG_RETURN_POINTER(pgrtn);
3942  }
3943  pixeltype = text_to_cstring(pixeltypetext);
3944  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: pixeltype %s", pixeltype);
3945  pixtype = rt_pixtype_index_from_name(pixeltype);
3946 
3947  /* nodata */
3948  tupv = GetAttributeByName(tup, "nodataval", &isnull);
3949  if (isnull) {
3950  nodataval = 0;
3951  hasnodata = FALSE;
3952  }
3953  else {
3954  nodataval = DatumGetFloat8(tupv);
3955  hasnodata = TRUE;
3956  }
3957  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: nodataval %f", nodataval);
3958  POSTGIS_RT_DEBUGF(3, "RASTER_reclass: hasnodata %d", hasnodata);
3959 
3960  /* do reclass */
3962  if (!band) {
3963  elog(NOTICE, "Could not find raster band of index %d. Returning original raster", nband);
3964  for (k = 0; k < j; k++) pfree(exprset[k]);
3965  pfree(exprset);
3966 
3967  pgrtn = rt_raster_serialize(raster);
3969  PG_FREE_IF_COPY(pgraster, 0);
3970  if (!pgrtn)
3971  PG_RETURN_NULL();
3972 
3973  SET_VARSIZE(pgrtn, pgrtn->size);
3974  PG_RETURN_POINTER(pgrtn);
3975  }
3976  newband = rt_band_reclass(band, pixtype, hasnodata, nodataval, exprset, j);
3977  if (!newband) {
3978  for (k = 0; k < j; k++) pfree(exprset[k]);
3979  pfree(exprset);
3980 
3982  PG_FREE_IF_COPY(pgraster, 0);
3983 
3984  elog(ERROR, "RASTER_reclass: Could not reclassify raster band of index %d", nband);
3985  PG_RETURN_NULL();
3986  }
3987 
3988  /* replace old band with new band */
3989  if (rt_raster_replace_band(raster, newband, nband - 1) == NULL) {
3990  for (k = 0; k < j; k++) pfree(exprset[k]);
3991  pfree(exprset);
3992 
3993  rt_band_destroy(newband);
3995  PG_FREE_IF_COPY(pgraster, 0);
3996 
3997  elog(ERROR, "RASTER_reclass: Could not replace raster band of index %d with reclassified band", nband);
3998  PG_RETURN_NULL();
3999  }
4000 
4001  /* old band is in the variable band */
4003 
4004  /* free exprset */
4005  for (k = 0; k < j; k++) pfree(exprset[k]);
4006  pfree(exprset);
4007  }
4008 
4009  pgrtn = rt_raster_serialize(raster);
4011  PG_FREE_IF_COPY(pgraster, 0);
4012  if (!pgrtn)
4013  PG_RETURN_NULL();
4014 
4015  POSTGIS_RT_DEBUG(3, "RASTER_reclass: Finished");
4016 
4017  SET_VARSIZE(pgrtn, pgrtn->size);
4018  PG_RETURN_POINTER(pgrtn);
4019 }
4020 
4021 /* ---------------------------------------------------------------- */
4022 /* apply colormap to specified band of a raster */
4023 /* ---------------------------------------------------------------- */
4024 
4028  int nband; /* 1-based */
4031 
4034 
4035  char **entry;
4036  uint32_t nentry;
4037  char **element;
4038  uint32_t nelement;
4039 };
4040 
4041 static rtpg_colormap_arg
4043  rtpg_colormap_arg arg = NULL;
4044 
4045  arg = palloc(sizeof(struct rtpg_colormap_arg_t));
4046  if (arg == NULL) {
4047  elog(ERROR, "rtpg_colormap_arg: Could not allocate memory for function arguments");
4048  return NULL;
4049  }
4050 
4051  arg->raster = NULL;
4052  arg->nband = 1;
4053  arg->band = NULL;
4054  arg->bandstats = NULL;
4055 
4056  arg->colormap = palloc(sizeof(struct rt_colormap_t));
4057  if (arg->colormap == NULL) {
4058  elog(ERROR, "rtpg_colormap_arg: Could not allocate memory for function arguments");
4059  return NULL;
4060  }
4061  arg->colormap->nentry = 0;
4062  arg->colormap->entry = NULL;
4063  arg->colormap->ncolor = 4; /* assume RGBA */
4064  arg->colormap->method = CM_INTERPOLATE;
4065  arg->nodataentry = -1;
4066 
4067  arg->entry = NULL;
4068  arg->nentry = 0;
4069  arg->element = NULL;
4070  arg->nelement = 0;
4071 
4072  return arg;
4073 }
4074 
4075 static void
4077  uint32_t i = 0;
4078  if (arg->raster != NULL)
4079  rt_raster_destroy(arg->raster);
4080 
4081  if (arg->bandstats != NULL)
4082  pfree(arg->bandstats);
4083 
4084  if (arg->colormap != NULL) {
4085  if (arg->colormap->entry != NULL)
4086  pfree(arg->colormap->entry);
4087  pfree(arg->colormap);
4088  }
4089 
4090  if (arg->nentry) {
4091  for (i = 0; i < arg->nentry; i++) {
4092  if (arg->entry[i] != NULL)
4093  pfree(arg->entry[i]);
4094  }
4095  pfree(arg->entry);
4096  }
4097 
4098  if (arg->nelement) {
4099  for (i = 0; i < arg->nelement; i++)
4100  pfree(arg->element[i]);
4101  pfree(arg->element);
4102  }
4103 
4104  pfree(arg);
4105  arg = NULL;
4106 }
4107 
4109 Datum RASTER_colorMap(PG_FUNCTION_ARGS)
4110 {
4111  rt_pgraster *pgraster = NULL;
4112  rtpg_colormap_arg arg = NULL;
4113  char *junk = NULL;
4114  rt_raster raster = NULL;
4115 
4116  POSTGIS_RT_DEBUG(3, "RASTER_colorMap: Starting");
4117 
4118  /* pgraster is NULL, return NULL */
4119  if (PG_ARGISNULL(0))
4120  PG_RETURN_NULL();
4121 
4122  /* init arg */
4123  arg = rtpg_colormap_arg_init();
4124  if (arg == NULL) {
4125  elog(ERROR, "RASTER_colorMap: Could not initialize argument structure");
4126  PG_RETURN_NULL();
4127  }
4128 
4129  /* raster (0) */
4130  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4131 
4132  /* raster object */
4133  arg->raster = rt_raster_deserialize(pgraster, FALSE);
4134  if (!arg->raster) {
4136  PG_FREE_IF_COPY(pgraster, 0);
4137  elog(ERROR, "RASTER_colorMap: Could not deserialize raster");
4138  PG_RETURN_NULL();
4139  }
4140 
4141  /* nband (1) */
4142  if (!PG_ARGISNULL(1))
4143  arg->nband = PG_GETARG_INT32(1);
4144  POSTGIS_RT_DEBUGF(4, "nband = %d", arg->nband);
4145 
4146  /* check that band works */
4147  if (!rt_raster_has_band(arg->raster, arg->nband - 1)) {
4148  elog(NOTICE, "Raster does not have band at index %d. Returning empty raster", arg->nband);
4149 
4150  raster = rt_raster_clone(arg->raster, 0);
4151  if (raster == NULL) {
4153  PG_FREE_IF_COPY(pgraster, 0);
4154  elog(ERROR, "RASTER_colorMap: Could not create empty raster");
4155  PG_RETURN_NULL();
4156  }
4157 
4159  PG_FREE_IF_COPY(pgraster, 0);
4160 
4161  pgraster = rt_raster_serialize(raster);
4163  if (pgraster == NULL)
4164  PG_RETURN_NULL();
4165 
4166  SET_VARSIZE(pgraster, ((rt_pgraster*) pgraster)->size);
4167  PG_RETURN_POINTER(pgraster);
4168  }
4169 
4170  /* get band */
4171  arg->band = rt_raster_get_band(arg->raster, arg->nband - 1);
4172  if (arg->band == NULL) {
4173  int nband = arg->nband;
4175  PG_FREE_IF_COPY(pgraster, 0);
4176  elog(ERROR, "RASTER_colorMap: Could not get band at index %d", nband);
4177  PG_RETURN_NULL();
4178  }
4179 
4180  /* method (3) */
4181  if (!PG_ARGISNULL(3)) {
4182  char *method = NULL;
4183  char *tmp = text_to_cstring(PG_GETARG_TEXT_P(3));
4184  POSTGIS_RT_DEBUGF(4, "raw method = %s", tmp);
4185 
4186  method = rtpg_trim(tmp);
4187  pfree(tmp);
4188  method = rtpg_strtoupper(method);
4189 
4190  if (strcmp(method, "INTERPOLATE") == 0)
4191  arg->colormap->method = CM_INTERPOLATE;
4192  else if (strcmp(method, "EXACT") == 0)
4193  arg->colormap->method = CM_EXACT;
4194  else if (strcmp(method, "NEAREST") == 0)
4195  arg->colormap->method = CM_NEAREST;
4196  else {
4197  elog(NOTICE, "Unknown value provided for method. Defaulting to INTERPOLATE");
4198  arg->colormap->method = CM_INTERPOLATE;
4199  }
4200  }
4201  /* default to INTERPOLATE */
4202  else
4203  arg->colormap->method = CM_INTERPOLATE;
4204  POSTGIS_RT_DEBUGF(4, "method = %d", arg->colormap->method);
4205 
4206  /* colormap (2) */
4207  if (PG_ARGISNULL(2)) {
4209  PG_FREE_IF_COPY(pgraster, 0);
4210  elog(ERROR, "RASTER_colorMap: Value must be provided for colormap");
4211  PG_RETURN_NULL();
4212  }
4213  else {
4214  char *tmp = NULL;
4215  char *colormap = text_to_cstring(PG_GETARG_TEXT_P(2));
4216  char *_entry;
4217  char *_element;
4218  uint32_t i = 0;
4219  uint32_t j = 0;
4220 
4221  POSTGIS_RT_DEBUGF(4, "colormap = %s", colormap);
4222 
4223  /* empty string */
4224  if (!strlen(colormap)) {
4226  PG_FREE_IF_COPY(pgraster, 0);
4227  elog(ERROR, "RASTER_colorMap: Value must be provided for colormap");
4228  PG_RETURN_NULL();
4229  }
4230 
4231  arg->entry = rtpg_strsplit(colormap, "\n", &(arg->nentry));
4232  pfree(colormap);
4233  if (arg->nentry < 1) {
4235  PG_FREE_IF_COPY(pgraster, 0);
4236  elog(ERROR, "RASTER_colorMap: Could not process the value provided for colormap");
4237  PG_RETURN_NULL();
4238  }
4239 
4240  /* allocate the max # of colormap entries */
4241  arg->colormap->entry = palloc(sizeof(struct rt_colormap_entry_t) * arg->nentry);
4242  if (arg->colormap->entry == NULL) {
4244  PG_FREE_IF_COPY(pgraster, 0);
4245  elog(ERROR, "RASTER_colorMap: Could not allocate memory for colormap entries");
4246  PG_RETURN_NULL();
4247  }
4248  memset(arg->colormap->entry, 0, sizeof(struct rt_colormap_entry_t) * arg->nentry);
4249 
4250  /* each entry */
4251  for (i = 0; i < arg->nentry; i++) {
4252  /* substitute space for other delimiters */
4253  tmp = rtpg_strreplace(arg->entry[i], ":", " ", NULL);
4254  _entry = rtpg_strreplace(tmp, ",", " ", NULL);
4255  pfree(tmp);
4256  tmp = rtpg_strreplace(_entry, "\t", " ", NULL);
4257  pfree(_entry);
4258  _entry = rtpg_trim(tmp);
4259  pfree(tmp);
4260 
4261  POSTGIS_RT_DEBUGF(4, "Processing entry[%d] = %s", i, arg->entry[i]);
4262  POSTGIS_RT_DEBUGF(4, "Cleaned entry[%d] = %s", i, _entry);
4263 
4264  /* empty entry, continue */
4265  if (!strlen(_entry)) {
4266  POSTGIS_RT_DEBUGF(3, "Skipping empty entry[%d]", i);
4267  pfree(_entry);
4268  continue;
4269  }
4270 
4271  arg->element = rtpg_strsplit(_entry, " ", &(arg->nelement));
4272  pfree(_entry);
4273  if (arg->nelement < 2) {
4275  PG_FREE_IF_COPY(pgraster, 0);
4276  elog(ERROR, "RASTER_colorMap: Could not process colormap entry %d", i + 1);
4277  PG_RETURN_NULL();
4278  }
4279  else if (arg->nelement > 5) {
4280  elog(NOTICE, "More than five elements in colormap entry %d. Using at most five elements", i + 1);
4281  arg->nelement = 5;
4282  }
4283 
4284  /* smallest # of colors */
4285  if (((int)arg->nelement - 1) < arg->colormap->ncolor)
4286  arg->colormap->ncolor = arg->nelement - 1;
4287 
4288  /* each element of entry */
4289  for (j = 0; j < arg->nelement; j++) {
4290 
4291  _element = rtpg_trim(arg->element[j]);
4292  _element = rtpg_strtoupper(_element);
4293  POSTGIS_RT_DEBUGF(4, "Processing entry[%d][%d] = %s", i, j, arg->element[j]);
4294  POSTGIS_RT_DEBUGF(4, "Cleaned entry[%d][%d] = %s", i, j, _element);
4295 
4296  /* first element is ALWAYS a band value, percentage OR "nv" string */
4297  if (j == 0) {
4298  char *percent = NULL;
4299 
4300  /* NODATA */
4301  if (
4302  strcmp(_element, "NV") == 0 ||
4303  strcmp(_element, "NULL") == 0 ||
4304  strcmp(_element, "NODATA") == 0
4305  ) {
4306  POSTGIS_RT_DEBUG(4, "Processing NODATA string");
4307 
4308  if (arg->nodataentry > -1) {
4309  elog(NOTICE, "More than one NODATA entry found. Using only the first one");
4310  }
4311  else {
4312  arg->colormap->entry[arg->colormap->nentry].isnodata = 1;
4313  /* no need to set value as value comes from band's NODATA */
4314  arg->colormap->entry[arg->colormap->nentry].value = 0;
4315  }
4316  }
4317  /* percent value */
4318  else if ((percent = strchr(_element, '%')) != NULL) {
4319  double value;
4320  POSTGIS_RT_DEBUG(4, "Processing percent string");
4321 
4322  /* get the band stats */
4323  if (arg->bandstats == NULL) {
4324  POSTGIS_RT_DEBUG(4, "Getting band stats");
4325 
4326  arg->bandstats = rt_band_get_summary_stats(arg->band, 1, 1, 0, NULL, NULL, NULL);
4327  if (arg->bandstats == NULL) {
4328  pfree(_element);
4330  PG_FREE_IF_COPY(pgraster, 0);
4331  elog(ERROR, "RASTER_colorMap: Could not get band's summary stats to process percentages");
4332  PG_RETURN_NULL();
4333  }
4334  }
4335 
4336  /* get the string before the percent char */
4337  tmp = palloc(sizeof(char) * (percent - _element + 1));
4338  if (tmp == NULL) {
4339  pfree(_element);
4341  PG_FREE_IF_COPY(pgraster, 0);
4342  elog(ERROR, "RASTER_colorMap: Could not allocate memory for value of percentage");
4343  PG_RETURN_NULL();
4344  }
4345 
4346  memcpy(tmp, _element, percent - _element);
4347  tmp[percent - _element] = '\0';
4348  POSTGIS_RT_DEBUGF(4, "Percent value = %s", tmp);
4349 
4350  /* get percentage value */
4351  errno = 0;
4352  value = strtod(tmp, NULL);
4353  pfree(tmp);
4354  if (errno != 0 || _element == junk) {
4355  pfree(_element);
4357  PG_FREE_IF_COPY(pgraster, 0);
4358  elog(ERROR, "RASTER_colorMap: Could not process percent string to value");
4359  PG_RETURN_NULL();
4360  }
4361 
4362  /* check percentage */
4363  if (value < 0.) {
4364  elog(NOTICE, "Percentage values cannot be less than zero. Defaulting to zero");
4365  value = 0.;
4366  }
4367  else if (value > 100.) {
4368  elog(NOTICE, "Percentage values cannot be greater than 100. Defaulting to 100");
4369  value = 100.;
4370  }
4371 
4372  /* get the true pixel value */
4373  /* TODO: should the percentage be quantile based? */
4374  arg->colormap->entry[arg->colormap->nentry].value = ((value / 100.) * (arg->bandstats->max - arg->bandstats->min)) + arg->bandstats->min;
4375  }
4376  /* straight value */
4377  else {
4378  errno = 0;
4379  arg->colormap->entry[arg->colormap->nentry].value = strtod(_element, &junk);
4380  if (errno != 0 || _element == junk) {
4381  pfree(_element);
4383  PG_FREE_IF_COPY(pgraster, 0);
4384  elog(ERROR, "RASTER_colorMap: Could not process string to value");
4385  PG_RETURN_NULL();
4386  }
4387  }
4388 
4389  }
4390  /* RGB values (0 - 255) */
4391  else {
4392  int value = 0;
4393 
4394  errno = 0;
4395  value = (int) 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  if (value > 255) {
4405  elog(NOTICE, "RGBA value cannot be greater than 255. Defaulting to 255");
4406  value = 255;
4407  }
4408  else if (value < 0) {
4409  elog(NOTICE, "RGBA value cannot be less than zero. Defaulting to zero");
4410  value = 0;
4411  }
4412  arg->colormap->entry[arg->colormap->nentry].color[j - 1] = value;
4413  }
4414 
4415  pfree(_element);
4416  }
4417 
4418  POSTGIS_RT_DEBUGF(4, "colormap->entry[%d] (isnodata, value, R, G, B, A) = (%d, %f, %d, %d, %d, %d)",
4419  arg->colormap->nentry,
4420  arg->colormap->entry[arg->colormap->nentry].isnodata,
4421  arg->colormap->entry[arg->colormap->nentry].value,
4422  arg->colormap->entry[arg->colormap->nentry].color[0],
4423  arg->colormap->entry[arg->colormap->nentry].color[1],
4424  arg->colormap->entry[arg->colormap->nentry].color[2],
4425  arg->colormap->entry[arg->colormap->nentry].color[3]
4426  );
4427 
4428  arg->colormap->nentry++;
4429  }
4430 
4431  POSTGIS_RT_DEBUGF(4, "colormap->nentry = %d", arg->colormap->nentry);
4432  POSTGIS_RT_DEBUGF(4, "colormap->ncolor = %d", arg->colormap->ncolor);
4433  }
4434 
4435  /* call colormap */
4436  raster = rt_raster_colormap(arg->raster, arg->nband - 1, arg->colormap);
4437  if (raster == NULL) {
4439  PG_FREE_IF_COPY(pgraster, 0);
4440  elog(ERROR, "RASTER_colorMap: Could not create new raster with applied colormap");
4441  PG_RETURN_NULL();
4442  }
4443 
4445  PG_FREE_IF_COPY(pgraster, 0);
4446  pgraster = rt_raster_serialize(raster);
4448 
4449  POSTGIS_RT_DEBUG(3, "RASTER_colorMap: Done");
4450 
4451  if (pgraster == NULL)
4452  PG_RETURN_NULL();
4453 
4454  SET_VARSIZE(pgraster, ((rt_pgraster*) pgraster)->size);
4455  PG_RETURN_POINTER(pgraster);
4456 }
4457 
4459 Datum RASTER_mapAlgebraExpr(PG_FUNCTION_ARGS)
4460 {
4461  rt_pgraster *pgraster = NULL;
4462  rt_pgraster *pgrtn = NULL;
4463  rt_raster raster = NULL;
4464  rt_raster newrast = NULL;
4465  rt_band band = NULL;
4466  rt_band newband = NULL;
4467  int x, y, nband, width, height;
4468  double r;
4469  double newnodatavalue = 0.0;
4470  double newinitialvalue = 0.0;
4471  double newval = 0.0;
4472  char *newexpr = NULL;
4473  char *initexpr = NULL;
4474  char *expression = NULL;
4475  int hasnodataval = 0;
4476  double nodataval = 0.;
4477  rt_pixtype newpixeltype;
4478  int skipcomputation = 0;
4479  int len = 0;
4480  const int argkwcount = 3;
4481  enum KEYWORDS { kVAL=0, kX=1, kY=2 };
4482  char *argkw[] = {"[rast]", "[rast.x]", "[rast.y]"};
4483  Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
4484  int argcount = 0;
4485  Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
4486  uint8_t argpos[3] = {0};
4487  char place[12];
4488  int idx = 0;
4489  int ret = -1;
4490  TupleDesc tupdesc;
4491  SPIPlanPtr spi_plan = NULL;
4492  SPITupleTable * tuptable = NULL;
4493  HeapTuple tuple;
4494  char * strFromText = NULL;
4495  Datum *values = NULL;
4496  Datum datum = (Datum)NULL;
4497  char *nulls = NULL;
4498  bool isnull = FALSE;
4499  int i = 0;
4500  int j = 0;
4501 
4502  POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraExpr: Starting...");
4503 
4504  /* Check raster */
4505  if (PG_ARGISNULL(0)) {
4506  elog(NOTICE, "Raster is NULL. Returning NULL");
4507  PG_RETURN_NULL();
4508  }
4509 
4510 
4511  /* Deserialize raster */
4512  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4513  raster = rt_raster_deserialize(pgraster, FALSE);
4514  if (NULL == raster) {
4515  PG_FREE_IF_COPY(pgraster, 0);
4516  elog(ERROR, "RASTER_mapAlgebraExpr: Could not deserialize raster");
4517  PG_RETURN_NULL();
4518  }
4519 
4520  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting arguments...");
4521 
4522  if (PG_ARGISNULL(1))
4523  nband = 1;
4524  else
4525  nband = PG_GETARG_INT32(1);
4526 
4527  if (nband < 1)
4528  nband = 1;
4529 
4530 
4531  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new empty raster...");
4532 
4537  width = rt_raster_get_width(raster);
4538  height = rt_raster_get_height(raster);
4539 
4540  newrast = rt_raster_new(width, height);
4541 
4542  if ( NULL == newrast ) {
4543  PG_FREE_IF_COPY(pgraster, 0);
4544  elog(ERROR, "RASTER_mapAlgebraExpr: Could not create a new raster");
4545  PG_RETURN_NULL();
4546  }
4547 
4548  rt_raster_set_scale(newrast,
4551 
4552  rt_raster_set_offsets(newrast,
4555 
4556  rt_raster_set_skews(newrast,
4559 
4561 
4562 
4567  if (rt_raster_is_empty(newrast))
4568  {
4569  elog(NOTICE, "Raster is empty. Returning an empty raster");
4571  PG_FREE_IF_COPY(pgraster, 0);
4572 
4573  pgrtn = rt_raster_serialize(newrast);
4574  rt_raster_destroy(newrast);
4575  if (NULL == pgrtn) {
4576 
4577  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4578  PG_RETURN_NULL();
4579  }
4580 
4581  SET_VARSIZE(pgrtn, pgrtn->size);
4582  PG_RETURN_POINTER(pgrtn);
4583  }
4584 
4585 
4586  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Getting raster band %d...", nband);
4587 
4592  if (!rt_raster_has_band(raster, nband - 1)) {
4593  elog(NOTICE, "Raster does not have the required band. Returning a raster "
4594  "without a band");
4596  PG_FREE_IF_COPY(pgraster, 0);
4597 
4598  pgrtn = rt_raster_serialize(newrast);
4599  rt_raster_destroy(newrast);
4600  if (NULL == pgrtn) {
4601  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4602  PG_RETURN_NULL();
4603  }
4604 
4605  SET_VARSIZE(pgrtn, pgrtn->size);
4606  PG_RETURN_POINTER(pgrtn);
4607  }
4608 
4609  /* Get the raster band */
4611  if ( NULL == band ) {
4612  elog(NOTICE, "Could not get the required band. Returning a raster "
4613  "without a band");
4615  PG_FREE_IF_COPY(pgraster, 0);
4616 
4617  pgrtn = rt_raster_serialize(newrast);
4618  rt_raster_destroy(newrast);
4619  if (NULL == pgrtn) {
4620  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4621  PG_RETURN_NULL();
4622  }
4623 
4624  SET_VARSIZE(pgrtn, pgrtn->size);
4625  PG_RETURN_POINTER(pgrtn);
4626  }
4627 
4628  /*
4629  * Get NODATA value
4630  */
4631  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting NODATA value for band...");
4632 
4634  rt_band_get_nodata(band, &newnodatavalue);
4635  }
4636 
4637  else {
4638  newnodatavalue = rt_band_get_min_value(band);
4639  }
4640 
4641  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: NODATA value for band: = %f",
4642  newnodatavalue);
4643 
4649  newinitialvalue = newnodatavalue;
4650 
4654  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Setting pixeltype...");
4655 
4656  if (PG_ARGISNULL(2)) {
4657  newpixeltype = rt_band_get_pixtype(band);
4658  }
4659 
4660  else {
4661  strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
4662  newpixeltype = rt_pixtype_index_from_name(strFromText);
4663  pfree(strFromText);
4664  if (newpixeltype == PT_END)
4665  newpixeltype = rt_band_get_pixtype(band);
4666  }
4667 
4668  if (newpixeltype == PT_END) {
4669  PG_FREE_IF_COPY(pgraster, 0);
4670  elog(ERROR, "RASTER_mapAlgebraExpr: Invalid pixeltype");
4671  PG_RETURN_NULL();
4672  }
4673 
4674  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Pixeltype set to %s",
4675  rt_pixtype_name(newpixeltype));
4676 
4677 
4678  /* Construct expression for raster values */
4679  if (!PG_ARGISNULL(3)) {
4680  expression = text_to_cstring(PG_GETARG_TEXT_P(3));
4681  len = strlen("SELECT (") + strlen(expression) + strlen(")::double precision");
4682  initexpr = (char *)palloc(len + 1);
4683 
4684  memcpy(initexpr, "SELECT (", strlen("SELECT ("));
4685  memcpy(initexpr + strlen("SELECT ("), expression, strlen(expression));
4686  memcpy(initexpr + strlen("SELECT (") + strlen(expression), ")::double precision", strlen(")::double precision"));
4687  initexpr[len] = '\0';
4688 
4689  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression is %s", initexpr);
4690 
4691  /* We don't need this memory */
4692  /*
4693  pfree(expression);
4694  expression = NULL;
4695  */
4696  }
4697 
4698 
4699 
4705  if (!PG_ARGISNULL(4)) {
4706  hasnodataval = 1;
4707  nodataval = PG_GETARG_FLOAT8(4);
4708  newinitialvalue = nodataval;
4709 
4710  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new initial value = %f",
4711  newinitialvalue);
4712  }
4713  else
4714  hasnodataval = 0;
4715 
4716 
4717 
4724 
4725  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Band is a nodata band, returning "
4726  "a raster filled with nodata");
4727 
4728  ret = rt_raster_generate_new_band(newrast, newpixeltype,
4729  newinitialvalue, TRUE, newnodatavalue, 0);
4730 
4731  /* Free memory */
4732  if (initexpr)
4733  pfree(initexpr);
4735  PG_FREE_IF_COPY(pgraster, 0);
4736 
4737  /* Serialize created raster */
4738  pgrtn = rt_raster_serialize(newrast);
4739  rt_raster_destroy(newrast);
4740  if (NULL == pgrtn) {
4741  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4742  PG_RETURN_NULL();
4743  }
4744 
4745  SET_VARSIZE(pgrtn, pgrtn->size);
4746  PG_RETURN_POINTER(pgrtn);
4747  }
4748 
4749 
4754  if (initexpr != NULL && ( !strcmp(initexpr, "SELECT [rast]") || !strcmp(initexpr, "SELECT [rast.val]") ) && !hasnodataval) {
4755 
4756  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression resumes to RAST. "
4757  "Returning raster with band %d from original raster", nband);
4758 
4759  POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster has %d bands",
4760  rt_raster_get_num_bands(newrast));
4761 
4762  rt_raster_copy_band(newrast, raster, nband - 1, 0);
4763 
4764  POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster now has %d bands",
4765  rt_raster_get_num_bands(newrast));
4766 
4767  pfree(initexpr);
4769  PG_FREE_IF_COPY(pgraster, 0);
4770 
4771  /* Serialize created raster */
4772  pgrtn = rt_raster_serialize(newrast);
4773  rt_raster_destroy(newrast);
4774  if (NULL == pgrtn) {
4775  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4776  PG_RETURN_NULL();
4777  }
4778 
4779  SET_VARSIZE(pgrtn, pgrtn->size);
4780  PG_RETURN_POINTER(pgrtn);
4781  }
4782 
4787  if (initexpr != NULL && strstr(initexpr, "[rast") == NULL) {
4788  ret = SPI_connect();
4789  if (ret != SPI_OK_CONNECT) {
4790  PG_FREE_IF_COPY(pgraster, 0);
4791  elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4792  PG_RETURN_NULL();
4793  };
4794 
4795  /* Execute the expresion into newval */
4796  ret = SPI_execute(initexpr, FALSE, 0);
4797 
4798  if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
4799 
4800  /* Free memory allocated out of the current context */
4801  if (SPI_tuptable)
4802  SPI_freetuptable(tuptable);
4803  PG_FREE_IF_COPY(pgraster, 0);
4804 
4805  SPI_finish();
4806  elog(ERROR, "RASTER_mapAlgebraExpr: Invalid construction for expression");
4807  PG_RETURN_NULL();
4808  }
4809 
4810  tupdesc = SPI_tuptable->tupdesc;
4811  tuptable = SPI_tuptable;
4812 
4813  tuple = tuptable->vals[0];
4814  newexpr = SPI_getvalue(tuple, tupdesc, 1);
4815  if ( ! newexpr ) {
4816  POSTGIS_RT_DEBUG(3, "Constant expression evaluated to NULL, keeping initvalue");
4817  newval = newinitialvalue;
4818  } else {
4819  newval = atof(newexpr);
4820  }
4821 
4822  SPI_freetuptable(tuptable);
4823 
4824  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: New raster value = %f",
4825  newval);
4826 
4827  SPI_finish();
4828 
4829  skipcomputation = 1;
4830 
4835  if (!hasnodataval) {
4836  newinitialvalue = newval;
4837  skipcomputation = 2;
4838  }
4839 
4840  /* Return the new raster as it will be before computing pixel by pixel */
4841  else if (FLT_NEQ(newval, newinitialvalue)) {
4842  skipcomputation = 2;
4843  }
4844  }
4845 
4850  ret = rt_raster_generate_new_band(newrast, newpixeltype,
4851  newinitialvalue, TRUE, newnodatavalue, 0);
4852 
4857  /*if (initexpr == NULL || skipcomputation == 2) {*/
4858  if (expression == NULL || skipcomputation == 2) {
4859 
4860  /* Free memory */
4861  if (initexpr)
4862  pfree(initexpr);
4864  PG_FREE_IF_COPY(pgraster, 0);
4865 
4866  /* Serialize created raster */
4867  pgrtn = rt_raster_serialize(newrast);
4868  rt_raster_destroy(newrast);
4869  if (NULL == pgrtn) {
4870  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4871  PG_RETURN_NULL();
4872  }
4873 
4874  SET_VARSIZE(pgrtn, pgrtn->size);
4875  PG_RETURN_POINTER(pgrtn);
4876  }
4877 
4878  RASTER_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new raster band...");
4879 
4880  /* Get the new raster band */
4881  newband = rt_raster_get_band(newrast, 0);
4882  if ( NULL == newband ) {
4883  elog(NOTICE, "Could not modify band for new raster. Returning new "
4884  "raster with the original band");
4885 
4886  if (initexpr)
4887  pfree(initexpr);
4889  PG_FREE_IF_COPY(pgraster, 0);
4890 
4891  /* Serialize created raster */
4892  pgrtn = rt_raster_serialize(newrast);
4893  rt_raster_destroy(newrast);
4894  if (NULL == pgrtn) {
4895  elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
4896  PG_RETURN_NULL();
4897  }
4898 
4899  SET_VARSIZE(pgrtn, pgrtn->size);
4900  PG_RETURN_POINTER(pgrtn);
4901  }
4902 
4903  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Main computing loop (%d x %d)",
4904  width, height);
4905 
4906  if (initexpr != NULL) {
4907  /* Convert [rast.val] to [rast] */
4908  newexpr = rtpg_strreplace(initexpr, "[rast.val]", "[rast]", NULL);
4909  pfree(initexpr); initexpr=newexpr;
4910 
4911  sprintf(place,"$1");
4912  for (i = 0, j = 1; i < argkwcount; i++) {
4913  len = 0;
4914  newexpr = rtpg_strreplace(initexpr, argkw[i], place, &len);
4915  pfree(initexpr); initexpr=newexpr;
4916  if (len > 0) {
4917  argtype[argcount] = argkwtypes[i];
4918  argcount++;
4919  argpos[i] = j++;
4920 
4921  sprintf(place, "$%d", j);
4922  }
4923  else {
4924  argpos[i] = 0;
4925  }
4926  }
4927 
4928  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: initexpr = %s", initexpr);
4929 
4930  /* define values */
4931  values = (Datum *) palloc(sizeof(Datum) * argcount);
4932  if (values == NULL) {
4933 
4934  SPI_finish();
4935 
4937  PG_FREE_IF_COPY(pgraster, 0);
4938  rt_raster_destroy(newrast);
4939 
4940  elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
4941  PG_RETURN_NULL();
4942  }
4943 
4944  /* define nulls */
4945  nulls = (char *)palloc(argcount);
4946  if (nulls == NULL) {
4947 
4948  SPI_finish();
4949 
4951  PG_FREE_IF_COPY(pgraster, 0);
4952  rt_raster_destroy(newrast);
4953 
4954  elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
4955  PG_RETURN_NULL();
4956  }
4957 
4958  /* Connect to SPI and prepare the expression */
4959  ret = SPI_connect();
4960  if (ret != SPI_OK_CONNECT) {
4961 
4962  if (initexpr)
4963  pfree(initexpr);
4965  PG_FREE_IF_COPY(pgraster, 0);
4966  rt_raster_destroy(newrast);
4967 
4968  elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4969  PG_RETURN_NULL();
4970  };
4971 
4972  /* Type of all arguments is FLOAT8OID */
4973  spi_plan = SPI_prepare(initexpr, argcount, argtype);
4974 
4975  if (spi_plan == NULL) {
4976 
4978  PG_FREE_IF_COPY(pgraster, 0);
4979  rt_raster_destroy(newrast);
4980 
4981  SPI_finish();
4982 
4983  pfree(initexpr);
4984 
4985  elog(ERROR, "RASTER_mapAlgebraExpr: Could not prepare expression");
4986  PG_RETURN_NULL();
4987  }
4988  }
4989 
4990  for (x = 0; x < width; x++) {
4991  for(y = 0; y < height; y++) {
4992  ret = rt_band_get_pixel(band, x, y, &r, NULL);
4993 
4998  if (ret == ES_NONE && FLT_NEQ(r, newnodatavalue)) {
4999  if (skipcomputation == 0) {
5000  if (initexpr != NULL) {
5001  /* Reset the null arg flags. */
5002  memset(nulls, 'n', argcount);
5003 
5004  for (i = 0; i < argkwcount; i++) {
5005  idx = argpos[i];
5006  if (idx < 1) continue;
5007  idx--;
5008 
5009  if (i == kX ) {
5010  /* x is 0 based index, but SQL expects 1 based index */
5011  values[idx] = Int32GetDatum(x+1);
5012  nulls[idx] = ' ';
5013  }
5014  else if (i == kY) {
5015  /* y is 0 based index, but SQL expects 1 based index */
5016  values[idx] = Int32GetDatum(y+1);
5017  nulls[idx] = ' ';
5018  }
5019  else if (i == kVAL ) {
5020  values[idx] = Float8GetDatum(r);
5021  nulls[idx] = ' ';
5022  }
5023 
5024  }
5025 
5026  ret = SPI_execute_plan(spi_plan, values, nulls, FALSE, 0);
5027  if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
5028  SPI_processed != 1) {
5029  if (SPI_tuptable)
5030  SPI_freetuptable(tuptable);
5031 
5032  SPI_freeplan(spi_plan);
5033  SPI_finish();
5034 
5035  pfree(values);
5036  pfree(nulls);
5037  pfree(initexpr);
5038 
5040  PG_FREE_IF_COPY(pgraster, 0);
5041  rt_raster_destroy(newrast);
5042 
5043  elog(ERROR, "RASTER_mapAlgebraExpr: Error executing prepared plan");
5044 
5045  PG_RETURN_NULL();
5046  }
5047 
5048  tupdesc = SPI_tuptable->tupdesc;
5049  tuptable = SPI_tuptable;
5050 
5051  tuple = tuptable->vals[0];
5052  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
5053  if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
5054  POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) errored, skip setting", x+1,y+1,r);
5055  newval = newinitialvalue;
5056  }
5057  else if ( isnull ) {
5058  POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting", x+1,y+1,r);
5059  newval = newinitialvalue;
5060  } else {
5061  newval = DatumGetFloat8(datum);
5062  }
5063 
5064  SPI_freetuptable(tuptable);
5065  }
5066 
5067  else
5068  newval = newinitialvalue;
5069 
5070  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new value = %f",
5071  newval);
5072  }
5073 
5074 
5075  rt_band_set_pixel(newband, x, y, newval, NULL);
5076  }
5077 
5078  }
5079  }
5080 
5081  if (initexpr != NULL) {
5082  SPI_freeplan(spi_plan);
5083  SPI_finish();
5084 
5085  pfree(values);
5086  pfree(nulls);
5087  pfree(initexpr);
5088  }
5089  else {
5090  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: no SPI cleanup");
5091  }
5092 
5093 
5094  /* The newrast band has been modified */
5095 
5096  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster modified, serializing it.");
5097  /* Serialize created raster */
5098 
5100  PG_FREE_IF_COPY(pgraster, 0);
5101 
5102  pgrtn = rt_raster_serialize(newrast);
5103  rt_raster_destroy(newrast);
5104  if (NULL == pgrtn)
5105  PG_RETURN_NULL();
5106 
5107  SET_VARSIZE(pgrtn, pgrtn->size);
5108 
5109  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster serialized");
5110 
5111 
5112  POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraExpr: returning raster");
5113 
5114 
5115  PG_RETURN_POINTER(pgrtn);
5116 }
5117 
5118 /*
5119  * One raster user function MapAlgebra.
5120  */
5122 Datum RASTER_mapAlgebraFct(PG_FUNCTION_ARGS)
5123 {
5124  rt_pgraster *pgraster = NULL;
5125  rt_pgraster *pgrtn = NULL;
5126  rt_raster raster = NULL;
5127  rt_raster newrast = NULL;
5128  rt_band band = NULL;
5129  rt_band newband = NULL;
5130  int x, y, nband, width, height;
5131  double r;
5132  double newnodatavalue = 0.0;
5133  double newinitialvalue = 0.0;
5134  double newval = 0.0;
5135  rt_pixtype newpixeltype;
5136  int ret = -1;
5137  Oid oid;
5138  FmgrInfo cbinfo;
5139 #if POSTGIS_PGSQL_VERSION < 120
5140  FunctionCallInfoData cbdata;
5141 #else
5142  LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS); /* Could be optimized */
5143 #endif
5144  Datum tmpnewval;
5145  char * strFromText = NULL;
5146  int k = 0;
5147 
5148  POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFct: STARTING...");
5149 
5150  /* Check raster */
5151  if (PG_ARGISNULL(0)) {
5152  elog(WARNING, "Raster is NULL. Returning NULL");
5153  PG_RETURN_NULL();
5154  }
5155 
5156 
5157  /* Deserialize raster */
5158  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5159  raster = rt_raster_deserialize(pgraster, FALSE);
5160  if (NULL == raster) {
5161  PG_FREE_IF_COPY(pgraster, 0);
5162  elog(ERROR, "RASTER_mapAlgebraFct: Could not deserialize raster");
5163  PG_RETURN_NULL();
5164  }
5165 
5166  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting arguments...");
5167 
5168  /* Get the rest of the arguments */
5169 
5170  if (PG_ARGISNULL(1))
5171  nband = 1;
5172  else
5173  nband = PG_GETARG_INT32(1);
5174 
5175  if (nband < 1)
5176  nband = 1;
5177 
5178  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Creating new empty raster...");
5179 
5184  width = rt_raster_get_width(raster);
5185  height = rt_raster_get_height(raster);
5186 
5187  newrast = rt_raster_new(width, height);
5188 
5189  if ( NULL == newrast ) {
5190 
5192  PG_FREE_IF_COPY(pgraster, 0);
5193 
5194  elog(ERROR, "RASTER_mapAlgebraFct: Could not create a new raster");
5195  PG_RETURN_NULL();
5196  }
5197 
5198  rt_raster_set_scale(newrast,
5201 
5202  rt_raster_set_offsets(newrast,
5205 
5206  rt_raster_set_skews(newrast,
5209 
5211 
5212 
5217  if (rt_raster_is_empty(newrast))
5218  {
5219  elog(NOTICE, "Raster is empty. Returning an empty raster");
5221  PG_FREE_IF_COPY(pgraster, 0);
5222 
5223  pgrtn = rt_raster_serialize(newrast);
5224  rt_raster_destroy(newrast);
5225  if (NULL == pgrtn) {
5226  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5227  PG_RETURN_NULL();
5228  }
5229 
5230  SET_VARSIZE(pgrtn, pgrtn->size);
5231  PG_RETURN_POINTER(pgrtn);
5232  }
5233 
5234  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Getting raster band %d...", nband);
5235 
5240  if (!rt_raster_has_band(raster, nband - 1)) {
5241  elog(NOTICE, "Raster does not have the required band. Returning a raster "
5242  "without a band");
5244  PG_FREE_IF_COPY(pgraster, 0);
5245 
5246  pgrtn = rt_raster_serialize(newrast);
5247  rt_raster_destroy(newrast);
5248  if (NULL == pgrtn) {
5249  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5250  PG_RETURN_NULL();
5251  }
5252 
5253  SET_VARSIZE(pgrtn, pgrtn->size);
5254  PG_RETURN_POINTER(pgrtn);
5255  }
5256 
5257  /* Get the raster band */
5259  if ( NULL == band ) {
5260  elog(NOTICE, "Could not get the required band. Returning a raster "
5261  "without a band");
5263  PG_FREE_IF_COPY(pgraster, 0);
5264 
5265  pgrtn = rt_raster_serialize(newrast);
5266  rt_raster_destroy(newrast);
5267  if (NULL == pgrtn) {
5268  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5269  PG_RETURN_NULL();
5270  }
5271 
5272  SET_VARSIZE(pgrtn, pgrtn->size);
5273  PG_RETURN_POINTER(pgrtn);
5274  }
5275 
5276  /*
5277  * Get NODATA value
5278  */
5279  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting NODATA value for band...");
5280 
5282  rt_band_get_nodata(band, &newnodatavalue);
5283  }
5284 
5285  else {
5286  newnodatavalue = rt_band_get_min_value(band);
5287  }
5288 
5289  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: NODATA value for band: %f",
5290  newnodatavalue);
5296  newinitialvalue = newnodatavalue;
5297 
5301  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Setting pixeltype...");
5302 
5303  if (PG_ARGISNULL(2)) {
5304  newpixeltype = rt_band_get_pixtype(band);
5305  }
5306 
5307  else {
5308  strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5309  newpixeltype = rt_pixtype_index_from_name(strFromText);
5310  pfree(strFromText);
5311  if (newpixeltype == PT_END)
5312  newpixeltype = rt_band_get_pixtype(band);
5313  }
5314 
5315  if (newpixeltype == PT_END) {
5316 
5318  PG_FREE_IF_COPY(pgraster, 0);
5319  rt_raster_destroy(newrast);
5320 
5321  elog(ERROR, "RASTER_mapAlgebraFct: Invalid pixeltype");
5322  PG_RETURN_NULL();
5323  }
5324 
5325  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Pixeltype set to %s",
5326  rt_pixtype_name(newpixeltype));
5327 
5328  /* Get the name of the callback user function for raster values */
5329  if (PG_ARGISNULL(3)) {
5330 
5332  PG_FREE_IF_COPY(pgraster, 0);
5333  rt_raster_destroy(newrast);
5334 
5335  elog(ERROR, "RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
5336  PG_RETURN_NULL();
5337  }
5338 
5339  oid = PG_GETARG_OID(3);
5340  if (oid == InvalidOid) {
5341 
5343  PG_FREE_IF_COPY(pgraster, 0);
5344  rt_raster_destroy(newrast);
5345 
5346  elog(ERROR, "RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
5347  PG_RETURN_NULL();
5348  }
5349 
5350  fmgr_info(oid, &cbinfo);
5351 
5352  /* function cannot return set */
5353  if (cbinfo.fn_retset) {
5354 
5356  PG_FREE_IF_COPY(pgraster, 0);
5357  rt_raster_destroy(newrast);
5358 
5359  elog(ERROR, "RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
5360  PG_RETURN_NULL();
5361  }
5362  /* function should have correct # of args */
5363  else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {
5364 
5366  PG_FREE_IF_COPY(pgraster, 0);
5367  rt_raster_destroy(newrast);
5368 
5369  elog(ERROR, "RASTER_mapAlgebraFct: Function does not have two or three input parameters");
5370  PG_RETURN_NULL();
5371  }
5372 
5373  if (cbinfo.fn_nargs == 2)
5374  k = 1;
5375  else
5376  k = 2;
5377 
5378  if (func_volatile(oid) == 'v') {
5379  elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5380  }
5381 
5382  /* prep function call data */
5383 #if POSTGIS_PGSQL_VERSION < 120
5384  InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5385 
5386  cbdata.argnull[0] = FALSE;
5387  cbdata.argnull[1] = FALSE;
5388  cbdata.argnull[2] = FALSE;
5389 #else
5390  InitFunctionCallInfoData(*cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5391 
5392  cbdata->args[0].isnull = FALSE;
5393  cbdata->args[1].isnull = FALSE;
5394  cbdata->args[2].isnull = FALSE;
5395 #endif
5396 
5397  /* check that the function isn't strict if the args are null. */
5398  if (PG_ARGISNULL(4)) {
5399  if (cbinfo.fn_strict) {
5400 
5402  PG_FREE_IF_COPY(pgraster, 0);
5403  rt_raster_destroy(newrast);
5404 
5405  elog(ERROR, "RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
5406  PG_RETURN_NULL();
5407  }
5408 
5409 #if POSTGIS_PGSQL_VERSION < 120
5410  cbdata.arg[k] = (Datum)NULL;
5411  cbdata.argnull[k] = TRUE;
5412 #else
5413  cbdata->args[k].value = (Datum)NULL;
5414  cbdata->args[k].isnull = TRUE;
5415 #endif
5416  }
5417  else {
5418 #if POSTGIS_PGSQL_VERSION < 120
5419  cbdata.arg[k] = PG_GETARG_DATUM(4);
5420 #else
5421  cbdata->args[k].value = PG_GETARG_DATUM(4);
5422 #endif
5423  }
5424 
5431 
5432  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Band is a nodata band, returning "
5433  "a raster filled with nodata");
5434 
5435  ret = rt_raster_generate_new_band(newrast, newpixeltype,
5436  newinitialvalue, TRUE, newnodatavalue, 0);
5437 
5439  PG_FREE_IF_COPY(pgraster, 0);
5440 
5441  /* Serialize created raster */
5442  pgrtn = rt_raster_serialize(newrast);
5443  rt_raster_destroy(newrast);
5444  if (NULL == pgrtn) {
5445  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5446  PG_RETURN_NULL();
5447  }
5448 
5449  SET_VARSIZE(pgrtn, pgrtn->size);
5450  PG_RETURN_POINTER(pgrtn);
5451  }
5452 
5453 
5458  ret = rt_raster_generate_new_band(newrast, newpixeltype,
5459  newinitialvalue, TRUE, newnodatavalue, 0);
5460 
5461  /* Get the new raster band */
5462  newband = rt_raster_get_band(newrast, 0);
5463  if ( NULL == newband ) {
5464  elog(NOTICE, "Could not modify band for new raster. Returning new "
5465  "raster with the original band");
5466 
5468  PG_FREE_IF_COPY(pgraster, 0);
5469 
5470  /* Serialize created raster */
5471  pgrtn = rt_raster_serialize(newrast);
5472  rt_raster_destroy(newrast);
5473  if (NULL == pgrtn) {
5474  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5475  PG_RETURN_NULL();
5476  }
5477 
5478  SET_VARSIZE(pgrtn, pgrtn->size);
5479  PG_RETURN_POINTER(pgrtn);
5480  }
5481 
5482  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Main computing loop (%d x %d)",
5483  width, height);
5484 
5485  for (x = 0; x < width; x++) {
5486  for(y = 0; y < height; y++) {
5487  ret = rt_band_get_pixel(band, x, y, &r, NULL);
5488 
5493  if (ret == ES_NONE) {
5494  if (FLT_EQ(r, newnodatavalue)) {
5495  if (cbinfo.fn_strict) {
5496  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
5497  continue;
5498  }
5499 #if POSTGIS_PGSQL_VERSION < 120
5500  cbdata.argnull[0] = TRUE;
5501  cbdata.arg[0] = (Datum)NULL;
5502 #else
5503  cbdata->args[0].isnull = TRUE;
5504  cbdata->args[0].value = (Datum)NULL;
5505 #endif
5506  }
5507  else {
5508 #if POSTGIS_PGSQL_VERSION < 120
5509  cbdata.argnull[0] = FALSE;
5510  cbdata.arg[0] = Float8GetDatum(r);
5511 #else
5512  cbdata->args[0].isnull = FALSE;
5513  cbdata->args[0].value = Float8GetDatum(r);
5514 #endif
5515  }
5516 
5517  /* Add pixel positions if callback has proper # of args */
5518  if (cbinfo.fn_nargs == 3) {
5519  Datum d[2];
5520  ArrayType *a;
5521 
5522  d[0] = Int32GetDatum(x+1);
5523  d[1] = Int32GetDatum(y+1);
5524 
5525  a = construct_array(d, 2, INT4OID, sizeof(int32), true, 'i');
5526 
5527 #if POSTGIS_PGSQL_VERSION < 120
5528  cbdata.argnull[1] = FALSE;
5529  cbdata.arg[1] = PointerGetDatum(a);
5530 #else
5531  cbdata->args[1].isnull = FALSE;
5532  cbdata->args[1].value = PointerGetDatum(a);
5533 #endif
5534  }
5535 
5536  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: (%dx%d), r = %f",
5537  x, y, r);
5538 
5539 #if POSTGIS_PGSQL_VERSION < 120
5540  tmpnewval = FunctionCallInvoke(&cbdata);
5541 
5542  if (cbdata.isnull) {
5543  newval = newnodatavalue;
5544  }
5545 #else
5546  tmpnewval = FunctionCallInvoke(cbdata);
5547 
5548  if (cbdata->isnull)
5549  {
5550  newval = newnodatavalue;
5551  }
5552 #endif
5553  else {
5554  newval = DatumGetFloat8(tmpnewval);
5555  }
5556 
5557  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: new value = %f",
5558  newval);
5559 
5560  rt_band_set_pixel(newband, x, y, newval, NULL);
5561  }
5562 
5563  }
5564  }
5565 
5566  /* The newrast band has been modified */
5567 
5568  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster modified, serializing it.");
5569  /* Serialize created raster */
5570 
5572  PG_FREE_IF_COPY(pgraster, 0);
5573 
5574  pgrtn = rt_raster_serialize(newrast);
5575  rt_raster_destroy(newrast);
5576  if (NULL == pgrtn)
5577  PG_RETURN_NULL();
5578 
5579  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster serialized");
5580 
5581  POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFct: returning raster");
5582 
5583  SET_VARSIZE(pgrtn, pgrtn->size);
5584  PG_RETURN_POINTER(pgrtn);
5585 }
5586 
5591 Datum RASTER_mapAlgebraFctNgb(PG_FUNCTION_ARGS)
5592 {
5593  rt_pgraster *pgraster = NULL;
5594  rt_pgraster *pgrtn = NULL;
5595  rt_raster raster = NULL;
5596  rt_raster newrast = NULL;
5597  rt_band band = NULL;
5598  rt_band newband = NULL;
5599  int x, y, nband, width, height, ngbwidth, ngbheight, winwidth, winheight, u, v, nIndex, nNullItems;
5600  double r, rpix;
5601  double newnodatavalue = 0.0;
5602  double newinitialvalue = 0.0;
5603  double newval = 0.0;
5604  rt_pixtype newpixeltype;
5605  int ret = -1;
5606  Oid oid;
5607  FmgrInfo cbinfo;
5608 #if POSTGIS_PGSQL_VERSION < 120
5609  FunctionCallInfoData cbdata;
5610 #else
5611  LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS); /* Could be optimized */
5612 #endif
5613  Datum tmpnewval;
5614  ArrayType * neighborDatum;
5615  char * strFromText = NULL;
5616  text * txtNodataMode = NULL;
5617  text * txtCallbackParam = NULL;
5618  int intReplace = 0;
5619  float fltReplace = 0;
5620  bool valuereplace = false, pixelreplace, nNodataOnly = true, nNullSkip = false;
5621  Datum * neighborData = NULL;
5622  bool * neighborNulls = NULL;
5623  int neighborDims[2];
5624  int neighborLbs[2];
5625  int16 typlen;
5626  bool typbyval;
5627  char typalign;
5628 
5629  POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFctNgb: STARTING...");
5630 
5631  /* Check raster */
5632  if (PG_ARGISNULL(0)) {
5633  elog(WARNING, "Raster is NULL. Returning NULL");
5634  PG_RETURN_NULL();
5635  }
5636 
5637 
5638  /* Deserialize raster */
5639  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5640  raster = rt_raster_deserialize(pgraster, FALSE);
5641  if (NULL == raster)
5642  {
5643  PG_FREE_IF_COPY(pgraster, 0);
5644  elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not deserialize raster");
5645  PG_RETURN_NULL();
5646  }
5647 
5648  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Getting arguments...");
5649 
5650  /* Get the rest of the arguments */
5651 
5652  if (PG_ARGISNULL(1))
5653  nband = 1;
5654  else
5655  nband = PG_GETARG_INT32(1);
5656 
5657  if (nband < 1)
5658  nband = 1;
5659 
5660  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Creating new empty raster...");
5661 
5666  width = rt_raster_get_width(raster);
5667  height = rt_raster_get_height(raster);
5668 
5669  newrast = rt_raster_new(width, height);
5670 
5671  if ( NULL == newrast ) {
5673  PG_FREE_IF_COPY(pgraster, 0);
5674  elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not create a new raster");
5675  PG_RETURN_NULL();
5676  }
5677 
5678  rt_raster_set_scale(newrast,
5681 
5682  rt_raster_set_offsets(newrast,
5685 
5686  rt_raster_set_skews(newrast,
5689 
5691 
5692 
5697  if (rt_raster_is_empty(newrast))
5698  {
5699  elog(NOTICE, "Raster is empty. Returning an empty raster");
5701  PG_FREE_IF_COPY(pgraster, 0);
5702 
5703  pgrtn = rt_raster_serialize(newrast);
5704  rt_raster_destroy(newrast);
5705  if (NULL == pgrtn) {
5706  elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5707  PG_RETURN_NULL();
5708  }
5709 
5710  SET_VARSIZE(pgrtn, pgrtn->size);
5711  PG_RETURN_POINTER(pgrtn);
5712  }
5713 
5714  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Getting raster band %d...", nband);
5715 
5720  if (!rt_raster_has_band(raster, nband - 1)) {
5721  elog(NOTICE, "Raster does not have the required band. Returning a raster "
5722  "without a band");
5724  PG_FREE_IF_COPY(pgraster, 0);
5725 
5726  pgrtn = rt_raster_serialize(newrast);
5727  rt_raster_destroy(newrast);
5728  if (NULL == pgrtn) {
5729  elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5730  PG_RETURN_NULL();
5731  }
5732 
5733  SET_VARSIZE(pgrtn, pgrtn->size);
5734  PG_RETURN_POINTER(pgrtn);
5735  }
5736 
5737  /* Get the raster band */
5739  if ( NULL == band ) {
5740  elog(NOTICE, "Could not get the required band. Returning a raster "
5741  "without a band");
5743  PG_FREE_IF_COPY(pgraster, 0);
5744 
5745  pgrtn = rt_raster_serialize(newrast);
5746  rt_raster_destroy(newrast);
5747  if (NULL == pgrtn) {
5748  elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5749  PG_RETURN_NULL();
5750  }
5751 
5752  SET_VARSIZE(pgrtn, pgrtn->size);
5753  PG_RETURN_POINTER(pgrtn);
5754  }
5755 
5756  /*
5757  * Get NODATA value
5758  */
5759  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Getting NODATA value for band...");
5760 
5762  rt_band_get_nodata(band, &newnodatavalue);
5763  }
5764 
5765  else {
5766  newnodatavalue = rt_band_get_min_value(band);
5767  }
5768 
5769  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: NODATA value for band: %f",
5770  newnodatavalue);
5776  newinitialvalue = newnodatavalue;
5777 
5781  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Setting pixeltype...");
5782 
5783  if (PG_ARGISNULL(2)) {
5784  newpixeltype = rt_band_get_pixtype(band);
5785  }
5786 
5787  else {
5788  strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5789  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Pixeltype parameter: %s", strFromText);
5790  newpixeltype = rt_pixtype_index_from_name(strFromText);
5791  pfree(strFromText);
5792  if (newpixeltype == PT_END)
5793  newpixeltype = rt_band_get_pixtype(band);
5794  }
5795 
5796  if (newpixeltype == PT_END) {
5797 
5799  PG_FREE_IF_COPY(pgraster, 0);
5800  rt_raster_destroy(newrast);
5801 
5802  elog(ERROR, "RASTER_mapAlgebraFctNgb: Invalid pixeltype");
5803  PG_RETURN_NULL();
5804  }
5805 
5806  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Pixeltype set to %s (%d)",
5807  rt_pixtype_name(newpixeltype), newpixeltype);
5808 
5809  /* Get the name of the callback userfunction */
5810  if (PG_ARGISNULL(5)) {
5811 
5813  PG_FREE_IF_COPY(pgraster, 0);
5814  rt_raster_destroy(newrast);
5815 
5816  elog(ERROR, "RASTER_mapAlgebraFctNgb: Required function is missing");
5817  PG_RETURN_NULL();
5818  }
5819 
5820  oid = PG_GETARG_OID(5);
5821  if (oid == InvalidOid) {
5822 
5824  PG_FREE_IF_COPY(pgraster, 0);
5825  rt_raster_destroy(newrast);
5826 
5827  elog(ERROR, "RASTER_mapAlgebraFctNgb: Got invalid function object id");
5828  PG_RETURN_NULL();
5829  }
5830 
5831  fmgr_info(oid, &cbinfo);
5832 
5833  /* function cannot return set */
5834  if (cbinfo.fn_retset) {
5835 
5837  PG_FREE_IF_COPY(pgraster, 0);
5838  rt_raster_destroy(newrast);
5839 
5840  elog(ERROR, "RASTER_mapAlgebraFctNgb: Function provided must return double precision not resultset");
5841  PG_RETURN_NULL();
5842  }
5843  /* function should have correct # of args */
5844  else if (cbinfo.fn_nargs != 3) {
5845 
5847  PG_FREE_IF_COPY(pgraster, 0);
5848  rt_raster_destroy(newrast);
5849 
5850  elog(ERROR, "RASTER_mapAlgebraFctNgb: Function does not have three input parameters");
5851  PG_RETURN_NULL();
5852  }
5853 
5854  if (func_volatile(oid) == 'v') {
5855  elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5856  }
5857 
5858  /* prep function call data */
5859 #if POSTGIS_PGSQL_VERSION < 120
5860  InitFunctionCallInfoData(cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5861  memset(cbdata.argnull, FALSE, sizeof(bool) * 3);
5862 #else
5863  InitFunctionCallInfoData(*cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5864  cbdata->args[0].isnull = FALSE;
5865  cbdata->args[1].isnull = FALSE;
5866  cbdata->args[2].isnull = FALSE;
5867 #endif
5868 
5869  /* check that the function isn't strict if the args are null. */
5870  if (PG_ARGISNULL(7)) {
5871  if (cbinfo.fn_strict) {
5872 
5874  PG_FREE_IF_COPY(pgraster, 0);
5875  rt_raster_destroy(newrast);
5876 
5877  elog(ERROR, "RASTER_mapAlgebraFctNgb: Strict callback functions cannot have NULL parameters");
5878  PG_RETURN_NULL();
5879  }
5880 
5881 #if POSTGIS_PGSQL_VERSION < 120
5882  cbdata.arg[2] = (Datum)NULL;
5883  cbdata.argnull[2] = TRUE;
5884 #else
5885  cbdata->args[2].value = (Datum)NULL;
5886  cbdata->args[2].isnull = TRUE;
5887 #endif
5888  }
5889  else {
5890 #if POSTGIS_PGSQL_VERSION < 120
5891  cbdata.arg[2] = PG_GETARG_DATUM(7);
5892 #else
5893  cbdata->args[2].value = PG_GETARG_DATUM(7);
5894 #endif
5895  }
5896 
5903 
5904  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Band is a nodata band, returning "
5905  "a raster filled with nodata");
5906 
5907  rt_raster_generate_new_band(newrast, newpixeltype,
5908  newinitialvalue, TRUE, newnodatavalue, 0);
5909 
5911  PG_FREE_IF_COPY(pgraster, 0);
5912 
5913  /* Serialize created raster */
5914  pgrtn = rt_raster_serialize(newrast);
5915  rt_raster_destroy(newrast);
5916  if (NULL == pgrtn) {
5917  elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5918  PG_RETURN_NULL();
5919  }
5920 
5921  SET_VARSIZE(pgrtn, pgrtn->size);
5922  PG_RETURN_POINTER(pgrtn);
5923  }
5924 
5925 
5930  rt_raster_generate_new_band(newrast, newpixeltype,
5931  newinitialvalue, TRUE, newnodatavalue, 0);
5932 
5933  /* Get the new raster band */
5934  newband = rt_raster_get_band(newrast, 0);
5935  if ( NULL == newband ) {
5936  elog(NOTICE, "Could not modify band for new raster. Returning new "
5937  "raster with the original band");
5938 
5940  PG_FREE_IF_COPY(pgraster, 0);
5941 
5942  /* Serialize created raster */
5943  pgrtn = rt_raster_serialize(newrast);
5944  rt_raster_destroy(newrast);
5945  if (NULL == pgrtn) {
5946  elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5947  PG_RETURN_NULL();
5948  }
5949 
5950  SET_VARSIZE(pgrtn, pgrtn->size);
5951  PG_RETURN_POINTER(pgrtn);
5952  }
5953 
5954  /* Get the width of the neighborhood */
5955  if (PG_ARGISNULL(3) || PG_GETARG_INT32(3) <= 0) {
5956  elog(NOTICE, "Neighborhood width is NULL or <= 0. Returning new "
5957  "raster with the original band");
5958 
5960  PG_FREE_IF_COPY(pgraster, 0);
5961 
5962  /* Serialize created raster */
5963  pgrtn = rt_raster_serialize(newrast);
5964  rt_raster_destroy(newrast);
5965  if (NULL == pgrtn) {
5966  elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5967  PG_RETURN_NULL();
5968  }
5969 
5970  SET_VARSIZE(pgrtn, pgrtn->size);
5971  PG_RETURN_POINTER(pgrtn);
5972  }
5973 
5974  ngbwidth = PG_GETARG_INT32(3);
5975  winwidth = ngbwidth * 2 + 1;
5976 
5977  /* Get the height of the neighborhood */
5978  if (PG_ARGISNULL(4) || PG_GETARG_INT32(4) <= 0) {
5979  elog(NOTICE, "Neighborhood height is NULL or <= 0. Returning new "
5980  "raster with the original band");
5981 
5983  PG_FREE_IF_COPY(pgraster, 0);
5984 
5985  /* Serialize created raster */
5986  pgrtn = rt_raster_serialize(newrast);
5987  rt_raster_destroy(newrast);
5988  if (NULL == pgrtn) {
5989  elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
5990  PG_RETURN_NULL();
5991  }
5992 
5993  SET_VARSIZE(pgrtn, pgrtn->size);
5994  PG_RETURN_POINTER(pgrtn);
5995  }
5996 
5997  ngbheight = PG_GETARG_INT32(4);
5998  winheight = ngbheight * 2 + 1;
5999 
6000  /* Get the type of NODATA behavior for the neighborhoods. */
6001  if (PG_ARGISNULL(6)) {
6002  elog(NOTICE, "Neighborhood NODATA behavior defaulting to 'ignore'");
6003  txtNodataMode = cstring_to_text("ignore");
6004  }
6005  else {
6006  txtNodataMode = PG_GETARG_TEXT_P(6);
6007  }
6008 
6009  txtCallbackParam = (text*)palloc(VARSIZE(txtNodataMode));
6010  SET_VARSIZE(txtCallbackParam, VARSIZE(txtNodataMode));
6011  memcpy((void *)VARDATA(txtCallbackParam), (void *)VARDATA(txtNodataMode), VARSIZE(txtNodataMode) - VARHDRSZ);
6012 
6013  /* pass the nodata mode into the user function */
6014 #if POSTGIS_PGSQL_VERSION < 120
6015  cbdata.arg[1] = PointerGetDatum(txtCallbackParam);
6016 #else
6017  cbdata->args[1].value = PointerGetDatum(txtCallbackParam);
6018 #endif
6019 
6020  strFromText = text_to_cstring(txtNodataMode);
6021  strFromText = rtpg_strtoupper(strFromText);
6022 
6023  if (strcmp(strFromText, "VALUE") == 0)
6024  valuereplace = true;
6025  else if (strcmp(strFromText, "IGNORE") != 0 && strcmp(strFromText, "NULL") != 0) {
6026  /* if the text is not "IGNORE" or "NULL", it may be a numerical value */
6027  if (sscanf(strFromText, "%d", &intReplace) <= 0 && sscanf(strFromText, "%f", &fltReplace) <= 0) {
6028  /* the value is NOT an integer NOR a floating point */
6029  elog(NOTICE, "Neighborhood NODATA mode is not recognized. Must be one of 'value', 'ignore', "
6030  "'NULL', or a numeric value. Returning new raster with the original band");
6031 
6032  /* clean up the nodatamode string */
6033  pfree(txtCallbackParam);
6034  pfree(strFromText);
6035 
6037  PG_FREE_IF_COPY(pgraster, 0);
6038 
6039  /* Serialize created raster */
6040  pgrtn = rt_raster_serialize(newrast);
6041  rt_raster_destroy(newrast);
6042  if (NULL == pgrtn) {
6043  elog(ERROR, "RASTER_mapAlgebraFctNgb: Could not serialize raster");
6044  PG_RETURN_NULL();
6045  }
6046 
6047  SET_VARSIZE(pgrtn, pgrtn->size);
6048  PG_RETURN_POINTER(pgrtn);
6049  }
6050  }
6051  else if (strcmp(strFromText, "NULL") == 0) {
6052  /* this setting means that the neighborhood should be skipped if any of the values are null */
6053  nNullSkip = true;
6054  }
6055 
6056  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Main computing loop (%d x %d)",
6057  width, height);
6058 
6059  /* Allocate room for the neighborhood. */
6060  neighborData = (Datum *)palloc(winwidth * winheight * sizeof(Datum));
6061  neighborNulls = (bool *)palloc(winwidth * winheight * sizeof(bool));
6062 
6063  /* The dimensions of the neighborhood array, for creating a multi-dimensional array. */
6064  neighborDims[0] = winwidth;
6065  neighborDims[1] = winheight;
6066 
6067  /* The lower bounds for the new multi-dimensional array. */
6068  neighborLbs[0] = 1;
6069  neighborLbs[1] = 1;
6070 
6071  /* Get information about the type of item in the multi-dimensional array (float8). */
6072  get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
6073 
6074  for (x = 0 + ngbwidth; x < width - ngbwidth; x++) {
6075  for(y = 0 + ngbheight; y < height - ngbheight; y++) {
6076  /* populate an array with the pixel values in the neighborhood */
6077  nIndex = 0;
6078  nNullItems = 0;
6079  nNodataOnly = true;
6080  pixelreplace = false;
6081  if (valuereplace) {
6082  ret = rt_band_get_pixel(band, x, y, &rpix, NULL);
6083  if (ret == ES_NONE && FLT_NEQ(rpix, newnodatavalue)) {
6084  pixelreplace = true;
6085  }
6086  }
6087  for (u = x - ngbwidth; u <= x + ngbwidth; u++) {
6088  for (v = y - ngbheight; v <= y + ngbheight; v++) {
6089  ret = rt_band_get_pixel(band, u, v, &r, NULL);
6090  if (ret == ES_NONE) {
6091  if (FLT_NEQ(r, newnodatavalue)) {
6092  /* If the pixel value for this neighbor cell is not NODATA */
6093  neighborData[nIndex] = Float8GetDatum((double)r);
6094  neighborNulls[nIndex] = false;
6095  nNodataOnly = false;
6096  }
6097  else {
6098  /* If the pixel value for this neighbor cell is NODATA */
6099  if (valuereplace && pixelreplace) {
6100  /* Replace the NODATA value with the currently processing pixel. */
6101  neighborData[nIndex] = Float8GetDatum((double)rpix);
6102  neighborNulls[nIndex] = false;
6103  /* do not increment nNullItems, since the user requested that the */
6104  /* neighborhood replace NODATA values with the central pixel value */
6105  }
6106  else {
6107  neighborData[nIndex] = PointerGetDatum(NULL);
6108  neighborNulls[nIndex] = true;
6109  nNullItems++;
6110  }
6111  }
6112  }
6113  else {
6114  /* Fill this will NULL if we can't read the raster pixel. */
6115  neighborData[nIndex] = PointerGetDatum(NULL);
6116  neighborNulls[nIndex] = true;
6117  nNullItems++;
6118  }
6119  /* Next neighbor position */
6120  nIndex++;
6121  }
6122  }
6123 
6128  if (!(nNodataOnly || /* neighborhood only contains NODATA -- OR -- */
6129  (nNullSkip && nNullItems > 0) || /* neighborhood should skip any NODATA cells, and a NODATA cell was detected -- OR -- */
6130  (valuereplace && nNullItems > 0))) { /* neighborhood should replace NODATA cells with the central pixel value, and a NODATA cell was detected */
6131  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: (%dx%d), %dx%d neighborhood",
6132  x, y, winwidth, winheight);
6133 
6134  neighborDatum = construct_md_array((void *)neighborData, neighborNulls, 2, neighborDims, neighborLbs,
6135  FLOAT8OID, typlen, typbyval, typalign);
6136 
6137 #if POSTGIS_PGSQL_VERSION < 120
6138  /* Assign the neighbor matrix as the first argument to the user function */
6139  cbdata.arg[0] = PointerGetDatum(neighborDatum);
6140 
6141  /* Invoke the user function */
6142  tmpnewval = FunctionCallInvoke(&cbdata);
6143 
6144  /* Get the return value of the user function */
6145  if (cbdata.isnull) {
6146  newval = newnodatavalue;
6147  }
6148 #else
6149  /* Assign the neighbor matrix as the first argument to the user function */
6150  cbdata->args[0].value = PointerGetDatum(neighborDatum);
6151 
6152  /* Invoke the user function */
6153  tmpnewval = FunctionCallInvoke(cbdata);
6154 
6155  /* Get the return value of the user function */
6156  if (cbdata->isnull)
6157  {
6158  newval = newnodatavalue;
6159  }
6160 #endif
6161  else {
6162  newval = DatumGetFloat8(tmpnewval);
6163  }
6164 
6165  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: new value = %f",
6166  newval);
6167 
6168  rt_band_set_pixel(newband, x, y, newval, NULL);
6169  }
6170 
6171  /* reset the number of null items in the neighborhood */
6172  nNullItems = 0;
6173  }
6174  }
6175 
6176 
6177  /* clean up */
6178  pfree(neighborNulls);
6179  pfree(neighborData);
6180  pfree(strFromText);
6181  pfree(txtCallbackParam);
6182 
6184  PG_FREE_IF_COPY(pgraster, 0);
6185 
6186  /* The newrast band has been modified */
6187 
6188  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: raster modified, serializing it.");
6189  /* Serialize created raster */
6190 
6191  pgrtn = rt_raster_serialize(newrast);
6192  rt_raster_destroy(newrast);
6193  if (NULL == pgrtn)
6194  PG_RETURN_NULL();
6195 
6196  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: raster serialized");
6197  POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFctNgb: returning raster");
6198 
6199  SET_VARSIZE(pgrtn, pgrtn->size);
6200  PG_RETURN_POINTER(pgrtn);
6201 }
6202 
6203 #define ARGKWCOUNT 8
6204 
6209 Datum RASTER_mapAlgebra2(PG_FUNCTION_ARGS)
6210 {
6211  const uint32_t set_count = 2;
6212  rt_pgraster *pgrast[2] = { NULL, NULL };
6213  int pgrastpos[2] = {-1, -1};
6214  rt_pgraster *pgrtn;
6215  rt_raster rast[2] = {NULL};
6216  int _isempty[2] = {0};
6217  uint32_t bandindex[2] = {0};
6218  rt_raster _rast[2] = {NULL};
6219  rt_band _band[2] = {NULL};
6220  int _hasnodata[2] = {0};
6221  double _nodataval[2] = {0};
6222  double _offset[4] = {0.};
6223  double _rastoffset[2][4] = {{0.}};
6224  int _haspixel[2] = {0};
6225  double _pixel[2] = {0};
6226  int _pos[2][2] = {{0}};
6227  uint16_t _dim[2][2] = {{0}};
6228 
6229  char *pixtypename = NULL;
6230  rt_pixtype pixtype = PT_END;
6231  char *extenttypename = NULL;
6232  rt_extenttype extenttype = ET_INTERSECTION;
6233 
6234  rt_raster raster = NULL;
6235  rt_band band = NULL;
6236  uint16_t dim[2] = {0};
6237  int haspixel = 0;
6238  double pixel = 0.;
6239  double nodataval = 0;
6240  double gt[6] = {0.};
6241 
6242  Oid calltype = InvalidOid;
6243 
6244  const uint32_t spi_count = 3;
6245  uint16_t spi_exprpos[3] = {4, 7, 8};
6246  uint32_t spi_argcount[3] = {0};
6247  char *expr = NULL;
6248  char *sql = NULL;
6249  SPIPlanPtr spi_plan[3] = {NULL};
6250  uint16_t spi_empty = 0;
6251  Oid *argtype = NULL;
6252  uint8_t argpos[3][8] = {{0}};
6253  char *argkw[] = {"[rast1.x]", "[rast1.y]", "[rast1.val]", "[rast1]", "[rast2.x]", "[rast2.y]", "[rast2.val]", "[rast2]"};
6254  Datum values[ARGKWCOUNT];
6255  char nulls[ARGKWCOUNT];
6256  TupleDesc tupdesc;
6257  SPITupleTable *tuptable = NULL;
6258  HeapTuple tuple;
6259  Datum datum;
6260  bool isnull = FALSE;
6261  int hasargval[3] = {0};
6262  double argval[3] = {0.};
6263  int hasnodatanodataval = 0;
6264  double nodatanodataval = 0;
6265  int isnodata = 0;
6266 
6267  Oid ufc_noid = InvalidOid;
6268  FmgrInfo ufl_info;
6269 #if POSTGIS_PGSQL_VERSION < 120
6270  FunctionCallInfoData ufc_info;
6271 #else
6272  LOCAL_FCINFO(ufc_info, FUNC_MAX_ARGS); /* Could be optimized */
6273 #endif
6274  int ufc_nullcount = 0;
6275 
6276  int idx = 0;
6277  uint32_t i = 0;
6278  uint32_t j = 0;
6279  uint32_t k = 0;
6280  uint32_t x = 0;
6281  uint32_t y = 0;
6282  int _x = 0;
6283  int _y = 0;
6284  int err;
6285  int aligned = 0;
6286  int len = 0;
6287 
6288  POSTGIS_RT_DEBUG(3, "Starting RASTER_mapAlgebra2");
6289 
6290  for (i = 0, j = 0; i < set_count; i++) {
6291  if (!PG_ARGISNULL(j)) {
6292  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
6293  pgrastpos[i] = j;
6294  j++;
6295 
6296  /* raster */
6297  rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
6298  if (!rast[i]) {
6299  for (k = 0; k <= i; k++) {
6300  if (k < i && rast[k] != NULL)
6301  rt_raster_destroy(rast[k]);
6302  if (pgrastpos[k] != -1)
6303  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6304  }
6305  elog(ERROR, "RASTER_mapAlgebra2: Could not deserialize the %s raster", i < 1 ? "first" : "second");
6306  PG_RETURN_NULL();
6307  }
6308 
6309  /* empty */
6310  _isempty[i] = rt_raster_is_empty(rast[i]);
6311 
6312  /* band index */
6313  if (!PG_ARGISNULL(j)) {
6314  bandindex[i] = PG_GETARG_INT32(j);
6315  }
6316  j++;
6317  }
6318  else {
6319  _isempty[i] = 1;
6320  j += 2;
6321  }
6322 
6323  POSTGIS_RT_DEBUGF(3, "_isempty[%d] = %d", i, _isempty[i]);
6324  }
6325 
6326  /* both rasters are NULL */
6327  if (rast[0] == NULL && rast[1] == NULL) {
6328  elog(NOTICE, "The two rasters provided are NULL. Returning NULL");
6329  for (k = 0; k < set_count; k++) {
6330  if (pgrastpos[k] != -1)
6331  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6332  }
6333  PG_RETURN_NULL();
6334  }
6335 
6336  /* both rasters are empty */
6337  if (_isempty[0] && _isempty[1]) {
6338  elog(NOTICE, "The two rasters provided are empty. Returning empty raster");
6339 
6340  raster = rt_raster_new(0, 0);
6341  if (raster == NULL) {
6342  for (k = 0; k < set_count; k++) {
6343  if (rast[k] != NULL)
6344  rt_raster_destroy(rast[k]);
6345  if (pgrastpos[k] != -1)
6346  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6347  }
6348  elog(ERROR, "RASTER_mapAlgebra2: Could not create empty raster");
6349  PG_RETURN_NULL();
6350  }
6351  rt_raster_set_scale(raster, 0, 0);
6352 
6353  pgrtn = rt_raster_serialize(raster);
6355  if (!pgrtn)
6356  PG_RETURN_NULL();
6357 
6358  SET_VARSIZE(pgrtn, pgrtn->size);
6359  PG_RETURN_POINTER(pgrtn);
6360  }
6361 
6362  /* replace the empty or NULL raster with one matching the other */
6363  if (
6364  (rast[0] == NULL || _isempty[0]) ||
6365  (rast[1] == NULL || _isempty[1])
6366  ) {
6367  /* first raster is empty */
6368  if (rast[0] == NULL || _isempty[0]) {
6369  i = 0;
6370  j = 1;
6371  }
6372  /* second raster is empty */
6373  else {
6374  i = 1;
6375  j = 0;
6376  }
6377 
6378  _rast[j] = rast[j];
6379 
6380  /* raster is empty, destroy it */
6381  if (_rast[i] != NULL)
6382  rt_raster_destroy(_rast[i]);
6383 
6384  _dim[i][0] = rt_raster_get_width(_rast[j]);
6385  _dim[i][1] = rt_raster_get_height(_rast[j]);
6386  _dim[j][0] = rt_raster_get_width(_rast[j]);
6387  _dim[j][1] = rt_raster_get_height(_rast[j]);
6388 
6389  _rast[i] = rt_raster_new(
6390  _dim[j][0],
6391  _dim[j][1]
6392  );
6393  if (_rast[i] == NULL) {
6394  rt_raster_destroy(_rast[j]);
6395  for (k = 0; k < set_count; k++) {
6396  if (pgrastpos[k] != -1)
6397  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6398  }
6399  elog(ERROR, "RASTER_mapAlgebra2: Could not create NODATA raster");
6400  PG_RETURN_NULL();
6401  }
6402  rt_raster_set_srid(_rast[i], rt_raster_get_srid(_rast[j]));
6403 
6406  }
6407  else {
6408  _rast[0] = rast[0];
6409  _dim[0][0] = rt_raster_get_width(_rast[0]);
6410  _dim[0][1] = rt_raster_get_height(_rast[0]);
6411 
6412  _rast[1] = rast[1];
6413  _dim[1][0] = rt_raster_get_width(_rast[1]);
6414  _dim[1][1] = rt_raster_get_height(_rast[1]);
6415  }
6416 
6417  /* SRID must match */
6418  /*
6419  if (rt_raster_get_srid(_rast[0]) != rt_raster_get_srid(_rast[1])) {
6420  elog(NOTICE, "The two rasters provided have different SRIDs. Returning NULL");
6421  for (k = 0; k < set_count; k++) {
6422  if (_rast[k] != NULL)
6423  rt_raster_destroy(_rast[k]);
6424  if (pgrastpos[k] != -1)
6425  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6426  }
6427  PG_RETURN_NULL();
6428  }
6429  */
6430 
6431  /* same alignment */
6432  if (rt_raster_same_alignment(_rast[0], _rast[1], &aligned, NULL) != ES_NONE) {
6433  for (k = 0; k < set_count; k++) {
6434  if (_rast[k] != NULL)
6435  rt_raster_destroy(_rast[k]);
6436  if (pgrastpos[k] != -1)
6437  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6438  }
6439  elog(ERROR, "RASTER_mapAlgebra2: Could not test for alignment on the two rasters");
6440  PG_RETURN_NULL();
6441  }
6442  if (!aligned) {
6443  elog(NOTICE, "The two rasters provided do not have the same alignment. Returning NULL");
6444  for (k = 0; k < set_count; k++) {
6445  if (_rast[k] != NULL)
6446  rt_raster_destroy(_rast[k]);
6447  if (pgrastpos[k] != -1)
6448  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6449  }
6450  PG_RETURN_NULL();
6451  }
6452 
6453  /* pixel type */
6454  if (!PG_ARGISNULL(5)) {
6455  pixtypename = text_to_cstring(PG_GETARG_TEXT_P(5));
6456  /* Get the pixel type index */
6457  pixtype = rt_pixtype_index_from_name(pixtypename);
6458  if (pixtype == PT_END ) {
6459  for (k = 0; k < set_count; k++) {
6460  if (_rast[k] != NULL)
6461  rt_raster_destroy(_rast[k]);
6462  if (pgrastpos[k] != -1)
6463  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6464  }
6465  elog(ERROR, "RASTER_mapAlgebra2: Invalid pixel type: %s", pixtypename);
6466  PG_RETURN_NULL();
6467  }
6468  }
6469 
6470  /* extent type */
6471  if (!PG_ARGISNULL(6)) {
6472  extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(6))));
6473  extenttype = rt_util_extent_type(extenttypename);
6474  }
6475  POSTGIS_RT_DEBUGF(3, "extenttype: %d %s", extenttype, extenttypename);
6476 
6477  /* computed raster from extent type */
6479  _rast[0], _rast[1],
6480  extenttype,
6481  &raster, _offset
6482  );
6483  if (err != ES_NONE) {
6484  for (k = 0; k < set_count; k++) {
6485  if (_rast[k] != NULL)
6486  rt_raster_destroy(_rast[k]);
6487  if (pgrastpos[k] != -1)
6488  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6489  }
6490  elog(ERROR, "RASTER_mapAlgebra2: Could not get output raster of correct extent");
6491  PG_RETURN_NULL();
6492  }
6493 
6494  /* copy offsets */
6495  _rastoffset[0][0] = _offset[0];
6496  _rastoffset[0][1] = _offset[1];
6497  _rastoffset[1][0] = _offset[2];
6498  _rastoffset[1][1] = _offset[3];
6499 
6500  /* get output raster dimensions */
6501  dim[0] = rt_raster_get_width(raster);
6502  dim[1] = rt_raster_get_height(raster);
6503 
6504  i = 2;
6505  /* handle special cases for extent */
6506  switch (extenttype) {
6507  case ET_FIRST:
6508  i = 0;
6509  /* fall through */
6510  case ET_SECOND:
6511  if (i > 1)
6512  i = 1;
6513 
6514  if (
6515  _isempty[i] && (
6516  (extenttype == ET_FIRST && i == 0) ||
6517  (extenttype == ET_SECOND && i == 1)
6518  )
6519  ) {
6520  elog(NOTICE, "The %s raster is NULL. Returning NULL", (i != 1 ? "FIRST" : "SECOND"));
6521  for (k = 0; k < set_count; k++) {
6522  if (_rast[k] != NULL)
6523  rt_raster_destroy(_rast[k]);
6524  if (pgrastpos[k] != -1)
6525  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6526  }
6528  PG_RETURN_NULL();
6529  }
6530 
6531  /* specified band not found */
6532  if (!rt_raster_has_band(_rast[i], bandindex[i] - 1)) {
6533  elog(NOTICE, "The %s raster does not have the band at index %d. Returning no band raster of correct extent",
6534  (i != 1 ? "FIRST" : "SECOND"), bandindex[i]
6535  );
6536 
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  }
6543 
6544  pgrtn = rt_raster_serialize(raster);
6546  if (!pgrtn) PG_RETURN_NULL();
6547 
6548  SET_VARSIZE(pgrtn, pgrtn->size);
6549  PG_RETURN_POINTER(pgrtn);
6550  }
6551  break;
6552  case ET_UNION:
6553  break;
6554  case ET_INTERSECTION:
6555  /* no intersection */
6556  if (
6557  _isempty[0] || _isempty[1] ||
6558  !dim[0] || !dim[1]
6559  ) {
6560  elog(NOTICE, "The two rasters provided have no intersection. Returning no band raster");
6561 
6562  /* raster has dimension, replace with no band raster */
6563  if (dim[0] || dim[1]) {
6565 
6566  raster = rt_raster_new(0, 0);
6567  if (raster == NULL) {
6568  for (k = 0; k < set_count; k++) {
6569  if (_rast[k] != NULL)
6570  rt_raster_destroy(_rast[k]);
6571  if (pgrastpos[k] != -1)
6572  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6573  }
6574  elog(ERROR, "RASTER_mapAlgebra2: Could not create no band raster");
6575  PG_RETURN_NULL();
6576  }
6577 
6578  rt_raster_set_scale(raster, 0, 0);
6580  }
6581 
6582  for (k = 0; k < set_count; k++) {
6583  if (_rast[k] != NULL)
6584  rt_raster_destroy(_rast[k]);
6585  if (pgrastpos[k] != -1)
6586  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6587  }
6588 
6589  pgrtn = rt_raster_serialize(raster);
6591  if (!pgrtn) PG_RETURN_NULL();
6592 
6593  SET_VARSIZE(pgrtn, pgrtn->size);
6594  PG_RETURN_POINTER(pgrtn);
6595  }
6596  break;
6597  case ET_LAST:
6598  case ET_CUSTOM:
6599  for (k = 0; k < set_count; k++) {
6600  if (_rast[k] != NULL)
6601  rt_raster_destroy(_rast[k]);
6602  if (pgrastpos[k] != -1)
6603  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6604  }
6605  elog(ERROR, "RASTER_mapAlgebra2: ET_LAST and ET_CUSTOM are not implemented");
6606  PG_RETURN_NULL();
6607  break;
6608  }
6609 
6610  /* both rasters do not have specified bands */
6611  if (
6612  (!_isempty[0] && !rt_raster_has_band(_rast[0], bandindex[0] - 1)) &&
6613  (!_isempty[1] && !rt_raster_has_band(_rast[1], bandindex[1] - 1))
6614  ) {
6615  elog(NOTICE, "The two rasters provided do not have the respectively specified band indices. Returning no band raster of correct extent");
6616 
6617  for (k = 0; k < set_count; k++) {
6618  if (_rast[k] != NULL)
6619  rt_raster_destroy(_rast[k]);
6620  if (pgrastpos[k] != -1)
6621  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6622  }
6623 
6624  pgrtn = rt_raster_serialize(raster);
6626  if (!pgrtn) PG_RETURN_NULL();
6627 
6628  SET_VARSIZE(pgrtn, pgrtn->size);
6629  PG_RETURN_POINTER(pgrtn);
6630  }
6631 
6632  /* get bands */
6633  for (i = 0; i < set_count; i++) {
6634  if (_isempty[i] || !rt_raster_has_band(_rast[i], bandindex[i] - 1)) {
6635  _hasnodata[i] = 1;
6636  _nodataval[i] = 0;
6637 
6638  continue;
6639  }
6640 
6641  _band[i] = rt_raster_get_band(_rast[i], bandindex[i] - 1);
6642  if (_band[i] == NULL) {
6643  for (k = 0; k < set_count; k++) {
6644  if (_rast[k] != NULL)
6645  rt_raster_destroy(_rast[k]);
6646  if (pgrastpos[k] != -1)
6647  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6648  }
6650  elog(ERROR, "RASTER_mapAlgebra2: Could not get band %d of the %s raster",
6651  bandindex[i],
6652  (i < 1 ? "FIRST" : "SECOND")
6653  );
6654  PG_RETURN_NULL();
6655  }
6656 
6657  _hasnodata[i] = rt_band_get_hasnodata_flag(_band[i]);
6658  if (_hasnodata[i])
6659  rt_band_get_nodata(_band[i], &(_nodataval[i]));
6660  }
6661 
6662  /* pixtype is PT_END, get pixtype based upon extent */
6663  if (pixtype == PT_END) {
6664  if ((extenttype == ET_SECOND && !_isempty[1]) || _isempty[0])
6665  pixtype = rt_band_get_pixtype(_band[1]);
6666  else
6667  pixtype = rt_band_get_pixtype(_band[0]);
6668  }
6669 
6670  /* nodata value for new band */
6671  if (extenttype == ET_SECOND && !_isempty[1] && _hasnodata[1]) {
6672  nodataval = _nodataval[1];
6673  }
6674  else if (!_isempty[0] && _hasnodata[0]) {
6675  nodataval = _nodataval[0];
6676  }
6677  else if (!_isempty[1] && _hasnodata[1]) {
6678  nodataval = _nodataval[1];
6679  }
6680  else {
6681  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));
6682  nodataval = rt_pixtype_get_min_value(pixtype);
6683  }
6684 
6685  /* add band to output raster */
6687  raster,
6688  pixtype,
6689  nodataval,
6690  1, nodataval,
6691  0
6692  ) < 0) {
6693  for (k = 0; k < set_count; k++) {
6694  if (_rast[k] != NULL)
6695  rt_raster_destroy(_rast[k]);
6696  if (pgrastpos[k] != -1)
6697  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6698  }
6700  elog(ERROR, "RASTER_mapAlgebra2: Could not add new band to output raster");
6701  PG_RETURN_NULL();
6702  }
6703 
6704  /* get output band */
6706  if (band == NULL) {
6707  for (k = 0; k < set_count; k++) {
6708  if (_rast[k] != NULL)
6709  rt_raster_destroy(_rast[k]);
6710  if (pgrastpos[k] != -1)
6711  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6712  }
6714  elog(ERROR, "RASTER_mapAlgebra2: Could not get newly added band of output raster");
6715  PG_RETURN_NULL();
6716  }
6717 
6718  POSTGIS_RT_DEBUGF(4, "offsets = (%d, %d, %d, %d)",
6719  (int) _rastoffset[0][0],
6720  (int) _rastoffset[0][1],
6721  (int) _rastoffset[1][0],
6722  (int) _rastoffset[1][1]
6723  );
6724 
6725  POSTGIS_RT_DEBUGF(4, "metadata = (%f, %f, %d, %d, %f, %f, %f, %f, %d)",
6735  );
6736 
6737  /*
6738  determine who called this function
6739  Arg 4 will either be text or regprocedure
6740  */
6741  POSTGIS_RT_DEBUG(3, "checking parameter type for arg 4");
6742  calltype = get_fn_expr_argtype(fcinfo->flinfo, 4);
6743 
6744  switch(calltype) {
6745  case TEXTOID: {
6746  POSTGIS_RT_DEBUG(3, "arg 4 is \"expression\"!");
6747 
6748  /* connect SPI */
6749  if (SPI_connect() != SPI_OK_CONNECT) {
6750  for (k = 0; k < set_count; k++) {
6751  if (_rast[k] != NULL)
6752  rt_raster_destroy(_rast[k]);
6753  if (pgrastpos[k] != -1)
6754  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6755  }
6757  elog(ERROR, "RASTER_mapAlgebra2: Could not connect to the SPI manager");
6758  PG_RETURN_NULL();
6759  }
6760 
6761  /* reset hasargval */
6762  memset(hasargval, 0, sizeof(int) * spi_count);
6763 
6764  /*
6765  process expressions
6766 
6767  spi_exprpos elements are:
6768  4 - expression => spi_plan[0]
6769  7 - nodata1expr => spi_plan[1]
6770  8 - nodata2expr => spi_plan[2]
6771  */
6772  for (i = 0; i < spi_count; i++) {
6773  if (!PG_ARGISNULL(spi_exprpos[i])) {
6774  char *tmp = NULL;
6775  char place[5] = "$1";
6776  expr = text_to_cstring(PG_GETARG_TEXT_P(spi_exprpos[i]));
6777  POSTGIS_RT_DEBUGF(3, "raw expr #%d: %s", i, expr);
6778 
6779  for (j = 0, k = 1; j < ARGKWCOUNT; j++) {
6780  /* attempt to replace keyword with placeholder */
6781  len = 0;
6782  tmp = rtpg_strreplace(expr, argkw[j], place, &len);
6783  pfree(expr);
6784  expr = tmp;
6785 
6786  if (len) {
6787  spi_argcount[i]++;
6788  argpos[i][j] = k++;
6789 
6790  sprintf(place, "$%d", k);
6791  }
6792  else
6793  argpos[i][j] = 0;
6794  }
6795 
6796  len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision");
6797  sql = (char *) palloc(len + 1);
6798  if (sql == NULL) {
6799 
6800  for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6801  SPI_finish();
6802 
6803  for (k = 0; k < set_count; k++) {
6804  if (_rast[k] != NULL)
6805  rt_raster_destroy(_rast[k]);
6806  if (pgrastpos[k] != -1)
6807  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6808  }
6810 
6811  elog(ERROR, "RASTER_mapAlgebra2: Could not allocate memory for expression parameter %d", spi_exprpos[i]);
6812  PG_RETURN_NULL();
6813  }
6814 
6815  memcpy(sql, "SELECT (", strlen("SELECT ("));
6816  memcpy(sql + strlen("SELECT ("), expr, strlen(expr));
6817  memcpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision"));
6818  sql[len] = '\0';
6819 
6820  POSTGIS_RT_DEBUGF(3, "sql #%d: %s", i, sql);
6821 
6822  /* create prepared plan */
6823  if (spi_argcount[i]) {
6824  argtype = (Oid *) palloc(spi_argcount[i] * sizeof(Oid));
6825  if (argtype == NULL) {
6826 
6827  pfree(sql);
6828  for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6829  SPI_finish();
6830 
6831  for (k = 0; k < set_count; k++) {
6832  if (_rast[k] != NULL)
6833  rt_raster_destroy(_rast[k]);
6834  if (pgrastpos[k] != -1)
6835  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6836  }
6838 
6839  elog(ERROR, "RASTER_mapAlgebra2: Could not allocate memory for prepared plan argtypes of expression parameter %d", spi_exprpos[i]);
6840  PG_RETURN_NULL();
6841  }
6842 
6843  /* specify datatypes of parameters */
6844  for (j = 0, k = 0; j < ARGKWCOUNT; j++) {
6845  if (argpos[i][j] < 1) continue;
6846 
6847  /* positions are INT4 */
6848  if (
6849  (strstr(argkw[j], "[rast1.x]") != NULL) ||
6850  (strstr(argkw[j], "[rast1.y]") != NULL) ||
6851  (strstr(argkw[j], "[rast2.x]") != NULL) ||
6852  (strstr(argkw[j], "[rast2.y]") != NULL)
6853  ) {
6854  argtype[k] = INT4OID;
6855  }
6856  /* everything else is FLOAT8 */
6857  else {
6858  argtype[k] = FLOAT8OID;
6859  }
6860 
6861  k++;
6862  }
6863 
6864  spi_plan[i] = SPI_prepare(sql, spi_argcount[i], argtype);
6865  pfree(argtype);
6866 
6867  if (spi_plan[i] == NULL) {
6868 
6869  pfree(sql);
6870  for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6871  SPI_finish();
6872 
6873  for (k = 0; k < set_count; k++) {
6874  if (_rast[k] != NULL)
6875  rt_raster_destroy(_rast[k]);
6876  if (pgrastpos[k] != -1)
6877  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6878  }
6880 
6881  elog(ERROR, "RASTER_mapAlgebra2: Could not create prepared plan of expression parameter %d", spi_exprpos[i]);
6882  PG_RETURN_NULL();
6883  }
6884  }
6885  /* no args, just execute query */
6886  else {
6887  err = SPI_execute(sql, TRUE, 0);
6888  if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
6889 
6890  pfree(sql);
6891  for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6892  SPI_finish();
6893 
6894  for (k = 0; k < set_count; k++) {
6895  if (_rast[k] != NULL)
6896  rt_raster_destroy(_rast[k]);
6897  if (pgrastpos[k] != -1)
6898  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6899  }
6901 
6902  elog(ERROR, "RASTER_mapAlgebra2: Could not evaluate expression parameter %d", spi_exprpos[i]);
6903  PG_RETURN_NULL();
6904  }
6905 
6906  /* get output of prepared plan */
6907  tupdesc = SPI_tuptable->tupdesc;
6908  tuptable = SPI_tuptable;
6909  tuple = tuptable->vals[0];
6910 
6911  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
6912  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
6913 
6914  pfree(sql);
6915  if (SPI_tuptable) SPI_freetuptable(tuptable);
6916  for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6917  SPI_finish();
6918 
6919  for (k = 0; k < set_count; k++) {
6920  if (_rast[k] != NULL)
6921  rt_raster_destroy(_rast[k]);
6922  if (pgrastpos[k] != -1)
6923  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6924  }
6926 
6927  elog(ERROR, "RASTER_mapAlgebra2: Could not get result of expression parameter %d", spi_exprpos[i]);
6928  PG_RETURN_NULL();
6929  }
6930 
6931  if (!isnull) {
6932  hasargval[i] = 1;
6933  argval[i] = DatumGetFloat8(datum);
6934  }
6935 
6936  if (SPI_tuptable) SPI_freetuptable(tuptable);
6937  }
6938 
6939  pfree(sql);
6940  }
6941  else
6942  spi_empty++;
6943  }
6944 
6945  /* nodatanodataval */
6946  if (!PG_ARGISNULL(9)) {
6947  hasnodatanodataval = 1;
6948  nodatanodataval = PG_GETARG_FLOAT8(9);
6949  }
6950  else
6951  hasnodatanodataval = 0;
6952  break;
6953  }
6954  case REGPROCEDUREOID: {
6955  POSTGIS_RT_DEBUG(3, "arg 4 is \"userfunction\"!");
6956  if (!PG_ARGISNULL(4)) {
6957 
6958  ufc_nullcount = 0;
6959  ufc_noid = PG_GETARG_OID(4);
6960 
6961  /* get function info */
6962  fmgr_info(ufc_noid, &ufl_info);
6963 
6964  /* function cannot return set */
6965  err = 0;
6966  if (ufl_info.fn_retset) {
6967  err = 1;
6968  }
6969  /* function should have correct # of args */
6970  else if (ufl_info.fn_nargs < 3 || ufl_info.fn_nargs > 4) {
6971  err = 2;
6972  }
6973 
6974  /*
6975  TODO: consider adding checks of the userfunction parameters
6976  should be able to use get_fn_expr_argtype() of fmgr.c
6977  */
6978 
6979  if (err > 0) {
6980  for (k = 0; k < set_count; k++) {
6981  if (_rast[k] != NULL)
6982  rt_raster_destroy(_rast[k]);
6983  if (pgrastpos[k] != -1)
6984  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6985  }
6987 
6988  if (err > 1)
6989  elog(ERROR, "RASTER_mapAlgebra2: Function provided must have three or four input parameters");
6990  else
6991  elog(ERROR, "RASTER_mapAlgebra2: Function provided must return double precision not resultset");
6992  PG_RETURN_NULL();
6993  }
6994 
6995  if (func_volatile(ufc_noid) == 'v') {
6996  elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
6997  }
6998 
6999  /* prep function call data */
7000 #if POSTGIS_PGSQL_VERSION < 120
7001  InitFunctionCallInfoData(ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7002  memset(ufc_info.argnull, FALSE, sizeof(bool) * ufl_info.fn_nargs);
7003 #else
7004  InitFunctionCallInfoData(
7005  *ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7006  ufc_info->args[0].isnull = FALSE;
7007  ufc_info->args[1].isnull = FALSE;
7008  ufc_info->args[2].isnull = FALSE;
7009  if (ufl_info.fn_nargs == 4)
7010  ufc_info->args[3].isnull = FALSE;
7011 #endif
7012 
7013  if (ufl_info.fn_nargs != 4)
7014  k = 2;
7015  else
7016  k = 3;
7017 #if POSTGIS_PGSQL_VERSION < 120
7018  if (!PG_ARGISNULL(7)) {
7019  ufc_info.arg[k] = PG_GETARG_DATUM(7);
7020  }
7021  else {
7022  ufc_info.arg[k] = (Datum) NULL;
7023  ufc_info.argnull[k] = TRUE;
7024  ufc_nullcount++;
7025  }
7026 #else
7027  if (!PG_ARGISNULL(7))
7028  {
7029  ufc_info->args[k].value = PG_GETARG_DATUM(7);
7030  }
7031  else
7032  {
7033  ufc_info->args[k].value = (Datum)NULL;
7034  ufc_info->args[k].isnull = TRUE;
7035  ufc_nullcount++;
7036  }
7037 #endif
7038  }
7039  break;
7040  }
7041  default:
7042  for (k = 0; k < set_count; k++) {
7043  if (_rast[k] != NULL)
7044  rt_raster_destroy(_rast[k]);
7045  if (pgrastpos[k] != -1)
7046  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7047  }
7049  elog(ERROR, "RASTER_mapAlgebra2: Invalid data type for expression or userfunction");
7050  PG_RETURN_NULL();
7051  break;
7052  }
7053 
7054  /* loop over pixels */
7055  /* if any expression present, run */
7056  if ((
7057  (calltype == TEXTOID) && (
7058  (spi_empty != spi_count) || hasnodatanodataval
7059  )
7060  ) || (
7061  (calltype == REGPROCEDUREOID) && (ufc_noid != InvalidOid)
7062  )) {
7063  for (x = 0; x < dim[0]; x++) {
7064  for (y = 0; y < dim[1]; y++) {
7065 
7066  /* get pixel from each raster */
7067  for (i = 0; i < set_count; i++) {
7068  _haspixel[i] = 0;
7069  _pixel[i] = 0;
7070 
7071  /* row/column */
7072  _x = (int)x - (int)_rastoffset[i][0];
7073  _y = (int)y - (int)_rastoffset[i][1];
7074 
7075  /* store _x and _y in 1-based */
7076  _pos[i][0] = _x + 1;
7077  _pos[i][1] = _y + 1;
7078 
7079  /* get pixel value */
7080  if (_band[i] == NULL) {
7081  if (!_hasnodata[i]) {
7082  _haspixel[i] = 1;
7083  _pixel[i] = _nodataval[i];
7084  }
7085  }
7086  else if (
7087  !_isempty[i] &&
7088  (_x >= 0 && _x < _dim[i][0]) &&
7089  (_y >= 0 && _y < _dim[i][1])
7090  ) {
7091  err = rt_band_get_pixel(_band[i], _x, _y, &(_pixel[i]), &isnodata);
7092  if (err != ES_NONE) {
7093 
7094  if (calltype == TEXTOID) {
7095  for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7096  SPI_finish();
7097  }
7098 
7099  for (k = 0; k < set_count; k++) {
7100  if (_rast[k] != NULL)
7101  rt_raster_destroy(_rast[k]);
7102  if (pgrastpos[k] != -1)
7103  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7104  }
7106 
7107  elog(ERROR, "RASTER_mapAlgebra2: Could not get pixel of %s raster", (i < 1 ? "FIRST" : "SECOND"));
7108  PG_RETURN_NULL();
7109  }
7110 
7111  if (!_hasnodata[i] || !isnodata)
7112  _haspixel[i] = 1;
7113  }
7114 
7115  POSTGIS_RT_DEBUGF(5, "pixel r%d(%d, %d) = %d, %f",
7116  i,
7117  _x, _y,
7118  _haspixel[i],
7119  _pixel[i]
7120  );
7121  }
7122 
7123  haspixel = 0;
7124 
7125  switch (calltype) {
7126  case TEXTOID: {
7127  /* which prepared plan to use? */
7128  /* !pixel0 && !pixel1 */
7129  /* use nodatanodataval */
7130  if (!_haspixel[0] && !_haspixel[1])
7131  i = 3;
7132  /* pixel0 && !pixel1 */
7133  /* run spi_plan[2] (nodata2expr) */
7134  else if (_haspixel[0] && !_haspixel[1])
7135  i = 2;
7136  /* !pixel0 && pixel1 */
7137  /* run spi_plan[1] (nodata1expr) */
7138  else if (!_haspixel[0] && _haspixel[1])
7139  i = 1;
7140  /* pixel0 && pixel1 */
7141  /* run spi_plan[0] (expression) */
7142  else
7143  i = 0;
7144 
7145  /* process values */
7146  if (i == 3) {
7147  if (hasnodatanodataval) {
7148  haspixel = 1;
7149  pixel = nodatanodataval;
7150  }
7151  }
7152  /* has an evaluated value */
7153  else if (hasargval[i]) {
7154  haspixel = 1;
7155  pixel = argval[i];
7156  }
7157  /* prepared plan exists */
7158  else if (spi_plan[i] != NULL) {
7159  POSTGIS_RT_DEBUGF(4, "Using prepared plan: %d", i);
7160 
7161  /* reset values to (Datum) NULL */
7162  memset(values, (Datum) NULL, sizeof(Datum) * ARGKWCOUNT);
7163  /* reset nulls to FALSE */
7164  memset(nulls, FALSE, sizeof(char) * ARGKWCOUNT);
7165 
7166  /* expression has argument(s) */
7167  if (spi_argcount[i]) {
7168  /* set values and nulls */
7169  for (j = 0; j < ARGKWCOUNT; j++) {
7170  idx = argpos[i][j];
7171  if (idx < 1) continue;
7172  idx--; /* 1-based becomes 0-based */
7173 
7174  if (strstr(argkw[j], "[rast1.x]") != NULL) {
7175  values[idx] = _pos[0][0];
7176  }
7177  else if (strstr(argkw[j], "[rast1.y]") != NULL) {
7178  values[idx] = _pos[0][1];
7179  }
7180  else if (
7181  (strstr(argkw[j], "[rast1.val]") != NULL) ||
7182  (strstr(argkw[j], "[rast1]") != NULL)
7183  ) {
7184  if (_isempty[0] || !_haspixel[0])
7185  nulls[idx] = TRUE;
7186  else
7187  values[idx] = Float8GetDatum(_pixel[0]);
7188  }
7189  else if (strstr(argkw[j], "[rast2.x]") != NULL) {
7190  values[idx] = _pos[1][0];
7191  }
7192  else if (strstr(argkw[j], "[rast2.y]") != NULL) {
7193  values[idx] = _pos[1][1];
7194  }
7195  else if (
7196  (strstr(argkw[j], "[rast2.val]") != NULL) ||
7197  (strstr(argkw[j], "[rast2]") != NULL)
7198  ) {
7199  if (_isempty[1] || !_haspixel[1])
7200  nulls[idx] = TRUE;
7201  else
7202  values[idx] = Float8GetDatum(_pixel[1]);
7203  }
7204  }
7205  }
7206 
7207  /* run prepared plan */
7208  err = SPI_execute_plan(spi_plan[i], values, nulls, TRUE, 1);
7209  if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
7210 
7211  for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7212  SPI_finish();
7213 
7214  for (k = 0; k < set_count; k++) {
7215  if (_rast[k] != NULL)
7216  rt_raster_destroy(_rast[k]);
7217  if (pgrastpos[k] != -1)
7218  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7219  }
7221 
7222  elog(ERROR, "RASTER_mapAlgebra2: Unexpected error when running prepared statement %d", i);
7223  PG_RETURN_NULL();
7224  }
7225 
7226  /* get output of prepared plan */
7227  tupdesc = SPI_tuptable->tupdesc;
7228  tuptable = SPI_tuptable;
7229  tuple = tuptable->vals[0];
7230 
7231  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
7232  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
7233 
7234  if (SPI_tuptable) SPI_freetuptable(tuptable);
7235  for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7236  SPI_finish();
7237 
7238  for (k = 0; k < set_count; k++) {
7239  if (_rast[k] != NULL)
7240  rt_raster_destroy(_rast[k]);
7241  if (pgrastpos[k] != -1)
7242  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7243  }
7245 
7246  elog(ERROR, "RASTER_mapAlgebra2: Could not get result of prepared statement %d", i);
7247  PG_RETURN_NULL();
7248  }
7249 
7250  if (!isnull) {
7251  haspixel = 1;
7252  pixel = DatumGetFloat8(datum);
7253  }
7254 
7255  if (SPI_tuptable) SPI_freetuptable(tuptable);
7256  }
7257  } break;
7258  case REGPROCEDUREOID: {
7259  Datum d[4];
7260  ArrayType *a;
7261 
7262  /* build fcnarg */
7263  for (i = 0; i < set_count; i++) {
7264 #if POSTGIS_PGSQL_VERSION < 120
7265  ufc_info.arg[i] = Float8GetDatum(_pixel[i]);
7266 #else
7267  ufc_info->args[i].value = Float8GetDatum(_pixel[i]);
7268 #endif
7269 
7270  if (_haspixel[i]) {
7271 #if POSTGIS_PGSQL_VERSION < 120
7272  ufc_info.argnull[i] = FALSE;
7273 #else
7274  ufc_info->args[i].isnull = FALSE;
7275 #endif
7276  ufc_nullcount--;
7277  }
7278  else {
7279 #if POSTGIS_PGSQL_VERSION < 120
7280  ufc_info.argnull[i] = TRUE;
7281 #else
7282  ufc_info->args[i].isnull = TRUE;
7283 #endif
7284  ufc_nullcount++;
7285  }
7286  }
7287 
7288  /* function is strict and null parameter is passed */
7289  /* http://archives.postgresql.org/pgsql-general/2011-11/msg00424.php */
7290  if (ufl_info.fn_strict && ufc_nullcount)
7291  break;
7292 
7293  /* 4 parameters, add position */
7294  if (ufl_info.fn_nargs == 4) {
7295  /* Datum of 4 element array */
7296  /* array is (x1, y1, x2, y2) */
7297  for (i = 0; i < set_count; i++) {
7298  if (i < 1) {
7299  d[0] = Int32GetDatum(_pos[i][0]);
7300  d[1] = Int32GetDatum(_pos[i][1]);
7301  }
7302  else {
7303  d[2] = Int32GetDatum(_pos[i][0]);
7304  d[3] = Int32GetDatum(_pos[i][1]);
7305  }
7306  }
7307 
7308  a = construct_array(d, 4, INT4OID, sizeof(int32), true, 'i');
7309 #if POSTGIS_PGSQL_VERSION < 120
7310  ufc_info.arg[2] = PointerGetDatum(a);
7311  ufc_info.argnull[2] = FALSE;
7312 #else
7313  ufc_info->args[2].value = PointerGetDatum(a);
7314  ufc_info->args[2].isnull = FALSE;
7315 #endif
7316  }
7317 
7318 #if POSTGIS_PGSQL_VERSION < 120
7319  datum = FunctionCallInvoke(&ufc_info);
7320 
7321  /* result is not null*/
7322  if (!ufc_info.isnull) {
7323  haspixel = 1;
7324  pixel = DatumGetFloat8(datum);
7325  }
7326 #else
7327  datum = FunctionCallInvoke(ufc_info);
7328 
7329  /* result is not null*/
7330  if (!ufc_info->isnull)
7331  {
7332  haspixel = 1;
7333  pixel = DatumGetFloat8(datum);
7334  }
7335 #endif
7336  } break;
7337  }
7338 
7339  /* burn pixel if haspixel != 0 */
7340  if (haspixel) {
7341  if (rt_band_set_pixel(band, x, y, pixel, NULL) != ES_NONE) {
7342 
7343  if (calltype == TEXTOID) {
7344  for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7345  SPI_finish();
7346  }
7347 
7348  for (k = 0; k < set_count; k++) {
7349  if (_rast[k] != NULL)
7350  rt_raster_destroy(_rast[k]);
7351  if (pgrastpos[k] != -1)
7352  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7353  }
7355 
7356  elog(ERROR, "RASTER_mapAlgebra2: Could not set pixel value of output raster");
7357  PG_RETURN_NULL();
7358  }
7359  }
7360 
7361  POSTGIS_RT_DEBUGF(5, "(x, y, val) = (%d, %d, %f)", x, y, haspixel ? pixel : nodataval);
7362 
7363  } /* y: height */
7364  } /* x: width */
7365  }
7366 
7367  /* CLEANUP */
7368  if (calltype == TEXTOID) {
7369  for (i = 0; i < spi_count; i++) {
7370  if (spi_plan[i] != NULL) SPI_freeplan(spi_plan[i]);
7371  }
7372  SPI_finish();
7373  }
7374 
7375  for (k = 0; k < set_count; k++) {
7376  if (_rast[k] != NULL)
7377  rt_raster_destroy(_rast[k]);
7378  if (pgrastpos[k] != -1)
7379  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7380  }
7381 
7382  pgrtn = rt_raster_serialize(raster);
7384  if (!pgrtn) PG_RETURN_NULL();
7385 
7386  POSTGIS_RT_DEBUG(3, "Finished RASTER_mapAlgebra2");
7387 
7388  SET_VARSIZE(pgrtn, pgrtn->size);
7389  PG_RETURN_POINTER(pgrtn);
7390 }
char * r
Definition: cu_in_wkt.c:24
#define TRUE
Definition: dbfopen.c:73
#define FALSE
Definition: dbfopen.c:72
int32_t gserialized_get_srid(const GSERIALIZED *g)
Extract the SRID from the serialized form (it is packed into three bytes so this is a handy function)...
Definition: gserialized.c:126
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
Definition: gserialized.c:239
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition: lwgeom.c:955
#define LWVARHDRSZ
Definition: liblwgeom.h:326
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1155
LWGEOM * lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
#define LWSIZE_GET(varsize)
Macro for reading the size from the GSERIALIZED size attribute.
Definition: liblwgeom.h:339
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:230
lwvarlena_t * lwgeom_to_wkb_varlena(const LWGEOM *geom, uint8_t variant)
Definition: lwout_wkb.c:851
#define WKB_SFSQL
Definition: liblwgeom.h:2158
int32_t clamp_srid(int32_t srid)
Return a valid SRID from an arbitrary integer Raises a notice if what comes out is different from wha...
Definition: lwutil.c:333
LWGEOM * lwgeom_force_2d(const LWGEOM *geom)
Strip out the Z/M components of an LWGEOM.
Definition: lwgeom.c:793
rt_band rt_band_reclass(rt_band srcband, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, rt_reclassexpr *exprset, int exprcount)
Returns new band with values reclassified.
Definition: rt_mapalgebra.c:50
#define FLT_NEQ(x, y)
Definition: librtcore.h:2379
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:302
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:360
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition: rt_raster.c:185
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
Definition: rt_raster.c:217
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
Definition: rt_raster.c:489
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster's geotransform using 6-element array.
Definition: rt_raster.c:731
rt_raster rt_raster_colormap(rt_raster raster, int nband, rt_colormap colormap)
Returns a new raster with up to four 8BUI bands (RGBA) from applying a colormap to the user-specified...
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_raster.c:141
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c: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:1376
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:86
rt_pixtype
Definition: librtcore.h:187
@ PT_32BUI
Definition: librtcore.h:196
@ PT_END
Definition: librtcore.h:199
@ PT_64BF
Definition: librtcore.h:198
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition: rt_raster.c:172
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:714
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:52
rt_extenttype rt_util_extent_type(const char *name)
Definition: rt_util.c:194
double rt_pixtype_get_min_value(rt_pixtype pixtype)
Return minimum value possible for pixel type.
Definition: rt_pixel.c:150
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
Definition: rt_raster.c:1375
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:2380
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:154
rt_band rt_raster_replace_band(rt_raster raster, rt_band band, int index)
Replace band at provided index with new band.
Definition: rt_raster.c:1526
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:1902
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:182
rt_errorstate rt_band_set_pixel_line(rt_band band, int x, int y, void *vals, uint32_t len)
Set values of multiple pixels.
Definition: rt_band.c: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:2618
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:376
rt_raster rt_raster_clone(rt_raster raster, uint8_t deep)
Clone an existing raster.
Definition: rt_raster.c:1568
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition: rt_geometry.c:838
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:133
rt_errorstate rt_raster_iterator(rt_iterator itrset, uint16_t itrcount, rt_extenttype extenttype, rt_raster customextent, rt_pixtype pixtype, uint8_t hasnodata, double nodataval, uint16_t distancex, uint16_t distancey, rt_mask mask, void *userarg, int(*callback)(rt_iterator_arg arg, void *userarg, double *value, int *nodata), rt_raster *rtnraster)
n-raster iterator.
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_raster.c:367
rt_extenttype
Definition: librtcore.h:202
@ ET_CUSTOM
Definition: librtcore.h:208
@ ET_LAST
Definition: librtcore.h:207
@ ET_INTERSECTION
Definition: librtcore.h:203
@ ET_UNION
Definition: librtcore.h:204
@ ET_SECOND
Definition: librtcore.h:206
@ ET_FIRST
Definition: librtcore.h:205
rt_raster rt_raster_from_band(rt_raster raster, uint32_t *bandNums, int count)
Construct a new rt_raster from an existing rt_raster and an array of band numbers.
Definition: rt_raster.c:1463
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1887
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:125
rt_errorstate rt_raster_from_two_rasters(rt_raster rast1, rt_raster rast2, rt_extenttype extenttype, rt_raster *rtnraster, double *offset)
Definition: rt_raster.c:3464
int rt_util_same_geotransform_matrix(double *gt1, double *gt2)
Definition: rt_util.c:535
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:1398
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:163
rt_errorstate rt_raster_same_alignment(rt_raster rast1, rt_raster rast2, int *aligned, char **reason)
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
Definition: rt_raster.c:194
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:710
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_raster.c:203
rt_bandstats rt_band_get_summary_stats(rt_band band, int exclude_nodata_value, double sample, int inc_vals, uint64_t *cK, double *cM, double *cQ)
Compute summary statistics for a band.
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:725
rt_errorstate rt_band_get_pixel_line(rt_band band, int x, int y, uint16_t len, void **vals, uint16_t *nvals)
Get values of multiple pixels.
Definition: rt_band.c:1137
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
Definition: rt_raster.c:226
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1362
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:385
if(!(yy_init))
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition: lwinline.h:203
int value
Definition: genraster.py:62
band
Definition: ovdump.py:58
pixel
Definition: pixval.py:91
nband
Definition: pixval.py:53
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
gt
Definition: window.py:78
char * rtpg_removespaces(char *str)
char * rtpg_strtoupper(char *str)
char * rtpg_chartrim(const char *input, char *remove)
char * rtpg_strrstr(const char *s1, const char *s2)
char * rtpg_trim(const char *input)
char ** rtpg_strsplit(const char *str, const char *delimiter, uint32_t *n)
char * rtpg_strreplace(const char *str, const char *oldstr, const char *newstr, int *count)
Definition: rtpg_internal.c:55
Datum RASTER_nMapAlgebra(PG_FUNCTION_ARGS)
#define ARGKWCOUNT
struct rtpg_colormap_arg_t * rtpg_colormap_arg
static rtpg_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)
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:69
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:73
unsigned int int32
Definition: shpopen.c:54
uint32_t size
Definition: liblwgeom.h:322
char data[]
Definition: liblwgeom.h:323
uint8_t color[4]
Definition: librtcore.h:2630
double value
Definition: librtcore.h:2629
int isnodata
Definition: librtcore.h:2628
Definition: librtcore.h:2627
enum rt_colormap_t::@10 method
rt_colormap_entry entry
Definition: librtcore.h:2642
uint16_t nentry
Definition: librtcore.h:2641
double *** values
Definition: librtcore.h:2605
uint32_t columns
Definition: librtcore.h:2601
uint16_t rasters
Definition: librtcore.h:2597
rt_raster raster
Definition: librtcore.h:2589
uint16_t nband
Definition: librtcore.h:2590
uint8_t nbnodata
Definition: librtcore.h:2591
double ** values
Definition: librtcore.h:2492
uint16_t dimy
Definition: librtcore.h:2491
uint16_t dimx
Definition: librtcore.h:2490
int ** nodata
Definition: librtcore.h:2493
int weighted
Definition: librtcore.h:2494
Struct definitions.
Definition: librtcore.h:2396
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
union rtpg_nmapalgebra_callback_arg::@21 ufc_info_data
FunctionCallInfoBaseData fcinfo
rtpg_nmapalgebraexpr_callback_arg callback
rtpg_nmapalgebra_arg bandarg
struct rtpg_nmapalgebraexpr_callback_arg::@24 kw
struct rtpg_nmapalgebraexpr_callback_arg::@23 nodatanodata
struct rtpg_nmapalgebraexpr_callback_arg::@22 expr[3]
rtpg_union_band_arg bandarg
rtpg_union_type uniontype