From 2ea489449eeb41e9ffe55e875111ef1da6b0627d Mon Sep 17 00:00:00 2001 From: Bborie Park Date: Mon, 27 Feb 2012 16:42:03 +0000 Subject: [PATCH] Correct handling of 8BSI pixel types when converting a raster to a GDAL MEM dataset. This should resolve the failures in #1617. git-svn-id: http://svn.osgeo.org/postgis/trunk@9313 b70326c6-7e19-0410-871a-916f4a2858ee --- raster/rt_core/rt_api.c | 221 ++++++++++++++++++++++++++++++------- raster/test/core/testapi.c | 145 ++++++++++++++++-------- 2 files changed, 281 insertions(+), 85 deletions(-) diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index 396571e93..53b46d417 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -473,6 +473,16 @@ rt_util_compute_skewed_extent( } rt_raster_set_geotransform_matrix(raster, _gt); + /* get inverse geotransform matrix */ + if (!GDALInvGeoTransform(_gt, _igt)) { + rterror("rt_util_compute_skewed_extent: Unable to compute inverse geotransform matrix"); + rt_raster_destroy(raster); + return 0; + } + RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f", + _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5] + ); + /* shift along axis */ for (i = 0; i < 2; i++) { covers = 0; @@ -489,16 +499,6 @@ rt_util_compute_skewed_extent( return 0; } - /* get inverse geotransform matrix */ - if (!GDALInvGeoTransform(_gt, _igt)) { - rterror("rt_util_compute_skewed_extent: Unable to compute inverse geotransform matrix"); - rt_raster_destroy(raster); - return 0; - } - RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f", - _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5] - ); - /* check the four corners that they are covered along the specific axis pixel column should be >= 0 @@ -539,6 +539,8 @@ rt_util_compute_skewed_extent( return 0; } + RASTER_DEBUGF(4, "Point %d at cell %d x %d", j, (int) _r[0], (int) _r[1]); + /* raster doesn't cover point */ if ((int) _r[i] < 0) { RASTER_DEBUGF(4, "Point outside of skewed extent: %d", j); @@ -590,6 +592,16 @@ rt_util_compute_skewed_extent( RASTER_DEBUGF(4, "Shifted geotransform: %f, %f, %f, %f, %f, %f", _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5] ); + + /* get inverse geotransform matrix */ + if (!GDALInvGeoTransform(_gt, _igt)) { + rterror("rt_util_compute_skewed_extent: Unable to compute inverse geotransform matrix"); + rt_raster_destroy(raster); + return 0; + } + RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f", + _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5] + ); } run++; @@ -614,6 +626,8 @@ rt_util_compute_skewed_extent( return 0; } + RASTER_DEBUGF(4, "geopoint %f x %f at cell %d x %d", extent.MaxX, extent.MinY, (int) _r[0], (int) _r[1]); + rtn = rt_raster_cell_to_geopoint( raster, _r[0] + 1, _r[1] + 1, @@ -631,9 +645,9 @@ rt_util_compute_skewed_extent( RASTER_DEBUGF(4, "Suggested skewed extent: %f, %f, %f, %f", skewedextent->MinX, - skewedextent->MaxY, + skewedextent->MinY, skewedextent->MaxX, - skewedextent->MinY + skewedextent->MaxY ); return 1; @@ -7366,6 +7380,8 @@ rt_raster_to_gdal_mem(rt_raster raster, const char *srs, int i; int numBands; + uint32_t width = 0; + uint32_t height = 0; rt_band rtband = NULL; rt_pixtype pt = PT_END; @@ -7381,8 +7397,13 @@ rt_raster_to_gdal_mem(rt_raster raster, const char *srs, } *rtn_drv = drv; - ds = GDALCreate(drv, "", rt_raster_get_width(raster), - rt_raster_get_height(raster), 0, GDT_Byte, NULL); + width = rt_raster_get_width(raster); + height = rt_raster_get_height(raster); + ds = GDALCreate( + drv, "", + width, height, + 0, GDT_Byte, NULL + ); if (NULL == ds) { rterror("rt_raster_to_gdal_mem: Could not create a GDALDataset to convert into"); return 0; @@ -7446,34 +7467,53 @@ rt_raster_to_gdal_mem(rt_raster raster, const char *srs, if (gdal_pt == GDT_Unknown) rtwarn("rt_raster_to_gdal_mem: Unknown pixel type for band"); - pVoid = rt_band_get_data(rtband); - RASTER_DEBUGF(4, "Band data is at pos %p", pVoid); + /* + For all pixel types other than PT_8BSI, set pointer to start of data + */ + if (pt != PT_8BSI) { + pVoid = rt_band_get_data(rtband); + RASTER_DEBUGF(4, "Band data is at pos %p", pVoid); - pszDataPointer = (char *) rtalloc(20 * sizeof (char)); - sprintf(pszDataPointer, "%p", pVoid); - RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p", - pszDataPointer); + pszDataPointer = (char *) rtalloc(20 * sizeof (char)); + sprintf(pszDataPointer, "%p", pVoid); + RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p", + pszDataPointer); - if (strnicmp(pszDataPointer, "0x", 2) == 0) - sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer); - else - sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer); + if (strnicmp(pszDataPointer, "0x", 2) == 0) + sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer); + else + sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer); - RASTER_DEBUG(3, "Storing info for GDAL MEM raster band"); + RASTER_DEBUG(3, "Storing info for GDAL MEM raster band"); - apszOptions[0] = szGDALOption; - apszOptions[1] = NULL; /* pixel offset, not needed */ - apszOptions[2] = NULL; /* line offset, not needed */ - apszOptions[3] = NULL; + apszOptions[0] = szGDALOption; + apszOptions[1] = NULL; /* pixel offset, not needed */ + apszOptions[2] = NULL; /* line offset, not needed */ + apszOptions[3] = NULL; - /* free */ - rtdealloc(pszDataPointer); + /* free */ + rtdealloc(pszDataPointer); - if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) { - rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band"); - if (allocBandNums) rtdealloc(bandNums); - GDALClose(ds); - return 0; + /* add band */ + if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) { + rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band"); + if (allocBandNums) rtdealloc(bandNums); + GDALClose(ds); + return 0; + } + } + /* + PT_8BSI is special as GDAL has no equivalent pixel type. + Must convert 8BSI to 16BSI so create basic band + */ + else { + /* add band */ + if (GDALAddBand(ds, gdal_pt, NULL) == CE_Failure) { + rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band"); + if (allocBandNums) rtdealloc(bandNums); + GDALClose(ds); + return 0; + } } /* check band count */ @@ -7484,7 +7524,7 @@ rt_raster_to_gdal_mem(rt_raster raster, const char *srs, return 0; } - /* get band to write data to */ + /* get new band */ band = NULL; band = GDALGetRasterBand(ds, i + 1); if (NULL == band) { @@ -7494,6 +7534,99 @@ rt_raster_to_gdal_mem(rt_raster raster, const char *srs, return 0; } + /* PT_8BSI requires manual setting of pixels */ + if (pt == PT_8BSI) { + int nXBlocks, nYBlocks; + int nXBlockSize, nYBlockSize; + int iXBlock, iYBlock; + int nXValid, nYValid; + int iX, iY; + int iXMax, iYMax; + + int x, y, z; + uint32_t valueslen = 0; + int16_t *values = NULL; + double value = 0.; + + /* this makes use of GDAL's "natural" blocks */ + GDALGetBlockSize(band, &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); + + /* length is for the desired pixel type */ + valueslen = rt_pixtype_size(PT_16BSI) * nXBlockSize * nYBlockSize; + values = rtalloc(valueslen); + if (NULL == values) { + rterror("rt_raster_to_gdal_mem: Could not allocate memory for GDAL band pixel values"); + if (allocBandNums) rtdealloc(bandNums); + GDALClose(ds); + return 0; + } + + for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) { + for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) { + memset(values, 0, valueslen); + + x = iXBlock * nXBlockSize; + y = iYBlock * nYBlockSize; + RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock); + RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y); + + /* valid block width */ + if ((iXBlock + 1) * nXBlockSize > width) + nXValid = width - (iXBlock * nXBlockSize); + else + nXValid = nXBlockSize; + + /* valid block height */ + if ((iYBlock + 1) * nYBlockSize > height) + nYValid = height - (iYBlock * nYBlockSize); + else + nYValid = nYBlockSize; + + RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid); + + /* convert 8BSI values to 16BSI */ + z = 0; + iYMax = y + nYValid; + iXMax = x + nXValid; + for (iY = y; iY < iYMax; iY++) { + for (iX = x; iX < iXMax; iX++) { + if (rt_band_get_pixel(rtband, iX, iY, &value) != 0) { + rterror("rt_raster_to_gdal_mem: Could not get pixel value to convert from 8BSI to 16BSI"); + rtdealloc(values); + if (allocBandNums) rtdealloc(bandNums); + GDALClose(ds); + return 0; + } + + values[z++] = rt_util_clamp_to_16BSI(value); + } + } + + /* burn values */ + if (GDALRasterIO( + band, GF_Write, + x, y, + nXValid, nYValid, + values, nXValid, nYValid, + gdal_pt, + 0, 0 + ) != CE_None) { + rterror("rt_raster_to_gdal_mem: Could not write converted 8BSI to 16BSI values to GDAL band"); + rtdealloc(values); + if (allocBandNums) rtdealloc(bandNums); + GDALClose(ds); + return 0; + } + } + } + + rtdealloc(values); + } + /* Add nodata value for band */ if (rt_band_get_hasnodata_flag(rtband) != FALSE) { nodata = rt_band_get_nodata(rtband); @@ -8842,7 +8975,7 @@ rt_raster_gdal_rasterize(const unsigned char *wkb, OGR_G_GetEnvelope(src_geom, &src_env); RASTER_DEBUGF(3, "Suggested extent: %f, %f, %f, %f", - src_env.MinX, src_env.MaxY, src_env.MaxX, src_env.MinY); + src_env.MinX, src_env.MinY, src_env.MaxX, src_env.MaxY); /* user-defined skew */ if (NULL != skew_x) @@ -9016,7 +9149,7 @@ rt_raster_gdal_rasterize(const unsigned char *wkb, #endif RASTER_DEBUGF(3, "Adjusted extent: %f, %f, %f, %f", - src_env.MinX, src_env.MaxY, src_env.MaxX, src_env.MinY); + src_env.MinX, src_env.MinY, src_env.MaxX, src_env.MaxY); } @@ -9054,7 +9187,7 @@ rt_raster_gdal_rasterize(const unsigned char *wkb, } RASTER_DEBUGF(3, "Suggested skewed extent: %f, %f, %f, %f", - skewedextent.MinX, skewedextent.MaxY, skewedextent.MaxX, skewedextent.MinY); + skewedextent.MinX, skewedextent.MinY, skewedextent.MaxX, skewedextent.MaxY); src_env.MinX = skewedextent.MinX; src_env.MinY = skewedextent.MinY; @@ -9069,6 +9202,8 @@ rt_raster_gdal_rasterize(const unsigned char *wkb, ) { ul_user = 1; + RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw); + /* raster dimensions */ if (!_dim[0]) _dim[0] = (int) fmax((fabs(src_env.MaxX - src_env.MinX) + (_scale[0] / 2.)) / _scale[0], 1); @@ -9180,6 +9315,8 @@ rt_raster_gdal_rasterize(const unsigned char *wkb, return NULL; } + RASTER_DEBUGF(4, "Aligning extent to user-specified grid: %f, %f", *grid_xw, *grid_yw); + do { double _d[2] = {0}; double _w[2] = {0}; @@ -9390,6 +9527,7 @@ rt_raster_gdal_rasterize(const unsigned char *wkb, (NULL != scale_x) && (*scale_x < 0.) ) { + RASTER_DEBUG(3, "Processing negative scale-x"); if (!rt_raster_cell_to_geopoint( rast, _dim[0], 0, @@ -9423,6 +9561,7 @@ rt_raster_gdal_rasterize(const unsigned char *wkb, (NULL != scale_y) && (*scale_y > 0) ) { + RASTER_DEBUG(3, "Processing positive scale-y"); if (!rt_raster_cell_to_geopoint( rast, 0, _dim[1], @@ -9457,7 +9596,7 @@ rt_raster_gdal_rasterize(const unsigned char *wkb, } RASTER_DEBUGF(3, "Applied extent: %f, %f, %f, %f", - src_env.MinX, src_env.MaxY, src_env.MaxX, src_env.MinY); + src_env.MinX, src_env.MinY, src_env.MaxX, src_env.MaxY); RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f", dst_gt[0], dst_gt[1], dst_gt[2], dst_gt[3], dst_gt[4], dst_gt[5]); RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d", diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c index 0cdb9f56a..524e070c0 100644 --- a/raster/test/core/testapi.c +++ b/raster/test/core/testapi.c @@ -1733,6 +1733,7 @@ static void testGDALToRaster() { const uint32_t ymax = 100; uint32_t x; uint32_t y; + int v; double values[xmax][ymax]; int rtn = 0; double value; @@ -1779,6 +1780,50 @@ static void testGDALToRaster() { deepRelease(rast); deepRelease(raster); + + raster = rt_raster_new(xmax, ymax); + assert(raster); /* or we're out of virtual memory */ + band = addBand(raster, PT_8BSI, 0, 0); + CHECK(band); + rt_band_set_nodata(band, 0); + + v = -127; + for (x = 0; x < xmax; x++) { + for (y = 0; y < ymax; y++) { + values[x][y] = v++; + rtn = rt_band_set_pixel(band, x, y, values[x][y]); + CHECK((rtn != -1)); + if (v == 128) + v = -127; + } + } + + gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, 0, &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); + CHECK((rt_band_get_pixtype(band) == PT_16BSI)); + + for (x = 0; x < xmax; x++) { + for (y = 0; y < ymax; y++) { + rtn = rt_band_get_pixel(band, x, y, &value); + CHECK((rtn != -1)); + CHECK(FLT_EQ(value, values[x][y])); + } + } + + GDALClose(gdds); + + deepRelease(rast); + deepRelease(raster); } static void testGDALWarp() { @@ -2712,116 +2757,128 @@ main() rt_raster_set_skews(raster, 0, 0); } - printf("Testing rt_raster_gdal_polygonize\n"); + printf("Testing rt_raster_gdal_polygonize... "); testGDALPolygonize(); - printf("Successfully tested rt_raster_gdal_polygonize\n"); + printf("OK\n"); - printf("Testing 1BB band\n"); + printf("Testing 1BB band... "); band_1BB = addBand(raster, PT_1BB, 0, 0); testBand1BB(band_1BB); + printf("OK\n"); - printf("Testing 2BB band\n"); + printf("Testing 2BB band... "); band_2BUI = addBand(raster, PT_2BUI, 0, 0); testBand2BUI(band_2BUI); + printf("OK\n"); - printf("Testing 4BUI band\n"); + printf("Testing 4BUI band... "); band_4BUI = addBand(raster, PT_4BUI, 0, 0); testBand4BUI(band_4BUI); + printf("OK\n"); - printf("Testing 8BUI band\n"); + printf("Testing 8BUI band... "); band_8BUI = addBand(raster, PT_8BUI, 0, 0); testBand8BUI(band_8BUI); + printf("OK\n"); - printf("Testing 8BSI band\n"); + printf("Testing 8BSI band... "); band_8BSI = addBand(raster, PT_8BSI, 0, 0); testBand8BSI(band_8BSI); + printf("OK\n"); - printf("Testing 16BSI band\n"); + printf("Testing 16BSI band... "); band_16BSI = addBand(raster, PT_16BSI, 0, 0); testBand16BSI(band_16BSI); + printf("OK\n"); - printf("Testing 16BUI band\n"); + printf("Testing 16BUI band... "); band_16BUI = addBand(raster, PT_16BUI, 0, 0); testBand16BUI(band_16BUI); + printf("OK\n"); - printf("Testing 32BUI band\n"); + printf("Testing 32BUI band... "); band_32BUI = addBand(raster, PT_32BUI, 0, 0); testBand32BUI(band_32BUI); + printf("OK\n"); - printf("Testing 32BSI band\n"); + printf("Testing 32BSI band... "); band_32BSI = addBand(raster, PT_32BSI, 0, 0); testBand32BSI(band_32BSI); + printf("OK\n"); - printf("Testing 32BF band\n"); + printf("Testing 32BF band... "); band_32BF = addBand(raster, PT_32BF, 0, 0); testBand32BF(band_32BF); + printf("OK\n"); - printf("Testing 64BF band\n"); + printf("Testing 64BF band... "); band_64BF = addBand(raster, PT_64BF, 0, 0); testBand64BF(band_64BF); + printf("OK\n"); - printf("Testing band hasnodata flag\n"); + printf("Testing band hasnodata flag... "); testBandHasNoData(band_64BF); + printf("OK\n"); - printf("Testing rt_raster_from_band\n"); + printf("Testing rt_raster_from_band... "); testRasterFromBand(); - printf("Successfully tested rt_raster_from_band\n"); + printf("OK\n"); - printf("Testing band stats\n"); + printf("Testing band stats... "); testBandStats(); - printf("Successfully tested band stats\n"); + printf("OK\n"); - printf("Testing rt_raster_replace_band\n"); + printf("Testing rt_raster_replace_band... "); testRasterReplaceBand(); - printf("Successfully tested rt_raster_replace_band\n"); + printf("OK\n"); - printf("Testing rt_band_reclass\n"); + printf("Testing rt_band_reclass... "); testBandReclass(); - printf("Successfully tested rt_band_reclass\n"); + printf("OK\n"); - printf("Testing rt_raster_to_gdal\n"); + printf("Testing rt_raster_to_gdal... "); testRasterToGDAL(); - printf("Successfully tested rt_raster_to_gdal\n"); + printf("OK\n"); - printf("Testing rt_raster_gdal_drivers\n"); + printf("Testing rt_raster_gdal_drivers... "); testGDALDrivers(); - printf("Successfully tested rt_raster_gdal_drivers\n"); + printf("OK\n"); - printf("Testing rt_band_get_value_count\n"); + printf("Testing rt_band_get_value_count... "); testValueCount(); - printf("Successfully tested rt_band_get_value_count\n"); + printf("OK\n"); - printf("Testing rt_raster_from_gdal_dataset\n"); + printf("Testing rt_raster_from_gdal_dataset... "); testGDALToRaster(); - printf("Successfully tested rt_raster_from_gdal_dataset\n"); + printf("OK\n"); - printf("Testing rt_util_compute_skewed_extent\n"); + printf("Testing rt_util_compute_skewed_extent... "); testComputeSkewedExtent(); - printf("Successfully tested rt_util_compute_skewed_extent\n"); + printf("OK\n"); - printf("Testing rt_raster_gdal_warp\n"); + printf("Testing rt_raster_gdal_warp... "); testGDALWarp(); - printf("Successfully tested rt_raster_gdal_warp\n"); + printf("OK\n"); - printf("Testing rt_raster_gdal_rasterize\n"); + printf("Testing rt_raster_gdal_rasterize... "); testGDALRasterize(); - printf("Successfully tested rt_raster_gdal_rasterize\n"); + printf("OK\n"); - printf("Testing rt_raster_intersects\n"); + printf("Testing rt_raster_intersects... "); testIntersects(); - printf("Successfully tested rt_raster_intersects\n"); + printf("OK\n"); - printf("Testing rt_raster_same_alignment\n"); + printf("Testing rt_raster_same_alignment... "); testAlignment(); - printf("Successfully tested rt_raster_same_alignment\n"); + printf("OK\n"); - printf("Testing rt_raster_from_two_rasters\n"); + printf("Testing rt_raster_from_two_rasters... "); testFromTwoRasters(); - printf("Successfully tested rt_raster_from_two_rasters\n"); + printf("OK\n"); - printf("Testing rt_raster_load_offline_band\n"); + printf("Testing rt_raster_load_offline_band... "); testLoadOfflineBand(); - printf("Successfully tested rt_raster_load_offline_band\n"); + printf("OK\n"); deepRelease(raster); -- 2.40.0