python raster.sample uses batch sampling, added tests

clients/python/tests/test_arcticdem.py

@@ -5,16 +5,17 @@
 import os.path
 import sliderule
 from sliderule import icesat2
+from sliderule import raster
 
 TESTDIR = Path(__file__).parent
 
 sigma = 1.0e-9
 
 vrtLon = -150.0
 vrtLat =  70.0
-vrtElevation = 116.250000000
+vrtElevation = 116.25
 vrtFile      = '/vsis3/pgc-opendata-dems/arcticdem/mosaics/v4.1/2m_dem_tiles.vrt'
-vrtFileTime  = 1358108640000.0
+vrtFileTime  = 1358108640.0
 
 @pytest.mark.network
 class TestMosaic:
@@ -26,6 +27,19 @@ def test_vrt(self, init):
         assert rsps["samples"][0][0]["file"] ==  vrtFile
         assert rsps["samples"][0][0]["time"] ==  vrtFileTime
 
+    def test_sample_api_serial(self, init):
+        gdf = raster.sample("arcticdem-mosaic", [[vrtLon,vrtLat]])
+        assert init
+        assert len(gdf) == 1
+        assert abs(gdf["value"].iat[0] - vrtElevation) < sigma
+        assert gdf["file"].iat[0] ==  vrtFile
+
+    def test_sample_api_batch(self, init):
+        gdf = raster.sample("arcticdem-mosaic", [[vrtLon,vrtLat],[vrtLon+0.01,vrtLat+0.01]])
+        assert init
+        assert len(gdf) == 2
+        assert abs(gdf["value"].iat[0] - vrtElevation) < sigma
+        assert gdf["file"].iat[0] ==  vrtFile
 
     def test_vrt_with_aoi(self, init):
         bbox = [-179, 50, -177, 52]
@@ -119,3 +133,13 @@ def test_indexed_raster(self, init):
             assert '/pgc-opendata-dems/arcticdem/strips/' in gdf.attrs['file_directory'][file_id]
             assert '_dem.tif' in gdf.attrs['file_directory'][file_id]  # only dems, no flags
 
+    def test_sample_api_serial(self, init):
+        gdf = raster.sample("arcticdem-strips", [[vrtLon,vrtLat]])
+        assert init
+        assert len(gdf) == 8
+
+    def test_sample_api_batch(self, init):
+        gdf = raster.sample("arcticdem-strips", [[vrtLon,vrtLat],[vrtLon+0.01,vrtLat+0.01]])
+        assert init
+        assert len(gdf) == 16
+

clients/python/tests/test_bluetopo.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
 
@@ -30,3 +32,11 @@ def test_samples(self, init):
         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
+
+    def test_sample_api(self, init):
+        gfp = gpd.GeoDataFrame(geometry=gpd.points_from_xy([lon, lon+0.001], [lat, lat+0.001]), crs='EPSG:4326')
+        points = [[x,y] for x,y in zip(gfp.geometry.x , gfp.geometry.y)]
+        gdf = raster.sample("bluetopo-bathy", points, {"bands": ["Elevation", "Uncertainty", "Contributor"]})
+        assert init
+        assert len(gdf) == 6
+

clients/python/tests/test_gebco.py

@@ -3,6 +3,7 @@
 import pytest
 from pathlib import Path
 import sliderule
+from sliderule import raster
 
 TESTDIR = Path(__file__).parent
 
@@ -31,7 +32,6 @@ def test_samples_2023(self, init):
         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
@@ -42,11 +42,17 @@ def test_samples_2024(self, init):
 
     # 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"] == expected_depth_2024
         assert rsps["samples"][0][0]["file"]  == expected_file_2024
         assert rsps["samples"][0][0]["flags"] == expected_tid
+
+    def test_sample_api(self, init):
+        gdf = raster.sample("gebco-bathy", [[lon,lat],[lon+0.01,lat+0.01]])
+        assert init
+        assert len(gdf) == 2
+        assert gdf["value"].iloc[0] == expected_depth_2024
+        assert gdf["file"].iloc[0]  == expected_file_2024

clients/python/tests/test_globalcanopy.py

@@ -3,6 +3,7 @@
 import pytest
 from pathlib import Path
 import sliderule
+from sliderule import raster
 
 TESTDIR = Path(__file__).parent
 
@@ -24,3 +25,18 @@ def test_samples(self, init):
         assert abs(rsps["samples"][0][0]["value"] - vrtValue) < sigma
         assert rsps["samples"][0][0]["file"] ==  vrtFile
         assert rsps["samples"][0][0]["time"] ==  vrtFileTime
+
+    def test_sample_api_serial(self, init):
+        gdf = raster.sample("meta-globalcanopy-1meter", [[vrtLon,vrtLat]])
+        assert init
+        assert len(gdf) == 1
+        assert abs(gdf["value"].iat[0] - vrtValue) < sigma
+        assert gdf["file"].iat[0] ==  vrtFile
+
+    def test_sample_api_batch(self, init):
+        gdf = raster.sample("meta-globalcanopy-1meter", [[vrtLon,vrtLat],[vrtLon+0.01,vrtLat+0.01]])
+        assert init
+        assert len(gdf) == 2
+        assert abs(gdf["value"].iat[0] - vrtValue) < sigma
+        assert gdf["file"].iat[0] ==  vrtFile
+

clients/python/tests/test_rema.py

@@ -4,7 +4,7 @@
 from pathlib import Path
 import os.path
 import sliderule
-from sliderule import icesat2
+from sliderule import raster
 
 TESTDIR = Path(__file__).parent
 
@@ -15,7 +15,7 @@
 
 vrtElevation = 328.0156250
 vrtFile      = '/vsis3/pgc-opendata-dems/rema/mosaics/v2.0/2m/2m_dem_tiles.vrt'
-vrtFileTime  = 1361299922000.0
+vrtFileTime  = 1361299922
 
 @pytest.mark.network
 class TestMosaic:
@@ -44,4 +44,17 @@ def test_vrt_with_proj_pipeline(self, init):
         assert init
         assert abs(rsps["samples"][0][0]["value"] - vrtElevation) < sigma
         assert rsps["samples"][0][0]["file"] ==  vrtFile
-        assert rsps["samples"][0][0]["time"] ==  vrtFileTime
+
+    def test_sample_api_serial(self, init):
+        gdf = raster.sample("rema-mosaic", [[vrtLon,vrtLat]])
+        assert init
+        assert len(gdf) == 1
+        assert abs(gdf["value"].iat[0] - vrtElevation) < sigma
+        assert gdf["file"].iat[0] ==  vrtFile
+
+    def test_sample_api_batch(self, init):
+        gdf = raster.sample("rema-mosaic", [[vrtLon,vrtLat],[vrtLon+0.01,vrtLat+0.01]])
+        assert init
+        assert len(gdf) == 2
+        assert abs(gdf["value"].iat[0] - vrtElevation) < sigma
+        assert gdf["file"].iat[0] ==  vrtFile

clients/python/tests/test_worldcover.py

@@ -29,12 +29,6 @@ def test_vrt(self, init):
         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)

datasets/pgc/package/ArcticDemMosaicRaster.h

@@ -63,7 +63,7 @@ class ArcticDemMosaicRaster: public GeoRaster
         ArcticDemMosaicRaster(lua_State* L, RequestFields* rqst_parms, const char* key):
          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),
+                  TimeLib::datetime2gps(2023, 01, 18, 20, 23, 42) / 1000,
                   1,                   /* elevationBandNum */
                   GdalRaster::NO_BAND, /* maskBandNum      */
                   &overrideTargetCRS) {}

datasets/pgc/package/RemaDemMosaicRaster.h

@@ -64,7 +64,7 @@ class RemaDemMosaicRaster: public GeoRaster
         RemaDemMosaicRaster(lua_State* L, RequestFields* rqst_parms, const char* key):
          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),
+                  TimeLib::datetime2gps(2023, 02, 24, 18, 51, 44) / 1000,
                   1,                   /* elevationBandNum */
                   GdalRaster::NO_BAND, /* flagsBandNum     */
                   &overrideTargetCRS) {}

scripts/endpoints/samples.lua

@@ -1,7 +1,7 @@
 --
 -- ENDPOINT:    /source/samples
 --
--- PURPOSE:     return a sample from a rastar dataset given a lat and lon
+-- PURPOSE:     return a sample from a raster dataset given a lat and lon
 --
 -- INPUT:       rqst
 --              {
@@ -27,14 +27,40 @@ local dem = geo.raster(geo.parms(rqst[geo.PARMS]))
 -- Build Table --
 local samples = {}
 local errors = {}
+
 if dem then
-    for _, position in ipairs(coord) do
-        local lon = position[1]
-        local lat = position[2]
-        local height = position[3] -- optional
+    if #coord == 1 then
+        -- Single coordinate: Use dem:sample
+        local lon = coord[1][1]
+        local lat = coord[1][2]
+        local height = coord[1][3] or 0 -- Default to 0 if not provided
         local sample, error = dem:sample(lon, lat, height)
         table.insert(samples, sample)
-        table.insert(errors, error)
+        table.insert(errors, error or 0)
+    else
+        -- Multiple coordinates: Use dem:batchsample
+        local lons = {}
+        local lats = {}
+        local heights = {}
+
+        for _, position in ipairs(coord) do
+            table.insert(lons, position[1])
+            table.insert(lats, position[2])
+            table.insert(heights, position[3] or 0) -- Default to 0 if not provided
+        end
+
+        local batchSamples, error = dem:batchsample(lons, lats, heights)
+
+        for i, sample in ipairs(batchSamples) do
+            table.insert(samples, sample)
+            table.insert(errors, error) -- Apply the same error code to all samples
+        end
+    end
+else
+    -- Handle Missing DEM
+    for _ in ipairs(coord) do
+        table.insert(samples, nil)
+        table.insert(errors, "DEM not available")
     end
 end