From af981b24bce152926a06814cde17d6beec3d91ac Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Thu, 7 Nov 2024 19:53:39 +0000 Subject: [PATCH 01/17] initial bluetopo, samples only one layer --- CMakeLists.txt | 5 + datasets/bluetopo/CMakeLists.txt | 23 ++ datasets/bluetopo/bluetopo.md | 1 + .../bluetopo/package/BlueTopoBathyRaster.cpp | 229 ++++++++++++++++++ .../bluetopo/package/BlueTopoBathyRaster.h | 98 ++++++++ datasets/bluetopo/package/bluetopo.cpp | 93 +++++++ datasets/bluetopo/package/bluetopo.h | 52 ++++ .../bluetopo/selftests/bluetopo_reader.lua | 56 +++++ packages/geo/GeoIndexedRaster.h | 2 +- platforms/linux/SlideRule.cpp | 12 + .../slideruleearth-aws/asset_directory.csv | 1 + 11 files changed, 571 insertions(+), 1 deletion(-) create mode 100644 datasets/bluetopo/CMakeLists.txt create mode 100644 datasets/bluetopo/bluetopo.md create mode 100644 datasets/bluetopo/package/BlueTopoBathyRaster.cpp create mode 100644 datasets/bluetopo/package/BlueTopoBathyRaster.h create mode 100644 datasets/bluetopo/package/bluetopo.cpp create mode 100644 datasets/bluetopo/package/bluetopo.h create mode 100644 datasets/bluetopo/selftests/bluetopo_reader.lua diff --git a/CMakeLists.txt b/CMakeLists.txt index 1e3cf1e9..cf6d733e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -183,6 +183,7 @@ option (ENABLE_TERMINAL "Instantiate terminal messages" ON) # Dataset Options # option (USE_BATHY_DATASET "Include the BATHY dataset" ON) +option (USE_BLUETOPO_DATASET "Include the BLUETOPO dataset" ON) option (USE_GEBCO_DATASET "Include the GEBCO dataset" ON) option (USE_GEDI_DATASET "Include the GEDI dataset" ON) option (USE_ICESAT2_DATASET "Include the ICESat2 dataset" ON) @@ -364,6 +365,10 @@ if(${USE_BATHY_DATASET}) add_subdirectory (datasets/bathy) endif() +if(${USE_BLUETOPO_DATASET}) + add_subdirectory (datasets/bluetopo) +endif() + if(${USE_GEBCO_DATASET}) add_subdirectory (datasets/gebco) endif() diff --git a/datasets/bluetopo/CMakeLists.txt b/datasets/bluetopo/CMakeLists.txt new file mode 100644 index 00000000..3c5186ab --- /dev/null +++ b/datasets/bluetopo/CMakeLists.txt @@ -0,0 +1,23 @@ +message (STATUS "Including bluetopo dataset") + +target_compile_definitions (slideruleLib PUBLIC __bluetopo__) + +target_sources(slideruleLib + PRIVATE + ${CMAKE_CURRENT_LIST_DIR}/package/bluetopo.cpp + ${CMAKE_CURRENT_LIST_DIR}/package/BlueTopoBathyRaster.cpp +) + +target_include_directories (slideruleLib + PUBLIC + $ + $/package +) + +install ( + FILES + ${CMAKE_CURRENT_LIST_DIR}/package/bluetopo.h + ${CMAKE_CURRENT_LIST_DIR}/package/BlueTopoBathyRaster.h + DESTINATION + ${INCDIR}/bluetopo +) diff --git a/datasets/bluetopo/bluetopo.md b/datasets/bluetopo/bluetopo.md new file mode 100644 index 00000000..10347982 --- /dev/null +++ b/datasets/bluetopo/bluetopo.md @@ -0,0 +1 @@ +Raster sampling support for BlueTopo (The compilation of the nation's best available bathymetric data.) \ No newline at end of file diff --git a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp new file mode 100644 index 00000000..4e493b42 --- /dev/null +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp @@ -0,0 +1,229 @@ +/* + * Copyright (c) 2021, University of Washington + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * 3. Neither the name of the University of Washington nor the names of its + * contributors may be used to endorse or promote products derived from this + * software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF WASHINGTON AND CONTRIBUTORS + * “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF WASHINGTON OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/****************************************************************************** + * INCLUDES + ******************************************************************************/ + +#include "BlueTopoBathyRaster.h" + +/****************************************************************************** + * PROTECTED METHODS + ******************************************************************************/ + +/*---------------------------------------------------------------------------- + * Constructor + *----------------------------------------------------------------------------*/ +BlueTopoBathyRaster::BlueTopoBathyRaster(lua_State* L, RequestFields* rqst_parms, const char* key): + GeoIndexedRaster(L, rqst_parms, key), + filePath(parms->asset.asset->getPath()), + indexBucket(parms->asset.asset->getIndex()) +{ + std::string bucketPath = filePath + indexBucket; + if(!findIndexFileInS3Bucket(bucketPath)) + { + throw RunTimeException(CRITICAL, RTE_ERROR, "Failed to create BlueTopoBathyRaster"); + } +} + +/*---------------------------------------------------------------------------- + * Destructor + *----------------------------------------------------------------------------*/ +BlueTopoBathyRaster::~BlueTopoBathyRaster(void) = default; + +/*---------------------------------------------------------------------------- + * getIndexFile + *----------------------------------------------------------------------------*/ +void BlueTopoBathyRaster::getIndexFile(const OGRGeometry* geo, std::string& file, const std::vector* points) +{ + static_cast(geo); + static_cast(points); + file = indexFile; + mlog(DEBUG, "Using %s", file.c_str()); +} + + +/*---------------------------------------------------------------------------- + * findRasters + *----------------------------------------------------------------------------*/ +bool BlueTopoBathyRaster::findRasters(raster_finder_t* finder) +{ + const std::vector* flist = finder->featuresList; + const OGRGeometry* geo = finder->geo; + const uint32_t start = 0; + const uint32_t end = flist->size(); + + try + { + const std::string token(".amazonaws.com/BlueTopo/"); + + for(uint32_t i = start; i < end; i++) + { + OGRFeature* feature = flist->at(i); + OGRGeometry *rastergeo = feature->GetGeometryRef(); + + if (!rastergeo->Intersects(geo)) continue; + + rasters_group_t* rgroup = new rasters_group_t; + rgroup->gpsTime = getGmtDate(feature, "Delivered_Date", rgroup->gmtDate); + + const char* dataFile = feature->GetFieldAsString("GeoTIFF_link"); + if(dataFile && (strlen(dataFile) > 0)) + { + std::string fullPath(dataFile); + size_t pos = fullPath.find(token); + if(pos != std::string::npos) + { + pos += token.length(); + + std::string rasterName = fullPath.substr(pos); + fullPath = filePath + rasterName; + } + else + { + mlog(WARNING, "Could not find token %s in %s", token.c_str(), dataFile); + continue; + } + + raster_info_t rinfo; + rinfo.dataIsElevation = true; + rinfo.tag = VALUE_TAG; + rinfo.fileId = finder->fileDict.add(fullPath); + rgroup->infovect.push_back(rinfo); + } + rgroup->infovect.shrink_to_fit(); + + mlog(DEBUG, "Added group with %ld rasters", rgroup->infovect.size()); + for(unsigned j = 0; j < rgroup->infovect.size(); j++) + mlog(DEBUG, " %s", finder->fileDict.get(rgroup->infovect[j].fileId)); + + // Add the group + finder->rasterGroups.push_back(rgroup); + } + finder->rasterGroups.shrink_to_fit(); + + mlog(DEBUG, "Found %ld raster groups", finder->rasterGroups.size()); + } + catch (const RunTimeException &e) + { + mlog(e.level(), "Error getting time from raster feature file: %s", e.what()); + } + + return (!finder->rasterGroups.empty()); +} + + +/*---------------------------------------------------------------------------- + * getGmtDate + *----------------------------------------------------------------------------*/ +double BlueTopoBathyRaster::getGmtDate(const OGRFeature* feature, const char* field, TimeLib::gmt_time_t& gmtDate) +{ + const int i = feature->GetFieldIndex(field); + if(i == -1) + { + mlog(ERROR, "Time field: %s not found, unable to get GMT date", field); + return 0; + } + + double gpstime = 0; + double seconds; + int year; + int month; + int day; + int hour; + int minute; + + /* + * The raster's 'Delivered_Date' in the GPKG index file is not in ISO8601 format. + * Instead, it uses the format "YYYY-MM-DD HH:MM:SS". + */ + if(const char* dateStr = feature->GetFieldAsString(i)) + { + if(sscanf(dateStr, "%04d-%02d-%02d %02d:%02d:%lf", + &year, &month, &day, &hour, &minute, &seconds) == 6) + { + gmtDate.year = year; + gmtDate.doy = TimeLib::dayofyear(year, month, day); + gmtDate.hour = hour; + gmtDate.minute = minute; + gmtDate.second = seconds; + gmtDate.millisecond = 0; + gpstime = TimeLib::gmt2gpstime(gmtDate); + // mlog(DEBUG, "%04d:%02d:%02d:%02d:%02d:%02d", year, month, day, hour, minute, (int)seconds); + } + else mlog(DEBUG, "Unable to parse date string [%s]", dateStr); + } + else mlog(DEBUG, "Date field is invalid"); + + return gpstime; +} + + +/****************************************************************************** + * PRIVATE METHODS + ******************************************************************************/ + +/*---------------------------------------------------------------------------- + * findIndexFileInS3Bucket + *----------------------------------------------------------------------------*/ +bool BlueTopoBathyRaster::findIndexFileInS3Bucket(const std::string& bucketPath) +{ + bool found = false; + char** fileList = VSIReadDir(bucketPath.c_str()); + + if(fileList) + { + /* Look for .gpkg file, assume there is only one in the bucket */ + std::string token = ".gpkg"; + + for(int i = 0; fileList[i] != nullptr; i++) + { + std::string fileName = fileList[i]; + + if(fileName.size() >= token.size() && + fileName.compare(fileName.size() - token.size(), token.size(), token) == 0) + { + indexFile = bucketPath + "/" + fileName; + found = true; + mlog(DEBUG, "Found index file: %s", indexFile.c_str()); + break; + } + } + CSLDestroy(fileList); + } + + if(!found) + { + mlog(CRITICAL, "Failed to find index file in bucket: %s", bucketPath.c_str()); + } + + return found; +} diff --git a/datasets/bluetopo/package/BlueTopoBathyRaster.h b/datasets/bluetopo/package/BlueTopoBathyRaster.h new file mode 100644 index 00000000..d25bf8e3 --- /dev/null +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.h @@ -0,0 +1,98 @@ +/* + * Copyright (c) 2021, University of Washington + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * 3. Neither the name of the University of Washington nor the names of its + * contributors may be used to endorse or promote products derived from this + * software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF WASHINGTON AND CONTRIBUTORS + * “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF WASHINGTON OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef __bluetopo_bathy_raster__ +#define __bluetopo_bathy_raster__ + +/****************************************************************************** + * INCLUDES + ******************************************************************************/ + +#include "GeoIndexedRaster.h" + +/****************************************************************************** + * GEBCO BATHY RASTER CLASS + ******************************************************************************/ + +class BlueTopoBathyRaster: public GeoIndexedRaster +{ + public: + + /*-------------------------------------------------------------------- + * Constants + *--------------------------------------------------------------------*/ + + /*-------------------------------------------------------------------- + * Typedefs + *--------------------------------------------------------------------*/ + + /*-------------------------------------------------------------------- + * Methods + *--------------------------------------------------------------------*/ + + static RasterObject* create(lua_State* L, RequestFields* rqst_parms, const char* key) + { return new BlueTopoBathyRaster(L, rqst_parms, key); } + + + protected: + + /*-------------------------------------------------------------------- + * Methods + *--------------------------------------------------------------------*/ + double getGmtDate (const OGRFeature* feature, const char* field, TimeLib::gmt_time_t& gmtDate) final; + + BlueTopoBathyRaster (lua_State* L, RequestFields* rqst_parms, const char* key); + ~BlueTopoBathyRaster (void) override; + + void getIndexFile (const OGRGeometry* geo, std::string& file, const std::vector* points) final; + bool findRasters (raster_finder_t* finder) final; + + /*-------------------------------------------------------------------- + * Data + *--------------------------------------------------------------------*/ + + private: + + /*-------------------------------------------------------------------- + * Methods + *--------------------------------------------------------------------*/ + + bool findIndexFileInS3Bucket(const std::string& bucketPath); + + /*-------------------------------------------------------------------- + * Data + *--------------------------------------------------------------------*/ + std::string filePath; + std::string indexBucket; + std::string indexFile; +}; + +#endif /* __gebco_bathy_raster__ */ diff --git a/datasets/bluetopo/package/bluetopo.cpp b/datasets/bluetopo/package/bluetopo.cpp new file mode 100644 index 00000000..52a6b50c --- /dev/null +++ b/datasets/bluetopo/package/bluetopo.cpp @@ -0,0 +1,93 @@ +/* + * Copyright (c) 2021, University of Washington + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * 3. Neither the name of the University of Washington nor the names of its + * contributors may be used to endorse or promote products derived from this + * software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF WASHINGTON AND CONTRIBUTORS + * “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF WASHINGTON OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/****************************************************************************** + *INCLUDES + ******************************************************************************/ + +#include "OsApi.h" +#include "bluetopo.h" + +/****************************************************************************** + * DEFINES + ******************************************************************************/ + +#define LUA_BLUETOPO_LIBNAME "bluetopo" +#define LUA_BLUETOPO_RASTER_NAME "bluetopo-bathy" + +/****************************************************************************** + * LOCAL FUNCTIONS + ******************************************************************************/ + +/*---------------------------------------------------------------------------- + * bluetopo_open + *----------------------------------------------------------------------------*/ +int bluetopo_open (lua_State *L) +{ + static const struct luaL_Reg bluetopo_functions[] = { + {NULL, NULL} + }; + + /* Set Library */ + luaL_newlib(L, bluetopo_functions); + + return 1; +} + +/****************************************************************************** + * EXPORTED FUNCTIONS + ******************************************************************************/ + +extern "C" { +void initbluetopo(void) +{ + /* Initialize Modules */ + BlueTopoBathyRaster::init(); + + /* Register Rasters */ + RasterObject::registerRaster(LUA_BLUETOPO_RASTER_NAME, BlueTopoBathyRaster::create); + + /* Extend Lua */ + LuaEngine::extend(LUA_BLUETOPO_LIBNAME, bluetopo_open); + + /* Indicate Presence of Package */ + LuaEngine::indicate(LUA_BLUETOPO_LIBNAME, LIBID); + + /* Display Status */ + print2term("%s package initialized (%s)\n", LUA_BLUETOPO_LIBNAME, LIBID); +} + +void deinitbluetopo(void) +{ + /* Uninitialize Modules */ + BlueTopoBathyRaster::deinit(); +} +} diff --git a/datasets/bluetopo/package/bluetopo.h b/datasets/bluetopo/package/bluetopo.h new file mode 100644 index 00000000..4073d494 --- /dev/null +++ b/datasets/bluetopo/package/bluetopo.h @@ -0,0 +1,52 @@ +/* + * Copyright (c) 2021, University of Washington + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * 3. Neither the name of the University of Washington nor the names of its + * contributors may be used to endorse or promote products derived from this + * software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF WASHINGTON AND CONTRIBUTORS + * “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF WASHINGTON OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef __bluetopo_plugin__ +#define __bluetopo_plugin__ + +/****************************************************************************** + * INCLUDES + ******************************************************************************/ + +#include "BlueTopoBathyRaster.h" + +/****************************************************************************** + * PROTOTYPES + ******************************************************************************/ + +extern "C" { +void initbluetopo(void); +void deinitbluetopo(void); +} + +#endif /* __bluetopo_plugin__ */ + + diff --git a/datasets/bluetopo/selftests/bluetopo_reader.lua b/datasets/bluetopo/selftests/bluetopo_reader.lua new file mode 100644 index 00000000..8422c0c0 --- /dev/null +++ b/datasets/bluetopo/selftests/bluetopo_reader.lua @@ -0,0 +1,56 @@ +local runner = require("test_executive") +local console = require("console") +local asset = require("asset") +local _,td = runner.srcscript() + +-- console.monitor:config(core.LOG, core.DEBUG) +-- sys.setlvl(core.LOG, core.DEBUG) + +-- Setup -- +local assets = asset.loaddir() + +-- Correct values test for different POIs + +local lons = {-80.87, -81.02, -89.66, -94.72} +local lats = { 32.06, 31.86, 29.99, 29.35} +height = 0 + +local expDepth = { -2.02, -14.03, -4.28, -17.18} +local expFlags = { 3.32} + +-- Tolerance for depth values is needed because of the different versions of PROJLIB +-- The values are different in the different versions of the library. +-- We add 1 meter of depth tolerance for the test to pass +local depth_tolerance = 1; + +print(string.format("\n--------------------------------\nTest: BlueTopo Correct Values\n--------------------------------")) +local dem = geo.raster(geo.parms({ asset = "bluetopo-bathy", algorithm = "NearestNeighbour", with_flags=true, sort_by_index = true })) +runner.check(dem ~= nil) + +for j, lon in ipairs(lons) do + lat = lats[j] + tbl, err = dem:sample(lon, lat, height) + runner.check(err == 0) + runner.check(tbl ~= nil) + + if err ~= 0 then + print(string.format("Point: %d, (%.3f, %.3f) ======> FAILED to read, err# %d",j, lon, lat, err)) + else + local el, fname + for k, v in ipairs(tbl) do + value = v["value"] + fname = v["file"] + flags = v["flags"] + print(string.format("(%02d) (%6.2f, %6.2f) %6.2fm %4.2f %s", k, lon, lat, value, flags, fname)) + + assert( math.abs(value - expDepth[j]) < depth_tolerance, string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) + -- assert(flags == expFlags[j], string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) + end + end +end + +-- Report Results -- + +runner.report() + + diff --git a/packages/geo/GeoIndexedRaster.h b/packages/geo/GeoIndexedRaster.h index cfeae801..158fbc53 100644 --- a/packages/geo/GeoIndexedRaster.h +++ b/packages/geo/GeoIndexedRaster.h @@ -228,7 +228,7 @@ class GeoIndexedRaster: public RasterObject virtual void getGroupSubsets (const rasters_group_t* rgroup, List& slist); uint32_t getGroupFlags (const rasters_group_t* rgroup); - static double getGmtDate (const OGRFeature* feature, const char* field, TimeLib::gmt_time_t& gmtDate); + virtual double getGmtDate (const OGRFeature* feature, const char* field, TimeLib::gmt_time_t& gmtDate); bool openGeoIndex (const OGRGeometry* geo, const std::vector* points); virtual bool getFeatureDate (const OGRFeature* feature, TimeLib::gmt_time_t& gmtDate); virtual void getIndexFile (const OGRGeometry* geo, std::string& file, const std::vector* points=NULL) = 0; diff --git a/platforms/linux/SlideRule.cpp b/platforms/linux/SlideRule.cpp index e1ac0a8b..99d8638c 100644 --- a/platforms/linux/SlideRule.cpp +++ b/platforms/linux/SlideRule.cpp @@ -75,6 +75,10 @@ #include "bathy.h" #endif +#ifdef __bluetopo__ +#include "bluetopo.h" +#endif + #ifdef __gebco__ #include "gebco.h" #endif @@ -368,6 +372,10 @@ int main (int argc, char* argv[]) initbathy(); #endif + #ifdef __bluetopo__ + initbluetopo(); + #endif + #ifdef __gebco__ initgebco(); #endif @@ -441,6 +449,10 @@ int main (int argc, char* argv[]) deinitbathy(); #endif + #ifdef __bluetopo__ + deinitbluetopo(); + #endif + #ifdef __gebco__ deinitgebco(); #endif diff --git a/targets/slideruleearth-aws/asset_directory.csv b/targets/slideruleearth-aws/asset_directory.csv index 3bae45c4..5c99bc51 100644 --- a/targets/slideruleearth-aws/asset_directory.csv +++ b/targets/slideruleearth-aws/asset_directory.csv @@ -20,6 +20,7 @@ usgs3dep-1meter-dem, nil, nil, /vsis3/prd-tnm, esa-worldcover-10meter, nil, nil, /vsis3/esa-worldcover/v200/2021/map, /vsis3/sliderule/data/WORLDCOVER/ESA_WorldCover_10m_2021_v200_Map.vrt, eu-central-1, https://s3.eu-central-1.amazonaws.com meta-globalcanopy-1meter, nil, nil, /vsis3/dataforgood-fb-data/forests/v1/alsgedi_global_v6_float/, /vsis3/sliderule/data/GLOBALCANOPY/META_GlobalCanopyHeight_1m_2024_v1.vrt, us-east-1, https://s3.us-east-1.amazonaws.com gebco-bathy, iam-role, s3, sliderule/data/GEBCO/2023, index.geojson, us-west-2, https://s3.us-west-2.amazonaws.com +bluetopo-bathy, iam-role, s3, /vsis3/noaa-ocs-nationalbathymetry-pds/BlueTopo/, _BlueTopo_Tile_Scheme, us-east-1, https://s3.us-east-1.amazonaws.com landsat-hls, lpdaac-cloud, nil, /vsis3/lp-prod-protected, nil, us-west-2, https://s3.us-west-2.amazonaws.com arcticdem-mosaic, nil, nil, /vsis3/pgc-opendata-dems/arcticdem/mosaics/v4.1, 2m_dem_tiles.vrt, us-west-2, https://s3.us-west-2.amazonaws.com arcticdem-strips, nil, nil, /vsis3/pgc-opendata-dems/arcticdem/strips/s2s041/2m, nil, us-west-2, https://s3.us-west-2.amazonaws.com From 5412497df06aab89fcee6d5e75a86ef3e25b1a00 Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Mon, 11 Nov 2024 17:05:19 +0000 Subject: [PATCH 02/17] added bandName to RasterSample, fixed landsat bands --- .../bluetopo/package/BlueTopoBathyRaster.cpp | 12 +++-- datasets/landsat/package/LandsatHlsRaster.cpp | 13 +++++ packages/arrow/ArrowSamplerImpl.cpp | 48 +++++++++++++++++-- packages/geo/GeoIndexedRaster.cpp | 10 ++-- packages/geo/GeoIndexedRaster.h | 6 +-- packages/geo/RasterObject.h | 5 ++ packages/geo/RasterSample.h | 2 + scripts/selftests/test_runner.lua | 10 +++- 8 files changed, 86 insertions(+), 20 deletions(-) diff --git a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp index 4e493b42..73b701b9 100644 --- a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp @@ -47,7 +47,7 @@ BlueTopoBathyRaster::BlueTopoBathyRaster(lua_State* L, RequestFields* rqst_parms filePath(parms->asset.asset->getPath()), indexBucket(parms->asset.asset->getIndex()) { - std::string bucketPath = filePath + indexBucket; + const std::string bucketPath = filePath + indexBucket; if(!findIndexFileInS3Bucket(bucketPath)) { throw RunTimeException(CRITICAL, RTE_ERROR, "Failed to create BlueTopoBathyRaster"); @@ -104,12 +104,13 @@ bool BlueTopoBathyRaster::findRasters(raster_finder_t* finder) { pos += token.length(); - std::string rasterName = fullPath.substr(pos); + const std::string rasterName = fullPath.substr(pos); fullPath = filePath + rasterName; } else { mlog(WARNING, "Could not find token %s in %s", token.c_str(), dataFile); + delete rgroup; continue; } @@ -202,16 +203,17 @@ bool BlueTopoBathyRaster::findIndexFileInS3Bucket(const std::string& bucketPath) if(fileList) { /* Look for .gpkg file, assume there is only one in the bucket */ - std::string token = ".gpkg"; + const std::string token = ".gpkg"; for(int i = 0; fileList[i] != nullptr; i++) { - std::string fileName = fileList[i]; + const std::string fileName = fileList[i]; if(fileName.size() >= token.size() && fileName.compare(fileName.size() - token.size(), token.size(), token) == 0) { - indexFile = bucketPath + "/" + fileName; + indexFile = bucketPath; + indexFile.append("/").append(fileName); found = true; mlog(DEBUG, "Found index file: %s", indexFile.c_str()); break; diff --git a/datasets/landsat/package/LandsatHlsRaster.cpp b/datasets/landsat/package/LandsatHlsRaster.cpp index c6284707..73875528 100644 --- a/datasets/landsat/package/LandsatHlsRaster.cpp +++ b/datasets/landsat/package/LandsatHlsRaster.cpp @@ -313,6 +313,7 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr const bool returnBandSample = it->second; if(returnBandSample) { + sample->bandName = bandName; // TODO: Remove when sample contains band name slist->add(sample); item->sample = NULL; } @@ -372,6 +373,9 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr s = new RasterSample(*ps.sample); } + /* Set band name for this sample */ + s->bandName = bandName; //TODO: remove when sample contains band name + /* Set time for this sample */ s->time = rgroup->gpsTime / 1000; @@ -417,6 +421,9 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr if((green != invalid) && (swir16 != invalid)) sample->value = (green - swir16) / (green + swir16); else sample->value = invalid; + + /* Set band name for this sample */ + sample->bandName = "NDSI"; slist->add(sample); } @@ -426,6 +433,9 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr if((red != invalid) && (nir08 != invalid)) sample->value = (nir08 - red) / (nir08 + red); else sample->value = invalid; + + /* Set band name for this sample */ + sample->bandName = "NDVI"; slist->add(sample); } @@ -435,6 +445,9 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr if((nir08 != invalid) && (swir16 != invalid)) sample->value = (nir08 - swir16) / (nir08 + swir16); else sample->value = invalid; + + /* Set band name for this sample */ + sample->bandName = "NDWI"; slist->add(sample); } diff --git a/packages/arrow/ArrowSamplerImpl.cpp b/packages/arrow/ArrowSamplerImpl.cpp index 11f47789..8c38929c 100644 --- a/packages/arrow/ArrowSamplerImpl.cpp +++ b/packages/arrow/ArrowSamplerImpl.cpp @@ -434,6 +434,9 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl const RasterObject* robj = sampler->robj; /* Create list builders for the new columns */ + arrow::ListBuilder band_list_builder(pool, std::make_shared()); + auto* band_builder = dynamic_cast(band_list_builder.value_builder()); + arrow::ListBuilder value_list_builder(pool, std::make_shared()); auto* value_builder = dynamic_cast(value_list_builder.value_builder()); @@ -468,7 +471,7 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl arrow::ListBuilder mad_list_builder(pool, std::make_shared()); auto* mad_builder = dynamic_cast(mad_list_builder.value_builder()); - std::shared_ptr value_list_array, time_list_array, fileid_list_array, flags_list_array; + std::shared_ptr band_list_array, value_list_array, time_list_array, fileid_list_array, flags_list_array; std::shared_ptr count_list_array, min_list_array, max_list_array, mean_list_array, median_list_array, stdev_list_array, mad_list_array; /* Iterate over each sample in a vector of lists of samples */ @@ -477,6 +480,10 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl sample_list_t* slist = sampler->samples[i]; /* Start new lists */ + if(robj->hasBands()) + { + PARQUET_THROW_NOT_OK(band_list_builder.Append()); + } PARQUET_THROW_NOT_OK(value_list_builder.Append()); PARQUET_THROW_NOT_OK(time_list_builder.Append()); if(robj->hasFlags()) @@ -503,7 +510,10 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl { const RasterSample* sample = slist->get(j); - /* Append the value to the value list */ + if(robj->hasBands()) + { + PARQUET_THROW_NOT_OK(band_builder->Append(sample->bandName)); + } PARQUET_THROW_NOT_OK(value_builder->Append(sample->value)); PARQUET_THROW_NOT_OK(time_builder->Append(sample->time)); if(robj->hasFlags()) @@ -525,6 +535,10 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl } /* Finish the list builders */ + if(robj->hasBands()) + { + PARQUET_THROW_NOT_OK(band_list_builder.Finish(&band_list_array)); + } PARQUET_THROW_NOT_OK(value_list_builder.Finish(&value_list_array)); PARQUET_THROW_NOT_OK(time_list_builder.Finish(&time_list_array)); if(robj->hasFlags()) @@ -546,6 +560,7 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl const std::string prefix = sampler->rkey; /* Create fields for the new columns */ + auto band_field = std::make_shared(prefix + ".band", arrow::list(arrow::utf8())); auto value_field = std::make_shared(prefix + ".value", arrow::list(arrow::float64())); auto time_field = std::make_shared(prefix + ".time", arrow::list(arrow::float64())); auto flags_field = std::make_shared(prefix + ".flags", arrow::list(arrow::uint32())); @@ -565,6 +580,10 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl mutex.lock(); { /* Add new columns fields */ + if(robj->hasBands()) + { + newFields.push_back(band_field); + } newFields.push_back(value_field); newFields.push_back(time_field); if(robj->hasFlags()) @@ -584,6 +603,10 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl } /* Add new columns */ + if(robj->hasBands()) + { + newColumns.push_back(std::make_shared(band_list_array)); + } newColumns.push_back(std::make_shared(value_list_array)); newColumns.push_back(std::make_shared(time_list_array)); if(robj->hasFlags()) @@ -616,6 +639,7 @@ bool ArrowSamplerImpl::makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* s const RasterObject* robj = sampler->robj; /* Create builders for the new columns */ + arrow::StringBuilder band_builder(pool); arrow::DoubleBuilder value_builder(pool); arrow::DoubleBuilder time_builder(pool); arrow::UInt32Builder flags_builder(pool); @@ -630,7 +654,7 @@ bool ArrowSamplerImpl::makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* s arrow::DoubleBuilder stdev_builder(pool); arrow::DoubleBuilder mad_builder(pool); - std::shared_ptr value_array, time_array, fileid_array, flags_array; + std::shared_ptr band_array, value_array, time_array, fileid_array, flags_array; std::shared_ptr count_array, min_array, max_array, mean_array, median_array, stdev_array, mad_array; RasterSample fakeSample(0.0, 0); @@ -652,7 +676,10 @@ bool ArrowSamplerImpl::makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* s */ sample = &fakeSample; } - + if(robj->hasBands()) + { + PARQUET_THROW_NOT_OK(band_builder.Append(sample->bandName)); + } PARQUET_THROW_NOT_OK(value_builder.Append(sample->value)); PARQUET_THROW_NOT_OK(time_builder.Append(sample->time)); if(robj->hasFlags()) @@ -673,6 +700,10 @@ bool ArrowSamplerImpl::makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* s } /* Finish the builders */ + if(robj->hasBands()) + { + PARQUET_THROW_NOT_OK(band_builder.Finish(&band_array)); + } PARQUET_THROW_NOT_OK(value_builder.Finish(&value_array)); PARQUET_THROW_NOT_OK(time_builder.Finish(&time_array)); if(robj->hasFlags()) @@ -694,6 +725,7 @@ bool ArrowSamplerImpl::makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* s const std::string prefix = sampler->rkey; /* Create fields for the new columns */ + auto band_field = std::make_shared(prefix + ".band", arrow::utf8()); auto value_field = std::make_shared(prefix + ".value", arrow::float64()); auto time_field = std::make_shared(prefix + ".time", arrow::float64()); auto flags_field = std::make_shared(prefix + ".flags", arrow::uint32()); @@ -713,6 +745,10 @@ bool ArrowSamplerImpl::makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* s mutex.lock(); { /* Add new columns fields */ + if(robj->hasBands()) + { + newFields.push_back(band_field); + } newFields.push_back(value_field); newFields.push_back(time_field); if(robj->hasFlags()) @@ -732,6 +768,10 @@ bool ArrowSamplerImpl::makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* s } /* Add new columns */ + if(robj->hasBands()) + { + newColumns.push_back(std::make_shared(band_array)); + } newColumns.push_back(std::make_shared(value_array)); newColumns.push_back(std::make_shared(time_array)); if(robj->hasFlags()) diff --git a/packages/geo/GeoIndexedRaster.cpp b/packages/geo/GeoIndexedRaster.cpp index 8f0e6570..ed9b2139 100644 --- a/packages/geo/GeoIndexedRaster.cpp +++ b/packages/geo/GeoIndexedRaster.cpp @@ -958,10 +958,10 @@ void* GeoIndexedRaster::batchReaderThread(void *param) { unique_raster_t* ur = breader->uraster; GdalRaster* raster = new GdalRaster(breader->obj->parms, - breader->obj->fileDict.get(ur->fileId), + breader->obj->fileDict.get(ur->rinfo->fileId), 0, /* Sample collecting code will set it to group's gpsTime */ - ur->fileId, - ur->dataIsElevation, + ur->rinfo->fileId, + ur->rinfo->dataIsElevation, breader->obj->crscb); /* Sample all points for this raster */ @@ -1569,7 +1569,7 @@ bool GeoIndexedRaster::findUniqueRasters(std::vector& uniqueRa else { /* Raster is not in the vector of unique rasters */ - unique_raster_t* ur = new unique_raster_t(rinfo.dataIsElevation, rinfo.fileId); + unique_raster_t* ur = new unique_raster_t(&rinfo); uniqueRasters.push_back(ur); /* Set pointer in rinfo to new unique raster */ @@ -1586,7 +1586,7 @@ bool GeoIndexedRaster::findUniqueRasters(std::vector& uniqueRa mlog(DEBUG, "Finding points for unique rasters"); for(unique_raster_t* ur : uniqueRasters) { - auto it = rasterToPointsMap.find(fileDict.get(ur->fileId)); + auto it = rasterToPointsMap.find(fileDict.get(ur->rinfo->fileId)); if(it != rasterToPointsMap.end()) { for(const uint32_t pointIndx : it->second) diff --git a/packages/geo/GeoIndexedRaster.h b/packages/geo/GeoIndexedRaster.h index 158fbc53..e1ecfe1b 100644 --- a/packages/geo/GeoIndexedRaster.h +++ b/packages/geo/GeoIndexedRaster.h @@ -114,11 +114,9 @@ class GeoIndexedRaster: public RasterObject /* Raster and associated points to sample, used by batch sampling */ typedef struct UniqueRaster { - bool dataIsElevation; - uint64_t fileId; // file dictionary id + const raster_info_t* rinfo; std::vector pointSamples; // vector of samples for each point in this raster - explicit UniqueRaster(bool _dataIsElevation, uint64_t _fileId): - dataIsElevation(_dataIsElevation), fileId(_fileId) {} + explicit UniqueRaster(const raster_info_t* _rinfo): rinfo(_rinfo) {} } unique_raster_t; typedef Ordering GroupOrdering; diff --git a/packages/geo/RasterObject.h b/packages/geo/RasterObject.h index 7e8e03c5..13803719 100644 --- a/packages/geo/RasterObject.h +++ b/packages/geo/RasterObject.h @@ -110,6 +110,11 @@ class RasterObject: public LuaObject virtual uint8_t* getPixels (uint32_t ulx, uint32_t uly, uint32_t xsize=0, uint32_t ysize=0, void* param=NULL); ~RasterObject (void) override; + bool hasBands (void) const + { + return parms->bands.length() > 0; + } + bool hasZonalStats (void) const { return parms->zonal_stats; diff --git a/packages/geo/RasterSample.h b/packages/geo/RasterSample.h index c6b09d87..9e92b9ac 100644 --- a/packages/geo/RasterSample.h +++ b/packages/geo/RasterSample.h @@ -50,6 +50,7 @@ class RasterSample double verticalShift; uint64_t fileId; uint32_t flags; + std::string bandName; struct zonal_t { @@ -82,6 +83,7 @@ class RasterSample { value = sample.value; flags = sample.flags; + bandName = sample.bandName; stats = sample.stats; } diff --git a/scripts/selftests/test_runner.lua b/scripts/selftests/test_runner.lua index 339da13c..1fe59b5b 100644 --- a/scripts/selftests/test_runner.lua +++ b/scripts/selftests/test_runner.lua @@ -72,7 +72,7 @@ if __icesat2__ and incloud then runner.script(icesat2_td .. "s3_driver.lua") end --- Run opendata Plugin Self Tests +-- Run OPENDATA Plugin Self Tests if __opendata__ and incloud then local opendata_td = td .. "../../datasets/opendata/selftests/" runner.script(opendata_td .. "worldcover_reader.lua") @@ -94,7 +94,7 @@ if __pgc__ and incloud then runner.script(pgc_td .. "subset_test.lua") end --- Run Landsat Plugin Self Tests +-- Run LANDSAT Plugin Self Tests if __landsat__ and incloud then local landsat_td = td .. "../../datasets/landsat/selftests/" runner.script(landsat_td .. "plugin_unittest.lua") @@ -114,6 +114,12 @@ if __gebco__ and incloud then runner.script(gebco_td .. "gebco_reader.lua") end +-- Run BLUETOPO Plugin Self Tests +if __bluetopo__ and incloud then + local bluetopo_td = td .. "../../datasets/bluetopo/selftests/" + runner.script(bluetopo_td.. "bluetopo_reader.lua") +end + -- Run Default Parameters Self Tests for all modules runner.script(td .. "parms_tojson.lua") From 4ec1d6c684c6a1e01213fa7769c28f5740a7743a Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Tue, 12 Nov 2024 22:00:33 +0000 Subject: [PATCH 03/17] work in progress, added banNum --- .../bluetopo/package/BlueTopoBathyRaster.cpp | 1 - datasets/landsat/package/LandsatHlsRaster.cpp | 4 +- datasets/landsat/package/LandsatHlsRaster.h | 2 + packages/geo/GdalRaster.cpp | 87 +++++++--- packages/geo/GdalRaster.h | 23 +-- packages/geo/GeoIndexedRaster.cpp | 150 +++++++++++++----- packages/geo/GeoIndexedRaster.h | 10 +- packages/geo/GeoJsonRaster.cpp | 9 +- packages/geo/GeoRaster.cpp | 10 +- packages/geo/GeoRaster.h | 2 +- packages/geo/RasterObject.cpp | 15 +- packages/geo/RasterObject.h | 5 +- 12 files changed, 223 insertions(+), 95 deletions(-) diff --git a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp index 73b701b9..344196dc 100644 --- a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp @@ -70,7 +70,6 @@ void BlueTopoBathyRaster::getIndexFile(const OGRGeometry* geo, std::string& file mlog(DEBUG, "Using %s", file.c_str()); } - /*---------------------------------------------------------------------------- * findRasters *----------------------------------------------------------------------------*/ diff --git a/datasets/landsat/package/LandsatHlsRaster.cpp b/datasets/landsat/package/LandsatHlsRaster.cpp index 73875528..c6fbe912 100644 --- a/datasets/landsat/package/LandsatHlsRaster.cpp +++ b/datasets/landsat/package/LandsatHlsRaster.cpp @@ -300,7 +300,7 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr cacheitem_t* item; if(cache.find(key, &item)) { - RasterSample* sample = item->sample; + RasterSample* sample = item->bandSample[0]; if(sample) { sample->flags = flags; @@ -315,7 +315,7 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr { sample->bandName = bandName; // TODO: Remove when sample contains band name slist->add(sample); - item->sample = NULL; + item->bandSample[0] = NULL; } } diff --git a/datasets/landsat/package/LandsatHlsRaster.h b/datasets/landsat/package/LandsatHlsRaster.h index 575ef729..bdeb25ff 100644 --- a/datasets/landsat/package/LandsatHlsRaster.h +++ b/datasets/landsat/package/LandsatHlsRaster.h @@ -100,6 +100,8 @@ class LandsatHlsRaster: public GeoIndexedRaster uint32_t getBatchGroupSamples(const rasters_group_t* rgroup, List* slist, uint32_t flags, uint32_t pointIndx) final { return _getGroupSamples(BATCH, rgroup, slist, flags, pointIndx);} + void getBands (std::vector& bands) final + { bands.clear(); } /* Landsat has bands in dedicated rasters */ private: diff --git a/packages/geo/GdalRaster.cpp b/packages/geo/GdalRaster.cpp index c2e718e1..3a8ad6ac 100644 --- a/packages/geo/GdalRaster.cpp +++ b/packages/geo/GdalRaster.cpp @@ -63,7 +63,6 @@ GdalRaster::GdalRaster(const GeoFields* _parms, const std::string& _fileName, do overrideCRS(cb), fileName (_fileName), dset (NULL), - band (NULL), dataIsElevation(_dataIsElevation), xsize (0), ysize (0), @@ -109,6 +108,28 @@ void GdalRaster::open(void) mlog(DEBUG, "Opened %s", fileName.c_str()); + const int bandCount = dset->GetRasterCount(); + if (bandCount == 0) + { + throw RunTimeException(CRITICAL, RTE_ERROR, "No bands found in raster: %s:", fileName.c_str()); + } + + /* Populate the mapping of band names to indices */ + for (int i = 1; i <= bandCount; ++i) + { + GDALRasterBand* band = dset->GetRasterBand(i); + CHECKPTR(band); + + const char* bandName = band->GetDescription(); + + /* Only store bands that have valid names */ + if (bandName && strlen(bandName) > 0) + { + bandMap[std::string(bandName)] = i; + mlog(DEBUG, "Band %d: %s", i, bandName); + } + } + /* Store information about raster */ xsize = dset->GetRasterXSize(); ysize = dset->GetRasterYSize(); @@ -140,9 +161,6 @@ void GdalRaster::open(void) parms->sampling_radius.value, MAX_SAMPLING_RADIUS_IN_PIXELS * static_cast(cellSize)); } - band = dset->GetRasterBand(1); - CHECKPTR(band); - /* Create coordinates transform for raster */ createTransform(); @@ -157,7 +175,7 @@ void GdalRaster::open(void) dset = NULL; OGRCoordinateTransformation::DestroyCT(transf); transf = NULL; - band = NULL; + bandMap.clear(); throw; } } @@ -165,7 +183,7 @@ void GdalRaster::open(void) /*---------------------------------------------------------------------------- * samplePOI *----------------------------------------------------------------------------*/ -RasterSample* GdalRaster::samplePOI(OGRPoint* poi) +RasterSample* GdalRaster::samplePOI(OGRPoint* poi, int bandNum) { RasterSample* sample = NULL; @@ -177,6 +195,9 @@ RasterSample* GdalRaster::samplePOI(OGRPoint* poi) if(dset == NULL) open(); + GDALRasterBand* band = dset->GetRasterBand(bandNum); + CHECKPTR(band); + const double z = poi->getZ(); // mlog(DEBUG, "Before transform x,y,z: (%.4lf, %.4lf, %.4lf)", poi->getX(), poi->getY(), poi->getZ()); if(poi->transform(transf) != OGRERR_NONE) @@ -192,12 +213,12 @@ RasterSample* GdalRaster::samplePOI(OGRPoint* poi) const double vertical_shift = z - poi->getZ(); sample = new RasterSample(gpsTime, fileId, vertical_shift); if(parms->sampling_algo == GRIORA_NearestNeighbour) - readPixel(poi, sample); + readPixel(poi, band, sample); else - resamplePixel(poi, sample); + resamplePixel(poi, band, sample); if(parms->zonal_stats) - computeZonalStats(poi, sample); + computeZonalStats(poi, band, sample); } else { @@ -218,7 +239,7 @@ RasterSample* GdalRaster::samplePOI(OGRPoint* poi) /*---------------------------------------------------------------------------- * subsetAOI *----------------------------------------------------------------------------*/ -RasterSubset* GdalRaster::subsetAOI(OGRPolygon* poly) +RasterSubset* GdalRaster::subsetAOI(OGRPolygon* poly, int bandNum) { /* * Notes on extent format: @@ -341,7 +362,7 @@ 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); } - subset = getSubset(ulx, uly, _xsize, _ysize); + subset = getSubset(ulx, uly, _xsize, _ysize, bandNum); } catch (const RunTimeException &e) { @@ -359,7 +380,7 @@ RasterSubset* GdalRaster::subsetAOI(OGRPolygon* poly) /*---------------------------------------------------------------------------- * getPixels *----------------------------------------------------------------------------*/ -uint8_t* GdalRaster::getPixels(uint32_t ulx, uint32_t uly, uint32_t _xsize, uint32_t _ysize) +uint8_t* GdalRaster::getPixels(uint32_t ulx, uint32_t uly, uint32_t _xsize, uint32_t _ysize, int bandNum) { /* Clear error status */ ssError = SS_NO_ERRORS; @@ -395,6 +416,8 @@ uint8_t* GdalRaster::getPixels(uint32_t ulx, uint32_t uly, uint32_t _xsize, uint if(uly + _ysize > ysize) throw RunTimeException(CRITICAL, RTE_ERROR, "rows out of bounds"); + GDALRasterBand* band = dset->GetRasterBand(bandNum); + CHECKPTR(band); const GDALDataType dtype = band->GetRasterDataType(); /* Make all uint64_t, with uint32_t got an overflow */ @@ -441,6 +464,20 @@ uint8_t* GdalRaster::getPixels(uint32_t ulx, uint32_t uly, uint32_t _xsize, uint } +/*---------------------------------------------------------------------------- + * getBandNumber + *----------------------------------------------------------------------------*/ +int GdalRaster::getBandNumber(const std::string& bandName) +{ + auto it = bandMap.find(bandName); + if(it == bandMap.end()) + { + mlog(WARNING, "Band name not found: %s", bandName.c_str()); + } + return it->second; +} + + /*---------------------------------------------------------------------------- * setCRSfromWkt *----------------------------------------------------------------------------*/ @@ -527,7 +564,7 @@ OGRPolygon GdalRaster::makeRectangle(double minx, double miny, double maxx, doub /*---------------------------------------------------------------------------- * readPixel *----------------------------------------------------------------------------*/ -void GdalRaster::readPixel(const OGRPoint* poi, RasterSample* sample) +void GdalRaster::readPixel(const OGRPoint* poi, GDALRasterBand* band, RasterSample* sample) { /* Use fast method recomended by GDAL docs to read individual pixel */ try @@ -540,6 +577,7 @@ void GdalRaster::readPixel(const OGRPoint* poi, RasterSample* sample) int xBlockSize = 0; int yBlockSize = 0; + band->GetBlockSize(&xBlockSize, &yBlockSize); /* Raster offsets to block of interest */ @@ -659,7 +697,7 @@ void GdalRaster::readPixel(const OGRPoint* poi, RasterSample* sample) /* Done reading, release block lock */ block->DropLock(); - if(nodataCheck(sample) && dataIsElevation) + if(nodataCheck(sample, band) && dataIsElevation) { sample->value += sample->verticalShift; } @@ -678,7 +716,7 @@ void GdalRaster::readPixel(const OGRPoint* poi, RasterSample* sample) /*---------------------------------------------------------------------------- * resamplePixel *----------------------------------------------------------------------------*/ -void GdalRaster::resamplePixel(const OGRPoint* poi, RasterSample* sample) +void GdalRaster::resamplePixel(const OGRPoint* poi, GDALRasterBand* band, RasterSample* sample) { try { @@ -727,8 +765,8 @@ void GdalRaster::resamplePixel(const OGRPoint* poi, RasterSample* sample) const bool validWindow = containsWindow(_x, _y, xsize, ysize, windowSize); if(validWindow) { - readWithRetry(_x, _y, windowSize, windowSize, &sample->value, 1, 1, &args); - if(nodataCheck(sample) && dataIsElevation) + readWithRetry(band, _x, _y, windowSize, windowSize, &sample->value, 1, 1, &args); + if(nodataCheck(sample, band) && dataIsElevation) { sample->value += sample->verticalShift; } @@ -736,7 +774,7 @@ void GdalRaster::resamplePixel(const OGRPoint* poi, RasterSample* sample) else { /* At least return pixel value if unable to resample raster */ - readPixel(poi, sample); + readPixel(poi, band, sample); } } catch (const RunTimeException &e) @@ -749,7 +787,7 @@ void GdalRaster::resamplePixel(const OGRPoint* poi, RasterSample* sample) /*---------------------------------------------------------------------------- * computeZonalStats *----------------------------------------------------------------------------*/ -void GdalRaster::computeZonalStats(const OGRPoint* poi, RasterSample* sample) +void GdalRaster::computeZonalStats(const OGRPoint* poi, GDALRasterBand* band, RasterSample* sample) { double *samplesArray = NULL; @@ -771,12 +809,13 @@ void GdalRaster::computeZonalStats(const OGRPoint* poi, RasterSample* sample) const bool validWindow = containsWindow(newx, newy, xsize, ysize, windowSize); if(validWindow) { - readWithRetry(newx, newy, windowSize, windowSize, samplesArray, windowSize, windowSize, &args); + readWithRetry(band, newx, newy, windowSize, windowSize, samplesArray, windowSize, windowSize, &args); /* One of the windows (raster or index data set) was valid. Compute zonal stats */ double min = std::numeric_limits::max(); double max = std::numeric_limits::min(); double sum = 0; + const double nodata = band->GetNoDataValue(); std::vector validSamples; /* @@ -872,7 +911,7 @@ void GdalRaster::computeZonalStats(const OGRPoint* poi, RasterSample* sample) /*---------------------------------------------------------------------------- * nodataCheck *----------------------------------------------------------------------------*/ -bool GdalRaster::nodataCheck(RasterSample* sample) +bool GdalRaster::nodataCheck(RasterSample* sample, GDALRasterBand* band) { /* * Replace nodata with NAN @@ -981,7 +1020,7 @@ bool GdalRaster::containsWindow(int x, int y, int maxx, int maxy, int windowSize /*---------------------------------------------------------------------------- * readWithRetry *----------------------------------------------------------------------------*/ -void GdalRaster::readWithRetry(int x, int y, int _xsize, int _ysize, void *data, int dataXsize, int dataYsize, GDALRasterIOExtraArg *args) +void GdalRaster::readWithRetry(GDALRasterBand* band, int x, int y, int _xsize, int _ysize, void *data, int dataXsize, int dataYsize, GDALRasterIOExtraArg *args) { /* * On AWS reading from S3 buckets may result in failed reads due to network issues/timeouts. @@ -1009,7 +1048,7 @@ void GdalRaster::readWithRetry(int x, int y, int _xsize, int _ysize, void *data, /*---------------------------------------------------------------------------- * getSubset *----------------------------------------------------------------------------*/ -RasterSubset* GdalRaster::getSubset(uint32_t ulx, uint32_t uly, uint32_t _xsize, uint32_t _ysize) +RasterSubset* GdalRaster::getSubset(uint32_t ulx, uint32_t uly, uint32_t _xsize, uint32_t _ysize, int bandNum) { RasterSubset* subset = NULL; char** options = NULL; @@ -1027,6 +1066,8 @@ RasterSubset* GdalRaster::getSubset(uint32_t ulx, uint32_t uly, uint32_t _xsize, } + GDALRasterBand* band = dset->GetRasterBand(bandNum); + CHECKPTR(band); const GDALDataType dtype = band->GetRasterDataType(); /* Calculate size of subset */ diff --git a/packages/geo/GdalRaster.h b/packages/geo/GdalRaster.h index 6136ce44..705dee1f 100644 --- a/packages/geo/GdalRaster.h +++ b/packages/geo/GdalRaster.h @@ -41,6 +41,7 @@ #include "RasterSubset.h" #include #include +#include /****************************************************************************** * Typedef and macros used by GDAL class @@ -104,9 +105,9 @@ class GdalRaster GdalRaster (const GeoFields* _parms, const std::string& _fileName, double _gpsTime, uint64_t _fileId, bool _dataIsElevation, overrideCRS_t cb, bbox_t* aoi_bbox_override=NULL); virtual ~GdalRaster (void); void open (void); - RasterSample* samplePOI (OGRPoint* poi); - RasterSubset* subsetAOI (OGRPolygon* poly); - uint8_t* getPixels (uint32_t ulx, uint32_t uly, uint32_t _xsize, uint32_t _ysize); + RasterSample* samplePOI (OGRPoint* poi, int bandNum); + RasterSubset* subsetAOI (OGRPolygon* poly, int bandNum); + uint8_t* getPixels (uint32_t ulx, uint32_t uly, uint32_t _xsize, uint32_t _ysize, int bandNum); const std::string& getFileName (void) const { return fileName;} int getRows (void) const { return ysize; } int getCols (void) const { return xsize; } @@ -116,6 +117,7 @@ class GdalRaster bool isElevation (void) const { return dataIsElevation; } overrideCRS_t getOverrideCRS (void) const { return overrideCRS; } double getGpsTime (void) const { return gpsTime; } + int getBandNumber (const std::string& bandName); /*-------------------------------------------------------------------- * Static Methods @@ -145,7 +147,6 @@ class GdalRaster std::string fileName; GDALDataset *dset; - GDALRasterBand* band; bool dataIsElevation; uint32_t xsize; uint32_t ysize; @@ -157,19 +158,21 @@ class GdalRaster double invGeoTransform[6]; uint32_t ssError; + std::unordered_map bandMap; /* Maps raster band names to band numbers */ + /*-------------------------------------------------------------------- * Methods *--------------------------------------------------------------------*/ - void readPixel (const OGRPoint* poi, RasterSample* sample); - void resamplePixel (const OGRPoint* poi, RasterSample* sample); - void computeZonalStats (const OGRPoint* poi, RasterSample* sample); - inline bool nodataCheck (RasterSample* sample); + void readPixel (const OGRPoint* poi, GDALRasterBand* band, RasterSample* sample); + void resamplePixel (const OGRPoint* poi, GDALRasterBand* band, RasterSample* sample); + void computeZonalStats (const OGRPoint* poi, GDALRasterBand* band, RasterSample* sample); + inline bool nodataCheck (RasterSample* sample, GDALRasterBand* band); void createTransform (void); int radius2pixels (int _radius) const; static inline bool containsWindow(int x, int y, int maxx, int maxy, int windowSize); - inline void readWithRetry (int x, int y, int xsize, int ysize, void* data, int dataXsize, int dataYsize, GDALRasterIOExtraArg* args); - RasterSubset* getSubset (uint32_t ulx, uint32_t uly, uint32_t _xsize, uint32_t _ysize); + inline void readWithRetry (GDALRasterBand* band, int x, int y, int xsize, int ysize, void* data, int dataXsize, int dataYsize, GDALRasterIOExtraArg* args); + RasterSubset* getSubset (uint32_t ulx, uint32_t uly, uint32_t _xsize, uint32_t _ysize, int bandNum); 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/packages/geo/GeoIndexedRaster.cpp b/packages/geo/GeoIndexedRaster.cpp index ed9b2139..b13bd7c4 100644 --- a/packages/geo/GeoIndexedRaster.cpp +++ b/packages/geo/GeoIndexedRaster.cpp @@ -183,15 +183,22 @@ uint32_t GeoIndexedRaster::getSamples(const MathLib::point_3d_t& point, int64_t const char* key = cache.first(&item); while(key != NULL) { - if(item->sample) + for(uint32_t i = 0; i < item->bandSample.size(); i++) { - delete item->sample; - item->sample = NULL; + if(item->bandSample[i]) + { + delete item->bandSample[i]; + item->bandSample[i] = NULL; + } } - if(item->subset) + + for(uint32_t i = 0; i < item->bandSubset.size(); i++) { - delete item->subset; - item->subset = NULL; + if(item->bandSubset[i]) + { + delete item->bandSubset[i]; + item->bandSubset[i] = NULL; + } } key = cache.next(&item); } @@ -480,12 +487,15 @@ void GeoIndexedRaster::getGroupSamples(const rasters_group_t* rgroup, Listsample; - if(_sample) + for(uint32_t i = 0; i < item->bandSample.size(); i++) { - _sample->flags = flags; - slist.add(_sample); - item->sample = NULL; + RasterSample* _sample = item->bandSample[i]; + if(_sample) + { + _sample->flags = flags; + slist.add(_sample); + item->bandSample[i] = NULL; + } } /* Get sampling/subset error status */ @@ -512,11 +522,14 @@ void GeoIndexedRaster::getGroupSubsets(const rasters_group_t* rgroup, Listsubset; - if(subset) + for(uint32_t i = 0; i < item->bandSubset.size(); i++) { - slist.add(subset); - item->subset = NULL; + RasterSubset* subset = item->bandSubset[i]; + if(subset) + { + slist.add(subset); + item->bandSubset[i] = NULL; + } } /* Get sampling/subset error status */ @@ -539,7 +552,10 @@ uint32_t GeoIndexedRaster::getGroupFlags(const rasters_group_t* rgroup) const char* key = fileDict.get(rinfo.fileId); if(cache.find(key, &item)) { - const RasterSample* _sample = item->sample; + /* TODO: rinfo needs to specify which band has the flags, + * for now just use the first + */ + const RasterSample* _sample = item->bandSample[0]; if(_sample) { return _sample->value; @@ -685,7 +701,7 @@ bool GeoIndexedRaster::openGeoIndex(const OGRGeometry* geo, const std::vectorentry; if(entry != NULL) { + std::vector bands; + reader->obj->getRasterBands(reader->obj, entry->raster, bands); + + GdalRaster* raster = entry->raster; if(GdalRaster::ispoint(reader->geo)) { - entry->sample = entry->raster->samplePOI(dynamic_cast(reader->geo)); + /* Sample raster bands */ + for(size_t i = 0; i < bands.size(); i++) + { + RasterSample* sample = raster->samplePOI(dynamic_cast(reader->geo), bands[i]); + entry->bandSample.push_back(sample); + } } else if(GdalRaster::ispoly(reader->geo)) { - entry->subset = entry->raster->subsetAOI(dynamic_cast(reader->geo)); - if(entry->subset) + /* Subset raster bands */ + for(size_t i = 0; i < bands.size(); i++) { - /* - * Create new GeoRaster object for subsetted raster - * Use NULL for LuaState, using parent's causes memory corruption - * NOTE: cannot use RasterObject::cppCreate(parms) here, it would create - * new GeoIndexRaster with the same file path as parent raster. - */ - entry->subset->robj = new GeoRaster(NULL, - reader->obj->rqstParms, - reader->obj->samplerKey, - entry->subset->rasterName, - entry->raster->getGpsTime(), - entry->raster->isElevation(), - entry->raster->getOverrideCRS()); - - /* RequestFields are shared with subsseted raster and other readers */ - GeoIndexedRaster::referenceLuaObject(reader->obj->rqstParms); + RasterSubset* subset = raster->subsetAOI(dynamic_cast(reader->geo), bands[i]); + if(subset) + { + /* + * Create new GeoRaster object for subsetted raster + * Use NULL for LuaState, using parent's causes memory corruption + * NOTE: cannot use RasterObject::cppCreate(parms) here, it would create + * new GeoIndexRaster with the same file path as parent raster. + */ + subset->robj = new GeoRaster(NULL, + reader->obj->rqstParms, + reader->obj->samplerKey, + subset->rasterName, + raster->getGpsTime(), + raster->isElevation(), + raster->getOverrideCRS()); + + entry->bandSubset.push_back(subset); + + /* RequestFields are shared with subsseted raster and other readers */ + GeoIndexedRaster::referenceLuaObject(reader->obj->rqstParms); + } } } entry->enabled = false; /* raster samples/subsetted */ @@ -964,11 +995,18 @@ void* GeoIndexedRaster::batchReaderThread(void *param) ur->rinfo->dataIsElevation, breader->obj->crscb); + std::vector bands; + breader->obj->getRasterBands(breader->obj, raster, bands); + /* Sample all points for this raster */ - for(point_sample_t& ps : ur->pointSamples) + for(point_sample_t& ps : ur->pointSamples) //TODO: fxi this, must be vector of samples { - ps.sample = raster->samplePOI(&ps.point); - ps.ssErrors |= raster->getSSerror(); + /* Sample all bands for this point */ + for(size_t i = 0; i < bands.size(); i++) + { + ps.sample = raster->samplePOI(&ps.point, bands[i]); + ps.ssErrors |= raster->getSSerror(); + } } delete raster; @@ -1214,12 +1252,14 @@ bool GeoIndexedRaster::updateCache(uint32_t& rasters2sample, const GroupOrdering static_cast(rgroup->gpsTime / 1000), rinfo.fileId, rinfo.dataIsElevation, crscb, &bbox); - item->sample = NULL; - item->subset = NULL; const bool status = cache.add(key, item); assert(status); (void)status; // cannot fail; prevents linter warnings } + /* Clear from previous run */ + item->bandSample.clear(); + item->bandSample.clear(); + /* Mark as Enabled */ item->enabled = true; rasters2sample++; @@ -1786,3 +1826,31 @@ bool GeoIndexedRaster::collectSamples(const std::vector& pointsG +/*---------------------------------------------------------------------------- + * getRasterBands + *----------------------------------------------------------------------------*/ +void GeoIndexedRaster::getRasterBands(RasterObject* robj, GdalRaster* raster, std::vector& bands) +{ + std::vector bandsNames; + robj->getBands(bandsNames); + + if(bandsNames.empty()) + { + /* Default to first band */ + bands.push_back(1); + } + else + { + for(const std::string& bname : bandsNames) + { + const int bandNum = raster->getBandNumber(bname); + if(bandNum > 0) + { + bands.push_back(bandNum); + } + } + } +} + + + diff --git a/packages/geo/GeoIndexedRaster.h b/packages/geo/GeoIndexedRaster.h index e1ecfe1b..dc7cb3a8 100644 --- a/packages/geo/GeoIndexedRaster.h +++ b/packages/geo/GeoIndexedRaster.h @@ -134,10 +134,10 @@ class GeoIndexedRaster: public RasterObject /* Cache used by serial sampling code */ typedef struct CacheItem { - bool enabled; - RasterSample* sample; - RasterSubset* subset; - GdalRaster* raster; + bool enabled; + std::vector bandSample; + std::vector bandSubset; + GdalRaster* raster; ~CacheItem(void) {delete raster;} } cacheitem_t; @@ -326,6 +326,8 @@ class GeoIndexedRaster: public RasterObject bool collectSamples (const std::vector& pointsGroups, List& sllist); + void getRasterBands (RasterObject* robj, GdalRaster* raster, std::vector& bands); + }; #endif /* __geo_indexed_raster__ */ diff --git a/packages/geo/GeoJsonRaster.cpp b/packages/geo/GeoJsonRaster.cpp index 1953f4ff..e17a5380 100644 --- a/packages/geo/GeoJsonRaster.cpp +++ b/packages/geo/GeoJsonRaster.cpp @@ -216,8 +216,7 @@ GeoJsonRaster::GeoJsonRaster(lua_State* L, RequestFields* rqst_parms, const char rasterDset->SetProjection(wkt); CPLFree(wkt); - const int bandInx = 1; /* Band index starts at 1, not 0 */ - GDALRasterBand *rb = rasterDset->GetRasterBand(bandInx); + GDALRasterBand *rb = rasterDset->GetRasterBand(1); CHECKPTR(rb); rb->SetNoDataValue(RASTER_NODATA_VALUE); @@ -228,7 +227,7 @@ GeoJsonRaster::GeoJsonRaster(lua_State* L, RequestFields* rqst_parms, const char const int BANDCNT = 1; int bandlist[BANDCNT]; - bandlist[0] = bandInx; + bandlist[0] = 1; OGRLayer *layers[BANDCNT]; layers[0] = srcLayer; @@ -280,8 +279,8 @@ GeoJsonRaster::GeoJsonRaster(lua_State* L, RequestFields* rqst_parms, const char /* Cleanup */ VSIUnlink(jsonFile.c_str()); - if(jsonDset) GDALClose(reinterpret_cast(jsonDset)); - if(rasterDset) GDALClose(reinterpret_cast(rasterDset)); + GDALClose(reinterpret_cast(jsonDset)); + GDALClose(reinterpret_cast(rasterDset)); if(!rasterCreated) throw RunTimeException(CRITICAL, RTE_ERROR, "GeoJsonRaster failed"); diff --git a/packages/geo/GeoRaster.cpp b/packages/geo/GeoRaster.cpp index 23ac319c..d9179171 100644 --- a/packages/geo/GeoRaster.cpp +++ b/packages/geo/GeoRaster.cpp @@ -87,7 +87,7 @@ uint32_t GeoRaster::getSamples(const MathLib::point_3d_t& point, int64_t gps, Li try { OGRPoint ogr_point(point.x, point.y, point.z); - RasterSample* sample = raster.samplePOI(&ogr_point); + RasterSample* sample = raster.samplePOI(&ogr_point, 1); //TODO: fix this if(sample) slist.add(sample); } catch (const RunTimeException &e) @@ -117,7 +117,7 @@ uint32_t GeoRaster::getSubsets(const MathLib::extent_t& extent, int64_t gps, Lis OGRPolygon poly = GdalRaster::makeRectangle(extent.ll.x, extent.ll.y, extent.ur.x, extent.ur.y); /* Get subset rasters, if none found, return */ - RasterSubset* subset = raster.subsetAOI(&poly); + RasterSubset* subset = raster.subsetAOI(&poly, 1); //TODO: fix this if(subset) { /* @@ -133,10 +133,10 @@ uint32_t GeoRaster::getSubsets(const MathLib::extent_t& extent, int64_t gps, Lis raster.getGpsTime(), raster.isElevation(), raster.getOverrideCRS()); + slist.add(subset); /* RequestFields are shared with subsseted raster */ referenceLuaObject(rqstParms); - slist.add(subset); } } catch (const RunTimeException &e) @@ -155,7 +155,7 @@ uint32_t GeoRaster::getSubsets(const MathLib::extent_t& extent, int64_t gps, Lis /*---------------------------------------------------------------------------- * getPixels *----------------------------------------------------------------------------*/ -uint8_t* GeoRaster::getPixels(uint32_t ulx, uint32_t uly, uint32_t xsize, uint32_t ysize, void* param) +uint8_t* GeoRaster::getPixels(uint32_t ulx, uint32_t uly, uint32_t xsize, uint32_t ysize, int bandNum, void* param) { static_cast(param); uint8_t* data = NULL; @@ -165,7 +165,7 @@ uint8_t* GeoRaster::getPixels(uint32_t ulx, uint32_t uly, uint32_t xsize, uint32 /* Enable multi-threaded decompression in Gtiff driver */ CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"); - data = raster.getPixels(ulx, uly, xsize, ysize); + data = raster.getPixels(ulx, uly, xsize, ysize, bandNum); /* Disable multi-threaded decompression in Gtiff driver */ CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "1"); diff --git a/packages/geo/GeoRaster.h b/packages/geo/GeoRaster.h index 8134f9d1..a6a930c2 100644 --- a/packages/geo/GeoRaster.h +++ b/packages/geo/GeoRaster.h @@ -65,7 +65,7 @@ class GeoRaster: public RasterObject ~GeoRaster (void) override; uint32_t getSamples (const MathLib::point_3d_t& point, int64_t gps, List& slist, void* param=NULL) final; uint32_t getSubsets (const MathLib::extent_t& extent, int64_t gps, List& slist, void* param=NULL) final; - uint8_t* getPixels (uint32_t ulx, uint32_t uly, uint32_t xsize=0, uint32_t ysize=0, void* param=NULL) override; + uint8_t* getPixels (uint32_t ulx, uint32_t uly, uint32_t xsize=0, uint32_t ysize=0, int bandNum=1, void* param=NULL) override; uint32_t getRows (void) const { return raster.getRows(); } uint32_t getCols (void) const { return raster.getCols(); } diff --git a/packages/geo/RasterObject.cpp b/packages/geo/RasterObject.cpp index be066279..67caba4f 100644 --- a/packages/geo/RasterObject.cpp +++ b/packages/geo/RasterObject.cpp @@ -252,16 +252,29 @@ uint32_t RasterObject::getSamples(const std::vector& points, List< /*---------------------------------------------------------------------------- * getPixels *----------------------------------------------------------------------------*/ -uint8_t* RasterObject::getPixels(uint32_t ulx, uint32_t uly, uint32_t xsize, uint32_t ysize, void* param) +uint8_t* RasterObject::getPixels(uint32_t ulx, uint32_t uly, uint32_t xsize, uint32_t ysize, int bandNum, void* param) { static_cast(ulx); static_cast(uly); static_cast(xsize); static_cast(ysize); + static_cast(bandNum); static_cast(param); return NULL; } +/*---------------------------------------------------------------------------- + * getBands + *----------------------------------------------------------------------------*/ +void RasterObject::getBands(std::vector& bands) +{ + for(long i = 0; i < parms->bands.length(); i++) + { + const std::string& bandName = parms->bands[i]; + bands.push_back(bandName); + } +} + /*---------------------------------------------------------------------------- * Destructor *----------------------------------------------------------------------------*/ diff --git a/packages/geo/RasterObject.h b/packages/geo/RasterObject.h index 13803719..0c890869 100644 --- a/packages/geo/RasterObject.h +++ b/packages/geo/RasterObject.h @@ -107,7 +107,8 @@ class RasterObject: public LuaObject virtual uint32_t getSamples (const MathLib::point_3d_t& point, int64_t gps, List& slist, void* param=NULL) = 0; virtual uint32_t getSamples (const std::vector& points, List& sllist, void* param=NULL); virtual uint32_t getSubsets (const MathLib::extent_t& extent, int64_t gps, List& slist, void* param=NULL) = 0; - virtual uint8_t* getPixels (uint32_t ulx, uint32_t uly, uint32_t xsize=0, uint32_t ysize=0, void* param=NULL); + virtual uint8_t* getPixels (uint32_t ulx, uint32_t uly, uint32_t xsize=0, uint32_t ysize=0, int bandNum=1, void* param=NULL); + virtual void getBands (std::vector& bands); ~RasterObject (void) override; bool hasBands (void) const @@ -180,7 +181,7 @@ class RasterObject: public LuaObject static void getThreadsRanges(std::vector& ranges, uint32_t num, uint32_t minPerThread, uint32_t maxNumThreads); - void fileDictSetSamples(List* slist); + void fileDictSetSamples(List* slist); /*-------------------------------------------------------------------- From cf1e99fe95a1f957ef2acc0072cf7dc42e3fa312 Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Wed, 13 Nov 2024 13:51:14 +0000 Subject: [PATCH 04/17] implemented band samples for index rasters --- datasets/landsat/package/LandsatHlsRaster.cpp | 84 ++++----- datasets/landsat/package/LandsatHlsRaster.h | 4 +- packages/geo/GdalRaster.h | 2 +- packages/geo/GeoIndexedRaster.cpp | 160 ++++++++++-------- packages/geo/GeoIndexedRaster.h | 24 +-- packages/geo/GeoRaster.cpp | 55 +++--- packages/geo/RasterObject.cpp | 36 ++++ packages/geo/RasterObject.h | 4 +- 8 files changed, 216 insertions(+), 153 deletions(-) diff --git a/datasets/landsat/package/LandsatHlsRaster.cpp b/datasets/landsat/package/LandsatHlsRaster.cpp index c6fbe912..48d86517 100644 --- a/datasets/landsat/package/LandsatHlsRaster.cpp +++ b/datasets/landsat/package/LandsatHlsRaster.cpp @@ -291,6 +291,9 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr double swir16; green = red = nir08 = swir16 = invalid; + /* Landsat rasters use only the first inner band */ + const int INNER_BAND_INDX = 0; + /* Collect samples for all rasters */ if(mode == SERIAL) { @@ -300,40 +303,41 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr cacheitem_t* item; if(cache.find(key, &item)) { - RasterSample* sample = item->bandSample[0]; - if(sample) - { - sample->flags = flags; + RasterSample* sample = item->bandSample[INNER_BAND_INDX]; - /* Is this band's sample to be returned to the user? */ - const char* bandName = rinfo.tag.c_str(); - auto it = bandsDict.find(bandName); - if(it != bandsDict.end()) - { - const bool returnBandSample = it->second; - if(returnBandSample) - { - sample->bandName = bandName; // TODO: Remove when sample contains band name - slist->add(sample); - item->bandSample[0] = NULL; - } - } + /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ + if(sample == NULL) continue;; - /* green and red bands are the same for L8 and S2 */ - if(rinfo.tag == "B03") green = sample->value; - if(rinfo.tag == "B04") red = sample->value; + sample->flags = flags; - if(isL8) - { - if(rinfo.tag == "B05") nir08 = sample->value; - if(rinfo.tag == "B06") swir16 = sample->value; - } - else /* Must be Sentinel2 */ + /* Is this band's sample to be returned to the user? */ + const char* bandName = rinfo.tag.c_str(); + auto it = bandsDict.find(bandName); + if(it != bandsDict.end()) + { + const bool returnBandSample = it->second; + if(returnBandSample) { - if(rinfo.tag == "B8A") nir08 = sample->value; - if(rinfo.tag == "B11") swir16 = sample->value; + sample->bandName = bandName; + slist->add(sample); + item->bandSample[INNER_BAND_INDX] = NULL; } } + + /* green and red bands are the same for L8 and S2 */ + if(rinfo.tag == "B03") green = sample->value; + if(rinfo.tag == "B04") red = sample->value; + + if(isL8) + { + if(rinfo.tag == "B05") nir08 = sample->value; + if(rinfo.tag == "B06") swir16 = sample->value; + } + else /* Must be Sentinel2 */ + { + if(rinfo.tag == "B8A") nir08 = sample->value; + if(rinfo.tag == "B11") swir16 = sample->value; + } } } } @@ -350,8 +354,10 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr { if(ps.pointIndex == pointIndx) { + RasterSample* sample = ps.bandSample[INNER_BAND_INDX]; + /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ - if(ps.sample == NULL) break; + if(sample == NULL) break; /* Is this band's sample to be returned to the user? */ const char* bandName = rinfo.tag.c_str(); @@ -363,18 +369,18 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr if(returnBandSample) { RasterSample* s; - if(!ps.sampleReturned.exchange(true)) + if(!ps.bandSampleReturned[INNER_BAND_INDX]->exchange(true)) { - s = ps.sample; + s = ps.bandSample[INNER_BAND_INDX]; } else { /* Sample has already been returned, must create a copy */ - s = new RasterSample(*ps.sample); + s = new RasterSample(*ps.bandSample[INNER_BAND_INDX]); } /* Set band name for this sample */ - s->bandName = bandName; //TODO: remove when sample contains band name + s->bandName = bandName; /* Set time for this sample */ s->time = rgroup->gpsTime / 1000; @@ -388,18 +394,18 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr errors |= ps.ssErrors; /* green and red bands are the same for L8 and S2 */ - if(rinfo.tag == "B03") green = ps.sample->value; - if(rinfo.tag == "B04") red = ps.sample->value; + if(rinfo.tag == "B03") green = sample->value; + if(rinfo.tag == "B04") red = sample->value; if(isL8) { - if(rinfo.tag == "B05") nir08 = ps.sample->value; - if(rinfo.tag == "B06") swir16 = ps.sample->value; + if(rinfo.tag == "B05") nir08 = sample->value; + if(rinfo.tag == "B06") swir16 = sample->value; } else /* Must be Sentinel2 */ { - if(rinfo.tag == "B8A") nir08 = ps.sample->value; - if(rinfo.tag == "B11") swir16 = ps.sample->value; + if(rinfo.tag == "B8A") nir08 = sample->value; + if(rinfo.tag == "B11") swir16 = sample->value; } break; } diff --git a/datasets/landsat/package/LandsatHlsRaster.h b/datasets/landsat/package/LandsatHlsRaster.h index bdeb25ff..5a55b522 100644 --- a/datasets/landsat/package/LandsatHlsRaster.h +++ b/datasets/landsat/package/LandsatHlsRaster.h @@ -100,8 +100,8 @@ class LandsatHlsRaster: public GeoIndexedRaster uint32_t getBatchGroupSamples(const rasters_group_t* rgroup, List* slist, uint32_t flags, uint32_t pointIndx) final { return _getGroupSamples(BATCH, rgroup, slist, flags, pointIndx);} - void getBands (std::vector& bands) final - { bands.clear(); } /* Landsat has bands in dedicated rasters */ + void getInnerBands (std::vector& bands) final + { bands.clear(); } /* Landsat bands are in seperate rasters */ private: diff --git a/packages/geo/GdalRaster.h b/packages/geo/GdalRaster.h index 705dee1f..cc7cbb97 100644 --- a/packages/geo/GdalRaster.h +++ b/packages/geo/GdalRaster.h @@ -167,7 +167,7 @@ class GdalRaster void readPixel (const OGRPoint* poi, GDALRasterBand* band, RasterSample* sample); void resamplePixel (const OGRPoint* poi, GDALRasterBand* band, RasterSample* sample); void computeZonalStats (const OGRPoint* poi, GDALRasterBand* band, RasterSample* sample); - inline bool nodataCheck (RasterSample* sample, GDALRasterBand* band); + static inline bool nodataCheck (RasterSample* sample, GDALRasterBand* band); void createTransform (void); int radius2pixels (int _radius) const; static inline bool containsWindow(int x, int y, int maxx, int maxy, int windowSize); diff --git a/packages/geo/GeoIndexedRaster.cpp b/packages/geo/GeoIndexedRaster.cpp index b13bd7c4..3c73d36e 100644 --- a/packages/geo/GeoIndexedRaster.cpp +++ b/packages/geo/GeoIndexedRaster.cpp @@ -58,6 +58,35 @@ const char* GeoIndexedRaster::DATE_TAG = "datetime"; * PUBLIC METHODS ******************************************************************************/ +/*---------------------------------------------------------------------------- + * PointSample Constructor + *----------------------------------------------------------------------------*/ +GeoIndexedRaster::PointSample::PointSample(const OGRPoint& _point, int64_t _pointIndex): + point(_point), pointIndex(_pointIndex), ssErrors(SS_NO_ERRORS) {} + + +/*---------------------------------------------------------------------------- + * PointSample Copy Constructor + *----------------------------------------------------------------------------*/ +GeoIndexedRaster::PointSample::PointSample(const PointSample& ps): + point(ps.point), pointIndex(ps.pointIndex), bandSample(ps.bandSample), ssErrors(ps.ssErrors) +{ + bandSampleReturned.resize(ps.bandSampleReturned.size()); + + for (size_t i = 0; i < ps.bandSampleReturned.size(); ++i) + { + if(ps.bandSampleReturned[i]) + { + /* Create a new atomic with the value loaded from the original */ + bandSampleReturned[i] = std::make_unique>(ps.bandSampleReturned[i]->load()); + } + else + { + bandSampleReturned[i] = NULL; + } + } +} + /*---------------------------------------------------------------------------- * Reader Constructor *----------------------------------------------------------------------------*/ @@ -282,9 +311,12 @@ uint32_t GeoIndexedRaster::getSamples(const std::vector& points, L for(const point_sample_t& ps : ur->pointSamples) { /* Delete samples which have not been returned (quality masks, etc) */ - if(!ps.sampleReturned) + for(size_t i = 0; i < ps.bandSample.size(); i++) { - delete ps.sample; + if(!ps.bandSampleReturned[i]->load()) + { + delete ps.bandSample[i]; + } } } @@ -409,31 +441,34 @@ uint32_t GeoIndexedRaster::getBatchGroupSamples(const rasters_group_t* rgroup, L { if(ps.pointIndex == pointIndx) { - /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ - if(ps.sample == NULL) break; - - RasterSample* s; - if(!ps.sampleReturned.exchange(true)) - { - s = ps.sample; - } - else + for(size_t i = 0; i < ps.bandSample.size(); i++) { - /* Sample has already been returned, must create a copy */ - s = new RasterSample(*ps.sample); - } + /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ + if(ps.bandSample[i] == NULL) continue;; - /* Set time for this sample */ - s->time = rgroup->gpsTime / 1000; + RasterSample* s; + if(!ps.bandSampleReturned[i]->exchange(true)) + { + s = ps.bandSample[i]; + } + else + { + /* Sample has already been returned, must create a copy */ + s = new RasterSample(*ps.bandSample[i]); + } + + /* Set time for this sample */ + s->time = rgroup->gpsTime / 1000; - /* Set flags for this sample, add it to the list */ - s->flags = flags; - slist->add(s); - errors |= ps.ssErrors; + /* Set flags for this sample, add it to the list */ + s->flags = flags; + slist->add(s); + errors |= ps.ssErrors; + } /* - * There is only one raster with VALUE_TAG and one with FLAGS_TAG in a group. - * If group has other type rasters the dataset must override the getGroupFlags method. + * This function assumes that there is only one raster with VALUE_TAG in a group. + * If group has other value rasters the dataset must override this function. */ return errors; } @@ -461,11 +496,17 @@ uint32_t GeoIndexedRaster::getBatchGroupFlags(const rasters_group_t* rgroup, uin { if(ps.pointIndex == pointIndx) { + /* + * This function assumes that there is only one raster with FLAGS_TAG in a group. + * The flags value must be in the first band. + * If these assumptions are not met the dataset must override this function. + */ + RasterSample* s = ps.bandSample[0]; + /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ - if(ps.sample) - { - return ps.sample->value; - } + if(s == NULL) break;; + + return s->value; } } } @@ -501,10 +542,10 @@ void GeoIndexedRaster::getGroupSamples(const rasters_group_t* rgroup, Listraster->getSSerror(); - /* - * There is only one raster with VALUE_TAG and one with FLAGS_TAG in a group. - * If group has other type rasters the dataset must override the getGroupFlags method. - */ + /* + * This function assumes that there is only one raster with VALUE_TAG in a group. + * If group has other value rasters the dataset must override this function. + */ break; } } @@ -552,8 +593,10 @@ uint32_t GeoIndexedRaster::getGroupFlags(const rasters_group_t* rgroup) const char* key = fileDict.get(rinfo.fileId); if(cache.find(key, &item)) { - /* TODO: rinfo needs to specify which band has the flags, - * for now just use the first + /* + * This function assumes that there is only one raster with FLAGS_TAG in a group. + * The flags value must be in the first band. + * If these assumptions are not met the dataset must override this function. */ const RasterSample* _sample = item->bandSample[0]; if(_sample) @@ -913,24 +956,24 @@ void* GeoIndexedRaster::readerThread(void *param) if(entry != NULL) { std::vector bands; - reader->obj->getRasterBands(reader->obj, entry->raster, bands); + reader->obj->getInnerBands(entry->raster, bands); GdalRaster* raster = entry->raster; if(GdalRaster::ispoint(reader->geo)) { /* Sample raster bands */ - for(size_t i = 0; i < bands.size(); i++) + for(const int bandNum : bands) { - RasterSample* sample = raster->samplePOI(dynamic_cast(reader->geo), bands[i]); + RasterSample* sample = raster->samplePOI(dynamic_cast(reader->geo), bandNum); entry->bandSample.push_back(sample); } } else if(GdalRaster::ispoly(reader->geo)) { /* Subset raster bands */ - for(size_t i = 0; i < bands.size(); i++) + for(const int bandNum : bands) { - RasterSubset* subset = raster->subsetAOI(dynamic_cast(reader->geo), bands[i]); + RasterSubset* subset = raster->subsetAOI(dynamic_cast(reader->geo), bandNum); if(subset) { /* @@ -996,15 +1039,17 @@ void* GeoIndexedRaster::batchReaderThread(void *param) breader->obj->crscb); std::vector bands; - breader->obj->getRasterBands(breader->obj, raster, bands); + breader->obj->getInnerBands(raster, bands); /* Sample all points for this raster */ - for(point_sample_t& ps : ur->pointSamples) //TODO: fxi this, must be vector of samples + for(point_sample_t& ps : ur->pointSamples) { /* Sample all bands for this point */ - for(size_t i = 0; i < bands.size(); i++) + for(const int bandNum : bands) { - ps.sample = raster->samplePOI(&ps.point, bands[i]); + RasterSample* sample = raster->samplePOI(&ps.point, bandNum); + ps.bandSample.push_back(sample); + ps.bandSampleReturned.emplace_back(std::make_unique>(false)); ps.ssErrors |= raster->getSSerror(); } } @@ -1822,35 +1867,4 @@ bool GeoIndexedRaster::collectSamples(const std::vector& pointsG mlog(INFO, "Populated sllist with %d lists of samples, time: %lf", sllist.length(), perfStats.collectSamplesTime); return true; -} - - - -/*---------------------------------------------------------------------------- - * getRasterBands - *----------------------------------------------------------------------------*/ -void GeoIndexedRaster::getRasterBands(RasterObject* robj, GdalRaster* raster, std::vector& bands) -{ - std::vector bandsNames; - robj->getBands(bandsNames); - - if(bandsNames.empty()) - { - /* Default to first band */ - bands.push_back(1); - } - else - { - for(const std::string& bname : bandsNames) - { - const int bandNum = raster->getBandNumber(bname); - if(bandNum > 0) - { - bands.push_back(bandNum); - } - } - } -} - - - +} \ No newline at end of file diff --git a/packages/geo/GeoIndexedRaster.h b/packages/geo/GeoIndexedRaster.h index dc7cb3a8..3dd17e85 100644 --- a/packages/geo/GeoIndexedRaster.h +++ b/packages/geo/GeoIndexedRaster.h @@ -74,18 +74,15 @@ class GeoIndexedRaster: public RasterObject /* Raster Sample used by batch sampling */ typedef struct PointSample { - OGRPoint point; - int64_t pointIndex; // index to the user provided list of points to sample - RasterSample* sample; // sample created by the batch reader thread - std::atomic sampleReturned;// multiple rasters may share the same sample, - // this flag is used to avoid returning the same sample, if set a copy of the sample is returned - uint32_t ssErrors; // sampling errors + OGRPoint point; + int64_t pointIndex; // index to the user provided list of points to sample + std::vector bandSample; // vector of samples for each band + std::vector>> bandSampleReturned; // multiple rasters may share the same sample, + // these flags are used to avoid returning the same sample, if set a copy of the sample is returned + uint32_t ssErrors; // sampling errors - PointSample(const OGRPoint& _point, int64_t _pointIndex): - point(_point), pointIndex(_pointIndex), sample(NULL), sampleReturned(false), ssErrors(SS_NO_ERRORS) {} - - PointSample(const PointSample& ps): - point(ps.point), pointIndex(ps.pointIndex), sample(ps.sample), sampleReturned(ps.sampleReturned.load()), ssErrors(ps.ssErrors) {} + PointSample(const OGRPoint& _point, int64_t _pointIndex); + PointSample(const PointSample& ps); } point_sample_t; @@ -94,7 +91,7 @@ class GeoIndexedRaster: public RasterObject /** Raster information needed for sampling */ typedef struct RasterInfo { bool dataIsElevation; - std::string tag; // "Dem", "Flags", "Date" + std::string tag; // value, fmask, etc uint64_t fileId; // file dictionary id UniqueRaster* uraster; // Pointer to the unique raster which contains the sample for this raster @@ -325,9 +322,6 @@ class GeoIndexedRaster: public RasterObject bool collectSamples (const std::vector& pointsGroups, List& sllist); - - void getRasterBands (RasterObject* robj, GdalRaster* raster, std::vector& bands); - }; #endif /* __geo_indexed_raster__ */ diff --git a/packages/geo/GeoRaster.cpp b/packages/geo/GeoRaster.cpp index d9179171..6aa1d778 100644 --- a/packages/geo/GeoRaster.cpp +++ b/packages/geo/GeoRaster.cpp @@ -87,8 +87,14 @@ uint32_t GeoRaster::getSamples(const MathLib::point_3d_t& point, int64_t gps, Li try { OGRPoint ogr_point(point.x, point.y, point.z); - RasterSample* sample = raster.samplePOI(&ogr_point, 1); //TODO: fix this - if(sample) slist.add(sample); + + std::vector bands; + getInnerBands(&raster, bands); + for(const int bandNum : bands) + { + RasterSample* sample = raster.samplePOI(&ogr_point, bandNum); + if(sample) slist.add(sample); + } } catch (const RunTimeException &e) { @@ -116,27 +122,32 @@ uint32_t GeoRaster::getSubsets(const MathLib::extent_t& extent, int64_t gps, Lis { OGRPolygon poly = GdalRaster::makeRectangle(extent.ll.x, extent.ll.y, extent.ur.x, extent.ur.y); - /* Get subset rasters, if none found, return */ - RasterSubset* subset = raster.subsetAOI(&poly, 1); //TODO: fix this - if(subset) + std::vector bands; + getInnerBands(&raster, bands); + for(const int bandNum : bands) { - /* - * Create new GeoRaster object for subsetted raster - * Use NULL for LuaState, using parent's causes memory corruption - * NOTE: cannot use RasterObject::cppCreate(parms) here, - * it would create subsetted raster with the same file path as parent raster. - */ - subset->robj = new GeoRaster(NULL, - rqstParms, - samplerKey, - subset->rasterName, - raster.getGpsTime(), - raster.isElevation(), - raster.getOverrideCRS()); - slist.add(subset); - - /* RequestFields are shared with subsseted raster */ - referenceLuaObject(rqstParms); + /* Get subset rasters, if none found, return */ + RasterSubset* subset = raster.subsetAOI(&poly, bandNum); + if(subset) + { + /* + * Create new GeoRaster object for subsetted raster + * Use NULL for LuaState, using parent's causes memory corruption + * NOTE: cannot use RasterObject::cppCreate(parms) here, + * it would create subsetted raster with the same file path as parent raster. + */ + subset->robj = new GeoRaster(NULL, + rqstParms, + samplerKey, + subset->rasterName, + raster.getGpsTime(), + raster.isElevation(), + raster.getOverrideCRS()); + slist.add(subset); + + /* RequestFields are shared with subsseted raster */ + referenceLuaObject(rqstParms); + } } } catch (const RunTimeException &e) diff --git a/packages/geo/RasterObject.cpp b/packages/geo/RasterObject.cpp index 67caba4f..8960b370 100644 --- a/packages/geo/RasterObject.cpp +++ b/packages/geo/RasterObject.cpp @@ -275,6 +275,42 @@ void RasterObject::getBands(std::vector& bands) } } +/*---------------------------------------------------------------------------- + * getInnerBands + *----------------------------------------------------------------------------*/ +void RasterObject::getInnerBands(std::vector& bands) +{ + return getBands(bands); +} + +/*---------------------------------------------------------------------------- + * getInnerBands + *----------------------------------------------------------------------------*/ +void RasterObject::getInnerBands(void* rptr, std::vector& bands) +{ + GdalRaster* raster = static_cast(rptr); + + std::vector bandsNames; + getInnerBands(bandsNames); + + if(bandsNames.empty()) + { + /* Default to first band */ + bands.push_back(1); + } + else + { + for(const std::string& bname : bandsNames) + { + const int bandNum = raster->getBandNumber(bname); + if(bandNum > 0) + { + bands.push_back(bandNum); + } + } + } +} + /*---------------------------------------------------------------------------- * Destructor *----------------------------------------------------------------------------*/ diff --git a/packages/geo/RasterObject.h b/packages/geo/RasterObject.h index 0c890869..98512c5a 100644 --- a/packages/geo/RasterObject.h +++ b/packages/geo/RasterObject.h @@ -108,7 +108,9 @@ class RasterObject: public LuaObject virtual uint32_t getSamples (const std::vector& points, List& sllist, void* param=NULL); virtual uint32_t getSubsets (const MathLib::extent_t& extent, int64_t gps, List& slist, void* param=NULL) = 0; virtual uint8_t* getPixels (uint32_t ulx, uint32_t uly, uint32_t xsize=0, uint32_t ysize=0, int bandNum=1, void* param=NULL); - virtual void getBands (std::vector& bands); + void getBands (std::vector& bands); + virtual void getInnerBands (std::vector& bands); + void getInnerBands (void* rptr, std::vector& bands); ~RasterObject (void) override; bool hasBands (void) const From 46c0cf0bbe7f33f3e4ee3a2bce828cb1a2cced51 Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Wed, 13 Nov 2024 21:17:00 +0000 Subject: [PATCH 05/17] implemented bands for bluetopo --- .../bluetopo/package/BlueTopoBathyRaster.cpp | 6 +- .../bluetopo/selftests/bluetopo_reader.lua | 25 ++- datasets/gebco/package/GebcoBathyRaster.cpp | 8 +- datasets/gedi/package/Gedi03Raster.h | 3 +- datasets/gedi/package/Gedi04bRaster.h | 3 +- datasets/landsat/package/LandsatHlsRaster.cpp | 2 +- .../package/EsaWorldCover10meterRaster.h | 3 +- .../package/MetaGlobalCanopy1meterRaster.h | 3 +- datasets/pgc/package/ArcticDemMosaicRaster.h | 3 +- datasets/pgc/package/PgcDemStripsRaster.cpp | 4 +- datasets/pgc/package/RemaDemMosaicRaster.h | 3 +- .../package/Usgs3dep1meterDemRaster.cpp | 2 +- packages/geo/GdalRaster.cpp | 176 ++++++++++-------- packages/geo/GdalRaster.h | 15 +- packages/geo/GeoIndexedRaster.cpp | 29 ++- packages/geo/GeoIndexedRaster.h | 11 +- packages/geo/GeoJsonRaster.cpp | 2 +- packages/geo/GeoRaster.cpp | 8 +- packages/geo/GeoRaster.h | 3 +- packages/geo/GeoUserRaster.cpp | 12 +- packages/geo/GeoUserRaster.h | 4 +- packages/geo/RasterObject.cpp | 1 + 22 files changed, 199 insertions(+), 127 deletions(-) diff --git a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp index 344196dc..e7f1a1b1 100644 --- a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp @@ -114,9 +114,9 @@ bool BlueTopoBathyRaster::findRasters(raster_finder_t* finder) } raster_info_t rinfo; - rinfo.dataIsElevation = true; - rinfo.tag = VALUE_TAG; - rinfo.fileId = finder->fileDict.add(fullPath); + rinfo.elevationBandNum = 1; + rinfo.tag = VALUE_TAG; + rinfo.fileId = finder->fileDict.add(fullPath); rgroup->infovect.push_back(rinfo); } rgroup->infovect.shrink_to_fit(); diff --git a/datasets/bluetopo/selftests/bluetopo_reader.lua b/datasets/bluetopo/selftests/bluetopo_reader.lua index 8422c0c0..04650864 100644 --- a/datasets/bluetopo/selftests/bluetopo_reader.lua +++ b/datasets/bluetopo/selftests/bluetopo_reader.lua @@ -15,16 +15,17 @@ local lons = {-80.87, -81.02, -89.66, -94.72} local lats = { 32.06, 31.86, 29.99, 29.35} height = 0 -local expDepth = { -2.02, -14.03, -4.28, -17.18} -local expFlags = { 3.32} +local expElevation = { -2.02, -14.03, -4.28, -17.18} +local expUncertainty = { 1.00, 4.37, 0.34, 1.32} +local expContributor = { 31156, 0, 24955, 45641} -- Tolerance for depth values is needed because of the different versions of PROJLIB -- The values are different in the different versions of the library. -- We add 1 meter of depth tolerance for the test to pass -local depth_tolerance = 1; +local elevation_tolerance = 0.1; print(string.format("\n--------------------------------\nTest: BlueTopo Correct Values\n--------------------------------")) -local dem = geo.raster(geo.parms({ asset = "bluetopo-bathy", algorithm = "NearestNeighbour", with_flags=true, sort_by_index = true })) +local dem = geo.raster(geo.parms({ asset = "bluetopo-bathy", algorithm = "NearestNeighbour", bands = {"Elevation", "Uncertainty", "Contributor"}, sort_by_index = true })) runner.check(dem ~= nil) for j, lon in ipairs(lons) do @@ -38,14 +39,20 @@ for j, lon in ipairs(lons) do else local el, fname for k, v in ipairs(tbl) do + band = v["band"] value = v["value"] fname = v["file"] - flags = v["flags"] - print(string.format("(%02d) (%6.2f, %6.2f) %6.2fm %4.2f %s", k, lon, lat, value, flags, fname)) - - assert( math.abs(value - expDepth[j]) < depth_tolerance, string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) - -- assert(flags == expFlags[j], string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) + print(string.format("(%02d) (%6.2f, %6.2f) %12s %10.2f %s", k, lon, lat, band, value, fname)) + + if band == "Elevation" then + runner.check(math.abs(value - expElevation[j]) < elevation_tolerance, string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) + elseif band == "Uncertainty" then + runner.check(math.abs(value - expUncertainty[j]) < elevation_tolerance, string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) + elseif band == "Contributor" then + runner.check(value == expContributor[j], string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) + end end + print("\n") end end diff --git a/datasets/gebco/package/GebcoBathyRaster.cpp b/datasets/gebco/package/GebcoBathyRaster.cpp index ac35cb7c..bfaf091d 100644 --- a/datasets/gebco/package/GebcoBathyRaster.cpp +++ b/datasets/gebco/package/GebcoBathyRaster.cpp @@ -92,9 +92,9 @@ bool GebcoBathyRaster::findRasters(raster_finder_t* finder) if(dataFile && (strlen(dataFile) > 0)) { raster_info_t rinfo; - rinfo.dataIsElevation = true; - rinfo.tag = VALUE_TAG; - rinfo.fileId = finder->fileDict.add(filePath + "/" + dataFile); + rinfo.elevationBandNum = 1; + rinfo.tag = VALUE_TAG; + rinfo.fileId = finder->fileDict.add(filePath + "/" + dataFile); rgroup->infovect.push_back(rinfo); } @@ -104,7 +104,7 @@ bool GebcoBathyRaster::findRasters(raster_finder_t* finder) if(flagsFile && (strlen(flagsFile) > 0)) { raster_info_t rinfo; - rinfo.dataIsElevation = false; + rinfo.flagsBandNum = 1; rinfo.tag = FLAGS_TAG; rinfo.fileId = finder->fileDict.add(filePath + "/" + flagsFile); rgroup->infovect.push_back(rinfo); diff --git a/datasets/gedi/package/Gedi03Raster.h b/datasets/gedi/package/Gedi03Raster.h index 8ddddade..e392f97b 100644 --- a/datasets/gedi/package/Gedi03Raster.h +++ b/datasets/gedi/package/Gedi03Raster.h @@ -66,7 +66,8 @@ class Gedi03Raster: public GeoRaster GeoRaster(L, rqst_parms, key, std::string(rqst_parms->geoFields(key)->asset.asset->getPath()).append("/").append(rqst_parms->geoFields(key)->asset.asset->getIndex()), TimeLib::datetime2gps(2022, 1, 19), - true /* Data is elevation */) {} + 1, /* elevationBandNum */ + GdalRaster::NO_BAND /* flagsBandNum */) {} }; #endif /* __gedi03_raster__ */ diff --git a/datasets/gedi/package/Gedi04bRaster.h b/datasets/gedi/package/Gedi04bRaster.h index ed026470..b273e537 100644 --- a/datasets/gedi/package/Gedi04bRaster.h +++ b/datasets/gedi/package/Gedi04bRaster.h @@ -65,7 +65,8 @@ class Gedi04bRaster: public GeoRaster GeoRaster(L, rqst_parms, key, std::string(rqst_parms->geoFields(key)->asset.asset->getPath()).append("/").append(rqst_parms->geoFields(key)->asset.asset->getIndex()), TimeLib::datetime2gps(2021, 8, 4), - true /* Data is elevation */) {} + 1, /* elevationBandNum */ + GdalRaster::NO_BAND /* flagsBandNum */) {} }; #endif /* __gedi04b_raster__ */ diff --git a/datasets/landsat/package/LandsatHlsRaster.cpp b/datasets/landsat/package/LandsatHlsRaster.cpp index 48d86517..5ee995e1 100644 --- a/datasets/landsat/package/LandsatHlsRaster.cpp +++ b/datasets/landsat/package/LandsatHlsRaster.cpp @@ -187,12 +187,12 @@ bool LandsatHlsRaster::findRasters(raster_finder_t* finder) const size_t pos = strlen(URL_str); raster_info_t rinfo; - rinfo.dataIsElevation = false; /* All bands are not elevation */ rinfo.fileId = finder->fileDict.add(filePath + fileName.substr(pos)); if(strcmp(bandName, "Fmask") == 0) { /* Use base class generic flags tag */ + rinfo.flagsBandNum = 1; rinfo.tag = FLAGS_TAG; if(parms->flags_file) { diff --git a/datasets/opendata/package/EsaWorldCover10meterRaster.h b/datasets/opendata/package/EsaWorldCover10meterRaster.h index c5d69a70..0c069a2b 100644 --- a/datasets/opendata/package/EsaWorldCover10meterRaster.h +++ b/datasets/opendata/package/EsaWorldCover10meterRaster.h @@ -63,8 +63,7 @@ class EsaWorldCover10meterRaster: public GeoRaster EsaWorldCover10meterRaster (lua_State* L, RequestFields* rqst_parms, const char* key): GeoRaster(L, rqst_parms, key, rqst_parms->geoFields(key)->asset.asset->getIndex(), - TimeLib::datetime2gps(2021, 06, 30, 0, 0, 0), /* Mid point for year data was collected */ - false /* Data is elevation */ ) {} + TimeLib::datetime2gps(2021, 06, 30, 0, 0, 0) /* Mid point for year data was collected */) {} private: diff --git a/datasets/opendata/package/MetaGlobalCanopy1meterRaster.h b/datasets/opendata/package/MetaGlobalCanopy1meterRaster.h index 4b129897..2f1b3b6c 100644 --- a/datasets/opendata/package/MetaGlobalCanopy1meterRaster.h +++ b/datasets/opendata/package/MetaGlobalCanopy1meterRaster.h @@ -63,8 +63,7 @@ class MetaGlobalCanopy1meterRaster: public GeoRaster MetaGlobalCanopy1meterRaster (lua_State* L, RequestFields* rqst_parms, const char* key): GeoRaster(L, rqst_parms, key, rqst_parms->geoFields(key)->asset.asset->getIndex(), - TimeLib::datetime2gps(2024, 4, 7, 0, 0, 0), - false /* Data is elevation */ ) {} + TimeLib::datetime2gps(2024, 4, 7, 0, 0, 0)) {} }; diff --git a/datasets/pgc/package/ArcticDemMosaicRaster.h b/datasets/pgc/package/ArcticDemMosaicRaster.h index 2822717b..5a772f7a 100644 --- a/datasets/pgc/package/ArcticDemMosaicRaster.h +++ b/datasets/pgc/package/ArcticDemMosaicRaster.h @@ -64,7 +64,8 @@ class ArcticDemMosaicRaster: public GeoRaster GeoRaster(L, rqst_parms, key, std::string(rqst_parms->geoFields(key)->asset.asset->getPath()).append("/").append(rqst_parms->geoFields(key)->asset.asset->getIndex()), TimeLib::datetime2gps(2023, 01, 18, 20, 23, 42), - true, /* Data is elevation */ + 1, /* elevationBandNum */ + GdalRaster::NO_BAND, /* maskBandNum */ &overrideTargetCRS) {} static OGRErr overrideTargetCRS(OGRSpatialReference& target) diff --git a/datasets/pgc/package/PgcDemStripsRaster.cpp b/datasets/pgc/package/PgcDemStripsRaster.cpp index c434c965..c60b9471 100644 --- a/datasets/pgc/package/PgcDemStripsRaster.cpp +++ b/datasets/pgc/package/PgcDemStripsRaster.cpp @@ -226,7 +226,7 @@ bool PgcDemStripsRaster::findRasters(raster_finder_t* finder) rasters_group_t* rgroup = new rasters_group_t; raster_info_t demRinfo; - demRinfo.dataIsElevation = true; + demRinfo.elevationBandNum = 1; demRinfo.tag = VALUE_TAG; demRinfo.fileId = finder->fileDict.add(fileName); @@ -245,7 +245,7 @@ bool PgcDemStripsRaster::findRasters(raster_finder_t* finder) if(!fileName.empty()) { raster_info_t flagsRinfo; - flagsRinfo.dataIsElevation = false; + flagsRinfo.flagsBandNum = 1; flagsRinfo.tag = FLAGS_TAG; flagsRinfo.fileId = finder->fileDict.add(fileName); rgroup->infovect.push_back(flagsRinfo); diff --git a/datasets/pgc/package/RemaDemMosaicRaster.h b/datasets/pgc/package/RemaDemMosaicRaster.h index 29e98d9d..0a6af124 100644 --- a/datasets/pgc/package/RemaDemMosaicRaster.h +++ b/datasets/pgc/package/RemaDemMosaicRaster.h @@ -65,7 +65,8 @@ class RemaDemMosaicRaster: public GeoRaster GeoRaster(L, rqst_parms, key, std::string(rqst_parms->geoFields(key)->asset.asset->getPath()).append("/").append(rqst_parms->geoFields(key)->asset.asset->getIndex()), TimeLib::datetime2gps(2023, 02, 24, 18, 51, 44), - true, /* Data is elevation */ + 1, /* elevationBandNum */ + GdalRaster::NO_BAND, /* flagsBandNum */ &overrideTargetCRS) {} static OGRErr overrideTargetCRS(OGRSpatialReference& target) diff --git a/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp b/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp index 3a65aef5..291eb7bd 100644 --- a/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp +++ b/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp @@ -121,7 +121,7 @@ bool Usgs3dep1meterDemRaster::findRasters(raster_finder_t* finder) const size_t pos = strlen(URL_str); raster_info_t rinfo; - rinfo.dataIsElevation = true; + rinfo.elevationBandNum = 1; rinfo.tag = VALUE_TAG; rinfo.fileId = finder->fileDict.add(filePath + fileName.substr(pos)); rgroup->infovect.push_back(rinfo); diff --git a/packages/geo/GdalRaster.cpp b/packages/geo/GdalRaster.cpp index 3a8ad6ac..141c5794 100644 --- a/packages/geo/GdalRaster.cpp +++ b/packages/geo/GdalRaster.cpp @@ -55,7 +55,10 @@ /*---------------------------------------------------------------------------- * Constructor *----------------------------------------------------------------------------*/ -GdalRaster::GdalRaster(const GeoFields* _parms, const std::string& _fileName, double _gpsTime, uint64_t _fileId, bool _dataIsElevation, overrideCRS_t cb, bbox_t* aoi_bbox_override): +GdalRaster::GdalRaster(const GeoFields* _parms, const std::string& _fileName, + double _gpsTime, uint64_t _fileId, + int _elevationBandNum, int _flagsBandNum, + overrideCRS_t cb, bbox_t* aoi_bbox_override): parms (_parms), gpsTime (_gpsTime), fileId (_fileId), @@ -63,7 +66,10 @@ GdalRaster::GdalRaster(const GeoFields* _parms, const std::string& _fileName, do overrideCRS(cb), fileName (_fileName), dset (NULL), - dataIsElevation(_dataIsElevation), + elevationBandNum(_elevationBandNum), + elevationBand(NULL), + flagsBandNum(_flagsBandNum), + flagsBand (NULL), xsize (0), ysize (0), cellSize (0), @@ -114,7 +120,7 @@ void GdalRaster::open(void) throw RunTimeException(CRITICAL, RTE_ERROR, "No bands found in raster: %s:", fileName.c_str()); } - /* Populate the mapping of band names to indices */ + /* Populate the mapping of band names to band numbers*/ for (int i = 1; i <= bandCount; ++i) { GDALRasterBand* band = dset->GetRasterBand(i); @@ -130,6 +136,20 @@ void GdalRaster::open(void) } } + /* Get elevation band */ + if(elevationBandNum > 0 && elevationBandNum <= bandCount) + { + elevationBand = dset->GetRasterBand(elevationBandNum); + CHECKPTR(elevationBand); + } + + /* Get flags band */ + if(flagsBandNum > 0 && flagsBandNum <= bandCount) + { + flagsBand = dset->GetRasterBand(flagsBandNum); + CHECKPTR(flagsBand); + } + /* Store information about raster */ xsize = dset->GetRasterXSize(); ysize = dset->GetRasterYSize(); @@ -176,6 +196,8 @@ void GdalRaster::open(void) OGRCoordinateTransformation::DestroyCT(transf); transf = NULL; bandMap.clear(); + elevationBand = NULL; + flagsBand = NULL; throw; } } @@ -199,10 +221,10 @@ RasterSample* GdalRaster::samplePOI(OGRPoint* poi, int bandNum) CHECKPTR(band); const double z = poi->getZ(); - // mlog(DEBUG, "Before transform x,y,z: (%.4lf, %.4lf, %.4lf)", poi->getX(), poi->getY(), poi->getZ()); + mlog(DEBUG, "Before transform x,y,z: (%.4lf, %.4lf, %.4lf)", poi->getX(), poi->getY(), poi->getZ()); if(poi->transform(transf) != OGRERR_NONE) throw RunTimeException(CRITICAL, RTE_ERROR, "Coordinates Transform failed for x,y,z (%lf, %lf, %lf)", poi->getX(), poi->getY(), poi->getZ()); - // mlog(DEBUG, "After transform x,y,z: (%.4lf, %.4lf, %.4lf)", poi->getX(), poi->getY(), poi->getZ()); + mlog(DEBUG, "After transform x,y,z: (%.4lf, %.4lf, %.4lf)", poi->getX(), poi->getY(), poi->getZ()); /* * Attempt to read raster only if it contains the point of interest. @@ -472,7 +494,8 @@ int GdalRaster::getBandNumber(const std::string& bandName) auto it = bandMap.find(bandName); if(it == bandMap.end()) { - mlog(WARNING, "Band name not found: %s", bandName.c_str()); + mlog(ERROR, "Band %s not found", bandName.c_str()); + return NO_BAND; } return it->second; } @@ -617,91 +640,93 @@ void GdalRaster::readPixel(const OGRPoint* poi, GDALRasterBand* band, RasterSamp /* Be carefull using offset based on the pixel data type */ switch(band->GetRasterDataType()) { - case GDT_Byte: - { - const uint8_t* p = static_cast(data); - sample->value = p[offset]; - } - break; + case GDT_Byte: + { + const uint8_t* p = static_cast(data); + sample->value = p[offset]; + } + break; - case GDT_Int8: - { - const int8_t* p = static_cast(data); - sample->value = p[offset]; - } - break; + case GDT_Int8: + { + const int8_t* p = static_cast(data); + sample->value = p[offset]; + } + break; - case GDT_UInt16: - { - const uint16_t* p = static_cast(data); - sample->value = p[offset]; - } - break; + case GDT_UInt16: + { + const uint16_t* p = static_cast(data); + sample->value = p[offset]; + } + break; - case GDT_Int16: - { - const int16_t* p = static_cast(data); - sample->value = p[offset]; - } - break; + case GDT_Int16: + { + const int16_t* p = static_cast(data); + sample->value = p[offset]; + } + break; - case GDT_UInt32: - { - const uint32_t* p = static_cast(data); - sample->value = p[offset]; - } - break; + case GDT_UInt32: + { + const uint32_t* p = static_cast(data); + sample->value = p[offset]; + } + break; - case GDT_Int32: - { - const int32_t* p = static_cast(data); - sample->value = p[offset]; - } - break; + case GDT_Int32: + { + const int32_t* p = static_cast(data); + sample->value = p[offset]; + } + break; - case GDT_Int64: - { - const int64_t* p = static_cast(data); - sample->value = p[offset]; - } - break; + case GDT_Int64: + { + const int64_t* p = static_cast(data); + sample->value = p[offset]; + } + break; - case GDT_UInt64: - { - const uint64_t* p = static_cast(data); - sample->value = p[offset]; - } - break; + case GDT_UInt64: + { + const uint64_t* p = static_cast(data); + sample->value = p[offset]; + } + break; - case GDT_Float32: - { - const float* p = static_cast(data); - sample->value = p[offset]; - } - break; + case GDT_Float32: + { + const float* p = static_cast(data); + sample->value = p[offset]; + } + break; - case GDT_Float64: - { - const double* p = static_cast(data); - sample->value = p[offset]; - } - break; + case GDT_Float64: + { + const double* p = static_cast(data); + sample->value = p[offset]; + } + break; - default: - /* - * Complex numbers are supported but not needed at this point. - */ - block->DropLock(); - throw RunTimeException(CRITICAL, RTE_ERROR, "Unsuported data type %d, in raster: %s:", band->GetRasterDataType(), fileName.c_str()); + default: + /* + * Complex numbers are supported but not needed at this point. + */ + block->DropLock(); + throw RunTimeException(CRITICAL, RTE_ERROR, "Unsuported data type %d, in raster: %s:", band->GetRasterDataType(), fileName.c_str()); } /* Done reading, release block lock */ block->DropLock(); - if(nodataCheck(sample, band) && dataIsElevation) + if(nodataCheck(sample, band) && band == elevationBand) { sample->value += sample->verticalShift; } + sample->bandName = band->GetDescription(); + // mlog(DEBUG, "Value: %.2lf, x: %u, y: %u, xblk: %u, yblk: %u, bcol: %u, brow: %u, offset: %u", // sample->value, x, y, xblk, yblk, _x, _y, offset); } @@ -766,10 +791,11 @@ void GdalRaster::resamplePixel(const OGRPoint* poi, GDALRasterBand* band, Raster if(validWindow) { readWithRetry(band, _x, _y, windowSize, windowSize, &sample->value, 1, 1, &args); - if(nodataCheck(sample, band) && dataIsElevation) + if(nodataCheck(sample, band) && band == elevationBand) { sample->value += sample->verticalShift; } + sample->bandName = band->GetDescription(); } else { @@ -832,7 +858,7 @@ void GdalRaster::computeZonalStats(const OGRPoint* poi, GDALRasterBand* band, Ra double value = samplesArray[_y*windowSize + _x]; if(value == nodata) continue; - if(dataIsElevation) + if(band == elevationBand) value += sample->verticalShift; const double x2 = _x + newx; /* Current pixel in buffer */ diff --git a/packages/geo/GdalRaster.h b/packages/geo/GdalRaster.h index cc7cbb97..d3819bf3 100644 --- a/packages/geo/GdalRaster.h +++ b/packages/geo/GdalRaster.h @@ -83,6 +83,7 @@ class GdalRaster static const int MAX_SAMPLING_RADIUS_IN_PIXELS = 50; static const int SLIDERULE_EPSG = 7912; + static const int NO_BAND = 0; /*-------------------------------------------------------------------- * Typedefs @@ -102,7 +103,11 @@ class GdalRaster * Methods *--------------------------------------------------------------------*/ - GdalRaster (const GeoFields* _parms, const std::string& _fileName, double _gpsTime, uint64_t _fileId, bool _dataIsElevation, overrideCRS_t cb, bbox_t* aoi_bbox_override=NULL); + GdalRaster (const GeoFields* _parms, const std::string& _fileName, + double _gpsTime, uint64_t _fileId, + int _elevationBandNum, int _flagsBandNum, + overrideCRS_t cb, bbox_t* aoi_bbox_override=NULL); + virtual ~GdalRaster (void); void open (void); RasterSample* samplePOI (OGRPoint* poi, int bandNum); @@ -114,7 +119,8 @@ class GdalRaster const bbox_t& getBbox (void) const { return bbox; } double getCellSize (void) const { return cellSize; } uint32_t getSSerror (void) const { return ssError; } - bool isElevation (void) const { return dataIsElevation; } + int getElevationBandNum (void) const { return elevationBandNum; } + int getFLagsBandNum (void) const { return flagsBandNum; } overrideCRS_t getOverrideCRS (void) const { return overrideCRS; } double getGpsTime (void) const { return gpsTime; } int getBandNumber (const std::string& bandName); @@ -147,7 +153,10 @@ class GdalRaster std::string fileName; GDALDataset *dset; - bool dataIsElevation; + int elevationBandNum; + GDALRasterBand* elevationBand; + int flagsBandNum; + GDALRasterBand* flagsBand; uint32_t xsize; uint32_t ysize; double cellSize; diff --git a/packages/geo/GeoIndexedRaster.cpp b/packages/geo/GeoIndexedRaster.cpp index 3c73d36e..b31710f3 100644 --- a/packages/geo/GeoIndexedRaster.cpp +++ b/packages/geo/GeoIndexedRaster.cpp @@ -955,17 +955,25 @@ void* GeoIndexedRaster::readerThread(void *param) cacheitem_t* entry = reader->entry; if(entry != NULL) { + GdalRaster* raster = entry->raster; + + /* Open raster so we can get inner bands from it */ + raster->open(); + std::vector bands; - reader->obj->getInnerBands(entry->raster, bands); + reader->obj->getInnerBands(raster, bands); - GdalRaster* raster = entry->raster; if(GdalRaster::ispoint(reader->geo)) { /* Sample raster bands */ for(const int bandNum : bands) { - RasterSample* sample = raster->samplePOI(dynamic_cast(reader->geo), bandNum); + /* Use local copy of point, it will be projected in samplePOI. We do not want project it again */ + OGRPoint point(*(dynamic_cast(reader->geo))); + + RasterSample* sample = raster->samplePOI(&point, bandNum); entry->bandSample.push_back(sample); + mlog(DEBUG, "Band: %d, %s", bandNum, sample->toString().c_str()); } } else if(GdalRaster::ispoly(reader->geo)) @@ -973,6 +981,7 @@ void* GeoIndexedRaster::readerThread(void *param) /* Subset raster bands */ for(const int bandNum : bands) { + /* No need to use local copy of polygon, subsetAOI will use it's envelope and not project it */ RasterSubset* subset = raster->subsetAOI(dynamic_cast(reader->geo), bandNum); if(subset) { @@ -987,7 +996,8 @@ void* GeoIndexedRaster::readerThread(void *param) reader->obj->samplerKey, subset->rasterName, raster->getGpsTime(), - raster->isElevation(), + raster->getElevationBandNum(), + raster->getFLagsBandNum(), raster->getOverrideCRS()); entry->bandSubset.push_back(subset); @@ -1035,7 +1045,8 @@ void* GeoIndexedRaster::batchReaderThread(void *param) breader->obj->fileDict.get(ur->rinfo->fileId), 0, /* Sample collecting code will set it to group's gpsTime */ ur->rinfo->fileId, - ur->rinfo->dataIsElevation, + ur->rinfo->elevationBandNum, + ur->rinfo->flagsBandNum, breader->obj->crscb); std::vector bands; @@ -1293,10 +1304,14 @@ bool GeoIndexedRaster::updateCache(uint32_t& rasters2sample, const GroupOrdering note use of bbox in construcutor - it limits area of interest to the extent of vector index file */ item = new cacheitem_t; - item->raster = new GdalRaster(parms, key, + item->raster = new GdalRaster(parms, + key, static_cast(rgroup->gpsTime / 1000), rinfo.fileId, - rinfo.dataIsElevation, crscb, &bbox); + rinfo.elevationBandNum, + rinfo.flagsBandNum, + crscb, + &bbox); const bool status = cache.add(key, item); assert(status); (void)status; // cannot fail; prevents linter warnings } diff --git a/packages/geo/GeoIndexedRaster.h b/packages/geo/GeoIndexedRaster.h index 3dd17e85..3423c37f 100644 --- a/packages/geo/GeoIndexedRaster.h +++ b/packages/geo/GeoIndexedRaster.h @@ -90,12 +90,13 @@ class GeoIndexedRaster: public RasterObject /** Raster information needed for sampling */ typedef struct RasterInfo { - bool dataIsElevation; - std::string tag; // value, fmask, etc - uint64_t fileId; // file dictionary id - UniqueRaster* uraster; // Pointer to the unique raster which contains the sample for this raster + int elevationBandNum; // band number for elevation + int flagsBandNum; // band number for fmask + std::string tag; // value, fmask, etc + uint64_t fileId; // file dictionary id + UniqueRaster* uraster; // Pointer to the unique raster which contains the sample for this raster - RasterInfo(void): dataIsElevation(false), uraster(NULL) {} + RasterInfo(void): elevationBandNum(GdalRaster::NO_BAND), flagsBandNum(GdalRaster::NO_BAND), uraster(NULL) {} } raster_info_t; /* Group of rasters belonging to the same geojson stac catalog feature */ diff --git a/packages/geo/GeoJsonRaster.cpp b/packages/geo/GeoJsonRaster.cpp index e17a5380..decfcd72 100644 --- a/packages/geo/GeoJsonRaster.cpp +++ b/packages/geo/GeoJsonRaster.cpp @@ -155,7 +155,7 @@ GeoJsonRaster::~GeoJsonRaster(void) * Constructor *----------------------------------------------------------------------------*/ GeoJsonRaster::GeoJsonRaster(lua_State* L, RequestFields* rqst_parms, const char* key, const char* _geojstr, double _cellsize): - GeoRaster(L, rqst_parms, key, std::string("/vsimem/" + GdalRaster::getUUID() + ".tif"), TimeLib::gpstime(), false /* not elevation*/), + GeoRaster(L, rqst_parms, key, std::string("/vsimem/" + GdalRaster::getUUID() + ".tif"), TimeLib::gpstime()), data(NULL), cellsize(_cellsize), cols(0), diff --git a/packages/geo/GeoRaster.cpp b/packages/geo/GeoRaster.cpp index 6aa1d778..ffd5f028 100644 --- a/packages/geo/GeoRaster.cpp +++ b/packages/geo/GeoRaster.cpp @@ -56,9 +56,10 @@ void GeoRaster::deinit (void) /*---------------------------------------------------------------------------- * Constructor *----------------------------------------------------------------------------*/ -GeoRaster::GeoRaster(lua_State *L, RequestFields* rqst_parms, const char* key, const std::string& _fileName, double _gpsTime, bool dataIsElevation, GdalRaster::overrideCRS_t cb): +GeoRaster::GeoRaster(lua_State *L, RequestFields* rqst_parms, const char* key, const std::string& _fileName, + double _gpsTime, int elevationBandNum, int flagsBandNum, GdalRaster::overrideCRS_t cb): RasterObject(L, rqst_parms, key), - raster(parms, _fileName, _gpsTime, fileDict.add(_fileName, true), dataIsElevation, cb) + raster(parms, _fileName, _gpsTime, fileDict.add(_fileName, true), elevationBandNum, flagsBandNum, cb) { /* Add Lua Functions */ LuaEngine::setAttrFunc(L, "dim", luaDimensions); @@ -141,7 +142,8 @@ uint32_t GeoRaster::getSubsets(const MathLib::extent_t& extent, int64_t gps, Lis samplerKey, subset->rasterName, raster.getGpsTime(), - raster.isElevation(), + raster.getElevationBandNum(), + raster.getFLagsBandNum(), raster.getOverrideCRS()); slist.add(subset); diff --git a/packages/geo/GeoRaster.h b/packages/geo/GeoRaster.h index a6a930c2..1b0f3a92 100644 --- a/packages/geo/GeoRaster.h +++ b/packages/geo/GeoRaster.h @@ -61,7 +61,8 @@ class GeoRaster: public RasterObject static void init (void); static void deinit (void); - GeoRaster (lua_State* L, RequestFields* rqst_parms, const char* key, const std::string& _fileName, double _gpsTime, bool dataIsElevation, GdalRaster::overrideCRS_t cb=NULL); + GeoRaster (lua_State* L, RequestFields* rqst_parms, const char* key, const std::string& _fileName, double _gpsTime, + int elevationBandNum=GdalRaster::NO_BAND, int flagsBandNum=GdalRaster::NO_BAND, GdalRaster::overrideCRS_t cb=NULL); ~GeoRaster (void) override; uint32_t getSamples (const MathLib::point_3d_t& point, int64_t gps, List& slist, void* param=NULL) final; uint32_t getSubsets (const MathLib::extent_t& extent, int64_t gps, List& slist, void* param=NULL) final; diff --git a/packages/geo/GeoUserRaster.cpp b/packages/geo/GeoUserRaster.cpp index 976f358e..2f864309 100644 --- a/packages/geo/GeoUserRaster.cpp +++ b/packages/geo/GeoUserRaster.cpp @@ -105,8 +105,12 @@ int GeoUserRaster::luaCreate (lua_State* L) if(rasterlength > maxSize) throw RunTimeException(CRITICAL, RTE_ERROR, "User raster too big, size is: %lu, max allowed: %u", rasterlength, maxSize); + /* If raster has elevation assume it is in the first band */ + const int elevationBandNum = iselevation ? 1 : GdalRaster::NO_BAND; + const int flagsBandNum = GdalRaster::NO_BAND; + /* Create GeoUserRaster */ - return createLuaObject(L, new GeoUserRaster(L, rqst_parms, GeoFields::DEFAULT_KEY, tiff.c_str(), tiff.size(), gps, iselevation)); + return createLuaObject(L, new GeoUserRaster(L, rqst_parms, GeoFields::DEFAULT_KEY, tiff.c_str(), tiff.size(), gps, elevationBandNum, flagsBandNum)); } catch(const RunTimeException& e) { @@ -132,8 +136,10 @@ GeoUserRaster::~GeoUserRaster(void) /*---------------------------------------------------------------------------- * Constructor *----------------------------------------------------------------------------*/ -GeoUserRaster::GeoUserRaster(lua_State *L, RequestFields* rqst_parms, const char* key, const char *file, long filelength, double gps, bool iselevation): - GeoRaster(L, rqst_parms, key, std::string("/vsimem/userraster/" + GdalRaster::getUUID() + ".tif"), gps, iselevation ), +GeoUserRaster::GeoUserRaster(lua_State *L, RequestFields* rqst_parms, const char* key, + const char *file, long filelength, double gps, + int elevationBandNum, int flagsBandNum) : + GeoRaster(L, rqst_parms, key, std::string("/vsimem/userraster/" + GdalRaster::getUUID() + ".tif"), gps, elevationBandNum, flagsBandNum), data(NULL) { if(file == NULL) diff --git a/packages/geo/GeoUserRaster.h b/packages/geo/GeoUserRaster.h index 85a26d3a..7b3fd785 100644 --- a/packages/geo/GeoUserRaster.h +++ b/packages/geo/GeoUserRaster.h @@ -71,7 +71,9 @@ class GeoUserRaster: public GeoRaster * Methods *--------------------------------------------------------------------*/ - GeoUserRaster(lua_State* L, RequestFields* rqst_parms, const char* key, const char* file, long filelength, double gps, bool iselevation); + GeoUserRaster(lua_State* L, RequestFields* rqst_parms, const char* key, + const char* file, long filelength, double gps, + int elevationBandNum, int flagsBandNum); private: diff --git a/packages/geo/RasterObject.cpp b/packages/geo/RasterObject.cpp index 8960b370..7542fc60 100644 --- a/packages/geo/RasterObject.cpp +++ b/packages/geo/RasterObject.cpp @@ -437,6 +437,7 @@ int RasterObject::luaSamples(lua_State *L) LuaEngine::setAttrInt(L, "fileid", sample->fileId); LuaEngine::setAttrNum(L, "time", sample->time); LuaEngine::setAttrNum(L, "value", sample->value); + LuaEngine::setAttrStr(L, "band", sample->bandName.c_str()); lua_rawseti(L, -2, i+1); } } else mlog(DEBUG, "No samples read for (%.2lf, %.2lf)", lon, lat); From 36e135e4f63ff233c8d862ce731ed1fee5cada74 Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Wed, 13 Nov 2024 22:50:52 +0000 Subject: [PATCH 06/17] fixed crash in debug msg --- packages/geo/GeoIndexedRaster.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/geo/GeoIndexedRaster.cpp b/packages/geo/GeoIndexedRaster.cpp index b31710f3..82d37cc7 100644 --- a/packages/geo/GeoIndexedRaster.cpp +++ b/packages/geo/GeoIndexedRaster.cpp @@ -973,7 +973,7 @@ void* GeoIndexedRaster::readerThread(void *param) RasterSample* sample = raster->samplePOI(&point, bandNum); entry->bandSample.push_back(sample); - mlog(DEBUG, "Band: %d, %s", bandNum, sample->toString().c_str()); + mlog(DEBUG, "Band: %d, %s", bandNum, sample ? sample->toString().c_str() : "NULL"); } } else if(GdalRaster::ispoly(reader->geo)) From a0f8ee5d9f1518434472b323c0a18b6dba4e5f5a Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Thu, 14 Nov 2024 16:40:16 +0000 Subject: [PATCH 07/17] added unit test for bluetopo, all working now --- .../bluetopo/package/BlueTopoBathyRaster.cpp | 39 ++++++++++ .../bluetopo/package/BlueTopoBathyRaster.h | 4 +- .../bluetopo/selftests/bluetopo_reader.lua | 8 +- .../bluetopo/selftests/plugin_unittest.lua | 39 ++++++++++ datasets/landsat/package/LandsatHlsRaster.cpp | 77 +++++++++++++++---- datasets/landsat/package/LandsatHlsRaster.h | 13 ++-- datasets/landsat/selftests/landsat_reader.lua | 2 +- packages/geo/GdalRaster.cpp | 2 +- packages/geo/GeoIndexedRaster.cpp | 49 +++++++++--- packages/geo/UT_RasterSample.cpp | 11 +++ scripts/selftests/test_runner.lua | 1 + 11 files changed, 204 insertions(+), 41 deletions(-) create mode 100644 datasets/bluetopo/selftests/plugin_unittest.lua diff --git a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp index e7f1a1b1..a836696f 100644 --- a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp @@ -35,6 +35,10 @@ #include "BlueTopoBathyRaster.h" +/* Valid Bands for BlueTopo Bathy Raster */ +const char* BlueTopoBathyRaster::validBands[] = {"Elevation", "Uncertainty", "Contributor"}; +#define validBandsCnt (sizeof (validBands) / sizeof (const char *)) + /****************************************************************************** * PROTECTED METHODS ******************************************************************************/ @@ -47,6 +51,11 @@ BlueTopoBathyRaster::BlueTopoBathyRaster(lua_State* L, RequestFields* rqst_parms filePath(parms->asset.asset->getPath()), indexBucket(parms->asset.asset->getIndex()) { + if(!validateBandNames()) + { + throw RunTimeException(DEBUG, RTE_ERROR, "Invalid band name specified"); + } + const std::string bucketPath = filePath + indexBucket; if(!findIndexFileInS3Bucket(bucketPath)) { @@ -191,6 +200,36 @@ double BlueTopoBathyRaster::getGmtDate(const OGRFeature* feature, const char* fi * PRIVATE METHODS ******************************************************************************/ +/*---------------------------------------------------------------------------- + * validateBandNames + *----------------------------------------------------------------------------*/ +bool BlueTopoBathyRaster::validateBandNames(void) +{ + for (int i = 0; i < parms->bands.length(); i++) + { + bool valid = false; + + for(size_t j = 0; j < validBandsCnt; j++) + { + /* Raster band names are case insensitive */ + if(strcmp(parms->bands[i].c_str(), validBands[j]) == 0) + { + valid = true; + break; + } + } + + if(!valid) + { + mlog(ERROR, "Invalid band name: %s", parms->bands[i].c_str()); + return false; + } + } + + return true; +} + + /*---------------------------------------------------------------------------- * findIndexFileInS3Bucket *----------------------------------------------------------------------------*/ diff --git a/datasets/bluetopo/package/BlueTopoBathyRaster.h b/datasets/bluetopo/package/BlueTopoBathyRaster.h index d25bf8e3..6a81daea 100644 --- a/datasets/bluetopo/package/BlueTopoBathyRaster.h +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.h @@ -49,6 +49,7 @@ class BlueTopoBathyRaster: public GeoIndexedRaster /*-------------------------------------------------------------------- * Constants *--------------------------------------------------------------------*/ + static const char* validBands[]; /*-------------------------------------------------------------------- * Typedefs @@ -85,7 +86,8 @@ class BlueTopoBathyRaster: public GeoIndexedRaster * Methods *--------------------------------------------------------------------*/ - bool findIndexFileInS3Bucket(const std::string& bucketPath); + bool validateBandNames (void); + bool findIndexFileInS3Bucket (const std::string& bucketPath); /*-------------------------------------------------------------------- * Data diff --git a/datasets/bluetopo/selftests/bluetopo_reader.lua b/datasets/bluetopo/selftests/bluetopo_reader.lua index 04650864..846fced5 100644 --- a/datasets/bluetopo/selftests/bluetopo_reader.lua +++ b/datasets/bluetopo/selftests/bluetopo_reader.lua @@ -19,10 +19,8 @@ local expElevation = { -2.02, -14.03, -4.28, -17.18} local expUncertainty = { 1.00, 4.37, 0.34, 1.32} local expContributor = { 31156, 0, 24955, 45641} --- Tolerance for depth values is needed because of the different versions of PROJLIB --- The values are different in the different versions of the library. --- We add 1 meter of depth tolerance for the test to pass -local elevation_tolerance = 0.1; + +local elevation_tolerance = 0.01; print(string.format("\n--------------------------------\nTest: BlueTopo Correct Values\n--------------------------------")) local dem = geo.raster(geo.parms({ asset = "bluetopo-bathy", algorithm = "NearestNeighbour", bands = {"Elevation", "Uncertainty", "Contributor"}, sort_by_index = true })) @@ -42,7 +40,7 @@ for j, lon in ipairs(lons) do band = v["band"] value = v["value"] fname = v["file"] - print(string.format("(%02d) (%6.2f, %6.2f) %12s %10.2f %s", k, lon, lat, band, value, fname)) + print(string.format("(%6.2f, %6.2f) Band: %11s %8.2f %s", lon, lat, band, value, fname)) if band == "Elevation" then runner.check(math.abs(value - expElevation[j]) < elevation_tolerance, string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) diff --git a/datasets/bluetopo/selftests/plugin_unittest.lua b/datasets/bluetopo/selftests/plugin_unittest.lua new file mode 100644 index 00000000..0e55c887 --- /dev/null +++ b/datasets/bluetopo/selftests/plugin_unittest.lua @@ -0,0 +1,39 @@ +local runner = require("test_executive") +local console = require("console") +local asset = require("asset") +local assets = asset.loaddir() + +-- console.monitor:config(core.LOG, core.DEBUG) +-- sys.setlvl(core.LOG, core.DEBUG) + +-- Check If Present -- +if not core.UNITTEST then return end + +-- Setup -- + +-- Unit Test -- + +lon = -80.87 +lat = 32.06 + +lon_incr = 0.0001 +lat_incr = 0.0001 +pointCount = 10 + +print('\n------------------\nTest: RasterSampler::bluetopo-bathy\n------------------') +local dem = geo.raster(geo.parms({ asset = "bluetopo-bathy", algorithm = "NearestNeighbour", bands = {"Elevation", "Uncertainty", "Contributor"}, sort_by_index = true })) +runner.check(dem ~= nil) + +ut = geo.ut_sample(dem) +runner.check(ut ~= nil) +status = ut:test(lon, lat, lon_incr, lat_incr, pointCount) +runner.check(status, "Failed sampling test") +ut = nil + + +-- Clean Up -- + +-- Report Results -- + +runner.report() + diff --git a/datasets/landsat/package/LandsatHlsRaster.cpp b/datasets/landsat/package/LandsatHlsRaster.cpp index 5ee995e1..c4977204 100644 --- a/datasets/landsat/package/LandsatHlsRaster.cpp +++ b/datasets/landsat/package/LandsatHlsRaster.cpp @@ -78,6 +78,13 @@ LandsatHlsRaster::LandsatHlsRaster(lua_State *L, RequestFields* rqst_parms, cons { ndsi = ndvi = ndwi = false; + if(!validateBandNames()) + { + VSIUnlink(indexFile.c_str()); + throw RunTimeException(DEBUG, RTE_ERROR, "Invalid band names specified"); + } + + /* Create in memory index file (geojson) */ VSILFILE* fp = VSIFileFromMemBuffer(indexFile.c_str(), const_cast(reinterpret_cast(parms->catalog.value.c_str())), // source bytes @@ -90,20 +97,18 @@ LandsatHlsRaster::LandsatHlsRaster(lua_State *L, RequestFields* rqst_parms, cons for(int i = 0; i < parms->bands.length(); i++) { const string& name = parms->bands[i]; - if( isValidL8Band(name.c_str()) || isValidS2Band(name.c_str()) || isValidAlgoName(name.c_str())) - { - /* Add band to dictionary of bands but don't override if already there */ - auto it = bandsDict.find(name); - if(it == bandsDict.end()) - bandsDict[name] = true; - if(strcasecmp(parms->bands[i].c_str(), "NDSI") == 0) ndsi = true; - if(strcasecmp(parms->bands[i].c_str(), "NDVI") == 0) ndvi = true; - if(strcasecmp(parms->bands[i].c_str(), "NDWI") == 0) ndwi = true; - } + /* Add band to dictionary of bands but don't override if already there */ + auto it = bandsDict.find(name); + if(it == bandsDict.end()) + bandsDict[name] = true; + + if(strcasecmp(name.c_str(), "NDSI") == 0) ndsi = true; + if(strcasecmp(name.c_str(), "NDVI") == 0) ndvi = true; + if(strcasecmp(name.c_str(), "NDWI") == 0) ndwi = true; } - /* If user specified only algortithm(s), add needed bands to dictionary of bands */ + /* If user specified algortithm(s), add needed bands to dictionary of bands */ if(ndsi || ndvi || ndwi) { for (unsigned i=0; ifirst.c_str(); /* skip algo names (NDIS, etc) */ - if(isValidAlgoName(bandName)) + if(validAlgoName(bandName)) continue; const char* fname = feature->GetFieldAsString(bandName); @@ -225,6 +230,28 @@ bool LandsatHlsRaster::findRasters(raster_finder_t* finder) * PRIVATE METHODS ******************************************************************************/ +/*---------------------------------------------------------------------------- + * validateBandNames + *----------------------------------------------------------------------------*/ +bool LandsatHlsRaster::validateBandNames(void) +{ + bool valid = false; + + for (int i = 0; i < parms->bands.length(); i++) + { + const std::string& bandName = parms->bands[i]; + valid = validL8Band(bandName.c_str()) || validS2Band(bandName.c_str()) || validAlgoName(bandName.c_str()); + + /* User specified invalid band name */ + if(!valid) mlog(ERROR, "Invalid band name: %s", bandName.c_str()); + } + + return valid; +} + +/*---------------------------------------------------------------------------- + * validateBand + *----------------------------------------------------------------------------*/ bool LandsatHlsRaster::validateBand(band_type_t type, const char* bandName) { const char** tags; @@ -294,6 +321,9 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr /* Landsat rasters use only the first inner band */ const int INNER_BAND_INDX = 0; + /* Samples to be returned to the user */ + std::vector sampleVect; + /* Collect samples for all rasters */ if(mode == SERIAL) { @@ -319,7 +349,7 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr if(returnBandSample) { sample->bandName = bandName; - slist->add(sample); + sampleVect.push_back(sample); item->bandSample[INNER_BAND_INDX] = NULL; } } @@ -387,7 +417,7 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr /* Set flags for this sample, add it to the list */ s->flags = flags; - slist->add(s); + sampleVect.push_back(sample); errors |= ps.ssErrors; } } @@ -430,7 +460,7 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr /* Set band name for this sample */ sample->bandName = "NDSI"; - slist->add(sample); + sampleVect.push_back(sample); } if(ndvi) @@ -442,7 +472,7 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr /* Set band name for this sample */ sample->bandName = "NDVI"; - slist->add(sample); + sampleVect.push_back(sample); } if(ndwi) @@ -454,7 +484,20 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr /* Set band name for this sample */ sample->bandName = "NDWI"; - slist->add(sample); + sampleVect.push_back(sample); + } + + /* Add samples to the slist in the order of bands specified by the user in parms */ + for(int i = 0; i < parms->bands.length(); i++) + { + for(RasterSample* sample : sampleVect) + { + if(strcasecmp(parms->bands[i].c_str(), sample->bandName.c_str()) == 0) + { + slist->add(sample); + break; + } + } } return errors; diff --git a/datasets/landsat/package/LandsatHlsRaster.h b/datasets/landsat/package/LandsatHlsRaster.h index 5a55b522..2b27e2be 100644 --- a/datasets/landsat/package/LandsatHlsRaster.h +++ b/datasets/landsat/package/LandsatHlsRaster.h @@ -109,12 +109,13 @@ class LandsatHlsRaster: public GeoIndexedRaster * Methods *--------------------------------------------------------------------*/ - static bool validateBand (band_type_t type, const char* bandName); + bool validateBandNames (void); + static bool validateBand (band_type_t type, const char* bandName); - static bool isValidL8Band (const char* bandName) {return validateBand(LANDSAT8, bandName);} - static bool isValidS2Band (const char* bandName) {return validateBand(SENTINEL2,bandName);} - static bool isValidAlgoBand (const char* bandName) {return validateBand(ALGOBAND, bandName);} - static bool isValidAlgoName (const char* bandName) {return validateBand(ALGONAME, bandName);} + static bool validL8Band (const char* bandName) {return validateBand(LANDSAT8, bandName);} + static bool validS2Band (const char* bandName) {return validateBand(SENTINEL2,bandName);} + static bool validAlgoBand (const char* bandName) {return validateBand(ALGOBAND, bandName);} + static bool validAlgoName (const char* bandName) {return validateBand(ALGONAME, bandName);} uint32_t _getGroupSamples(sample_mode_t mode, const rasters_group_t* rgroup, List* slist, uint32_t flags, uint32_t pointIndx=0); @@ -125,7 +126,7 @@ class LandsatHlsRaster: public GeoIndexedRaster std::string filePath; std::string indexFile; - std::unordered_map bandsDict; + std::unordered_map bandsDict; /* Bands (rasters) to sample */ bool ndsi; bool ndvi; diff --git a/datasets/landsat/selftests/landsat_reader.lua b/datasets/landsat/selftests/landsat_reader.lua index 6b1a7414..bb5fa62e 100644 --- a/datasets/landsat/selftests/landsat_reader.lua +++ b/datasets/landsat/selftests/landsat_reader.lua @@ -270,7 +270,7 @@ expResults = {{ 523.0, 1293577339.0, 0xc2, '/vsis3/lp-prod-protected/HLSS30.020 print(string.format("\n-------------------------------------------------\nLandsat Sample test (BO3 and Fmask raster)\n-------------------------------------------------")) -dem = geo.raster(geo.parms({ asset = demType, algorithm = "NearestNeighbour", radius = 0, with_flags=true, bands = {"Fmask","B03"}, catalog = contents, sort_by_index = true })) +dem = geo.raster(geo.parms({ asset = demType, algorithm = "NearestNeighbour", radius = 0, with_flags=true, bands = {"B03", "Fmask"}, catalog = contents, sort_by_index = true })) fmaskCnt = 0 sampleCnt = 0 b03 = 0 diff --git a/packages/geo/GdalRaster.cpp b/packages/geo/GdalRaster.cpp index 141c5794..47dc72a9 100644 --- a/packages/geo/GdalRaster.cpp +++ b/packages/geo/GdalRaster.cpp @@ -494,7 +494,7 @@ int GdalRaster::getBandNumber(const std::string& bandName) auto it = bandMap.find(bandName); if(it == bandMap.end()) { - mlog(ERROR, "Band %s not found", bandName.c_str()); + mlog(ERROR, "Band \"%s\" not found", bandName.c_str()); return NO_BAND; } return it->second; diff --git a/packages/geo/GeoIndexedRaster.cpp b/packages/geo/GeoIndexedRaster.cpp index 82d37cc7..39fc1dfd 100644 --- a/packages/geo/GeoIndexedRaster.cpp +++ b/packages/geo/GeoIndexedRaster.cpp @@ -966,14 +966,26 @@ void* GeoIndexedRaster::readerThread(void *param) if(GdalRaster::ispoint(reader->geo)) { /* Sample raster bands */ - for(const int bandNum : bands) + const bool oneBand = bands.size() == 1; + if(oneBand) { - /* Use local copy of point, it will be projected in samplePOI. We do not want project it again */ - OGRPoint point(*(dynamic_cast(reader->geo))); - - RasterSample* sample = raster->samplePOI(&point, bandNum); + OGRPoint* p = (dynamic_cast(reader->geo)); + RasterSample* sample = raster->samplePOI(p, bands[0]); entry->bandSample.push_back(sample); - mlog(DEBUG, "Band: %d, %s", bandNum, sample ? sample->toString().c_str() : "NULL"); + } + else + { + /* Multiple bands */ + for(const int bandNum : bands) + { + /* Use local copy of point, it will be projected in samplePOI. We do not want to project it again */ + OGRPoint* p = (dynamic_cast(reader->geo)); + OGRPoint point(*p); + + RasterSample* sample = raster->samplePOI(&point, bandNum); + entry->bandSample.push_back(sample); + mlog(DEBUG, "Band: %d, %s", bandNum, sample ? sample->toString().c_str() : "NULL"); + } } } else if(GdalRaster::ispoly(reader->geo)) @@ -1049,19 +1061,36 @@ void* GeoIndexedRaster::batchReaderThread(void *param) ur->rinfo->flagsBandNum, breader->obj->crscb); + /* Open raster so we can get inner bands from it */ + raster->open(); + std::vector bands; breader->obj->getInnerBands(raster, bands); /* Sample all points for this raster */ for(point_sample_t& ps : ur->pointSamples) { - /* Sample all bands for this point */ - for(const int bandNum : bands) + /* Sample raster bands */ + const bool oneBand = bands.size() == 1; + if(oneBand) { - RasterSample* sample = raster->samplePOI(&ps.point, bandNum); + RasterSample* sample = raster->samplePOI(&ps.point, bands[0]); ps.bandSample.push_back(sample); ps.bandSampleReturned.emplace_back(std::make_unique>(false)); - ps.ssErrors |= raster->getSSerror(); + } + else + { + /* Multiple bands */ + for(const int bandNum : bands) + { + /* Use local copy of point, it will be projected in samplePOI. We do not want to project it again */ + OGRPoint point(ps.point); + + RasterSample* sample = raster->samplePOI(&point, bandNum); + ps.bandSample.push_back(sample); + ps.bandSampleReturned.emplace_back(std::make_unique>(false)); + ps.ssErrors |= raster->getSSerror(); + } } } diff --git a/packages/geo/UT_RasterSample.cpp b/packages/geo/UT_RasterSample.cpp index 2b393e05..89874282 100644 --- a/packages/geo/UT_RasterSample.cpp +++ b/packages/geo/UT_RasterSample.cpp @@ -367,6 +367,17 @@ int UT_RasterSample::luaSampleTest(lua_State* L) errors++; } + const std::string& serialBand = serial->bandName; + const std::string& batchBand = batch->bandName; + + if (serialBand != batchBand) + { + print2term("Bands differ:\n"); + print2term("Serial: %s\n", serialBand.c_str()); + print2term("Batch: %s\n", batchBand.c_str()); + errors++; + } + /* Compare as whole seconds */ if(static_cast(serial->time) != static_cast(batch->time)) { diff --git a/scripts/selftests/test_runner.lua b/scripts/selftests/test_runner.lua index 1fe59b5b..ce413e4c 100644 --- a/scripts/selftests/test_runner.lua +++ b/scripts/selftests/test_runner.lua @@ -117,6 +117,7 @@ end -- Run BLUETOPO Plugin Self Tests if __bluetopo__ and incloud then local bluetopo_td = td .. "../../datasets/bluetopo/selftests/" + runner.script(bluetopo_td.. "plugin_unittest.lua") runner.script(bluetopo_td.. "bluetopo_reader.lua") end From 6fdfff8aed8e6e3b25e0340c376946ccebfc84ef Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Thu, 14 Nov 2024 19:12:58 +0000 Subject: [PATCH 08/17] added bluetopo python test --- clients/python/tests/test_bluetopo.py | 33 ++++ clients/python/tests/test_landsat.py | 5 +- .../bluetopo/package/BlueTopoBathyRaster.cpp | 7 +- datasets/landsat/package/LandsatHlsRaster.cpp | 17 +- packages/geo/GeoIndexedRaster.cpp | 183 ++++++++++-------- 5 files changed, 156 insertions(+), 89 deletions(-) create mode 100644 clients/python/tests/test_bluetopo.py diff --git a/clients/python/tests/test_bluetopo.py b/clients/python/tests/test_bluetopo.py new file mode 100644 index 00000000..9b1595a6 --- /dev/null +++ b/clients/python/tests/test_bluetopo.py @@ -0,0 +1,33 @@ +"""Tests for sliderule bluetopo raster support.""" + +import pytest +from pathlib import Path +import sliderule + +TESTDIR = Path(__file__).parent + +sigma = 1.0e-9 + +lon = -80.87 +lat = 32.06 + +file = '/vsis3/noaa-ocs-nationalbathymetry-pds/BlueTopo/BH4SX59B/BlueTopo_BH4SX59B_20241021.tiff' + +expElevation = -2.02 +expUncertainty = 1.00 +expContributor = 31156 + +@pytest.mark.network +class TestBlueTopo: + def test_samples(self, init): + rqst = {"samples": {"asset": "bluetopo-bathy", "bands": ["Elevation", "Uncertainty", "Contributor"]}, "coordinates": [[lon,lat]]} + rsps = sliderule.source("samples", rqst) + assert init + print(rsps) + assert rsps["samples"][0][0]["file"] == file + assert rsps["samples"][0][0]["band"] == "Elevation" + assert rsps['samples'][0][0]['value'] == pytest.approx(expElevation, rel=1e-3) + assert rsps["samples"][0][1]["band"] == "Uncertainty" + assert rsps['samples'][0][1]['value'] == pytest.approx(expUncertainty, rel=1e-3) + assert rsps["samples"][0][2]["band"] == "Contributor" + assert rsps['samples'][0][2]['value'] == expContributor diff --git a/clients/python/tests/test_landsat.py b/clients/python/tests/test_landsat.py index a374f005..819d3a28 100644 --- a/clients/python/tests/test_landsat.py +++ b/clients/python/tests/test_landsat.py @@ -18,10 +18,13 @@ def test_samples(self, init): {"lon": -177.0000000001, "lat": 49.0000000001}, {"lon": -177.0000000001, "lat": 51.0000000001} ] catalog = earthdata.stac(short_name="HLS", polygon=polygon, time_start=time_start, time_end=time_end, as_str=True) - rqst = {"samples": {"asset": "landsat-hls", "catalog": catalog, "bands": ["B02"]}, "coordinates": [[-178.0, 50.7]]} + rqst = {"samples": {"asset": "landsat-hls", "catalog": catalog, "closest_time": "2020-01-05T00:00:00Z", "bands": ["NDVI"]}, "coordinates": [[-179.0, 51.0]]} rsps = sliderule.source("samples", rqst) assert init + print(rsps) assert len(rsps) > 0 + assert rsps['samples'][0][0]['band'] == "NDVI" + assert rsps['samples'][0][0]['value'] == pytest.approx(-0.259439707674, rel=1e-6) def test_cmr_stac(self, init): time_start = "2000-01-01T00:00:00Z" diff --git a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp index a836696f..92313399 100644 --- a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp @@ -205,10 +205,15 @@ double BlueTopoBathyRaster::getGmtDate(const OGRFeature* feature, const char* fi *----------------------------------------------------------------------------*/ bool BlueTopoBathyRaster::validateBandNames(void) { + if(parms->bands.length() == 0) + { + mlog(ERROR, "No bands specified"); + return false; + } + for (int i = 0; i < parms->bands.length(); i++) { bool valid = false; - for(size_t j = 0; j < validBandsCnt; j++) { /* Raster band names are case insensitive */ diff --git a/datasets/landsat/package/LandsatHlsRaster.cpp b/datasets/landsat/package/LandsatHlsRaster.cpp index c4977204..bae88c57 100644 --- a/datasets/landsat/package/LandsatHlsRaster.cpp +++ b/datasets/landsat/package/LandsatHlsRaster.cpp @@ -235,18 +235,23 @@ bool LandsatHlsRaster::findRasters(raster_finder_t* finder) *----------------------------------------------------------------------------*/ bool LandsatHlsRaster::validateBandNames(void) { - bool valid = false; + if(parms->bands.length() == 0) + { + mlog(ERROR, "No bands specified"); + return false; + } for (int i = 0; i < parms->bands.length(); i++) { const std::string& bandName = parms->bands[i]; - valid = validL8Band(bandName.c_str()) || validS2Band(bandName.c_str()) || validAlgoName(bandName.c_str()); - - /* User specified invalid band name */ - if(!valid) mlog(ERROR, "Invalid band name: %s", bandName.c_str()); + if(!validL8Band(bandName.c_str()) && !validS2Band(bandName.c_str()) && !validAlgoName(bandName.c_str())) + { + mlog(ERROR, "Invalid band name: %s", bandName.c_str()); + return false; + } } - return valid; + return true; } /*---------------------------------------------------------------------------- diff --git a/packages/geo/GeoIndexedRaster.cpp b/packages/geo/GeoIndexedRaster.cpp index 39fc1dfd..b04588de 100644 --- a/packages/geo/GeoIndexedRaster.cpp +++ b/packages/geo/GeoIndexedRaster.cpp @@ -957,68 +957,78 @@ void* GeoIndexedRaster::readerThread(void *param) { GdalRaster* raster = entry->raster; - /* Open raster so we can get inner bands from it */ - raster->open(); + try + { + CHECKPTR(raster); - std::vector bands; - reader->obj->getInnerBands(raster, bands); + /* Open raster so we can get inner bands from it */ + raster->open(); - if(GdalRaster::ispoint(reader->geo)) - { - /* Sample raster bands */ - const bool oneBand = bands.size() == 1; - if(oneBand) - { - OGRPoint* p = (dynamic_cast(reader->geo)); - RasterSample* sample = raster->samplePOI(p, bands[0]); - entry->bandSample.push_back(sample); - } - else + std::vector bands; + reader->obj->getInnerBands(raster, bands); + + if(GdalRaster::ispoint(reader->geo)) { - /* Multiple bands */ - for(const int bandNum : bands) + /* Sample raster bands */ + const bool oneBand = bands.size() == 1; + if(oneBand) { - /* Use local copy of point, it will be projected in samplePOI. We do not want to project it again */ OGRPoint* p = (dynamic_cast(reader->geo)); - OGRPoint point(*p); - - RasterSample* sample = raster->samplePOI(&point, bandNum); + RasterSample* sample = raster->samplePOI(p, bands[0]); entry->bandSample.push_back(sample); - mlog(DEBUG, "Band: %d, %s", bandNum, sample ? sample->toString().c_str() : "NULL"); + } + else + { + /* Multiple bands */ + for(const int bandNum : bands) + { + /* Use local copy of point, it will be projected in samplePOI. We do not want to project it again */ + OGRPoint* p = (dynamic_cast(reader->geo)); + OGRPoint point(*p); + + RasterSample* sample = raster->samplePOI(&point, bandNum); + entry->bandSample.push_back(sample); + mlog(DEBUG, "Band: %d, %s", bandNum, sample ? sample->toString().c_str() : "NULL"); + } } } - } - else if(GdalRaster::ispoly(reader->geo)) - { - /* Subset raster bands */ - for(const int bandNum : bands) + else if(GdalRaster::ispoly(reader->geo)) { - /* No need to use local copy of polygon, subsetAOI will use it's envelope and not project it */ - RasterSubset* subset = raster->subsetAOI(dynamic_cast(reader->geo), bandNum); - if(subset) + /* Subset raster bands */ + for(const int bandNum : bands) { - /* - * Create new GeoRaster object for subsetted raster - * Use NULL for LuaState, using parent's causes memory corruption - * NOTE: cannot use RasterObject::cppCreate(parms) here, it would create - * new GeoIndexRaster with the same file path as parent raster. - */ - subset->robj = new GeoRaster(NULL, - reader->obj->rqstParms, - reader->obj->samplerKey, - subset->rasterName, - raster->getGpsTime(), - raster->getElevationBandNum(), - raster->getFLagsBandNum(), - raster->getOverrideCRS()); - - entry->bandSubset.push_back(subset); - - /* RequestFields are shared with subsseted raster and other readers */ - GeoIndexedRaster::referenceLuaObject(reader->obj->rqstParms); + /* No need to use local copy of polygon, subsetAOI will use it's envelope and not project it */ + RasterSubset* subset = raster->subsetAOI(dynamic_cast(reader->geo), bandNum); + if(subset) + { + /* + * Create new GeoRaster object for subsetted raster + * Use NULL for LuaState, using parent's causes memory corruption + * NOTE: cannot use RasterObject::cppCreate(parms) here, it would create + * new GeoIndexRaster with the same file path as parent raster. + */ + subset->robj = new GeoRaster(NULL, + reader->obj->rqstParms, + reader->obj->samplerKey, + subset->rasterName, + raster->getGpsTime(), + raster->getElevationBandNum(), + raster->getFLagsBandNum(), + raster->getOverrideCRS()); + + entry->bandSubset.push_back(subset); + + /* RequestFields are shared with subsseted raster and other readers */ + GeoIndexedRaster::referenceLuaObject(reader->obj->rqstParms); + } } } } + catch(const RunTimeException& e) + { + mlog(e.level(), "%s", e.what()); + } + entry->enabled = false; /* raster samples/subsetted */ reader->sync.lock(); @@ -1053,46 +1063,57 @@ void* GeoIndexedRaster::batchReaderThread(void *param) if(breader->uraster != NULL) { unique_raster_t* ur = breader->uraster; - GdalRaster* raster = new GdalRaster(breader->obj->parms, - breader->obj->fileDict.get(ur->rinfo->fileId), - 0, /* Sample collecting code will set it to group's gpsTime */ - ur->rinfo->fileId, - ur->rinfo->elevationBandNum, - ur->rinfo->flagsBandNum, - breader->obj->crscb); - - /* Open raster so we can get inner bands from it */ - raster->open(); - - std::vector bands; - breader->obj->getInnerBands(raster, bands); - - /* Sample all points for this raster */ - for(point_sample_t& ps : ur->pointSamples) + GdalRaster* raster = NULL; + + try { - /* Sample raster bands */ - const bool oneBand = bands.size() == 1; - if(oneBand) - { - RasterSample* sample = raster->samplePOI(&ps.point, bands[0]); - ps.bandSample.push_back(sample); - ps.bandSampleReturned.emplace_back(std::make_unique>(false)); - } - else + raster = new GdalRaster(breader->obj->parms, + breader->obj->fileDict.get(ur->rinfo->fileId), + 0, /* Sample collecting code will set it to group's gpsTime */ + ur->rinfo->fileId, + ur->rinfo->elevationBandNum, + ur->rinfo->flagsBandNum, + breader->obj->crscb); + + CHECKPTR(raster); + + /* Open raster so we can get inner bands from it */ + raster->open(); + + std::vector bands; + breader->obj->getInnerBands(raster, bands); + + /* Sample all points for this raster */ + for(point_sample_t& ps : ur->pointSamples) { - /* Multiple bands */ - for(const int bandNum : bands) + /* Sample raster bands */ + const bool oneBand = bands.size() == 1; + if(oneBand) { - /* Use local copy of point, it will be projected in samplePOI. We do not want to project it again */ - OGRPoint point(ps.point); - - RasterSample* sample = raster->samplePOI(&point, bandNum); + RasterSample* sample = raster->samplePOI(&ps.point, bands[0]); ps.bandSample.push_back(sample); ps.bandSampleReturned.emplace_back(std::make_unique>(false)); - ps.ssErrors |= raster->getSSerror(); + } + else + { + /* Multiple bands */ + for(const int bandNum : bands) + { + /* Use local copy of point, it will be projected in samplePOI. We do not want to project it again */ + OGRPoint point(ps.point); + + RasterSample* sample = raster->samplePOI(&point, bandNum); + ps.bandSample.push_back(sample); + ps.bandSampleReturned.emplace_back(std::make_unique>(false)); + ps.ssErrors |= raster->getSSerror(); + } } } } + catch(const RunTimeException& e) + { + mlog(e.level(), "%s", e.what()); + } delete raster; breader->sync.lock(); From 3a46c5ec0ad0cfcc6cf47c9bc3fe09e5121a90bc Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Thu, 14 Nov 2024 21:08:47 +0000 Subject: [PATCH 09/17] fixed bug in GeoRaster --- packages/geo/GeoRaster.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/packages/geo/GeoRaster.cpp b/packages/geo/GeoRaster.cpp index ffd5f028..c73e3487 100644 --- a/packages/geo/GeoRaster.cpp +++ b/packages/geo/GeoRaster.cpp @@ -87,12 +87,11 @@ uint32_t GeoRaster::getSamples(const MathLib::point_3d_t& point, int64_t gps, Li lockSampling(); try { - OGRPoint ogr_point(point.x, point.y, point.z); - std::vector bands; getInnerBands(&raster, bands); for(const int bandNum : bands) { + OGRPoint ogr_point(point.x, point.y, point.z); RasterSample* sample = raster.samplePOI(&ogr_point, bandNum); if(sample) slist.add(sample); } From 03941f8513e806d3a1d3b0f20ac276f6fbf9b487 Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Thu, 14 Nov 2024 21:14:53 +0000 Subject: [PATCH 10/17] added comment to bug fix --- packages/geo/GeoRaster.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/packages/geo/GeoRaster.cpp b/packages/geo/GeoRaster.cpp index c73e3487..2ff9193c 100644 --- a/packages/geo/GeoRaster.cpp +++ b/packages/geo/GeoRaster.cpp @@ -91,7 +91,9 @@ uint32_t GeoRaster::getSamples(const MathLib::point_3d_t& point, int64_t gps, Li getInnerBands(&raster, bands); for(const int bandNum : bands) { + /* Must create OGRPoint for each bandNum, samplePOI projects it to raster CRS */ OGRPoint ogr_point(point.x, point.y, point.z); + RasterSample* sample = raster.samplePOI(&ogr_point, bandNum); if(sample) slist.add(sample); } From d06829dca35d35d61a824734ff8c0ee2ceccf09b Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Fri, 15 Nov 2024 10:07:15 +0000 Subject: [PATCH 11/17] skip resampling and zonal stats for fmask band --- packages/geo/GdalRaster.cpp | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/packages/geo/GdalRaster.cpp b/packages/geo/GdalRaster.cpp index 47dc72a9..bac1c50a 100644 --- a/packages/geo/GdalRaster.cpp +++ b/packages/geo/GdalRaster.cpp @@ -234,13 +234,22 @@ RasterSample* GdalRaster::samplePOI(OGRPoint* poi, int bandNum) { const double vertical_shift = z - poi->getZ(); sample = new RasterSample(gpsTime, fileId, vertical_shift); - if(parms->sampling_algo == GRIORA_NearestNeighbour) + + if(band == flagsBand) + { + /* Skip resampling and zonal stats for quality mask band (value is bitmask) */ readPixel(poi, band, sample); + } else - resamplePixel(poi, band, sample); + { + if(parms->sampling_algo == GRIORA_NearestNeighbour) + readPixel(poi, band, sample); + else + resamplePixel(poi, band, sample); - if(parms->zonal_stats) - computeZonalStats(poi, band, sample); + if(parms->zonal_stats) + computeZonalStats(poi, band, sample); + } } else { From 0f380ac136faecfdf399dc485ddaad4c06df1a9b Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Fri, 15 Nov 2024 18:50:12 +0000 Subject: [PATCH 12/17] Updated code to support GEBCO 2024 dataset --- datasets/gebco/package/GebcoBathyRaster.cpp | 28 ++++++++- datasets/gebco/selftests/gebco_reader.lua | 61 ++++++++++++++++++- datasets/gebco/utils/README.txt | 47 ++++++++++++++ datasets/gebco/utils/cog_maker.py | 59 +++++++++++++++--- datasets/gebco/utils/index_maker.py | 49 ++++++++++----- .../slideruleearth-aws/asset_directory.csv | 2 +- 6 files changed, 217 insertions(+), 29 deletions(-) create mode 100644 datasets/gebco/utils/README.txt diff --git a/datasets/gebco/package/GebcoBathyRaster.cpp b/datasets/gebco/package/GebcoBathyRaster.cpp index bfaf091d..c752f0a5 100644 --- a/datasets/gebco/package/GebcoBathyRaster.cpp +++ b/datasets/gebco/package/GebcoBathyRaster.cpp @@ -47,6 +47,32 @@ GebcoBathyRaster::GebcoBathyRaster(lua_State* L, RequestFields* rqst_parms, cons filePath("/vsis3/" + std::string(parms->asset.asset->getPath())), indexFile(parms->asset.asset->getIndex()) { + /* + * To support datasets from different years, we use the 'bands' parameter to identify which year's data to sample. + * This band parameter must be cleared in 'parms' since it does not correspond to an actual band. + */ + std::string year("2024"); + if(parms->bands.length() == 0) + { + mlog(INFO, "Using latest GEBCO data from %s", year.c_str()); + } + else if(parms->bands.length() == 1) + { + if(parms->bands[0] == "2023" || parms->bands[0] == "2024") + { + year = parms->bands[0]; + mlog(INFO, "Using GEBCO data from %s", year.c_str()); + + /* Clear params->bands */ + const_cast&>(parms->bands).clear(); + } + else throw RunTimeException(CRITICAL, RTE_ERROR, "Invalid band name specified"); + } + else throw RunTimeException(CRITICAL, RTE_ERROR, "Invalid number of bands specified"); + + /* Build data path on s3 */ + filePath += "/" + year; + mlog(DEBUG, "Using data path: %s", filePath.c_str()); } /*---------------------------------------------------------------------------- @@ -62,7 +88,7 @@ void GebcoBathyRaster::getIndexFile(const OGRGeometry* geo, std::string& file, c static_cast(geo); static_cast(points); file = filePath + "/" + indexFile; - mlog(DEBUG, "Using %s", file.c_str()); + mlog(DEBUG, "Using index file: %s", file.c_str()); } diff --git a/datasets/gebco/selftests/gebco_reader.lua b/datasets/gebco/selftests/gebco_reader.lua index 8e24d124..14846ff6 100644 --- a/datasets/gebco/selftests/gebco_reader.lua +++ b/datasets/gebco/selftests/gebco_reader.lua @@ -15,7 +15,8 @@ local lons = {-40, -20, -120} local lats = { 70, -20, 20} height = 0 -local expDepth = { -64, -4933, -4072} +-- Updataed results for GEBCO 2024 +local expDepth = { -64, -4940, -4072} local expFlags = { 70, 40, 44} -- Tolerance for depth values is needed because of the different versions of PROJLIB @@ -23,7 +24,7 @@ local expFlags = { 70, 40, 44} -- We add 1 meter of depth tolerance for the test to pass local depth_tolerance = 1; -print(string.format("\n--------------------------------\nTest: GEBCO Correct Values\n--------------------------------")) +print(string.format("\n--------------------------------\nTest: GEBCO Correct Values default with no band specified\n--------------------------------")) local dem = geo.raster(geo.parms({ asset = "gebco-bathy", algorithm = "NearestNeighbour", with_flags=true, sort_by_index = true })) runner.check(dem ~= nil) @@ -49,6 +50,62 @@ for j, lon in ipairs(lons) do end end +print(string.format("\n--------------------------------\nTest: GEBCO Correct Values band=2024\n--------------------------------")) +dem = geo.raster(geo.parms({ asset = "gebco-bathy", bands = {"2024"}, algorithm = "NearestNeighbour", with_flags=true, sort_by_index = true })) +runner.check(dem ~= nil) + +for j, lon in ipairs(lons) do + lat = lats[j] + tbl, err = dem:sample(lon, lat, height) + runner.check(err == 0) + runner.check(tbl ~= nil) + + if err ~= 0 then + print(string.format("Point: %d, (%.3f, %.3f) ======> FAILED to read, err# %d",j, lon, lat, err)) + else + local el, fname + for k, v in ipairs(tbl) do + value = v["value"] + fname = v["file"] + flags = v["flags"] + print(string.format("(%02d) (%6.1f, %5.1f) %8.1fm %02d %s", k, lon, lat, value, flags, fname)) + + assert( math.abs(value - expDepth[j]) < depth_tolerance, string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) + assert(flags == expFlags[j], string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) + end + end +end + + +-- Different results for GEBCO 2023 +expDepth = { -64, -4933, -4072} + +print(string.format("\n--------------------------------\nTest: GEBCO Correct Values band=2023\n--------------------------------")) +dem = geo.raster(geo.parms({ asset = "gebco-bathy", bands = {"2023"}, algorithm = "NearestNeighbour", with_flags=true, sort_by_index = true })) +runner.check(dem ~= nil) + +for j, lon in ipairs(lons) do + lat = lats[j] + tbl, err = dem:sample(lon, lat, height) + runner.check(err == 0) + runner.check(tbl ~= nil) + + if err ~= 0 then + print(string.format("Point: %d, (%.3f, %.3f) ======> FAILED to read, err# %d",j, lon, lat, err)) + else + local el, fname + for k, v in ipairs(tbl) do + value = v["value"] + fname = v["file"] + flags = v["flags"] + print(string.format("(%02d) (%6.1f, %5.1f) %8.1fm %02d %s", k, lon, lat, value, flags, fname)) + + assert( math.abs(value - expDepth[j]) < depth_tolerance, string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) + assert(flags == expFlags[j], string.format("Point: %d, (%.3f, %.3f) ======> FAILED",j, lon, lat)) + end + end +end + -- Report Results -- runner.report() diff --git a/datasets/gebco/utils/README.txt b/datasets/gebco/utils/README.txt new file mode 100644 index 00000000..7f4ed0c3 --- /dev/null +++ b/datasets/gebco/utils/README.txt @@ -0,0 +1,47 @@ +GEBCO Data Processing Guide + +The GEBCO data is not hosted on AWS. Users must download the compressed files and run the provided scripts to create COGs and a GeoJSON index file. +Below are the steps to perform this process. In these examples, the year "2024" is used, but the Python scripts support the "year" parameter for other years as well. + +*** Create Directory Structure *** +Create the following directories with the specified year in the path: + +/data/GEBCO/2024/TIFS +/data/GEBCO/2024/COGS + +*** Download the GEBCO Data *** +Go to: https://www.gebco.net/data_and_products/gridded_bathymetry_data/#a1 +Find the table entry "GEBCO_2024 Grid (sub-ice topo/bathy)" download DataGeoTiff +It should be named something like: gebco_2024_sub_ice_topo_geotiff.zip + +Find the table entry "GEBCO_2024 TID Grid" download DataGeoTiff +It should be named something like: gebco_2024_tid_geotiff.zip + +Uncompress both files into: /data/GEBCO/2024/TIFS +>> unzip gebco_2024_sub_ice_topo_geotiff.zip -d /data/GEBCO/2024/TIFS +>> unzip gebco_2024_tid_geotiff.zip -d /data/GEBCO/2024/TIFS + + +*** Run the Python Scripts *** +From the utils directory where the scripts are located, run the following scripts in the specified order. +Make sure to include the year as a parameter. + +Step 1: Run the script cog_maker.py with the year as a parameter +>> python cog_maker.py 2024 + +Step 2: Run the script index_maker.py with the year as a parameter +>> python index_maker.py 2024 + + +*** Upload to Sliderule S3 Bucket *** +In the Sliderule public bucket (s3://sliderule/data/GEBCO/), create a new directory for the year, e.g., 2024. + +Upload all COG tif files and index.geojso from /data/GEBCO/2024/COGS to: s3://sliderule/data/GEBCO/2024 +>> aws s3 cp . s3://sliderule/data/GEBCO/2024/ --recursive + + + +Upload the original downloaded data files to: s3://sliderule/data/GEBCO/2024 +>> aws s3 cp gebco_2024_sub_ice_topo_geotiff.zip s3://sliderule/data/GEBCO/2024/gebco_2024_sub_ice_topo_geotiff.zip +>> aws s3 cp gebco_2024_tid_geotiff.zip s3://sliderule/data/GEBCO/2024/gebco_2024_tid_geotiff.zip + diff --git a/datasets/gebco/utils/cog_maker.py b/datasets/gebco/utils/cog_maker.py index 85fe0331..b6e369d1 100644 --- a/datasets/gebco/utils/cog_maker.py +++ b/datasets/gebco/utils/cog_maker.py @@ -1,25 +1,37 @@ import os -from osgeo import gdal, org +import sys +from osgeo import gdal + +# Enable GDAL exceptions to avoid future warnings +gdal.UseExceptions() # Set GDAL configurations for large cache and multithreading gdal.SetConfigOption("GDAL_CACHEMAX", "1024") # 1024 MB cache gdal.SetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS") # Use all available CPU cores -def convert_geotiffs_to_cogs(tiffs_dir, cogs_directory): +def convert_geotiffs_to_cogs(year, tiffs_dir, cogs_dir): # Empty the COGS directory if it exists - if os.path.exists(cogs_directory): - for file in os.listdir(cogs_directory): - os.remove(os.path.join(cogs_directory, file)) + if os.path.exists(cogs_dir): + for file in os.listdir(cogs_dir): + os.remove(os.path.join(cogs_dir, file)) else: - os.makedirs(cogs_directory) + os.makedirs(cogs_dir) # List of GeoTIFF files to convert to COGs geotiff_files = [f for f in os.listdir(tiffs_dir) if f.endswith(".tif")] + # tif files have year in their name, example of 2023 tif file: + # gebco_2023_sub_ice_n0.0_s-90.0_w-180.0_e-90.0.tif + # Check if the year is part of the file name + for geotiff_file in geotiff_files: + if year not in geotiff_file: + print(f"Error: File '{geotiff_file}' does not contain the year '{year}' in its name.") + sys.exit(1) + # Convert each GeoTIFF to a Cloud-Optimized GeoTIFF for geotiff_file in geotiff_files: input_path = os.path.join(tiffs_dir, geotiff_file) - output_cog_path = os.path.join(cogs_directory, f"cog_{geotiff_file}") + output_cog_path = os.path.join(cogs_dir, f"cog_{geotiff_file}") gdal.Translate( output_cog_path, input_path, @@ -38,7 +50,34 @@ def convert_geotiffs_to_cogs(tiffs_dir, cogs_directory): # Convert original tif files to COGs # https://www.gebco.net/data_and_products/gridded_bathymetry_data/#area -tifs_directory = "/data/GEBCO/2023/TIFS" -cogs_directory = "/data/GEBCO/2023/COGS" -convert_geotiffs_to_cogs(tifs_directory, cogs_directory) +if __name__ == "__main__": + # List of supported years + supported_years = ["2023", "2024"] + + if len(sys.argv) != 2: + print("Usage: python cog_maker.py ") + sys.exit(1) + + # Get the year from the command line argument + year = sys.argv[1] + + # Ensure the provided year is valid + if year not in supported_years: + print(f"Error: Invalid year. Supported years are: {', '.join(supported_years)}") + sys.exit(1) + + # Compose directory paths based on the year + tifs_dir = f"/data/GEBCO/{year}/TIFS" + cogs_dir = f"/data/GEBCO/{year}/COGS" + + # Check if the directories exist + if not os.path.exists(tifs_dir): + print(f"Error: Directory '{tifs_dir}' does not exist.") + sys.exit(1) + + if not os.path.exists(cogs_dir): + print(f"Error: Directory '{cogs_dir}' does not exist.") + sys.exit(1) + # Call the conversion function + convert_geotiffs_to_cogs(year, tifs_dir, cogs_dir) \ No newline at end of file diff --git a/datasets/gebco/utils/index_maker.py b/datasets/gebco/utils/index_maker.py index 4175b973..78431b50 100644 --- a/datasets/gebco/utils/index_maker.py +++ b/datasets/gebco/utils/index_maker.py @@ -1,19 +1,22 @@ import os +import sys from osgeo import gdal, ogr +# Enable GDAL exceptions to avoid future warnings +gdal.UseExceptions() -def create_geojson(directory_path, geojson_name): - if not os.path.exists(directory_path): - raise ValueError(f"The directory '{directory_path}' does not exist.") +def create_geojson(dir_path, index_file): + if not os.path.exists(dir_path): + raise ValueError(f"The directory '{dir_path}' does not exist.") # Find all .tif files in the directory - tif_files = [os.path.join(directory_path, f) for f in os.listdir(directory_path) if f.endswith(".tif")] + tif_files = [os.path.join(dir_path, f) for f in os.listdir(dir_path) if f.endswith(".tif")] if not tif_files: - raise ValueError("No .vrt files found in the specified directory.") + raise ValueError(f"No .tif files found in directory: {dir_path}") # Set the output path for the geojson file in the same directory - geojson_output_path = os.path.join(directory_path, geojson_name) + geojson_output_path = os.path.join(dir_path, index_file) # Remove if it already exists if os.path.exists(geojson_output_path): @@ -48,7 +51,7 @@ def create_geojson(directory_path, geojson_name): if '_tid_' in tif_file: continue - print(f"TIF File: {tif_file}") + # print(f"Processed: {tif_file}") # Get the geotransform and size of the raster geotransform = raster_ds.GetGeoTransform() @@ -80,9 +83,9 @@ def create_geojson(directory_path, geojson_name): feature.SetField("id", id) id += 1 - # Set the "datetime" attribute to the middle of 2023 - middle_of_2023 = "2023-07-02T00:00:00Z" - feature.SetField("datetime", middle_of_2023) + # Set the "datetime" attribute to the middle of the specified year + middle_of_year = f"{year}-07-02T00:00:00Z" + feature.SetField("datetime", middle_of_year) # Remove the directory path from the tif file name tif_file = os.path.basename(tif_file) @@ -101,11 +104,27 @@ def create_geojson(directory_path, geojson_name): geojson_ds.Destroy() # Close and save the GeoJSON file - print(f"GeoJSON created at: {geojson_output_path}") + print(f"GeoJSON index file created at: {geojson_output_path}") -# Make index geojson for tif files -cogs_directory = "/data/GEBCO/2023/COGS" -geojson_name = "index.geojson" -create_geojson(cogs_directory, geojson_name) +if __name__ == "__main__": + # List of supported years + supported_years = ["2023", "2024"] + if len(sys.argv) != 2: + print("Usage: python index_maker.py ") + sys.exit(1) + + # Get the year from the command line argument + year = sys.argv[1] + + # Ensure the provided year is valid + if year not in supported_years: + print(f"Error: Invalid year. Supported years are: {', '.join(supported_years)}") + sys.exit(1) + + # Compose directory paths based on the year + cogs_dir = f"/data/GEBCO/{year}/COGS" + index_file = "index.geojson" + + create_geojson(cogs_dir, index_file) diff --git a/targets/slideruleearth-aws/asset_directory.csv b/targets/slideruleearth-aws/asset_directory.csv index 5c99bc51..998c1ba8 100644 --- a/targets/slideruleearth-aws/asset_directory.csv +++ b/targets/slideruleearth-aws/asset_directory.csv @@ -19,7 +19,7 @@ swot-sim-glorys, podaac-cloud, s3, podaac-ops-cumulus-prote usgs3dep-1meter-dem, nil, nil, /vsis3/prd-tnm, nil, us-west-2, https://s3.us-west-2.amazonaws.com esa-worldcover-10meter, nil, nil, /vsis3/esa-worldcover/v200/2021/map, /vsis3/sliderule/data/WORLDCOVER/ESA_WorldCover_10m_2021_v200_Map.vrt, eu-central-1, https://s3.eu-central-1.amazonaws.com meta-globalcanopy-1meter, nil, nil, /vsis3/dataforgood-fb-data/forests/v1/alsgedi_global_v6_float/, /vsis3/sliderule/data/GLOBALCANOPY/META_GlobalCanopyHeight_1m_2024_v1.vrt, us-east-1, https://s3.us-east-1.amazonaws.com -gebco-bathy, iam-role, s3, sliderule/data/GEBCO/2023, index.geojson, us-west-2, https://s3.us-west-2.amazonaws.com +gebco-bathy, iam-role, s3, sliderule/data/GEBCO, index.geojson, us-west-2, https://s3.us-west-2.amazonaws.com bluetopo-bathy, iam-role, s3, /vsis3/noaa-ocs-nationalbathymetry-pds/BlueTopo/, _BlueTopo_Tile_Scheme, us-east-1, https://s3.us-east-1.amazonaws.com landsat-hls, lpdaac-cloud, nil, /vsis3/lp-prod-protected, nil, us-west-2, https://s3.us-west-2.amazonaws.com arcticdem-mosaic, nil, nil, /vsis3/pgc-opendata-dems/arcticdem/mosaics/v4.1, 2m_dem_tiles.vrt, us-west-2, https://s3.us-west-2.amazonaws.com From 25125fb7f50febdf346457229ff410788cba6f12 Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Fri, 15 Nov 2024 19:50:05 +0000 Subject: [PATCH 13/17] updated gebco python test for 2023 and 2024 data --- clients/python/tests/test_gebco.py | 43 ++++++++++++++++++++++++------ 1 file changed, 35 insertions(+), 8 deletions(-) diff --git a/clients/python/tests/test_gebco.py b/clients/python/tests/test_gebco.py index 5398f66f..314bb8e4 100644 --- a/clients/python/tests/test_gebco.py +++ b/clients/python/tests/test_gebco.py @@ -8,18 +8,45 @@ sigma = 1.0e-9 -lon = -40 -lat = 70 +lon = -20 +lat = -20 -value = -64 # meters -file = '/vsis3/sliderule/data/GEBCO/2023/cog_gebco_2023_sub_ice_n90.0_s0.0_w-90.0_e0.0.tif' +expected_tid = 40 # fmask for type identifier from TID grid + +expected_depth_2023 = -4933 # meters +expected_file_2023 = '/vsis3/sliderule/data/GEBCO/2023/cog_gebco_2023_sub_ice_n0.0_s-90.0_w-90.0_e0.0.tif' + +expected_depth_2024 = -4940 # different from 2023 +expected_file_2024 = '/vsis3/sliderule/data/GEBCO/2024/cog_gebco_2024_sub_ice_n0.0_s-90.0_w-90.0_e0.0.tif' @pytest.mark.network class TestGebco: - def test_samples(self, init): - rqst = {"samples": {"asset": "gebco-bathy"}, "coordinates": [[lon,lat]]} + def test_samples_2023(self, init): + rqst = {"samples": {"asset": "gebco-bathy", "with_flags": True, "bands": ["2023"]}, "coordinates": [[lon,lat]]} + rsps = sliderule.source("samples", rqst) + assert init + print(rsps) + assert rsps["samples"][0][0]["value"] == expected_depth_2023 + assert rsps["samples"][0][0]["file"] == expected_file_2023 + assert rsps["samples"][0][0]["flags"] == expected_tid + + def test_samples_2024(self, init): + value = -64 # meters + rqst = {"samples": {"asset": "gebco-bathy", "with_flags": True, "bands": ["2024"]}, "coordinates": [[lon,lat]]} + rsps = sliderule.source("samples", rqst) + assert init + print(rsps) + assert rsps["samples"][0][0]["value"] == expected_depth_2024 + assert rsps["samples"][0][0]["file"] == expected_file_2024 + assert rsps["samples"][0][0]["flags"] == expected_tid + + # Test default data set (no bands parameter) + def test_samples_default(self, init): + value = -64 # meters + rqst = {"samples": {"asset": "gebco-bathy", "with_flags": True}, "coordinates": [[lon,lat]]} rsps = sliderule.source("samples", rqst) assert init print(rsps) - assert rsps["samples"][0][0]["value"] == value - assert rsps["samples"][0][0]["file"] == file + assert rsps["samples"][0][0]["value"] == expected_depth_2024 + assert rsps["samples"][0][0]["file"] == expected_file_2024 + assert rsps["samples"][0][0]["flags"] == expected_tid From ad9ed05a743285e53468eb5bd75d02ef59064eb7 Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Mon, 18 Nov 2024 17:58:40 +0000 Subject: [PATCH 14/17] fixed int64 time overlfow, fixes #445 --- clients/python/sliderule/raster.py | 5 ++++- clients/python/tests/test_worldcover.py | 25 +++++++++++++++++++++++++ 2 files changed, 29 insertions(+), 1 deletion(-) diff --git a/clients/python/sliderule/raster.py b/clients/python/sliderule/raster.py index ca70ae9f..899a5542 100644 --- a/clients/python/sliderule/raster.py +++ b/clients/python/sliderule/raster.py @@ -134,6 +134,9 @@ def sample(asset, coordinates, parms={}): **columns } + # GPS Offset between unix epoch and gps epoch in milliseconds + gps_offset_msecs = 315964800.0 * 1000 + # Populate Columns i = 0 per_coord_index = 0 @@ -142,7 +145,7 @@ def sample(asset, coordinates, parms={}): latitude = coordinates[per_coord_index][1] per_coord_index += 1 for raster_sample in input_coord_response: - columns['time'][i] = numpy.int64((raster_sample['time'] + 315964800.0) * 1000000000) + columns['time'][i] = numpy.int64((raster_sample['time'] + gps_offset_msecs) * 1e6) # int64 with gps nanoseconds will overflow in year 2272 columns['longitude'][i] = numpy.double(longitude) columns['latitude'][i] = numpy.double(latitude) columns['file'] += raster_sample['file'], diff --git a/clients/python/tests/test_worldcover.py b/clients/python/tests/test_worldcover.py index 6a36e0df..01fb5c62 100644 --- a/clients/python/tests/test_worldcover.py +++ b/clients/python/tests/test_worldcover.py @@ -3,6 +3,8 @@ import pytest from pathlib import Path import sliderule +from sliderule import raster +import geopandas as gpd TESTDIR = Path(__file__).parent @@ -24,3 +26,26 @@ def test_vrt(self, init): assert abs(rsps["samples"][0][0]["value"] - vrtValue) < sigma assert rsps["samples"][0][0]["file"] == vrtFile assert rsps["samples"][0][0]["time"] == vrtFileTime + + def test_time_overflow(self, init): + region = [ {"lon": -108.34, "lat": 38.89}, + {"lon": -107.76, "lat": 38.90}, + {"lon": -107.78, "lat": 39.26}, + {"lon": -108.36, "lat": 39.25}, + {"lon": -108.34, "lat": 38.89} ] + + gfp = gpd.GeoDataFrame(geometry=gpd.points_from_xy([-108.1, -108.3], [39.1, 39.2]), crs='EPSG:4326') + points = [[x,y] for x,y in zip(gfp.geometry.x , gfp.geometry.y)] + gdf = raster.sample("esa-worldcover-10meter", points) + + expected_time = "2021-06-30 00:00:18" + expected_value = 10.0 + + # Check if all timestamps match the expected_time (used as index) + assert all(str(index) == expected_time for index in gdf.index), \ + f"Time mismatch. Expected: {expected_time}, Found: {gdf.index.tolist()}" + + # Check if all values match the expected value + assert all(gdf["value"] == expected_value), \ + f"Value mismatch. Expected: {expected_value}, Found: {gdf['value'].tolist()}" + From 9500edf4fae6931348c8037c050fc76135de6dcc Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Tue, 19 Nov 2024 21:11:05 +0000 Subject: [PATCH 15/17] cleanup and update of gpsTime in RasterSample --- clients/python/sliderule/raster.py | 5 +-- clients/python/tests/test_bluetopo.py | 2 +- clients/python/tests/test_globalcanopy.py | 2 +- clients/python/tests/test_worldcover.py | 29 +++++++------ .../bluetopo/package/BlueTopoBathyRaster.cpp | 2 +- .../bluetopo/selftests/bluetopo_reader.lua | 6 +-- datasets/gebco/package/GebcoBathyRaster.cpp | 2 +- datasets/landsat/package/LandsatHlsRaster.cpp | 22 ++++++---- .../package/EsaWorldCover10meterRaster.h | 2 +- .../package/MetaGlobalCanopy1meterRaster.h | 2 +- .../selftests/globalcanopy_reader.lua | 2 +- .../opendata/selftests/worldcover_reader.lua | 2 +- datasets/pgc/package/PgcDemStripsRaster.cpp | 12 +++--- .../package/Usgs3dep1meterDemRaster.cpp | 2 +- packages/geo/GeoIndexedRaster.cpp | 41 +++++++++++-------- packages/geo/GeoIndexedRaster.h | 6 +-- packages/geo/RasterObject.cpp | 2 +- 17 files changed, 78 insertions(+), 63 deletions(-) diff --git a/clients/python/sliderule/raster.py b/clients/python/sliderule/raster.py index 899a5542..e8779b22 100644 --- a/clients/python/sliderule/raster.py +++ b/clients/python/sliderule/raster.py @@ -134,9 +134,6 @@ def sample(asset, coordinates, parms={}): **columns } - # GPS Offset between unix epoch and gps epoch in milliseconds - gps_offset_msecs = 315964800.0 * 1000 - # Populate Columns i = 0 per_coord_index = 0 @@ -145,7 +142,7 @@ def sample(asset, coordinates, parms={}): latitude = coordinates[per_coord_index][1] per_coord_index += 1 for raster_sample in input_coord_response: - columns['time'][i] = numpy.int64((raster_sample['time'] + gps_offset_msecs) * 1e6) # int64 with gps nanoseconds will overflow in year 2272 + columns['time'][i] = numpy.int64((raster_sample['time'] + 315564800.0) * 1e9) columns['longitude'][i] = numpy.double(longitude) columns['latitude'][i] = numpy.double(latitude) columns['file'] += raster_sample['file'], diff --git a/clients/python/tests/test_bluetopo.py b/clients/python/tests/test_bluetopo.py index 9b1595a6..119564e1 100644 --- a/clients/python/tests/test_bluetopo.py +++ b/clients/python/tests/test_bluetopo.py @@ -11,7 +11,7 @@ lon = -80.87 lat = 32.06 -file = '/vsis3/noaa-ocs-nationalbathymetry-pds/BlueTopo/BH4SX59B/BlueTopo_BH4SX59B_20241021.tiff' +file = '/vsis3/noaa-ocs-nationalbathymetry-pds/BlueTopo/BH4SX59B/BlueTopo_BH4SX59B_20241115.tiff' expElevation = -2.02 expUncertainty = 1.00 diff --git a/clients/python/tests/test_globalcanopy.py b/clients/python/tests/test_globalcanopy.py index 4378fec3..8854e863 100644 --- a/clients/python/tests/test_globalcanopy.py +++ b/clients/python/tests/test_globalcanopy.py @@ -13,7 +13,7 @@ vrtValue = 20 vrtFile = '/vsis3/sliderule/data/GLOBALCANOPY/META_GlobalCanopyHeight_1m_2024_v1.vrt' -vrtFileTime = 1396483218000 +vrtFileTime = 1396483218 @pytest.mark.network class TestGlobalCanopy: diff --git a/clients/python/tests/test_worldcover.py b/clients/python/tests/test_worldcover.py index 01fb5c62..d0fb247a 100644 --- a/clients/python/tests/test_worldcover.py +++ b/clients/python/tests/test_worldcover.py @@ -13,9 +13,9 @@ vrtLon = -108.1 vrtLat = 39.1 -vrtValue = 10 +expValue = 10 vrtFile = '/vsis3/sliderule/data/WORLDCOVER/ESA_WorldCover_10m_2021_v200_Map.vrt' -vrtFileTime = 1309046418000 +vrtFileTime = 1309046418 # gpsTime in seconds 2021-06-30 @pytest.mark.network class TestMosaic: @@ -23,9 +23,9 @@ def test_vrt(self, init): rqst = {"samples": {"asset": "esa-worldcover-10meter"}, "coordinates": [[vrtLon,vrtLat]]} rsps = sliderule.source("samples", rqst) assert init - assert abs(rsps["samples"][0][0]["value"] - vrtValue) < sigma + assert abs(rsps["samples"][0][0]["value"] - expValue) < sigma assert rsps["samples"][0][0]["file"] == vrtFile - assert rsps["samples"][0][0]["time"] == vrtFileTime + assert rsps["samples"][0][0]["time"] == vrtFileTime # datetime in seconds def test_time_overflow(self, init): region = [ {"lon": -108.34, "lat": 38.89}, @@ -38,14 +38,19 @@ def test_time_overflow(self, init): points = [[x,y] for x,y in zip(gfp.geometry.x , gfp.geometry.y)] gdf = raster.sample("esa-worldcover-10meter", points) - expected_time = "2021-06-30 00:00:18" - expected_value = 10.0 + # Reset the index to make 'time' a regular column + gdf_reset = gdf.reset_index() - # Check if all timestamps match the expected_time (used as index) - assert all(str(index) == expected_time for index in gdf.index), \ - f"Time mismatch. Expected: {expected_time}, Found: {gdf.index.tolist()}" + # Access raw time column values as NumPy array + raw_time_values = gdf_reset['time'].astype('int64') - # Check if all values match the expected value - assert all(gdf["value"] == expected_value), \ - f"Value mismatch. Expected: {expected_value}, Found: {gdf['value'].tolist()}" + # Expected gps epoch time in nanoseconds + expected_time = (vrtFileTime + 315564800) * 1e9 + print("Raw time column values (NumPy array):", raw_time_values) + print("expected_time:", expected_time) + + assert (raw_time_values == expected_time).all() + + values = gdf_reset['value'] + assert (values == expValue).all() diff --git a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp index 92313399..bfee12bb 100644 --- a/datasets/bluetopo/package/BlueTopoBathyRaster.cpp +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp @@ -101,7 +101,7 @@ bool BlueTopoBathyRaster::findRasters(raster_finder_t* finder) if (!rastergeo->Intersects(geo)) continue; rasters_group_t* rgroup = new rasters_group_t; - rgroup->gpsTime = getGmtDate(feature, "Delivered_Date", rgroup->gmtDate); + rgroup->gpsTime = getGmtDate(feature, "Delivered_Date", rgroup->gmtDate) / 1000.0; const char* dataFile = feature->GetFieldAsString("GeoTIFF_link"); if(dataFile && (strlen(dataFile) > 0)) diff --git a/datasets/bluetopo/selftests/bluetopo_reader.lua b/datasets/bluetopo/selftests/bluetopo_reader.lua index 846fced5..d1663d5b 100644 --- a/datasets/bluetopo/selftests/bluetopo_reader.lua +++ b/datasets/bluetopo/selftests/bluetopo_reader.lua @@ -15,9 +15,9 @@ local lons = {-80.87, -81.02, -89.66, -94.72} local lats = { 32.06, 31.86, 29.99, 29.35} height = 0 -local expElevation = { -2.02, -14.03, -4.28, -17.18} -local expUncertainty = { 1.00, 4.37, 0.34, 1.32} -local expContributor = { 31156, 0, 24955, 45641} +local expElevation = { -2.02, -14.10, -4.28, -17.18} +local expUncertainty = { 1.00, 2.58, 0.34, 1.32} +local expContributor = { 31156, 63846, 24955, 45641} local elevation_tolerance = 0.01; diff --git a/datasets/gebco/package/GebcoBathyRaster.cpp b/datasets/gebco/package/GebcoBathyRaster.cpp index c752f0a5..a06e6c41 100644 --- a/datasets/gebco/package/GebcoBathyRaster.cpp +++ b/datasets/gebco/package/GebcoBathyRaster.cpp @@ -112,7 +112,7 @@ bool GebcoBathyRaster::findRasters(raster_finder_t* finder) if (!rastergeo->Intersects(geo)) continue; rasters_group_t* rgroup = new rasters_group_t; - rgroup->gpsTime = getGmtDate(feature, DATE_TAG, rgroup->gmtDate); + rgroup->gpsTime = getGmtDate(feature, DATE_TAG, rgroup->gmtDate) / 1000.0; const char* dataFile = feature->GetFieldAsString("data_raster"); if(dataFile && (strlen(dataFile) > 0)) diff --git a/datasets/landsat/package/LandsatHlsRaster.cpp b/datasets/landsat/package/LandsatHlsRaster.cpp index bae88c57..78de2718 100644 --- a/datasets/landsat/package/LandsatHlsRaster.cpp +++ b/datasets/landsat/package/LandsatHlsRaster.cpp @@ -174,7 +174,7 @@ bool LandsatHlsRaster::findRasters(raster_finder_t* finder) /* Set raster group time and group featureId */ rasters_group_t* rgroup = new rasters_group_t; rgroup->featureId = StringLib::duplicate(feature->GetFieldAsString("id")); - rgroup->gpsTime = getGmtDate(feature, DATE_TAG, rgroup->gmtDate); + rgroup->gpsTime = getGmtDate(feature, DATE_TAG, rgroup->gmtDate) / 1000.0; /* Find each requested band in the index file */ for(auto it = bandsDict.begin(); it != bandsDict.end(); it++) @@ -336,7 +336,7 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr { const char* key = fileDictGet(rinfo.fileId); cacheitem_t* item; - if(cache.find(key, &item)) + if(cache.find(key, &item) && !item->bandSample.empty()) { RasterSample* sample = item->bandSample[INNER_BAND_INDX]; @@ -389,7 +389,13 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr { if(ps.pointIndex == pointIndx) { - RasterSample* sample = ps.bandSample[INNER_BAND_INDX]; + RasterSample* sample = NULL; + + /* bandSample can be empty if raster failed to open */ + if(!ps.bandSample.empty()) + { + sample = ps.bandSample[INNER_BAND_INDX]; + } /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ if(sample == NULL) break; @@ -406,23 +412,23 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr RasterSample* s; if(!ps.bandSampleReturned[INNER_BAND_INDX]->exchange(true)) { - s = ps.bandSample[INNER_BAND_INDX]; + s = sample; } else { /* Sample has already been returned, must create a copy */ - s = new RasterSample(*ps.bandSample[INNER_BAND_INDX]); + s = new RasterSample(*sample); } /* Set band name for this sample */ s->bandName = bandName; /* Set time for this sample */ - s->time = rgroup->gpsTime / 1000; + s->time = rgroup->gpsTime; /* Set flags for this sample, add it to the list */ s->flags = flags; - sampleVect.push_back(sample); + sampleVect.push_back(s); errors |= ps.ssErrors; } } @@ -452,7 +458,7 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr throw RunTimeException(DEBUG, RTE_ERROR, "Invalid sample mode"); } - const double groupTime = rgroup->gpsTime / 1000; + const double groupTime = rgroup->gpsTime; const std::string groupName = featureId + " {\"algo\": \""; /* Calculate algos - make sure that all the necessary bands were read */ diff --git a/datasets/opendata/package/EsaWorldCover10meterRaster.h b/datasets/opendata/package/EsaWorldCover10meterRaster.h index 0c069a2b..02eba9de 100644 --- a/datasets/opendata/package/EsaWorldCover10meterRaster.h +++ b/datasets/opendata/package/EsaWorldCover10meterRaster.h @@ -63,7 +63,7 @@ class EsaWorldCover10meterRaster: public GeoRaster EsaWorldCover10meterRaster (lua_State* L, RequestFields* rqst_parms, const char* key): GeoRaster(L, rqst_parms, key, rqst_parms->geoFields(key)->asset.asset->getIndex(), - TimeLib::datetime2gps(2021, 06, 30, 0, 0, 0) /* Mid point for year data was collected */) {} + TimeLib::datetime2gps(2021, 06, 30, 0, 0, 0) / 1000 /* Mid point for year data was collected */) {} private: diff --git a/datasets/opendata/package/MetaGlobalCanopy1meterRaster.h b/datasets/opendata/package/MetaGlobalCanopy1meterRaster.h index 2f1b3b6c..36233562 100644 --- a/datasets/opendata/package/MetaGlobalCanopy1meterRaster.h +++ b/datasets/opendata/package/MetaGlobalCanopy1meterRaster.h @@ -63,7 +63,7 @@ class MetaGlobalCanopy1meterRaster: public GeoRaster MetaGlobalCanopy1meterRaster (lua_State* L, RequestFields* rqst_parms, const char* key): GeoRaster(L, rqst_parms, key, rqst_parms->geoFields(key)->asset.asset->getIndex(), - TimeLib::datetime2gps(2024, 4, 7, 0, 0, 0)) {} + TimeLib::datetime2gps(2024, 4, 7, 0, 0, 0) / 1000) {} }; diff --git a/datasets/opendata/selftests/globalcanopy_reader.lua b/datasets/opendata/selftests/globalcanopy_reader.lua index 5c9bd46c..ac7d941e 100644 --- a/datasets/opendata/selftests/globalcanopy_reader.lua +++ b/datasets/opendata/selftests/globalcanopy_reader.lua @@ -21,7 +21,7 @@ local height = 0.0 print(string.format("\n-------------------------------------------------\nMeta (\u{221E}) global canopy height sample POI\n-------------------------------------------------")) -local expResults = {{20.0, 1396483218000, '/vsis3/sliderule/data/GLOBALCANOPY/META_GlobalCanopyHeight_1m_2024_v1.vrt'}} +local expResults = {{20.0, 1396483218, '/vsis3/sliderule/data/GLOBALCANOPY/META_GlobalCanopyHeight_1m_2024_v1.vrt'}} local demType = "meta-globalcanopy-1meter" local dem = geo.raster(geo.parms({ asset = demType, algorithm = "NearestNeighbour", radius = 0})) diff --git a/datasets/opendata/selftests/worldcover_reader.lua b/datasets/opendata/selftests/worldcover_reader.lua index 13d19ab3..695c60fa 100644 --- a/datasets/opendata/selftests/worldcover_reader.lua +++ b/datasets/opendata/selftests/worldcover_reader.lua @@ -20,7 +20,7 @@ local height = 0.0 print(string.format("\n-------------------------------------------------\nesa worldcover 10meter sample POI\n-------------------------------------------------")) -local expResults = {{10.0, 1309046418000, '/vsis3/sliderule/data/WORLDCOVER/ESA_WorldCover_10m_2021_v200_Map.vrt'}} +local expResults = {{10.0, 1309046418, '/vsis3/sliderule/data/WORLDCOVER/ESA_WorldCover_10m_2021_v200_Map.vrt'}} local demType = "esa-worldcover-10meter" local dem = geo.raster(geo.parms({ asset = demType, algorithm = "NearestNeighbour", radius = 0})) diff --git a/datasets/pgc/package/PgcDemStripsRaster.cpp b/datasets/pgc/package/PgcDemStripsRaster.cpp index c60b9471..68583c55 100644 --- a/datasets/pgc/package/PgcDemStripsRaster.cpp +++ b/datasets/pgc/package/PgcDemStripsRaster.cpp @@ -245,24 +245,24 @@ bool PgcDemStripsRaster::findRasters(raster_finder_t* finder) if(!fileName.empty()) { raster_info_t flagsRinfo; - flagsRinfo.flagsBandNum = 1; + flagsRinfo.flagsBandNum = 1; flagsRinfo.tag = FLAGS_TAG; flagsRinfo.fileId = finder->fileDict.add(fileName); rgroup->infovect.push_back(flagsRinfo); } } - double gps = 0; + double gps_msecs = 0; for(const auto& s: dates) { TimeLib::gmt_time_t gmt; - gps += getGmtDate(feature, s, gmt); + gps_msecs += getGmtDate(feature, s, gmt); } - gps = gps/dates.size(); + gps_msecs = gps_msecs/dates.size(); /* Set raster group time and group id */ - rgroup->gmtDate = TimeLib::gps2gmttime(static_cast(gps)); - rgroup->gpsTime = static_cast(gps); + rgroup->gmtDate = TimeLib::gps2gmttime(static_cast(gps_msecs)); + rgroup->gpsTime = static_cast(gps_msecs / 1000.0); rgroup->infovect.push_back(demRinfo); rgroup->infovect.shrink_to_fit(); finder->rasterGroups.push_back(rgroup); diff --git a/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp b/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp index 291eb7bd..7aaa5770 100644 --- a/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp +++ b/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp @@ -112,7 +112,7 @@ bool Usgs3dep1meterDemRaster::findRasters(raster_finder_t* finder) if (!rastergeo->Intersects(geo)) continue; rasters_group_t* rgroup = new rasters_group_t; - rgroup->gpsTime = getGmtDate(feature, DATE_TAG, rgroup->gmtDate); + rgroup->gpsTime = getGmtDate(feature, DATE_TAG, rgroup->gmtDate) / 1000.0; const char* fname = feature->GetFieldAsString("url"); if(fname && strlen(fname) > 0) diff --git a/packages/geo/GeoIndexedRaster.cpp b/packages/geo/GeoIndexedRaster.cpp index b04588de..b2ace1d4 100644 --- a/packages/geo/GeoIndexedRaster.cpp +++ b/packages/geo/GeoIndexedRaster.cpp @@ -446,23 +446,23 @@ uint32_t GeoIndexedRaster::getBatchGroupSamples(const rasters_group_t* rgroup, L /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ if(ps.bandSample[i] == NULL) continue;; - RasterSample* s; + RasterSample* sample; if(!ps.bandSampleReturned[i]->exchange(true)) { - s = ps.bandSample[i]; + sample = ps.bandSample[i]; } else { /* Sample has already been returned, must create a copy */ - s = new RasterSample(*ps.bandSample[i]); + sample = new RasterSample(*ps.bandSample[i]); } /* Set time for this sample */ - s->time = rgroup->gpsTime / 1000; + sample->time = rgroup->gpsTime; /* Set flags for this sample, add it to the list */ - s->flags = flags; - slist->add(s); + sample->flags = flags; + slist->add(sample); errors |= ps.ssErrors; } @@ -501,12 +501,19 @@ uint32_t GeoIndexedRaster::getBatchGroupFlags(const rasters_group_t* rgroup, uin * The flags value must be in the first band. * If these assumptions are not met the dataset must override this function. */ - RasterSample* s = ps.bandSample[0]; + RasterSample* sample = NULL; + + /* bandSample can be empty if raster failed to open */ + if(!ps.bandSample.empty()) + { + sample = ps.bandSample[0]; + } + /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ - if(s == NULL) break;; + if(sample == NULL) break;; - return s->value; + return sample->value; } } } @@ -591,7 +598,7 @@ uint32_t GeoIndexedRaster::getGroupFlags(const rasters_group_t* rgroup) cacheitem_t* item; const char* key = fileDict.get(rinfo.fileId); - if(cache.find(key, &item)) + if(cache.find(key, &item) && !item->bandSample.empty()) { /* * This function assumes that there is only one raster with FLAGS_TAG in a group. @@ -800,7 +807,7 @@ void GeoIndexedRaster::sampleRasters(OGRGeometry* geo) /*---------------------------------------------------------------------------- * sample *----------------------------------------------------------------------------*/ -bool GeoIndexedRaster::sample(OGRGeometry* geo, int64_t gps, GroupOrdering* groupList) +bool GeoIndexedRaster::sample(OGRGeometry* geo, int64_t gps_secs, GroupOrdering* groupList) { /* Open the index file, if not already open */ if(!openGeoIndex(geo, NULL)) @@ -823,7 +830,7 @@ bool GeoIndexedRaster::sample(OGRGeometry* geo, int64_t gps, GroupOrdering* grou groupList->add(groupList->length(), rgroup); } - if(!filterRasters(gps, groupList, fileDict)) + if(!filterRasters(gps_secs, groupList, fileDict)) return false; uint32_t rasters2sample = 0; @@ -1356,7 +1363,7 @@ bool GeoIndexedRaster::updateCache(uint32_t& rasters2sample, const GroupOrdering item = new cacheitem_t; item->raster = new GdalRaster(parms, key, - static_cast(rgroup->gpsTime / 1000), + static_cast(rgroup->gpsTime), rinfo.fileId, rinfo.elevationBandNum, rinfo.flagsBandNum, @@ -1414,7 +1421,7 @@ bool GeoIndexedRaster::updateCache(uint32_t& rasters2sample, const GroupOrdering /*---------------------------------------------------------------------------- * filterRasters *----------------------------------------------------------------------------*/ -bool GeoIndexedRaster::filterRasters(int64_t gps, GroupOrdering* groupList, RasterFileDictionary& dict) +bool GeoIndexedRaster::filterRasters(int64_t gps_secs, GroupOrdering* groupList, RasterFileDictionary& dict) { /* NOTE: temporal filter is applied in openGeoIndex() */ if(!parms->url_substring.value.empty() || parms->filter_doy_range) @@ -1470,15 +1477,15 @@ bool GeoIndexedRaster::filterRasters(int64_t gps, GroupOrdering* groupList, Rast /* Closest time filter - using raster group time, not individual reaster time */ int64_t closestGps = 0; - if(gps > 0) + if(gps_secs > 0) { /* Caller provided gps time, use it insead of time from params */ - closestGps = gps; + closestGps = gps_secs; } else if (parms->filter_closest_time) { /* Params provided closest time */ - closestGps = TimeLib::gmt2gpstime(parms->closest_time); + closestGps = TimeLib::gmt2gpstime(parms->closest_time) / 1000; } if(closestGps > 0) diff --git a/packages/geo/GeoIndexedRaster.h b/packages/geo/GeoIndexedRaster.h index 3423c37f..81a27e9c 100644 --- a/packages/geo/GeoIndexedRaster.h +++ b/packages/geo/GeoIndexedRaster.h @@ -104,7 +104,7 @@ class GeoIndexedRaster: public RasterObject char* featureId; // stac catalog feature id std::vector infovect; // vector of rasters belonging to the same stac catalog feature TimeLib::gmt_time_t gmtDate; // feature date (can be computed from start/end dates) - int64_t gpsTime; // feature gps time + int64_t gpsTime; // feature gps time in seconds RaserGroup(void): featureId(NULL), gmtDate{0,0,0,0,0,0}, gpsTime(0) {} ~RaserGroup(void) { delete[] featureId; } @@ -200,7 +200,7 @@ class GeoIndexedRaster: public RasterObject static void init (void); static void deinit (void); - uint32_t getSamples (const MathLib::point_3d_t& point, int64_t gps, List& slist, void* param=NULL) final; + uint32_t getSamples (const MathLib::point_3d_t& point, int64_t gps_secs, List& slist, void* param=NULL) final; uint32_t getSamples (const std::vector& points, List& sllist, void* param=NULL) final; uint32_t getSubsets (const MathLib::extent_t& extent, int64_t gps, List& slist, void* param=NULL) final; @@ -307,7 +307,7 @@ class GeoIndexedRaster: public RasterObject bool createBatchReaderThreads(uint32_t rasters2sample); bool updateCache (uint32_t& rasters2sample, const GroupOrdering* groupList); - bool filterRasters (int64_t gps, GroupOrdering* groupList, RasterFileDictionary& dict); + bool filterRasters (int64_t gps_secs, GroupOrdering* groupList, RasterFileDictionary& dict); static OGRGeometry* getConvexHull (const std::vector* points); void applySpatialFilter (OGRLayer* layer, const std::vector* points); diff --git a/packages/geo/RasterObject.cpp b/packages/geo/RasterObject.cpp index 7542fc60..35bc758c 100644 --- a/packages/geo/RasterObject.cpp +++ b/packages/geo/RasterObject.cpp @@ -383,7 +383,7 @@ int RasterObject::luaSamples(lua_State *L) int64_t gps = 0; if(closest_time_str != NULL) { - gps = TimeLib::str2gpstime(closest_time_str); + gps = TimeLib::str2gpstime(closest_time_str) / 1000; } /* Get samples */ From d22e4b9156b646330d32c51059f634a84163036f Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Wed, 20 Nov 2024 20:17:09 +0000 Subject: [PATCH 16/17] fixed time bug, updated arrowsampler point time --- clients/python/sliderule/raster.py | 2 +- clients/python/tests/test_worldcover.py | 19 ++++++++----------- packages/arrow/ArrowSamplerImpl.cpp | 18 ++++++++++++++---- 3 files changed, 23 insertions(+), 16 deletions(-) diff --git a/clients/python/sliderule/raster.py b/clients/python/sliderule/raster.py index e8779b22..f09652d5 100644 --- a/clients/python/sliderule/raster.py +++ b/clients/python/sliderule/raster.py @@ -142,7 +142,7 @@ def sample(asset, coordinates, parms={}): latitude = coordinates[per_coord_index][1] per_coord_index += 1 for raster_sample in input_coord_response: - columns['time'][i] = numpy.int64((raster_sample['time'] + 315564800.0) * 1e9) + columns['time'][i] = numpy.int64((raster_sample['time'] + 315964800.0) * 1e9) columns['longitude'][i] = numpy.double(longitude) columns['latitude'][i] = numpy.double(latitude) columns['file'] += raster_sample['file'], diff --git a/clients/python/tests/test_worldcover.py b/clients/python/tests/test_worldcover.py index d0fb247a..23032448 100644 --- a/clients/python/tests/test_worldcover.py +++ b/clients/python/tests/test_worldcover.py @@ -15,7 +15,8 @@ expValue = 10 vrtFile = '/vsis3/sliderule/data/WORLDCOVER/ESA_WorldCover_10m_2021_v200_Map.vrt' -vrtFileTime = 1309046418 # gpsTime in seconds 2021-06-30 +timeAsDateStr = '2021-06-30 00:00:18' +timeAsUnixSecs = 1309046418 # includes 18 leap seconds @pytest.mark.network class TestMosaic: @@ -25,7 +26,7 @@ def test_vrt(self, init): assert init assert abs(rsps["samples"][0][0]["value"] - expValue) < sigma assert rsps["samples"][0][0]["file"] == vrtFile - assert rsps["samples"][0][0]["time"] == vrtFileTime # datetime in seconds + assert rsps["samples"][0][0]["time"] == timeAsUnixSecs # datetime in seconds def test_time_overflow(self, init): region = [ {"lon": -108.34, "lat": 38.89}, @@ -41,16 +42,12 @@ def test_time_overflow(self, init): # Reset the index to make 'time' a regular column gdf_reset = gdf.reset_index() - # Access raw time column values as NumPy array - raw_time_values = gdf_reset['time'].astype('int64') + # Access the 'time' column + time_values = gdf_reset['time'].tolist() + print(time_values) - # Expected gps epoch time in nanoseconds - expected_time = (vrtFileTime + 315564800) * 1e9 - - print("Raw time column values (NumPy array):", raw_time_values) - print("expected_time:", expected_time) - - assert (raw_time_values == expected_time).all() + assert [str(t) for t in time_values] == [timeAsDateStr, timeAsDateStr], \ + f"Time values do not match. Expected: {[timeAsDateStr, timeAsDateStr]}, Found: {[str(t) for t in time_values]}" values = gdf_reset['value'] assert (values == expValue).all() diff --git a/packages/arrow/ArrowSamplerImpl.cpp b/packages/arrow/ArrowSamplerImpl.cpp index 8c38929c..617edce3 100644 --- a/packages/arrow/ArrowSamplerImpl.cpp +++ b/packages/arrow/ArrowSamplerImpl.cpp @@ -281,13 +281,23 @@ void ArrowSamplerImpl::getPoints(std::vector& points) if(time_column_index > -1) { auto time_column = std::static_pointer_cast(table->column(time_column_index)->chunk(0)); - mlog(DEBUG, "Time column elements: %ld", time_column->length()); + auto timestamp_type = std::static_pointer_cast(table->column(time_column_index)->type()); - /* Update gps time for each point */ + if(timestamp_type->unit() != arrow::TimeUnit::NANO) + { + mlog(ERROR, "Time column must be in nanoseconds."); + points.clear(); + return; + } + + /* Convert unix nanoseconds to gps seconds */ for(int64_t i = 0; i < time_column->length(); i++) - points[i].gps = time_column->Value(i); + { + double unix_nsecs = time_column->Value(i); + points[i].gps = TimeLib::sys2gpstime(unix_nsecs / 1000.0) * 1000.0; + } } - else mlog(DEBUG, "Time column not found."); + else mlog(ERROR, "Time column not found."); } /*---------------------------------------------------------------------------- From 456d222e9858470dc17c81b96675ef680594f3f4 Mon Sep 17 00:00:00 2001 From: Eric Lidwa Date: Wed, 20 Nov 2024 20:37:46 +0000 Subject: [PATCH 17/17] minor fix which stopped the build --- packages/arrow/ArrowSamplerImpl.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/arrow/ArrowSamplerImpl.cpp b/packages/arrow/ArrowSamplerImpl.cpp index 617edce3..93b413ee 100644 --- a/packages/arrow/ArrowSamplerImpl.cpp +++ b/packages/arrow/ArrowSamplerImpl.cpp @@ -293,7 +293,7 @@ void ArrowSamplerImpl::getPoints(std::vector& points) /* Convert unix nanoseconds to gps seconds */ for(int64_t i = 0; i < time_column->length(); i++) { - double unix_nsecs = time_column->Value(i); + const double unix_nsecs = time_column->Value(i); points[i].gps = TimeLib::sys2gpstime(unix_nsecs / 1000.0) * 1000.0; } }