PostGIS 3.0.6dev-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

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 guarentee 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 guarentee 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 /* initialize GEOS */
437 initGEOS(rtinfo, lwgeom_geos_error);
438
439 /* use gdal polygonize */
440 gv = rt_raster_gdal_polygonize(raster, nband, 1, &gvcount);
441 /* no polygons returned */
442 if (gvcount < 1) {
443 RASTER_DEBUG(3, "All pixels of band are NODATA. Returning NULL");
444 if (gv != NULL) rtdealloc(gv);
445 return ES_NONE;
446 }
447 /* more than 1 polygon */
448 else if (gvcount > 1) {
449 /* convert LWPOLY to GEOSGeometry */
450 geomscount = gvcount;
451 geoms = rtalloc(sizeof(GEOSGeometry *) * geomscount);
452 if (geoms == NULL) {
453 rterror("rt_raster_surface: Could not allocate memory for pixel polygons as GEOSGeometry");
454 for (i = 0; i < gvcount; i++) lwpoly_free(gv[i].geom);
455 rtdealloc(gv);
456 return ES_ERROR;
457 }
458 for (i = 0; i < gvcount; i++) {
459#if POSTGIS_DEBUG_LEVEL > 3
460 {
461 char *wkt = lwgeom_to_wkt(lwpoly_as_lwgeom(gv[i].geom), WKT_ISO, DBL_DIG, NULL);
462 RASTER_DEBUGF(4, "geom %d = %s", i, wkt);
463 rtdealloc(wkt);
464 }
465#endif
466
467 geoms[i] = LWGEOM2GEOS(lwpoly_as_lwgeom(gv[i].geom), 0);
468 lwpoly_free(gv[i].geom);
469 }
470 rtdealloc(gv);
471
472 /* create geometry collection */
473 gc = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, geomscount);
474
475 if (gc == NULL) {
476 rterror("rt_raster_surface: Could not create GEOS GEOMETRYCOLLECTION from set of pixel polygons");
477
478 for (i = 0; i < geomscount; i++)
479 GEOSGeom_destroy(geoms[i]);
480 rtdealloc(geoms);
481 return ES_ERROR;
482 }
483
484 /* run the union */
485 gunion = GEOSUnaryUnion(gc);
486
487 GEOSGeom_destroy(gc);
488 rtdealloc(geoms);
489
490 if (gunion == NULL) {
491 rterror("rt_raster_surface: Could not union the pixel polygons using GEOSUnaryUnion()");
492 return ES_ERROR;
493 }
494
495 /* convert union result from GEOSGeometry to LWGEOM */
496 mpoly = GEOS2LWGEOM(gunion, 0);
497
498 /*
499 is geometry valid?
500 if not, try to make valid
501 */
502 do {
503 LWGEOM *mpolyValid = NULL;
504
505 if (GEOSisValid(gunion))
506 break;
507
508 /* make geometry valid */
509 mpolyValid = lwgeom_make_valid(mpoly);
510 if (mpolyValid == NULL) {
511 rtwarn("Cannot fix invalid geometry");
512 break;
513 }
514
515 lwgeom_free(mpoly);
516 mpoly = mpolyValid;
517 }
518 while (0);
519
520 GEOSGeom_destroy(gunion);
521 }
522 else {
523 mpoly = lwpoly_as_lwgeom(gv[0].geom);
524 rtdealloc(gv);
525
526#if POSTGIS_DEBUG_LEVEL > 3
527 {
528 char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
529 RASTER_DEBUGF(4, "geom 0 = %s", wkt);
530 rtdealloc(wkt);
531 }
532#endif
533 }
534
535 /* specify SRID */
536 lwgeom_set_srid(mpoly, rt_raster_get_srid(raster));
537
538 if (mpoly != NULL) {
539 /* convert to multi */
540 if (!lwgeom_is_collection(mpoly)) {
541 tmp = mpoly;
542
543#if POSTGIS_DEBUG_LEVEL > 3
544 {
545 char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
546 RASTER_DEBUGF(4, "before multi = %s", wkt);
547 rtdealloc(wkt);
548 }
549#endif
550
551 RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
552
553 /*
554 lwgeom_as_multi() only does a shallow clone internally
555 so input and output geometries may share memory
556 hence the deep clone of the output geometry for returning
557 is the only way to guarentee the memory isn't shared
558 */
559 mpoly = lwgeom_as_multi(tmp);
560 clone = lwgeom_clone_deep(mpoly);
561 lwgeom_free(tmp);
562 lwgeom_free(mpoly);
563 mpoly = clone;
564
565 RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
566
567#if POSTGIS_DEBUG_LEVEL > 3
568 {
569 char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
570 RASTER_DEBUGF(4, "after multi = %s", wkt);
571 rtdealloc(wkt);
572 }
573#endif
574 }
575
576#if POSTGIS_DEBUG_LEVEL > 3
577 {
578 char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
579 RASTER_DEBUGF(4, "returning geometry = %s", wkt);
580 rtdealloc(wkt);
581 }
582#endif
583
584 *surface = lwgeom_as_lwmpoly(mpoly);
585 return ES_NONE;
586 }
587
588 return ES_NONE;
589}
LWGEOM2GEOS
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
Definition liblwgeom/lwgeom_geos.c:390
GEOS2LWGEOM
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, uint8_t want3d)
Definition liblwgeom/lwgeom_geos.c:166
lwgeom_geos_error
void lwgeom_geos_error(const char *fmt,...)
Definition liblwgeom/lwgeom_geos.c:43
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:1530
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition lwgeom.c:1138
lwgeom_as_multi
LWGEOM * lwgeom_as_multi(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate MULTI* type.
Definition lwgeom.c:362
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:676
lwgeom_is_collection
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition lwgeom.c:1079
lwgeom_as_lwmpoly
LWMPOLY * lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
Definition lwgeom.c:242
WKT_ISO
#define WKT_ISO
Definition liblwgeom.h:2130
lwpoly_free
void lwpoly_free(LWPOLY *poly)
Definition lwpoly.c:175
lwpoly_as_lwgeom
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition lwgeom.c:311
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:854
lwgeom_clone_deep
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition lwgeom.c:511
rterror
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition rt_context.c:199
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition rt_context.c:171
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition librtcore.h:295
rt_raster_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition rt_raster.c:356
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition librtcore.h:299
rtinfo
void rtinfo(const char *fmt,...)
Definition rt_context.c:211
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:674
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition rt_band.c:714
rtwarn
void rtwarn(const char *fmt,...)
Definition rt_context.c:224
ES_NONE
@ ES_NONE
Definition librtcore.h:180
ES_ERROR
@ ES_ERROR
Definition librtcore.h:181
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:372
rtdealloc
void rtdealloc(void *mem)
Definition rt_context.c:186
rt_raster_is_empty
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition rt_raster.c:1329
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:381
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:940
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:803
rt_band_t
Definition librtcore.h:2314
rt_geomval_t
Definition librtcore.h:2353

References ES_ERROR, ES_NONE, GEOS2LWGEOM(), LWGEOM2GEOS(), lwgeom_as_lwmpoly(), lwgeom_as_multi(), lwgeom_clone_deep(), lwgeom_free(), lwgeom_geos_error(), 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(), rtinfo(), 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: