From dea913f1232dcb3a2b94e7704bb48e32ec7f65c8 Mon Sep 17 00:00:00 2001 From: Bborie Park Date: Thu, 23 Feb 2012 21:49:51 +0000 Subject: [PATCH] Renamed rt_raster_dump_as_wktpolygons() to rt_raster_gdal_polygonize(). Refactored rt_raster_gdal_polygonize() to use LWPOLY objects instead of WKT strings. Also added cleanup code to make invalid dumped geometries valid. Associated tickets are #1586 and #637. git-svn-id: http://svn.osgeo.org/postgis/trunk@9280 b70326c6-7e19-0410-871a-916f4a2858ee --- raster/rt_core/rt_api.c | 468 +++++++++++++++----------- raster/rt_core/rt_api.h | 23 +- raster/rt_pg/rt_pg.c | 255 +++++++------- raster/rt_pg/rtpostgis.sql.in.c | 24 +- raster/rt_pg/rtpostgis_drop.sql.in.c | 5 + raster/test/core/testapi.c | 483 +++++++++++++++------------ 6 files changed, 679 insertions(+), 579 deletions(-) diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index 50d68fb44..87c6bfc92 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -1266,6 +1266,7 @@ rt_band_new_inline( band->width = width; band->height = height; band->hasnodata = hasnodata; + band->nodataval = 0; band->data.mem = data; band->ownsData = 0; band->isnodata = FALSE; @@ -1326,6 +1327,7 @@ rt_band_new_offline( band->width = width; band->height = height; band->hasnodata = hasnodata; + band->nodataval = 0; band->isnodata = FALSE; band->raster = NULL; @@ -5178,8 +5180,8 @@ rt_raster_geopoint_to_cell(rt_raster raster, * 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 the wkt - * representation of a geometry (one for each group of pixel sharing + * 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. * * @param raster: the raster to get info from. @@ -5187,258 +5189,324 @@ rt_raster_geopoint_to_cell(rt_raster raster, * * @return A set of "geomval" values, one for each group of pixels * sharing the same value for the provided band. The returned values are - * WKT geometries, not real PostGIS geometries (this may change in the - * future, and the function returns real geometries) + * LWPOLY geometries. */ rt_geomval -rt_raster_dump_as_wktpolygons(rt_raster raster, int nband, int * pnElements) -{ - char * pszQuery; - long j; - OGRSFDriverH ogr_drv = NULL; - GDALDriverH gdal_drv = NULL; - 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; - char * pszSrcText = NULL; - int nFeatureCount = 0; - rt_band band = NULL; - int iPixVal = -1; - double dValue = 0.0; - int iBandHasNodataValue = FALSE; - double dBandNoData = 0.0; - - uint32_t bandNums[1] = {nband}; +rt_raster_gdal_polygonize( + rt_raster raster, int nband, + int *pnElements +) { + CPLErr cplerr = CE_None; + char *pszQuery; + long j; + OGRSFDriverH ogr_drv = NULL; + GDALDriverH gdal_drv = NULL; + 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; - /* Checkings */ + /* for checking that a geometry is valid */ + GEOSGeom ggeom = NULL; + int msgValid = 0; + int isValid; + LWGEOM *lwgeomValid = NULL; + uint32_t bandNums[1] = {nband}; - assert(NULL != raster); - assert(nband >= 0 && nband < rt_raster_get_num_bands(raster)); + /* checks */ + assert(NULL != raster); + assert(nband >= 0 && nband < rt_raster_get_num_bands(raster)); - RASTER_DEBUG(2, "In rt_raster_dump_as_polygons"); + RASTER_DEBUG(2, "In rt_raster_gdal_polygonize"); - /******************************* - * Get band - *******************************/ - band = rt_raster_get_band(raster, nband); - if (NULL == band) { - rterror("rt_raster_dump_as_wktpolygons: Error getting band %d from raster", nband); - return 0; - } - - iBandHasNodataValue = rt_band_get_hasnodata_flag(band); - if (iBandHasNodataValue) dBandNoData = rt_band_get_nodata(band); + /******************************* + * 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; + } - /***************************************************** - * Convert raster to GDAL MEM dataset - *****************************************************/ - memdataset = rt_raster_to_gdal_mem(raster, NULL, bandNums, 1, &gdal_drv); - if (NULL == memdataset) { - rterror("rt_raster_dump_as_wktpolygons: Couldn't convert raster to GDAL MEM dataset"); + iBandHasNodataValue = rt_band_get_hasnodata_flag(band); + if (iBandHasNodataValue) dBandNoData = rt_band_get_nodata(band); - return 0; - } + /***************************************************** + * Convert raster to GDAL MEM dataset + *****************************************************/ + memdataset = rt_raster_to_gdal_mem(raster, NULL, bandNums, 1, &gdal_drv); + if (NULL == memdataset) { + rterror("rt_raster_gdal_polygonize: Couldn't convert raster to GDAL MEM dataset"); + return NULL; + } - /***************************** - * Register ogr mem driver - *****************************/ - OGRRegisterAll(); + /***************************** + * Register ogr mem driver + *****************************/ + OGRRegisterAll(); - RASTER_DEBUG(3, "creating OGR MEM vector"); + 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); + return NULL; + } - /***************************************************** - * Create an OGR in-memory vector for layers - *****************************************************/ - ogr_drv = OGRGetDriverByName("Memory"); - memdatasource = OGR_Dr_CreateDataSource(ogr_drv, "", 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"); - if (NULL == memdatasource) { - rterror("rt_raster_dump_as_wktpolygons: Couldn't create a OGR Datasource to store pols"); - GDALClose(memdataset); - return 0; - } + /* xxx jorgearevalo: what should we do now? */ + GDALClose(memdataset); + OGRReleaseDataSource(memdatasource); - /** - * Can MEM driver create new layers? - **/ - if (!OGR_DS_TestCapability(memdatasource, ODsCCreateLayer)) { - rterror("rt_raster_dump_as_wktpolygons: MEM driver can't create new layers, aborting"); - /* xxx jorgearevalo: what should we do now? */ - GDALClose(memdataset); - OGRReleaseDataSource(memdatasource); - return 0; - } + return NULL; + } - RASTER_DEBUG(3, "polygonizying GDAL MEM raster band"); + RASTER_DEBUG(3, "polygonizying GDAL MEM raster band"); - /***************************** - * Polygonize the raster band - *****************************/ + /***************************** + * Polygonize the raster band + *****************************/ - /** - * 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" - **/ - hLayer = OGR_DS_CreateLayer(memdatasource, "PolygonizedLayer", NULL, - wkbPolygon, NULL); + /** + * 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" + **/ + hLayer = OGR_DS_CreateLayer(memdatasource, "PolygonizedLayer", NULL, wkbPolygon, NULL); - if (NULL == hLayer) { - rterror("rt_raster_dump_as_wktpolygons: Couldn't create layer to store polygons"); - GDALClose(memdataset); - OGRReleaseDataSource(memdatasource); - return 0; - } + if (NULL == hLayer) { + rterror("rt_raster_gdal_polygonize: Couldn't create layer to store polygons"); - /** - * Create a new field in the layer, to store the px value - */ + GDALClose(memdataset); + OGRReleaseDataSource(memdatasource); - /* First, create a field definition to create the field */ - hFldDfn = OGR_Fld_Create("PixelValue", OFTReal); + return NULL; + } - /* 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; - } + /** + * Create a new field in the layer, to store the px value + */ + /* First, create a field definition to create the field */ + hFldDfn = OGR_Fld_Create("PixelValue", OFTReal); - /* Get GDAL raster band */ - gdal_band = GDALGetRasterBand(memdataset, 1); - if (NULL == gdal_band) { - rterror("rt_raster_dump_as_wktpolygons: Couldn't get GDAL band to polygonize"); - GDALClose(memdataset); + /* 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; + } - OGR_Fld_Destroy(hFldDfn); - OGR_DS_DeleteLayer(memdatasource, 0); - OGRReleaseDataSource(memdatasource); + /* 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"); - return 0; - } + GDALClose(memdataset); + 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. - **/ + /** + * We don't need a raster mask band. Each band has a nodata value. + **/ #ifdef GDALFPOLYGONIZE - GDALFPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL); + cplerr = GDALFPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL); #else - GDALPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL); + cplerr = GDALPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL); #endif - /** - * Optimization: Apply a OGR SQL filter to the layer to select the - * features different from NODATA value. - * - * Thanks to David Zwarg. - **/ - 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"); - } - } + if (cplerr != CE_None) { + rterror("rt_raster_gdal_polygonize: Could not polygonize GDAL band"); - else { - pszQuery = NULL; - } + GDALClose(memdataset); + OGR_Fld_Destroy(hFldDfn); + OGR_DS_DeleteLayer(memdatasource, 0); + OGRReleaseDataSource(memdatasource); + return NULL; + } + /** + * Optimization: Apply a OGR SQL filter to the layer to select the + * features different from NODATA value. + * + * Thanks to David Zwarg. + **/ + 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 in WKT 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. - * TODO: modify GDALPolygonize to directly emit polygons in WKT format? - *********************************************************************/ - nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE); + /********************************************************************* + * 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)); + /* Allocate memory for pols */ + pols = (rt_geomval) rtalloc(nFeatureCount * sizeof(struct rt_geomval_t)); - if (NULL == pols) { - rterror("rt_raster_dump_as_wktpolygons: Couldn't allocate memory for " - "geomval structure"); - GDALClose(memdataset); + if (NULL == pols) { + rterror("rt_raster_gdal_polygonize: Could not allocate memory for geomval set"); - OGR_Fld_Destroy(hFldDfn); - OGR_DS_DeleteLayer(memdatasource, 0); - if (NULL != pszQuery) rtdealloc(pszQuery); - OGRReleaseDataSource(memdatasource); + GDALClose(memdataset); + OGR_Fld_Destroy(hFldDfn); + OGR_DS_DeleteLayer(memdatasource, 0); + if (NULL != pszQuery) + rtdealloc(pszQuery); + OGRReleaseDataSource(memdatasource); - return 0; - } + return NULL; + } - RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount); + /* initialize GEOS */ + initGEOS(lwnotice, lwgeom_geos_error); - if (pnElements) - *pnElements = 0; + RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount); - /* Reset feature reading to start in the first feature */ - OGR_L_ResetReading(hLayer); + if (pnElements) + *pnElements = 0; - for (j = 0; j < nFeatureCount; j++) { + /* Reset feature reading to start in the first feature */ + OGR_L_ResetReading(hLayer); - hFeature = OGR_L_GetNextFeature(hLayer); - dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal); + for (j = 0; j < nFeatureCount; j++) { + hFeature = OGR_L_GetNextFeature(hLayer); + dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal); - hGeom = OGR_F_GetGeometryRef(hFeature); - OGR_G_ExportToWkt(hGeom, &pszSrcText); + hGeom = OGR_F_GetGeometryRef(hFeature); + wkbsize = OGR_G_WkbSize(hGeom); - pols[j].val = dValue; - pols[j].srid = rt_raster_get_srid(raster); - pols[j].geom = (char *) rtalloc((1 + strlen(pszSrcText)) - * sizeof (char)); - strcpy(pols[j].geom, pszSrcText); + /* allocate wkb buffer */ + wkb = rtalloc(sizeof(unsigned char) * wkbsize); + if (wkb == NULL) { + rterror("rt_raster_gdal_polygonize: Could not allocate memory for WKB buffer"); - RASTER_DEBUGF(4, "Feature %d, Polygon: %s", j, pols[j].geom); - RASTER_DEBUGF(4, "Feature %d, value: %f", j, pols[j].val); - RASTER_DEBUGF(4, "Feature %d, srid: %d", j, pols[j].srid); + OGR_F_Destroy(hFeature); + GDALClose(memdataset); + OGR_Fld_Destroy(hFldDfn); + OGR_DS_DeleteLayer(memdatasource, 0); + if (NULL != pszQuery) + rtdealloc(pszQuery); + OGRReleaseDataSource(memdatasource); - /** - * We can't use deallocator from rt_context, because it comes from - * postgresql backend, that uses pfree. This function needs a - * postgresql memory context to work with, and the memory created - * for pszSrcText is created outside this context. - **/ - /*rtdealloc(pszSrcText);*/ - free(pszSrcText); - pszSrcText = NULL; + return NULL; + } - OGR_F_Destroy(hFeature); - } + /* export WKB with LSB byte order */ + OGR_G_ExportToWkb(hGeom, wkbNDR, wkb); - if (pnElements) - *pnElements = nFeatureCount; + /* convert WKB to LWGEOM */ + lwgeom = lwgeom_from_wkb(wkb, wkbsize, LW_PARSER_CHECK_NONE); - RASTER_DEBUG(3, "destroying GDAL MEM raster"); + /* 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 { +#if POSTGIS_GEOS_VERSION < 33 + if (!msgValid) { + rtwarn("Skipping check for invalid geometry. GEOS-3.3.0 or up is required to fix an invalid geometry"); + msgValid = 1; + } + break; +#endif + + ggeom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom); + if (ggeom == NULL) { + rtwarn("Cannot test geometry for validity"); + break; + } + + isValid = GEOSisValid(ggeom); + + GEOSGeom_destroy(ggeom); + ggeom = NULL; + + /* geometry is valid */ + if (isValid) + break; + + /* 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); + + /* set pixel value */ + pols[j].val = dValue; + } - GDALClose(memdataset); + if (pnElements) + *pnElements = nFeatureCount; - RASTER_DEBUG(3, "destroying OGR MEM vector"); + RASTER_DEBUG(3, "destroying GDAL MEM raster"); + GDALClose(memdataset); - OGR_Fld_Destroy(hFldDfn); - OGR_DS_DeleteLayer(memdatasource, 0); - if (NULL != pszQuery) rtdealloc(pszQuery); - OGRReleaseDataSource(memdatasource); + 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; + return pols; } LWPOLY* @@ -9886,7 +9954,7 @@ rt_raster_intersects( rtn = 0; break; } - ghull[i] = (GEOSGeometry *) LWGEOM2GEOS((LWGEOM *) hull[i]); + ghull[i] = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(hull[i])); if (NULL == ghull[i]) { for (j = 0; j < i; j++) { GEOSGeom_destroy(ghull[j]); diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index f57ba5529..ac26d11bb 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -1041,13 +1041,12 @@ LWPOLY* rt_raster_get_convex_hull(rt_raster raster); */ LWPOLY* rt_raster_pixel_as_polygon(rt_raster raster, int x, int y); - /** * 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 the wkt - * representation of a geometry (one for each group of pixel sharing + * 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. * * @param raster: the raster to get info from. @@ -1055,13 +1054,13 @@ LWPOLY* rt_raster_pixel_as_polygon(rt_raster raster, int x, int y); * * @return A set of "geomval" values, one for each group of pixels * sharing the same value for the provided band. The returned values are - * WKT geometries, not real PostGIS geometries (this may change in the - * future, and the function returns real geometries) + * LWPOLY geometries. */ rt_geomval -rt_raster_dump_as_wktpolygons(rt_raster raster, int nband, - int * pnElements); - +rt_raster_gdal_polygonize( + rt_raster raster, int nband, + int * pnElements +); /** * Return this raster in serialized form. @@ -1568,12 +1567,10 @@ struct rt_band_t { }; -/* WKT string representing each polygon in WKT format accompanied by its -corresponding value */ +/* polygon as LWPOLY with associated value */ struct rt_geomval_t { - int srid; - double val; - char * geom; + LWPOLY *geom; + double val; }; /* summary stats of specified band */ diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index 7799c92fd..36348fbd1 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -127,7 +127,7 @@ static char *rtpg_getSR(int srid); * - Use SRF funcapi, and storing the data at multi_call_memory_ctx (by pointing * the data you want to keep with funcctx->user_fctx. funcctx is created by * funcctx = SPI_FIRSTCALL_INIT()). Recommended way in functions returning rows, - * like RASTER_dumpAsWKTPolygons (see section 34.9.9 at + * like RASTER_dumpAsPolygons (see section 34.9.9 at * http://www.postgresql.org/docs/8.4/static/xfunc-c.html). * * But raster code follows the same philosophy than the rest of PostGIS: keep @@ -153,7 +153,7 @@ Datum RASTER_to_binary(PG_FUNCTION_ARGS); /* Raster as geometry operations */ Datum RASTER_convex_hull(PG_FUNCTION_ARGS); -Datum RASTER_dumpAsWKTPolygons(PG_FUNCTION_ARGS); +Datum RASTER_dumpAsPolygons(PG_FUNCTION_ARGS); /* Get all the properties of a raster */ Datum RASTER_getSRID(PG_FUNCTION_ARGS); @@ -809,155 +809,150 @@ Datum RASTER_convex_hull(PG_FUNCTION_ARGS) } -PG_FUNCTION_INFO_V1(RASTER_dumpAsWKTPolygons); -Datum RASTER_dumpAsWKTPolygons(PG_FUNCTION_ARGS) -{ - rt_pgraster *pgraster = NULL; - rt_raster raster = NULL; - FuncCallContext *funcctx; - TupleDesc tupdesc; - int nband; - rt_geomval geomval; - rt_geomval geomval2; - int call_cntr; - int max_calls; - int nElements; - MemoryContext oldcontext; - - /* stuff done only on the first call of the function */ - if (SRF_IS_FIRSTCALL()) - { - int numbands; - - POSTGIS_RT_DEBUG(2, "RASTER_dumpAsWKTPolygons first call"); - - /* create a function context for cross-call persistence */ - funcctx = SRF_FIRSTCALL_INIT(); - - /* switch to memory context appropriate for multiple function calls */ - oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx); - - /* Get input arguments */ - if (PG_ARGISNULL(0)) SRF_RETURN_DONE(funcctx); - pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - raster = rt_raster_deserialize(pgraster, FALSE); - if ( ! raster ) - { - ereport(ERROR, - (errcode(ERRCODE_OUT_OF_MEMORY), - errmsg("Could not deserialize raster"))); - MemoryContextSwitchTo(oldcontext); - SRF_RETURN_DONE(funcctx); - } +PG_FUNCTION_INFO_V1(RASTER_dumpAsPolygons); +Datum RASTER_dumpAsPolygons(PG_FUNCTION_ARGS) { + FuncCallContext *funcctx; + TupleDesc tupdesc; + rt_geomval geomval; + rt_geomval geomval2; + int call_cntr; + int max_calls; - if (PG_NARGS() == 2) - nband = PG_GETARG_UINT32(1); - else - nband = 1; /* By default, first band */ + /* stuff done only on the first call of the function */ + if (SRF_IS_FIRSTCALL()) { + MemoryContext oldcontext; + int numbands; + rt_pgraster *pgraster = NULL; + rt_raster raster = NULL; + int nband; + int nElements; - POSTGIS_RT_DEBUGF(3, "band %d", nband); + POSTGIS_RT_DEBUG(2, "RASTER_dumpAsPolygons first call"); - numbands = rt_raster_get_num_bands(raster); - if (nband < 1 || nband > numbands) { - elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL"); - rt_raster_destroy(raster); - MemoryContextSwitchTo(oldcontext); - SRF_RETURN_DONE(funcctx); - } + /* create a function context for cross-call persistence */ + funcctx = SRF_FIRSTCALL_INIT(); - /* Polygonize raster */ + /* switch to memory context appropriate for multiple function calls */ + oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx); - /** - * Dump raster - */ - geomval = rt_raster_dump_as_wktpolygons(raster, nband - 1, &nElements); - rt_raster_destroy(raster); - if (NULL == geomval) - { - ereport(ERROR, - (errcode(ERRCODE_NO_DATA_FOUND), - errmsg("Could not polygonize raster"))); - rt_raster_destroy(raster); - MemoryContextSwitchTo(oldcontext); - SRF_RETURN_DONE(funcctx); - } + /* Get input arguments */ + if (PG_ARGISNULL(0)) { + MemoryContextSwitchTo(oldcontext); + SRF_RETURN_DONE(funcctx); + } + pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + raster = rt_raster_deserialize(pgraster, FALSE); + if (!raster) { + ereport(ERROR, ( + errcode(ERRCODE_OUT_OF_MEMORY), + errmsg("Could not deserialize raster") + )); + MemoryContextSwitchTo(oldcontext); + SRF_RETURN_DONE(funcctx); + } + + if (PG_NARGS() == 2) + nband = PG_GETARG_UINT32(1); + else + nband = 1; /* By default, first band */ + + POSTGIS_RT_DEBUGF(3, "band %d", nband); + numbands = rt_raster_get_num_bands(raster); + + if (nband < 1 || nband > numbands) { + elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL"); + rt_raster_destroy(raster); + MemoryContextSwitchTo(oldcontext); + SRF_RETURN_DONE(funcctx); + } - POSTGIS_RT_DEBUGF(3, "raster dump, %d elements returned", nElements); + /* Polygonize raster */ - /* Store needed information */ - funcctx->user_fctx = geomval; + /** + * Dump raster + */ + geomval = rt_raster_gdal_polygonize(raster, nband - 1, &nElements); + rt_raster_destroy(raster); + if (NULL == geomval) { + ereport(ERROR, ( + errcode(ERRCODE_NO_DATA_FOUND), + errmsg("Could not polygonize raster") + )); + MemoryContextSwitchTo(oldcontext); + SRF_RETURN_DONE(funcctx); + } - /* total number of tuples to be returned */ - funcctx->max_calls = nElements; + POSTGIS_RT_DEBUGF(3, "raster dump, %d elements returned", nElements); - /* Build a tuple descriptor for our result type */ - if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) - ereport(ERROR, - (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), - errmsg("function returning record called in context " - "that cannot accept type record"))); + /* Store needed information */ + funcctx->user_fctx = geomval; - BlessTupleDesc(tupdesc); - funcctx->tuple_desc = tupdesc; + /* total number of tuples to be returned */ + funcctx->max_calls = nElements; - MemoryContextSwitchTo(oldcontext); - } + /* Build a tuple descriptor for our result type */ + if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) { + ereport(ERROR, ( + errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("function returning record called in context that cannot accept type record") + )); + } - /* stuff done on every call of the function */ - funcctx = SRF_PERCALL_SETUP(); + BlessTupleDesc(tupdesc); + funcctx->tuple_desc = tupdesc; - call_cntr = funcctx->call_cntr; - max_calls = funcctx->max_calls; - tupdesc = funcctx->tuple_desc; - geomval2 = funcctx->user_fctx; + MemoryContextSwitchTo(oldcontext); + } - if (call_cntr < max_calls) /* do when there is more left to send */ - { - bool *nulls = NULL; - int values_length = 3; - Datum values[values_length]; - HeapTuple tuple; - Datum result; + /* stuff done on every call of the function */ + funcctx = SRF_PERCALL_SETUP(); - POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr); + call_cntr = funcctx->call_cntr; + max_calls = funcctx->max_calls; + tupdesc = funcctx->tuple_desc; + geomval2 = funcctx->user_fctx; - nulls = palloc(sizeof(bool) * values_length); - memset(nulls, FALSE, values_length); + /* do when there is more left to send */ + if (call_cntr < max_calls) { + bool *nulls = NULL; + int values_length = 2; + Datum values[values_length]; + HeapTuple tuple; + Datum result; - values[0] = CStringGetTextDatum(geomval2[call_cntr].geom); - values[1] = Float8GetDatum(geomval2[call_cntr].val); - values[2] = Int32GetDatum(geomval2[call_cntr].srid); + GSERIALIZED *gser = NULL; + size_t gser_size = 0; - POSTGIS_RT_DEBUGF(4, "Result %d, Polygon %s", call_cntr, geomval2[call_cntr].geom); - POSTGIS_RT_DEBUGF(4, "Result %d, val %f", call_cntr, geomval2[call_cntr].val); - POSTGIS_RT_DEBUGF(4, "Result %d, srid %d", call_cntr, geomval2[call_cntr].srid); + POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr); - /** - * Free resources. - */ - pfree(geomval2[call_cntr].geom); + nulls = palloc(sizeof(bool) * values_length); + memset(nulls, FALSE, values_length); - /* build a tuple */ - tuple = heap_form_tuple(tupdesc, values, nulls); + /* convert LWGEOM to GSERIALIZED */ + gser = gserialized_from_lwgeom(lwpoly_as_lwgeom(geomval2[call_cntr].geom), 0, &gser_size); + lwgeom_free(lwpoly_as_lwgeom(geomval2[call_cntr].geom)); - /* make the tuple into a datum */ - result = HeapTupleGetDatum(tuple); + values[0] = PointerGetDatum(gser); + values[1] = Float8GetDatum(geomval2[call_cntr].val); - /* clean up (this is not really necessary) */ - pfree(nulls); + /* build a tuple */ + tuple = heap_form_tuple(tupdesc, values, nulls); + /* make the tuple into a datum */ + result = HeapTupleGetDatum(tuple); - SRF_RETURN_NEXT(funcctx, result); - } - else /* do when there is no more left */ - { - pfree(geomval2); - SRF_RETURN_DONE(funcctx); - } + /* clean up (this is not really necessary) */ + pfree(nulls); + SRF_RETURN_NEXT(funcctx, result); + } + /* do when there is no more left */ + else { + pfree(geomval2); + SRF_RETURN_DONE(funcctx); + } } - /** * rt_MakeEmptyRaster( , , , , * , , @@ -1613,8 +1608,10 @@ Datum RASTER_getGeotransform(PG_FUNCTION_ARGS) double jmag; double theta_i; double theta_ij; + /* double xoffset; double yoffset; + */ TupleDesc result_tuple; /* for returning a composite */ bool *nulls = NULL; @@ -7168,7 +7165,7 @@ Datum RASTER_getGDALDrivers(PG_FUNCTION_ARGS) PG_FUNCTION_INFO_V1(RASTER_asRaster); Datum RASTER_asRaster(PG_FUNCTION_ARGS) { - GSERIALIZED *pggeom = NULL; + GSERIALIZED *gser = NULL; LWGEOM *geom = NULL; rt_raster rast = NULL; @@ -7242,8 +7239,8 @@ Datum RASTER_asRaster(PG_FUNCTION_ARGS) if (PG_ARGISNULL(0)) PG_RETURN_NULL(); - pggeom = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - geom = lwgeom_from_gserialized(pggeom); + gser = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + geom = lwgeom_from_gserialized(gser); /* Get a 2D version of the geometry if necessary */ if (lwgeom_ndims(geom) > 2) { @@ -7655,7 +7652,7 @@ Datum RASTER_asRaster(PG_FUNCTION_ARGS) } /* get geometry's srid */ - srid = gserialized_get_srid(pggeom); + srid = gserialized_get_srid(gser); POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: srid = %d", srid); if (clamp_srid(srid) != SRID_UNKNOWN) { @@ -8447,7 +8444,7 @@ Datum RASTER_intersects(PG_FUNCTION_ARGS) rtn = 0; break; } - ghull[i] = (GEOSGeometry *) LWGEOM2GEOS((LWGEOM *) hull[i]); + ghull[i] = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(hull[i])); if (NULL == ghull[i]) { for (j = 0; j < i; j++) { GEOSGeom_destroy(ghull[j]); diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index 0d787ba28..4b3f7655f 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -2425,28 +2425,14 @@ CREATE OR REPLACE FUNCTION st_setvalue(rast raster, pt geometry, newvalue float8 ----------------------------------------------------------------------- CREATE TYPE geomval AS ( - geom geometry, - val double precision -); - -CREATE TYPE wktgeomval AS ( - wktgeom text, - val double precision, - srid int + geom geometry, + val double precision ); -CREATE OR REPLACE FUNCTION _st_dumpaswktpolygons(rast raster, band integer) - RETURNS SETOF wktgeomval - AS 'MODULE_PATHNAME','RASTER_dumpAsWKTPolygons' - LANGUAGE 'C' IMMUTABLE STRICT; - CREATE OR REPLACE FUNCTION st_dumpaspolygons(rast raster, band integer DEFAULT 1) - RETURNS SETOF geomval AS - $$ - SELECT st_geomfromtext(wktgeomval.wktgeom, wktgeomval.srid), wktgeomval.val - FROM _st_dumpaswktpolygons($1, $2) AS wktgeomval; - $$ - LANGUAGE 'SQL' IMMUTABLE STRICT; + RETURNS SETOF geomval + AS 'MODULE_PATHNAME','RASTER_dumpAsPolygons' + LANGUAGE 'C' IMMUTABLE STRICT; CREATE OR REPLACE FUNCTION st_polygon(rast raster, band integer DEFAULT 1) RETURNS geometry AS diff --git a/raster/rt_pg/rtpostgis_drop.sql.in.c b/raster/rt_pg/rtpostgis_drop.sql.in.c index c5e257096..7f41bc0f2 100644 --- a/raster/rt_pg/rtpostgis_drop.sql.in.c +++ b/raster/rt_pg/rtpostgis_drop.sql.in.c @@ -310,3 +310,8 @@ DROP FUNCTION IF EXISTS st_intersection(raster, int, geometry, text, regprocedur DROP FUNCTION IF EXISTS st_intersection(raster, int, geometry, regprocedure); DROP FUNCTION IF EXISTS st_intersection(raster, geometry, text, regprocedure); DROP FUNCTION IF EXISTS st_intersection(raster, geometry, regprocedure); + +-- refactoring of functions +DROP TYPE IF EXISTS wktgeomval; +DROP FUNCTION IF EXISTS _st_dumpaswktpolygons(rast raster, band integer); + diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c index dc0479628..d12997290 100644 --- a/raster/test/core/testapi.c +++ b/raster/test/core/testapi.c @@ -7,10 +7,15 @@ #include "rt_api.h" #include "check.h" -static rt_band addBand(rt_raster raster, rt_pixtype pixtype, int hasnodata, double nodataval); -static void deepRelease(rt_raster raster); -static void testBand1BB(rt_band band); -static rt_raster fillRasterToPolygonize(int hasnodata, double nodatavalue); +static char * +lwgeom_to_text(const LWGEOM *lwgeom) { + char *wkt; + size_t wkt_size; + + wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, &wkt_size); + + return wkt; +} static rt_band addBand(rt_raster raster, rt_pixtype pixtype, int hasnodata, double nodataval) @@ -159,6 +164,259 @@ fillRasterToPolygonize(int hasnodata, double nodatavalue) return raster; } +static void testGDALPolygonize() { + int i; + rt_raster rt; + int nPols = 0; + rt_geomval gv = NULL; + char *wkt = NULL; + + rt = fillRasterToPolygonize(1, -1.0); + CHECK(!rt_raster_has_no_band(rt, 0)); + + nPols = 0; + gv = rt_raster_gdal_polygonize(rt, 0, &nPols); + + /* + for (i = 0; i < nPols; i++) { + wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom); + printf("(i, val, geom) = (%d, %f, '%s')\n", i, gv[i].val, wkt); + rtdealloc(wkt); + } + */ + +#ifdef GDALFPOLYGONIZE + CHECK(FLT_EQ(gv[0].val, 1.8)); +#else + CHECK(FLT_EQ(gv[0].val, 2.0)); +#endif + + wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom); + CHECK(!strcmp(wkt, "POLYGON((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))")); + rtdealloc(wkt); + + CHECK_EQUALS_DOUBLE(gv[1].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom); + CHECK(!strcmp(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))")); + rtdealloc(wkt); + +#ifdef GDALFPOLYGONIZE + CHECK(FLT_EQ(gv[2].val, 2.8)); +#else + CHECK(FLT_EQ(gv[2].val, 3.0)); +#endif + + wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom); + CHECK(!strcmp(wkt, "POLYGON((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))")); + rtdealloc(wkt); + + CHECK_EQUALS_DOUBLE(gv[3].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom); + CHECK(!strcmp(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); + rtdealloc(wkt); + + for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom); + rtdealloc(gv); + deepRelease(rt); + + /* Second test: NODATA value = 1.8 */ +#ifdef GDALFPOLYGONIZE + rt = fillRasterToPolygonize(1, 1.8); +#else + rt = fillRasterToPolygonize(1, 2.0); +#endif + + /* We can check rt_raster_has_no_band here too */ + CHECK(!rt_raster_has_no_band(rt, 0)); + + nPols = 0; + gv = rt_raster_gdal_polygonize(rt, 0, &nPols); + + /* + for (i = 0; i < nPols; i++) { + wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom); + printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt); + rtdealloc(wkt); + } + */ + +#ifdef GDALFPOLYGONIZE + CHECK_EQUALS_DOUBLE(gv[1].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom); + CHECK(!strcmp(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))")); + rtdealloc(wkt); + + CHECK(FLT_EQ(gv[2].val, 2.8)); + wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom); + CHECK(!strcmp(wkt, "POLYGON((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))")); + rtdealloc(wkt); + + CHECK_EQUALS_DOUBLE(gv[3].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom); + CHECK(!strcmp(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); + rtdealloc(wkt); +#else + CHECK_EQUALS_DOUBLE(gv[0].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom); + CHECK(!strcmp(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))")); + rtdealloc(wkt); + + CHECK(FLT_EQ(gv[1].val, 3.0)); + wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom); + CHECK(!strcmp(wkt, "POLYGON((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))")); + rtdealloc(wkt); + + CHECK_EQUALS_DOUBLE(gv[2].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom); + CHECK(!strcmp(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); + rtdealloc(wkt); +#endif + + for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom); + rtdealloc(gv); + deepRelease(rt); + + /* Third test: NODATA value = 2.8 */ +#ifdef GDALFPOLYGONIZE + rt = fillRasterToPolygonize(1, 2.8); +#else + rt = fillRasterToPolygonize(1, 3.0); +#endif + + /* We can check rt_raster_has_no_band here too */ + CHECK(!rt_raster_has_no_band(rt, 0)); + + nPols = 0; + gv = rt_raster_gdal_polygonize(rt, 0, &nPols); + + /* + for (i = 0; i < nPols; i++) { + wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom); + printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt); + rtdealloc(wkt); + } + */ + +#ifdef GDALFPOLYGONIZE + CHECK(FLT_EQ(gv[0].val, 1.8)); + + CHECK_EQUALS_DOUBLE(gv[3].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom); + CHECK(!strcmp(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); + rtdealloc(wkt); +#else + CHECK(FLT_EQ(gv[0].val, 2.0)); + + CHECK_EQUALS_DOUBLE(gv[2].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom); + CHECK(!strcmp(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); + rtdealloc(wkt); +#endif + + wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom); + CHECK(!strcmp(wkt, "POLYGON((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))")); + rtdealloc(wkt); + + CHECK_EQUALS_DOUBLE(gv[1].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom); + CHECK(!strcmp(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))")); + rtdealloc(wkt); + + for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom); + rtdealloc(gv); + deepRelease(rt); + + /* Fourth test: NODATA value = 0 */ + rt = fillRasterToPolygonize(1, 0.0); + /* We can check rt_raster_has_no_band here too */ + CHECK(!rt_raster_has_no_band(rt, 0)); + + nPols = 0; + gv = rt_raster_gdal_polygonize(rt, 0, &nPols); + + /* + for (i = 0; i < nPols; i++) { + wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom); + printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt); + rtdealloc(wkt); + } + */ + +#ifdef GDALFPOLYGONIZE + CHECK(FLT_EQ(gv[0].val, 1.8)); +#else + CHECK(FLT_EQ(gv[0].val, 2.0)); +#endif + + wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom); + CHECK(!strcmp(wkt, "POLYGON((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))")); + rtdealloc(wkt); + +#ifdef GDALFPOLYGONIZE + CHECK(FLT_EQ(gv[1].val, 2.8)); +#else + CHECK(FLT_EQ(gv[1].val, 3.0)); +#endif + + wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom); + CHECK(!strcmp(wkt, "POLYGON((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))")); + rtdealloc(wkt); + + for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom); + rtdealloc(gv); + deepRelease(rt); + + /* Last test: There is no NODATA value (all values are valid) */ + rt = fillRasterToPolygonize(0, 0.0); + /* We can check rt_raster_has_no_band here too */ + CHECK(!rt_raster_has_no_band(rt, 0)); + + nPols = 0; + gv = rt_raster_gdal_polygonize(rt, 0, &nPols); + + /* + for (i = 0; i < nPols; i++) { + wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom); + printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt); + rtdealloc(wkt); + } + */ + +#ifdef GDALFPOLYGONIZE + CHECK(FLT_EQ(gv[0].val, 1.8)); +#else + CHECK(FLT_EQ(gv[0].val, 2.0)); +#endif + + wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom); + CHECK(!strcmp(wkt, "POLYGON((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))")); + rtdealloc(wkt); + + CHECK_EQUALS_DOUBLE(gv[1].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom); + CHECK(!strcmp(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))")); + rtdealloc(wkt); + +#ifdef GDALFPOLYGONIZE + CHECK(FLT_EQ(gv[2].val, 2.8)); +#else + CHECK(FLT_EQ(gv[2].val, 3.0)); +#endif + + wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom); + CHECK(!strcmp(wkt, "POLYGON((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))")); + rtdealloc(wkt); + + CHECK_EQUALS_DOUBLE(gv[3].val, 0.0); + wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom); + CHECK(!strcmp(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); + rtdealloc(wkt); + + for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom); + rtdealloc(gv); + deepRelease(rt); +} + static void testBand1BB(rt_band band) { int failure; @@ -2298,7 +2556,6 @@ static void testLoadOfflineBand() { int main() { - int i; rt_raster raster; rt_band band_1BB, band_2BUI, band_4BUI, band_8BSI, band_8BUI, band_16BSI, band_16BUI, @@ -2431,219 +2688,9 @@ main() rt_raster_set_skews(raster, 0, 0); } - { /* Check ST_AsPolygon */ - printf("Testing polygonize function\n"); - - /* First test: NODATA value = -1 */ - rt_raster rt = fillRasterToPolygonize(1, -1.0); - - /* We can check rt_raster_has_no_band here too */ - CHECK(!rt_raster_has_no_band(rt, 0)); - - int nPols = 0; - - rt_geomval gv = (rt_geomval) rt_raster_dump_as_wktpolygons(rt, 0, &nPols); - - /* - for (i = 0; i < nPols; i++) { - printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, gv[i].geom); - } - */ - - -#ifdef GDALFPOLYGONIZE - CHECK(FLT_EQ(gv[0].val, 1.8)); -#else - CHECK(FLT_EQ(gv[0].val, 2.0)); -#endif - - CHECK(!strcmp(gv[0].geom, "POLYGON ((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))")); - - CHECK_EQUALS_DOUBLE(gv[1].val, 0.0); - CHECK(!strcmp(gv[1].geom, "POLYGON ((3 3,3 6,6 6,6 3,3 3))")); - -#ifdef GDALFPOLYGONIZE - CHECK(FLT_EQ(gv[2].val, 2.8)); -#else - CHECK(FLT_EQ(gv[2].val, 3.0)); -#endif - - CHECK(!strcmp(gv[2].geom, "POLYGON ((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))")); - - CHECK_EQUALS_DOUBLE(gv[3].val, 0.0); - CHECK(!strcmp(gv[3].geom, "POLYGON ((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); - - - for (i = 0; i < nPols; i++) rtdealloc(gv[i].geom); - rtdealloc(gv); - deepRelease(rt); - - - /* Second test: NODATA value = 1.8 */ -#ifdef GDALFPOLYGONIZE - rt = fillRasterToPolygonize(1, 1.8); -#else - rt = fillRasterToPolygonize(1, 2.0); -#endif - - - /* We can check rt_raster_has_no_band here too */ - CHECK(!rt_raster_has_no_band(rt, 0)); - - nPols = 0; - - gv = (rt_geomval) rt_raster_dump_as_wktpolygons(rt, 0, &nPols); - - /* - for (i = 0; i < nPols; i++) { - printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, gv[i].geom); - } - */ - -#ifdef GDALFPOLYGONIZE - CHECK_EQUALS_DOUBLE(gv[1].val, 0.0); - CHECK(!strcmp(gv[1].geom, "POLYGON ((3 3,3 6,6 6,6 3,3 3))")); - - CHECK(FLT_EQ(gv[2].val, 2.8)); - CHECK(!strcmp(gv[2].geom, "POLYGON ((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))")); - - CHECK_EQUALS_DOUBLE(gv[3].val, 0.0); - CHECK(!strcmp(gv[3].geom, "POLYGON ((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); -#else - CHECK_EQUALS_DOUBLE(gv[0].val, 0.0); - CHECK(!strcmp(gv[0].geom, "POLYGON ((3 3,3 6,6 6,6 3,3 3))")); - - CHECK(FLT_EQ(gv[1].val, 3.0)); - CHECK(!strcmp(gv[1].geom, "POLYGON ((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))")); - - CHECK_EQUALS_DOUBLE(gv[2].val, 0.0); - CHECK(!strcmp(gv[2].geom, "POLYGON ((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); -#endif - - for (i = 0; i < nPols; i++) rtdealloc(gv[i].geom); - rtdealloc(gv); - deepRelease(rt); - - /* Third test: NODATA value = 2.8 */ -#ifdef GDALFPOLYGONIZE - rt = fillRasterToPolygonize(1, 2.8); -#else - rt = fillRasterToPolygonize(1, 3.0); -#endif - - /* We can check rt_raster_has_no_band here too */ - CHECK(!rt_raster_has_no_band(rt, 0)); - - nPols = 0; - - gv = (rt_geomval) rt_raster_dump_as_wktpolygons(rt, 0, &nPols); - - /* - for (i = 0; i < nPols; i++) { - printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, gv[i].geom); - } - */ - -#ifdef GDALFPOLYGONIZE - CHECK(FLT_EQ(gv[0].val, 1.8)); - - CHECK_EQUALS_DOUBLE(gv[3].val, 0.0); - CHECK(!strcmp(gv[3].geom, "POLYGON ((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); -#else - CHECK(FLT_EQ(gv[0].val, 2.0)); - - CHECK_EQUALS_DOUBLE(gv[2].val, 0.0); - CHECK(!strcmp(gv[2].geom, "POLYGON ((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); -#endif - - CHECK(!strcmp(gv[0].geom, "POLYGON ((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))")); - - CHECK_EQUALS_DOUBLE(gv[1].val, 0.0); - CHECK(!strcmp(gv[1].geom, "POLYGON ((3 3,3 6,6 6,6 3,3 3))")); - - for (i = 0; i < nPols; i++) rtdealloc(gv[i].geom); - rtdealloc(gv); - deepRelease(rt); - - /* Fourth test: NODATA value = 0 */ - rt = fillRasterToPolygonize(1, 0.0); - - /* We can check rt_raster_has_no_band here too */ - CHECK(!rt_raster_has_no_band(rt, 0)); - - nPols = 0; - - gv = (rt_geomval) rt_raster_dump_as_wktpolygons(rt, 0, &nPols); - - /* - for (i = 0; i < nPols; i++) { - printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, gv[i].geom); - } - */ - -#ifdef GDALFPOLYGONIZE - CHECK(FLT_EQ(gv[0].val, 1.8)); -#else - CHECK(FLT_EQ(gv[0].val, 2.0)); -#endif - - CHECK(!strcmp(gv[0].geom, "POLYGON ((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))")); - -#ifdef GDALFPOLYGONIZE - CHECK(FLT_EQ(gv[1].val, 2.8)); -#else - CHECK(FLT_EQ(gv[1].val, 3.0)); -#endif - - CHECK(!strcmp(gv[1].geom, "POLYGON ((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))")); - - for (i = 0; i < nPols; i++) rtdealloc(gv[i].geom); - rtdealloc(gv); - deepRelease(rt); - - /* Last test: There is no NODATA value (all values are valid) */ - rt = fillRasterToPolygonize(0, 0.0); - - /* We can check rt_raster_has_no_band here too */ - CHECK(!rt_raster_has_no_band(rt, 0)); - - nPols = 0; - - gv = (rt_geomval) rt_raster_dump_as_wktpolygons(rt, 0, &nPols); - - /* - for (i = 0; i < nPols; i++) { - printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, gv[i].geom); - } - */ - -#ifdef GDALFPOLYGONIZE - CHECK(FLT_EQ(gv[0].val, 1.8)); -#else - CHECK(FLT_EQ(gv[0].val, 2.0)); -#endif - - CHECK(!strcmp(gv[0].geom, "POLYGON ((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))")); - - CHECK_EQUALS_DOUBLE(gv[1].val, 0.0); - CHECK(!strcmp(gv[1].geom, "POLYGON ((3 3,3 6,6 6,6 3,3 3))")); - -#ifdef GDALFPOLYGONIZE - CHECK(FLT_EQ(gv[2].val, 2.8)); -#else - CHECK(FLT_EQ(gv[2].val, 3.0)); -#endif - - CHECK(!strcmp(gv[2].geom, "POLYGON ((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))")); - - CHECK_EQUALS_DOUBLE(gv[3].val, 0.0); - CHECK(!strcmp(gv[3].geom, "POLYGON ((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))")); - - for (i = 0; i < nPols; i++) rtdealloc(gv[i].geom); - rtdealloc(gv); - deepRelease(rt); - - } + printf("Testing rt_raster_gdal_polygonize\n"); + testGDALPolygonize(); + printf("Successfully tested rt_raster_gdal_polygonize\n"); printf("Testing 1BB band\n"); band_1BB = addBand(raster, PT_1BB, 0, 0); -- 2.40.0