From 5f7e7742bbe514cc3874493cf033e83e32f7400c Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Fri, 27 Oct 2023 19:12:11 +0000 Subject: [PATCH] code cleanup - getRasterSubset added --- packages/geo/GdalRaster.cpp | 203 +++++++++++--------------- packages/geo/GdalRaster.h | 1 + scripts/systests/subset_perf_test.lua | 4 +- 3 files changed, 88 insertions(+), 120 deletions(-) diff --git a/packages/geo/GdalRaster.cpp b/packages/geo/GdalRaster.cpp index bc663f3a2..6e033f5bb 100644 --- a/packages/geo/GdalRaster.cpp +++ b/packages/geo/GdalRaster.cpp @@ -215,7 +215,6 @@ RasterSubset* GdalRaster::subsetAOI(OGRPolygon* poly) #if SUBSET_DEBUG_TRACE static Mutex mutex; mutex.lock(); - #warning Runing with serialized subsetAOI (very slow!!!) mlog(CRITICAL, "Runing with serialized subsetAOI (very slow!!!)"); #endif @@ -301,8 +300,8 @@ RasterSubset* GdalRaster::subsetAOI(OGRPolygon* poly) map2pixel(aoi_maxx, aoi_miny, lrx, lry); if(SUBSET_DEBUG_TRACE) mlog(DEBUG, "pixel aoi: (%13d, %13d) (%13d, %13d)", ulx, uly, lrx, lry); - int64_t cols2read = lrx - ulx; - int64_t rows2read = lry - uly; + uint32_t _xsize = lrx - ulx; + uint32_t _ysize = lry - uly; /* Sanity check for GCC optimizer 'bug'. Raster's top left corner pixel must be (0, 0) */ int raster_ulx; @@ -321,62 +320,10 @@ RasterSubset* GdalRaster::subsetAOI(OGRPolygon* poly) throw RunTimeException(CRITICAL, RTE_ERROR, "AOI upleft pixel (%d, %d) < raster upleft pixel (%d, %d)", ulx, uly, raster_ulx, raster_uly); } - GDALDataType dtype = band->GetRasterDataType(); - RecordObject::fieldType_t datatype = RecordObject::INVALID_FIELD; - switch (dtype) { - case GDT_Byte: datatype = RecordObject::UINT8; break; - case GDT_Int8: datatype = RecordObject::UINT8; break; - case GDT_UInt16: datatype = RecordObject::UINT16; break; - case GDT_Int16: datatype = RecordObject::INT16; break; - case GDT_UInt32: datatype = RecordObject::UINT32; break; - case GDT_Int32: datatype = RecordObject::INT32; break; - case GDT_UInt64: datatype = RecordObject::UINT64; break; - case GDT_Int64: datatype = RecordObject::INT64; break; - case GDT_Float32: datatype = RecordObject::FLOAT; break; - case GDT_Float64: datatype = RecordObject::DOUBLE; break; - default: throw RunTimeException(CRITICAL, RTE_ERROR, "Unsupported GDT Datatype: %d", dtype); - } - - char* wkt; - targetCRS.exportToWkt(&wkt); - subset = new RasterSubset(cols2read, rows2read, datatype, aoi_minx, aoi_maxy, cellSize, bbox, wkt, gpsTime, fileId); - CPLFree(wkt); - - if(subset->data) - { - int cnt = 1; - OGRErr err = CE_None; - do - { - err = band->RasterIO(GF_Read, ulx, uly, subset->cols, subset->rows, subset->data, subset->cols, subset->rows, dtype, 0, 0, NULL); - } - while(err != CE_None && cnt--); - - if(err != CE_None) - { - ssError |= SS_AOI_FAILED_TO_READ_ERROR; - throw RunTimeException(CRITICAL, RTE_ERROR, "RasterIO call failed: %d", err); - } - - if(SUBSET_DEBUG_TRACE) mlog(DEBUG, "read %ld bytes (%.1fMB), pixel_ulx: %d, pixel_uly: %d, map_ulx: %.2lf, map_uly: %.2lf, cols2read: %ld, rows2read: %ld, datatype %s\n", - subset->size, (float)subset->size/(1024*1024), ulx, uly, subset->map_ulx, subset->map_uly, subset->cols, subset->rows, GDALGetDataTypeName(dtype)); - } - else - { - ssError |= SS_MEMPOOL_ERROR; - uint64_t size = cols2read * rows2read * GDALGetDataTypeSizeBytes(dtype); - mlog(ERROR, "RasterSubset requested memory: %lu MB, available: %lu MB, max: %lu MB", size / (1024*1024), - RasterSubset::poolsize / (1024*1024), - RasterSubset::MAX_SIZE / (1024*1024)); - } + subset = getRasterSubset(ulx, uly, aoi_minx, aoi_maxy, _xsize, _ysize); } catch (const RunTimeException &e) { - if(subset) - { - delete subset; - subset = NULL; - } mlog(e.level(), "Error subsetting: %s", e.what()); } @@ -427,70 +374,14 @@ RasterSubset* GdalRaster::subsetAOI(uint32_t ulx, uint32_t uly, uint32_t _xsize, if(uly + _ysize > ysize) throw RunTimeException(CRITICAL, RTE_ERROR, "rows out of bounds"); + double map_ulx, map_uly; + pixel2map(ulx, uly, map_ulx, map_uly); + // mlog(DEBUG, "upleft pixel: (%13d, %13d) (%.6lf, %.6lf)", ulx, uly, map_ulx, map_uly); - int64_t cols2read = _xsize; - int64_t rows2read = _ysize; - - GDALDataType dtype = band->GetRasterDataType(); - RecordObject::fieldType_t datatype = RecordObject::INVALID_FIELD; - switch (dtype) { - case GDT_Byte: datatype = RecordObject::UINT8; break; - case GDT_Int8: datatype = RecordObject::UINT8; break; - case GDT_UInt16: datatype = RecordObject::UINT16; break; - case GDT_Int16: datatype = RecordObject::INT16; break; - case GDT_UInt32: datatype = RecordObject::UINT32; break; - case GDT_Int32: datatype = RecordObject::INT32; break; - case GDT_UInt64: datatype = RecordObject::UINT64; break; - case GDT_Int64: datatype = RecordObject::INT64; break; - case GDT_Float32: datatype = RecordObject::FLOAT; break; - case GDT_Float64: datatype = RecordObject::DOUBLE; break; - default: throw RunTimeException(CRITICAL, RTE_ERROR, "Unsupported GDT Datatype: %d", dtype); - } - - char* wkt; - targetCRS.exportToWkt(&wkt); - - double mapx, mapy; - pixel2map(ulx, uly, mapx, mapy); - - subset = new RasterSubset(cols2read, rows2read, datatype, mapx, mapy, cellSize, bbox, wkt, gpsTime, fileId); - CPLFree(wkt); - - if(subset->data) - { - int cnt = 1; - OGRErr err = CE_None; - do - { - err = band->RasterIO(GF_Read, ulx, uly, subset->cols, subset->rows, subset->data, subset->cols, subset->rows, dtype, 0, 0, NULL); - } - while(err != CE_None && cnt--); - - if(err != CE_None) - { - ssError |= SS_AOI_FAILED_TO_READ_ERROR; - throw RunTimeException(CRITICAL, RTE_ERROR, "RasterIO call failed: %d", err); - } - - mlog(DEBUG, "read %ld bytes (%.1fMB), pixel_ulx: %d, pixel_uly: %d, map_ulx: %.2lf, map_uly: %.2lf, cols2read: %ld, rows2read: %ld, datatype %s\n", - subset->size, (float)subset->size/(1024*1024), ulx, uly, subset->map_ulx, subset->map_uly, subset->cols, subset->rows, GDALGetDataTypeName(dtype)); - } - else - { - ssError |= SS_MEMPOOL_ERROR; - uint64_t size = cols2read * rows2read * GDALGetDataTypeSizeBytes(dtype); - mlog(ERROR, "RasterSubset requested memory: %lu MB, available: %lu MB, max: %lu MB", size / (1024*1024), - RasterSubset::poolsize / (1024*1024), - RasterSubset::MAX_SIZE / (1024*1024)); - } + subset = getRasterSubset(ulx, uly, map_ulx, map_uly, _xsize, _ysize); } catch (const RunTimeException &e) { - if(subset) - { - delete subset; - subset = NULL; - } mlog(e.level(), "Error subsetting: %s", e.what()); } @@ -1033,6 +924,79 @@ void GdalRaster::readRasterWithRetry(int x, int y, int _xsize, int _ysize, void } +/*---------------------------------------------------------------------------- + * getRasterSubset + *----------------------------------------------------------------------------*/ +RasterSubset* GdalRaster::getRasterSubset(uint32_t ulx, uint32_t uly, double map_ulx, double map_uly, uint32_t _xsize, uint32_t _ysize) +{ + RasterSubset* subset = NULL; + + try + { + GDALDataType dtype = band->GetRasterDataType(); + RecordObject::fieldType_t datatype = RecordObject::INVALID_FIELD; + switch (dtype) { + case GDT_Byte: datatype = RecordObject::UINT8; break; + case GDT_Int8: datatype = RecordObject::UINT8; break; + case GDT_UInt16: datatype = RecordObject::UINT16; break; + case GDT_Int16: datatype = RecordObject::INT16; break; + case GDT_UInt32: datatype = RecordObject::UINT32; break; + case GDT_Int32: datatype = RecordObject::INT32; break; + case GDT_UInt64: datatype = RecordObject::UINT64; break; + case GDT_Int64: datatype = RecordObject::INT64; break; + case GDT_Float32: datatype = RecordObject::FLOAT; break; + case GDT_Float64: datatype = RecordObject::DOUBLE; break; + default: throw RunTimeException(CRITICAL, RTE_ERROR, "Unsupported GDT Datatype: %d", dtype); + } + + char* wkt; + targetCRS.exportToWkt(&wkt); + + subset = new RasterSubset(_xsize, _ysize, datatype, map_ulx, map_uly, cellSize, bbox, wkt, gpsTime, fileId); + CPLFree(wkt); + + if(subset->data) + { + int cnt = 1; + OGRErr err = CE_None; + do + { + err = band->RasterIO(GF_Read, ulx, uly, subset->cols, subset->rows, subset->data, subset->cols, subset->rows, dtype, 0, 0, NULL); + } + while(err != CE_None && cnt--); + + if(err != CE_None) + { + ssError |= SS_AOI_FAILED_TO_READ_ERROR; + throw RunTimeException(CRITICAL, RTE_ERROR, "RasterIO call failed: %d", err); + } + + // mlog(DEBUG, "read %ld bytes (%.1fMB), pixel_ulx: %d, pixel_uly: %d, map_ulx: %.6lf, map_uly: %.6lf, cols2read: %ld, rows2read: %ld, datatype %s\n", + // subset->size, (float)subset->size/(1024*1024), ulx, uly, subset->map_ulx, subset->map_uly, subset->cols, subset->rows, GDALGetDataTypeName(dtype)); + } + else + { + ssError |= SS_MEMPOOL_ERROR; + uint64_t size = _xsize * _ysize * GDALGetDataTypeSizeBytes(dtype); + mlog(ERROR, "RasterSubset requested memory: %lu MB, available: %lu MB, max: %lu MB", size / (1024*1024), + RasterSubset::poolsize / (1024*1024), + RasterSubset::MAX_SIZE / (1024*1024)); + } + } + catch (const RunTimeException &e) + { + if(subset) + { + delete subset; + subset = NULL; + } + mlog(e.level(), "Error subsetting: %s", e.what()); + } + + return subset; +} + + /*---------------------------------------------------------------------------- * maptopixel *----------------------------------------------------------------------------*/ @@ -1048,7 +1012,10 @@ void GdalRaster::map2pixel(double mapx, double mapy, int& x, int& y) *----------------------------------------------------------------------------*/ void GdalRaster::pixel2map(int x, int y, double& mapx, double& mapy) { - mapx = (geoTransform[0] + ((geoTransform[1] * x) + (geoTransform[2] * y))); - mapy = (geoTransform[3] + ((geoTransform[4] * y) + (geoTransform[5] * y))); + double _x = static_cast(x) + 0.5; + double _y = static_cast(y) + 0.5; + + mapx = (geoTransform[0] + ((geoTransform[1] * _x) + (geoTransform[2] * _y))); + mapy = (geoTransform[3] + ((geoTransform[4] * _y) + (geoTransform[5] * _y))); } diff --git a/packages/geo/GdalRaster.h b/packages/geo/GdalRaster.h index a0f55acd0..8f7dedb4a 100644 --- a/packages/geo/GdalRaster.h +++ b/packages/geo/GdalRaster.h @@ -165,6 +165,7 @@ class GdalRaster int radius2pixels (int _radius) const; static inline bool containsWindow (int x, int y, int maxx, int maxy, int windowSize); inline void readRasterWithRetry (int x, int y, int xsize, int ysize, void* data, int dataXsize, int dataYsize, GDALRasterIOExtraArg* args); + RasterSubset* getRasterSubset (uint32_t ulx, uint32_t uly, double map_ulx, double map_uly, uint32_t _xsize, uint32_t _ysize); void map2pixel (double mapx, double mapy, int& x, int& y); void map2pixel (const OGRPoint* poi, int& x, int& y) { map2pixel(poi->getX(), poi->getY(), x, y); } diff --git a/scripts/systests/subset_perf_test.lua b/scripts/systests/subset_perf_test.lua index 38206b425..85fc6a9a8 100644 --- a/scripts/systests/subset_perf_test.lua +++ b/scripts/systests/subset_perf_test.lua @@ -396,8 +396,8 @@ for i, v in ipairs(tbl) do runner.check(cols == 7344) runner.check(rows == 4464) - runner.check(math.abs(ulx - -108.34120000) < sigma) - runner.check(math.abs(uly - 39.19560000) < sigma) + runner.check(math.abs(ulx - gm_llx) < sigma) + runner.check(math.abs(uly - gm_ury) < sigma) runner.check(msg.datatype(datatype) == "UINT8") runner.check(math.abs(cellsize - 0.000083333) < sigma) runner.check(wkt ~= "")