From: Bborie Park Date: Fri, 13 Jul 2012 23:28:32 +0000 (+0000) Subject: Initial commit of rt_raster_surface() with base regression tests X-Git-Tag: 2.1.0beta2~806 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=fb817af897f1d3c14544a1d8065f1c9efa12e5a0;p=postgis Initial commit of rt_raster_surface() with base regression tests git-svn-id: http://svn.osgeo.org/postgis/trunk@10054 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index 6603622fd..9a4b001f1 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -2585,7 +2585,6 @@ rt_band_get_min_value(rt_band band) { return rt_pixtype_get_min_value(band->pixtype); } - int rt_band_check_is_nodata(rt_band band) { int i, j, err; @@ -11879,7 +11878,20 @@ rt_raster_from_two_rasters( return raster; } - +/** + * Get a raster pixel as a polygon. + * + * The pixel shape is a 4 vertices (5 to be closed) single + * ring polygon bearing the raster's rotation + * and using projection coordinates + * + * @param raster : the raster to get pixel from + * @param x : the column number + * @param y : the row number + * + * @return the pixel polygon, or NULL on error. + * + */ LWPOLY* rt_raster_pixel_as_polygon(rt_raster rast, int x, int y) { @@ -11891,6 +11903,8 @@ rt_raster_pixel_as_polygon(rt_raster rast, int x, int y) POINT4D p, p0; LWPOLY *poly; + assert(rast != NULL); + scale_x = rt_raster_get_x_scale(rast); scale_y = rt_raster_get_y_scale(rast); skew_x = rt_raster_get_x_skew(rast); @@ -11925,3 +11939,177 @@ rt_raster_pixel_as_polygon(rt_raster rast, int x, int y) return poly; } + +/** + * 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. + * + * @param raster: the raster to convert to a multipolygon + * @param nband : the 0-based band of raster rast to use + * if value is less than zero, bands are ignored. + * + * @return the raster surface or NULL on error + */ +LWMPOLY* rt_raster_surface(rt_raster raster, int nband) { + rt_band band = NULL; + LWGEOM *mpoly = NULL; + LWGEOM *tmp = NULL; + LWGEOM *poly = NULL; + uint16_t width = 0; + uint16_t height = 0; + int x = 0; + int y = 0; + int err = 0; + double val = 0; + double nodata = 0; + + /* raster is empty, return NULL */ + if (rt_raster_is_empty(raster)) + return NULL; + + /* if nband < 0, return the convex hull as a multipolygon */ + if (nband < 0) { + return lwgeom_as_lwmpoly( + lwgeom_as_multi( + lwpoly_as_lwgeom( + rt_raster_get_convex_hull(raster) + ) + ) + ); + } + /* 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 NULL; + } + + /* get band */ + band = rt_raster_get_band(raster, nband); + if (band == NULL) { + rterror("rt_raster_surface: Error getting band %d from raster", nband); + return NULL; + } + + /* band is NODATA or does not have a NODATA flag, return convex hull */ + if ( + rt_band_get_isnodata_flag(band) || + !rt_band_get_hasnodata_flag(band) + ) { + return lwgeom_as_lwmpoly( + lwgeom_as_multi( + lwpoly_as_lwgeom( + rt_raster_get_convex_hull(raster) + ) + ) + ); + } + + /* NODATA value */ + nodata = rt_band_get_nodata(band); + + /* process each pixel of band */ + width = rt_raster_get_width(raster); + height = rt_raster_get_height(raster); + for (y = 0; y < height; y++) { + for (x = 0; x < width; x++) { + err = rt_band_get_pixel(band, x, y, &val); + if (err != 0) { + rterror("rt_raster_surface: Unable to get pixel value"); + return NULL; + } + + /* val is NODATA, skip */ + if ( + FLT_EQ(val, nodata) || + rt_band_clamped_value_is_nodata(band, val) != 0 + ) { + continue; + } + + /* convert pixel to polygon */ + poly = lwpoly_as_lwgeom(rt_raster_pixel_as_polygon(raster, x, y)); + if (poly == NULL) { + rterror("rt_raster_surface: Unable to create polygon of pixel"); + if (mpoly != NULL) lwgeom_free(mpoly); + return NULL; + } + + /* union polygon */ + if (mpoly == NULL) + mpoly = poly; + else { + tmp = mpoly; + mpoly = lwgeom_union(tmp, poly); + +#if POSTGIS_DEBUG_LEVEL > 3 + { + char *wkt = NULL; + + wkt = lwgeom_to_wkt(poly, WKT_ISO, DBL_DIG, NULL); + RASTER_DEBUGF(4, "poly = %s", wkt); + rtdealloc(wkt); + + wkt = lwgeom_to_wkt(tmp, WKT_ISO, DBL_DIG, NULL); + RASTER_DEBUGF(4, "tmp = %s", wkt); + rtdealloc(wkt); + + wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL); + RASTER_DEBUGF(4, "mpoly = %s", wkt); + rtdealloc(wkt); + } +#endif + + lwgeom_free(tmp); + lwgeom_free(poly); + + if (mpoly == NULL) { + rterror("rt_raster_surface: Unable to union pixel polygons"); + return NULL; + } + } + } + } + + if (mpoly != NULL) { + /* convert to multi */ + if (!lwgeom_is_collection(mpoly)) { + tmp = mpoly; + +#if POSTGIS_DEBUG_LEVEL > 3 + { + char *wkt = NULL; + wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL); + RASTER_DEBUGF(4, "before multi = %s", wkt); + rtdealloc(wkt); + } +#endif + + mpoly = lwgeom_clone_deep(lwgeom_as_multi(tmp)); + lwgeom_free(tmp); + +#if POSTGIS_DEBUG_LEVEL > 3 + { + char *wkt = NULL; + 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 = NULL; + wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL); + RASTER_DEBUGF(4, "returning geometry = %s", wkt); + rtdealloc(wkt); + } +#endif + + return lwgeom_as_lwmpoly(mpoly); + } + + return NULL; +} diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index 450843983..e4740d1cc 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -605,13 +605,6 @@ int rt_band_get_pixel_of_value( */ double rt_band_get_min_value(rt_band band); -/** - * Returns TRUE if the band is only nodata values - * @param band: the band to get info from - * @return TRUE if the band is only nodata values, FALSE otherwise - */ -int rt_band_is_nodata(rt_band band); - /** * Returns TRUE if the band is only nodata values * @param band: the band to get info from @@ -1178,6 +1171,19 @@ rt_raster_compute_skewed_raster( */ LWPOLY* rt_raster_pixel_as_polygon(rt_raster raster, int x, int y); +/** + * 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. + * + * @param raster: the raster to convert to a multipolygon + * @param nband : the 0-based band of raster rast to use + * if value is less than zero, bands are ignored. + * + * @return the raster surface or NULL on error + */ +LWMPOLY* rt_raster_surface(rt_raster raster, int nband); + /** * Returns a set of "geomval" value, one for each group of pixel * sharing the same value for the provided band. diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c index 83a831776..682ede3d6 100644 --- a/raster/test/core/testapi.c +++ b/raster/test/core/testapi.c @@ -2894,6 +2894,63 @@ static void testPixelOfValue() { deepRelease(rast); } +static void testRasterSurface() { + rt_raster rast; + rt_band band; + const int maxX = 5; + const int maxY = 5; + int x, y; + char *wkt = NULL; + LWMPOLY *mpoly = NULL; + + rast = rt_raster_new(maxX, maxY); + assert(rast); + + rt_raster_set_offsets(rast, 0, 0); + rt_raster_set_scale(rast, 1, -1); + + band = addBand(rast, PT_32BUI, 1, 0); + CHECK(band); + + for (y = 0; y < maxY; y++) { + for (x = 0; x < maxX; x++) { + rt_band_set_pixel(band, x, y, 1); + } + } + + mpoly = rt_raster_surface(rast, 0); + CHECK((mpoly != NULL)); + wkt = lwgeom_to_wkt((const LWGEOM *) lwmpoly_as_lwgeom(mpoly), WKT_ISO, DBL_DIG, NULL); + RASTER_DEBUGF(4, "wkt = %s", wkt); + CHECK(!strcmp(wkt, "MULTIPOLYGON(((0 0,1 0,2 0,3 0,4 0,5 0,5 -1,5 -2,5 -3,5 -4,5 -5,4 -5,3 -5,2 -5,1 -5,0 -5,0 -4,0 -3,0 -2,0 -1,0 0)))")); + rtdealloc(wkt); + lwmpoly_free(mpoly); + + /* 0,0 NODATA */ + rt_band_set_pixel(band, 0, 0, 0); + + mpoly = rt_raster_surface(rast, 0); + CHECK((mpoly != NULL)); + wkt = lwgeom_to_wkt((const LWGEOM *) lwmpoly_as_lwgeom(mpoly), WKT_ISO, DBL_DIG, NULL); + RASTER_DEBUGF(4, "wkt = %s", wkt); + CHECK(!strcmp(wkt, "MULTIPOLYGON(((1 0,2 0,3 0,4 0,5 0,5 -1,5 -2,5 -3,5 -4,5 -5,4 -5,3 -5,2 -5,1 -5,0 -5,0 -4,0 -3,0 -2,0 -1,1 -1,1 0)))")); + rtdealloc(wkt); + lwmpoly_free(mpoly); + + /* plus 1,1 NODATA */ + rt_band_set_pixel(band, 1, 1, 0); + + mpoly = rt_raster_surface(rast, 0); + CHECK((mpoly != NULL)); + wkt = lwgeom_to_wkt((const LWGEOM *) lwmpoly_as_lwgeom(mpoly), WKT_ISO, DBL_DIG, NULL); + RASTER_DEBUGF(4, "wkt = %s", wkt); + CHECK(!strcmp(wkt, "MULTIPOLYGON(((1 0,2 0,3 0,4 0,5 0,5 -1,5 -2,5 -3,5 -4,5 -5,4 -5,3 -5,2 -5,1 -5,0 -5,0 -4,0 -3,0 -2,0 -1,1 -1,1 0),(1 -1,1 -2,2 -2,2 -1,1 -1)))")); + rtdealloc(wkt); + lwmpoly_free(mpoly); + + deepRelease(rast); +} + int main() { @@ -3168,6 +3225,10 @@ main() testPixelOfValue(); printf("OK\n"); + printf("Testing rt_raster_surface... "); + testRasterSurface(); + printf("OK\n"); + deepRelease(raster); return EXIT_SUCCESS;