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

  {
        CPLErr cplerr = CE_None;
        char *pszQuery;
        long j;
        OGRSFDriverH ogr_drv = NULL;
        GDALDriverH gdal_drv = NULL;
        int destroy_gdal_drv = 0;
        GDALDatasetH memdataset = NULL;
        GDALRasterBandH gdal_band = NULL;
        OGRDataSourceH memdatasource = NULL;
        rt_geomval pols = NULL;
        OGRLayerH hLayer = NULL;
        OGRFeatureH hFeature = NULL;
        OGRGeometryH hGeom = NULL;
        OGRFieldDefnH hFldDfn = NULL;
        unsigned char *wkb = NULL;
        int wkbsize = 0;
        LWGEOM *lwgeom = NULL;
        int nFeatureCount = 0;
        rt_band band = NULL;
        int iPixVal = -1;
        double dValue = 0.0;
        int iBandHasNodataValue = FALSE;
        double dBandNoData = 0.0;

        /* for checking that a geometry is valid */
        GEOSGeometry *ggeom = NULL;
        int isValid;
        LWGEOM *lwgeomValid = NULL;

        uint32_t bandNums[1] = {nband};
        int excludeNodataValues[1] = {exclude_nodata_value};

        /* checks */
        assert(NULL != raster);
        assert(NULL != pnElements);

        RASTER_DEBUG(2, "In rt_raster_gdal_polygonize");

        *pnElements = 0;

        /*******************************
         * Get band
         *******************************/
        band = rt_raster_get_band(raster, nband);
        if (NULL == band) {
                rterror("rt_raster_gdal_polygonize: Error getting band %d from raster", nband);
                return NULL;
        }

        if (exclude_nodata_value) {

                /* band is NODATA */
                if (rt_band_get_isnodata_flag(band)) {
                        RASTER_DEBUG(3, "Band is NODATA.  Returning null");
                        *pnElements = 0;
                        return NULL;
                }

                iBandHasNodataValue = rt_band_get_hasnodata_flag(band);
                if (iBandHasNodataValue)
                        rt_band_get_nodata(band, &dBandNoData);
                else
                        exclude_nodata_value = FALSE;
        }

        /*****************************************************
         * Convert raster to GDAL MEM dataset
         *****************************************************/
        memdataset = rt_raster_to_gdal_mem(raster, NULL, bandNums, excludeNodataValues, 1, &gdal_drv, &destroy_gdal_drv);
        if (NULL == memdataset) {
                rterror("rt_raster_gdal_polygonize: Couldn't convert raster to GDAL MEM dataset");
                return NULL;
        }

        /*****************************
         * Register ogr mem driver
         *****************************/
#ifdef GDAL_DCAP_RASTER
        /* in GDAL 2.0, OGRRegisterAll() is an alias to GDALAllRegister() */
        rt_util_gdal_register_all(0);
#else
        OGRRegisterAll();
#endif

        RASTER_DEBUG(3, "creating OGR MEM vector");

        /*****************************************************
         * Create an OGR in-memory vector for layers
         *****************************************************/
        ogr_drv = OGRGetDriverByName("Memory");
        memdatasource = OGR_Dr_CreateDataSource(ogr_drv, "", NULL);
        if (NULL == memdatasource) {
                rterror("rt_raster_gdal_polygonize: Couldn't create a OGR Datasource to store pols");
                GDALClose(memdataset);
                if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
                return NULL;
        }

        /* Can MEM driver create new layers? */
        if (!OGR_DS_TestCapability(memdatasource, ODsCCreateLayer)) {
                rterror("rt_raster_gdal_polygonize: MEM driver can't create new layers, aborting");

                /* xxx jorgearevalo: what should we do now? */
                GDALClose(memdataset);
                if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
                OGRReleaseDataSource(memdatasource);

                return NULL;
        }

        RASTER_DEBUG(3, "polygonizying GDAL MEM raster band");

        /*****************************
         * Polygonize the raster band
         *****************************/

        hLayer = OGR_DS_CreateLayer(memdatasource, "PolygonizedLayer", NULL, wkbPolygon, NULL);

        if (NULL == hLayer) {
                rterror("rt_raster_gdal_polygonize: Couldn't create layer to store polygons");

                GDALClose(memdataset);
                if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
                OGRReleaseDataSource(memdatasource);

                return NULL;
        }

        /* First, create a field definition to create the field */
        hFldDfn = OGR_Fld_Create("PixelValue", OFTReal);

        /* Second, create the field */
        if (OGR_L_CreateField(hLayer, hFldDfn, TRUE) != OGRERR_NONE) {
                rtwarn("Couldn't create a field in OGR Layer. The polygons generated won't be able to store the pixel value");
                iPixVal = -1;
        }
        else {
                /* Index to the new field created in the layer */
                iPixVal = 0;
        }

        /* Get GDAL raster band */
        gdal_band = GDALGetRasterBand(memdataset, 1);
        if (NULL == gdal_band) {
                rterror("rt_raster_gdal_polygonize: Couldn't get GDAL band to polygonize");

                GDALClose(memdataset);
                if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
                OGR_Fld_Destroy(hFldDfn);
                OGR_DS_DeleteLayer(memdatasource, 0);
                OGRReleaseDataSource(memdatasource);

                return NULL;
        }

        /* We don't need a raster mask band. Each band has a nodata value. */
        cplerr = GDALFPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);

        if (cplerr != CE_None) {
                rterror("rt_raster_gdal_polygonize: Could not polygonize GDAL band");

                GDALClose(memdataset);
                if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
                OGR_Fld_Destroy(hFldDfn);
                OGR_DS_DeleteLayer(memdatasource, 0);
                OGRReleaseDataSource(memdatasource);

                return NULL;
        }

        if (iBandHasNodataValue) {
                pszQuery = (char *) rtalloc(50 * sizeof (char));
                sprintf(pszQuery, "PixelValue != %f", dBandNoData );
                OGRErr e = OGR_L_SetAttributeFilter(hLayer, pszQuery);
                if (e != OGRERR_NONE) {
                        rtwarn("Error filtering NODATA values for band. All values will be treated as data values");
                }
        }
        else {
                pszQuery = NULL;
        }

        /*********************************************************************
         * Transform OGR layers to WKB polygons
         * XXX jorgearevalo: GDALPolygonize does not set the coordinate system
         * on the output layer. Application code should do this when the layer
         * is created, presumably matching the raster coordinate system.
         *********************************************************************/
        nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE);

        /* Allocate memory for pols */
        pols = (rt_geomval) rtalloc(nFeatureCount * sizeof(struct rt_geomval_t));

        if (NULL == pols) {
                rterror("rt_raster_gdal_polygonize: Could not allocate memory for geomval set");

                GDALClose(memdataset);
                if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
                OGR_Fld_Destroy(hFldDfn);
                OGR_DS_DeleteLayer(memdatasource, 0);
                if (NULL != pszQuery)
                        rtdealloc(pszQuery);
                OGRReleaseDataSource(memdatasource);

                return NULL;
        }

        /* initialize GEOS */
        initGEOS(rtinfo, lwgeom_geos_error);

        RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount);

        /* Reset feature reading to start in the first feature */
        OGR_L_ResetReading(hLayer);

        for (j = 0; j < nFeatureCount; j++) {
                hFeature = OGR_L_GetNextFeature(hLayer);
                dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal);

                hGeom = OGR_F_GetGeometryRef(hFeature);
                wkbsize = OGR_G_WkbSize(hGeom);

                /* allocate wkb buffer */
                wkb = rtalloc(sizeof(unsigned char) * wkbsize);
                if (wkb == NULL) {
                        rterror("rt_raster_gdal_polygonize: Could not allocate memory for WKB buffer");

                        OGR_F_Destroy(hFeature);
                        GDALClose(memdataset);
                        if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
                        OGR_Fld_Destroy(hFldDfn);
                        OGR_DS_DeleteLayer(memdatasource, 0);
                        if (NULL != pszQuery)
                                rtdealloc(pszQuery);
                        OGRReleaseDataSource(memdatasource);

                        return NULL;
                }

                /* export WKB with LSB byte order */
                OGR_G_ExportToWkb(hGeom, wkbNDR, wkb);

                /* convert WKB to LWGEOM */
                lwgeom = lwgeom_from_wkb(wkb, wkbsize, LW_PARSER_CHECK_NONE);

#if POSTGIS_DEBUG_LEVEL > 3
                {
                        char *wkt = NULL;
                        OGR_G_ExportToWkt(hGeom, &wkt);
                        RASTER_DEBUGF(4, "GDAL wkt = %s", wkt);
                        CPLFree(wkt);

                        d_print_binary_hex("GDAL wkb", wkb, wkbsize);
                }
#endif

                /* cleanup unnecessary stuff */
                rtdealloc(wkb);
                wkb = NULL;
                wkbsize = 0;

                OGR_F_Destroy(hFeature);

                /* specify SRID */
                lwgeom_set_srid(lwgeom, rt_raster_get_srid(raster));

                /*
                        is geometry valid?
                        if not, try to make valid
                */
                do {
                        ggeom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom, 0);
                        if (ggeom == NULL) {
                                rtwarn("Cannot test geometry for validity");
                                break;
                        }

                        isValid = GEOSisValid(ggeom);

                        GEOSGeom_destroy(ggeom);
                        ggeom = NULL;

                        /* geometry is valid */
                        if (isValid)
                                break;

                        RASTER_DEBUG(3, "fixing invalid geometry");

                        /* make geometry valid */
                        lwgeomValid = lwgeom_make_valid(lwgeom);
                        if (lwgeomValid == NULL) {
                                rtwarn("Cannot fix invalid geometry");
                                break;
                        }

                        lwgeom_free(lwgeom);
                        lwgeom = lwgeomValid;
                }
                while (0);

                /* save lwgeom */
                pols[j].geom = lwgeom_as_lwpoly(lwgeom);

#if POSTGIS_DEBUG_LEVEL > 3
                {
                        char *wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, NULL);
                        RASTER_DEBUGF(4, "LWGEOM wkt = %s", wkt);
                        rtdealloc(wkt);

                        size_t lwwkbsize = 0;
                        uint8_t *lwwkb = lwgeom_to_wkb(lwgeom, WKB_ISO | WKB_NDR, &lwwkbsize);
                        if (lwwkbsize) {
                                d_print_binary_hex("LWGEOM wkb", lwwkb, lwwkbsize);
                                rtdealloc(lwwkb);
                        }
                }
#endif

                /* set pixel value */
                pols[j].val = dValue;
        }

        *pnElements = nFeatureCount;

        RASTER_DEBUG(3, "destroying GDAL MEM raster");
        GDALClose(memdataset);
        if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);

        RASTER_DEBUG(3, "destroying OGR MEM vector");
        OGR_Fld_Destroy(hFldDfn);
        OGR_DS_DeleteLayer(memdatasource, 0);
        if (NULL != pszQuery) rtdealloc(pszQuery);
        OGRReleaseDataSource(memdatasource);

        return pols;
}

Here is the call graph for this function:

Here is the caller graph for this function: