PostGIS  3.4.0dev-r@@SVN_REVISION@@
rtpg_pixel.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  *
14  * This program is free software; you can redistribute it and/or
15  * modify it under the terms of the GNU General Public License
16  * as published by the Free Software Foundation; either version 2
17  * of the License, or (at your option) any later version.
18  *
19  * This program is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with this program; if not, write to the Free Software Foundation,
26  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
27  *
28  */
29 
30 #include <postgres.h>
31 #include <fmgr.h>
32 #include "utils/lsyscache.h" /* for get_typlenbyvalalign */
33 #include <funcapi.h>
34 #include "utils/array.h" /* for ArrayType */
35 #include "utils/builtins.h" /* for text_to_cstring */
36 #include "utils/formatting.h" /* for asc_tolower */
37 #include "catalog/pg_type.h" /* for INT2OID, INT4OID, FLOAT4OID, FLOAT8OID and TEXTOID */
38 
39 #include "../../postgis_config.h"
40 #include "lwgeom_pg.h"
41 
42 
43 
44 #include "access/htup_details.h" /* for heap_form_tuple() */
45 
46 
47 #include "rtpostgis.h"
48 
49 /* Get pixel value */
50 Datum RASTER_getPixelValue(PG_FUNCTION_ARGS);
51 Datum RASTER_getPixelValueResample(PG_FUNCTION_ARGS);
52 Datum RASTER_dumpValues(PG_FUNCTION_ARGS);
53 
54 /* Set pixel value(s) */
55 Datum RASTER_setPixelValue(PG_FUNCTION_ARGS);
56 Datum RASTER_setPixelValuesArray(PG_FUNCTION_ARGS);
57 Datum RASTER_setPixelValuesGeomval(PG_FUNCTION_ARGS);
58 Datum RASTER_getGeometryValues(PG_FUNCTION_ARGS);
59 
60 /* Get pixels of value */
61 Datum RASTER_pixelOfValue(PG_FUNCTION_ARGS);
62 
63 /* Get nearest value to a point */
64 Datum RASTER_nearestValue(PG_FUNCTION_ARGS);
65 
66 /* Get the neighborhood around a pixel */
67 Datum RASTER_neighborhood(PG_FUNCTION_ARGS);
68 
76 PG_FUNCTION_INFO_V1(RASTER_getPixelValue);
77 Datum RASTER_getPixelValue(PG_FUNCTION_ARGS)
78 {
79  rt_pgraster *pgraster = NULL;
80  rt_raster raster = NULL;
81  rt_band band = NULL;
82  double pixvalue = 0;
83  int32_t bandindex = 0;
84  int32_t x = 0;
85  int32_t y = 0;
86  int result = 0;
87  bool exclude_nodata_value = TRUE;
88  int isnodata = 0;
89 
90  /* Index is 1-based */
91  bandindex = PG_GETARG_INT32(1);
92  if ( bandindex < 1 ) {
93  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
94  PG_RETURN_NULL();
95  }
96 
97  x = PG_GETARG_INT32(2);
98 
99  y = PG_GETARG_INT32(3);
100 
101  exclude_nodata_value = PG_GETARG_BOOL(4);
102 
103  POSTGIS_RT_DEBUGF(3, "Pixel coordinates (%d, %d)", x, y);
104 
105  /* Deserialize raster */
106  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
107  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
108 
109  raster = rt_raster_deserialize(pgraster, FALSE);
110  if (!raster) {
111  PG_FREE_IF_COPY(pgraster, 0);
112  elog(ERROR, "RASTER_getPixelValue: Could not deserialize raster");
113  PG_RETURN_NULL();
114  }
115 
116  /* Fetch Nth band using 0-based internal index */
117  band = rt_raster_get_band(raster, bandindex - 1);
118  if (! band) {
119  elog(NOTICE, "Could not find raster band of index %d when getting pixel "
120  "value. Returning NULL", bandindex);
121  rt_raster_destroy(raster);
122  PG_FREE_IF_COPY(pgraster, 0);
123  PG_RETURN_NULL();
124  }
125  /* Fetch pixel using 0-based coordinates */
126  result = rt_band_get_pixel(band, x - 1, y - 1, &pixvalue, &isnodata);
127 
128  /* If the result is -1 or the value is nodata and we take nodata into account
129  * then return nodata = NULL */
130  if (result != ES_NONE || (exclude_nodata_value && isnodata)) {
131  rt_raster_destroy(raster);
132  PG_FREE_IF_COPY(pgraster, 0);
133  PG_RETURN_NULL();
134  }
135 
136  rt_raster_destroy(raster);
137  PG_FREE_IF_COPY(pgraster, 0);
138 
139  PG_RETURN_FLOAT8(pixvalue);
140 }
141 
142 static rt_resample_type resample_text_to_type(text *txt)
143 {
144  char *resample = asc_tolower(VARDATA(txt), VARSIZE_ANY_EXHDR(txt));
145  if (strncmp(resample, "bilinear", 8) == 0)
146  return RT_BILINEAR;
147  else if (strncmp(resample, "nearest", 7) == 0)
148  return RT_NEAREST;
149  else {
150  elog(ERROR, "Unknown resample type '%s' requested", resample);
151  }
152  pfree(resample);
153  return RT_NEAREST;
154 }
155 
156 /*
157 * ST_Value(
158 * rast raster,
159 * band integer,
160 * pt geometry,
161 * exclude_nodata_value boolean DEFAULT TRUE,
162 * resample text DEFAULT 'nearest'
163 */
164 PG_FUNCTION_INFO_V1(RASTER_getPixelValueResample);
165 Datum RASTER_getPixelValueResample(PG_FUNCTION_ARGS)
166 {
167  rt_raster raster = NULL;
168  rt_band band = NULL;
169  rt_pgraster *pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
170  int32_t bandnum = PG_GETARG_INT32(1);
171  GSERIALIZED *gser;
172  LWPOINT *lwpoint;
173  LWGEOM *lwgeom;
174  bool exclude_nodata_value = PG_GETARG_BOOL(3);
175  rt_resample_type resample_type = RT_NEAREST;
176  double x, y, xr, yr;
177  double pixvalue = 0.0;
178  int isnodata = 0;
179  rt_errorstate err;
180 
181  /* Index is 1-based */
182  if (bandnum < 1) {
183  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
184  PG_RETURN_NULL();
185  }
186 
187  gser = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(2));
188  if (gserialized_get_type(gser) != POINTTYPE || gserialized_is_empty(gser)) {
189  elog(ERROR, "Attempting to get the value of a pixel with a non-point geometry");
190  PG_RETURN_NULL();
191  }
192 
193  raster = rt_raster_deserialize(pgraster, FALSE);
194  if (!raster) {
195  elog(ERROR, "RASTER_getPixelValue: Could not deserialize raster");
196  PG_RETURN_NULL();
197  }
198 
199  if (gserialized_get_srid(gser) != rt_raster_get_srid(raster)) {
200  elog(ERROR, "Raster and geometry do not have the same SRID");
201  PG_RETURN_NULL();
202  }
203 
204  if (PG_NARGS() > 4) {
205  text *resample_text = PG_GETARG_TEXT_P(4);
206  resample_type = resample_text_to_type(resample_text);
207  }
208 
209  /* Fetch Nth band using 0-based internal index */
210  band = rt_raster_get_band(raster, bandnum - 1);
211  if (!band) {
212  elog(ERROR, "Could not find raster band of index %d when getting pixel "
213  "value. Returning NULL", bandnum);
214  PG_RETURN_NULL();
215  }
216 
217  /* Get the X/Y coordinates */
218  lwgeom = lwgeom_from_gserialized(gser);
219  lwpoint = lwgeom_as_lwpoint(lwgeom);
220  x = lwpoint_get_x(lwpoint);
221  y = lwpoint_get_y(lwpoint);
222 
223  /* Convert X/Y world coordinates into raster coordinates */
224  err = rt_raster_geopoint_to_rasterpoint(raster, x, y, &xr, &yr, NULL);
225  if (err != ES_NONE) {
226  elog(ERROR, "Could not convert world coordinate to raster coordinate");
227  PG_RETURN_NULL();
228  }
229 
230  /* Use appropriate resample algorithm */
231  err = rt_band_get_pixel_resample(
232  band,
233  xr, yr,
234  resample_type,
235  &pixvalue, &isnodata
236  );
237 
238  /* If the result is -1 or the value is nodata and we take nodata into account
239  * then return nodata = NULL */
240  rt_raster_destroy(raster);
241  lwgeom_free(lwgeom);
242  if (err != ES_NONE || (exclude_nodata_value && isnodata)) {
243  PG_RETURN_NULL();
244  }
245  PG_RETURN_FLOAT8(pixvalue);
246 }
247 
248 
249 /*
250 * ST_SetZ(
251 * rast raster,
252 * pt geometry,
253 * resample text DEFAULT 'nearest' ),
254 * band integer default 1,
255 * dimension text default 'z'
256 */
257 PG_FUNCTION_INFO_V1(RASTER_getGeometryValues);
258 Datum RASTER_getGeometryValues(PG_FUNCTION_ARGS)
259 {
260  rt_raster raster = NULL;
261  rt_pgraster *pgraster = NULL;
262  GSERIALIZED *gser;
263  LWGEOM *lwgeom_in, *lwgeom_out;
264  rt_resample_type resample_type = RT_NEAREST;
265  rt_errorstate err;
266  char dimension;
267  const char *func_name;
268  uint16_t num_bands;
269  int32_t band;
270 
271  text *resample_text = PG_GETARG_TEXT_P(2);
272 
273  /* Dimension depends on the name of calling SQL function */
274  /* ST_SetZ()? or ST_SetM()? */
275  func_name = get_func_name(fcinfo->flinfo->fn_oid);
276  if (strcmp(func_name, "st_setz") == 0)
277  dimension = 'z';
278  else if (strcmp(func_name, "st_setm") == 0)
279  dimension = 'm';
280  else
281  elog(ERROR, "%s called from unexpected SQL signature", __func__);
282 
283  /* Geometry */
284  gser = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
285  if (gserialized_is_empty(gser)) {
286  elog(ERROR, "Cannot copy value into an empty geometry");
287  PG_RETURN_NULL();
288  }
289 
290  /* Raster */
291  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
292  raster = rt_raster_deserialize(pgraster, FALSE);
293  num_bands = rt_raster_get_num_bands(raster);
294  if (!raster) {
295  elog(ERROR, "Could not deserialize raster");
296  PG_RETURN_NULL();
297  }
298 
299  /* Bandnidex is 1-based */
300  band = PG_GETARG_INT32(3);
301  if (band < 1 || band > num_bands) {
302  elog(NOTICE, "Invalid band index %d. Must be between 1 and %u", band, num_bands);
303  PG_RETURN_NULL();
304  }
305 
306  /* SRID consistency */
307  if (gserialized_get_srid(gser) != rt_raster_get_srid(raster)) {
308  elog(ERROR, "Raster and geometry do not have the same SRID");
309  PG_RETURN_NULL();
310  }
311 
312  /* Process arguments */
313  resample_type = resample_text_to_type(resample_text);
314 
315  /* Get the geometry */
316  lwgeom_in = lwgeom_from_gserialized(gser);
317 
318  /* Run the sample */
319  err = rt_raster_copy_to_geometry(
320  raster,
321  band - 1, /* rtcore uses 0-based band number */
322  dimension, /* 'z' or 'm' */
323  resample_type, /* bilinear or nearest */
324  lwgeom_in,
325  &lwgeom_out
326  );
327 
328  rt_raster_destroy(raster);
329  lwgeom_free(lwgeom_in);
330  if (err != ES_NONE || !lwgeom_out) {
331  PG_RETURN_NULL();
332  }
333  PG_RETURN_POINTER(gserialized_from_lwgeom(lwgeom_out, NULL));
334 }
335 
336 
337 /* ---------------------------------------------------------------- */
338 /* ST_DumpValues function */
339 /* ---------------------------------------------------------------- */
340 
341 typedef struct rtpg_dumpvalues_arg_t *rtpg_dumpvalues_arg;
342 struct rtpg_dumpvalues_arg_t {
343  int numbands;
344  int rows;
345  int columns;
346 
347  int *nbands; /* 0-based */
348  Datum **values;
349  bool **nodata;
350 };
351 
352 static rtpg_dumpvalues_arg rtpg_dumpvalues_arg_init() {
353  rtpg_dumpvalues_arg arg = NULL;
354 
355  arg = palloc(sizeof(struct rtpg_dumpvalues_arg_t));
356  if (arg == NULL) {
357  elog(ERROR, "rtpg_dumpvalues_arg_init: Could not allocate memory for arguments");
358  return NULL;
359  }
360 
361  arg->numbands = 0;
362  arg->rows = 0;
363  arg->columns = 0;
364 
365  arg->nbands = NULL;
366  arg->values = NULL;
367  arg->nodata = NULL;
368 
369  return arg;
370 }
371 
372 static void rtpg_dumpvalues_arg_destroy(rtpg_dumpvalues_arg arg) {
373  int i = 0;
374 
375  if (arg->numbands > 0) {
376  if (arg->nbands != NULL)
377  pfree(arg->nbands);
378 
379  if (arg->values != NULL) {
380  for (i = 0; i < arg->numbands; i++) {
381 
382  if (arg->values[i] != NULL)
383  pfree(arg->values[i]);
384 
385  if (arg->nodata[i] != NULL)
386  pfree(arg->nodata[i]);
387  }
388 
389  pfree(arg->values);
390  }
391 
392  if (arg->nodata != NULL)
393  pfree(arg->nodata);
394  }
395 
396  pfree(arg);
397 }
398 
399 #define VALUES_LENGTH 2
400 
401 PG_FUNCTION_INFO_V1(RASTER_dumpValues);
402 Datum RASTER_dumpValues(PG_FUNCTION_ARGS)
403 {
404  FuncCallContext *funcctx;
405  TupleDesc tupdesc;
406  int call_cntr;
407  int max_calls;
408  int i = 0;
409  int x = 0;
410  int y = 0;
411  int z = 0;
412 
413  int16 typlen;
414  bool typbyval;
415  char typalign;
416 
417  rtpg_dumpvalues_arg arg1 = NULL;
418  rtpg_dumpvalues_arg arg2 = NULL;
419 
420  /* stuff done only on the first call of the function */
421  if (SRF_IS_FIRSTCALL()) {
422  MemoryContext oldcontext;
423  rt_pgraster *pgraster = NULL;
424  rt_raster raster = NULL;
425  rt_band band = NULL;
426  int numbands = 0;
427  int j = 0;
428  bool exclude_nodata_value = TRUE;
429 
430  ArrayType *array;
431  Oid etype;
432  Datum *e;
433  bool *nulls;
434 
435  double val = 0;
436  int isnodata = 0;
437 
438  POSTGIS_RT_DEBUG(2, "RASTER_dumpValues first call");
439 
440  /* create a function context for cross-call persistence */
441  funcctx = SRF_FIRSTCALL_INIT();
442 
443  /* switch to memory context appropriate for multiple function calls */
444  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
445 
446  /* Get input arguments */
447  if (PG_ARGISNULL(0)) {
448  MemoryContextSwitchTo(oldcontext);
449  SRF_RETURN_DONE(funcctx);
450  }
451  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
452 
453  raster = rt_raster_deserialize(pgraster, FALSE);
454  if (!raster) {
455  PG_FREE_IF_COPY(pgraster, 0);
456  ereport(ERROR, (
457  errcode(ERRCODE_OUT_OF_MEMORY),
458  errmsg("Could not deserialize raster")
459  ));
460  MemoryContextSwitchTo(oldcontext);
461  SRF_RETURN_DONE(funcctx);
462  }
463 
464  /* check that raster is not empty */
465  /*
466  if (rt_raster_is_empty(raster)) {
467  elog(NOTICE, "Raster provided is empty");
468  rt_raster_destroy(raster);
469  PG_FREE_IF_COPY(pgraster, 0);
470  MemoryContextSwitchTo(oldcontext);
471  SRF_RETURN_DONE(funcctx);
472  }
473  */
474 
475  /* raster has bands */
476  numbands = rt_raster_get_num_bands(raster);
477  if (!numbands) {
478  elog(NOTICE, "Raster provided has no bands");
479  rt_raster_destroy(raster);
480  PG_FREE_IF_COPY(pgraster, 0);
481  MemoryContextSwitchTo(oldcontext);
482  SRF_RETURN_DONE(funcctx);
483  }
484 
485  /* initialize arg1 */
486  arg1 = rtpg_dumpvalues_arg_init();
487  if (arg1 == NULL) {
488  rt_raster_destroy(raster);
489  PG_FREE_IF_COPY(pgraster, 0);
490  MemoryContextSwitchTo(oldcontext);
491  elog(ERROR, "RASTER_dumpValues: Could not initialize argument structure");
492  SRF_RETURN_DONE(funcctx);
493  }
494 
495  /* nband, array */
496  if (!PG_ARGISNULL(1)) {
497  array = PG_GETARG_ARRAYTYPE_P(1);
498  etype = ARR_ELEMTYPE(array);
499  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
500 
501  switch (etype) {
502  case INT2OID:
503  case INT4OID:
504  break;
505  default:
506  rtpg_dumpvalues_arg_destroy(arg1);
507  rt_raster_destroy(raster);
508  PG_FREE_IF_COPY(pgraster, 0);
509  MemoryContextSwitchTo(oldcontext);
510  elog(ERROR, "RASTER_dumpValues: Invalid data type for band indexes");
511  SRF_RETURN_DONE(funcctx);
512  break;
513  }
514 
515  deconstruct_array(array, etype, typlen, typbyval, typalign, &e, &nulls, &(arg1->numbands));
516 
517  arg1->nbands = palloc(sizeof(int) * arg1->numbands);
518  if (arg1->nbands == NULL) {
519  rtpg_dumpvalues_arg_destroy(arg1);
520  rt_raster_destroy(raster);
521  PG_FREE_IF_COPY(pgraster, 0);
522  MemoryContextSwitchTo(oldcontext);
523  elog(ERROR, "RASTER_dumpValues: Could not allocate memory for band indexes");
524  SRF_RETURN_DONE(funcctx);
525  }
526 
527  for (i = 0, j = 0; i < arg1->numbands; i++) {
528  if (nulls[i]) continue;
529 
530  switch (etype) {
531  case INT2OID:
532  arg1->nbands[j] = DatumGetInt16(e[i]) - 1;
533  break;
534  case INT4OID:
535  arg1->nbands[j] = DatumGetInt32(e[i]) - 1;
536  break;
537  }
538 
539  j++;
540  }
541 
542  if (j < arg1->numbands) {
543  arg1->nbands = repalloc(arg1->nbands, sizeof(int) * j);
544  if (arg1->nbands == NULL) {
545  rtpg_dumpvalues_arg_destroy(arg1);
546  rt_raster_destroy(raster);
547  PG_FREE_IF_COPY(pgraster, 0);
548  MemoryContextSwitchTo(oldcontext);
549  elog(ERROR, "RASTER_dumpValues: Could not reallocate memory for band indexes");
550  SRF_RETURN_DONE(funcctx);
551  }
552 
553  arg1->numbands = j;
554  }
555 
556  /* validate nbands */
557  for (i = 0; i < arg1->numbands; i++) {
558  if (!rt_raster_has_band(raster, arg1->nbands[i])) {
559  elog(NOTICE, "Band at index %d not found in raster", arg1->nbands[i] + 1);
560  rtpg_dumpvalues_arg_destroy(arg1);
561  rt_raster_destroy(raster);
562  PG_FREE_IF_COPY(pgraster, 0);
563  MemoryContextSwitchTo(oldcontext);
564  SRF_RETURN_DONE(funcctx);
565  }
566  }
567 
568  }
569  /* no bands specified, return all bands */
570  else {
571  arg1->numbands = numbands;
572  arg1->nbands = palloc(sizeof(int) * arg1->numbands);
573 
574  if (arg1->nbands == NULL) {
575  rtpg_dumpvalues_arg_destroy(arg1);
576  rt_raster_destroy(raster);
577  PG_FREE_IF_COPY(pgraster, 0);
578  MemoryContextSwitchTo(oldcontext);
579  elog(ERROR, "RASTER_dumpValues: Could not allocate memory for band indexes");
580  SRF_RETURN_DONE(funcctx);
581  }
582 
583  for (i = 0; i < arg1->numbands; i++) {
584  arg1->nbands[i] = i;
585  POSTGIS_RT_DEBUGF(4, "arg1->nbands[%d] = %d", arg1->nbands[i], i);
586  }
587  }
588 
589  arg1->rows = rt_raster_get_height(raster);
590  arg1->columns = rt_raster_get_width(raster);
591 
592  /* exclude_nodata_value */
593  if (!PG_ARGISNULL(2))
594  exclude_nodata_value = PG_GETARG_BOOL(2);
595  POSTGIS_RT_DEBUGF(4, "exclude_nodata_value = %d", exclude_nodata_value);
596 
597  /* allocate memory for each band's values and nodata flags */
598  arg1->values = palloc(sizeof(Datum *) * arg1->numbands);
599  arg1->nodata = palloc(sizeof(bool *) * arg1->numbands);
600  if (arg1->values == NULL || arg1->nodata == NULL) {
601  rtpg_dumpvalues_arg_destroy(arg1);
602  rt_raster_destroy(raster);
603  PG_FREE_IF_COPY(pgraster, 0);
604  MemoryContextSwitchTo(oldcontext);
605  elog(ERROR, "RASTER_dumpValues: Could not allocate memory for pixel values");
606  SRF_RETURN_DONE(funcctx);
607  }
608  memset(arg1->values, 0, sizeof(Datum *) * arg1->numbands);
609  memset(arg1->nodata, 0, sizeof(bool *) * arg1->numbands);
610 
611  /* get each band and dump data */
612  for (z = 0; z < arg1->numbands; z++) {
613  /* shortcut if raster is empty */
614  if (rt_raster_is_empty(raster))
615  break;
616 
617  band = rt_raster_get_band(raster, arg1->nbands[z]);
618  if (!band) {
619  int nband = arg1->nbands[z] + 1;
620  rtpg_dumpvalues_arg_destroy(arg1);
621  rt_raster_destroy(raster);
622  PG_FREE_IF_COPY(pgraster, 0);
623  MemoryContextSwitchTo(oldcontext);
624  elog(ERROR, "RASTER_dumpValues: Could not get band at index %d", nband);
625  SRF_RETURN_DONE(funcctx);
626  }
627 
628  /* allocate memory for values and nodata flags */
629  arg1->values[z] = palloc(sizeof(Datum) * arg1->rows * arg1->columns);
630  arg1->nodata[z] = palloc(sizeof(bool) * arg1->rows * arg1->columns);
631  if (arg1->values[z] == NULL || arg1->nodata[z] == NULL) {
632  rtpg_dumpvalues_arg_destroy(arg1);
633  rt_raster_destroy(raster);
634  PG_FREE_IF_COPY(pgraster, 0);
635  MemoryContextSwitchTo(oldcontext);
636  elog(ERROR, "RASTER_dumpValues: Could not allocate memory for pixel values");
637  SRF_RETURN_DONE(funcctx);
638  }
639  memset(arg1->values[z], 0, sizeof(Datum) * arg1->rows * arg1->columns);
640  memset(arg1->nodata[z], 0, sizeof(bool) * arg1->rows * arg1->columns);
641 
642  i = 0;
643 
644  /* shortcut if band is NODATA */
645  if (rt_band_get_isnodata_flag(band)) {
646  for (i = (arg1->rows * arg1->columns) - 1; i >= 0; i--)
647  arg1->nodata[z][i] = TRUE;
648  continue;
649  }
650 
651  for (y = 0; y < arg1->rows; y++) {
652  for (x = 0; x < arg1->columns; x++) {
653  /* get pixel */
654  if (rt_band_get_pixel(band, x, y, &val, &isnodata) != ES_NONE) {
655  int nband = arg1->nbands[z] + 1;
656  rtpg_dumpvalues_arg_destroy(arg1);
657  rt_raster_destroy(raster);
658  PG_FREE_IF_COPY(pgraster, 0);
659  MemoryContextSwitchTo(oldcontext);
660  elog(ERROR, "RASTER_dumpValues: Could not pixel (%d, %d) of band %d", x, y, nband);
661  SRF_RETURN_DONE(funcctx);
662  }
663 
664  arg1->values[z][i] = Float8GetDatum(val);
665  POSTGIS_RT_DEBUGF(5, "arg1->values[z][i] = %f", DatumGetFloat8(arg1->values[z][i]));
666  POSTGIS_RT_DEBUGF(5, "clamped is?: %d", rt_band_clamped_value_is_nodata(band, val));
667 
668  if (exclude_nodata_value && isnodata) {
669  arg1->nodata[z][i] = TRUE;
670  POSTGIS_RT_DEBUG(5, "nodata = 1");
671  }
672  else
673  POSTGIS_RT_DEBUG(5, "nodata = 0");
674 
675  i++;
676  }
677  }
678  }
679 
680  /* cleanup */
681  rt_raster_destroy(raster);
682  PG_FREE_IF_COPY(pgraster, 0);
683 
684  /* Store needed information */
685  funcctx->user_fctx = arg1;
686 
687  /* total number of tuples to be returned */
688  funcctx->max_calls = arg1->numbands;
689 
690  /* Build a tuple descriptor for our result type */
691  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
692  MemoryContextSwitchTo(oldcontext);
693  ereport(ERROR, (
694  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
695  errmsg(
696  "function returning record called in context "
697  "that cannot accept type record"
698  )
699  ));
700  }
701 
702  BlessTupleDesc(tupdesc);
703  funcctx->tuple_desc = tupdesc;
704 
705  MemoryContextSwitchTo(oldcontext);
706  }
707 
708  /* stuff done on every call of the function */
709  funcctx = SRF_PERCALL_SETUP();
710 
711  call_cntr = funcctx->call_cntr;
712  max_calls = funcctx->max_calls;
713  tupdesc = funcctx->tuple_desc;
714  arg2 = funcctx->user_fctx;
715 
716  /* do when there is more left to send */
717  if (call_cntr < max_calls) {
718  Datum values[VALUES_LENGTH];
719  bool nulls[VALUES_LENGTH];
720  HeapTuple tuple;
721  Datum result;
722  ArrayType *mdValues = NULL;
723  int ndim = 2;
724  int dim[2] = {arg2->rows, arg2->columns};
725  int lbound[2] = {1, 1};
726 
727  POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr);
728  POSTGIS_RT_DEBUGF(4, "dim = %d, %d", dim[0], dim[1]);
729 
730  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
731 
732  values[0] = Int32GetDatum(arg2->nbands[call_cntr] + 1);
733 
734  /* info about the type of item in the multi-dimensional array (float8). */
735  get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
736 
737  /* if values is NULL, return empty array */
738  if (arg2->values[call_cntr] == NULL)
739  ndim = 0;
740 
741  /* assemble 3-dimension array of values */
742  mdValues = construct_md_array(
743  arg2->values[call_cntr], arg2->nodata[call_cntr],
744  ndim, dim, lbound,
745  FLOAT8OID,
746  typlen, typbyval, typalign
747  );
748  values[1] = PointerGetDatum(mdValues);
749 
750  /* build a tuple and datum */
751  tuple = heap_form_tuple(tupdesc, values, nulls);
752  result = HeapTupleGetDatum(tuple);
753 
754  SRF_RETURN_NEXT(funcctx, result);
755  }
756  /* do when there is no more left */
757  else {
758  rtpg_dumpvalues_arg_destroy(arg2);
759  SRF_RETURN_DONE(funcctx);
760  }
761 }
762 
766 PG_FUNCTION_INFO_V1(RASTER_setPixelValue);
767 Datum RASTER_setPixelValue(PG_FUNCTION_ARGS)
768 {
769  rt_pgraster *pgraster = NULL;
770  rt_pgraster *pgrtn = NULL;
771  rt_raster raster = NULL;
772  rt_band band = NULL;
773  double pixvalue = 0;
774  int32_t bandindex = 0;
775  int32_t x = 0;
776  int32_t y = 0;
777  bool skipset = FALSE;
778 
779  if (PG_ARGISNULL(0))
780  PG_RETURN_NULL();
781 
782  /* Check index is not NULL or < 1 */
783  if (PG_ARGISNULL(1))
784  bandindex = -1;
785  else
786  bandindex = PG_GETARG_INT32(1);
787 
788  if (bandindex < 1) {
789  elog(NOTICE, "Invalid band index (must use 1-based). Value not set. Returning original raster");
790  skipset = TRUE;
791  }
792 
793  /* Validate pixel coordinates are not null */
794  if (PG_ARGISNULL(2)) {
795  elog(NOTICE, "X coordinate can not be NULL when setting pixel value. Value not set. Returning original raster");
796  skipset = TRUE;
797  }
798  else
799  x = PG_GETARG_INT32(2);
800 
801  if (PG_ARGISNULL(3)) {
802  elog(NOTICE, "Y coordinate can not be NULL when setting pixel value. Value not set. Returning original raster");
803  skipset = TRUE;
804  }
805  else
806  y = PG_GETARG_INT32(3);
807 
808  POSTGIS_RT_DEBUGF(3, "Pixel coordinates (%d, %d)", x, y);
809 
810  /* Deserialize raster */
811  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
812 
813  raster = rt_raster_deserialize(pgraster, FALSE);
814  if (!raster) {
815  PG_FREE_IF_COPY(pgraster, 0);
816  elog(ERROR, "RASTER_setPixelValue: Could not deserialize raster");
817  PG_RETURN_NULL();
818  }
819 
820  if (!skipset) {
821  /* Fetch requested band */
822  band = rt_raster_get_band(raster, bandindex - 1);
823  if (!band) {
824  elog(NOTICE, "Could not find raster band of index %d when setting "
825  "pixel value. Value not set. Returning original raster",
826  bandindex);
827  PG_RETURN_POINTER(pgraster);
828  }
829  else {
830  /* Set the pixel value */
831  if (PG_ARGISNULL(4)) {
832  if (!rt_band_get_hasnodata_flag(band)) {
833  elog(NOTICE, "Raster do not have a nodata value defined. "
834  "Set band nodata value first. Nodata value not set. "
835  "Returning original raster");
836  PG_RETURN_POINTER(pgraster);
837  }
838  else {
839  rt_band_get_nodata(band, &pixvalue);
840  rt_band_set_pixel(band, x - 1, y - 1, pixvalue, NULL);
841  }
842  }
843  else {
844  pixvalue = PG_GETARG_FLOAT8(4);
845  rt_band_set_pixel(band, x - 1, y - 1, pixvalue, NULL);
846  }
847  }
848  }
849 
850  pgrtn = rt_raster_serialize(raster);
851  rt_raster_destroy(raster);
852  PG_FREE_IF_COPY(pgraster, 0);
853  if (!pgrtn)
854  PG_RETURN_NULL();
855 
856  SET_VARSIZE(pgrtn, pgrtn->size);
857  PG_RETURN_POINTER(pgrtn);
858 }
859 
863 PG_FUNCTION_INFO_V1(RASTER_setPixelValuesArray);
864 Datum RASTER_setPixelValuesArray(PG_FUNCTION_ARGS)
865 {
866  rt_pgraster *pgraster = NULL;
867  rt_pgraster *pgrtn = NULL;
868  rt_raster raster = NULL;
869  rt_band band = NULL;
870  int numbands = 0;
871 
872  int nband = 0;
873  int width = 0;
874  int height = 0;
875 
876  ArrayType *array;
877  Oid etype;
878  Datum *elements;
879  bool *nulls;
880  int16 typlen;
881  bool typbyval;
882  char typalign;
883  int ndims = 1;
884  int *dims;
885  int num = 0;
886 
887  int ul[2] = {0};
888  struct pixelvalue {
889  int x;
890  int y;
891 
892  bool noset;
893  bool nodata;
894  double value;
895  };
896  struct pixelvalue *pixval = NULL;
897  int numpixval = 0;
898  int dimpixval[2] = {1, 1};
899  int dimnoset[2] = {1, 1};
900  int hasnodata = FALSE;
901  double nodataval = 0;
902  bool keepnodata = FALSE;
903  bool hasnosetval = FALSE;
904  bool nosetvalisnull = FALSE;
905  double nosetval = 0;
906 
907  int rtn = 0;
908  double val = 0;
909  int isnodata = 0;
910 
911  int i = 0;
912  int j = 0;
913  int x = 0;
914  int y = 0;
915 
916  /* pgraster is null, return null */
917  if (PG_ARGISNULL(0))
918  PG_RETURN_NULL();
919  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
920 
921  /* raster */
922  raster = rt_raster_deserialize(pgraster, FALSE);
923  if (!raster) {
924  PG_FREE_IF_COPY(pgraster, 0);
925  elog(ERROR, "RASTER_setPixelValuesArray: Could not deserialize raster");
926  PG_RETURN_NULL();
927  }
928 
929  /* raster attributes */
930  numbands = rt_raster_get_num_bands(raster);
931  width = rt_raster_get_width(raster);
932  height = rt_raster_get_height(raster);
933 
934  /* nband */
935  if (PG_ARGISNULL(1)) {
936  elog(NOTICE, "Band index cannot be NULL. Value must be 1-based. Returning original raster");
937  rt_raster_destroy(raster);
938  PG_RETURN_POINTER(pgraster);
939  }
940 
941  nband = PG_GETARG_INT32(1);
942  if (nband < 1 || nband > numbands) {
943  elog(NOTICE, "Band index is invalid. Value must be 1-based. Returning original raster");
944  rt_raster_destroy(raster);
945  PG_RETURN_POINTER(pgraster);
946  }
947 
948  /* x, y */
949  for (i = 2, j = 0; i < 4; i++, j++) {
950  if (PG_ARGISNULL(i)) {
951  elog(NOTICE, "%s cannot be NULL. Value must be 1-based. Returning original raster", j < 1 ? "X" : "Y");
952  rt_raster_destroy(raster);
953  PG_RETURN_POINTER(pgraster);
954  }
955 
956  ul[j] = PG_GETARG_INT32(i);
957  if (
958  (ul[j] < 1) || (
959  (j < 1 && ul[j] > width) ||
960  (j > 0 && ul[j] > height)
961  )
962  ) {
963  elog(NOTICE, "%s is invalid. Value must be 1-based. Returning original raster", j < 1 ? "X" : "Y");
964  rt_raster_destroy(raster);
965  PG_RETURN_POINTER(pgraster);
966  }
967 
968  /* force 0-based from 1-based */
969  ul[j] -= 1;
970  }
971 
972  /* new value set */
973  if (PG_ARGISNULL(4)) {
974  elog(NOTICE, "No values to set. Returning original raster");
975  rt_raster_destroy(raster);
976  PG_RETURN_POINTER(pgraster);
977  }
978 
979  array = PG_GETARG_ARRAYTYPE_P(4);
980  etype = ARR_ELEMTYPE(array);
981  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
982 
983  switch (etype) {
984  case FLOAT4OID:
985  case FLOAT8OID:
986  break;
987  default:
988  rt_raster_destroy(raster);
989  PG_FREE_IF_COPY(pgraster, 0);
990  elog(ERROR, "RASTER_setPixelValuesArray: Invalid data type for new values");
991  PG_RETURN_NULL();
992  break;
993  }
994 
995  ndims = ARR_NDIM(array);
996  dims = ARR_DIMS(array);
997  POSTGIS_RT_DEBUGF(4, "ndims = %d", ndims);
998 
999  if (ndims < 1 || ndims > 2) {
1000  elog(NOTICE, "New values array must be of 1 or 2 dimensions. Returning original raster");
1001  rt_raster_destroy(raster);
1002  PG_RETURN_POINTER(pgraster);
1003  }
1004  /* outer element, then inner element */
1005  /* i = 0, y */
1006  /* i = 1, x */
1007  if (ndims != 2)
1008  dimpixval[1] = dims[0];
1009  else {
1010  dimpixval[0] = dims[0];
1011  dimpixval[1] = dims[1];
1012  }
1013  POSTGIS_RT_DEBUGF(4, "dimpixval = (%d, %d)", dimpixval[0], dimpixval[1]);
1014 
1015  deconstruct_array(
1016  array,
1017  etype,
1018  typlen, typbyval, typalign,
1019  &elements, &nulls, &num
1020  );
1021 
1022  /* # of elements doesn't match dims */
1023  if (num < 1 || num != (dimpixval[0] * dimpixval[1])) {
1024  if (num) {
1025  pfree(elements);
1026  pfree(nulls);
1027  }
1028  rt_raster_destroy(raster);
1029  PG_FREE_IF_COPY(pgraster, 0);
1030  elog(ERROR, "RASTER_setPixelValuesArray: Could not deconstruct new values array");
1031  PG_RETURN_NULL();
1032  }
1033 
1034  /* allocate memory for pixval */
1035  numpixval = num;
1036  pixval = palloc(sizeof(struct pixelvalue) * numpixval);
1037  if (pixval == NULL) {
1038  pfree(elements);
1039  pfree(nulls);
1040  rt_raster_destroy(raster);
1041  PG_FREE_IF_COPY(pgraster, 0);
1042  elog(ERROR, "RASTER_setPixelValuesArray: Could not allocate memory for new pixel values");
1043  PG_RETURN_NULL();
1044  }
1045 
1046  /* load new values into pixval */
1047  i = 0;
1048  for (y = 0; y < dimpixval[0]; y++) {
1049  for (x = 0; x < dimpixval[1]; x++) {
1050  /* 0-based */
1051  pixval[i].x = ul[0] + x;
1052  pixval[i].y = ul[1] + y;
1053 
1054  pixval[i].noset = FALSE;
1055  pixval[i].nodata = FALSE;
1056  pixval[i].value = 0;
1057 
1058  if (nulls[i])
1059  pixval[i].nodata = TRUE;
1060  else {
1061  switch (etype) {
1062  case FLOAT4OID:
1063  pixval[i].value = DatumGetFloat4(elements[i]);
1064  break;
1065  case FLOAT8OID:
1066  pixval[i].value = DatumGetFloat8(elements[i]);
1067  break;
1068  }
1069  }
1070 
1071  i++;
1072  }
1073  }
1074 
1075  pfree(elements);
1076  pfree(nulls);
1077 
1078  /* now load noset flags */
1079  if (!PG_ARGISNULL(5)) {
1080  array = PG_GETARG_ARRAYTYPE_P(5);
1081  etype = ARR_ELEMTYPE(array);
1082  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1083 
1084  switch (etype) {
1085  case BOOLOID:
1086  break;
1087  default:
1088  pfree(pixval);
1089  rt_raster_destroy(raster);
1090  PG_FREE_IF_COPY(pgraster, 0);
1091  elog(ERROR, "RASTER_setPixelValuesArray: Invalid data type for noset flags");
1092  PG_RETURN_NULL();
1093  break;
1094  }
1095 
1096  ndims = ARR_NDIM(array);
1097  dims = ARR_DIMS(array);
1098  POSTGIS_RT_DEBUGF(4, "ndims = %d", ndims);
1099 
1100  if (ndims < 1 || ndims > 2) {
1101  elog(NOTICE, "Noset flags array must be of 1 or 2 dimensions. Returning original raster");
1102  pfree(pixval);
1103  rt_raster_destroy(raster);
1104  PG_RETURN_POINTER(pgraster);
1105  }
1106  /* outer element, then inner element */
1107  /* i = 0, y */
1108  /* i = 1, x */
1109  if (ndims != 2)
1110  dimnoset[1] = dims[0];
1111  else {
1112  dimnoset[0] = dims[0];
1113  dimnoset[1] = dims[1];
1114  }
1115  POSTGIS_RT_DEBUGF(4, "dimnoset = (%d, %d)", dimnoset[0], dimnoset[1]);
1116 
1117  deconstruct_array(
1118  array,
1119  etype,
1120  typlen, typbyval, typalign,
1121  &elements, &nulls, &num
1122  );
1123 
1124  /* # of elements doesn't match dims */
1125  if (num < 1 || num != (dimnoset[0] * dimnoset[1])) {
1126  pfree(pixval);
1127  if (num) {
1128  pfree(elements);
1129  pfree(nulls);
1130  }
1131  rt_raster_destroy(raster);
1132  PG_FREE_IF_COPY(pgraster, 0);
1133  elog(ERROR, "RASTER_setPixelValuesArray: Could not deconstruct noset flags array");
1134  PG_RETURN_NULL();
1135  }
1136 
1137  i = 0;
1138  j = 0;
1139  for (y = 0; y < dimnoset[0]; y++) {
1140  if (y >= dimpixval[0]) break;
1141 
1142  for (x = 0; x < dimnoset[1]; x++) {
1143  /* fast forward noset elements */
1144  if (x >= dimpixval[1]) {
1145  i += (dimnoset[1] - dimpixval[1]);
1146  break;
1147  }
1148 
1149  if (!nulls[i] && DatumGetBool(elements[i]))
1150  pixval[j].noset = TRUE;
1151 
1152  i++;
1153  j++;
1154  }
1155 
1156  /* fast forward pixval */
1157  if (x < dimpixval[1])
1158  j += (dimpixval[1] - dimnoset[1]);
1159  }
1160 
1161  pfree(elements);
1162  pfree(nulls);
1163  }
1164  /* hasnosetvalue and nosetvalue */
1165  else if (!PG_ARGISNULL(6) && PG_GETARG_BOOL(6)) {
1166  hasnosetval = TRUE;
1167  if (PG_ARGISNULL(7))
1168  nosetvalisnull = TRUE;
1169  else
1170  nosetval = PG_GETARG_FLOAT8(7);
1171  }
1172 
1173 #if POSTGIS_DEBUG_LEVEL > 0
1174  for (i = 0; i < numpixval; i++) {
1175  POSTGIS_RT_DEBUGF(4, "pixval[%d](x, y, noset, nodata, value) = (%d, %d, %d, %d, %f)",
1176  i,
1177  pixval[i].x,
1178  pixval[i].y,
1179  pixval[i].noset,
1180  pixval[i].nodata,
1181  pixval[i].value
1182  );
1183  }
1184 #endif
1185 
1186  /* keepnodata flag */
1187  if (!PG_ARGISNULL(8))
1188  keepnodata = PG_GETARG_BOOL(8);
1189 
1190  /* get band */
1191  band = rt_raster_get_band(raster, nband - 1);
1192  if (!band) {
1193  elog(NOTICE, "Could not find band at index %d. Returning original raster", nband);
1194  pfree(pixval);
1195  rt_raster_destroy(raster);
1196  PG_RETURN_POINTER(pgraster);
1197  }
1198 
1199  /* get band nodata info */
1200  /* has NODATA, use NODATA */
1201  hasnodata = rt_band_get_hasnodata_flag(band);
1202  if (hasnodata)
1203  rt_band_get_nodata(band, &nodataval);
1204  /* no NODATA, use min possible value */
1205  else
1206  nodataval = rt_band_get_min_value(band);
1207 
1208  /* set pixels */
1209  for (i = 0; i < numpixval; i++) {
1210  /* noset = true, skip */
1211  if (pixval[i].noset)
1212  continue;
1213  /* check against nosetval */
1214  else if (hasnosetval) {
1215  /* pixel = NULL AND nosetval = NULL */
1216  if (pixval[i].nodata && nosetvalisnull)
1217  continue;
1218  /* pixel value = nosetval */
1219  else if (!pixval[i].nodata && !nosetvalisnull && FLT_EQ(pixval[i].value, nosetval))
1220  continue;
1221  }
1222 
1223  /* if pixel is outside bounds, skip */
1224  if (
1225  (pixval[i].x < 0 || pixval[i].x >= width) ||
1226  (pixval[i].y < 0 || pixval[i].y >= height)
1227  ) {
1228  elog(NOTICE, "Cannot set value for pixel (%d, %d) outside raster bounds: %d x %d",
1229  pixval[i].x + 1, pixval[i].y + 1,
1230  width, height
1231  );
1232  continue;
1233  }
1234 
1235  /* if hasnodata = TRUE and keepnodata = TRUE, inspect pixel value */
1236  if (hasnodata && keepnodata) {
1237  rtn = rt_band_get_pixel(band, pixval[i].x, pixval[i].y, &val, &isnodata);
1238  if (rtn != ES_NONE) {
1239  pfree(pixval);
1240  rt_raster_destroy(raster);
1241  PG_FREE_IF_COPY(pgraster, 0);
1242  elog(ERROR, "Cannot get value of pixel");
1243  PG_RETURN_NULL();
1244  }
1245 
1246  /* pixel value = NODATA, skip */
1247  if (isnodata) {
1248  continue;
1249  }
1250  }
1251 
1252  if (pixval[i].nodata)
1253  rt_band_set_pixel(band, pixval[i].x, pixval[i].y, nodataval, NULL);
1254  else
1255  rt_band_set_pixel(band, pixval[i].x, pixval[i].y, pixval[i].value, NULL);
1256  }
1257 
1258  pfree(pixval);
1259 
1260  /* serialize new raster */
1261  pgrtn = rt_raster_serialize(raster);
1262  rt_raster_destroy(raster);
1263  PG_FREE_IF_COPY(pgraster, 0);
1264  if (!pgrtn)
1265  PG_RETURN_NULL();
1266 
1267  SET_VARSIZE(pgrtn, pgrtn->size);
1268  PG_RETURN_POINTER(pgrtn);
1269 }
1270 
1271 /* ---------------------------------------------------------------- */
1272 /* ST_SetValues using geomval array */
1273 /* ---------------------------------------------------------------- */
1274 
1275 typedef struct rtpg_setvaluesgv_arg_t *rtpg_setvaluesgv_arg;
1276 typedef struct rtpg_setvaluesgv_geomval_t *rtpg_setvaluesgv_geomval;
1277 
1278 struct rtpg_setvaluesgv_arg_t {
1279  int ngv;
1280  rtpg_setvaluesgv_geomval gv;
1281 
1282  bool keepnodata;
1283 };
1284 
1285 struct rtpg_setvaluesgv_geomval_t {
1286  struct {
1287  int nodata;
1288  double value;
1289  } pixval;
1290 
1291  LWGEOM *geom;
1292  rt_raster mask;
1293 };
1294 
1295 static rtpg_setvaluesgv_arg rtpg_setvaluesgv_arg_init() {
1296  rtpg_setvaluesgv_arg arg = palloc(sizeof(struct rtpg_setvaluesgv_arg_t));
1297  if (arg == NULL) {
1298  elog(ERROR, "rtpg_setvaluesgv_arg_init: Could not allocate memory for function arguments");
1299  return NULL;
1300  }
1301 
1302  arg->ngv = 0;
1303  arg->gv = NULL;
1304  arg->keepnodata = 0;
1305 
1306  return arg;
1307 }
1308 
1309 static void rtpg_setvaluesgv_arg_destroy(rtpg_setvaluesgv_arg arg) {
1310  int i = 0;
1311 
1312  if (arg->gv != NULL) {
1313  for (i = 0; i < arg->ngv; i++) {
1314  if (arg->gv[i].geom != NULL)
1315  lwgeom_free(arg->gv[i].geom);
1316  if (arg->gv[i].mask != NULL)
1317  rt_raster_destroy(arg->gv[i].mask);
1318  }
1319 
1320  pfree(arg->gv);
1321  }
1322 
1323  pfree(arg);
1324 }
1325 
1326 static int rtpg_setvalues_geomval_callback(
1327  rt_iterator_arg arg, void *userarg,
1328  double *value, int *nodata
1329 ) {
1330  rtpg_setvaluesgv_arg funcarg = (rtpg_setvaluesgv_arg) userarg;
1331  int i = 0;
1332  int j = 0;
1333 
1334  *value = 0;
1335  *nodata = 0;
1336 
1337  POSTGIS_RT_DEBUGF(4, "keepnodata = %d", funcarg->keepnodata);
1338 
1339  /* keepnodata = TRUE */
1340  if (funcarg->keepnodata && arg->nodata[0][0][0]) {
1341  POSTGIS_RT_DEBUG(3, "keepnodata = 1 AND source raster pixel is NODATA");
1342  *nodata = 1;
1343  return 1;
1344  }
1345 
1346  for (i = arg->rasters - 1, j = funcarg->ngv - 1; i > 0; i--, j--) {
1347  POSTGIS_RT_DEBUGF(4, "checking raster %d", i);
1348 
1349  /* mask is NODATA */
1350  if (arg->nodata[i][0][0])
1351  continue;
1352  /* mask is NOT NODATA */
1353  else {
1354  POSTGIS_RT_DEBUGF(4, "Using information from geometry %d", j);
1355 
1356  if (funcarg->gv[j].pixval.nodata)
1357  *nodata = 1;
1358  else
1359  *value = funcarg->gv[j].pixval.value;
1360 
1361  return 1;
1362  }
1363  }
1364 
1365  POSTGIS_RT_DEBUG(4, "Using information from source raster");
1366 
1367  /* still here */
1368  if (arg->nodata[0][0][0])
1369  *nodata = 1;
1370  else
1371  *value = arg->values[0][0][0];
1372 
1373  return 1;
1374 }
1375 
1376 PG_FUNCTION_INFO_V1(RASTER_setPixelValuesGeomval);
1377 Datum RASTER_setPixelValuesGeomval(PG_FUNCTION_ARGS)
1378 {
1379  rt_pgraster *pgraster = NULL;
1380  rt_pgraster *pgrtn = NULL;
1381  rt_raster raster = NULL;
1382  rt_band band = NULL;
1383  rt_raster _raster = NULL;
1384  rt_band _band = NULL;
1385  int nband = 0; /* 1-based */
1386 
1387  int numbands = 0;
1388  int width = 0;
1389  int height = 0;
1390  int32_t srid = 0;
1391  double gt[6] = {0};
1392 
1393  rt_pixtype pixtype = PT_END;
1394  int hasnodata = 0;
1395  double nodataval = 0;
1396 
1397  rtpg_setvaluesgv_arg arg = NULL;
1398  int allpoint = 0;
1399 
1400  ArrayType *array;
1401  Oid etype;
1402  Datum *e;
1403  bool *nulls;
1404  int16 typlen;
1405  bool typbyval;
1406  char typalign;
1407  int n = 0;
1408 
1409  HeapTupleHeader tup;
1410  bool isnull;
1411  Datum tupv;
1412 
1413  GSERIALIZED *gser = NULL;
1414  uint8_t gtype;
1415  lwvarlena_t *wkb = NULL;
1416 
1417  int i = 0;
1418  uint32_t j = 0;
1419  int noerr = 1;
1420 
1421  /* pgraster is null, return null */
1422  if (PG_ARGISNULL(0))
1423  PG_RETURN_NULL();
1424  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
1425 
1426  /* raster */
1427  raster = rt_raster_deserialize(pgraster, FALSE);
1428  if (!raster) {
1429  PG_FREE_IF_COPY(pgraster, 0);
1430  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not deserialize raster");
1431  PG_RETURN_NULL();
1432  }
1433 
1434  /* raster attributes */
1435  numbands = rt_raster_get_num_bands(raster);
1436  width = rt_raster_get_width(raster);
1437  height = rt_raster_get_height(raster);
1438  srid = clamp_srid(rt_raster_get_srid(raster));
1439  rt_raster_get_geotransform_matrix(raster, gt);
1440 
1441  /* nband */
1442  if (PG_ARGISNULL(1)) {
1443  elog(NOTICE, "Band index cannot be NULL. Value must be 1-based. Returning original raster");
1444  rt_raster_destroy(raster);
1445  PG_RETURN_POINTER(pgraster);
1446  }
1447 
1448  nband = PG_GETARG_INT32(1);
1449  if (nband < 1 || nband > numbands) {
1450  elog(NOTICE, "Band index is invalid. Value must be 1-based. Returning original raster");
1451  rt_raster_destroy(raster);
1452  PG_RETURN_POINTER(pgraster);
1453  }
1454 
1455  /* get band attributes */
1456  band = rt_raster_get_band(raster, nband - 1);
1457  pixtype = rt_band_get_pixtype(band);
1458  hasnodata = rt_band_get_hasnodata_flag(band);
1459  if (hasnodata)
1460  rt_band_get_nodata(band, &nodataval);
1461 
1462  /* array of geomval (2) */
1463  if (PG_ARGISNULL(2)) {
1464  elog(NOTICE, "No values to set. Returning original raster");
1465  rt_raster_destroy(raster);
1466  PG_RETURN_POINTER(pgraster);
1467  }
1468 
1469  array = PG_GETARG_ARRAYTYPE_P(2);
1470  etype = ARR_ELEMTYPE(array);
1471  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1472 
1473  deconstruct_array(
1474  array,
1475  etype,
1476  typlen, typbyval, typalign,
1477  &e, &nulls, &n
1478  );
1479 
1480  if (!n) {
1481  elog(NOTICE, "No values to set. Returning original raster");
1482  rt_raster_destroy(raster);
1483  PG_RETURN_POINTER(pgraster);
1484  }
1485 
1486  /* init arg */
1487  arg = rtpg_setvaluesgv_arg_init();
1488  if (arg == NULL) {
1489  rt_raster_destroy(raster);
1490  PG_FREE_IF_COPY(pgraster, 0);
1491  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not intialize argument structure");
1492  PG_RETURN_NULL();
1493  }
1494 
1495  arg->gv = palloc(sizeof(struct rtpg_setvaluesgv_geomval_t) * n);
1496  if (arg->gv == NULL) {
1497  rtpg_setvaluesgv_arg_destroy(arg);
1498  rt_raster_destroy(raster);
1499  PG_FREE_IF_COPY(pgraster, 0);
1500  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not allocate memory for geomval array");
1501  PG_RETURN_NULL();
1502  }
1503 
1504  /* process each element */
1505  arg->ngv = 0;
1506  for (i = 0; i < n; i++) {
1507  if (nulls[i])
1508  continue;
1509 
1510  arg->gv[arg->ngv].pixval.nodata = 0;
1511  arg->gv[arg->ngv].pixval.value = 0;
1512  arg->gv[arg->ngv].geom = NULL;
1513  arg->gv[arg->ngv].mask = NULL;
1514 
1515  /* each element is a tuple */
1516  tup = (HeapTupleHeader) DatumGetPointer(e[i]);
1517  if (NULL == tup) {
1518  rtpg_setvaluesgv_arg_destroy(arg);
1519  rt_raster_destroy(raster);
1520  PG_FREE_IF_COPY(pgraster, 0);
1521  elog(ERROR, "RASTER_setPixelValuesGeomval: Invalid argument for geomval at index %d", i);
1522  PG_RETURN_NULL();
1523  }
1524 
1525  /* first element, geometry */
1526  POSTGIS_RT_DEBUG(4, "Processing first element (geometry)");
1527  tupv = GetAttributeByName(tup, "geom", &isnull);
1528  if (isnull) {
1529  elog(NOTICE, "First argument (geom) of geomval at index %d is NULL. Skipping", i);
1530  continue;
1531  }
1532 
1533  gser = (GSERIALIZED *) PG_DETOAST_DATUM(tupv);
1534  arg->gv[arg->ngv].geom = lwgeom_from_gserialized(gser);
1535  if (arg->gv[arg->ngv].geom == NULL) {
1536  rtpg_setvaluesgv_arg_destroy(arg);
1537  rt_raster_destroy(raster);
1538  PG_FREE_IF_COPY(pgraster, 0);
1539  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not deserialize geometry of geomval at index %d", i);
1540  PG_RETURN_NULL();
1541  }
1542 
1543  /* empty geometry */
1544  if (lwgeom_is_empty(arg->gv[arg->ngv].geom)) {
1545  elog(NOTICE, "First argument (geom) of geomval at index %d is an empty geometry. Skipping", i);
1546  continue;
1547  }
1548 
1549  /* check SRID */
1550  if (clamp_srid(gserialized_get_srid(gser)) != srid) {
1551  elog(NOTICE, "Geometry provided for geomval at index %d does not have the same SRID as the raster: %d. Returning original raster", i, srid);
1552  rtpg_setvaluesgv_arg_destroy(arg);
1553  rt_raster_destroy(raster);
1554  PG_RETURN_POINTER(pgraster);
1555  }
1556 
1557  /* Get a 2D version of the geometry if necessary */
1558  if (lwgeom_ndims(arg->gv[arg->ngv].geom) > 2) {
1559  LWGEOM *geom2d = lwgeom_force_2d(arg->gv[arg->ngv].geom);
1560  lwgeom_free(arg->gv[arg->ngv].geom);
1561  arg->gv[arg->ngv].geom = geom2d;
1562  }
1563 
1564  /* filter for types */
1565  gtype = gserialized_get_type(gser);
1566 
1567  /* shortcuts for POINT and MULTIPOINT */
1568  if (gtype == POINTTYPE || gtype == MULTIPOINTTYPE)
1569  allpoint++;
1570 
1571  /* get wkb of geometry */
1572  POSTGIS_RT_DEBUG(3, "getting wkb of geometry");
1573  wkb = lwgeom_to_wkb_varlena(arg->gv[arg->ngv].geom, WKB_SFSQL);
1574 
1575  /* rasterize geometry */
1576  arg->gv[arg->ngv].mask = rt_raster_gdal_rasterize((unsigned char *)wkb->data,
1577  LWSIZE_GET(wkb->size) - LWVARHDRSZ,
1578  NULL,
1579  0,
1580  NULL,
1581  NULL,
1582  NULL,
1583  NULL,
1584  NULL,
1585  NULL,
1586  NULL,
1587  &(gt[1]),
1588  &(gt[5]),
1589  NULL,
1590  NULL,
1591  &(gt[0]),
1592  &(gt[3]),
1593  &(gt[2]),
1594  &(gt[4]),
1595  NULL);
1596 
1597  pfree(wkb);
1598  if (gtype != POINTTYPE && gtype != MULTIPOINTTYPE) {
1599  lwgeom_free(arg->gv[arg->ngv].geom);
1600  arg->gv[arg->ngv].geom = NULL;
1601  }
1602 
1603  if (arg->gv[arg->ngv].mask == NULL) {
1604  rtpg_setvaluesgv_arg_destroy(arg);
1605  rt_raster_destroy(raster);
1606  PG_FREE_IF_COPY(pgraster, 0);
1607  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not rasterize geometry of geomval at index %d", i);
1608  PG_RETURN_NULL();
1609  }
1610 
1611  /* set SRID */
1612  rt_raster_set_srid(arg->gv[arg->ngv].mask, srid);
1613 
1614  /* second element, value */
1615  POSTGIS_RT_DEBUG(4, "Processing second element (val)");
1616  tupv = GetAttributeByName(tup, "val", &isnull);
1617  if (isnull) {
1618  elog(NOTICE, "Second argument (val) of geomval at index %d is NULL. Treating as NODATA", i);
1619  arg->gv[arg->ngv].pixval.nodata = 1;
1620  }
1621  else
1622  arg->gv[arg->ngv].pixval.value = DatumGetFloat8(tupv);
1623 
1624  (arg->ngv)++;
1625  }
1626 
1627  /* redim arg->gv if needed */
1628  if (arg->ngv < n) {
1629  arg->gv = repalloc(arg->gv, sizeof(struct rtpg_setvaluesgv_geomval_t) * arg->ngv);
1630  if (arg->gv == NULL) {
1631  rtpg_setvaluesgv_arg_destroy(arg);
1632  rt_raster_destroy(raster);
1633  PG_FREE_IF_COPY(pgraster, 0);
1634  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not reallocate memory for geomval array");
1635  PG_RETURN_NULL();
1636  }
1637  }
1638 
1639  /* keepnodata */
1640  if (!PG_ARGISNULL(3))
1641  arg->keepnodata = PG_GETARG_BOOL(3);
1642  POSTGIS_RT_DEBUGF(3, "keepnodata = %d", arg->keepnodata);
1643 
1644  /* keepnodata = TRUE and band is NODATA */
1645  if (arg->keepnodata && rt_band_get_isnodata_flag(band)) {
1646  POSTGIS_RT_DEBUG(3, "keepnodata = TRUE and band is NODATA. Not doing anything");
1647  }
1648  /* all elements are points */
1649  else if (allpoint == arg->ngv) {
1650  double igt[6] = {0};
1651  double xy[2] = {0};
1652  double value = 0;
1653  int isnodata = 0;
1654 
1655  LWCOLLECTION *coll = NULL;
1656  LWPOINT *point = NULL;
1657  POINT2D p;
1658 
1659  POSTGIS_RT_DEBUG(3, "all geometries are points, using direct to pixel method");
1660 
1661  /* cache inverse gretransform matrix */
1662  rt_raster_get_inverse_geotransform_matrix(NULL, gt, igt);
1663 
1664  /* process each geometry */
1665  for (i = 0; i < arg->ngv; i++) {
1666  /* convert geometry to collection */
1667  coll = lwgeom_as_lwcollection(lwgeom_as_multi(arg->gv[i].geom));
1668 
1669  /* process each point in collection */
1670  for (j = 0; j < coll->ngeoms; j++) {
1671  point = lwgeom_as_lwpoint(coll->geoms[j]);
1672  getPoint2d_p(point->point, 0, &p);
1673 
1674  if (rt_raster_geopoint_to_cell(raster, p.x, p.y, &(xy[0]), &(xy[1]), igt) != ES_NONE) {
1675  rtpg_setvaluesgv_arg_destroy(arg);
1676  rt_raster_destroy(raster);
1677  PG_FREE_IF_COPY(pgraster, 0);
1678  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not process coordinates of point");
1679  PG_RETURN_NULL();
1680  }
1681 
1682  /* skip point if outside raster */
1683  if (
1684  (xy[0] < 0 || xy[0] >= width) ||
1685  (xy[1] < 0 || xy[1] >= height)
1686  ) {
1687  elog(NOTICE, "Point is outside raster extent. Skipping");
1688  continue;
1689  }
1690 
1691  /* get pixel value */
1692  if (rt_band_get_pixel(band, xy[0], xy[1], &value, &isnodata) != ES_NONE) {
1693  rtpg_setvaluesgv_arg_destroy(arg);
1694  rt_raster_destroy(raster);
1695  PG_FREE_IF_COPY(pgraster, 0);
1696  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not get pixel value");
1697  PG_RETURN_NULL();
1698  }
1699 
1700  /* keepnodata = TRUE AND pixel value is NODATA */
1701  if (arg->keepnodata && isnodata)
1702  continue;
1703 
1704  /* set pixel */
1705  if (arg->gv[i].pixval.nodata)
1706  noerr = rt_band_set_pixel(band, xy[0], xy[1], nodataval, NULL);
1707  else
1708  noerr = rt_band_set_pixel(band, xy[0], xy[1], arg->gv[i].pixval.value, NULL);
1709 
1710  if (noerr != ES_NONE) {
1711  rtpg_setvaluesgv_arg_destroy(arg);
1712  rt_raster_destroy(raster);
1713  PG_FREE_IF_COPY(pgraster, 0);
1714  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not set pixel value");
1715  PG_RETURN_NULL();
1716  }
1717  }
1718  }
1719  }
1720  /* run iterator otherwise */
1721  else {
1722  rt_iterator itrset;
1723 
1724  POSTGIS_RT_DEBUG(3, "a mix of geometries, using iterator method");
1725 
1726  /* init itrset */
1727  itrset = palloc(sizeof(struct rt_iterator_t) * (arg->ngv + 1));
1728  if (itrset == NULL) {
1729  rtpg_setvaluesgv_arg_destroy(arg);
1730  rt_raster_destroy(raster);
1731  PG_FREE_IF_COPY(pgraster, 0);
1732  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not allocate memory for iterator arguments");
1733  PG_RETURN_NULL();
1734  }
1735 
1736  /* set first raster's info */
1737  itrset[0].raster = raster;
1738  itrset[0].nband = nband - 1;
1739  itrset[0].nbnodata = 1;
1740 
1741  /* set other raster's info */
1742  for (i = 0, j = 1; i < arg->ngv; i++, j++) {
1743  itrset[j].raster = arg->gv[i].mask;
1744  itrset[j].nband = 0;
1745  itrset[j].nbnodata = 1;
1746  }
1747 
1748  /* pass to iterator */
1749  noerr = rt_raster_iterator(
1750  itrset, arg->ngv + 1,
1751  ET_FIRST, NULL,
1752  pixtype,
1753  hasnodata, nodataval,
1754  0, 0,
1755  NULL,
1756  arg,
1757  rtpg_setvalues_geomval_callback,
1758  &_raster
1759  );
1760  pfree(itrset);
1761 
1762  if (noerr != ES_NONE) {
1763  rtpg_setvaluesgv_arg_destroy(arg);
1764  rt_raster_destroy(raster);
1765  PG_FREE_IF_COPY(pgraster, 0);
1766  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not run raster iterator function");
1767  PG_RETURN_NULL();
1768  }
1769 
1770  /* copy band from _raster to raster */
1771  _band = rt_raster_get_band(_raster, 0);
1772  if (_band == NULL) {
1773  rtpg_setvaluesgv_arg_destroy(arg);
1774  rt_raster_destroy(_raster);
1775  rt_raster_destroy(raster);
1776  PG_FREE_IF_COPY(pgraster, 0);
1777  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not get band from working raster");
1778  PG_RETURN_NULL();
1779  }
1780 
1781  _band = rt_raster_replace_band(raster, _band, nband - 1);
1782  if (_band == NULL) {
1783  rtpg_setvaluesgv_arg_destroy(arg);
1784  rt_raster_destroy(_raster);
1785  rt_raster_destroy(raster);
1786  PG_FREE_IF_COPY(pgraster, 0);
1787  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not replace band in output raster");
1788  PG_RETURN_NULL();
1789  }
1790 
1791  rt_band_destroy(_band);
1792  rt_raster_destroy(_raster);
1793  }
1794 
1795  rtpg_setvaluesgv_arg_destroy(arg);
1796 
1797  pgrtn = rt_raster_serialize(raster);
1798  rt_raster_destroy(raster);
1799  PG_FREE_IF_COPY(pgraster, 0);
1800 
1801  POSTGIS_RT_DEBUG(3, "Finished");
1802 
1803  if (!pgrtn)
1804  PG_RETURN_NULL();
1805 
1806  SET_VARSIZE(pgrtn, pgrtn->size);
1807  PG_RETURN_POINTER(pgrtn);
1808 }
1809 
1810 #undef VALUES_LENGTH
1811 #define VALUES_LENGTH 3
1812 
1816 PG_FUNCTION_INFO_V1(RASTER_pixelOfValue);
1817 Datum RASTER_pixelOfValue(PG_FUNCTION_ARGS)
1818 {
1819  FuncCallContext *funcctx;
1820  TupleDesc tupdesc;
1821 
1822  rt_pixel pixels = NULL;
1823  rt_pixel pixels2 = NULL;
1824  int count = 0;
1825  int i = 0;
1826  int n = 0;
1827  int call_cntr;
1828  int max_calls;
1829 
1830  if (SRF_IS_FIRSTCALL()) {
1831  MemoryContext oldcontext;
1832 
1833  rt_pgraster *pgraster = NULL;
1834  rt_raster raster = NULL;
1835  rt_band band = NULL;
1836  int nband = 1;
1837  int num_bands = 0;
1838  double *search = NULL;
1839  int nsearch = 0;
1840  double val;
1841  bool exclude_nodata_value = TRUE;
1842 
1843  ArrayType *array;
1844  Oid etype;
1845  Datum *e;
1846  bool *nulls;
1847  int16 typlen;
1848  bool typbyval;
1849  char typalign;
1850 
1851  /* create a function context for cross-call persistence */
1852  funcctx = SRF_FIRSTCALL_INIT();
1853 
1854  /* switch to memory context appropriate for multiple function calls */
1855  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1856 
1857  if (PG_ARGISNULL(0)) {
1858  MemoryContextSwitchTo(oldcontext);
1859  SRF_RETURN_DONE(funcctx);
1860  }
1861  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1862  raster = rt_raster_deserialize(pgraster, FALSE);
1863  if (!raster) {
1864  PG_FREE_IF_COPY(pgraster, 0);
1865  MemoryContextSwitchTo(oldcontext);
1866  elog(ERROR, "RASTER_pixelOfValue: Could not deserialize raster");
1867  SRF_RETURN_DONE(funcctx);
1868  }
1869 
1870  /* num_bands */
1871  num_bands = rt_raster_get_num_bands(raster);
1872  if (num_bands < 1) {
1873  elog(NOTICE, "Raster provided has no bands");
1874  rt_raster_destroy(raster);
1875  PG_FREE_IF_COPY(pgraster, 0);
1876  MemoryContextSwitchTo(oldcontext);
1877  SRF_RETURN_DONE(funcctx);
1878  }
1879 
1880  /* band index is 1-based */
1881  if (!PG_ARGISNULL(1))
1882  nband = PG_GETARG_INT32(1);
1883  if (nband < 1 || nband > num_bands) {
1884  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1885  rt_raster_destroy(raster);
1886  PG_FREE_IF_COPY(pgraster, 0);
1887  MemoryContextSwitchTo(oldcontext);
1888  SRF_RETURN_DONE(funcctx);
1889  }
1890 
1891  /* search values */
1892  array = PG_GETARG_ARRAYTYPE_P(2);
1893  etype = ARR_ELEMTYPE(array);
1894  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1895 
1896  switch (etype) {
1897  case FLOAT4OID:
1898  case FLOAT8OID:
1899  break;
1900  default:
1901  rt_raster_destroy(raster);
1902  PG_FREE_IF_COPY(pgraster, 0);
1903  MemoryContextSwitchTo(oldcontext);
1904  elog(ERROR, "RASTER_pixelOfValue: Invalid data type for pixel values");
1905  SRF_RETURN_DONE(funcctx);
1906  break;
1907  }
1908 
1909  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1910  &nulls, &n);
1911 
1912  search = palloc(sizeof(double) * n);
1913  for (i = 0, nsearch = 0; i < n; i++) {
1914  if (nulls[i]) continue;
1915 
1916  switch (etype) {
1917  case FLOAT4OID:
1918  val = (double) DatumGetFloat4(e[i]);
1919  break;
1920  case FLOAT8OID:
1921  val = (double) DatumGetFloat8(e[i]);
1922  break;
1923  }
1924 
1925  search[nsearch] = val;
1926  POSTGIS_RT_DEBUGF(3, "search[%d] = %f", nsearch, search[nsearch]);
1927  nsearch++;
1928  }
1929 
1930  /* not searching for anything */
1931  if (nsearch < 1) {
1932  elog(NOTICE, "No search values provided. Returning NULL");
1933  pfree(search);
1934  rt_raster_destroy(raster);
1935  PG_FREE_IF_COPY(pgraster, 0);
1936  MemoryContextSwitchTo(oldcontext);
1937  SRF_RETURN_DONE(funcctx);
1938  }
1939  else if (nsearch < n)
1940  search = repalloc(search, sizeof(double) * nsearch);
1941 
1942  /* exclude_nodata_value flag */
1943  if (!PG_ARGISNULL(3))
1944  exclude_nodata_value = PG_GETARG_BOOL(3);
1945 
1946  /* get band */
1947  band = rt_raster_get_band(raster, nband - 1);
1948  if (!band) {
1949  elog(NOTICE, "Could not find band at index %d. Returning NULL", nband);
1950  rt_raster_destroy(raster);
1951  PG_FREE_IF_COPY(pgraster, 0);
1952  MemoryContextSwitchTo(oldcontext);
1953  SRF_RETURN_DONE(funcctx);
1954  }
1955 
1956  /* get pixels of values */
1957  count = rt_band_get_pixel_of_value(
1958  band, exclude_nodata_value,
1959  search, nsearch,
1960  &pixels
1961  );
1962  pfree(search);
1963  rt_band_destroy(band);
1964  rt_raster_destroy(raster);
1965  PG_FREE_IF_COPY(pgraster, 0);
1966  if (count < 1) {
1967  /* error */
1968  if (count < 0)
1969  elog(NOTICE, "Could not get the pixels of search values for band at index %d", nband);
1970  /* no nearest pixel */
1971  else
1972  elog(NOTICE, "No pixels of search values found for band at index %d", nband);
1973 
1974  MemoryContextSwitchTo(oldcontext);
1975  SRF_RETURN_DONE(funcctx);
1976  }
1977 
1978  /* Store needed information */
1979  funcctx->user_fctx = pixels;
1980 
1981  /* total number of tuples to be returned */
1982  funcctx->max_calls = count;
1983 
1984  /* Build a tuple descriptor for our result type */
1985  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1986  ereport(ERROR, (
1987  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1988  errmsg(
1989  "function returning record called in context "
1990  "that cannot accept type record"
1991  )
1992  ));
1993  }
1994 
1995  BlessTupleDesc(tupdesc);
1996  funcctx->tuple_desc = tupdesc;
1997 
1998  MemoryContextSwitchTo(oldcontext);
1999  }
2000 
2001  /* stuff done on every call of the function */
2002  funcctx = SRF_PERCALL_SETUP();
2003 
2004  call_cntr = funcctx->call_cntr;
2005  max_calls = funcctx->max_calls;
2006  tupdesc = funcctx->tuple_desc;
2007  pixels2 = funcctx->user_fctx;
2008 
2009  /* do when there is more left to send */
2010  if (call_cntr < max_calls) {
2011  Datum values[VALUES_LENGTH];
2012  bool nulls[VALUES_LENGTH];
2013  HeapTuple tuple;
2014  Datum result;
2015 
2016  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
2017 
2018  /* 0-based to 1-based */
2019  pixels2[call_cntr].x += 1;
2020  pixels2[call_cntr].y += 1;
2021 
2022  values[0] = Float8GetDatum(pixels2[call_cntr].value);
2023  values[1] = Int32GetDatum(pixels2[call_cntr].x);
2024  values[2] = Int32GetDatum(pixels2[call_cntr].y);
2025 
2026  /* build a tuple */
2027  tuple = heap_form_tuple(tupdesc, values, nulls);
2028 
2029  /* make the tuple into a datum */
2030  result = HeapTupleGetDatum(tuple);
2031 
2032  SRF_RETURN_NEXT(funcctx, result);
2033  }
2034  else {
2035  pfree(pixels2);
2036  SRF_RETURN_DONE(funcctx);
2037  }
2038 }
2039 
2043 PG_FUNCTION_INFO_V1(RASTER_nearestValue);
2044 Datum RASTER_nearestValue(PG_FUNCTION_ARGS)
2045 {
2046  rt_pgraster *pgraster = NULL;
2047  rt_raster raster = NULL;
2048  rt_band band = NULL;
2049  int bandindex = 1;
2050  int num_bands = 0;
2051  GSERIALIZED *geom;
2052  bool exclude_nodata_value = TRUE;
2053  LWGEOM *lwgeom;
2054  LWPOINT *point = NULL;
2055  POINT2D p;
2056 
2057  double x;
2058  double y;
2059  int count;
2060  rt_pixel npixels = NULL;
2061  double value = 0;
2062  int hasvalue = 0;
2063  int isnodata = 0;
2064 
2065  if (PG_ARGISNULL(0))
2066  PG_RETURN_NULL();
2067  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2068  raster = rt_raster_deserialize(pgraster, FALSE);
2069  if (!raster) {
2070  PG_FREE_IF_COPY(pgraster, 0);
2071  elog(ERROR, "RASTER_nearestValue: Could not deserialize raster");
2072  PG_RETURN_NULL();
2073  }
2074 
2075  /* band index is 1-based */
2076  if (!PG_ARGISNULL(1))
2077  bandindex = PG_GETARG_INT32(1);
2078  num_bands = rt_raster_get_num_bands(raster);
2079  if (bandindex < 1 || bandindex > num_bands) {
2080  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2081  rt_raster_destroy(raster);
2082  PG_FREE_IF_COPY(pgraster, 0);
2083  PG_RETURN_NULL();
2084  }
2085 
2086  /* point */
2087  geom = PG_GETARG_GSERIALIZED_P(2);
2088  if (gserialized_get_type(geom) != POINTTYPE) {
2089  elog(NOTICE, "Geometry provided must be a point");
2090  rt_raster_destroy(raster);
2091  PG_FREE_IF_COPY(pgraster, 0);
2092  PG_FREE_IF_COPY(geom, 2);
2093  PG_RETURN_NULL();
2094  }
2095 
2096  /* exclude_nodata_value flag */
2097  if (!PG_ARGISNULL(3))
2098  exclude_nodata_value = PG_GETARG_BOOL(3);
2099 
2100  /* SRIDs of raster and geometry must match */
2101  if (clamp_srid(gserialized_get_srid(geom)) != clamp_srid(rt_raster_get_srid(raster))) {
2102  elog(NOTICE, "SRIDs of geometry and raster do not match");
2103  rt_raster_destroy(raster);
2104  PG_FREE_IF_COPY(pgraster, 0);
2105  PG_FREE_IF_COPY(geom, 2);
2106  PG_RETURN_NULL();
2107  }
2108 
2109  /* get band */
2110  band = rt_raster_get_band(raster, bandindex - 1);
2111  if (!band) {
2112  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
2113  rt_raster_destroy(raster);
2114  PG_FREE_IF_COPY(pgraster, 0);
2115  PG_FREE_IF_COPY(geom, 2);
2116  PG_RETURN_NULL();
2117  }
2118 
2119  /* process geometry */
2120  lwgeom = lwgeom_from_gserialized(geom);
2121 
2122  if (lwgeom_is_empty(lwgeom)) {
2123  elog(NOTICE, "Geometry provided cannot be empty");
2124  rt_raster_destroy(raster);
2125  PG_FREE_IF_COPY(pgraster, 0);
2126  PG_FREE_IF_COPY(geom, 2);
2127  PG_RETURN_NULL();
2128  }
2129 
2130  /* Get a 2D version of the geometry if necessary */
2131  if (lwgeom_ndims(lwgeom) > 2) {
2132  LWGEOM *lwgeom2d = lwgeom_force_2d(lwgeom);
2133  lwgeom_free(lwgeom);
2134  lwgeom = lwgeom2d;
2135  }
2136 
2137  point = lwgeom_as_lwpoint(lwgeom);
2138  getPoint2d_p(point->point, 0, &p);
2139 
2140  if (rt_raster_geopoint_to_cell(
2141  raster,
2142  p.x, p.y,
2143  &x, &y,
2144  NULL
2145  ) != ES_NONE) {
2146  rt_raster_destroy(raster);
2147  PG_FREE_IF_COPY(pgraster, 0);
2148  lwgeom_free(lwgeom);
2149  PG_FREE_IF_COPY(geom, 2);
2150  elog(ERROR, "RASTER_nearestValue: Could not compute pixel coordinates from spatial coordinates");
2151  PG_RETURN_NULL();
2152  }
2153 
2154  /* get value at point */
2155  if (
2156  (x >= 0 && x < rt_raster_get_width(raster)) &&
2157  (y >= 0 && y < rt_raster_get_height(raster))
2158  ) {
2159  if (rt_band_get_pixel(band, x, y, &value, &isnodata) != ES_NONE) {
2160  rt_raster_destroy(raster);
2161  PG_FREE_IF_COPY(pgraster, 0);
2162  lwgeom_free(lwgeom);
2163  PG_FREE_IF_COPY(geom, 2);
2164  elog(ERROR, "RASTER_nearestValue: Could not get pixel value for band at index %d", bandindex);
2165  PG_RETURN_NULL();
2166  }
2167 
2168  /* value at point, return value */
2169  if (!exclude_nodata_value || !isnodata) {
2170  rt_raster_destroy(raster);
2171  PG_FREE_IF_COPY(pgraster, 0);
2172  lwgeom_free(lwgeom);
2173  PG_FREE_IF_COPY(geom, 2);
2174 
2175  PG_RETURN_FLOAT8(value);
2176  }
2177  }
2178 
2179  /* get neighborhood */
2180  count = rt_band_get_nearest_pixel(
2181  band,
2182  x, y,
2183  0, 0,
2184  exclude_nodata_value,
2185  &npixels
2186  );
2187  rt_band_destroy(band);
2188  /* error or no neighbors */
2189  if (count < 1) {
2190  /* error */
2191  if (count < 0)
2192  elog(NOTICE, "Could not get the nearest value for band at index %d", bandindex);
2193  /* no nearest pixel */
2194  else
2195  elog(NOTICE, "No nearest value found for band at index %d", bandindex);
2196 
2197  lwgeom_free(lwgeom);
2198  PG_FREE_IF_COPY(geom, 2);
2199  rt_raster_destroy(raster);
2200  PG_FREE_IF_COPY(pgraster, 0);
2201  PG_RETURN_NULL();
2202  }
2203 
2204  /* more than one nearest value, see which one is closest */
2205  if (count > 1) {
2206  int i = 0;
2207  LWPOLY *poly = NULL;
2208  double lastdist = -1;
2209  double dist;
2210 
2211  for (i = 0; i < count; i++) {
2212  /* convex-hull of pixel */
2213  poly = rt_raster_pixel_as_polygon(raster, npixels[i].x, npixels[i].y);
2214  if (!poly) {
2215  lwgeom_free(lwgeom);
2216  PG_FREE_IF_COPY(geom, 2);
2217  rt_raster_destroy(raster);
2218  PG_FREE_IF_COPY(pgraster, 0);
2219  elog(ERROR, "RASTER_nearestValue: Could not get polygon of neighboring pixel");
2220  PG_RETURN_NULL();
2221  }
2222 
2223  /* distance between convex-hull and point */
2224  dist = lwgeom_mindistance2d(lwpoly_as_lwgeom(poly), lwgeom);
2225  if (lastdist < 0 || dist < lastdist) {
2226  value = npixels[i].value;
2227  hasvalue = 1;
2228  }
2229  lastdist = dist;
2230 
2231  lwpoly_free(poly);
2232  }
2233  }
2234  else {
2235  value = npixels[0].value;
2236  hasvalue = 1;
2237  }
2238 
2239  pfree(npixels);
2240  lwgeom_free(lwgeom);
2241  PG_FREE_IF_COPY(geom, 2);
2242  rt_raster_destroy(raster);
2243  PG_FREE_IF_COPY(pgraster, 0);
2244 
2245  if (hasvalue)
2246  PG_RETURN_FLOAT8(value);
2247  else
2248  PG_RETURN_NULL();
2249 }
2250 
2254 PG_FUNCTION_INFO_V1(RASTER_neighborhood);
2255 Datum RASTER_neighborhood(PG_FUNCTION_ARGS)
2256 {
2257  rt_pgraster *pgraster = NULL;
2258  rt_raster raster = NULL;
2259  rt_band band = NULL;
2260  int bandindex = 1;
2261  int num_bands = 0;
2262  int x = 0;
2263  int y = 0;
2264  int _x = 0;
2265  int _y = 0;
2266  int distance[2] = {0};
2267  bool exclude_nodata_value = TRUE;
2268  double pixval;
2269  int isnodata = 0;
2270 
2271  rt_pixel npixels = NULL;
2272  int count;
2273  double **value2D = NULL;
2274  int **nodata2D = NULL;
2275 
2276  int i = 0;
2277  int j = 0;
2278  int k = 0;
2279  Datum *value1D = NULL;
2280  bool *nodata1D = NULL;
2281  int dim[2] = {0};
2282  int lbound[2] = {1, 1};
2283  ArrayType *mdArray = NULL;
2284 
2285  int16 typlen;
2286  bool typbyval;
2287  char typalign;
2288 
2289  /* pgraster is null, return nothing */
2290  if (PG_ARGISNULL(0))
2291  PG_RETURN_NULL();
2292  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2293 
2294  raster = rt_raster_deserialize(pgraster, FALSE);
2295  if (!raster) {
2296  PG_FREE_IF_COPY(pgraster, 0);
2297  elog(ERROR, "RASTER_neighborhood: Could not deserialize raster");
2298  PG_RETURN_NULL();
2299  }
2300 
2301  /* band index is 1-based */
2302  if (!PG_ARGISNULL(1))
2303  bandindex = PG_GETARG_INT32(1);
2304  num_bands = rt_raster_get_num_bands(raster);
2305  if (bandindex < 1 || bandindex > num_bands) {
2306  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2307  rt_raster_destroy(raster);
2308  PG_FREE_IF_COPY(pgraster, 0);
2309  PG_RETURN_NULL();
2310  }
2311 
2312  /* pixel column, 1-based */
2313  x = PG_GETARG_INT32(2);
2314  _x = x - 1;
2315 
2316  /* pixel row, 1-based */
2317  y = PG_GETARG_INT32(3);
2318  _y = y - 1;
2319 
2320  /* distance X axis */
2321  distance[0] = PG_GETARG_INT32(4);
2322  if (distance[0] < 0) {
2323  elog(NOTICE, "Invalid value for distancex (must be >= zero). Returning NULL");
2324  rt_raster_destroy(raster);
2325  PG_FREE_IF_COPY(pgraster, 0);
2326  PG_RETURN_NULL();
2327  }
2328  distance[0] = (uint16_t) distance[0];
2329 
2330  /* distance Y axis */
2331  distance[1] = PG_GETARG_INT32(5);
2332  if (distance[1] < 0) {
2333  elog(NOTICE, "Invalid value for distancey (must be >= zero). Returning NULL");
2334  rt_raster_destroy(raster);
2335  PG_FREE_IF_COPY(pgraster, 0);
2336  PG_RETURN_NULL();
2337  }
2338  distance[1] = (uint16_t) distance[1];
2339 
2340  /* exclude_nodata_value flag */
2341  if (!PG_ARGISNULL(6))
2342  exclude_nodata_value = PG_GETARG_BOOL(6);
2343 
2344  /* get band */
2345  band = rt_raster_get_band(raster, bandindex - 1);
2346  if (!band) {
2347  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
2348  rt_raster_destroy(raster);
2349  PG_FREE_IF_COPY(pgraster, 0);
2350  PG_RETURN_NULL();
2351  }
2352 
2353  /* get neighborhood */
2354  count = 0;
2355  npixels = NULL;
2356  if (distance[0] > 0 || distance[1] > 0) {
2357  count = rt_band_get_nearest_pixel(
2358  band,
2359  _x, _y,
2360  distance[0], distance[1],
2361  exclude_nodata_value,
2362  &npixels
2363  );
2364  /* error */
2365  if (count < 0) {
2366  elog(NOTICE, "Could not get the pixel's neighborhood for band at index %d", bandindex);
2367 
2368  rt_band_destroy(band);
2369  rt_raster_destroy(raster);
2370  PG_FREE_IF_COPY(pgraster, 0);
2371 
2372  PG_RETURN_NULL();
2373  }
2374  }
2375 
2376  /* get pixel's value */
2377  if (
2378  (_x >= 0 && _x < rt_band_get_width(band)) &&
2379  (_y >= 0 && _y < rt_band_get_height(band))
2380  ) {
2381  if (rt_band_get_pixel(
2382  band,
2383  _x, _y,
2384  &pixval,
2385  &isnodata
2386  ) != ES_NONE) {
2387  elog(NOTICE, "Could not get the pixel of band at index %d. Returning NULL", bandindex);
2388  rt_band_destroy(band);
2389  rt_raster_destroy(raster);
2390  PG_FREE_IF_COPY(pgraster, 0);
2391  PG_RETURN_NULL();
2392  }
2393  }
2394  /* outside band extent, set to NODATA */
2395  else {
2396  /* has NODATA, use NODATA */
2397  if (rt_band_get_hasnodata_flag(band))
2398  rt_band_get_nodata(band, &pixval);
2399  /* no NODATA, use min possible value */
2400  else
2401  pixval = rt_band_get_min_value(band);
2402  isnodata = 1;
2403  }
2404  POSTGIS_RT_DEBUGF(4, "pixval: %f", pixval);
2405 
2406 
2407  /* add pixel to neighborhood */
2408  count++;
2409  if (count > 1)
2410  npixels = (rt_pixel) repalloc(npixels, sizeof(struct rt_pixel_t) * count);
2411  else
2412  npixels = (rt_pixel) palloc(sizeof(struct rt_pixel_t));
2413  if (npixels == NULL) {
2414 
2415  rt_band_destroy(band);
2416  rt_raster_destroy(raster);
2417  PG_FREE_IF_COPY(pgraster, 0);
2418 
2419  elog(ERROR, "RASTER_neighborhood: Could not reallocate memory for neighborhood");
2420  PG_RETURN_NULL();
2421  }
2422  npixels[count - 1].x = _x;
2423  npixels[count - 1].y = _y;
2424  npixels[count - 1].nodata = 1;
2425  npixels[count - 1].value = pixval;
2426 
2427  /* set NODATA */
2428  if (!exclude_nodata_value || !isnodata) {
2429  npixels[count - 1].nodata = 0;
2430  }
2431 
2432  /* free unnecessary stuff */
2433  rt_band_destroy(band);
2434  rt_raster_destroy(raster);
2435  PG_FREE_IF_COPY(pgraster, 0);
2436 
2437  /* convert set of rt_pixel to 2D array */
2438  /* dim is passed with element 0 being Y-axis and element 1 being X-axis */
2439  count = rt_pixel_set_to_array(
2440  npixels, count, NULL,
2441  _x, _y,
2442  distance[0], distance[1],
2443  &value2D,
2444  &nodata2D,
2445  &(dim[1]), &(dim[0])
2446  );
2447  pfree(npixels);
2448  if (count != ES_NONE) {
2449  elog(NOTICE, "Could not create 2D array of neighborhood");
2450  PG_RETURN_NULL();
2451  }
2452 
2453  /* 1D arrays for values and nodata from 2D arrays */
2454  value1D = palloc(sizeof(Datum) * dim[0] * dim[1]);
2455  nodata1D = palloc(sizeof(bool) * dim[0] * dim[1]);
2456 
2457  if (value1D == NULL || nodata1D == NULL) {
2458 
2459  for (i = 0; i < dim[0]; i++) {
2460  pfree(value2D[i]);
2461  pfree(nodata2D[i]);
2462  }
2463  pfree(value2D);
2464  pfree(nodata2D);
2465 
2466  elog(ERROR, "RASTER_neighborhood: Could not allocate memory for return 2D array");
2467  PG_RETURN_NULL();
2468  }
2469 
2470  /* copy values from 2D arrays to 1D arrays */
2471  k = 0;
2472  /* Y-axis */
2473  for (i = 0; i < dim[0]; i++) {
2474  /* X-axis */
2475  for (j = 0; j < dim[1]; j++) {
2476  nodata1D[k] = (bool) nodata2D[i][j];
2477  if (!nodata1D[k])
2478  value1D[k] = Float8GetDatum(value2D[i][j]);
2479  else
2480  value1D[k] = PointerGetDatum(NULL);
2481 
2482  k++;
2483  }
2484  }
2485 
2486  /* no more need for 2D arrays */
2487  for (i = 0; i < dim[0]; i++) {
2488  pfree(value2D[i]);
2489  pfree(nodata2D[i]);
2490  }
2491  pfree(value2D);
2492  pfree(nodata2D);
2493 
2494  /* info about the type of item in the multi-dimensional array (float8). */
2495  get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
2496 
2497  mdArray = construct_md_array(
2498  value1D, nodata1D,
2499  2, dim, lbound,
2500  FLOAT8OID,
2501  typlen, typbyval, typalign
2502  );
2503 
2504  pfree(value1D);
2505  pfree(nodata1D);
2506 
2507  PG_RETURN_ARRAYTYPE_P(mdArray);
2508 }
2509 
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition: cu_print.c:262
TRUE
#define TRUE
Definition: dbfopen.c:73
FALSE
#define FALSE
Definition: dbfopen.c:72
gserialized_get_srid
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_from_gserialized
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
Definition: gserialized.c:239
gserialized_is_empty
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition: gserialized.c:152
gserialized_get_type
uint32_t gserialized_get_type(const GSERIALIZED *g)
Extract the geometry type from the serialized form (it hides in the anonymous data area,...
Definition: gserialized.c:89
gserialized_from_lwgeom
GSERIALIZED * gserialized_from_lwgeom(LWGEOM *geom, size_t *size)
Allocate a new GSERIALIZED from an LWGEOM.
Definition: gserialized.c:222
lwgeom_ndims
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition: lwgeom.c:955
LWVARHDRSZ
#define LWVARHDRSZ
Definition: liblwgeom.h:311
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1155
lwgeom_as_multi
LWGEOM * lwgeom_as_multi(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate MULTI* type.
Definition: lwgeom.c:380
LWSIZE_GET
#define LWSIZE_GET(varsize)
Macro for reading the size from the GSERIALIZED size attribute.
Definition: liblwgeom.h:324
lwpoly_as_lwgeom
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:329
MULTIPOINTTYPE
#define MULTIPOINTTYPE
Definition: liblwgeom.h:105
lwpoint_get_x
double lwpoint_get_x(const LWPOINT *point)
Definition: lwpoint.c:63
getPoint2d_p
int getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point)
Definition: lwgeom_api.c:342
POINTTYPE
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:102
lwgeom_mindistance2d
double lwgeom_mindistance2d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initializing min distance calculation.
Definition: measures.c:197
lwgeom_as_lwcollection
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:233
lwpoly_free
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:175
lwgeom_to_wkb_varlena
lwvarlena_t * lwgeom_to_wkb_varlena(const LWGEOM *geom, uint8_t variant)
Definition: lwout_wkb.c:851
WKB_SFSQL
#define WKB_SFSQL
Definition: liblwgeom.h:2176
clamp_srid
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
lwpoint_get_y
double lwpoint_get_y(const LWPOINT *point)
Definition: lwpoint.c:76
lwgeom_force_2d
LWGEOM * lwgeom_force_2d(const LWGEOM *geom)
Strip out the Z/M components of an LWGEOM.
Definition: lwgeom.c:793
rt_band_get_pixel_of_value
int rt_band_get_pixel_of_value(rt_band band, int exclude_nodata_value, double *searchset, int searchcount, rt_pixel *pixels)
Search band for pixel(s) with search values.
Definition: rt_band.c:1809
rt_band_get_width
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition: rt_band.c:640
rt_raster_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:360
rt_raster_get_inverse_geotransform_matrix
rt_errorstate rt_raster_get_inverse_geotransform_matrix(rt_raster raster, double *gt, double *igt)
Get 6-element array of raster inverse geotransform matrix.
Definition: rt_raster.c:680
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:674
rt_raster_geopoint_to_cell
rt_errorstate rt_raster_geopoint_to_cell(rt_raster raster, double xw, double yw, double *xr, double *yr, double *igt)
Convert an xw, yw map point to a xr, yr raster point.
Definition: rt_raster.c:808
rt_band_get_pixel
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
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:86
rt_raster_pixel_as_polygon
LWPOLY * rt_raster_pixel_as_polygon(rt_raster raster, int x, int y)
Get a raster pixel as a polygon.
Definition: rt_geometry.c:607
rt_pixtype
rt_pixtype
Definition: librtcore.h:187
PT_END
@ PT_END
Definition: librtcore.h:199
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:714
rt_raster_geopoint_to_rasterpoint
rt_errorstate rt_raster_geopoint_to_rasterpoint(rt_raster raster, double xw, double yw, double *xr, double *yr, double *igt)
Convert an xw,yw map point to a xr,yr raster point.
Definition: rt_raster.c:852
rt_raster_has_band
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
rt_pixel
struct rt_pixel_t * rt_pixel
Definition: librtcore.h:149
rt_raster_serialize
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition: rt_serialize.c:521
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2387
rt_resample_type
rt_resample_type
Definition: librtcore.h:1450
RT_BILINEAR
@ RT_BILINEAR
Definition: librtcore.h:1452
RT_NEAREST
@ RT_NEAREST
Definition: librtcore.h:1451
rt_raster_replace_band
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
rt_band_get_min_value
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_band_set_pixel
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
rt_errorstate
rt_errorstate
Enum definitions.
Definition: librtcore.h:181
ES_NONE
@ ES_NONE
Definition: librtcore.h:182
rt_raster_gdal_rasterize
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
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:340
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:376
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:133
rt_raster_iterator
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.
Definition: rt_mapalgebra.c:812
rt_raster_set_srid
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_raster.c:367
ET_FIRST
@ ET_FIRST
Definition: librtcore.h:205
rt_band_get_pixel_resample
rt_errorstate rt_band_get_pixel_resample(rt_band band, double xr, double yr, rt_resample_type resample, double *r_value, int *r_nodata)
Retrieve a point value from the raster using a world coordinate and selected resampling method.
Definition: rt_band.c:1221
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1887
rt_band_get_pixtype
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:631
rt_pixel_set_to_array
rt_errorstate rt_pixel_set_to_array(rt_pixel npixel, uint32_t count, rt_mask mask, int x, int y, uint16_t distancex, uint16_t distancey, double ***value, int ***nodata, int *dimx, int *dimy)
Definition: rt_pixel.c:288
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:125
rt_band_get_nearest_pixel
uint32_t rt_band_get_nearest_pixel(rt_band band, int x, int y, uint16_t distancex, uint16_t distancey, int exclude_nodata_value, rt_pixel *npixels)
Get nearest pixel(s) with value (not NODATA) to specified pixel.
Definition: rt_band.c:1529
rt_band_clamped_value_is_nodata
int rt_band_clamped_value_is_nodata(rt_band band, double val)
Compare clamped value to band's clamped NODATA value.
Definition: rt_band.c:1955
rt_raster_get_geotransform_matrix
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:710
rt_raster_copy_to_geometry
rt_errorstate rt_raster_copy_to_geometry(rt_raster raster, uint32_t bandnum, char dim, rt_resample_type resample, const LWGEOM *lwgeom_in, LWGEOM **lwgeom_out)
Copy values from a raster to the points on a geometry using the requested interpolation type.
Definition: rt_raster.c:1614
rt_band_get_height
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition: rt_band.c:649
rt_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:725
rt_raster_is_empty
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1362
rt_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:385
lwgeom_is_empty
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
lwgeom_as_lwpoint
static LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwinline.h:131
distance
static double distance(double x1, double y1, double x2, double y2)
Definition: lwtree.c:1032
genraster.value
int value
Definition: genraster.py:62
genraster.count
int count
Definition: genraster.py:57
ovdump.band
band
Definition: ovdump.py:58
pixval.pixval
pixval
Definition: pixval.py:94
pixval.nband
nband
Definition: pixval.py:53
pixval
Definition: pixval.py:1
rtrowdump.raster
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
window.gt
gt
Definition: window.py:78
rtpg_dumpvalues_arg
struct rtpg_dumpvalues_arg_t * rtpg_dumpvalues_arg
Definition: rtpg_pixel.c:341
RASTER_getPixelValueResample
Datum RASTER_getPixelValueResample(PG_FUNCTION_ARGS)
Definition: rtpg_pixel.c:165
RASTER_setPixelValuesArray
Datum RASTER_setPixelValuesArray(PG_FUNCTION_ARGS)
Definition: rtpg_pixel.c:864
rtpg_dumpvalues_arg_destroy
static void rtpg_dumpvalues_arg_destroy(rtpg_dumpvalues_arg arg)
Definition: rtpg_pixel.c:372
rtpg_setvaluesgv_arg_destroy
static void rtpg_setvaluesgv_arg_destroy(rtpg_setvaluesgv_arg arg)
Definition: rtpg_pixel.c:1309
RASTER_dumpValues
Datum RASTER_dumpValues(PG_FUNCTION_ARGS)
Definition: rtpg_pixel.c:402
RASTER_getPixelValue
Datum RASTER_getPixelValue(PG_FUNCTION_ARGS)
Definition: rtpg_pixel.c:77
resample_text_to_type
static rt_resample_type resample_text_to_type(text *txt)
Definition: rtpg_pixel.c:142
rtpg_dumpvalues_arg_init
static rtpg_dumpvalues_arg rtpg_dumpvalues_arg_init()
Definition: rtpg_pixel.c:352
RASTER_nearestValue
Datum RASTER_nearestValue(PG_FUNCTION_ARGS)
Definition: rtpg_pixel.c:2044
RASTER_getGeometryValues
Datum RASTER_getGeometryValues(PG_FUNCTION_ARGS)
Definition: rtpg_pixel.c:258
RASTER_pixelOfValue
Datum RASTER_pixelOfValue(PG_FUNCTION_ARGS)
Definition: rtpg_pixel.c:1817
RASTER_neighborhood
Datum RASTER_neighborhood(PG_FUNCTION_ARGS)
Definition: rtpg_pixel.c:2255
rtpg_setvaluesgv_geomval
struct rtpg_setvaluesgv_geomval_t * rtpg_setvaluesgv_geomval
Definition: rtpg_pixel.c:1276
rtpg_setvalues_geomval_callback
static int rtpg_setvalues_geomval_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
Definition: rtpg_pixel.c:1326
rtpg_setvaluesgv_arg_init
static rtpg_setvaluesgv_arg rtpg_setvaluesgv_arg_init()
Definition: rtpg_pixel.c:1295
VALUES_LENGTH
#define VALUES_LENGTH
Definition: rtpg_pixel.c:1811
rtpg_setvaluesgv_arg
struct rtpg_setvaluesgv_arg_t * rtpg_setvaluesgv_arg
Definition: rtpg_pixel.c:1275
PG_FUNCTION_INFO_V1
PG_FUNCTION_INFO_V1(RASTER_getPixelValue)
Return value of a single pixel.
RASTER_setPixelValuesGeomval
Datum RASTER_setPixelValuesGeomval(PG_FUNCTION_ARGS)
Definition: rtpg_pixel.c:1377
RASTER_setPixelValue
Datum RASTER_setPixelValue(PG_FUNCTION_ARGS)
Definition: rtpg_pixel.c:767
POSTGIS_RT_DEBUG
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rtpostgis.h:65
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:69
LWCOLLECTION::ngeoms
uint32_t ngeoms
Definition: liblwgeom.h:580
LWCOLLECTION::geoms
LWGEOM ** geoms
Definition: liblwgeom.h:575
LWPOINT::point
POINTARRAY * point
Definition: liblwgeom.h:471
POINT2D::y
double y
Definition: liblwgeom.h:390
POINT2D::x
double x
Definition: liblwgeom.h:390
lwvarlena_t::size
uint32_t size
Definition: liblwgeom.h:307
lwvarlena_t::data
char data[]
Definition: liblwgeom.h:308
lwvarlena_t
Definition: liblwgeom.h:306
rt_band_t
Definition: librtcore.h:2466
rt_iterator_arg_t::values
double *** values
Definition: librtcore.h:2612
rt_iterator_arg_t::nodata
int *** nodata
Definition: librtcore.h:2614
rt_iterator_arg_t::rasters
uint16_t rasters
Definition: librtcore.h:2604
rt_iterator_arg_t
Definition: librtcore.h:2602
rt_iterator_t::raster
rt_raster raster
Definition: librtcore.h:2596
rt_iterator_t::nband
uint16_t nband
Definition: librtcore.h:2597
rt_iterator_t::nbnodata
uint8_t nbnodata
Definition: librtcore.h:2598
rt_iterator_t
Definition: librtcore.h:2595
rt_pixel_t::y
int y
Definition: librtcore.h:2488
rt_pixel_t::value
double value
Definition: librtcore.h:2491
rt_pixel_t::nodata
uint8_t nodata
Definition: librtcore.h:2490
rt_pixel_t::x
int x
Definition: librtcore.h:2487
rt_pixel_t
Definition: librtcore.h:2486
rt_raster_serialized_t::size
uint32_t size
Definition: librtcore.h:2405
rt_raster_serialized_t
Struct definitions.
Definition: librtcore.h:2403
rt_raster_t
Definition: librtcore.h:2437
rtpg_dumpvalues_arg_t::columns
int columns
Definition: rtpg_pixel.c:345
rtpg_dumpvalues_arg_t::rows
int rows
Definition: rtpg_pixel.c:344
rtpg_dumpvalues_arg_t::values
Datum ** values
Definition: rtpg_pixel.c:348
rtpg_dumpvalues_arg_t::nbands
int * nbands
Definition: rtpg_pixel.c:347
rtpg_dumpvalues_arg_t::nodata
bool ** nodata
Definition: rtpg_pixel.c:349
rtpg_dumpvalues_arg_t::numbands
int numbands
Definition: rtpg_pixel.c:343
rtpg_dumpvalues_arg_t
Definition: rtpg_pixel.c:342
rtpg_setvaluesgv_arg_t::ngv
int ngv
Definition: rtpg_pixel.c:1279
rtpg_setvaluesgv_arg_t::gv
rtpg_setvaluesgv_geomval gv
Definition: rtpg_pixel.c:1280
rtpg_setvaluesgv_arg_t::keepnodata
bool keepnodata
Definition: rtpg_pixel.c:1282
rtpg_setvaluesgv_arg_t
Definition: rtpg_pixel.c:1278
rtpg_setvaluesgv_geomval_t::mask
rt_raster mask
Definition: rtpg_pixel.c:1292
rtpg_setvaluesgv_geomval_t::value
double value
Definition: rtpg_pixel.c:1288
rtpg_setvaluesgv_geomval_t::nodata
int nodata
Definition: rtpg_pixel.c:1287
rtpg_setvaluesgv_geomval_t::geom
LWGEOM * geom
Definition: rtpg_pixel.c:1291
rtpg_setvaluesgv_geomval_t::pixval
struct rtpg_setvaluesgv_geomval_t::@27 pixval
rtpg_setvaluesgv_geomval_t
Definition: rtpg_pixel.c:1285