PostGIS 3.6.2dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ rt_raster_surface()

rt_errorstate rt_raster_surface ( rt_raster  raster,
int  nband,
LWMPOLY **  surface 
)

Get a raster as a surface (multipolygon).

If a band is specified, those pixels with value (not NODATA) contribute to the area of the output multipolygon.

Parameters
raster: the raster to convert to a multipolygon
nband: the 0-based band of raster rast to use if value is less than zero, bands are ignored.
**surface: raster as a surface (multipolygon). if all pixels are NODATA, NULL is set
Returns
ES_NONE on success, ES_ERROR on error

If a band is specified, those pixels with value (not NODATA) contribute to the area of the output multipolygon.

Parameters
raster: the raster to convert to a multipolygon
nband: the 0-based band of raster rast to use if value is less than zero, bands are ignored.
*surface: raster as a surface (multipolygon). if all pixels are NODATA, NULL is set
Returns
ES_NONE on success, ES_ERROR on error

Definition at line 355 of file rt_geometry.c.

355 {
356 rt_band band = NULL;
357 LWGEOM *mpoly = NULL;
358 LWGEOM *tmp = NULL;
359 LWGEOM *clone = NULL;
360 rt_geomval gv = NULL;
361 int gvcount = 0;
362 GEOSGeometry *gc = NULL;
363 GEOSGeometry *gunion = NULL;
364 GEOSGeometry **geoms = NULL;
365 int geomscount = 0;
366 int i = 0;
367
368 assert(surface != NULL);
369
370 /* init *surface to NULL */
371 *surface = NULL;
372
373 /* raster is empty, surface = NULL */
374 if (rt_raster_is_empty(raster))
375 return ES_NONE;
376
377 /* if nband < 0, return the convex hull as a multipolygon */
378 if (nband < 0) {
379 /*
380 lwgeom_as_multi() only does a shallow clone internally
381 so input and output geometries may share memory
382 hence the deep clone of the output geometry for returning
383 is the only way to guarantee the memory isn't shared
384 */
385 if (rt_raster_get_convex_hull(raster, &tmp) != ES_NONE) {
386 rterror("rt_raster_surface: Could not get convex hull of raster");
387 return ES_ERROR;
388 }
389 mpoly = lwgeom_as_multi(tmp);
390 clone = lwgeom_clone_deep(mpoly);
391 lwgeom_free(tmp);
392 lwgeom_free(mpoly);
393
394 *surface = lwgeom_as_lwmpoly(clone);
395 return ES_NONE;
396 }
397 /* check that nband is valid */
398 else if (nband >= rt_raster_get_num_bands(raster)) {
399 rterror("rt_raster_surface: The band index %d is invalid", nband);
400 return ES_ERROR;
401 }
402
403 /* get band */
404 band = rt_raster_get_band(raster, nband);
405 if (band == NULL) {
406 rterror("rt_raster_surface: Error getting band %d from raster", nband);
407 return ES_ERROR;
408 }
409
410 /* band does not have a NODATA flag, return convex hull */
411 if (!rt_band_get_hasnodata_flag(band)) {
412 /*
413 lwgeom_as_multi() only does a shallow clone internally
414 so input and output geometries may share memory
415 hence the deep clone of the output geometry for returning
416 is the only way to guarantee the memory isn't shared
417 */
418 if (rt_raster_get_convex_hull(raster, &tmp) != ES_NONE) {
419 rterror("rt_raster_surface: Could not get convex hull of raster");
420 return ES_ERROR;
421 }
422 mpoly = lwgeom_as_multi(tmp);
423 clone = lwgeom_clone_deep(mpoly);
424 lwgeom_free(tmp);
425 lwgeom_free(mpoly);
426
427 *surface = lwgeom_as_lwmpoly(clone);
428 return ES_NONE;
429 }
430 /* band is NODATA, return NULL */
431 else if (rt_band_get_isnodata_flag(band)) {
432 RASTER_DEBUG(3, "Band is NODATA. Returning NULL");
433 return ES_NONE;
434 }
435
436 /* use gdal polygonize */
437 gv = rt_raster_gdal_polygonize(raster, nband, 1, &gvcount);
438 /* no polygons returned */
439 if (gvcount < 1) {
440 RASTER_DEBUG(3, "All pixels of band are NODATA. Returning NULL");
441 if (gv != NULL) rtdealloc(gv);
442 return ES_NONE;
443 }
444 /* more than 1 polygon */
445 else if (gvcount > 1) {
446 /* convert LWPOLY to GEOSGeometry */
447 geomscount = gvcount;
448 geoms = rtalloc(sizeof(GEOSGeometry *) * geomscount);
449 if (geoms == NULL) {
450 rterror("rt_raster_surface: Could not allocate memory for pixel polygons as GEOSGeometry");
451 for (i = 0; i < gvcount; i++) lwpoly_free(gv[i].geom);
452 rtdealloc(gv);
453 return ES_ERROR;
454 }
455 for (i = 0; i < gvcount; i++) {
456#if POSTGIS_DEBUG_LEVEL > 3
457 {
458 char *wkt = lwgeom_to_wkt(lwpoly_as_lwgeom(gv[i].geom), WKT_ISO, DBL_DIG, NULL);
459 RASTER_DEBUGF(4, "geom %d = %s", i, wkt);
460 rtdealloc(wkt);
461 }
462#endif
463
464 geoms[i] = LWGEOM2GEOS(lwpoly_as_lwgeom(gv[i].geom), 0);
465 lwpoly_free(gv[i].geom);
466 }
467 rtdealloc(gv);
468
469 /* create geometry collection */
470 gc = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, geomscount);
471
472 if (gc == NULL) {
473 rterror("rt_raster_surface: Could not create GEOS GEOMETRYCOLLECTION from set of pixel polygons");
474
475 for (i = 0; i < geomscount; i++)
476 GEOSGeom_destroy(geoms[i]);
477 rtdealloc(geoms);
478 return ES_ERROR;
479 }
480
481 /* run the union */
482 gunion = GEOSUnaryUnion(gc);
483
484 GEOSGeom_destroy(gc);
485 rtdealloc(geoms);
486
487 if (gunion == NULL) {
488 rterror("rt_raster_surface: Could not union the pixel polygons using GEOSUnaryUnion()");
489 return ES_ERROR;
490 }
491
492 /* convert union result from GEOSGeometry to LWGEOM */
493 mpoly = GEOS2LWGEOM(gunion, 0);
494
495 /*
496 is geometry valid?
497 if not, try to make valid
498 */
499 do {
500 LWGEOM *mpolyValid = NULL;
501
502 if (GEOSisValid(gunion))
503 break;
504
505 /* make geometry valid */
506 mpolyValid = lwgeom_make_valid(mpoly);
507 if (mpolyValid == NULL) {
508 rtwarn("Cannot fix invalid geometry");
509 break;
510 }
511
512 lwgeom_free(mpoly);
513 mpoly = mpolyValid;
514 }
515 while (0);
516
517 GEOSGeom_destroy(gunion);
518 }
519 else {
520 mpoly = lwpoly_as_lwgeom(gv[0].geom);
521 rtdealloc(gv);
522
523#if POSTGIS_DEBUG_LEVEL > 3
524 {
525 char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
526 RASTER_DEBUGF(4, "geom 0 = %s", wkt);
527 rtdealloc(wkt);
528 }
529#endif
530 }
531
532 /* specify SRID */
533 lwgeom_set_srid(mpoly, rt_raster_get_srid(raster));
534
535 if (mpoly != NULL) {
536 /* convert to multi */
537 if (!lwgeom_is_collection(mpoly)) {
538 tmp = mpoly;
539
540#if POSTGIS_DEBUG_LEVEL > 3
541 {
542 char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
543 RASTER_DEBUGF(4, "before multi = %s", wkt);
544 rtdealloc(wkt);
545 }
546#endif
547
548 RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
549
550 /*
551 lwgeom_as_multi() only does a shallow clone internally
552 so input and output geometries may share memory
553 hence the deep clone of the output geometry for returning
554 is the only way to guarantee the memory isn't shared
555 */
556 mpoly = lwgeom_as_multi(tmp);
557 clone = lwgeom_clone_deep(mpoly);
558 lwgeom_free(tmp);
559 lwgeom_free(mpoly);
560 mpoly = clone;
561
562 RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
563
564#if POSTGIS_DEBUG_LEVEL > 3
565 {
566 char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
567 RASTER_DEBUGF(4, "after multi = %s", wkt);
568 rtdealloc(wkt);
569 }
570#endif
571 }
572
573#if POSTGIS_DEBUG_LEVEL > 3
574 {
575 char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
576 RASTER_DEBUGF(4, "returning geometry = %s", wkt);
577 rtdealloc(wkt);
578 }
579#endif
580
581 *surface = lwgeom_as_lwmpoly(mpoly);
582 return ES_NONE;
583 }
584
585 return ES_NONE;
586}
LWGEOM2GEOS
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
Definition liblwgeom/lwgeom_geos.c:400
GEOS2LWGEOM
void(*) LWGEOM GEOS2LWGEOM)(const GEOSGeometry *geom, uint8_t want3d)
Definition liblwgeom/lwgeom_geos.h:39
lwgeom_set_srid
void lwgeom_set_srid(LWGEOM *geom, int32_t srid)
Set the SRID on an LWGEOM For collections, only the parent gets an SRID, all the children get SRID_UN...
Definition lwgeom.c:1610
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition lwgeom.c:1218
lwgeom_as_multi
LWGEOM * lwgeom_as_multi(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate MULTI* type.
Definition lwgeom.c:380
lwgeom_to_wkt
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition lwout_wkt.c:708
lwgeom_is_collection
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM contains sub-geometries or not This basically just checks that the struct ...
Definition lwgeom.c:1097
lwgeom_as_lwmpoly
LWMPOLY * lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
Definition lwgeom.c:260
WKT_ISO
#define WKT_ISO
Definition liblwgeom.h:2216
lwpoly_free
void lwpoly_free(LWPOLY *poly)
Definition lwpoly.c:175
lwpoly_as_lwgeom
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition lwgeom.c:329
lwgeom_make_valid
LWGEOM * lwgeom_make_valid(LWGEOM *geom)
Attempts to make an invalid geometries valid w/out losing points.
Definition liblwgeom/lwgeom_geos_clean.c:286
lwgeom_clone_deep
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition lwgeom.c:529
rterror
void rterror(const char *fmt,...) __attribute__((format(printf
Wrappers used for reporting errors and info.
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition rt_context.c:191
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition librtcore.h:302
rt_raster_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition rt_raster.c:360
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition librtcore.h:306
rtwarn
void void void rtwarn(const char *fmt,...) __attribute__((format(printf
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:825
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition rt_band.c:865
ES_NONE
@ ES_NONE
Definition librtcore.h:182
ES_ERROR
@ ES_ERROR
Definition librtcore.h:183
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:376
rtdealloc
void rtdealloc(void *mem)
Definition rt_context.c:206
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
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:975
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:838
rt_band_t
Definition librtcore.h:2503
rt_geomval_t
Definition librtcore.h:2542

References ES_ERROR, ES_NONE, LWGEOM2GEOS(), lwgeom_as_lwmpoly(), lwgeom_as_multi(), lwgeom_clone_deep(), lwgeom_free(), lwgeom_is_collection(), lwgeom_make_valid(), lwgeom_set_srid(), lwgeom_to_wkt(), lwpoly_as_lwgeom(), lwpoly_free(), RASTER_DEBUG, RASTER_DEBUGF, rt_band_get_hasnodata_flag(), rt_band_get_isnodata_flag(), rt_raster_gdal_polygonize(), rt_raster_get_band(), rt_raster_get_convex_hull(), rt_raster_get_num_bands(), rt_raster_get_srid(), rt_raster_is_empty(), rtalloc(), rtdealloc(), rterror(), rtwarn(), and WKT_ISO.

Referenced by RASTER_getPolygon(), rt_raster_fully_within_distance(), rt_raster_geos_spatial_relationship(), rt_raster_within_distance(), and test_raster_surface().

Here is the call graph for this function:
Here is the caller graph for this function: