PostGIS 3.6.2dev-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 975 of file rt_geometry.c.

979 {
980 CPLErr cplerr = CE_None;
981 char *pszQuery;
982 long j;
983 OGRSFDriverH ogr_drv = NULL;
984 GDALDriverH gdal_drv = NULL;
985 int destroy_gdal_drv = 0;
986 GDALDatasetH memdataset = NULL;
987 GDALRasterBandH gdal_band = NULL;
988 OGRDataSourceH memdatasource = NULL;
989 rt_geomval pols = NULL;
990 OGRLayerH hLayer = NULL;
991 OGRFeatureH hFeature = NULL;
992 OGRGeometryH hGeom = NULL;
993 OGRFieldDefnH hFldDfn = NULL;
994 unsigned char *wkb = NULL;
995 int wkbsize = 0;
996 LWGEOM *lwgeom = NULL;
997 int nFeatureCount = 0;
998 rt_band band = NULL;
999 int iPixVal = -1;
1000 double dValue = 0.0;
1001 int iBandHasNodataValue = FALSE;
1002 double dBandNoData = 0.0;
1003
1004 uint32_t bandNums[1] = {nband};
1005 int excludeNodataValues[1] = {exclude_nodata_value};
1006
1007 /* checks */
1008 assert(NULL != raster);
1009 assert(NULL != pnElements);
1010
1011 RASTER_DEBUG(2, "In rt_raster_gdal_polygonize");
1012
1013 *pnElements = 0;
1014
1015 /*******************************
1016 * Get band
1017 *******************************/
1018 band = rt_raster_get_band(raster, nband);
1019 if (NULL == band) {
1020 rterror("rt_raster_gdal_polygonize: Error getting band %d from raster", nband);
1021 return NULL;
1022 }
1023
1024 if (exclude_nodata_value) {
1025
1026 /* band is NODATA */
1027 if (rt_band_get_isnodata_flag(band)) {
1028 RASTER_DEBUG(3, "Band is NODATA. Returning null");
1029 *pnElements = 0;
1030 return NULL;
1031 }
1032
1033 iBandHasNodataValue = rt_band_get_hasnodata_flag(band);
1034 if (iBandHasNodataValue)
1035 rt_band_get_nodata(band, &dBandNoData);
1036 else
1037 exclude_nodata_value = FALSE;
1038 }
1039
1040 /*****************************************************
1041 * Convert raster to GDAL MEM dataset
1042 *****************************************************/
1043 memdataset = rt_raster_to_gdal_mem(raster, NULL, bandNums, excludeNodataValues, 1, &gdal_drv, &destroy_gdal_drv);
1044 if (NULL == memdataset) {
1045 rterror("rt_raster_gdal_polygonize: Couldn't convert raster to GDAL MEM dataset");
1046 return NULL;
1047 }
1048
1049 /*****************************
1050 * Register ogr mem driver
1051 *****************************/
1052#ifdef GDAL_DCAP_RASTER
1053 /* in GDAL 2.0, OGRRegisterAll() is an alias to GDALAllRegister() */
1054 rt_util_gdal_register_all(0);
1055#else
1056 OGRRegisterAll();
1057#endif
1058
1059 RASTER_DEBUG(3, "creating OGR MEM vector");
1060
1061 /*****************************************************
1062 * Create an OGR in-memory vector for layers
1063 *****************************************************/
1064 ogr_drv = OGRGetDriverByName("Memory");
1065 memdatasource = OGR_Dr_CreateDataSource(ogr_drv, "", NULL);
1066 if (NULL == memdatasource) {
1067 rterror("rt_raster_gdal_polygonize: Couldn't create a OGR Datasource to store pols");
1068 GDALClose(memdataset);
1069 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1070 return NULL;
1071 }
1072
1073 /* Can MEM driver create new layers? */
1074 if (!OGR_DS_TestCapability(memdatasource, ODsCCreateLayer)) {
1075 rterror("rt_raster_gdal_polygonize: MEM driver can't create new layers, aborting");
1076
1077 /* xxx jorgearevalo: what should we do now? */
1078 GDALClose(memdataset);
1079 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1080 OGRReleaseDataSource(memdatasource);
1081
1082 return NULL;
1083 }
1084
1085 RASTER_DEBUG(3, "polygonizying GDAL MEM raster band");
1086
1087 /*****************************
1088 * Polygonize the raster band
1089 *****************************/
1090
1096 hLayer = OGR_DS_CreateLayer(memdatasource, "PolygonizedLayer", NULL, wkbPolygon, NULL);
1097
1098 if (NULL == hLayer) {
1099 rterror("rt_raster_gdal_polygonize: Couldn't create layer to store polygons");
1100
1101 GDALClose(memdataset);
1102 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1103 OGRReleaseDataSource(memdatasource);
1104
1105 return NULL;
1106 }
1107
1112 /* First, create a field definition to create the field */
1113 hFldDfn = OGR_Fld_Create("PixelValue", OFTReal);
1114
1115 /* Second, create the field */
1116 if (OGR_L_CreateField(hLayer, hFldDfn, TRUE) != OGRERR_NONE) {
1117 rtwarn("Couldn't create a field in OGR Layer. The polygons generated won't be able to store the pixel value");
1118 iPixVal = -1;
1119 }
1120 else {
1121 /* Index to the new field created in the layer */
1122 iPixVal = 0;
1123 }
1124
1125 /* Get GDAL raster band */
1126 gdal_band = GDALGetRasterBand(memdataset, 1);
1127 if (NULL == gdal_band) {
1128 rterror("rt_raster_gdal_polygonize: Couldn't get GDAL band to polygonize");
1129
1130 GDALClose(memdataset);
1131 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1132 OGR_Fld_Destroy(hFldDfn);
1133 OGR_DS_DeleteLayer(memdatasource, 0);
1134 OGRReleaseDataSource(memdatasource);
1135
1136 return NULL;
1137 }
1138
1139 /* We don't need a raster mask band. Each band has a nodata value. */
1140 cplerr = GDALFPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
1141
1142 if (cplerr != CE_None) {
1143 rterror("rt_raster_gdal_polygonize: Could not polygonize GDAL band");
1144
1145 GDALClose(memdataset);
1146 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1147 OGR_Fld_Destroy(hFldDfn);
1148 OGR_DS_DeleteLayer(memdatasource, 0);
1149 OGRReleaseDataSource(memdatasource);
1150
1151 return NULL;
1152 }
1153
1160 if (iBandHasNodataValue) {
1161 size_t sz = 50 * sizeof (char);
1162 pszQuery = (char *) rtalloc(sz);
1163 snprintf(pszQuery, sz, "PixelValue != %f", dBandNoData );
1164 OGRErr e = OGR_L_SetAttributeFilter(hLayer, pszQuery);
1165 if (e != OGRERR_NONE) {
1166 rtwarn("Error filtering NODATA values for band. All values will be treated as data values");
1167 }
1168 }
1169 else {
1170 pszQuery = NULL;
1171 }
1172
1173 /*********************************************************************
1174 * Transform OGR layers to WKB polygons
1175 * XXX jorgearevalo: GDALPolygonize does not set the coordinate system
1176 * on the output layer. Application code should do this when the layer
1177 * is created, presumably matching the raster coordinate system.
1178 *********************************************************************/
1179 nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE);
1180
1181 /* Allocate memory for pols */
1182 pols = (rt_geomval) rtalloc(nFeatureCount * sizeof(struct rt_geomval_t));
1183
1184 if (NULL == pols) {
1185 rterror("rt_raster_gdal_polygonize: Could not allocate memory for geomval set");
1186
1187 GDALClose(memdataset);
1188 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1189 OGR_Fld_Destroy(hFldDfn);
1190 OGR_DS_DeleteLayer(memdatasource, 0);
1191 if (NULL != pszQuery)
1192 rtdealloc(pszQuery);
1193 OGRReleaseDataSource(memdatasource);
1194
1195 return NULL;
1196 }
1197
1198 /* initialize GEOS */
1199 initGEOS(rtinfo, lwgeom_geos_error);
1200
1201 RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount);
1202
1203 /* Reset feature reading to start in the first feature */
1204 OGR_L_ResetReading(hLayer);
1205
1206 for (j = 0; j < nFeatureCount; j++) {
1207 hFeature = OGR_L_GetNextFeature(hLayer);
1208 dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal);
1209
1210 hGeom = OGR_F_GetGeometryRef(hFeature);
1211 wkbsize = OGR_G_WkbSize(hGeom);
1212
1213 /* allocate wkb buffer */
1214 wkb = rtalloc(sizeof(unsigned char) * wkbsize);
1215 if (wkb == NULL) {
1216 rterror("rt_raster_gdal_polygonize: Could not allocate memory for WKB buffer");
1217
1218 OGR_F_Destroy(hFeature);
1219 GDALClose(memdataset);
1220 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1221 OGR_Fld_Destroy(hFldDfn);
1222 OGR_DS_DeleteLayer(memdatasource, 0);
1223 if (NULL != pszQuery)
1224 rtdealloc(pszQuery);
1225 OGRReleaseDataSource(memdatasource);
1226
1227 return NULL;
1228 }
1229
1230 /* export WKB with LSB byte order */
1231 OGR_G_ExportToWkb(hGeom, wkbNDR, wkb);
1232
1233 /* convert WKB to LWGEOM */
1234 lwgeom = lwgeom_from_wkb(wkb, wkbsize, LW_PARSER_CHECK_NONE);
1235 if (!lwgeom) rterror("%s: invalid wkb", __func__);
1236
1237#if POSTGIS_DEBUG_LEVEL > 3
1238 {
1239 char *wkt = NULL;
1240 OGR_G_ExportToWkt(hGeom, &wkt);
1241 RASTER_DEBUGF(4, "GDAL wkt = %s", wkt);
1242 CPLFree(wkt);
1243
1244 d_print_binary_hex("GDAL wkb", wkb, wkbsize);
1245 }
1246#endif
1247
1248 /* cleanup unnecessary stuff */
1249 rtdealloc(wkb);
1250 wkb = NULL;
1251 wkbsize = 0;
1252
1253 OGR_F_Destroy(hFeature);
1254
1255 /* specify SRID */
1256 lwgeom_set_srid(lwgeom, rt_raster_get_srid(raster));
1257
1258 /* save lwgeom */
1259 pols[j].geom = lwgeom_as_lwpoly(lwgeom);
1260
1261#if POSTGIS_DEBUG_LEVEL > 3
1262 {
1263 char *wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, NULL);
1264 RASTER_DEBUGF(4, "LWGEOM wkt = %s", wkt);
1265 rtdealloc(wkt);
1266
1267 lwvarlena_t *lwwkb = lwgeom_to_wkb_varlena(lwgeom, WKB_ISO | WKB_NDR);
1268 size_t lwwkbsize = LWSIZE_GET(lwwkb->size) - LWVARHDRSZ;
1269 if (lwwkbsize) {
1270 d_print_binary_hex("LWGEOM wkb", (const uint8_t *)lwwkb->data, lwwkbsize);
1271 rtdealloc(lwwkb);
1272 }
1273 }
1274#endif
1275
1276 /* set pixel value */
1277 pols[j].val = dValue;
1278 }
1279
1280 *pnElements = nFeatureCount;
1281
1282 RASTER_DEBUG(3, "destroying GDAL MEM raster");
1283 GDALClose(memdataset);
1284 if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1285
1286 RASTER_DEBUG(3, "destroying OGR MEM vector");
1287 OGR_Fld_Destroy(hFldDfn);
1288 OGR_DS_DeleteLayer(memdatasource, 0);
1289 if (NULL != pszQuery) rtdealloc(pszQuery);
1290 OGRReleaseDataSource(memdatasource);
1291
1292 return pols;
1293}
TRUE
#define TRUE
Definition dbfopen.c:73
FALSE
#define FALSE
Definition dbfopen.c:72
lwgeom_geos_error
void lwgeom_geos_error(const char *fmt,...)
Definition liblwgeom/lwgeom_geos.c:57
WKB_ISO
#define WKB_ISO
Definition liblwgeom.h:2207
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:1610
LWVARHDRSZ
#define LWVARHDRSZ
Definition liblwgeom.h:311
LW_PARSER_CHECK_NONE
#define LW_PARSER_CHECK_NONE
Definition liblwgeom.h:2146
LWSIZE_GET
#define LWSIZE_GET(varsize)
Macro for reading the size from the GSERIALIZED size attribute.
Definition liblwgeom.h:324
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:708
lwgeom_to_wkb_varlena
lwvarlena_t * lwgeom_to_wkb_varlena(const LWGEOM *geom, uint8_t variant)
Definition lwout_wkb.c:851
lwgeom_as_lwpoly
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition lwgeom.c:215
WKT_ISO
#define WKT_ISO
Definition liblwgeom.h:2216
WKB_NDR
#define WKB_NDR
Definition liblwgeom.h:2210
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:842
rterror
void rterror(const char *fmt,...) __attribute__((format(printf
Wrappers used for reporting errors and info.
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition rt_context.c:191
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition librtcore.h:302
rt_raster_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition rt_raster.c:360
rt_util_gdal_register_all
int rt_util_gdal_register_all(int force_register_all)
Definition rt_util.c:347
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition librtcore.h:306
rtinfo
void void rtinfo(const char *fmt,...) __attribute__((format(printf
rtwarn
void void void rtwarn(const char *fmt,...) __attribute__((format(printf
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:825
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition rt_band.c:865
rt_geomval
struct rt_geomval_t * rt_geomval
Definition librtcore.h:150
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition rt_band.c:2038
rtdealloc
void rtdealloc(void *mem)
Definition rt_context.c:206
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:1813
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:385
pixval.nband
nband
Definition pixval.py:54
lwvarlena_t::size
uint32_t size
Definition liblwgeom.h:307
lwvarlena_t::data
char data[]
Definition liblwgeom.h:308
lwvarlena_t
Definition liblwgeom.h:306
rt_band_t
Definition librtcore.h:2503
rt_geomval_t::val
double val
Definition librtcore.h:2544
rt_geomval_t::geom
LWPOLY * geom
Definition librtcore.h:2543
rt_geomval_t
Definition librtcore.h:2542

References lwvarlena_t::data, FALSE, rt_geomval_t::geom, LW_PARSER_CHECK_NONE, lwgeom_as_lwpoly(), lwgeom_from_wkb(), lwgeom_geos_error(), lwgeom_set_srid(), lwgeom_to_wkb_varlena(), lwgeom_to_wkt(), LWSIZE_GET, LWVARHDRSZ, 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(), lwvarlena_t::size, 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: