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