From: Bborie Park Date: Fri, 3 Jun 2011 17:02:57 +0000 (+0000) Subject: Added rt_raster_from_gdal_dataset function to rt_core/rt_api.c. This function provid... X-Git-Tag: 2.0.0alpha1~1497 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=8e1a828f0f8ae4fc33abc1786fc2c1fc0c948b8d;p=postgis Added rt_raster_from_gdal_dataset function to rt_core/rt_api.c. This function provides the functionality needed to convert a GDAL dataset into a PostGIS Raster object. git-svn-id: http://svn.osgeo.org/postgis/trunk@7317 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index ea23f4456..27cd1d0f2 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -5455,13 +5455,13 @@ rt_raster_gdal_drivers(uint32_t *drv_count) { } /** - * Return GDAL datasource using GDAL MEM driver from raster + * Return GDAL dataset using GDAL MEM driver from raster * * @param raster : raster to convert to GDAL MEM * @param srs : the raster's coordinate system in OGC WKT or PROJ.4 * @param rtn_drv : is set to the GDAL driver object * - * @return GDAL datasource using GDAL MEM driver + * @return GDAL dataset using GDAL MEM driver */ GDALDatasetH rt_raster_to_gdal_mem(rt_raster raster, char *srs, @@ -5651,3 +5651,258 @@ rt_raster_to_gdal_mem(rt_raster raster, char *srs, return ds; } + +/** + * Return a raster from a GDAL dataset + * + * @param ds : the GDAL dataset to convert to a raster + * + * @return raster + */ +rt_raster +rt_raster_from_gdal_dataset(GDALDatasetH ds) { + rt_raster rast = NULL; + double gt[6] = {0}; + CPLErr cplerr; + uint32_t width = 0; + uint32_t height = 0; + uint32_t numBands = 0; + int i = 0; + int *status = NULL; + + GDALRasterBandH gdband = NULL; + GDALDataType gdpixtype = GDT_Unknown; + rt_band band; + int32_t idx; + rt_pixtype pt = PT_END; + uint32_t hasnodata = 0; + double nodataval; + + int x; + int y; + double value; + + int nXBlocks, nYBlocks; + int nXBlockSize, nYBlockSize; + int iXBlock, iYBlock; + int nXValid, nYValid; + int iY, iX; + double *values = NULL; + int byblock = 0; + + assert(NULL != ds); + + /* raster size */ + width = GDALGetRasterXSize(ds); + height = GDALGetRasterYSize(ds); + RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d", width, height); + + /* create new raster */ + RASTER_DEBUG(3, "Creating new raster"); + rast = rt_raster_new(width, height); + if (NULL == rast) { + rterror("rt_raster_from_gdal_dataset: Out of memory allocating new raster\n"); + return NULL; + } + + /* get raster attributes */ + cplerr = GDALGetGeoTransform(ds, gt); + if (cplerr != CE_None) { + rterror("rt_raster_from_gdal_dataset: Unable to get geotransformation\n"); + rt_raster_destroy(rast); + return NULL; + } + RASTER_DEBUGF(3, "Raster geotransform (%f, %f, %f, %f, %f, %f)", + gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]); + + /* apply raster attributes */ + /* scale */ + rt_raster_set_scale(rast, gt[1], gt[5]); + /* offset */ + rt_raster_set_offsets(rast, gt[0], gt[3]); + /* skew */ + rt_raster_set_skews(rast, gt[2], gt[4]); + /* srid */ + /* CANNOT SET srid AS srid IS NOT AVAILABLE IN gdal dataset */ + + /* copy bands */ + numBands = GDALGetRasterCount(ds); + for (i = 1; i <= numBands; i++) { + RASTER_DEBUGF(3, "Processing band %d of %d", i, numBands); + gdband = NULL; + gdband = GDALGetRasterBand(ds, i); + if (NULL == gdband) { + rterror("rt_raster_from_gdal_dataset: Unable to get transformed band\n"); + rt_raster_destroy(rast); + return NULL; + } + + /* pixtype */ + gdpixtype = GDALGetRasterDataType(gdband); + switch (gdpixtype) { + case GDT_Byte: + pt = PT_8BUI; + RASTER_DEBUG(3, "Pixel type is GDT_Byte"); + break; + case GDT_UInt16: + pt = PT_16BUI; + RASTER_DEBUG(3, "Pixel type is GDT_UInt16"); + break; + case GDT_Int16: + pt = PT_16BSI; + RASTER_DEBUG(3, "Pixel type is GDT_Int16"); + break; + case GDT_UInt32: + pt = PT_32BUI; + RASTER_DEBUG(3, "Pixel type is GDT_UInt32"); + break; + case GDT_Int32: + pt = PT_32BSI; + RASTER_DEBUG(3, "Pixel type is GDT_Int32"); + break; + case GDT_Float32: + pt = PT_32BF; + RASTER_DEBUG(3, "Pixel type is GDT_Float32"); + break; + case GDT_Float64: + pt = PT_64BF; + RASTER_DEBUG(3, "Pixel type is GDT_Float64"); + break; + default: + rterror("rt_raster_from_gdal_dataset: Unknown pixel type for transformed raster\n"); + rt_raster_destroy(rast); + return NULL; + break; + } + + /* size: width and height */ + width = GDALGetRasterBandXSize(gdband); + height = GDALGetRasterBandYSize(gdband); + RASTER_DEBUGF(3, "Band dimensions (width x height): %d x %d", width, height); + + /* nodata */ + nodataval = GDALGetRasterNoDataValue(gdband, status); + if (NULL == status || !*status) { + nodataval = 0; + hasnodata = 0; + } + else + hasnodata = 1; + RASTER_DEBUGF(3, "(hasnodata, nodataval) = (%d, %f)", hasnodata, nodataval); + + idx = rt_raster_generate_new_band( + rast, pt, + (hasnodata ? nodataval : 0), + hasnodata, nodataval, rt_raster_get_num_bands(rast) + ); + if (idx < 0) { + rterror("rt_raster_from_gdal_dataset: Could not allocate memory for band\n"); + rt_raster_destroy(rast); + return NULL; + } + band = rt_raster_get_band(rast, idx); + + /* copy data from gdband to band */ + /* try to duplicate entire band data at once */ + do { + values = rtalloc(width * height * sizeof(double)); + if (NULL == values) { + byblock = 1; + break; + } + + cplerr = GDALRasterIO( + gdband, GF_Read, + 0, 0, + width, height, + values, width, height, + GDT_Float64, + 0, 0 + ); + + if (cplerr != CE_None) { + byblock = 1; + rtdealloc(values); + break; + } + + for (y = 0; y < height; y++) { + for (x = 0; x < width; x++) { + value = values[x + y * width]; + + RASTER_DEBUGF(5, "(x, y, value) = (%d, %d, %f)", x, y, value); + + if (rt_band_set_pixel(band, x, y, value) < 0) { + rterror("rt_raster_from_gdal_dataset: Unable to save data from transformed raster\n"); + rt_raster_destroy(rast); + return NULL; + } + } + } + + rtdealloc(values); + } + while (0); + + /* this makes use of GDAL's "natural" blocks */ + if (byblock) { + GDALGetBlockSize(gdband, &nXBlockSize, &nYBlockSize); + nXBlocks = (width + nXBlockSize - 1) / nXBlockSize; + nYBlocks = (height + nYBlockSize - 1) / nYBlockSize; + RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize); + RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks); + + values = rtalloc(nXBlockSize * nYBlockSize * sizeof(double)); + for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) { + for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) { + x = iXBlock * nXBlockSize; + y = iYBlock * nYBlockSize; + + cplerr = GDALRasterIO( + gdband, GF_Read, + x, y, + nXBlockSize, nYBlockSize, + values, nXBlockSize, nYBlockSize, + GDT_Float64, + 0, 0 + ); + if (cplerr != CE_None) { + rterror("rt_raster_from_gdal_dataset: Unable to get data from transformed raster\n"); + rt_raster_destroy(rast); + return NULL; + } + + if ((iXBlock + 1) * nXBlockSize > width) + nXValid = width - (iXBlock * nXBlockSize); + else + nXValid = nXBlockSize; + + if ((iYBlock + 1) * nYBlockSize > height) + nYValid = height - (iYBlock * nYBlockSize); + else + nYValid = nYBlockSize; + + for (iY = 0; iY < nYValid; iY++) { + for (iX = 0; iX < nXValid; iX++) { + x = iX + (nXBlockSize * iXBlock); + y = iY + (nYBlockSize * iYBlock); + value = values[iX + iY * nXBlockSize]; + + RASTER_DEBUGF(5, "(x, y, value) = (%d, %d, %f)", x, y, value); + + if (rt_band_set_pixel(band, x, y, value) < 0) { + rterror("rt_raster_from_gdal_dataset: Unable to save data from transformed raster\n"); + rt_raster_destroy(rast); + return NULL; + } + } + } + } + } + + rtdealloc(values); + } + } + + return rast; +} diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index 9036a2364..3e44f2a0b 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -95,6 +95,7 @@ #include "ogr_api.h" #include "ogr_srs_api.h" #include "cpl_vsi.h" +#include "cpl_conv.h" #include "../../postgis_config.h" /** @@ -892,17 +893,26 @@ typedef struct rt_gdaldriver_t* rt_gdaldriver; rt_gdaldriver rt_raster_gdal_drivers(uint32_t *drv_count); /** - * Return GDAL datasource using GDAL MEM driver from raster + * Return GDAL dataset using GDAL MEM driver from raster * * @param raster : raster to convert to GDAL MEM * @param srs : the raster's coordinate system in OGC WKT or PROJ.4 * @param rtn_drv : is set to the GDAL driver object * - * @return GDAL datasource using GDAL MEM driver + * @return GDAL dataset using GDAL MEM driver */ GDALDatasetH rt_raster_to_gdal_mem(rt_raster raster, char *srs, GDALDriverH *rtn_drv); +/** + * Return a raster from a GDAL dataset + * + * @param ds : the GDAL dataset to convert to a raster + * + * @return raster + */ +rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds); + /*- utilities -------------------------------------------------------*/ /* diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c index 10746447a..e27b5b847 100644 --- a/raster/test/core/testapi.c +++ b/raster/test/core/testapi.c @@ -1282,28 +1282,43 @@ static void testGDALDrivers() { } } -static void testRasterToGDAL(rt_raster raster) { - uint32_t bandNums[] = {1}; - int lenBandNums = 1; +static void testRasterToGDAL() { + rt_raster raster; + rt_band band; + uint32_t x; + uint32_t xmax = 100; + uint32_t y; + uint32_t ymax = 100; + int rtn = 0; + uint64_t gdalSize; - rt_raster rast; - uint8_t *rtn = NULL; + uint8_t *gdal = NULL; - rast = rt_raster_from_band(raster, bandNums, lenBandNums); - assert(rast); + raster = rt_raster_new(xmax, ymax); + assert(raster); /* or we're out of virtual memory */ + band = addBand(raster, PT_32BF, 0, 0); + CHECK(band); + rt_band_set_nodata(band, 0); - rtn = rt_raster_to_gdal(rast, NULL, "GTiff", NULL, &gdalSize); + for (x = 0; x < xmax; x++) { + for (y = 0; y < ymax; y++) { + rtn = rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1)); + CHECK((rtn != -1)); + } + } + + gdal = rt_raster_to_gdal(raster, NULL, "GTiff", NULL, &gdalSize); /*printf("gdalSize: %d\n", (int) gdalSize);*/ CHECK(gdalSize); /* FILE *fh = NULL; fh = fopen("/tmp/out.tif", "w"); - fwrite(rtn, sizeof(uint8_t), gdalSize, fh); + fwrite(gdal, sizeof(uint8_t), gdalSize, fh); fclose(fh); */ - rt_raster_destroy(rast); + deepRelease(raster); } struct rt_valuecount_t { @@ -1369,6 +1384,66 @@ static void testValueCount() { deepRelease(raster); } +static void testGDALToRaster() { + rt_raster raster; + rt_raster rast; + rt_band band; + uint32_t x; + uint32_t xmax = 100; + uint32_t y; + uint32_t ymax = 100; + int rtn = 0; + double value; + + GDALDriverH gddrv = NULL; + GDALDatasetH gdds = NULL; + + raster = rt_raster_new(xmax, ymax); + assert(raster); /* or we're out of virtual memory */ + band = addBand(raster, PT_32BF, 0, 0); + CHECK(band); + rt_band_set_nodata(band, 0); + + for (x = 0; x < xmax; x++) { + for (y = 0; y < ymax; y++) { + rtn = rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1)); + CHECK((rtn != -1)); + } + } + + gdds = rt_raster_to_gdal_mem(raster, NULL, &gddrv); + CHECK(gddrv); + CHECK(gdds); + CHECK((GDALGetRasterXSize(gdds) == xmax)); + CHECK((GDALGetRasterYSize(gdds) == ymax)); + + rast = rt_raster_from_gdal_dataset(gdds); + CHECK(rast); + CHECK((rt_raster_get_num_bands(rast) == 1)); + + band = rt_raster_get_band(rast, 0); + CHECK(band); + + rtn = rt_band_get_pixel(band, 0, 3, &value); + CHECK((rtn != -1)); + CHECK((fabs(value - 0.75) < FLT_EPSILON)); + + rtn = rt_band_get_pixel(band, 99, 0, &value); + CHECK((rtn != -1)); + CHECK((fabs(value - 1.98) < FLT_EPSILON)); + + rtn = rt_band_get_pixel(band, 95, 4, &value); + CHECK((rtn != -1)); + CHECK((fabs(value - 9.54) < FLT_EPSILON)); + + GDALClose(gdds); + GDALDeregisterDriver(gddrv); + GDALDestroyDriver(gddrv); + + deepRelease(rast); + deepRelease(raster); +} + int main() { @@ -1740,7 +1815,7 @@ main() printf("Successfully tested rt_band_reclass\n"); printf("Testing rt_raster_to_gdal\n"); - testRasterToGDAL(raster); + testRasterToGDAL(); printf("Successfully tested rt_raster_to_gdal\n"); printf("Testing rt_raster_gdal_drivers\n"); @@ -1751,6 +1826,10 @@ main() testValueCount(); printf("Successfully tested rt_band_get_value_count\n"); + printf("Testing rt_raster_from_gdal_dataset\n"); + testGDALToRaster(); + printf("Successfully tested rt_raster_from_gdal_dataset\n"); + deepRelease(raster); return EXIT_SUCCESS;