rt_raster_surface()

rt_errorstate rt_raster_surface	(	rt_raster	raster,
		int	nband,
		LWMPOLY **	surface
	)

Get a raster as a surface (multipolygon).

If a band is specified, those pixels with value (not NODATA) contribute to the area of the output multipolygon.

Parameters

raster	: the raster to convert to a multipolygon
nband	: the 0-based band of raster rast to use if value is less than zero, bands are ignored.
*surface	: raster as a surface (multipolygon). if all pixels are NODATA, NULL is set

Returns

ES_NONE on success, ES_ERROR on error

Definition at line 355 of file rt_geometry.c.

                                                                                {
        rt_band band = NULL;
        LWGEOM *mpoly = NULL;
        LWGEOM *tmp = NULL;
        LWGEOM *clone = NULL;
        rt_geomval gv = NULL;
        int gvcount = 0;
        GEOSGeometry *gc = NULL;
        GEOSGeometry *gunion = NULL;
        GEOSGeometry **geoms = NULL;
        int geomscount = 0;
        int i = 0;
 
        assert(surface != NULL);
 
        /* init *surface to NULL */
        *surface = NULL;
 
        /* raster is empty, surface = NULL */
        if (rt_raster_is_empty(raster))
                return ES_NONE;
 
        /* if nband < 0, return the convex hull as a multipolygon */
        if (nband < 0) {
                /*
                        lwgeom_as_multi() only does a shallow clone internally
                        so input and output geometries may share memory
                        hence the deep clone of the output geometry for returning
                        is the only way to guarentee the memory isn't shared
                */
                if (rt_raster_get_convex_hull(raster, &tmp) != ES_NONE) {
                        rterror("rt_raster_surface: Could not get convex hull of raster");
                        return ES_ERROR;
                }
                mpoly = lwgeom_as_multi(tmp);
                clone = lwgeom_clone_deep(mpoly);
                lwgeom_free(tmp);
                lwgeom_free(mpoly);
 
                *surface = lwgeom_as_lwmpoly(clone);
                return ES_NONE;
        }
        /* check that nband is valid */
        else if (nband >= rt_raster_get_num_bands(raster)) {
                rterror("rt_raster_surface: The band index %d is invalid", nband);
                return ES_ERROR;
        }
 
        /* get band */
        band = rt_raster_get_band(raster, nband);
        if (band == NULL) {
                rterror("rt_raster_surface: Error getting band %d from raster", nband);
                return ES_ERROR;
        }
 
        /* band does not have a NODATA flag, return convex hull */
        if (!rt_band_get_hasnodata_flag(band)) {
                /*
                        lwgeom_as_multi() only does a shallow clone internally
                        so input and output geometries may share memory
                        hence the deep clone of the output geometry for returning
                        is the only way to guarentee the memory isn't shared
                */
                if (rt_raster_get_convex_hull(raster, &tmp) != ES_NONE) {
                        rterror("rt_raster_surface: Could not get convex hull of raster");
                        return ES_ERROR;
                }
                mpoly = lwgeom_as_multi(tmp);
                clone = lwgeom_clone_deep(mpoly);
                lwgeom_free(tmp);
                lwgeom_free(mpoly);
 
                *surface = lwgeom_as_lwmpoly(clone);
                return ES_NONE;
        }
        /* band is NODATA, return NULL */
        else if (rt_band_get_isnodata_flag(band)) {
                RASTER_DEBUG(3, "Band is NODATA.  Returning NULL");
                return ES_NONE;
        }
 
        /* use gdal polygonize */
        gv = rt_raster_gdal_polygonize(raster, nband, 1, &gvcount);
        /* no polygons returned */
        if (gvcount < 1) {
                RASTER_DEBUG(3, "All pixels of band are NODATA.  Returning NULL");
                if (gv != NULL) rtdealloc(gv);
                return ES_NONE;
        }
        /* more than 1 polygon */
        else if (gvcount > 1) {
                /* convert LWPOLY to GEOSGeometry */
                geomscount = gvcount;
                geoms = rtalloc(sizeof(GEOSGeometry *) * geomscount);
                if (geoms == NULL) {
                        rterror("rt_raster_surface: Could not allocate memory for pixel polygons as GEOSGeometry");
                        for (i = 0; i < gvcount; i++) lwpoly_free(gv[i].geom);
                        rtdealloc(gv);
                        return ES_ERROR;
                }
                for (i = 0; i < gvcount; i++) {
#if POSTGIS_DEBUG_LEVEL > 3
                        {
                                char *wkt = lwgeom_to_wkt(lwpoly_as_lwgeom(gv[i].geom), WKT_ISO, DBL_DIG, NULL);
                                RASTER_DEBUGF(4, "geom %d = %s", i, wkt);
                                rtdealloc(wkt);
                        }
#endif
 
                        geoms[i] = LWGEOM2GEOS(lwpoly_as_lwgeom(gv[i].geom), 0);
                        lwpoly_free(gv[i].geom);
                }
                rtdealloc(gv);
 
                /* create geometry collection */
                gc = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, geomscount);
 
                if (gc == NULL) {
                        rterror("rt_raster_surface: Could not create GEOS GEOMETRYCOLLECTION from set of pixel polygons");
 
                        for (i = 0; i < geomscount; i++)
                                GEOSGeom_destroy(geoms[i]);
                        rtdealloc(geoms);
                        return ES_ERROR;
                }
 
                /* run the union */
                gunion = GEOSUnaryUnion(gc);
 
                GEOSGeom_destroy(gc);
                rtdealloc(geoms);
 
                if (gunion == NULL) {
                        rterror("rt_raster_surface: Could not union the pixel polygons using GEOSUnaryUnion()");
                        return ES_ERROR;
                }
 
                /* convert union result from GEOSGeometry to LWGEOM */
                mpoly = GEOS2LWGEOM(gunion, 0);
 
                /*
                        is geometry valid?
                        if not, try to make valid
                */
                do {
                        LWGEOM *mpolyValid = NULL;
 
                        if (GEOSisValid(gunion))
                                break;
 
                        /* make geometry valid */
                        mpolyValid = lwgeom_make_valid(mpoly);
                        if (mpolyValid == NULL) {
                                rtwarn("Cannot fix invalid geometry");
                                break;
                        }
 
                        lwgeom_free(mpoly);
                        mpoly = mpolyValid;
                }
                while (0);
 
                GEOSGeom_destroy(gunion);
        }
        else {
                mpoly = lwpoly_as_lwgeom(gv[0].geom);
                rtdealloc(gv);
 
#if POSTGIS_DEBUG_LEVEL > 3
                        {
                                char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
                                RASTER_DEBUGF(4, "geom 0 = %s", wkt);
                                rtdealloc(wkt);
                        }
#endif
        }
 
        /* specify SRID */
        lwgeom_set_srid(mpoly, rt_raster_get_srid(raster));
 
        if (mpoly != NULL) {
                /* convert to multi */
                if (!lwgeom_is_collection(mpoly)) {
                        tmp = mpoly;
 
#if POSTGIS_DEBUG_LEVEL > 3
                        {
                                char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
                                RASTER_DEBUGF(4, "before multi = %s", wkt);
                                rtdealloc(wkt);
                        }
#endif
 
                        RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
 
                        /*
                                lwgeom_as_multi() only does a shallow clone internally
                                so input and output geometries may share memory
                                hence the deep clone of the output geometry for returning
                                is the only way to guarentee the memory isn't shared
                        */
                        mpoly = lwgeom_as_multi(tmp);
                        clone = lwgeom_clone_deep(mpoly);
                        lwgeom_free(tmp);
                        lwgeom_free(mpoly);
                        mpoly = clone;
 
                        RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
 
#if POSTGIS_DEBUG_LEVEL > 3
                        {
                                char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
                                RASTER_DEBUGF(4, "after multi = %s", wkt);
                                rtdealloc(wkt);
                        }
#endif
                }
 
#if POSTGIS_DEBUG_LEVEL > 3
                {
                        char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
                        RASTER_DEBUGF(4, "returning geometry = %s", wkt);
                        rtdealloc(wkt);
                }
#endif
 
                *surface = lwgeom_as_lwmpoly(mpoly);
                return ES_NONE;
        }
 
        return ES_NONE;
}

References ovdump::band, ES_ERROR, ES_NONE, GEOS2LWGEOM(), LWGEOM2GEOS(), lwgeom_as_lwmpoly(), lwgeom_as_multi(), lwgeom_clone_deep(), lwgeom_free(), lwgeom_is_collection(), lwgeom_make_valid(), lwgeom_set_srid(), lwgeom_to_wkt(), lwpoly_as_lwgeom(), lwpoly_free(), pixval::nband, rtrowdump::raster, RASTER_DEBUG, RASTER_DEBUGF, rt_band_get_hasnodata_flag(), rt_band_get_isnodata_flag(), rt_raster_gdal_polygonize(), rt_raster_get_band(), rt_raster_get_convex_hull(), rt_raster_get_num_bands(), rt_raster_get_srid(), rt_raster_is_empty(), rtalloc(), rtdealloc(), rterror(), rtwarn(), and WKT_ISO.

Referenced by RASTER_getPolygon(), rt_raster_fully_within_distance(), rt_raster_geos_spatial_relationship(), rt_raster_within_distance(), and test_raster_surface().

Here is the call graph for this function:

Here is the caller graph for this function: