PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches
rtpg_geometry.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 <funcapi.h>
33#include <utils/lsyscache.h> /* for get_typlenbyvalalign */
34#include <utils/array.h> /* for ArrayType */
35#include <catalog/pg_type.h> /* for INT2OID, INT4OID, FLOAT4OID, FLOAT8OID and TEXTOID */
36#include <utils/builtins.h> /* for text_to_cstring() */
37
38#include "../../postgis_config.h"
39
40#include "access/htup_details.h" /* for heap_form_tuple() */
41#include "lwgeom_pg.h"
42
43#include "rtpostgis.h"
44#include "rtpg_internal.h"
45
46Datum RASTER_envelope(PG_FUNCTION_ARGS);
47Datum RASTER_convex_hull(PG_FUNCTION_ARGS);
48Datum RASTER_dumpAsPolygons(PG_FUNCTION_ARGS);
49
50/* Get pixel geographical shape */
51Datum RASTER_getPixelPolygons(PG_FUNCTION_ARGS);
52
53/* Get pixel centroid points */
54Datum RASTER_getPixelCentroids(PG_FUNCTION_ARGS);
55
56/* Get raster band's polygon */
57Datum RASTER_getPolygon(PG_FUNCTION_ARGS);
58
59/* rasterize a geometry */
60Datum RASTER_asRaster(PG_FUNCTION_ARGS);
61
62/* ---------------------------------------------------------------- */
63/* Raster envelope */
64/* ---------------------------------------------------------------- */
65PG_FUNCTION_INFO_V1(RASTER_envelope);
66Datum RASTER_envelope(PG_FUNCTION_ARGS)
67{
68 rt_pgraster *pgraster;
69 rt_raster raster;
70 LWGEOM *geom = NULL;
71 GSERIALIZED* gser = NULL;
72 size_t gser_size;
73 int err = ES_NONE;
74
75 if (PG_ARGISNULL(0))
76 PG_RETURN_NULL();
77
78 pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(
79 PG_GETARG_DATUM(0),
80 0,
81 sizeof(struct rt_raster_serialized_t)
82 );
83 raster = rt_raster_deserialize(pgraster, TRUE);
84
85 if (!raster) {
86 PG_FREE_IF_COPY(pgraster, 0);
87 elog(ERROR, "RASTER_envelope: Could not deserialize raster");
88 PG_RETURN_NULL();
89 }
90
91 err = rt_raster_get_envelope_geom(raster, &geom);
92
93 rt_raster_destroy(raster);
94 PG_FREE_IF_COPY(pgraster, 0);
95
96 if (err != ES_NONE) {
97 elog(ERROR, "RASTER_envelope: Could not get raster's envelope");
98 PG_RETURN_NULL();
99 }
100 else if (geom == NULL) {
101 elog(NOTICE, "Raster's envelope is NULL");
102 PG_RETURN_NULL();
103 }
104
105 gser = gserialized_from_lwgeom(geom, &gser_size);
106 lwgeom_free(geom);
107
108 SET_VARSIZE(gser, gser_size);
109 PG_RETURN_POINTER(gser);
110}
111
115/* ---------------------------------------------------------------- */
116/* Raster convex hull */
117/* ---------------------------------------------------------------- */
118PG_FUNCTION_INFO_V1(RASTER_convex_hull);
119Datum RASTER_convex_hull(PG_FUNCTION_ARGS)
120{
121 rt_pgraster *pgraster;
122 rt_raster raster;
123 LWGEOM *geom = NULL;
124 GSERIALIZED* gser = NULL;
125 size_t gser_size;
126 int err = ES_NONE;
127
128 bool minhull = FALSE;
129
130 if (PG_ARGISNULL(0))
131 PG_RETURN_NULL();
132
133 /* # of args */
134 if (PG_NARGS() > 1)
135 minhull = TRUE;
136
137 if (!minhull) {
138 pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
139 raster = rt_raster_deserialize(pgraster, TRUE);
140 }
141 else {
142 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
143 raster = rt_raster_deserialize(pgraster, FALSE);
144 }
145
146 if (!raster) {
147 PG_FREE_IF_COPY(pgraster, 0);
148 elog(ERROR, "RASTER_convex_hull: Could not deserialize raster");
149 PG_RETURN_NULL();
150 }
151
152 if (!minhull)
153 err = rt_raster_get_convex_hull(raster, &geom);
154 else {
155 int nband = -1;
156
157 /* get arg 1 */
158 if (!PG_ARGISNULL(1)) {
159 nband = PG_GETARG_INT32(1);
160 if (!rt_raster_has_band(raster, nband - 1)) {
161 elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
162 rt_raster_destroy(raster);
163 PG_FREE_IF_COPY(pgraster, 0);
164 PG_RETURN_NULL();
165 }
166 nband = nband - 1;
167 }
168
169 err = rt_raster_get_perimeter(raster, nband, &geom);
170 }
171
172 rt_raster_destroy(raster);
173 PG_FREE_IF_COPY(pgraster, 0);
174
175 if (err != ES_NONE) {
176 elog(ERROR, "RASTER_convex_hull: Could not get raster's convex hull");
177 PG_RETURN_NULL();
178 }
179 else if (geom == NULL) {
180 elog(NOTICE, "Raster's convex hull is NULL");
181 PG_RETURN_NULL();
182 }
183
184 gser = gserialized_from_lwgeom(geom, &gser_size);
185 lwgeom_free(geom);
186
187 SET_VARSIZE(gser, gser_size);
188 PG_RETURN_POINTER(gser);
189}
190
191#define VALUES_LENGTH 2
192
193PG_FUNCTION_INFO_V1(RASTER_dumpAsPolygons);
194Datum RASTER_dumpAsPolygons(PG_FUNCTION_ARGS) {
195 FuncCallContext *funcctx;
196 TupleDesc tupdesc;
197 rt_geomval geomval;
198 rt_geomval geomval2;
199 int call_cntr;
200 int max_calls;
201
202 /* stuff done only on the first call of the function */
203 if (SRF_IS_FIRSTCALL()) {
204 MemoryContext oldcontext;
205 int numbands;
206 rt_pgraster *pgraster = NULL;
207 rt_raster raster = NULL;
208 int nband;
209 bool exclude_nodata_value = TRUE;
210 int nElements;
211
212 POSTGIS_RT_DEBUG(2, "RASTER_dumpAsPolygons first call");
213
214 /* create a function context for cross-call persistence */
215 funcctx = SRF_FIRSTCALL_INIT();
216
217 /* switch to memory context appropriate for multiple function calls */
218 oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
219
220 /* Get input arguments */
221 if (PG_ARGISNULL(0)) {
222 MemoryContextSwitchTo(oldcontext);
223 SRF_RETURN_DONE(funcctx);
224 }
225 pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
226 raster = rt_raster_deserialize(pgraster, FALSE);
227 if (!raster) {
228 PG_FREE_IF_COPY(pgraster, 0);
229 ereport(ERROR, (
230 errcode(ERRCODE_OUT_OF_MEMORY),
231 errmsg("Could not deserialize raster")
232 ));
233 MemoryContextSwitchTo(oldcontext);
234 SRF_RETURN_DONE(funcctx);
235 }
236
237 if (!PG_ARGISNULL(1))
238 nband = PG_GETARG_UINT32(1);
239 else
240 nband = 1; /* By default, first band */
241
242 POSTGIS_RT_DEBUGF(3, "band %d", nband);
243 numbands = rt_raster_get_num_bands(raster);
244
245 if (nband < 1 || nband > numbands) {
246 elog(NOTICE, "Invalid band index (must use 1-based). Returning empty set");
247 rt_raster_destroy(raster);
248 PG_FREE_IF_COPY(pgraster, 0);
249 MemoryContextSwitchTo(oldcontext);
250 SRF_RETURN_DONE(funcctx);
251 }
252
253 if (!PG_ARGISNULL(2))
254 exclude_nodata_value = PG_GETARG_BOOL(2);
255
256 /* see if band is NODATA */
257 if (rt_band_get_isnodata_flag(rt_raster_get_band(raster, nband - 1))) {
258 POSTGIS_RT_DEBUGF(3, "Band at index %d is NODATA. Returning empty set", nband);
259 rt_raster_destroy(raster);
260 PG_FREE_IF_COPY(pgraster, 0);
261 MemoryContextSwitchTo(oldcontext);
262 SRF_RETURN_DONE(funcctx);
263 }
264
265 /* Polygonize raster */
266
270 geomval = rt_raster_gdal_polygonize(raster, nband - 1, exclude_nodata_value, &nElements);
271 rt_raster_destroy(raster);
272 PG_FREE_IF_COPY(pgraster, 0);
273 if (NULL == geomval) {
274 ereport(ERROR, (
275 errcode(ERRCODE_NO_DATA_FOUND),
276 errmsg("Could not polygonize raster")
277 ));
278 MemoryContextSwitchTo(oldcontext);
279 SRF_RETURN_DONE(funcctx);
280 }
281
282 POSTGIS_RT_DEBUGF(3, "raster dump, %d elements returned", nElements);
283
284 /* Store needed information */
285 funcctx->user_fctx = geomval;
286
287 /* total number of tuples to be returned */
288 funcctx->max_calls = nElements;
289
290 /* Build a tuple descriptor for our result type */
291 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
292 ereport(ERROR, (
293 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
294 errmsg("function returning record called in context that cannot accept type record")
295 ));
296 }
297
298 BlessTupleDesc(tupdesc);
299 funcctx->tuple_desc = tupdesc;
300
301 MemoryContextSwitchTo(oldcontext);
302 }
303
304 /* stuff done on every call of the function */
305 funcctx = SRF_PERCALL_SETUP();
306
307 call_cntr = funcctx->call_cntr;
308 max_calls = funcctx->max_calls;
309 tupdesc = funcctx->tuple_desc;
310 geomval2 = funcctx->user_fctx;
311
312 /* do when there is more left to send */
313 if (call_cntr < max_calls) {
314 Datum values[VALUES_LENGTH];
315 bool nulls[VALUES_LENGTH];
316 HeapTuple tuple;
317 Datum result;
318
319 GSERIALIZED *gser = NULL;
320 size_t gser_size = 0;
321
322 POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr);
323
324 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
325
326 /* convert LWGEOM to GSERIALIZED */
327 gser = gserialized_from_lwgeom(lwpoly_as_lwgeom(geomval2[call_cntr].geom), &gser_size);
328 lwgeom_free(lwpoly_as_lwgeom(geomval2[call_cntr].geom));
329
330 values[0] = PointerGetDatum(gser);
331 values[1] = Float8GetDatum(geomval2[call_cntr].val);
332
333 /* build a tuple */
334 tuple = heap_form_tuple(tupdesc, values, nulls);
335
336 /* make the tuple into a datum */
337 result = HeapTupleGetDatum(tuple);
338
339 SRF_RETURN_NEXT(funcctx, result);
340 }
341 /* do when there is no more left */
342 else {
343 pfree(geomval2);
344 SRF_RETURN_DONE(funcctx);
345 }
346}
347
348#undef VALUES_LENGTH
349#define VALUES_LENGTH 4
350
354PG_FUNCTION_INFO_V1(RASTER_getPixelPolygons);
355Datum RASTER_getPixelPolygons(PG_FUNCTION_ARGS)
356{
357 FuncCallContext *funcctx;
358 TupleDesc tupdesc;
359 rt_pixel pix = NULL;
360 rt_pixel pix2;
361 int call_cntr;
362 int max_calls;
363 int i = 0;
364
365 /* stuff done only on the first call of the function */
366 if (SRF_IS_FIRSTCALL()) {
367 MemoryContext oldcontext;
368 rt_pgraster *pgraster = NULL;
369 rt_raster raster = NULL;
370 rt_band band = NULL;
371 int nband = 1;
372 int numbands;
373 bool hasband = TRUE;
374 bool exclude_nodata_value = TRUE;
375 bool nocolumnx = FALSE;
376 bool norowy = FALSE;
377 int x = 0;
378 int y = 0;
379 int bounds[4] = {0};
380 int pixcount = 0;
381 double value = 0;
382 int isnodata = 0;
383
384 LWPOLY *poly;
385
386 POSTGIS_RT_DEBUG(3, "RASTER_getPixelPolygons first call");
387
388 /* create a function context for cross-call persistence */
389 funcctx = SRF_FIRSTCALL_INIT();
390
391 /* switch to memory context appropriate for multiple function calls */
392 oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
393
394 if (PG_ARGISNULL(0)) {
395 MemoryContextSwitchTo(oldcontext);
396 SRF_RETURN_DONE(funcctx);
397 }
398 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
399
400 /* band */
401 if (PG_ARGISNULL(1))
402 hasband = FALSE;
403 else {
404 nband = PG_GETARG_INT32(1);
405 hasband = TRUE;
406 }
407
408 /* column */
409 if (PG_ARGISNULL(2))
410 nocolumnx = TRUE;
411 else {
412 bounds[0] = PG_GETARG_INT32(2);
413 bounds[1] = bounds[0];
414 }
415
416 /* row */
417 if (PG_ARGISNULL(3))
418 norowy = TRUE;
419 else {
420 bounds[2] = PG_GETARG_INT32(3);
421 bounds[3] = bounds[2];
422 }
423
424 /* exclude NODATA */
425 if (!PG_ARGISNULL(4))
426 exclude_nodata_value = PG_GETARG_BOOL(4);
427
428 raster = rt_raster_deserialize(pgraster, FALSE);
429 if (!raster) {
430 PG_FREE_IF_COPY(pgraster, 0);
431 ereport(ERROR, (
432 errcode(ERRCODE_OUT_OF_MEMORY),
433 errmsg("Could not deserialize raster")
434 ));
435 MemoryContextSwitchTo(oldcontext);
436 SRF_RETURN_DONE(funcctx);
437 }
438
439 /* raster empty, return empty set */
440 if (rt_raster_is_empty(raster)) {
441 elog(NOTICE, "Raster is empty. Returning empty set");
442 rt_raster_destroy(raster);
443 PG_FREE_IF_COPY(pgraster, 0);
444 MemoryContextSwitchTo(oldcontext);
445 SRF_RETURN_DONE(funcctx);
446 }
447
448 /* band specified, load band and info */
449 if (hasband) {
450 numbands = rt_raster_get_num_bands(raster);
451 POSTGIS_RT_DEBUGF(3, "band %d", nband);
452 POSTGIS_RT_DEBUGF(3, "# of bands %d", numbands);
453
454 if (nband < 1 || nband > numbands) {
455 elog(NOTICE, "Invalid band index (must use 1-based). Returning empty set");
456 rt_raster_destroy(raster);
457 PG_FREE_IF_COPY(pgraster, 0);
458 MemoryContextSwitchTo(oldcontext);
459 SRF_RETURN_DONE(funcctx);
460 }
461
462 band = rt_raster_get_band(raster, nband - 1);
463 if (!band) {
464 elog(NOTICE, "Could not find band at index %d. Returning empty set", nband);
465 rt_raster_destroy(raster);
466 PG_FREE_IF_COPY(pgraster, 0);
467 MemoryContextSwitchTo(oldcontext);
468 SRF_RETURN_DONE(funcctx);
469 }
470
471 if (!rt_band_get_hasnodata_flag(band))
472 exclude_nodata_value = FALSE;
473 }
474
475 /* set bounds if columnx, rowy not set */
476 if (nocolumnx) {
477 bounds[0] = 1;
478 bounds[1] = rt_raster_get_width(raster);
479 }
480 if (norowy) {
481 bounds[2] = 1;
482 bounds[3] = rt_raster_get_height(raster);
483 }
484 POSTGIS_RT_DEBUGF(3, "bounds (min x, max x, min y, max y) = (%d, %d, %d, %d)",
485 bounds[0], bounds[1], bounds[2], bounds[3]);
486
487 /* rowy */
488 pixcount = 0;
489 for (y = bounds[2]; y <= bounds[3]; y++) {
490 /* columnx */
491 for (x = bounds[0]; x <= bounds[1]; x++) {
492
493 value = 0;
494 isnodata = TRUE;
495
496 if (hasband) {
497 if (rt_band_get_pixel(band, x - 1, y - 1, &value, &isnodata) != ES_NONE) {
498
499 for (i = 0; i < pixcount; i++)
500 lwgeom_free(pix[i].geom);
501 if (pixcount) pfree(pix);
502
503 rt_band_destroy(band);
504 rt_raster_destroy(raster);
505 PG_FREE_IF_COPY(pgraster, 0);
506
507 MemoryContextSwitchTo(oldcontext);
508 elog(ERROR, "RASTER_getPixelPolygons: Could not get pixel value");
509 SRF_RETURN_DONE(funcctx);
510 }
511
512 /* don't continue if pixel is NODATA and to exclude NODATA */
513 if (isnodata && exclude_nodata_value) {
514 POSTGIS_RT_DEBUG(5, "pixel value is NODATA and exclude_nodata_value = TRUE");
515 continue;
516 }
517 }
518
519 /* geometry */
520 poly = rt_raster_pixel_as_polygon(raster, x - 1, y - 1);
521 if (!poly) {
522 for (i = 0; i < pixcount; i++)
523 lwgeom_free(pix[i].geom);
524 if (pixcount) pfree(pix);
525
526 if (hasband) rt_band_destroy(band);
527 rt_raster_destroy(raster);
528 PG_FREE_IF_COPY(pgraster, 0);
529
530 MemoryContextSwitchTo(oldcontext);
531 elog(ERROR, "RASTER_getPixelPolygons: Could not get pixel polygon");
532 SRF_RETURN_DONE(funcctx);
533 }
534
535 if (!pixcount)
536 pix = palloc(sizeof(struct rt_pixel_t) * (pixcount + 1));
537 else
538 pix = repalloc(pix, sizeof(struct rt_pixel_t) * (pixcount + 1));
539 if (pix == NULL) {
540
541 lwpoly_free(poly);
542 if (hasband) rt_band_destroy(band);
543 rt_raster_destroy(raster);
544 PG_FREE_IF_COPY(pgraster, 0);
545
546 MemoryContextSwitchTo(oldcontext);
547 elog(ERROR, "RASTER_getPixelPolygons: Could not allocate memory for storing pixel polygons");
548 SRF_RETURN_DONE(funcctx);
549 }
550 pix[pixcount].geom = (LWGEOM *) poly;
551 POSTGIS_RT_DEBUGF(5, "poly @ %p", poly);
552 POSTGIS_RT_DEBUGF(5, "geom @ %p", pix[pixcount].geom);
553
554 /* x, y */
555 pix[pixcount].x = x;
556 pix[pixcount].y = y;
557
558 /* value */
559 pix[pixcount].value = value;
560
561 /* NODATA */
562 if (hasband) {
563 if (exclude_nodata_value)
564 pix[pixcount].nodata = isnodata;
565 else
566 pix[pixcount].nodata = FALSE;
567 }
568 else {
569 pix[pixcount].nodata = isnodata;
570 }
571
572 pixcount++;
573 }
574 }
575
576 if (hasband) rt_band_destroy(band);
577 rt_raster_destroy(raster);
578 PG_FREE_IF_COPY(pgraster, 0);
579
580 /* shortcut if no pixcount */
581 if (pixcount < 1) {
582 elog(NOTICE, "No pixels found for band %d", nband);
583 MemoryContextSwitchTo(oldcontext);
584 SRF_RETURN_DONE(funcctx);
585 }
586
587 /* Store needed information */
588 funcctx->user_fctx = pix;
589
590 /* total number of tuples to be returned */
591 funcctx->max_calls = pixcount;
592 POSTGIS_RT_DEBUGF(3, "pixcount = %d", pixcount);
593
594 /* Build a tuple descriptor for our result type */
595 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
596 ereport(ERROR, (
597 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
598 errmsg("function returning record called in context that cannot accept type record")
599 ));
600 }
601
602 BlessTupleDesc(tupdesc);
603 funcctx->tuple_desc = tupdesc;
604
605 MemoryContextSwitchTo(oldcontext);
606 }
607
608 /* stuff done on every call of the function */
609 funcctx = SRF_PERCALL_SETUP();
610
611 call_cntr = funcctx->call_cntr;
612 max_calls = funcctx->max_calls;
613 tupdesc = funcctx->tuple_desc;
614 pix2 = funcctx->user_fctx;
615
616 /* do when there is more left to send */
617 if (call_cntr < max_calls) {
618 Datum values[VALUES_LENGTH];
619 bool nulls[VALUES_LENGTH];
620 HeapTuple tuple;
621 Datum result;
622
623 GSERIALIZED *gser = NULL;
624 size_t gser_size = 0;
625
626 POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr);
627
628 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
629
630 /* convert LWGEOM to GSERIALIZED */
631 gser = gserialized_from_lwgeom(pix2[call_cntr].geom, &gser_size);
632 lwgeom_free(pix2[call_cntr].geom);
633
634 values[0] = PointerGetDatum(gser);
635 if (pix2[call_cntr].nodata)
636 nulls[1] = TRUE;
637 else
638 values[1] = Float8GetDatum(pix2[call_cntr].value);
639 values[2] = Int32GetDatum(pix2[call_cntr].x);
640 values[3] = Int32GetDatum(pix2[call_cntr].y);
641
642 /* build a tuple */
643 tuple = heap_form_tuple(tupdesc, values, nulls);
644
645 /* make the tuple into a datum */
646 result = HeapTupleGetDatum(tuple);
647
648 SRF_RETURN_NEXT(funcctx, result);
649 }
650 /* do when there is no more left */
651 else {
652 pfree(pix2);
653 SRF_RETURN_DONE(funcctx);
654 }
655}
656
657#undef VALUES_LENGTH
658#define VALUES_LENGTH 4
659
663PG_FUNCTION_INFO_V1(RASTER_getPixelCentroids);
664Datum RASTER_getPixelCentroids(PG_FUNCTION_ARGS)
665{
666 FuncCallContext *funcctx;
667 TupleDesc tupdesc;
668 rt_pixel pix = NULL;
669 rt_pixel pix2;
670 int call_cntr;
671 int max_calls;
672 int i = 0;
673
674 if (SRF_IS_FIRSTCALL()) {
675 MemoryContext oldcontext;
676 rt_pgraster *pgraster = NULL;
677 rt_raster raster = NULL;
678 rt_band band = NULL;
679 int nband = 1;
680 int numbands;
681 bool hasband = FALSE;
682 bool exclude_nodata_value = TRUE;
683 bool nocolumnx = TRUE;
684 bool norowy = TRUE;
685 int x = 0;
686 int y = 0;
687 int bounds[4] = {0};
688 int pixcount = 0;
689 double value = 0;
690 int isnodata = 0;
691
692 LWPOINT* point = NULL;
693
694 POSTGIS_RT_DEBUG(3, "RASTER_getPixelCentroids first call");
695
696 /* create a function context for cross-call persistence */
697 funcctx = SRF_FIRSTCALL_INIT();
698
699 /* switch to memory context appropriate for multiple function calls */
700 oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
701
702 /* raster */
703 if (PG_ARGISNULL(0)) {
704 MemoryContextSwitchTo(oldcontext);
705 SRF_RETURN_DONE(funcctx);
706 }
707 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
708
709 /* band */
710 if (!PG_ARGISNULL(1)) {
711 nband = PG_GETARG_INT32(1);
712 hasband = TRUE;
713 }
714
715 /* column */
716 if (!PG_ARGISNULL(2)) {
717 bounds[0] = PG_GETARG_INT32(2);
718 bounds[1] = bounds[0];
719 nocolumnx = FALSE;
720 }
721
722 /* row */
723 if (!PG_ARGISNULL(3)) {
724 bounds[2] = PG_GETARG_INT32(3);
725 bounds[3] = bounds[2];
726 norowy = FALSE;
727 }
728
729 /* exclude NODATA */
730 if (!PG_ARGISNULL(4))
731 exclude_nodata_value = PG_GETARG_BOOL(4);
732
733 /* deserialize raster */
734 raster = rt_raster_deserialize(pgraster, FALSE);
735 if (!raster) {
736 PG_FREE_IF_COPY(pgraster, 0);
737 ereport(ERROR, (
738 errcode(ERRCODE_OUT_OF_MEMORY),
739 errmsg("Could not deserialize raster")
740 ));
741 MemoryContextSwitchTo(oldcontext);
742 SRF_RETURN_DONE(funcctx);
743 }
744
745 /* raster empty, return empty set */
746 if (rt_raster_is_empty(raster)) {
747 elog(NOTICE, "Raster is empty. Returning empty set");
748 rt_raster_destroy(raster);
749 PG_FREE_IF_COPY(pgraster, 0);
750 MemoryContextSwitchTo(oldcontext);
751 SRF_RETURN_DONE(funcctx);
752 }
753
754 /* band specified, load band and info */
755 if (hasband) {
756 numbands = rt_raster_get_num_bands(raster);
757 POSTGIS_RT_DEBUGF(3, "band %d", nband);
758 POSTGIS_RT_DEBUGF(3, "# of bands %d", numbands);
759
760 /* check band index */
761 if (nband < 1 || nband > numbands) {
762 elog(NOTICE, "Invalid band index (must use 1-based). Returning empty set");
763 rt_raster_destroy(raster);
764 PG_FREE_IF_COPY(pgraster, 0);
765 MemoryContextSwitchTo(oldcontext);
766 SRF_RETURN_DONE(funcctx);
767 }
768
769 /* get band */
770 band = rt_raster_get_band(raster, nband - 1);
771 if (!band) {
772 elog(NOTICE, "Could not find band at index %d. Returning empty set", nband);
773
774 rt_raster_destroy(raster);
775 PG_FREE_IF_COPY(pgraster, 0);
776
777 MemoryContextSwitchTo(oldcontext);
778 SRF_RETURN_DONE(funcctx);
779 }
780
781 if (!rt_band_get_hasnodata_flag(band))
782 exclude_nodata_value = FALSE;
783 }
784
785 /* set bounds if columnx, rowy not set */
786 if (nocolumnx) {
787 bounds[0] = 1;
788 bounds[1] = rt_raster_get_width(raster);
789 }
790 if (norowy) {
791 bounds[2] = 1;
792 bounds[3] = rt_raster_get_height(raster);
793 }
794 POSTGIS_RT_DEBUGF(3, "bounds (min x, max x, min y, max y) = (%d, %d, %d, %d)",
795 bounds[0], bounds[1], bounds[2], bounds[3]);
796
797 /* rowy */
798 pixcount = 0;
799 for (y = bounds[2]; y <= bounds[3]; y++) {
800 /* columnx */
801 for (x = bounds[0]; x <= bounds[1]; x++) {
802
803 value = 0;
804 isnodata = TRUE;
805
806 if (hasband) {
807 if (rt_band_get_pixel(band, x - 1, y - 1, &value, &isnodata) != ES_NONE) {
808 /* free already created centroid points */
809 for (i = 0; i < pixcount; i++)
810 lwgeom_free(pix[i].geom);
811 if (pixcount) pfree(pix);
812
813 rt_band_destroy(band);
814 rt_raster_destroy(raster);
815 PG_FREE_IF_COPY(pgraster, 0);
816
817 MemoryContextSwitchTo(oldcontext);
818 elog(ERROR, "RASTER_getPixelCentroids: Could not get pixel value");
819 SRF_RETURN_DONE(funcctx);
820 }
821
822 /* don't continue if pixel is NODATA and to exclude NODATA */
823 if (isnodata && exclude_nodata_value) {
824 POSTGIS_RT_DEBUG(5, "pixel value is NODATA and exclude_nodata_value = TRUE");
825 continue;
826 }
827 }
828
829 /* geometry */
830 point = rt_raster_pixel_as_centroid_point(raster, x - 1, y - 1);
831 if (!point) {
832 /* free already created centroid points */
833 for (i = 0; i < pixcount; i++)
834 lwgeom_free(pix[i].geom);
835 if (pixcount) pfree(pix);
836
837 if (hasband) rt_band_destroy(band);
838 rt_raster_destroy(raster);
839 PG_FREE_IF_COPY(pgraster, 0);
840
841 MemoryContextSwitchTo(oldcontext);
842 elog(ERROR, "RASTER_getPixelCentroids: Could not get pixel centroid");
843 SRF_RETURN_DONE(funcctx);
844 }
845
846 /* allocate space for new point */
847 if (!pixcount) {
848 pix = palloc(sizeof(struct rt_pixel_t) * (pixcount + 1));
849 }
850 else {
851 pix = repalloc(pix, sizeof(struct rt_pixel_t) * (pixcount + 1));
852 }
853 if (pix == NULL) {
854 lwpoint_free(point);
855
856 if (hasband) rt_band_destroy(band);
857 rt_raster_destroy(raster);
858 PG_FREE_IF_COPY(pgraster, 0);
859
860 MemoryContextSwitchTo(oldcontext);
861 elog(ERROR, "RASTER_getPixelCentroids: Could not allocate memory for storing pixel centroids");
862 SRF_RETURN_DONE(funcctx);
863 }
864
865 /* set pixel geometry */
866 pix[pixcount].geom = (LWGEOM *) point;
867 POSTGIS_RT_DEBUGF(5, "point @ %p", point);
868 POSTGIS_RT_DEBUGF(5, "geom @ %p", pix[pixcount].geom);
869
870 /* set pixel coordinates */
871 pix[pixcount].x = x;
872 pix[pixcount].y = y;
873
874 /* set pixel value */
875 pix[pixcount].value = value;
876
877 /* NODATA */
878 if (hasband) {
879 if (exclude_nodata_value)
880 pix[pixcount].nodata = isnodata;
881 else
882 pix[pixcount].nodata = FALSE;
883 }
884 else {
885 pix[pixcount].nodata = isnodata;
886 }
887
888 /* update pixcount*/
889 pixcount++;
890 }
891 }
892
893 /* cleanup */
894 if (hasband) rt_band_destroy(band);
895 rt_raster_destroy(raster);
896 PG_FREE_IF_COPY(pgraster, 0);
897
898 /* shortcut if no pixcount */
899 if (pixcount < 1) {
900 elog(NOTICE, "No pixels found for band %d", nband);
901 MemoryContextSwitchTo(oldcontext);
902 SRF_RETURN_DONE(funcctx);
903 }
904
905 /* store pixel data */
906 funcctx->user_fctx = pix;
907
908 /* set total number of tuples to be returned */
909 funcctx->max_calls = pixcount;
910 POSTGIS_RT_DEBUGF(3, "pixcount = %d", pixcount);
911
912 /* build a tuple descriptor for our result type */
913 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
914 ereport(ERROR, (
915 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
916 errmsg("function returning record called in context that cannot accept type record")
917 ));
918 }
919 /* set tuple descriptor */
920 BlessTupleDesc(tupdesc);
921 funcctx->tuple_desc = tupdesc;
922
923 MemoryContextSwitchTo(oldcontext);
924 }
925
926 /* stuff done on every call of the function */
927 funcctx = SRF_PERCALL_SETUP();
928
929 call_cntr = funcctx->call_cntr;
930 max_calls = funcctx->max_calls;
931 tupdesc = funcctx->tuple_desc;
932 pix2 = funcctx->user_fctx;
933
934 /* do when there is more left to send */
935 if (call_cntr < max_calls) {
936 Datum values[VALUES_LENGTH];
937 bool nulls[VALUES_LENGTH];
938 HeapTuple tuple;
939 Datum result;
940
941 GSERIALIZED *gser = NULL;
942 size_t gser_size = 0;
943
944 POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr);
945
946 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
947
948 /* convert LWGEOM go GSERIALIZED */
949 gser = gserialized_from_lwgeom(pix2[call_cntr].geom, &gser_size);
950 lwgeom_free(pix2[call_cntr].geom); /* free geometry object */
951
952 /* set result tuple data */
953 values[0] = PointerGetDatum(gser);
954 if (pix2[call_cntr].nodata)
955 nulls[1] = TRUE;
956 else
957 values[1] = Float8GetDatum(pix2[call_cntr].value);
958 values[2] = Int32GetDatum(pix2[call_cntr].x);
959 values[3] = Int32GetDatum(pix2[call_cntr].y);
960
961 /* build a tuple */
962 tuple = heap_form_tuple(tupdesc, values, nulls);
963
964 /* make the tuple into a datum */
965 result = HeapTupleGetDatum(tuple);
966
967 SRF_RETURN_NEXT(funcctx, result);
968 }
969 /* do when there is no more left */
970 else {
971 pfree(pix2);
972 SRF_RETURN_DONE(funcctx);
973 }
974}
975
979PG_FUNCTION_INFO_V1(RASTER_getPolygon);
980Datum RASTER_getPolygon(PG_FUNCTION_ARGS)
981{
982 rt_pgraster *pgraster = NULL;
983 rt_raster raster = NULL;
984 int num_bands = 0;
985 int nband = 1;
986 int err;
987 LWMPOLY *surface = NULL;
988 GSERIALIZED *rtn = NULL;
989
990 /* raster */
991 if (PG_ARGISNULL(0))
992 PG_RETURN_NULL();
993 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
994
995 raster = rt_raster_deserialize(pgraster, FALSE);
996 if (!raster) {
997 PG_FREE_IF_COPY(pgraster, 0);
998 elog(ERROR, "RASTER_getPolygon: Could not deserialize raster");
999 PG_RETURN_NULL();
1000 }
1001
1002 /* num_bands */
1003 num_bands = rt_raster_get_num_bands(raster);
1004 if (num_bands < 1) {
1005 elog(NOTICE, "Raster provided has no bands");
1006 rt_raster_destroy(raster);
1007 PG_FREE_IF_COPY(pgraster, 0);
1008 PG_RETURN_NULL();
1009 }
1010
1011 /* band index is 1-based */
1012 if (!PG_ARGISNULL(1))
1013 nband = PG_GETARG_INT32(1);
1014 if (nband < 1 || nband > num_bands) {
1015 elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1016 rt_raster_destroy(raster);
1017 PG_FREE_IF_COPY(pgraster, 0);
1018 PG_RETURN_NULL();
1019 }
1020
1021 /* get band surface */
1022 err = rt_raster_surface(raster, nband - 1, &surface);
1023 rt_raster_destroy(raster);
1024 PG_FREE_IF_COPY(pgraster, 0);
1025
1026 if (err != ES_NONE) {
1027 elog(ERROR, "RASTER_getPolygon: Could not get raster band's surface");
1028 PG_RETURN_NULL();
1029 }
1030 else if (surface == NULL) {
1031 elog(NOTICE, "Raster is empty or all pixels of band are NODATA. Returning NULL");
1032 PG_RETURN_NULL();
1033 }
1034
1035 rtn = geometry_serialize(lwmpoly_as_lwgeom(surface));
1036 lwmpoly_free(surface);
1037
1038 PG_RETURN_POINTER(rtn);
1039}
1040
1044PG_FUNCTION_INFO_V1(RASTER_asRaster);
1045Datum RASTER_asRaster(PG_FUNCTION_ARGS)
1046{
1047 GSERIALIZED *gser = NULL;
1048
1049 LWGEOM *geom = NULL;
1050 rt_raster rast = NULL;
1051 rt_pgraster *pgrast = NULL;
1052
1053 lwvarlena_t *wkb;
1054 unsigned char variant = WKB_SFSQL;
1055
1056 double scale[2] = {0};
1057 double *scale_x = NULL;
1058 double *scale_y = NULL;
1059
1060 int dim[2] = {0};
1061 int *dim_x = NULL;
1062 int *dim_y = NULL;
1063
1064 ArrayType *array;
1065 Oid etype;
1066 Datum *e;
1067 bool *nulls;
1068 int16 typlen;
1069 bool typbyval;
1070 char typalign;
1071 int n = 0;
1072 int i = 0;
1073 int j = 0;
1074 int haserr = 0;
1075
1076 text *pixeltypetext = NULL;
1077 char *pixeltype = NULL;
1078 rt_pixtype pixtype = PT_END;
1079 rt_pixtype *pixtypes = NULL;
1080 uint32_t pixtypes_len = 0;
1081
1082 double *values = NULL;
1083 uint32_t values_len = 0;
1084
1085 uint8_t *hasnodatas = NULL;
1086 double *nodatavals = NULL;
1087 uint32_t nodatavals_len = 0;
1088
1089 double ulw[2] = {0};
1090 double *ul_xw = NULL;
1091 double *ul_yw = NULL;
1092
1093 double gridw[2] = {0};
1094 double *grid_xw = NULL;
1095 double *grid_yw = NULL;
1096
1097 double skew[2] = {0};
1098 double *skew_x = NULL;
1099 double *skew_y = NULL;
1100
1101 char **options = NULL;
1102 int options_len = 0;
1103
1104 uint32_t num_bands = 0;
1105
1106 int32_t srid = SRID_UNKNOWN;
1107 char *srs = NULL;
1108
1109 POSTGIS_RT_DEBUG(3, "RASTER_asRaster: Starting");
1110
1111 /* based upon LWGEOM_asBinary function in postgis/lwgeom_ogc.c */
1112
1113 /* Get the geometry */
1114 if (PG_ARGISNULL(0))
1115 PG_RETURN_NULL();
1116
1117 gser = PG_GETARG_GSERIALIZED_P(0);
1118 geom = lwgeom_from_gserialized(gser);
1119
1120 /* Get a 2D version of the geometry if necessary */
1121 if (lwgeom_ndims(geom) > 2) {
1122 LWGEOM *geom2d = lwgeom_force_2d(geom);
1123 lwgeom_free(geom);
1124 geom = geom2d;
1125 }
1126
1127 /* empty geometry, return empty raster */
1128 if (lwgeom_is_empty(geom)) {
1129 POSTGIS_RT_DEBUG(3, "Input geometry is empty. Returning empty raster");
1130 lwgeom_free(geom);
1131 PG_FREE_IF_COPY(gser, 0);
1132
1133 rast = rt_raster_new(0, 0);
1134 if (rast == NULL)
1135 PG_RETURN_NULL();
1136
1137 pgrast = rt_raster_serialize(rast);
1138 rt_raster_destroy(rast);
1139
1140 if (NULL == pgrast)
1141 PG_RETURN_NULL();
1142
1143 SET_VARSIZE(pgrast, pgrast->size);
1144 PG_RETURN_POINTER(pgrast);
1145 }
1146
1147 /* scale x */
1148 if (!PG_ARGISNULL(1)) {
1149 scale[0] = PG_GETARG_FLOAT8(1);
1150 if (FLT_NEQ(scale[0], 0.0))
1151 scale_x = &scale[0];
1152 }
1153
1154 /* scale y */
1155 if (!PG_ARGISNULL(2)) {
1156 scale[1] = PG_GETARG_FLOAT8(2);
1157 if (FLT_NEQ(scale[1], 0.0))
1158 scale_y = &scale[1];
1159 }
1160 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: scale (x, y) = %f, %f", scale[0], scale[1]);
1161
1162 /* width */
1163 if (!PG_ARGISNULL(3)) {
1164 dim[0] = PG_GETARG_INT32(3);
1165 if (dim[0] < 0) dim[0] = 0;
1166 if (dim[0] != 0) dim_x = &dim[0];
1167 }
1168
1169 /* height */
1170 if (!PG_ARGISNULL(4)) {
1171 dim[1] = PG_GETARG_INT32(4);
1172 if (dim[1] < 0) dim[1] = 0;
1173 if (dim[1] != 0) dim_y = &dim[1];
1174 }
1175 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: dim (x, y) = %d, %d", dim[0], dim[1]);
1176
1177 /* pixeltype */
1178 if (!PG_ARGISNULL(5)) {
1179 array = PG_GETARG_ARRAYTYPE_P(5);
1180 etype = ARR_ELEMTYPE(array);
1181 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1182
1183 switch (etype) {
1184 case TEXTOID:
1185 break;
1186 default:
1187
1188 lwgeom_free(geom);
1189 PG_FREE_IF_COPY(gser, 0);
1190
1191 elog(ERROR, "RASTER_asRaster: Invalid data type for pixeltype");
1192 PG_RETURN_NULL();
1193 break;
1194 }
1195
1196 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1197 &nulls, &n);
1198
1199 if (n) {
1200 pixtypes = (rt_pixtype *) palloc(sizeof(rt_pixtype) * n);
1201 /* clean each pixeltype */
1202 for (i = 0, j = 0; i < n; i++) {
1203 if (nulls[i]) {
1204 pixtypes[j++] = PT_64BF;
1205 continue;
1206 }
1207
1208 pixeltype = NULL;
1209 switch (etype) {
1210 case TEXTOID:
1211 pixeltypetext = (text *) DatumGetPointer(e[i]);
1212 if (NULL == pixeltypetext) break;
1213 pixeltype = text_to_cstring(pixeltypetext);
1214
1215 /* trim string */
1216 pixeltype = rtpg_trim(pixeltype);
1217 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: pixeltype is '%s'", pixeltype);
1218 break;
1219 }
1220
1221 if (strlen(pixeltype)) {
1222 pixtype = rt_pixtype_index_from_name(pixeltype);
1223 if (pixtype == PT_END) {
1224
1225 pfree(pixtypes);
1226
1227 lwgeom_free(geom);
1228 PG_FREE_IF_COPY(gser, 0);
1229
1230 elog(ERROR, "RASTER_asRaster: Invalid pixel type provided: %s", pixeltype);
1231 PG_RETURN_NULL();
1232 }
1233
1234 pixtypes[j] = pixtype;
1235 j++;
1236 }
1237 }
1238
1239 if (j > 0) {
1240 /* trim allocation */
1241 pixtypes = repalloc(pixtypes, j * sizeof(rt_pixtype));
1242 pixtypes_len = j;
1243 }
1244 else {
1245 pfree(pixtypes);
1246 pixtypes = NULL;
1247 pixtypes_len = 0;
1248 }
1249 }
1250 }
1251#if POSTGIS_DEBUG_LEVEL > 0
1252 for (uint32_t u = 0; u < pixtypes_len; u++)
1253 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: pixtypes[%u] = %u", i, pixtypes[i]);
1254#endif
1255
1256 /* value */
1257 if (!PG_ARGISNULL(6)) {
1258 array = PG_GETARG_ARRAYTYPE_P(6);
1259 etype = ARR_ELEMTYPE(array);
1260 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1261
1262 switch (etype) {
1263 case FLOAT4OID:
1264 case FLOAT8OID:
1265 break;
1266 default:
1267
1268 if (pixtypes_len) pfree(pixtypes);
1269
1270 lwgeom_free(geom);
1271 PG_FREE_IF_COPY(gser, 0);
1272
1273 elog(ERROR, "RASTER_asRaster: Invalid data type for value");
1274 PG_RETURN_NULL();
1275 break;
1276 }
1277
1278 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1279 &nulls, &n);
1280
1281 if (n) {
1282 values = (double *) palloc(sizeof(double) * n);
1283 for (i = 0, j = 0; i < n; i++) {
1284 if (nulls[i]) {
1285 values[j++] = 1;
1286 continue;
1287 }
1288
1289 switch (etype) {
1290 case FLOAT4OID:
1291 values[j] = (double) DatumGetFloat4(e[i]);
1292 break;
1293 case FLOAT8OID:
1294 values[j] = (double) DatumGetFloat8(e[i]);
1295 break;
1296 }
1297 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: values[%d] = %f", j, values[j]);
1298
1299 j++;
1300 }
1301
1302 if (j > 0) {
1303 /* trim allocation */
1304 values = repalloc(values, j * sizeof(double));
1305 values_len = j;
1306 }
1307 else {
1308 pfree(values);
1309 values = NULL;
1310 values_len = 0;
1311 }
1312 }
1313 }
1314#if POSTGIS_DEBUG_LEVEL > 0
1315 for (uint32_t u = 0; u < values_len; u++)
1316 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: values[%u] = %f", i, values[i]);
1317#endif
1318
1319 /* nodataval */
1320 if (!PG_ARGISNULL(7)) {
1321 array = PG_GETARG_ARRAYTYPE_P(7);
1322 etype = ARR_ELEMTYPE(array);
1323 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1324
1325 switch (etype) {
1326 case FLOAT4OID:
1327 case FLOAT8OID:
1328 break;
1329 default:
1330
1331 if (pixtypes_len) pfree(pixtypes);
1332 if (values_len) pfree(values);
1333
1334 lwgeom_free(geom);
1335 PG_FREE_IF_COPY(gser, 0);
1336
1337 elog(ERROR, "RASTER_asRaster: Invalid data type for nodataval");
1338 PG_RETURN_NULL();
1339 break;
1340 }
1341
1342 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1343 &nulls, &n);
1344
1345 if (n) {
1346 nodatavals = (double *) palloc(sizeof(double) * n);
1347 hasnodatas = (uint8_t *) palloc(sizeof(uint8_t) * n);
1348 for (i = 0, j = 0; i < n; i++) {
1349 if (nulls[i]) {
1350 hasnodatas[j] = 0;
1351 nodatavals[j] = 0;
1352 j++;
1353 continue;
1354 }
1355
1356 hasnodatas[j] = 1;
1357 switch (etype) {
1358 case FLOAT4OID:
1359 nodatavals[j] = (double) DatumGetFloat4(e[i]);
1360 break;
1361 case FLOAT8OID:
1362 nodatavals[j] = (double) DatumGetFloat8(e[i]);
1363 break;
1364 }
1365 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: hasnodatas[%d] = %d", j, hasnodatas[j]);
1366 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: nodatavals[%d] = %f", j, nodatavals[j]);
1367
1368 j++;
1369 }
1370
1371 if (j > 0) {
1372 /* trim allocation */
1373 nodatavals = repalloc(nodatavals, j * sizeof(double));
1374 hasnodatas = repalloc(hasnodatas, j * sizeof(uint8_t));
1375 nodatavals_len = j;
1376 }
1377 else {
1378 pfree(nodatavals);
1379 pfree(hasnodatas);
1380 nodatavals = NULL;
1381 hasnodatas = NULL;
1382 nodatavals_len = 0;
1383 }
1384 }
1385 }
1386#if POSTGIS_DEBUG_LEVEL > 0
1387 for (uint32_t u = 0; u < nodatavals_len; u++)
1388 {
1389 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: hasnodatas[%u] = %d", u, hasnodatas[u]);
1390 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: nodatavals[%u] = %f", u, nodatavals[u]);
1391 }
1392#endif
1393
1394 /* upperleftx */
1395 if (!PG_ARGISNULL(8)) {
1396 ulw[0] = PG_GETARG_FLOAT8(8);
1397 ul_xw = &ulw[0];
1398 }
1399
1400 /* upperlefty */
1401 if (!PG_ARGISNULL(9)) {
1402 ulw[1] = PG_GETARG_FLOAT8(9);
1403 ul_yw = &ulw[1];
1404 }
1405 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: upperleft (x, y) = %f, %f", ulw[0], ulw[1]);
1406
1407 /* gridx */
1408 if (!PG_ARGISNULL(10)) {
1409 gridw[0] = PG_GETARG_FLOAT8(10);
1410 grid_xw = &gridw[0];
1411 }
1412
1413 /* gridy */
1414 if (!PG_ARGISNULL(11)) {
1415 gridw[1] = PG_GETARG_FLOAT8(11);
1416 grid_yw = &gridw[1];
1417 }
1418 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: grid (x, y) = %f, %f", gridw[0], gridw[1]);
1419
1420 /* check dependent variables */
1421 haserr = 0;
1422 do {
1423 /* only part of scale provided */
1424 if (
1425 (scale_x == NULL && scale_y != NULL) ||
1426 (scale_x != NULL && scale_y == NULL)
1427 ) {
1428 elog(NOTICE, "Values must be provided for both X and Y of scale if one is specified");
1429 haserr = 1;
1430 break;
1431 }
1432
1433 /* only part of dimension provided */
1434 if (
1435 (dim_x == NULL && dim_y != NULL) ||
1436 (dim_x != NULL && dim_y == NULL)
1437 ) {
1438 elog(NOTICE, "Values must be provided for both width and height if one is specified");
1439 haserr = 1;
1440 break;
1441 }
1442
1443 /* scale and dimension provided */
1444 if (
1445 (scale_x != NULL && scale_y != NULL) &&
1446 (dim_x != NULL && dim_y != NULL)
1447 ) {
1448 elog(NOTICE, "Values provided for X and Y of scale and width and height. Using the width and height");
1449 scale_x = NULL;
1450 scale_y = NULL;
1451 break;
1452 }
1453
1454 /* neither scale or dimension provided */
1455 if (
1456 (scale_x == NULL && scale_y == NULL) &&
1457 (dim_x == NULL && dim_y == NULL)
1458 ) {
1459 elog(NOTICE, "Values must be provided for X and Y of scale or width and height");
1460 haserr = 1;
1461 break;
1462 }
1463
1464 /* only part of upper-left provided */
1465 if (
1466 (ul_xw == NULL && ul_yw != NULL) ||
1467 (ul_xw != NULL && ul_yw == NULL)
1468 ) {
1469 elog(NOTICE, "Values must be provided for both X and Y when specifying the upper-left corner");
1470 haserr = 1;
1471 break;
1472 }
1473
1474 /* only part of alignment provided */
1475 if (
1476 (grid_xw == NULL && grid_yw != NULL) ||
1477 (grid_xw != NULL && grid_yw == NULL)
1478 ) {
1479 elog(NOTICE, "Values must be provided for both X and Y when specifying the alignment");
1480 haserr = 1;
1481 break;
1482 }
1483
1484 /* upper-left and alignment provided */
1485 if (
1486 (ul_xw != NULL && ul_yw != NULL) &&
1487 (grid_xw != NULL && grid_yw != NULL)
1488 ) {
1489 elog(NOTICE, "Values provided for both X and Y of upper-left corner and alignment. Using the values of upper-left corner");
1490 grid_xw = NULL;
1491 grid_yw = NULL;
1492 break;
1493 }
1494 }
1495 while (0);
1496
1497 if (haserr) {
1498 if (pixtypes_len) pfree(pixtypes);
1499 if (values_len) pfree(values);
1500 if (nodatavals_len) {
1501 pfree(nodatavals);
1502 pfree(hasnodatas);
1503 }
1504
1505 lwgeom_free(geom);
1506 PG_FREE_IF_COPY(gser, 0);
1507
1508 PG_RETURN_NULL();
1509 }
1510
1511 /* skewx */
1512 if (!PG_ARGISNULL(12)) {
1513 skew[0] = PG_GETARG_FLOAT8(12);
1514 if (FLT_NEQ(skew[0], 0.0))
1515 skew_x = &skew[0];
1516 }
1517
1518 /* skewy */
1519 if (!PG_ARGISNULL(13)) {
1520 skew[1] = PG_GETARG_FLOAT8(13);
1521 if (FLT_NEQ(skew[1], 0.0))
1522 skew_y = &skew[1];
1523 }
1524 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: skew (x, y) = %f, %f", skew[0], skew[1]);
1525
1526 /* all touched */
1527 if (!PG_ARGISNULL(14) && PG_GETARG_BOOL(14) == TRUE) {
1528 if (options_len == 0) {
1529 options_len = 1;
1530 options = (char **) palloc(sizeof(char *) * options_len);
1531 }
1532 else {
1533 options_len++;
1534 options = (char **) repalloc(options, sizeof(char *) * options_len);
1535 }
1536
1537 options[options_len - 1] = palloc(sizeof(char*) * (strlen("ALL_TOUCHED=TRUE") + 1));
1538 strcpy(options[options_len - 1], "ALL_TOUCHED=TRUE");
1539 }
1540
1541 if (options_len) {
1542 options_len++;
1543 options = (char **) repalloc(options, sizeof(char *) * options_len);
1544 options[options_len - 1] = NULL;
1545 }
1546
1547 /* get geometry's srid */
1548 srid = gserialized_get_srid(gser);
1549
1550 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: srid = %d", srid);
1551 if (clamp_srid(srid) != SRID_UNKNOWN) {
1552 srs = rtpg_getSR(srid);
1553 if (NULL == srs) {
1554
1555 if (pixtypes_len) pfree(pixtypes);
1556 if (values_len) pfree(values);
1557 if (nodatavals_len) {
1558 pfree(hasnodatas);
1559 pfree(nodatavals);
1560 }
1561 if (options_len) pfree(options);
1562
1563 lwgeom_free(geom);
1564 PG_FREE_IF_COPY(gser, 0);
1565
1566 elog(ERROR, "RASTER_asRaster: Could not find srtext for SRID (%d)", srid);
1567 PG_RETURN_NULL();
1568 }
1569 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: srs is %s", srs);
1570 }
1571 else
1572 srs = NULL;
1573
1574 /* determine number of bands */
1575 /* MIN macro is from GDAL's cpl_port.h */
1576 num_bands = MIN(pixtypes_len, values_len);
1577 num_bands = MIN(num_bands, nodatavals_len);
1578 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: pixtypes_len = %d", pixtypes_len);
1579 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: values_len = %d", values_len);
1580 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: nodatavals_len = %d", nodatavals_len);
1581 POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: num_bands = %d", num_bands);
1582
1583 /* warn of imbalanced number of band elements */
1584 if (!(
1585 (pixtypes_len == values_len) &&
1586 (values_len == nodatavals_len)
1587 )) {
1588 elog(
1589 NOTICE,
1590 "Imbalanced number of values provided for pixeltype (%d), value (%d) and nodataval (%d). Using the first %d values of each parameter",
1591 pixtypes_len,
1592 values_len,
1593 nodatavals_len,
1594 num_bands
1595 );
1596 }
1597
1598 /* get wkb of geometry */
1599 POSTGIS_RT_DEBUG(3, "RASTER_asRaster: getting wkb of geometry");
1600 wkb = lwgeom_to_wkb_varlena(geom, variant);
1601 lwgeom_free(geom);
1602 PG_FREE_IF_COPY(gser, 0);
1603
1604 /* rasterize geometry */
1605 POSTGIS_RT_DEBUG(3, "RASTER_asRaster: rasterizing geometry");
1606 /* use nodatavals for the init parameter */
1607 rast = rt_raster_gdal_rasterize((unsigned char *)wkb->data,
1608 LWSIZE_GET(wkb->size) - LWVARHDRSZ,
1609 srs,
1610 num_bands,
1611 pixtypes,
1612 nodatavals,
1613 values,
1614 nodatavals,
1615 hasnodatas,
1616 dim_x,
1617 dim_y,
1618 scale_x,
1619 scale_y,
1620 ul_xw,
1621 ul_yw,
1622 grid_xw,
1623 grid_yw,
1624 skew_x,
1625 skew_y,
1626 options);
1627
1628 if (pixtypes_len) pfree(pixtypes);
1629 if (values_len) pfree(values);
1630 if (nodatavals_len) {
1631 pfree(hasnodatas);
1632 pfree(nodatavals);
1633 }
1634 if (options_len) pfree(options);
1635
1636 if (!rast) {
1637 elog(ERROR, "RASTER_asRaster: Could not rasterize geometry");
1638 PG_RETURN_NULL();
1639 }
1640
1641 /* add target srid */
1642 rt_raster_set_srid(rast, srid);
1643
1644 pgrast = rt_raster_serialize(rast);
1645 rt_raster_destroy(rast);
1646
1647 if (NULL == pgrast) PG_RETURN_NULL();
1648
1649 POSTGIS_RT_DEBUG(3, "RASTER_asRaster: done");
1650
1651 SET_VARSIZE(pgrast, pgrast->size);
1652 PG_RETURN_POINTER(pgrast);
1653}
variant
static uint8_t variant
Definition cu_in_twkb.c:26
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
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
lwpoint_free
void lwpoint_free(LWPOINT *pt)
Definition lwpoint.c:213
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition lwgeom.c:1246
lwmpoly_free
void lwmpoly_free(LWMPOLY *mpoly)
Definition lwmpoly.c:53
LWSIZE_GET
#define LWSIZE_GET(varsize)
Macro for reading the size from the GSERIALIZED size attribute.
Definition liblwgeom.h:324
lwgeom_to_wkb_varlena
lwvarlena_t * lwgeom_to_wkb_varlena(const LWGEOM *geom, uint8_t variant)
Definition lwout_wkb.c:851
lwgeom_force_2d
LWGEOM * lwgeom_force_2d(const LWGEOM *geom)
Strip out the Z/M components of an LWGEOM.
Definition lwgeom.c:821
lwpoly_free
void lwpoly_free(LWPOLY *poly)
Definition lwpoly.c:175
SRID_UNKNOWN
#define SRID_UNKNOWN
Unknown SRID value.
Definition liblwgeom.h:215
lwmpoly_as_lwgeom
LWGEOM * lwmpoly_as_lwgeom(const LWMPOLY *obj)
Definition lwgeom.c:322
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
FLT_NEQ
#define FLT_NEQ(x, y)
Definition librtcore.h:2435
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_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:833
rt_pixtype_index_from_name
rt_pixtype rt_pixtype_index_from_name(const char *pixname)
Definition rt_pixel.c:82
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
PT_64BF
@ PT_64BF
Definition librtcore.h:200
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_surface
rt_errorstate rt_raster_surface(rt_raster raster, int nband, LWMPOLY **surface)
Get a raster as a surface (multipolygon).
Definition rt_geometry.c:356
rt_raster_new
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition rt_raster.c:52
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_raster_gdal_polygonize
rt_geomval rt_raster_gdal_polygonize(rt_raster raster, int nband, int exclude_nodata_value, int *pnElements)
Returns a set of "geomval" value, one for each group of pixel sharing the same value for the provided...
Definition rt_geometry.c:976
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_convex_hull
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition rt_geometry.c:839
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition rt_raster.c:133
rt_raster_set_srid
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition rt_raster.c:367
rt_raster_pixel_as_centroid_point
LWPOINT * rt_raster_pixel_as_centroid_point(rt_raster rast, int x, int y)
Get a raster pixel centroid point.
Definition rt_geometry.c:669
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition rt_raster.c:125
rt_raster_get_envelope_geom
rt_errorstate rt_raster_get_envelope_geom(rt_raster raster, LWGEOM **env)
Get raster's envelope as a geometry.
Definition rt_geometry.c:706
rt_raster_serialize
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition rt_serialize.c:522
rt_raster_get_perimeter
rt_errorstate rt_raster_get_perimeter(rt_raster raster, int nband, LWGEOM **perimeter)
Get raster perimeter.
Definition rt_geometry.c:183
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_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
RASTER_dumpAsPolygons
Datum RASTER_dumpAsPolygons(PG_FUNCTION_ARGS)
Definition rtpg_geometry.c:194
RASTER_getPolygon
Datum RASTER_getPolygon(PG_FUNCTION_ARGS)
Definition rtpg_geometry.c:980
RASTER_convex_hull
Datum RASTER_convex_hull(PG_FUNCTION_ARGS)
Definition rtpg_geometry.c:119
RASTER_envelope
Datum RASTER_envelope(PG_FUNCTION_ARGS)
Definition rtpg_geometry.c:66
RASTER_getPixelCentroids
Datum RASTER_getPixelCentroids(PG_FUNCTION_ARGS)
Definition rtpg_geometry.c:664
RASTER_asRaster
Datum RASTER_asRaster(PG_FUNCTION_ARGS)
Definition rtpg_geometry.c:1045
PG_FUNCTION_INFO_V1
PG_FUNCTION_INFO_V1(RASTER_envelope)
VALUES_LENGTH
#define VALUES_LENGTH
Definition rtpg_geometry.c:191
RASTER_getPixelPolygons
Datum RASTER_getPixelPolygons(PG_FUNCTION_ARGS)
Definition rtpg_geometry.c:355
rtpg_getSR
char * rtpg_getSR(int32_t srid)
Definition rtpg_internal.c:282
rtpg_trim
char * rtpg_trim(const char *input)
Definition rtpg_internal.c:226
rtpg_internal.h
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
MIN
#define MIN(a, b)
Definition shpopen.c:64
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_geomval_t
Definition librtcore.h:2554
rt_pixel_t::y
int y
Definition librtcore.h:2537
rt_pixel_t::geom
LWGEOM * geom
Definition librtcore.h:2542
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::srid
int32_t srid
Definition librtcore.h:2477
rt_raster_serialized_t
Struct definitions.
Definition librtcore.h:2452
rt_raster_t
Definition librtcore.h:2486