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