PostGIS 3.0.6dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

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

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() */
1024 rt_util_gdal_register_all(0);
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 size_t sz = 50 * sizeof (char);
1132 pszQuery = (char *) rtalloc(sz);
1133 snprintf(pszQuery, sz, "PixelValue != %f", dBandNoData );
1134 OGRErr e = OGR_L_SetAttributeFilter(hLayer, pszQuery);
1135 if (e != OGRERR_NONE) {
1136 rtwarn("Error filtering NODATA values for band. All values will be treated as data values");
1137 }
1138 }
1139 else {
1140 pszQuery = NULL;
1141 }
1142
1143 /*********************************************************************
1144 * Transform OGR layers to WKB polygons
1145 * XXX jorgearevalo: GDALPolygonize does not set the coordinate system
1146 * on the output layer. Application code should do this when the layer
1147 * is created, presumably matching the raster coordinate system.
1148 *********************************************************************/
1149 nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE);
1150
1151 /* Allocate memory for pols */
1152 pols = (rt_geomval) rtalloc(nFeatureCount * sizeof(struct rt_geomval_t));
1153
1154 if (NULL == pols) {
1155 rterror("rt_raster_gdal_polygonize: Could not allocate memory for geomval set");
1156
1157 GDALClose(memdataset);
1158 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1159 OGR_Fld_Destroy(hFldDfn);
1160 OGR_DS_DeleteLayer(memdatasource, 0);
1161 if (NULL != pszQuery)
1162 rtdealloc(pszQuery);
1163 OGRReleaseDataSource(memdatasource);
1164
1165 return NULL;
1166 }
1167
1168 /* initialize GEOS */
1169 initGEOS(rtinfo, lwgeom_geos_error);
1170
1171 RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount);
1172
1173 /* Reset feature reading to start in the first feature */
1174 OGR_L_ResetReading(hLayer);
1175
1176 for (j = 0; j < nFeatureCount; j++) {
1177 hFeature = OGR_L_GetNextFeature(hLayer);
1178 dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal);
1179
1180 hGeom = OGR_F_GetGeometryRef(hFeature);
1181 wkbsize = OGR_G_WkbSize(hGeom);
1182
1183 /* allocate wkb buffer */
1184 wkb = rtalloc(sizeof(unsigned char) * wkbsize);
1185 if (wkb == NULL) {
1186 rterror("rt_raster_gdal_polygonize: Could not allocate memory for WKB buffer");
1187
1188 OGR_F_Destroy(hFeature);
1189 GDALClose(memdataset);
1190 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1191 OGR_Fld_Destroy(hFldDfn);
1192 OGR_DS_DeleteLayer(memdatasource, 0);
1193 if (NULL != pszQuery)
1194 rtdealloc(pszQuery);
1195 OGRReleaseDataSource(memdatasource);
1196
1197 return NULL;
1198 }
1199
1200 /* export WKB with LSB byte order */
1201 OGR_G_ExportToWkb(hGeom, wkbNDR, wkb);
1202
1203 /* convert WKB to LWGEOM */
1204 lwgeom = lwgeom_from_wkb(wkb, wkbsize, LW_PARSER_CHECK_NONE);
1205 if (!lwgeom) rterror("%s: invalid wkb", __func__);
1206
1207#if POSTGIS_DEBUG_LEVEL > 3
1208 {
1209 char *wkt = NULL;
1210 OGR_G_ExportToWkt(hGeom, &wkt);
1211 RASTER_DEBUGF(4, "GDAL wkt = %s", wkt);
1212 CPLFree(wkt);
1213
1214 d_print_binary_hex("GDAL wkb", wkb, wkbsize);
1215 }
1216#endif
1217
1218 /* cleanup unnecessary stuff */
1219 rtdealloc(wkb);
1220 wkb = NULL;
1221 wkbsize = 0;
1222
1223 OGR_F_Destroy(hFeature);
1224
1225 /* specify SRID */
1226 lwgeom_set_srid(lwgeom, rt_raster_get_srid(raster));
1227
1228 /*
1229 is geometry valid?
1230 if not, try to make valid
1231 */
1232 do {
1233 ggeom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom, 0);
1234 if (ggeom == NULL) {
1235 rtwarn("Cannot test geometry for validity");
1236 break;
1237 }
1238
1239 isValid = GEOSisValid(ggeom);
1240
1241 GEOSGeom_destroy(ggeom);
1242 ggeom = NULL;
1243
1244 /* geometry is valid */
1245 if (isValid)
1246 break;
1247
1248 RASTER_DEBUG(3, "fixing invalid geometry");
1249
1250 /* make geometry valid */
1251 lwgeomValid = lwgeom_make_valid(lwgeom);
1252 if (lwgeomValid == NULL) {
1253 rtwarn("Cannot fix invalid geometry");
1254 break;
1255 }
1256
1257 lwgeom_free(lwgeom);
1258 lwgeom = lwgeomValid;
1259 }
1260 while (0);
1261
1262 /* save lwgeom */
1263 pols[j].geom = lwgeom_as_lwpoly(lwgeom);
1264
1265#if POSTGIS_DEBUG_LEVEL > 3
1266 {
1267 char *wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, NULL);
1268 RASTER_DEBUGF(4, "LWGEOM wkt = %s", wkt);
1269 rtdealloc(wkt);
1270
1271 size_t lwwkbsize = 0;
1272 uint8_t *lwwkb = lwgeom_to_wkb(lwgeom, WKB_ISO | WKB_NDR, &lwwkbsize);
1273 if (lwwkbsize) {
1274 d_print_binary_hex("LWGEOM wkb", lwwkb, lwwkbsize);
1275 rtdealloc(lwwkb);
1276 }
1277 }
1278#endif
1279
1280 /* set pixel value */
1281 pols[j].val = dValue;
1282 }
1283
1284 *pnElements = nFeatureCount;
1285
1286 RASTER_DEBUG(3, "destroying GDAL MEM raster");
1287 GDALClose(memdataset);
1288 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1289
1290 RASTER_DEBUG(3, "destroying OGR MEM vector");
1291 OGR_Fld_Destroy(hFldDfn);
1292 OGR_DS_DeleteLayer(memdatasource, 0);
1293 if (NULL != pszQuery) rtdealloc(pszQuery);
1294 OGRReleaseDataSource(memdatasource);
1295
1296 return pols;
1297}
TRUE
#define TRUE
Definition dbfopen.c:169
FALSE
#define FALSE
Definition dbfopen.c:168
LWGEOM2GEOS
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
Definition liblwgeom/lwgeom_geos.c:390
lwgeom_geos_error
void lwgeom_geos_error(const char *fmt,...)
Definition liblwgeom/lwgeom_geos.c:43
WKB_ISO
#define WKB_ISO
Definition liblwgeom.h:2121
lwgeom_set_srid
void lwgeom_set_srid(LWGEOM *geom, int32_t srid)
Set the SRID on an LWGEOM For collections, only the parent gets an SRID, all the children get SRID_UN...
Definition lwgeom.c:1530
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition lwgeom.c:1138
LW_PARSER_CHECK_NONE
#define LW_PARSER_CHECK_NONE
Definition liblwgeom.h:2060
lwgeom_to_wkt
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition lwout_wkt.c:676
lwgeom_as_lwpoly
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition lwgeom.c:197
WKT_ISO
#define WKT_ISO
Definition liblwgeom.h:2130
WKB_NDR
#define WKB_NDR
Definition liblwgeom.h:2124
lwgeom_to_wkb
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:790
lwgeom_from_wkb
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:833
lwgeom_make_valid
LWGEOM * lwgeom_make_valid(LWGEOM *geom)
Attempts to make an invalid geometries valid w/out losing points.
Definition liblwgeom/lwgeom_geos_clean.c:854
rterror
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition rt_context.c:199
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition rt_context.c:171
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition librtcore.h:295
rt_raster_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition rt_raster.c:356
rt_util_gdal_register_all
int rt_util_gdal_register_all(int force_register_all)
Definition rt_util.c:338
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition librtcore.h:299
rtinfo
void rtinfo(const char *fmt,...)
Definition rt_context.c:211
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:674
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition rt_band.c:714
rtwarn
void rtwarn(const char *fmt,...)
Definition rt_context.c:224
rt_geomval
struct rt_geomval_t * rt_geomval
Definition librtcore.h:149
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition rt_band.c:1730
rtdealloc
void rtdealloc(void *mem)
Definition rt_context.c:186
rt_raster_to_gdal_mem
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:1821
rt_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition rt_raster.c:381
pixval.nband
nband
Definition pixval.py:53
rt_band_t
Definition librtcore.h:2314
rt_geomval_t::val
double val
Definition librtcore.h:2355
rt_geomval_t::geom
LWPOLY * geom
Definition librtcore.h:2354
rt_geomval_t
Definition librtcore.h:2353

References 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().

Here is the call graph for this function:
Here is the caller graph for this function: