PostGIS  2.3.8dev-r@@SVN_REVISION@@

◆ 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.

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  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 }
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:1004
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:669
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1063
band
Definition: ovdump.py:57
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
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:446
LWGEOM * lwgeom_as_multi(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate MULTI* type.
Definition: lwgeom.c:297
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:252
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:183
nband
Definition: pixval.py:52
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster&#39;s convex hull.
Definition: rt_geometry.c:803
#define WKT_ISO
Definition: liblwgeom.h:2055
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:152
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_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&#39;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: