PostGIS  2.2.7dev-r@@SVN_REVISION@@
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.

References ovdump::band, 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().

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 #if POSTGIS_GEOS_VERSION >= 33
474  gc = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, geomscount);
475 #else
476  gc = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, geomscount);
477 #endif
478 
479  if (gc == NULL) {
480 #if POSTGIS_GEOS_VERSION >= 33
481  rterror("rt_raster_surface: Could not create GEOS GEOMETRYCOLLECTION from set of pixel polygons");
482 #else
483  rterror("rt_raster_surface: Could not create GEOS MULTIPOLYGON from set of pixel polygons");
484 #endif
485 
486  for (i = 0; i < geomscount; i++)
487  GEOSGeom_destroy(geoms[i]);
488  rtdealloc(geoms);
489  return ES_ERROR;
490  }
491 
492  /* run the union */
493 #if POSTGIS_GEOS_VERSION >= 33
494  gunion = GEOSUnaryUnion(gc);
495 #else
496  gunion = GEOSUnionCascaded(gc);
497 #endif
498  GEOSGeom_destroy(gc);
499  rtdealloc(geoms);
500 
501  if (gunion == NULL) {
502 #if POSTGIS_GEOS_VERSION >= 33
503  rterror("rt_raster_surface: Could not union the pixel polygons using GEOSUnaryUnion()");
504 #else
505  rterror("rt_raster_surface: Could not union the pixel polygons using GEOSUnionCascaded()");
506 #endif
507  return ES_ERROR;
508  }
509 
510  /* convert union result from GEOSGeometry to LWGEOM */
511  mpoly = GEOS2LWGEOM(gunion, 0);
512 
513  /*
514  is geometry valid?
515  if not, try to make valid
516  */
517  do {
518  LWGEOM *mpolyValid = NULL;
519 
520 #if POSTGIS_GEOS_VERSION < 33
521  break;
522 #endif
523 
524  if (GEOSisValid(gunion))
525  break;
526 
527  /* make geometry valid */
528  mpolyValid = lwgeom_make_valid(mpoly);
529  if (mpolyValid == NULL) {
530  rtwarn("Cannot fix invalid geometry");
531  break;
532  }
533 
534  lwgeom_free(mpoly);
535  mpoly = mpolyValid;
536  }
537  while (0);
538 
539  GEOSGeom_destroy(gunion);
540  }
541  else {
542  mpoly = lwpoly_as_lwgeom(gv[0].geom);
543  rtdealloc(gv);
544 
545 #if POSTGIS_DEBUG_LEVEL > 3
546  {
547  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
548  RASTER_DEBUGF(4, "geom 0 = %s", wkt);
549  rtdealloc(wkt);
550  }
551 #endif
552  }
553 
554  /* specify SRID */
555  lwgeom_set_srid(mpoly, rt_raster_get_srid(raster));
556 
557  if (mpoly != NULL) {
558  /* convert to multi */
559  if (!lwgeom_is_collection(mpoly)) {
560  tmp = mpoly;
561 
562 #if POSTGIS_DEBUG_LEVEL > 3
563  {
564  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
565  RASTER_DEBUGF(4, "before multi = %s", wkt);
566  rtdealloc(wkt);
567  }
568 #endif
569 
570  RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
571 
572  /*
573  lwgeom_as_multi() only does a shallow clone internally
574  so input and output geometries may share memory
575  hence the deep clone of the output geometry for returning
576  is the only way to guarentee the memory isn't shared
577  */
578  mpoly = lwgeom_as_multi(tmp);
579  clone = lwgeom_clone_deep(mpoly);
580  lwgeom_free(tmp);
581  lwgeom_free(mpoly);
582  mpoly = clone;
583 
584  RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
585 
586 #if POSTGIS_DEBUG_LEVEL > 3
587  {
588  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
589  RASTER_DEBUGF(4, "after multi = %s", wkt);
590  rtdealloc(wkt);
591  }
592 #endif
593  }
594 
595 #if POSTGIS_DEBUG_LEVEL > 3
596  {
597  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
598  RASTER_DEBUGF(4, "returning geometry = %s", wkt);
599  rtdealloc(wkt);
600  }
601 #endif
602 
603  *surface = lwgeom_as_lwmpoly(mpoly);
604  return ES_NONE;
605  }
606 
607  return ES_NONE;
608 }
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:991
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:655
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1050
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
tuple band
Definition: ovdump.py:57
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:433
LWGEOM * lwgeom_as_multi(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate MULTI* type.
Definition: lwgeom.c:284
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:239
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1338
LWMPOLY * lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:170
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition: rt_geometry.c:822
#define WKT_ISO
Definition: liblwgeom.h:1939
LWGEOM * lwgeom_make_valid(LWGEOM *geom)
Attempts to make an invalid geometries valid w/out losing points.
void lwgeom_geos_error(const char *fmt,...)
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:79
tuple nband
Definition: pixval.py:52
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:959
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
void rtinfo(const char *fmt,...)
Definition: rt_context.c:211
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:311
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:356
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
void rtdealloc(void *mem)
Definition: rt_context.c:186
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:307
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:581
void lwgeom_set_srid(LWGEOM *geom, int srid)
Set the SRID on an LWGEOM For collections, only the parent gets an SRID, all the children get SRID_UN...

Here is the call graph for this function:

Here is the caller graph for this function: