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