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

A "geomval" value is a complex type composed of a geometry in LWPOLY representation (one for each group of pixel sharing the same value) and the value associated with this geometry.

Parameters
raster: the raster to get info from.
nband: the band to polygonize. 0-based
exclude_nodata_value: if non-zero, ignore nodata values to check for pixels with value
Returns
A set of "geomval" values, one for each group of pixels sharing the same value for the provided band. The returned values are LWPOLY geometries.

From GDALPolygonize function header: "Polygon features will be created on the output layer, with polygon geometries representing the polygons". So,the WKB geometry type should be "wkbPolygon"

Create a new field in the layer, to store the px value

We don't need a raster mask band. Each band has a nodata value.

Optimization: Apply a OGR SQL filter to the layer to select the features different from NODATA value.

Thanks to David Zwarg.

Definition at line 6175 of file rt_api.c.

References ovdump::band, FALSE, rt_geomval_t::geom, LW_PARSER_CHECK_NONE, LWGEOM2GEOS(), lwgeom_as_lwpoly(), lwgeom_free(), lwgeom_from_wkb(), lwgeom_geos_error(), lwgeom_make_valid(), lwgeom_set_srid(), lwgeom_to_wkb(), lwgeom_to_wkt(), lwnotice(), RASTER_DEBUG, RASTER_DEBUGF, rt_band_get_hasnodata_flag(), rt_band_get_isnodata_flag(), rt_band_get_nodata(), rt_raster_get_band(), rt_raster_get_srid(), rt_raster_to_gdal_mem(), rt_util_gdal_register_all(), rtalloc(), rtdealloc(), rterror(), rtwarn(), TRUE, rt_geomval_t::val, WKB_ISO, WKB_NDR, and WKT_ISO.

Referenced by RASTER_dumpAsPolygons(), rt_raster_surface(), and test_gdal_polygonize().

6179  {
6180  CPLErr cplerr = CE_None;
6181  char *pszQuery;
6182  long j;
6183  OGRSFDriverH ogr_drv = NULL;
6184  GDALDriverH gdal_drv = NULL;
6185  int destroy_gdal_drv = 0;
6186  GDALDatasetH memdataset = NULL;
6187  GDALRasterBandH gdal_band = NULL;
6188  OGRDataSourceH memdatasource = NULL;
6189  rt_geomval pols = NULL;
6190  OGRLayerH hLayer = NULL;
6191  OGRFeatureH hFeature = NULL;
6192  OGRGeometryH hGeom = NULL;
6193  OGRFieldDefnH hFldDfn = NULL;
6194  unsigned char *wkb = NULL;
6195  int wkbsize = 0;
6196  LWGEOM *lwgeom = NULL;
6197  int nFeatureCount = 0;
6198  rt_band band = NULL;
6199  int iPixVal = -1;
6200  double dValue = 0.0;
6201  int iBandHasNodataValue = FALSE;
6202  double dBandNoData = 0.0;
6203 
6204  /* for checking that a geometry is valid */
6205  GEOSGeometry *ggeom = NULL;
6206  int isValid;
6207  LWGEOM *lwgeomValid = NULL;
6208 
6209  uint32_t bandNums[1] = {nband};
6210  int excludeNodataValues[1] = {exclude_nodata_value};
6211 
6212  /* checks */
6213  assert(NULL != raster);
6214  assert(NULL != pnElements);
6215 
6216  RASTER_DEBUG(2, "In rt_raster_gdal_polygonize");
6217 
6218  *pnElements = 0;
6219 
6220  /*******************************
6221  * Get band
6222  *******************************/
6223  band = rt_raster_get_band(raster, nband);
6224  if (NULL == band) {
6225  rterror("rt_raster_gdal_polygonize: Error getting band %d from raster", nband);
6226  return NULL;
6227  }
6228 
6229  if (exclude_nodata_value) {
6230 
6231  /* band is NODATA */
6232  if (rt_band_get_isnodata_flag(band)) {
6233  RASTER_DEBUG(3, "Band is NODATA. Returning null");
6234  *pnElements = 0;
6235  return NULL;
6236  }
6237 
6238  iBandHasNodataValue = rt_band_get_hasnodata_flag(band);
6239  if (iBandHasNodataValue)
6240  rt_band_get_nodata(band, &dBandNoData);
6241  else
6242  exclude_nodata_value = FALSE;
6243  }
6244 
6245  /*****************************************************
6246  * Convert raster to GDAL MEM dataset
6247  *****************************************************/
6248  memdataset = rt_raster_to_gdal_mem(raster, NULL, bandNums, excludeNodataValues, 1, &gdal_drv, &destroy_gdal_drv);
6249  if (NULL == memdataset) {
6250  rterror("rt_raster_gdal_polygonize: Couldn't convert raster to GDAL MEM dataset");
6251  return NULL;
6252  }
6253 
6254  /*****************************
6255  * Register ogr mem driver
6256  *****************************/
6257 #ifdef GDAL_DCAP_RASTER
6258  /* As for GDAL 2.0, OGRRegisterAll() is an alias for GDALAllRegister() */
6260 #else
6261  OGRRegisterAll();
6262 #endif
6263 
6264  RASTER_DEBUG(3, "creating OGR MEM vector");
6265 
6266  /*****************************************************
6267  * Create an OGR in-memory vector for layers
6268  *****************************************************/
6269  ogr_drv = OGRGetDriverByName("Memory");
6270  memdatasource = OGR_Dr_CreateDataSource(ogr_drv, "", NULL);
6271  if (NULL == memdatasource) {
6272  rterror("rt_raster_gdal_polygonize: Couldn't create a OGR Datasource to store pols");
6273  GDALClose(memdataset);
6274  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
6275  return NULL;
6276  }
6277 
6278  /* Can MEM driver create new layers? */
6279  if (!OGR_DS_TestCapability(memdatasource, ODsCCreateLayer)) {
6280  rterror("rt_raster_gdal_polygonize: MEM driver can't create new layers, aborting");
6281 
6282  /* xxx jorgearevalo: what should we do now? */
6283  GDALClose(memdataset);
6284  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
6285  OGRReleaseDataSource(memdatasource);
6286 
6287  return NULL;
6288  }
6289 
6290  RASTER_DEBUG(3, "polygonizying GDAL MEM raster band");
6291 
6292  /*****************************
6293  * Polygonize the raster band
6294  *****************************/
6295 
6301  hLayer = OGR_DS_CreateLayer(memdatasource, "PolygonizedLayer", NULL, wkbPolygon, NULL);
6302 
6303  if (NULL == hLayer) {
6304  rterror("rt_raster_gdal_polygonize: Couldn't create layer to store polygons");
6305 
6306  GDALClose(memdataset);
6307  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
6308  OGRReleaseDataSource(memdatasource);
6309 
6310  return NULL;
6311  }
6312 
6317  /* First, create a field definition to create the field */
6318  hFldDfn = OGR_Fld_Create("PixelValue", OFTReal);
6319 
6320  /* Second, create the field */
6321  if (OGR_L_CreateField(hLayer, hFldDfn, TRUE) != OGRERR_NONE) {
6322  rtwarn("Couldn't create a field in OGR Layer. The polygons generated won't be able to store the pixel value");
6323  iPixVal = -1;
6324  }
6325  else {
6326  /* Index to the new field created in the layer */
6327  iPixVal = 0;
6328  }
6329 
6330  /* Get GDAL raster band */
6331  gdal_band = GDALGetRasterBand(memdataset, 1);
6332  if (NULL == gdal_band) {
6333  rterror("rt_raster_gdal_polygonize: Couldn't get GDAL band to polygonize");
6334 
6335  GDALClose(memdataset);
6336  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
6337  OGR_Fld_Destroy(hFldDfn);
6338  OGR_DS_DeleteLayer(memdatasource, 0);
6339  OGRReleaseDataSource(memdatasource);
6340 
6341  return NULL;
6342  }
6343 
6347 #ifdef GDALFPOLYGONIZE
6348  cplerr = GDALFPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
6349 #else
6350  cplerr = GDALPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
6351 #endif
6352 
6353  if (cplerr != CE_None) {
6354  rterror("rt_raster_gdal_polygonize: Could not polygonize GDAL band");
6355 
6356  GDALClose(memdataset);
6357  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
6358  OGR_Fld_Destroy(hFldDfn);
6359  OGR_DS_DeleteLayer(memdatasource, 0);
6360  OGRReleaseDataSource(memdatasource);
6361 
6362  return NULL;
6363  }
6364 
6371  if (iBandHasNodataValue) {
6372  pszQuery = (char *) rtalloc(50 * sizeof (char));
6373  sprintf(pszQuery, "PixelValue != %f", dBandNoData );
6374  OGRErr e = OGR_L_SetAttributeFilter(hLayer, pszQuery);
6375  if (e != OGRERR_NONE) {
6376  rtwarn("Error filtering NODATA values for band. All values will be treated as data values");
6377  }
6378  }
6379  else {
6380  pszQuery = NULL;
6381  }
6382 
6383  /*********************************************************************
6384  * Transform OGR layers to WKB polygons
6385  * XXX jorgearevalo: GDALPolygonize does not set the coordinate system
6386  * on the output layer. Application code should do this when the layer
6387  * is created, presumably matching the raster coordinate system.
6388  *********************************************************************/
6389  nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE);
6390 
6391  /* Allocate memory for pols */
6392  pols = (rt_geomval) rtalloc(nFeatureCount * sizeof(struct rt_geomval_t));
6393 
6394  if (NULL == pols) {
6395  rterror("rt_raster_gdal_polygonize: Could not allocate memory for geomval set");
6396 
6397  GDALClose(memdataset);
6398  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
6399  OGR_Fld_Destroy(hFldDfn);
6400  OGR_DS_DeleteLayer(memdatasource, 0);
6401  if (NULL != pszQuery)
6402  rtdealloc(pszQuery);
6403  OGRReleaseDataSource(memdatasource);
6404 
6405  return NULL;
6406  }
6407 
6408  /* initialize GEOS */
6409  initGEOS(lwnotice, lwgeom_geos_error);
6410 
6411  RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount);
6412 
6413  /* Reset feature reading to start in the first feature */
6414  OGR_L_ResetReading(hLayer);
6415 
6416  for (j = 0; j < nFeatureCount; j++) {
6417  hFeature = OGR_L_GetNextFeature(hLayer);
6418  dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal);
6419 
6420  hGeom = OGR_F_GetGeometryRef(hFeature);
6421  wkbsize = OGR_G_WkbSize(hGeom);
6422 
6423  /* allocate wkb buffer */
6424  wkb = rtalloc(sizeof(unsigned char) * wkbsize);
6425  if (wkb == NULL) {
6426  rterror("rt_raster_gdal_polygonize: Could not allocate memory for WKB buffer");
6427 
6428  OGR_F_Destroy(hFeature);
6429  GDALClose(memdataset);
6430  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
6431  OGR_Fld_Destroy(hFldDfn);
6432  OGR_DS_DeleteLayer(memdatasource, 0);
6433  if (NULL != pszQuery)
6434  rtdealloc(pszQuery);
6435  OGRReleaseDataSource(memdatasource);
6436 
6437  return NULL;
6438  }
6439 
6440  /* export WKB with LSB byte order */
6441  OGR_G_ExportToWkb(hGeom, wkbNDR, wkb);
6442 
6443  /* convert WKB to LWGEOM */
6444  lwgeom = lwgeom_from_wkb(wkb, wkbsize, LW_PARSER_CHECK_NONE);
6445 
6446 #if POSTGIS_DEBUG_LEVEL > 3
6447  {
6448  char *wkt = NULL;
6449  OGR_G_ExportToWkt(hGeom, &wkt);
6450  RASTER_DEBUGF(4, "GDAL wkt = %s", wkt);
6451  CPLFree(wkt);
6452 
6453  d_print_binary_hex("GDAL wkb", wkb, wkbsize);
6454  }
6455 #endif
6456 
6457  /* cleanup unnecessary stuff */
6458  rtdealloc(wkb);
6459  wkb = NULL;
6460  wkbsize = 0;
6461 
6462  OGR_F_Destroy(hFeature);
6463 
6464  /* specify SRID */
6465  lwgeom_set_srid(lwgeom, rt_raster_get_srid(raster));
6466 
6467  /*
6468  is geometry valid?
6469  if not, try to make valid
6470  */
6471  do {
6472 #if POSTGIS_GEOS_VERSION < 33
6473  /* nothing can be done if the geometry was invalid if GEOS < 3.3 */
6474  break;
6475 #endif
6476 
6477  ggeom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom);
6478  if (ggeom == NULL) {
6479  rtwarn("Cannot test geometry for validity");
6480  break;
6481  }
6482 
6483  isValid = GEOSisValid(ggeom);
6484 
6485  GEOSGeom_destroy(ggeom);
6486  ggeom = NULL;
6487 
6488  /* geometry is valid */
6489  if (isValid)
6490  break;
6491 
6492  RASTER_DEBUG(3, "fixing invalid geometry");
6493 
6494  /* make geometry valid */
6495  lwgeomValid = lwgeom_make_valid(lwgeom);
6496  if (lwgeomValid == NULL) {
6497  rtwarn("Cannot fix invalid geometry");
6498  break;
6499  }
6500 
6501  lwgeom_free(lwgeom);
6502  lwgeom = lwgeomValid;
6503  }
6504  while (0);
6505 
6506  /* save lwgeom */
6507  pols[j].geom = lwgeom_as_lwpoly(lwgeom);
6508 
6509 #if POSTGIS_DEBUG_LEVEL > 3
6510  {
6511  char *wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, NULL);
6512  RASTER_DEBUGF(4, "LWGEOM wkt = %s", wkt);
6513  rtdealloc(wkt);
6514 
6515  size_t lwwkbsize = 0;
6516  uint8_t *lwwkb = lwgeom_to_wkb(lwgeom, WKB_ISO | WKB_NDR, &lwwkbsize);
6517  if (lwwkbsize) {
6518  d_print_binary_hex("LWGEOM wkb", lwwkb, lwwkbsize);
6519  rtdealloc(lwwkb);
6520  }
6521  }
6522 #endif
6523 
6524  /* set pixel value */
6525  pols[j].val = dValue;
6526  }
6527 
6528  *pnElements = nFeatureCount;
6529 
6530  RASTER_DEBUG(3, "destroying GDAL MEM raster");
6531  GDALClose(memdataset);
6532  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
6533 
6534  RASTER_DEBUG(3, "destroying OGR MEM vector");
6535  OGR_Fld_Destroy(hFldDfn);
6536  OGR_DS_DeleteLayer(memdatasource, 0);
6537  if (NULL != pszQuery) rtdealloc(pszQuery);
6538  OGRReleaseDataSource(memdatasource);
6539 
6540  return pols;
6541 }
int rt_util_gdal_register_all(int force_register_all)
Definition: rt_api.c:445
void rtdealloc(void *mem)
Definition: rt_api.c:882
#define WKB_NDR
Definition: liblwgeom.h:1770
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
GDALDatasetH rt_raster_to_gdal_mem(rt_raster raster, const char *srs, uint32_t *bandNums, int *excludeNodataValues, int count, GDALDriverH *rtn_drv, int *destroy_rtn_drv)
Return GDAL dataset using GDAL MEM driver from raster.
Definition: rt_api.c:8989
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:125
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom)
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
LWPOLY * geom
Definition: rt_api.h:2270
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:1706
#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
struct rt_geomval_t * rt_geomval
Definition: rt_api.h:136
void rtwarn(const char *fmt,...)
Definition: rt_api.c:920
double val
Definition: rt_api.h:2271
#define WKT_ISO
Definition: liblwgeom.h:1776
uint8_t * lwgeom_to_wkb(const LWGEOM *geom, uint8_t variant, size_t *size_out)
Convert LWGEOM to a char* in WKB format.
Definition: lwout_wkb.c:692
#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,...)
tuple nband
Definition: pixval.py:52
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_api.c:3058
#define WKB_ISO
Definition: liblwgeom.h:1767
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
#define FALSE
Definition: dbfopen.c:169
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...
LWGEOM * lwgeom_from_wkb(const uint8_t *wkb, const size_t wkb_size, const char check)
WKB inputs must have a declared size, to prevent malformed WKB from reading off the end of the memory...
Definition: lwin_wkb.c:728
#define TRUE
Definition: dbfopen.c:170
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: