PostGIS  2.2.8dev-r@@SVN_REVISION@@

◆ 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 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 959 of file rt_geometry.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(), 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(), rtinfo(), rtwarn(), TRUE, rt_geomval_t::val, WKB_ISO, WKB_NDR, and WKT_ISO.

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

963  {
964  CPLErr cplerr = CE_None;
965  char *pszQuery;
966  long j;
967  OGRSFDriverH ogr_drv = NULL;
968  GDALDriverH gdal_drv = NULL;
969  int destroy_gdal_drv = 0;
970  GDALDatasetH memdataset = NULL;
971  GDALRasterBandH gdal_band = NULL;
972  OGRDataSourceH memdatasource = NULL;
973  rt_geomval pols = NULL;
974  OGRLayerH hLayer = NULL;
975  OGRFeatureH hFeature = NULL;
976  OGRGeometryH hGeom = NULL;
977  OGRFieldDefnH hFldDfn = NULL;
978  unsigned char *wkb = NULL;
979  int wkbsize = 0;
980  LWGEOM *lwgeom = NULL;
981  int nFeatureCount = 0;
982  rt_band band = NULL;
983  int iPixVal = -1;
984  double dValue = 0.0;
985  int iBandHasNodataValue = FALSE;
986  double dBandNoData = 0.0;
987 
988  /* for checking that a geometry is valid */
989  GEOSGeometry *ggeom = NULL;
990  int isValid;
991  LWGEOM *lwgeomValid = NULL;
992 
993  uint32_t bandNums[1] = {nband};
994  int excludeNodataValues[1] = {exclude_nodata_value};
995 
996  /* checks */
997  assert(NULL != raster);
998  assert(NULL != pnElements);
999 
1000  RASTER_DEBUG(2, "In rt_raster_gdal_polygonize");
1001 
1002  *pnElements = 0;
1003 
1004  /*******************************
1005  * Get band
1006  *******************************/
1007  band = rt_raster_get_band(raster, nband);
1008  if (NULL == band) {
1009  rterror("rt_raster_gdal_polygonize: Error getting band %d from raster", nband);
1010  return NULL;
1011  }
1012 
1013  if (exclude_nodata_value) {
1014 
1015  /* band is NODATA */
1016  if (rt_band_get_isnodata_flag(band)) {
1017  RASTER_DEBUG(3, "Band is NODATA. Returning null");
1018  *pnElements = 0;
1019  return NULL;
1020  }
1021 
1022  iBandHasNodataValue = rt_band_get_hasnodata_flag(band);
1023  if (iBandHasNodataValue)
1024  rt_band_get_nodata(band, &dBandNoData);
1025  else
1026  exclude_nodata_value = FALSE;
1027  }
1028 
1029  /*****************************************************
1030  * Convert raster to GDAL MEM dataset
1031  *****************************************************/
1032  memdataset = rt_raster_to_gdal_mem(raster, NULL, bandNums, excludeNodataValues, 1, &gdal_drv, &destroy_gdal_drv);
1033  if (NULL == memdataset) {
1034  rterror("rt_raster_gdal_polygonize: Couldn't convert raster to GDAL MEM dataset");
1035  return NULL;
1036  }
1037 
1038  /*****************************
1039  * Register ogr mem driver
1040  *****************************/
1041 #ifdef GDAL_DCAP_RASTER
1042  /* in GDAL 2.0, OGRRegisterAll() is an alias to GDALAllRegister() */
1044 #else
1045  OGRRegisterAll();
1046 #endif
1047 
1048  RASTER_DEBUG(3, "creating OGR MEM vector");
1049 
1050  /*****************************************************
1051  * Create an OGR in-memory vector for layers
1052  *****************************************************/
1053  ogr_drv = OGRGetDriverByName("Memory");
1054  memdatasource = OGR_Dr_CreateDataSource(ogr_drv, "", NULL);
1055  if (NULL == memdatasource) {
1056  rterror("rt_raster_gdal_polygonize: Couldn't create a OGR Datasource to store pols");
1057  GDALClose(memdataset);
1058  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1059  return NULL;
1060  }
1061 
1062  /* Can MEM driver create new layers? */
1063  if (!OGR_DS_TestCapability(memdatasource, ODsCCreateLayer)) {
1064  rterror("rt_raster_gdal_polygonize: MEM driver can't create new layers, aborting");
1065 
1066  /* xxx jorgearevalo: what should we do now? */
1067  GDALClose(memdataset);
1068  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1069  OGRReleaseDataSource(memdatasource);
1070 
1071  return NULL;
1072  }
1073 
1074  RASTER_DEBUG(3, "polygonizying GDAL MEM raster band");
1075 
1076  /*****************************
1077  * Polygonize the raster band
1078  *****************************/
1079 
1085  hLayer = OGR_DS_CreateLayer(memdatasource, "PolygonizedLayer", NULL, wkbPolygon, NULL);
1086 
1087  if (NULL == hLayer) {
1088  rterror("rt_raster_gdal_polygonize: Couldn't create layer to store polygons");
1089 
1090  GDALClose(memdataset);
1091  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1092  OGRReleaseDataSource(memdatasource);
1093 
1094  return NULL;
1095  }
1096 
1101  /* First, create a field definition to create the field */
1102  hFldDfn = OGR_Fld_Create("PixelValue", OFTReal);
1103 
1104  /* Second, create the field */
1105  if (OGR_L_CreateField(hLayer, hFldDfn, TRUE) != OGRERR_NONE) {
1106  rtwarn("Couldn't create a field in OGR Layer. The polygons generated won't be able to store the pixel value");
1107  iPixVal = -1;
1108  }
1109  else {
1110  /* Index to the new field created in the layer */
1111  iPixVal = 0;
1112  }
1113 
1114  /* Get GDAL raster band */
1115  gdal_band = GDALGetRasterBand(memdataset, 1);
1116  if (NULL == gdal_band) {
1117  rterror("rt_raster_gdal_polygonize: Couldn't get GDAL band to polygonize");
1118 
1119  GDALClose(memdataset);
1120  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1121  OGR_Fld_Destroy(hFldDfn);
1122  OGR_DS_DeleteLayer(memdatasource, 0);
1123  OGRReleaseDataSource(memdatasource);
1124 
1125  return NULL;
1126  }
1127 
1131 #ifdef GDALFPOLYGONIZE
1132  cplerr = GDALFPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
1133 #else
1134  cplerr = GDALPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
1135 #endif
1136 
1137  if (cplerr != CE_None) {
1138  rterror("rt_raster_gdal_polygonize: Could not polygonize GDAL band");
1139 
1140  GDALClose(memdataset);
1141  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1142  OGR_Fld_Destroy(hFldDfn);
1143  OGR_DS_DeleteLayer(memdatasource, 0);
1144  OGRReleaseDataSource(memdatasource);
1145 
1146  return NULL;
1147  }
1148 
1155  if (iBandHasNodataValue) {
1156  pszQuery = (char *) rtalloc(50 * sizeof (char));
1157  sprintf(pszQuery, "PixelValue != %f", dBandNoData );
1158  OGRErr e = OGR_L_SetAttributeFilter(hLayer, pszQuery);
1159  if (e != OGRERR_NONE) {
1160  rtwarn("Error filtering NODATA values for band. All values will be treated as data values");
1161  }
1162  }
1163  else {
1164  pszQuery = NULL;
1165  }
1166 
1167  /*********************************************************************
1168  * Transform OGR layers to WKB polygons
1169  * XXX jorgearevalo: GDALPolygonize does not set the coordinate system
1170  * on the output layer. Application code should do this when the layer
1171  * is created, presumably matching the raster coordinate system.
1172  *********************************************************************/
1173  nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE);
1174 
1175  /* Allocate memory for pols */
1176  pols = (rt_geomval) rtalloc(nFeatureCount * sizeof(struct rt_geomval_t));
1177 
1178  if (NULL == pols) {
1179  rterror("rt_raster_gdal_polygonize: Could not allocate memory for geomval set");
1180 
1181  GDALClose(memdataset);
1182  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1183  OGR_Fld_Destroy(hFldDfn);
1184  OGR_DS_DeleteLayer(memdatasource, 0);
1185  if (NULL != pszQuery)
1186  rtdealloc(pszQuery);
1187  OGRReleaseDataSource(memdatasource);
1188 
1189  return NULL;
1190  }
1191 
1192  /* initialize GEOS */
1193  initGEOS(rtinfo, lwgeom_geos_error);
1194 
1195  RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount);
1196 
1197  /* Reset feature reading to start in the first feature */
1198  OGR_L_ResetReading(hLayer);
1199 
1200  for (j = 0; j < nFeatureCount; j++) {
1201  hFeature = OGR_L_GetNextFeature(hLayer);
1202  dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal);
1203 
1204  hGeom = OGR_F_GetGeometryRef(hFeature);
1205  wkbsize = OGR_G_WkbSize(hGeom);
1206 
1207  /* allocate wkb buffer */
1208  wkb = rtalloc(sizeof(unsigned char) * wkbsize);
1209  if (wkb == NULL) {
1210  rterror("rt_raster_gdal_polygonize: Could not allocate memory for WKB buffer");
1211 
1212  OGR_F_Destroy(hFeature);
1213  GDALClose(memdataset);
1214  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1215  OGR_Fld_Destroy(hFldDfn);
1216  OGR_DS_DeleteLayer(memdatasource, 0);
1217  if (NULL != pszQuery)
1218  rtdealloc(pszQuery);
1219  OGRReleaseDataSource(memdatasource);
1220 
1221  return NULL;
1222  }
1223 
1224  /* export WKB with LSB byte order */
1225  OGR_G_ExportToWkb(hGeom, wkbNDR, wkb);
1226 
1227  /* convert WKB to LWGEOM */
1228  lwgeom = lwgeom_from_wkb(wkb, wkbsize, LW_PARSER_CHECK_NONE);
1229 
1230 #if POSTGIS_DEBUG_LEVEL > 3
1231  {
1232  char *wkt = NULL;
1233  OGR_G_ExportToWkt(hGeom, &wkt);
1234  RASTER_DEBUGF(4, "GDAL wkt = %s", wkt);
1235  CPLFree(wkt);
1236 
1237  d_print_binary_hex("GDAL wkb", wkb, wkbsize);
1238  }
1239 #endif
1240 
1241  /* cleanup unnecessary stuff */
1242  rtdealloc(wkb);
1243  wkb = NULL;
1244  wkbsize = 0;
1245 
1246  OGR_F_Destroy(hFeature);
1247 
1248  /* specify SRID */
1249  lwgeom_set_srid(lwgeom, rt_raster_get_srid(raster));
1250 
1251  /*
1252  is geometry valid?
1253  if not, try to make valid
1254  */
1255  do {
1256 #if POSTGIS_GEOS_VERSION < 33
1257  /* nothing can be done if the geometry was invalid if GEOS < 3.3 */
1258  break;
1259 #endif
1260 
1261  ggeom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom, 0);
1262  if (ggeom == NULL) {
1263  rtwarn("Cannot test geometry for validity");
1264  break;
1265  }
1266 
1267  isValid = GEOSisValid(ggeom);
1268 
1269  GEOSGeom_destroy(ggeom);
1270  ggeom = NULL;
1271 
1272  /* geometry is valid */
1273  if (isValid)
1274  break;
1275 
1276  RASTER_DEBUG(3, "fixing invalid geometry");
1277 
1278  /* make geometry valid */
1279  lwgeomValid = lwgeom_make_valid(lwgeom);
1280  if (lwgeomValid == NULL) {
1281  rtwarn("Cannot fix invalid geometry");
1282  break;
1283  }
1284 
1285  lwgeom_free(lwgeom);
1286  lwgeom = lwgeomValid;
1287  }
1288  while (0);
1289 
1290  /* save lwgeom */
1291  pols[j].geom = lwgeom_as_lwpoly(lwgeom);
1292 
1293 #if POSTGIS_DEBUG_LEVEL > 3
1294  {
1295  char *wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, NULL);
1296  RASTER_DEBUGF(4, "LWGEOM wkt = %s", wkt);
1297  rtdealloc(wkt);
1298 
1299  size_t lwwkbsize = 0;
1300  uint8_t *lwwkb = lwgeom_to_wkb(lwgeom, WKB_ISO | WKB_NDR, &lwwkbsize);
1301  if (lwwkbsize) {
1302  d_print_binary_hex("LWGEOM wkb", lwwkb, lwwkbsize);
1303  rtdealloc(lwwkb);
1304  }
1305  }
1306 #endif
1307 
1308  /* set pixel value */
1309  pols[j].val = dValue;
1310  }
1311 
1312  *pnElements = nFeatureCount;
1313 
1314  RASTER_DEBUG(3, "destroying GDAL MEM raster");
1315  GDALClose(memdataset);
1316  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1317 
1318  RASTER_DEBUG(3, "destroying OGR MEM vector");
1319  OGR_Fld_Destroy(hFldDfn);
1320  OGR_DS_DeleteLayer(memdatasource, 0);
1321  if (NULL != pszQuery) rtdealloc(pszQuery);
1322  OGRReleaseDataSource(memdatasource);
1323 
1324  return pols;
1325 }
int rt_util_gdal_register_all(int force_register_all)
Definition: rt_util.c:334
#define WKB_NDR
Definition: liblwgeom.h:1933
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
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_raster.c:1809
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
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:125
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
LWPOLY * geom
Definition: librtcore.h:2316
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:1869
nband
Definition: pixval.py:52
double val
Definition: librtcore.h:2317
#define WKT_ISO
Definition: liblwgeom.h:1939
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:750
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
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)
#define WKB_ISO
Definition: liblwgeom.h:1930
struct rt_geomval_t * rt_geomval
Definition: librtcore.h:161
void rtdealloc(void *mem)
Definition: rt_context.c:186
#define FALSE
Definition: dbfopen.c:168
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:307
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...
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:754
#define TRUE
Definition: dbfopen.c:169
Here is the call graph for this function:
Here is the caller graph for this function: