added luaBatchSample, updated bluetopo tests

clients/python/tests/test_bluetopo.py

@@ -8,14 +8,13 @@
 
 sigma = 1.0e-9
 
-lon = -80.87
-lat =  32.06
+lon =  -81.02
+lat =   31.86
 
-file  = '/vsis3/noaa-ocs-nationalbathymetry-pds/BlueTopo/BH4SX59B/BlueTopo_BH4SX59B_20241203.tiff'
-
-expElevation   = -4.14
-expUncertainty =  1.03
-expContributor = 63824
+file  = '/vsis3/noaa-ocs-nationalbathymetry-pds/BlueTopo/BH4SV597/BlueTopo_BH4SV597_20241219.tiff'
+expElevation   = -14.10
+expUncertainty =   2.58
+expContributor =  63846
 
 @pytest.mark.network
 class TestBlueTopo:

datasets/bluetopo/selftests/bluetopo_reader.lua

@@ -11,13 +11,13 @@ 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}
+local lons = {-81.02, -89.66, -94.72}
+local lats = { 31.86,  29.99,  29.35}
 height = 0
 
-local expElevation   = { -4.14, -14.10, -4.28, -17.18}
-local expUncertainty = {  1.03,   2.58,  0.34,  1.32}
-local expContributor = { 63824,  63846, 24955, 45641}
+local expElevation   = {-14.10, -4.28, -17.18}
+local expUncertainty = {  2.58,  0.34,  1.32}
+local expContributor = { 63846, 24955, 45641}
 
 local elevation_tolerance = 0.01;
 

datasets/icesat2/package/MeritRaster.cpp

@@ -112,26 +112,25 @@ MeritRaster::MeritRaster(lua_State *L, RequestFields* rqst_parms, const char* ke
 /*----------------------------------------------------------------------------
  * getSamples
  *----------------------------------------------------------------------------*/
-uint32_t MeritRaster::getSamples (const MathLib::point_3d_t& point, int64_t gps, List<RasterSample*>& slist, void* param)
+uint32_t MeritRaster::getSamples (const point_info_t& pinfo, sample_list_t& slist, void* param)
 {
     (void)param;
-    (void)gps;
 
     lockSampling();
 
     /* Determine Upper Left Coordinates */
-    int left_lon = ((int)floor(point.x / 5.0)) * 5;
-    int upper_lat = ((int)ceil(point.y / 5.0)) * 5;
+    int left_lon = ((int)floor(pinfo.point3d.x / 5.0)) * 5;
+    int upper_lat = ((int)ceil(pinfo.point3d.y / 5.0)) * 5;
 
     /* Calculate Pixel Location */
-    const int x_offset = (int)(((double)point.x - left_lon) / X_SCALE);
-    const int y_offset = (int)(((double)point.y - upper_lat) / Y_SCALE);
+    const int x_offset = (int)(((double)pinfo.point3d.x - left_lon)  / X_SCALE);
+    const int y_offset = (int)(((double)pinfo.point3d.y - upper_lat) / Y_SCALE);
 
     /* Check Pixel Location */
     if( x_offset < 0 || x_offset >= X_MAX ||
         y_offset < 0 || y_offset >= Y_MAX )
     {
-        mlog(ERROR, "Invalid pixel location for MERIT DEM at %lf, %lf: %d, %d\n", point.x, point.y, x_offset, y_offset);
+        mlog(ERROR, "Invalid pixel location for MERIT DEM at %lf, %lf: %d, %d\n", pinfo.point3d.x, pinfo.point3d.y, x_offset, y_offset);
         return SS_OUT_OF_BOUNDS_ERROR;
     }
 

datasets/icesat2/package/MeritRaster.h

@@ -77,7 +77,7 @@ class MeritRaster: public RasterObject
          *--------------------------------------------------------------------*/
 
                  MeritRaster (lua_State *L, RequestFields* rqst_parms, const char* key);
-        uint32_t getSamples (const MathLib::point_3d_t& point, int64_t gps, List<RasterSample*>& slist, void* param=NULL) final;
+        uint32_t getSamples (const point_info_t& pinfo, sample_list_t& slist, void* param=NULL) final;
         uint32_t getSubsets (const MathLib::extent_t&  extent, int64_t gps, List<RasterSubset*>& slist, void* param=NULL) final;
     private:
 

datasets/pgc/package/PgcDemStripsRaster.cpp

@@ -180,10 +180,10 @@ void PgcDemStripsRaster::getIndexFile(const std::vector<point_info_t>* points, s
 
     /* Get all geojson files for all points */
     std::vector<std::string> files;
-    for(const auto& p : *points)
+    for(const auto& pinfo : *points)
     {
         std::string newFile;
-        _getIndexFile(p.point.x, p.point.y, newFile);
+        _getIndexFile(pinfo.point3d.x, pinfo.point3d.y, newFile);
         if(!newFile.empty())
         {
             files.push_back(newFile);

datasets/pgc/selftests/arcticdem_reader_batch.lua

@@ -0,0 +1,81 @@
+local runner = require("test_executive")
+local console = require("console")
+local asset = require("asset")
+local csv = require("csv")
+local json = require("json")
+
+-- console.monitor:config(core.LOG, core.DEBUG)
+-- sys.setlvl(core.LOG, core.DEBUG)
+
+local assets = asset.loaddir()
+
+-- Unit Test --
+
+local demTypes = {"arcticdem-mosaic", "arcticdem-strips"}
+
+-- Correct values test for different POIs
+
+local sigma = 1.0e-9
+
+local lons = {-150, -40, 100, 150}
+local lats = {  70,  70,  70,  75}
+local heights = {0,   0,   0,   0}
+
+local expResultsMosaic = { 116.250000000,  2968.085937500, 475.742187500, 19.890625000}
+local expResultsStrips = { 119.554687500,  2968.015625000, 474.156250000, 10.296875000}  -- Only first strip samples for each lon/lat strip group
+
+local expSamplesCnt = {8, 8, 4, 14}
+
+for i = 1, #demTypes do
+
+    local demType = demTypes[i];
+    print(string.format("\n--------------------------------\nTest: %s Correct Values\n--------------------------------", demType))
+    local dem = geo.raster(geo.parms({ asset = demType, algorithm = "NearestNeighbour", sort_by_index=true}))
+    local tbl, err = dem:batchsample(lons, lats, heights)
+
+    if err ~= 0 then
+        print("FAILED to batchsample")
+        runner.check(false)
+    else
+        for j, pointSamples in ipairs(tbl) do
+            local lon, lat = lons[j], lats[j]
+            local sampleCnt = #pointSamples
+            print(string.format("Point: %d, (%.3f, %.3f) Sample Count: %d", j, lon, lat, sampleCnt))
+
+            if sampleCnt > 0 then
+                for k, sample in ipairs(pointSamples) do
+                    local el = sample["value"]
+                    local fname = sample["file"]
+                    print(string.format("    Sample (%02d): %16.9fm   %s", k, el, fname))
+
+                    if k == 1 then -- Check all mosaics but only first strip for each POI
+                        local expElevation
+                        if demType == "arcticdem-mosaic" then
+                            expElevation = expResultsMosaic[j]
+                        else
+                            expElevation = expResultsStrips[j]
+                        end
+                        runner.check(math.abs(el - expElevation) < sigma)
+                    else
+                        runner.check(el > 10)  -- Check all other samples
+                    end
+                end
+            else
+                print(string.format("Point: %d, (%.3f, %.3f) ======> FAILED to read samples", j, lon, lat))
+            end
+
+            local expectedSamplesCnt
+            if demType == "arcticdem-mosaic" then
+                expectedSamplesCnt = 1
+            else
+                expectedSamplesCnt = expSamplesCnt[j]
+            end
+            runner.check(sampleCnt == expectedSamplesCnt)
+        end
+    end
+end
+
+-- Report Results --
+
+runner.report()
+

packages/arrow/ArrowSamplerImpl.cpp

@@ -323,7 +323,7 @@ void ArrowSamplerImpl::getXYPoints(std::vector<point_info_t>& points)
         const double x = x_column->Value(i);
         const double y = y_column->Value(i);
 
-        points.emplace_back(point_info_t({{x, y, 0.0}, 0.0}));
+        points.emplace_back(point_info_t({{x, y, 0.0}, 0}));
     }
     mlog(DEBUG, "Read %zu points from file", points.size());
 }
@@ -352,7 +352,7 @@ void ArrowSamplerImpl::getGeoPoints(std::vector<point_info_t>& points)
         const std::string wkb_data = binary_array->GetString(i);     /* Get WKB data as string (binary data) */
         const ArrowCommon::wkbpoint_t point = convertWKBToPoint(wkb_data);
 
-        points.emplace_back(point_info_t({{point.x, point.y, 0.0}, 0.0}));
+        points.emplace_back(point_info_t({{point.x, point.y, 0.0}, 0}));
     }
     mlog(INFO, "Read %zu geo points from file", points.size());
 }

packages/geo/GeoIndexedRaster.h

@@ -200,7 +200,7 @@ class GeoIndexedRaster: public RasterObject
 
         static void     init              (void);
         static void     deinit            (void);
-        uint32_t        getSamples        (const MathLib::point_3d_t& point, int64_t gps_secs, List<RasterSample*>& slist, void* param=NULL) final;
+        uint32_t        getSamples        (const point_info_t& pinfo, sample_list_t& slist, void* param=NULL) final;
         uint32_t        getSamples        (const std::vector<point_info_t>& points, List<sample_list_t*>& sllist, void* param=NULL) final;
         uint32_t        getSubsets        (const MathLib::extent_t&  extent, int64_t gps, List<RasterSubset*>& slist, void* param=NULL) final;
 

packages/geo/GeoIndexedRasterBatch.cpp

@@ -446,7 +446,7 @@ void* GeoIndexedRaster::groupsFinderThread(void *param)
         }
 
         const point_info_t& pinfo = gf->points->at(i);
-        const OGRPoint ogrPoint(pinfo.point.x, pinfo.point.y, pinfo.point.z);
+        const OGRPoint ogrPoint(pinfo.point3d.x, pinfo.point3d.y, pinfo.point3d.z);
 
         /* Query the R-tree with the OGRPoint and get the result features */
         std::vector<OGRFeature*> foundFeatures;
@@ -955,10 +955,10 @@ OGRGeometry* GeoIndexedRaster::getConvexHull(const std::vector<point_info_t>* po
     mlog(INFO, "Creating convex hull from %zu points", points->size());
 
     /* Collect all points into a geometry collection */
-    for(const point_info_t& point_info : *points)
+    for(const point_info_t& pinfo : *points)
     {
-        const double lon = point_info.point.x;
-        const double lat = point_info.point.y;
+        const double lon = pinfo.point3d.x;
+        const double lat = pinfo.point3d.y;
 
         OGRPoint* point = new OGRPoint(lon, lat);
         geometryCollection.addGeometryDirectly(point);

packages/geo/GeoIndexedRasterSerial.cpp

@@ -78,20 +78,20 @@ GeoIndexedRaster::Reader::~Reader (void)
 /*----------------------------------------------------------------------------
  * getSamples - serial sampling
  *----------------------------------------------------------------------------*/
-uint32_t GeoIndexedRaster::getSamples(const MathLib::point_3d_t& point, int64_t gps, List<RasterSample*>& slist, void* param)
+uint32_t GeoIndexedRaster::getSamples(const point_info_t& pinfo, sample_list_t& slist, void* param)
 {
     static_cast<void>(param);
 
     lockSampling();
     try
     {
         GroupOrdering groupList;
-        OGRPoint      ogrPoint(point.x, point.y, point.z);
+        OGRPoint      ogrPoint(pinfo.point3d.x, pinfo.point3d.y, pinfo.point3d.z);
 
         ssErrors = SS_NO_ERRORS;
 
         /* Sample Rasters */
-        if(serialSample(&ogrPoint, gps, &groupList))
+        if(serialSample(&ogrPoint, pinfo.gps, &groupList))
         {
             /* Populate Return List of Samples (slist) */
             const GroupOrdering::Iterator iter(groupList);

packages/geo/GeoRaster.cpp

@@ -79,9 +79,8 @@ GeoRaster::~GeoRaster(void) = default;
 /*----------------------------------------------------------------------------
  * getSamples
  *----------------------------------------------------------------------------*/
-uint32_t GeoRaster::getSamples(const MathLib::point_3d_t& point, int64_t gps, List<RasterSample*>& slist, void* param)
+uint32_t GeoRaster::getSamples(const point_info_t& pinfo, sample_list_t& slist, void* param)
 {
-    static_cast<void>(gps);
     static_cast<void>(param);
 
     lockSampling();
@@ -92,7 +91,7 @@ uint32_t GeoRaster::getSamples(const MathLib::point_3d_t& point, int64_t gps, Li
         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);
+            OGRPoint ogr_point(pinfo.point3d.x, pinfo.point3d.y, pinfo.point3d.z);
 
             RasterSample* sample = raster.samplePOI(&ogr_point, bandNum);
             if(sample) slist.add(sample);

packages/geo/GeoRaster.h

@@ -64,7 +64,7 @@ class GeoRaster: public RasterObject
                       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<RasterSample*>& slist, void* param=NULL) final;
+        uint32_t      getSamples (const point_info_t& pinfo, sample_list_t& slist, void* param=NULL) final;
         uint32_t      getSubsets (const MathLib::extent_t&  extent, int64_t gps, List<RasterSubset*>& slist, void* param=NULL) final;
         uint8_t*      getPixels  (uint32_t ulx, uint32_t uly, uint32_t xsize=0, uint32_t ysize=0, int bandNum=1, void* param=NULL) override;