gdallocationinfo: Allow querying raster corner points (OSGeo#12087)

* fix read out of bounds in GDALInterpExtractValuesWindow

* InterpolateAtPoint: exclude points just out of range with NearestNeighbour

* InterpolateAtPoint: allow points right at the bottom or right edge with NearestNeighbour

* InterpolateAtPoint: add a small margin to avoid numeric precision issues

* gdallocationinfo: Allow querying raster corner points

* InterpolateAtGeoLocation: add test for corner points

---------

Co-authored-by: Javier Jimenez Shaw <j1@jimenezshaw.com>
Co-authored-by: Chris Toney <jctoney@gmail.com>
Co-authored-by: Even Rouault <even.rouault@spatialys.com>

Author: 4 people

alg/gdal_interpolateatpoint.cpp

@@ -124,13 +124,20 @@ bool GDALInterpExtractValuesWindow(GDALRasterBand *pBand,
             // Compose the cached block to the final buffer
             for (int j = 0; j < nLinesToCopy; j++)
             {
-                memcpy(padfAsDouble + ((nFirstLineInOutput + j) * nWidth +
-                                       nFirstColInOutput) *
-                                          nTypeFactor,
-                       poValue->data() +
-                           ((nFirstLineInCachedBlock + j) * nReqXSize +
-                            nFirstColInCachedBlock) *
-                               nTypeFactor,
+                const size_t dstOffset =
+                    (static_cast<size_t>(nFirstLineInOutput + j) * nWidth +
+                     nFirstColInOutput) *
+                    nTypeFactor;
+                const size_t srcOffset =
+                    (static_cast<size_t>(nFirstLineInCachedBlock + j) *
+                         nReqXSize +
+                     nFirstColInCachedBlock) *
+                    nTypeFactor;
+                if (srcOffset + nColsToCopy * nTypeFactor > poValue->size())
+                {
+                    return false;
+                }
+                memcpy(padfAsDouble + dstOffset, poValue->data() + srcOffset,
                        nColsToCopy * sizeof(T));
             }
         }
@@ -158,9 +165,20 @@ template <typename T>
 bool GDALInterpolateAtPointImpl(GDALRasterBand *pBand,
                                 GDALRIOResampleAlg eResampleAlg,
                                 std::unique_ptr<DoublePointsCache> &cache,
-                                const double dfXIn, const double dfYIn, T &out)
+                                double dfXIn, double dfYIn, T &out)
 {
     const gdal::Vector2i rasterSize{pBand->GetXSize(), pBand->GetYSize()};
+
+    if (eResampleAlg == GRIORA_NearestNeighbour)
+    {
+        // Allow input coordinates right at the bottom or right edge
+        // with GRIORA_NearestNeighbour.
+        // "introduce" them in the pixel of the image.
+        if (dfXIn >= rasterSize.x() && dfXIn <= rasterSize.x() + 1e-5)
+            dfXIn -= 0.25;
+        if (dfYIn >= rasterSize.y() && dfYIn <= rasterSize.y() + 1e-5)
+            dfYIn -= 0.25;
+    }
     const gdal::Vector2d inLoc{dfXIn, dfYIn};
 
     int bGotNoDataValue = FALSE;

apps/gdallocationinfo.cpp

@@ -439,8 +439,9 @@ MAIN_START(argc, argv)
 
         bool bPixelReport = true;
 
-        if (iPixel < 0 || iLine < 0 || iPixel >= GDALGetRasterXSize(hSrcDS) ||
-            iLine >= GDALGetRasterYSize(hSrcDS))
+        if (iPixel < 0 || iLine < 0 ||
+            dfPixel > static_cast<double>(GDALGetRasterXSize(hSrcDS) + 1e-5) ||
+            dfLine > static_cast<double>(GDALGetRasterYSize(hSrcDS) + 1e-5))
         {
             if (bAsXML)
                 osXML += "<Alert>Location is off this file! No further details "

autotest/gcore/interpolateatpoint.py

@@ -44,6 +44,25 @@ def test_interpolateatpoint_outofrange():
     assert res is None
     res = ds.GetRasterBand(1).InterpolateAtPoint(10, 1200, gdal.GRIORA_Bilinear)
     assert res is None
+    res = ds.GetRasterBand(1).InterpolateAtPoint(-1, 0, gdal.GRIORA_NearestNeighbour)
+    assert res is None
+    res = ds.GetRasterBand(1).InterpolateAtPoint(0, -0.5, gdal.GRIORA_NearestNeighbour)
+    assert res is None
+
+
+def test_interpolateatpoint_right_at_border():
+
+    ds = gdal.Open("data/byte.tif")
+    res = ds.GetRasterBand(1).InterpolateAtPoint(20, 20, gdal.GRIORA_NearestNeighbour)
+    assert res == pytest.approx(107.0, 1e-5)
+    res = ds.GetRasterBand(1).InterpolateAtPoint(18, 20, gdal.GRIORA_NearestNeighbour)
+    assert res == pytest.approx(99.0, 1e-5)
+    res = ds.GetRasterBand(1).InterpolateAtPoint(20, 18, gdal.GRIORA_NearestNeighbour)
+    assert res == pytest.approx(123.0, 1e-5)
+    res = ds.GetRasterBand(1).InterpolateAtPoint(20, 20, gdal.GRIORA_Bilinear)
+    assert res == pytest.approx(107.0, 1e-5)
+    res = ds.GetRasterBand(1).InterpolateAtPoint(20.1, 20.1, gdal.GRIORA_Bilinear)
+    assert res is None
 
 
 @gdaltest.disable_exceptions()
@@ -444,3 +463,23 @@ def test_interpolateatgeolocation_1():
         assert ungeoreferenced_ds.GetRasterBand(1).InterpolateAtGeolocation(
             0, 0, None, gdal.GRIORA_NearestNeighbour
         )
+
+
+def test_interpolateatgeolocation_corners():
+
+    gdaltest.importorskip_gdal_array()
+
+    ds = gdal.Open("data/byte.tif")
+    band = ds.GetRasterBand(1)
+    data = band.ReadAsArray()
+
+    values = {}
+    for loc, coord in gdal.Info(ds, format="json")["cornerCoordinates"].items():
+        values[loc] = band.InterpolateAtGeolocation(
+            coord[0], coord[1], None, gdal.GRIORA_NearestNeighbour
+        )
+
+    assert values["upperLeft"] == data[0, 0]
+    assert values["upperRight"] == data[0, -1]
+    assert values["lowerLeft"] == data[-1, 0]
+    assert values["lowerRight"] == data[-1, -1]

autotest/utilities/test_gdallocationinfo.py

@@ -101,6 +101,20 @@ def test_gdallocationinfo_4(gdallocationinfo_path):
     assert ret.startswith(expected_ret)
 
 
+# Test -geoloc at lower right corner
+def test_gdallocationinfo_lr(gdallocationinfo_path):
+
+    ret = gdaltest.runexternal(
+        gdallocationinfo_path + " -geoloc ../gcore/data/byte.tif 441920.000 3750120.000"
+    )
+    ret = ret.replace("\r\n", "\n")
+    expected_ret = """Report:
+  Location: (20P,20L)
+  Band 1:
+    Value: 107"""
+    assert ret.startswith(expected_ret)
+
+
 ###############################################################################
 # Test -lifonly