PostGIS  2.1.10dev-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 13242 of file rt_api.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(), lwnotice(), 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().

13242  {
13243  rt_band band = NULL;
13244  LWGEOM *mpoly = NULL;
13245  LWGEOM *tmp = NULL;
13246  LWGEOM *clone = NULL;
13247  rt_geomval gv = NULL;
13248  int gvcount = 0;
13249  GEOSGeometry *gc = NULL;
13250  GEOSGeometry *gunion = NULL;
13251  GEOSGeometry **geoms = NULL;
13252  int geomscount = 0;
13253  int i = 0;
13254 
13255  assert(surface != NULL);
13256 
13257  /* init *surface to NULL */
13258  *surface = NULL;
13259 
13260  /* raster is empty, surface = NULL */
13261  if (rt_raster_is_empty(raster))
13262  return ES_NONE;
13263 
13264  /* if nband < 0, return the convex hull as a multipolygon */
13265  if (nband < 0) {
13266  /*
13267  lwgeom_as_multi() only does a shallow clone internally
13268  so input and output geometries may share memory
13269  hence the deep clone of the output geometry for returning
13270  is the only way to guarentee the memory isn't shared
13271  */
13272  if (rt_raster_get_convex_hull(raster, &tmp) != ES_NONE) {
13273  rterror("rt_raster_surface: Could not get convex hull of raster");
13274  return ES_ERROR;
13275  }
13276  mpoly = lwgeom_as_multi(tmp);
13277  clone = lwgeom_clone_deep(mpoly);
13278  lwgeom_free(tmp);
13279  lwgeom_free(mpoly);
13280 
13281  *surface = lwgeom_as_lwmpoly(clone);
13282  return ES_NONE;
13283  }
13284  /* check that nband is valid */
13285  else if (nband >= rt_raster_get_num_bands(raster)) {
13286  rterror("rt_raster_surface: The band index %d is invalid", nband);
13287  return ES_ERROR;
13288  }
13289 
13290  /* get band */
13291  band = rt_raster_get_band(raster, nband);
13292  if (band == NULL) {
13293  rterror("rt_raster_surface: Error getting band %d from raster", nband);
13294  return ES_ERROR;
13295  }
13296 
13297  /* band does not have a NODATA flag, return convex hull */
13298  if (!rt_band_get_hasnodata_flag(band)) {
13299  /*
13300  lwgeom_as_multi() only does a shallow clone internally
13301  so input and output geometries may share memory
13302  hence the deep clone of the output geometry for returning
13303  is the only way to guarentee the memory isn't shared
13304  */
13305  if (rt_raster_get_convex_hull(raster, &tmp) != ES_NONE) {
13306  rterror("rt_raster_surface: Could not get convex hull of raster");
13307  return ES_ERROR;
13308  }
13309  mpoly = lwgeom_as_multi(tmp);
13310  clone = lwgeom_clone_deep(mpoly);
13311  lwgeom_free(tmp);
13312  lwgeom_free(mpoly);
13313 
13314  *surface = lwgeom_as_lwmpoly(clone);
13315  return ES_NONE;
13316  }
13317  /* band is NODATA, return NULL */
13318  else if (rt_band_get_isnodata_flag(band)) {
13319  RASTER_DEBUG(3, "Band is NODATA. Returning NULL");
13320  return ES_NONE;
13321  }
13322 
13323  /* initialize GEOS */
13324  initGEOS(lwnotice, lwgeom_geos_error);
13325 
13326  /* use gdal polygonize */
13327  gv = rt_raster_gdal_polygonize(raster, nband, 1, &gvcount);
13328  /* no polygons returned */
13329  if (gvcount < 1) {
13330  RASTER_DEBUG(3, "All pixels of band are NODATA. Returning NULL");
13331  if (gv != NULL) rtdealloc(gv);
13332  return ES_NONE;
13333  }
13334  /* more than 1 polygon */
13335  else if (gvcount > 1) {
13336  /* convert LWPOLY to GEOSGeometry */
13337  geomscount = gvcount;
13338  geoms = rtalloc(sizeof(GEOSGeometry *) * geomscount);
13339  if (geoms == NULL) {
13340  rterror("rt_raster_surface: Could not allocate memory for pixel polygons as GEOSGeometry");
13341  for (i = 0; i < gvcount; i++) lwpoly_free(gv[i].geom);
13342  rtdealloc(gv);
13343  return ES_ERROR;
13344  }
13345  for (i = 0; i < gvcount; i++) {
13346 #if POSTGIS_DEBUG_LEVEL > 3
13347  {
13348  char *wkt = lwgeom_to_wkt(lwpoly_as_lwgeom(gv[i].geom), WKT_ISO, DBL_DIG, NULL);
13349  RASTER_DEBUGF(4, "geom %d = %s", i, wkt);
13350  rtdealloc(wkt);
13351  }
13352 #endif
13353 
13354  geoms[i] = LWGEOM2GEOS(lwpoly_as_lwgeom(gv[i].geom));
13355  lwpoly_free(gv[i].geom);
13356  }
13357  rtdealloc(gv);
13358 
13359  /* create geometry collection */
13360 #if POSTGIS_GEOS_VERSION >= 33
13361  gc = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, geomscount);
13362 #else
13363  gc = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, geomscount);
13364 #endif
13365 
13366  if (gc == NULL) {
13367 #if POSTGIS_GEOS_VERSION >= 33
13368  rterror("rt_raster_surface: Could not create GEOS GEOMETRYCOLLECTION from set of pixel polygons");
13369 #else
13370  rterror("rt_raster_surface: Could not create GEOS MULTIPOLYGON from set of pixel polygons");
13371 #endif
13372 
13373  for (i = 0; i < geomscount; i++)
13374  GEOSGeom_destroy(geoms[i]);
13375  rtdealloc(geoms);
13376  return ES_ERROR;
13377  }
13378 
13379  /* run the union */
13380 #if POSTGIS_GEOS_VERSION >= 33
13381  gunion = GEOSUnaryUnion(gc);
13382 #else
13383  gunion = GEOSUnionCascaded(gc);
13384 #endif
13385  GEOSGeom_destroy(gc);
13386  rtdealloc(geoms);
13387 
13388  if (gunion == NULL) {
13389 #if POSTGIS_GEOS_VERSION >= 33
13390  rterror("rt_raster_surface: Could not union the pixel polygons using GEOSUnaryUnion()");
13391 #else
13392  rterror("rt_raster_surface: Could not union the pixel polygons using GEOSUnionCascaded()");
13393 #endif
13394  return ES_ERROR;
13395  }
13396 
13397  /* convert union result from GEOSGeometry to LWGEOM */
13398  mpoly = GEOS2LWGEOM(gunion, 0);
13399 
13400  /*
13401  is geometry valid?
13402  if not, try to make valid
13403  */
13404  do {
13405  LWGEOM *mpolyValid = NULL;
13406 
13407 #if POSTGIS_GEOS_VERSION < 33
13408  break;
13409 #endif
13410 
13411  if (GEOSisValid(gunion))
13412  break;
13413 
13414  /* make geometry valid */
13415  mpolyValid = lwgeom_make_valid(mpoly);
13416  if (mpolyValid == NULL) {
13417  rtwarn("Cannot fix invalid geometry");
13418  break;
13419  }
13420 
13421  lwgeom_free(mpoly);
13422  mpoly = mpolyValid;
13423  }
13424  while (0);
13425 
13426  GEOSGeom_destroy(gunion);
13427  }
13428  else {
13429  mpoly = lwpoly_as_lwgeom(gv[0].geom);
13430  rtdealloc(gv);
13431 
13432 #if POSTGIS_DEBUG_LEVEL > 3
13433  {
13434  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
13435  RASTER_DEBUGF(4, "geom 0 = %s", wkt);
13436  rtdealloc(wkt);
13437  }
13438 #endif
13439  }
13440 
13441  /* specify SRID */
13442  lwgeom_set_srid(mpoly, rt_raster_get_srid(raster));
13443 
13444  if (mpoly != NULL) {
13445  /* convert to multi */
13446  if (!lwgeom_is_collection(mpoly)) {
13447  tmp = mpoly;
13448 
13449 #if POSTGIS_DEBUG_LEVEL > 3
13450  {
13451  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
13452  RASTER_DEBUGF(4, "before multi = %s", wkt);
13453  rtdealloc(wkt);
13454  }
13455 #endif
13456 
13457  RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
13458 
13459  /*
13460  lwgeom_as_multi() only does a shallow clone internally
13461  so input and output geometries may share memory
13462  hence the deep clone of the output geometry for returning
13463  is the only way to guarentee the memory isn't shared
13464  */
13465  mpoly = lwgeom_as_multi(tmp);
13466  clone = lwgeom_clone_deep(mpoly);
13467  lwgeom_free(tmp);
13468  lwgeom_free(mpoly);
13469  mpoly = clone;
13470 
13471  RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
13472 
13473 #if POSTGIS_DEBUG_LEVEL > 3
13474  {
13475  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
13476  RASTER_DEBUGF(4, "after multi = %s", wkt);
13477  rtdealloc(wkt);
13478  }
13479 #endif
13480  }
13481 
13482 #if POSTGIS_DEBUG_LEVEL > 3
13483  {
13484  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
13485  RASTER_DEBUGF(4, "returning geometry = %s", wkt);
13486  rtdealloc(wkt);
13487  }
13488 #endif
13489 
13490  *surface = lwgeom_as_lwmpoly(mpoly);
13491  return ES_NONE;
13492  }
13493 
13494  return ES_NONE;
13495 }
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_api.c:5677
void rtdealloc(void *mem)
Definition: rt_api.c:882
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:947
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:1006
tuple band
Definition: ovdump.py:57
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:389
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom)
LWGEOM * lwgeom_as_multi(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate MULTI* type.
Definition: lwgeom.c:284
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_api.c:5661
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_api.c:2042
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:239
LWMPOLY * lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:170
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:54
void rtwarn(const char *fmt,...)
Definition: rt_api.c:920
#define WKT_ISO
Definition: liblwgeom.h:1776
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
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 lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:79
tuple nband
Definition: pixval.py:52
void * rtalloc(size_t size)
Raster core memory management functions.
Definition: rt_api.c:867
void rterror(const char *fmt,...)
Raster core error and info handlers.
Definition: rt_api.c:895
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_api.c:5686
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_api.c:8550
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
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...
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition: rt_api.c:6556
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_api.c:6175
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_api.c:2002

Here is the call graph for this function:

Here is the caller graph for this function: