diff --git a/CMakeLists.txt b/CMakeLists.txt index 1e3cf1e9f..cf6d733e7 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/clients/python/sliderule/raster.py b/clients/python/sliderule/raster.py index ca70ae9fd..f09652d57 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'] + 315964800.0) * 1000000000) + 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_bluetopo.py b/clients/python/tests/test_bluetopo.py new file mode 100644 index 000000000..119564e19 --- /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_20241115.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_gebco.py b/clients/python/tests/test_gebco.py index 5398f66f7..314bb8e46 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 diff --git a/clients/python/tests/test_globalcanopy.py b/clients/python/tests/test_globalcanopy.py index 4378fec35..8854e863b 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_landsat.py b/clients/python/tests/test_landsat.py index a374f0055..819d3a288 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/clients/python/tests/test_worldcover.py b/clients/python/tests/test_worldcover.py index 6a36e0df6..230324484 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 @@ -11,9 +13,10 @@ vrtLon = -108.1 vrtLat = 39.1 -vrtValue = 10 +expValue = 10 vrtFile = '/vsis3/sliderule/data/WORLDCOVER/ESA_WorldCover_10m_2021_v200_Map.vrt' -vrtFileTime = 1309046418000 +timeAsDateStr = '2021-06-30 00:00:18' +timeAsUnixSecs = 1309046418 # includes 18 leap seconds @pytest.mark.network class TestMosaic: @@ -21,6 +24,30 @@ 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"] == timeAsUnixSecs # datetime in seconds + + 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) + + # Reset the index to make 'time' a regular column + gdf_reset = gdf.reset_index() + + # Access the 'time' column + time_values = gdf_reset['time'].tolist() + print(time_values) + + 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/datasets/bluetopo/CMakeLists.txt b/datasets/bluetopo/CMakeLists.txt new file mode 100644 index 000000000..3c5186abe --- /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 000000000..10347982f --- /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 000000000..bfee12bbf --- /dev/null +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.cpp @@ -0,0 +1,274 @@ +/* + * 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" + +/* Valid Bands for BlueTopo Bathy Raster */ +const char* BlueTopoBathyRaster::validBands[] = {"Elevation", "Uncertainty", "Contributor"}; +#define validBandsCnt (sizeof (validBands) / sizeof (const char *)) + +/****************************************************************************** + * 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()) +{ + if(!validateBandNames()) + { + throw RunTimeException(DEBUG, RTE_ERROR, "Invalid band name specified"); + } + + const 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) / 1000.0; + + 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(); + + 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; + } + + raster_info_t rinfo; + rinfo.elevationBandNum = 1; + 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 + ******************************************************************************/ + +/*---------------------------------------------------------------------------- + * validateBandNames + *----------------------------------------------------------------------------*/ +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 */ + 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 + *----------------------------------------------------------------------------*/ +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 */ + const std::string token = ".gpkg"; + + for(int i = 0; fileList[i] != nullptr; i++) + { + const std::string fileName = fileList[i]; + + if(fileName.size() >= token.size() && + fileName.compare(fileName.size() - token.size(), token.size(), token) == 0) + { + indexFile = bucketPath; + indexFile.append("/").append(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 000000000..6a81daea5 --- /dev/null +++ b/datasets/bluetopo/package/BlueTopoBathyRaster.h @@ -0,0 +1,100 @@ +/* + * 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 + *--------------------------------------------------------------------*/ + static const char* validBands[]; + + /*-------------------------------------------------------------------- + * 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 validateBandNames (void); + 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 000000000..52a6b50ca --- /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 000000000..4073d4943 --- /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 000000000..d1663d5b0 --- /dev/null +++ b/datasets/bluetopo/selftests/bluetopo_reader.lua @@ -0,0 +1,61 @@ +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 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; + +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 })) +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 + band = v["band"] + value = v["value"] + fname = v["file"] + 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)) + 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 + +-- Report Results -- + +runner.report() + + diff --git a/datasets/bluetopo/selftests/plugin_unittest.lua b/datasets/bluetopo/selftests/plugin_unittest.lua new file mode 100644 index 000000000..0e55c8871 --- /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/gebco/package/GebcoBathyRaster.cpp b/datasets/gebco/package/GebcoBathyRaster.cpp index ac35cb7c3..a06e6c41d 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()); } @@ -86,15 +112,15 @@ 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)) { 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 +130,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/gebco/selftests/gebco_reader.lua b/datasets/gebco/selftests/gebco_reader.lua index 8e24d124d..14846ff63 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 000000000..7f4ed0c32 --- /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 85fe03313..b6e369d1c 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 4175b9737..78431b501 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/datasets/gedi/package/Gedi03Raster.h b/datasets/gedi/package/Gedi03Raster.h index 8ddddadef..e392f97ba 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 ed0264705..b273e5379 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 c62847077..78de2718b 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; ifeatureId = 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++) @@ -177,7 +182,7 @@ bool LandsatHlsRaster::findRasters(raster_finder_t* finder) const char* bandName = it->first.c_str(); /* skip algo names (NDIS, etc) */ - if(isValidAlgoName(bandName)) + if(validAlgoName(bandName)) continue; const char* fname = feature->GetFieldAsString(bandName); @@ -187,12 +192,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) { @@ -225,6 +230,33 @@ bool LandsatHlsRaster::findRasters(raster_finder_t* finder) * PRIVATE METHODS ******************************************************************************/ +/*---------------------------------------------------------------------------- + * validateBandNames + *----------------------------------------------------------------------------*/ +bool LandsatHlsRaster::validateBandNames(void) +{ + 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]; + 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 true; +} + +/*---------------------------------------------------------------------------- + * validateBand + *----------------------------------------------------------------------------*/ bool LandsatHlsRaster::validateBand(band_type_t type, const char* bandName) { const char** tags; @@ -291,6 +323,12 @@ 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; + + /* Samples to be returned to the user */ + std::vector sampleVect; + /* Collect samples for all rasters */ if(mode == SERIAL) { @@ -298,41 +336,43 @@ 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->sample; - 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) - { - slist->add(sample); - item->sample = 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; + sampleVect.push_back(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; + } } } } @@ -349,8 +389,16 @@ uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_gr { if(ps.pointIndex == pointIndx) { + 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(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(); @@ -362,40 +410,43 @@ 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 = sample; } else { /* Sample has already been returned, must create a copy */ - s = new RasterSample(*ps.sample); + 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; - slist->add(s); + sampleVect.push_back(s); errors |= ps.ssErrors; } } 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; } @@ -407,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 */ @@ -417,7 +468,10 @@ 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; - slist->add(sample); + + /* Set band name for this sample */ + sample->bandName = "NDSI"; + sampleVect.push_back(sample); } if(ndvi) @@ -426,7 +480,10 @@ 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; - slist->add(sample); + + /* Set band name for this sample */ + sample->bandName = "NDVI"; + sampleVect.push_back(sample); } if(ndwi) @@ -435,7 +492,23 @@ 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; - slist->add(sample); + + /* Set band name for this sample */ + sample->bandName = "NDWI"; + 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 575ef7299..2b27e2be2 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 getInnerBands (std::vector& bands) final + { bands.clear(); } /* Landsat bands are in seperate rasters */ private: @@ -107,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); @@ -123,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 6b1a74146..bb5fa62e5 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/datasets/opendata/package/EsaWorldCover10meterRaster.h b/datasets/opendata/package/EsaWorldCover10meterRaster.h index c5d69a702..02eba9ded 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) / 1000 /* Mid point for year data was collected */) {} private: diff --git a/datasets/opendata/package/MetaGlobalCanopy1meterRaster.h b/datasets/opendata/package/MetaGlobalCanopy1meterRaster.h index 4b1298975..36233562b 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) / 1000) {} }; diff --git a/datasets/opendata/selftests/globalcanopy_reader.lua b/datasets/opendata/selftests/globalcanopy_reader.lua index 5c9bd46cc..ac7d941e3 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 13d19ab30..695c60fae 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/ArcticDemMosaicRaster.h b/datasets/pgc/package/ArcticDemMosaicRaster.h index 2822717b6..5a772f7a8 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 c434c9651..68583c551 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,24 +245,24 @@ 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); } } - 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/pgc/package/RemaDemMosaicRaster.h b/datasets/pgc/package/RemaDemMosaicRaster.h index 29e98d9d5..0a6af1249 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 3a65aef58..7aaa5770d 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) @@ -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/arrow/ArrowSamplerImpl.cpp b/packages/arrow/ArrowSamplerImpl.cpp index 11f47789e..93b413ee0 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); + { + const 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."); } /*---------------------------------------------------------------------------- @@ -434,6 +444,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 +481,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 +490,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 +520,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 +545,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 +570,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 +590,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 +613,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 +649,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 +664,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 +686,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 +710,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 +735,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 +755,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 +778,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/GdalRaster.cpp b/packages/geo/GdalRaster.cpp index c2e718e1f..bac1c50ae 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,8 +66,10 @@ GdalRaster::GdalRaster(const GeoFields* _parms, const std::string& _fileName, do overrideCRS(cb), fileName (_fileName), dset (NULL), - band (NULL), - dataIsElevation(_dataIsElevation), + elevationBandNum(_elevationBandNum), + elevationBand(NULL), + flagsBandNum(_flagsBandNum), + flagsBand (NULL), xsize (0), ysize (0), cellSize (0), @@ -109,6 +114,42 @@ 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 band numbers*/ + 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); + } + } + + /* 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(); @@ -140,9 +181,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 +195,9 @@ void GdalRaster::open(void) dset = NULL; OGRCoordinateTransformation::DestroyCT(transf); transf = NULL; - band = NULL; + bandMap.clear(); + elevationBand = NULL; + flagsBand = NULL; throw; } } @@ -165,7 +205,7 @@ void GdalRaster::open(void) /*---------------------------------------------------------------------------- * samplePOI *----------------------------------------------------------------------------*/ -RasterSample* GdalRaster::samplePOI(OGRPoint* poi) +RasterSample* GdalRaster::samplePOI(OGRPoint* poi, int bandNum) { RasterSample* sample = NULL; @@ -177,11 +217,14 @@ 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()); + 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. @@ -191,13 +234,22 @@ 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); + + if(band == flagsBand) + { + /* Skip resampling and zonal stats for quality mask band (value is bitmask) */ + readPixel(poi, band, sample); + } else - resamplePixel(poi, sample); + { + if(parms->sampling_algo == GRIORA_NearestNeighbour) + readPixel(poi, band, sample); + else + resamplePixel(poi, band, sample); - if(parms->zonal_stats) - computeZonalStats(poi, sample); + if(parms->zonal_stats) + computeZonalStats(poi, band, sample); + } } else { @@ -218,7 +270,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 +393,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 +411,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 +447,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 +495,21 @@ 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(ERROR, "Band \"%s\" not found", bandName.c_str()); + return NO_BAND; + } + return it->second; +} + + /*---------------------------------------------------------------------------- * setCRSfromWkt *----------------------------------------------------------------------------*/ @@ -527,7 +596,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 +609,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 */ @@ -579,91 +649,93 @@ void GdalRaster::readPixel(const OGRPoint* poi, RasterSample* sample) /* 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) && 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); } @@ -678,7 +750,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,16 +799,17 @@ 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) && band == elevationBand) { sample->value += sample->verticalShift; } + sample->bandName = band->GetDescription(); } else { /* At least return pixel value if unable to resample raster */ - readPixel(poi, sample); + readPixel(poi, band, sample); } } catch (const RunTimeException &e) @@ -749,7 +822,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 +844,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; /* @@ -793,7 +867,7 @@ void GdalRaster::computeZonalStats(const OGRPoint* poi, RasterSample* sample) 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 */ @@ -872,7 +946,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 +1055,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 +1083,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 +1101,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 6136ce442..d3819bf3b 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 @@ -82,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 @@ -101,21 +103,27 @@ 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); - 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; } 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); /*-------------------------------------------------------------------- * Static Methods @@ -145,8 +153,10 @@ class GdalRaster std::string fileName; GDALDataset *dset; - GDALRasterBand* band; - bool dataIsElevation; + int elevationBandNum; + GDALRasterBand* elevationBand; + int flagsBandNum; + GDALRasterBand* flagsBand; uint32_t xsize; uint32_t ysize; double cellSize; @@ -157,19 +167,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); + 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); - 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 8f0e6570a..b2ace1d41 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 *----------------------------------------------------------------------------*/ @@ -183,15 +212,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); } @@ -275,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]; + } } } @@ -402,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)) + for(size_t i = 0; i < ps.bandSample.size(); i++) { - s = ps.sample; - } - else - { - /* 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* sample; + if(!ps.bandSampleReturned[i]->exchange(true)) + { + sample = ps.bandSample[i]; + } + else + { + /* Sample has already been returned, must create a copy */ + sample = new RasterSample(*ps.bandSample[i]); + } + + /* Set time for this sample */ + sample->time = rgroup->gpsTime; - /* 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 */ + sample->flags = flags; + slist->add(sample); + 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; } @@ -454,11 +496,24 @@ uint32_t GeoIndexedRaster::getBatchGroupFlags(const rasters_group_t* rgroup, uin { if(ps.pointIndex == pointIndx) { - /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ - if(ps.sample) + /* + * 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* sample = NULL; + + /* bandSample can be empty if raster failed to open */ + if(!ps.bandSample.empty()) { - return ps.sample->value; + sample = ps.bandSample[0]; } + + + /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ + if(sample == NULL) break;; + + return sample->value; } } } @@ -480,21 +535,24 @@ 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 */ ssErrors |= item->raster->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; } } @@ -512,11 +570,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 */ @@ -537,9 +598,14 @@ 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()) { - const RasterSample* _sample = item->sample; + /* + * 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) { return _sample->value; @@ -685,7 +751,7 @@ bool GeoIndexedRaster::openGeoIndex(const OGRGeometry* geo, const std::vectoradd(groupList->length(), rgroup); } - if(!filterRasters(gps, groupList, fileDict)) + if(!filterRasters(gps_secs, groupList, fileDict)) return false; uint32_t rasters2sample = 0; @@ -896,33 +962,80 @@ void* GeoIndexedRaster::readerThread(void *param) cacheitem_t* entry = reader->entry; if(entry != NULL) { - if(GdalRaster::ispoint(reader->geo)) - { - entry->sample = entry->raster->samplePOI(dynamic_cast(reader->geo)); - } - else if(GdalRaster::ispoly(reader->geo)) + GdalRaster* raster = entry->raster; + + try { - entry->subset = entry->raster->subsetAOI(dynamic_cast(reader->geo)); - if(entry->subset) + CHECKPTR(raster); + + /* Open raster so we can get inner bands from it */ + raster->open(); + + std::vector bands; + reader->obj->getInnerBands(raster, bands); + + if(GdalRaster::ispoint(reader->geo)) { - /* - * 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, + /* 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 + { + /* 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) + { + /* 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, - entry->subset->rasterName, - entry->raster->getGpsTime(), - entry->raster->isElevation(), - entry->raster->getOverrideCRS()); + subset->rasterName, + raster->getGpsTime(), + raster->getElevationBandNum(), + raster->getFLagsBandNum(), + raster->getOverrideCRS()); - /* RequestFields are shared with subsseted raster and other readers */ - GeoIndexedRaster::referenceLuaObject(reader->obj->rqstParms); + 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(); @@ -957,18 +1070,56 @@ 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->fileId), - 0, /* Sample collecting code will set it to group's gpsTime */ - ur->fileId, - ur->dataIsElevation, - breader->obj->crscb); - - /* Sample all points for this raster */ - for(point_sample_t& ps : ur->pointSamples) + GdalRaster* raster = NULL; + + try { - ps.sample = raster->samplePOI(&ps.point); - ps.ssErrors |= raster->getSSerror(); + 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) + { + /* 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 + { + /* 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; @@ -1210,16 +1361,22 @@ 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, - static_cast(rgroup->gpsTime / 1000), + item->raster = new GdalRaster(parms, + key, + static_cast(rgroup->gpsTime), rinfo.fileId, - rinfo.dataIsElevation, crscb, &bbox); - item->sample = NULL; - item->subset = NULL; + rinfo.elevationBandNum, + rinfo.flagsBandNum, + crscb, + &bbox); 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++; @@ -1264,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) @@ -1320,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) @@ -1569,7 +1726,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 +1743,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) @@ -1782,7 +1939,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; -} - - - +} \ No newline at end of file diff --git a/packages/geo/GeoIndexedRaster.h b/packages/geo/GeoIndexedRaster.h index cfeae8017..81a27e9c1 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; @@ -93,12 +90,13 @@ class GeoIndexedRaster: public RasterObject /** Raster information needed for sampling */ typedef struct RasterInfo { - bool dataIsElevation; - std::string tag; // "Dem", "Flags", "Date" - 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 */ @@ -106,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; } @@ -114,11 +112,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; @@ -136,10 +132,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; @@ -204,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; @@ -228,7 +224,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; @@ -311,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); @@ -327,7 +323,6 @@ class GeoIndexedRaster: public RasterObject bool collectSamples (const std::vector& pointsGroups, List& sllist); - }; #endif /* __geo_indexed_raster__ */ diff --git a/packages/geo/GeoJsonRaster.cpp b/packages/geo/GeoJsonRaster.cpp index 1953f4ff8..decfcd724 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), @@ -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 23ac319c1..2ff9193c8 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); @@ -86,9 +87,16 @@ 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); - RasterSample* sample = raster.samplePOI(&ogr_point); - if(sample) slist.add(sample); + std::vector bands; + 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); + } } catch (const RunTimeException &e) { @@ -116,27 +124,33 @@ 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); - 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()); - - /* RequestFields are shared with subsseted raster */ - referenceLuaObject(rqstParms); - slist.add(subset); + /* 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.getElevationBandNum(), + raster.getFLagsBandNum(), + raster.getOverrideCRS()); + slist.add(subset); + + /* RequestFields are shared with subsseted raster */ + referenceLuaObject(rqstParms); + } } } catch (const RunTimeException &e) @@ -155,7 +169,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 +179,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 8134f9d16..1b0f3a92e 100644 --- a/packages/geo/GeoRaster.h +++ b/packages/geo/GeoRaster.h @@ -61,11 +61,12 @@ 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; - 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/GeoUserRaster.cpp b/packages/geo/GeoUserRaster.cpp index 976f358e8..2f8643093 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 85a26d3af..7b3fd7854 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 be066279d..35bc758cb 100644 --- a/packages/geo/RasterObject.cpp +++ b/packages/geo/RasterObject.cpp @@ -252,16 +252,65 @@ 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); + } +} + +/*---------------------------------------------------------------------------- + * 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 *----------------------------------------------------------------------------*/ @@ -334,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 */ @@ -388,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); diff --git a/packages/geo/RasterObject.h b/packages/geo/RasterObject.h index 7e8e03c58..98512c5ae 100644 --- a/packages/geo/RasterObject.h +++ b/packages/geo/RasterObject.h @@ -107,9 +107,17 @@ 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); + 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 + { + return parms->bands.length() > 0; + } + bool hasZonalStats (void) const { return parms->zonal_stats; @@ -175,7 +183,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); /*-------------------------------------------------------------------- diff --git a/packages/geo/RasterSample.h b/packages/geo/RasterSample.h index c6b09d876..9e92b9ac5 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/packages/geo/UT_RasterSample.cpp b/packages/geo/UT_RasterSample.cpp index 2b393e055..898742827 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/platforms/linux/SlideRule.cpp b/platforms/linux/SlideRule.cpp index e1ac0a8b3..99d8638c5 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/scripts/selftests/test_runner.lua b/scripts/selftests/test_runner.lua index 339da13c0..ce413e4cc 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,13 @@ 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.. "plugin_unittest.lua") + runner.script(bluetopo_td.. "bluetopo_reader.lua") +end + -- Run Default Parameters Self Tests for all modules runner.script(td .. "parms_tojson.lua") diff --git a/targets/slideruleearth-aws/asset_directory.csv b/targets/slideruleearth-aws/asset_directory.csv index 3bae45c4d..998c1ba83 100644 --- a/targets/slideruleearth-aws/asset_directory.csv +++ b/targets/slideruleearth-aws/asset_directory.csv @@ -19,7 +19,8 @@ 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 arcticdem-strips, nil, nil, /vsis3/pgc-opendata-dems/arcticdem/strips/s2s041/2m, nil, us-west-2, https://s3.us-west-2.amazonaws.com