Bluetopo support (#447)

Implemented bluetopo, updated gebco, fixed bugs

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()

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'],

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

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

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:

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"

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,16 +13,41 @@
 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:
     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()

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
+        $<INSTALL_INTERFACE:${INCDIR}>
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>/package
+)
+
+install (
+    FILES
+        ${CMAKE_CURRENT_LIST_DIR}/package/bluetopo.h
+        ${CMAKE_CURRENT_LIST_DIR}/package/BlueTopoBathyRaster.h
+    DESTINATION
+        ${INCDIR}/bluetopo
+)

datasets/bluetopo/bluetopo.md

@@ -0,0 +1 @@
+Raster sampling support for BlueTopo (The compilation of the nation's best available bathymetric data.)