PostGIS  2.4.9dev-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

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 
1109  /* We don't need a raster mask band. Each band has a nodata value. */
1110  cplerr = GDALFPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
1111 
1112  if (cplerr != CE_None) {
1113  rterror("rt_raster_gdal_polygonize: Could not polygonize GDAL band");
1114 
1115  GDALClose(memdataset);
1116  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1117  OGR_Fld_Destroy(hFldDfn);
1118  OGR_DS_DeleteLayer(memdatasource, 0);
1119  OGRReleaseDataSource(memdatasource);
1120 
1121  return NULL;
1122  }
1123 
1130  if (iBandHasNodataValue) {
1131  pszQuery = (char *) rtalloc(50 * sizeof (char));
1132  sprintf(pszQuery, "PixelValue != %f", dBandNoData );
1133  OGRErr e = OGR_L_SetAttributeFilter(hLayer, pszQuery);
1134  if (e != OGRERR_NONE) {
1135  rtwarn("Error filtering NODATA values for band. All values will be treated as data values");
1136  }
1137  }
1138  else {
1139  pszQuery = NULL;
1140  }
1141 
1142  /*********************************************************************
1143  * Transform OGR layers to WKB polygons
1144  * XXX jorgearevalo: GDALPolygonize does not set the coordinate system
1145  * on the output layer. Application code should do this when the layer
1146  * is created, presumably matching the raster coordinate system.
1147  *********************************************************************/
1148  nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE);
1149 
1150  /* Allocate memory for pols */
1151  pols = (rt_geomval) rtalloc(nFeatureCount * sizeof(struct rt_geomval_t));
1152 
1153  if (NULL == pols) {
1154  rterror("rt_raster_gdal_polygonize: Could not allocate memory for geomval set");
1155 
1156  GDALClose(memdataset);
1157  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1158  OGR_Fld_Destroy(hFldDfn);
1159  OGR_DS_DeleteLayer(memdatasource, 0);
1160  if (NULL != pszQuery)
1161  rtdealloc(pszQuery);
1162  OGRReleaseDataSource(memdatasource);
1163 
1164  return NULL;
1165  }
1166 
1167  /* initialize GEOS */
1168  initGEOS(rtinfo, lwgeom_geos_error);
1169 
1170  RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount);
1171 
1172  /* Reset feature reading to start in the first feature */
1173  OGR_L_ResetReading(hLayer);
1174 
1175  for (j = 0; j < nFeatureCount; j++) {
1176  hFeature = OGR_L_GetNextFeature(hLayer);
1177  dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal);
1178 
1179  hGeom = OGR_F_GetGeometryRef(hFeature);
1180  wkbsize = OGR_G_WkbSize(hGeom);
1181 
1182  /* allocate wkb buffer */
1183  wkb = rtalloc(sizeof(unsigned char) * wkbsize);
1184  if (wkb == NULL) {
1185  rterror("rt_raster_gdal_polygonize: Could not allocate memory for WKB buffer");
1186 
1187  OGR_F_Destroy(hFeature);
1188  GDALClose(memdataset);
1189  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1190  OGR_Fld_Destroy(hFldDfn);
1191  OGR_DS_DeleteLayer(memdatasource, 0);
1192  if (NULL != pszQuery)
1193  rtdealloc(pszQuery);
1194  OGRReleaseDataSource(memdatasource);
1195 
1196  return NULL;
1197  }
1198 
1199  /* export WKB with LSB byte order */
1200  OGR_G_ExportToWkb(hGeom, wkbNDR, wkb);
1201 
1202  /* convert WKB to LWGEOM */
1203  lwgeom = lwgeom_from_wkb(wkb, wkbsize, LW_PARSER_CHECK_NONE);
1204 
1205 #if POSTGIS_DEBUG_LEVEL > 3
1206  {
1207  char *wkt = NULL;
1208  OGR_G_ExportToWkt(hGeom, &wkt);
1209  RASTER_DEBUGF(4, "GDAL wkt = %s", wkt);
1210  CPLFree(wkt);
1211 
1212  d_print_binary_hex("GDAL wkb", wkb, wkbsize);
1213  }
1214 #endif
1215 
1216  /* cleanup unnecessary stuff */
1217  rtdealloc(wkb);
1218  wkb = NULL;
1219  wkbsize = 0;
1220 
1221  OGR_F_Destroy(hFeature);
1222 
1223  /* specify SRID */
1224  lwgeom_set_srid(lwgeom, rt_raster_get_srid(raster));
1225 
1226  /*
1227  is geometry valid?
1228  if not, try to make valid
1229  */
1230  do {
1231  ggeom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom, 0);
1232  if (ggeom == NULL) {
1233  rtwarn("Cannot test geometry for validity");
1234  break;
1235  }
1236 
1237  isValid = GEOSisValid(ggeom);
1238 
1239  GEOSGeom_destroy(ggeom);
1240  ggeom = NULL;
1241 
1242  /* geometry is valid */
1243  if (isValid)
1244  break;
1245 
1246  RASTER_DEBUG(3, "fixing invalid geometry");
1247 
1248  /* make geometry valid */
1249  lwgeomValid = lwgeom_make_valid(lwgeom);
1250  if (lwgeomValid == NULL) {
1251  rtwarn("Cannot fix invalid geometry");
1252  break;
1253  }
1254 
1255  lwgeom_free(lwgeom);
1256  lwgeom = lwgeomValid;
1257  }
1258  while (0);
1259 
1260  /* save lwgeom */
1261  pols[j].geom = lwgeom_as_lwpoly(lwgeom);
1262 
1263 #if POSTGIS_DEBUG_LEVEL > 3
1264  {
1265  char *wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, NULL);
1266  RASTER_DEBUGF(4, "LWGEOM wkt = %s", wkt);
1267  rtdealloc(wkt);
1268 
1269  size_t lwwkbsize = 0;
1270  uint8_t *lwwkb = lwgeom_to_wkb(lwgeom, WKB_ISO | WKB_NDR, &lwwkbsize);
1271  if (lwwkbsize) {
1272  d_print_binary_hex("LWGEOM wkb", lwwkb, lwwkbsize);
1273  rtdealloc(lwwkb);
1274  }
1275  }
1276 #endif
1277 
1278  /* set pixel value */
1279  pols[j].val = dValue;
1280  }
1281 
1282  *pnElements = nFeatureCount;
1283 
1284  RASTER_DEBUG(3, "destroying GDAL MEM raster");
1285  GDALClose(memdataset);
1286  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1287 
1288  RASTER_DEBUG(3, "destroying OGR MEM vector");
1289  OGR_Fld_Destroy(hFldDfn);
1290  OGR_DS_DeleteLayer(memdatasource, 0);
1291  if (NULL != pszQuery) rtdealloc(pszQuery);
1292  OGRReleaseDataSource(memdatasource);
1293 
1294  return pols;
1295 }
int rt_util_gdal_register_all(int force_register_all)
Definition: rt_util.c:334
#define WKB_NDR
Definition: liblwgeom.h:2077
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:1099
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:174
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
unsigned int uint32_t
Definition: uthash.h:78
LWPOLY * geom
Definition: librtcore.h:2304
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:2013
nband
Definition: pixval.py:52
double val
Definition: librtcore.h:2305
#define WKT_ISO
Definition: liblwgeom.h:2083
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:541
#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
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
#define WKB_ISO
Definition: liblwgeom.h:2074
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: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:772
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: