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