PostGIS 3.7.0dev-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 356 of file rt_geometry.c.

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

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: