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