PostGIS  3.0.0dev-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 940 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().

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