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