Add convenience GDALDataset::GeolocationToPixelLine() and GDALRasterB…

…and::InterpolateAtGeolocation()

and corresponding C functions GDALDatasetGeolocationToPixelLine() and
GDALRasterInterpolateAtGeolocation()

and Python Band.InterpolateAtGeolocation()

autotest/gcore/interpolateatpoint.py

@@ -353,3 +353,44 @@ def test_interpolateatpoint_big_complex():
     assert res == pytest.approx((182 + 122j), 1e-6)
     res = ds.GetRasterBand(1).InterpolateAtPoint(255.5, 63.5, gdal.GRIORA_Cubic)
     assert res == pytest.approx((182 + 122j), 1e-6)
+
+
+###############################################################################
+# Test InterpolateAtGeolocation
+
+
+def test_interpolateatgeolocation_1():
+
+    ds = gdal.Open("data/byte.tif")
+
+    gt = ds.GetGeoTransform()
+    X = gt[0] + 10 * gt[1]
+    Y = gt[3] + 12 * gt[5]
+
+    got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation(
+        X, Y, gdal.GRIORA_NearestNeighbour
+    )
+    assert got_nearest == pytest.approx(173, 1e-6)
+    got_bilinear = ds.GetRasterBand(1).InterpolateAtGeolocation(
+        X, Y, gdal.GRIORA_Bilinear
+    )
+    assert got_bilinear == pytest.approx(139.75, 1e-6)
+    got_cubic = ds.GetRasterBand(1).InterpolateAtGeolocation(
+        X, Y, gdal.GRIORA_CubicSpline
+    )
+    assert got_cubic == pytest.approx(138.02, 1e-2)
+    got_cubic = ds.GetRasterBand(1).InterpolateAtGeolocation(X, Y, gdal.GRIORA_Cubic)
+    assert got_cubic == pytest.approx(145.57, 1e-2)
+
+    X = gt[0] + -1 * gt[1]
+    Y = gt[3] + -1 * gt[5]
+    assert (
+        ds.GetRasterBand(1).InterpolateAtGeolocation(X, Y, gdal.GRIORA_NearestNeighbour)
+        is None
+    )
+
+    ungeoreferenced_ds = gdal.GetDriverByName("MEM").Create("", 1, 1)
+    with pytest.raises(Exception, match="Unable to compute a transformation"):
+        assert ungeoreferenced_ds.GetRasterBand(1).InterpolateAtGeolocation(
+            0, 0, gdal.GRIORA_NearestNeighbour
+        )

gcore/gdal.h

@@ -1247,6 +1247,10 @@ CPLErr CPL_DLL GDALSetSpatialRef(GDALDatasetH, OGRSpatialReferenceH);
 CPLErr CPL_DLL CPL_STDCALL GDALGetGeoTransform(GDALDatasetH, double *);
 CPLErr CPL_DLL CPL_STDCALL GDALSetGeoTransform(GDALDatasetH, double *);
 
+CPLErr CPL_DLL GDALDatasetGeolocationToPixelLine(
+    GDALDatasetH, double dfGeolocX, double dfGeolocY, double *pdfPixel,
+    double *pdfLine, CSLConstList papszTransformerOptions);
+
 int CPL_DLL CPL_STDCALL GDALGetGCPCount(GDALDatasetH);
 const char CPL_DLL *CPL_STDCALL GDALGetGCPProjection(GDALDatasetH);
 OGRSpatialReferenceH CPL_DLL GDALGetGCPSpatialRef(GDALDatasetH);
@@ -1703,6 +1707,11 @@ CPLErr CPL_DLL GDALRasterInterpolateAtPoint(GDALRasterBandH hBand,
                                             double *pdfRealValue,
                                             double *pdfImagValue);
 
+CPLErr CPL_DLL GDALRasterInterpolateAtGeolocation(
+    GDALRasterBandH hBand, double dfGeolocX, double dfGeolocY,
+    GDALRIOResampleAlg eInterpolation, double *pdfRealValue,
+    double *pdfImagValue, CSLConstList papszTransformerOptions);
+
 /** Generic pointer for the working structure of VRTProcessedDataset
  * function. */
 typedef void *VRTPDWorkingDataPtr;

gcore/gdal_priv.h

@@ -715,6 +715,10 @@ class CPL_DLL GDALDataset : public GDALMajorObject
     virtual CPLErr GetGeoTransform(double *padfTransform);
     virtual CPLErr SetGeoTransform(double *padfTransform);
 
+    CPLErr GeolocationToPixelLine(
+        double dfGeolocX, double dfGeolocY, double *pdfPixel, double *pdfLine,
+        CSLConstList papszTransformerOptions = nullptr) const;
+
     virtual CPLErr AddBand(GDALDataType eType, char **papszOptions = nullptr);
 
     virtual void *GetInternalHandle(const char *pszHandleName);
@@ -1874,6 +1878,11 @@ class CPL_DLL GDALRasterBand : public GDALMajorObject
 
     std::shared_ptr<GDALMDArray> AsMDArray() const;
 
+    CPLErr InterpolateAtGeolocation(
+        double dfGeolocX, double dfGeolocY, GDALRIOResampleAlg eInterpolation,
+        double *pdfRealValue, double *pdfImagValue = nullptr,
+        CSLConstList papszTransfomerOptions = nullptr) const;
+
     virtual CPLErr InterpolateAtPoint(double dfPixel, double dfLine,
                                       GDALRIOResampleAlg eInterpolation,
                                       double *pdfRealValue,

gcore/gdaldataset.cpp

@@ -36,6 +36,7 @@
 #include "cpl_string.h"
 #include "cpl_vsi.h"
 #include "cpl_vsi_error.h"
+#include "gdal_alg.h"
 #include "ogr_api.h"
 #include "ogr_attrind.h"
 #include "ogr_core.h"
@@ -10301,3 +10302,81 @@ GDALDataset::Clone(int nScopeFlags, [[maybe_unused]] bool bCanShareState) const
 }
 
 //! @endcond
+
+/************************************************************************/
+/*                    GeolocationToPixelLine()                          */
+/************************************************************************/
+
+/** Transform georeferenced coordinates to pixel/line coordinates.
+ *
+ * Those georeferenced coordinates must be in the SRS of the dataset,
+ * that is the one returned by GetSpatialRef() if there is a geotransform,
+ * GetGCPSpatialRef() if there are GCPs, WGS 84 if there are RPC coefficients,
+ * or the SRS of the geolocation array (generally WGS 84) if there is a
+ * geolocation array.
+ *
+ * This method uses GDALCreateGenImgProjTransformer2() underneath.
+ *
+ * @param dfGeolocX X coordinate of the position (longitude or easting), where interpolation should be done.
+ * @param dfGeolocY Y coordinate of the position (latitude or northing), where interpolation should be done.
+ * @param[out] pdfPixel Pointer to the variable where to the store the pixel/column coordinate.
+ * @param[out] pdfLine Pointer to the variable where to the store the line coordinate.
+ * @param papszTransformerOptions Options accepted by GDALCreateGenImgProjTransformer2(), or nullptr.
+ *
+ * @return CE_None on success, or an error code on failure.
+ * @since GDAL 3.11
+ */
+
+CPLErr
+GDALDataset::GeolocationToPixelLine(double dfGeolocX, double dfGeolocY,
+                                    double *pdfPixel, double *pdfLine,
+                                    CSLConstList papszTransformerOptions) const
+{
+    auto hTransformer = GDALCreateGenImgProjTransformer2(
+        GDALDataset::ToHandle(const_cast<GDALDataset *>(this)), nullptr,
+        const_cast<char **>(papszTransformerOptions));
+    if (hTransformer == nullptr)
+    {
+        return CE_Failure;
+    }
+
+    double z = 0;
+    int bSuccess = 0;
+    GDALGenImgProjTransform(hTransformer, TRUE, 1, &dfGeolocX, &dfGeolocY, &z,
+                            &bSuccess);
+    GDALDestroyTransformer(hTransformer);
+    if (bSuccess)
+    {
+        if (pdfPixel)
+            *pdfPixel = dfGeolocX;
+        if (pdfLine)
+            *pdfLine = dfGeolocY;
+        return CE_None;
+    }
+    else
+    {
+        return CE_Failure;
+    }
+}
+
+/************************************************************************/
+/*                  GDALDatasetGeolocationToPixelLine()                 */
+/************************************************************************/
+
+/** Transform georeferenced coordinates to pixel/line coordinates.
+ *
+ * @see GDALDataset::GeolocationToPixelLine()
+ * @since GDAL 3.11
+ */
+
+CPLErr GDALDatasetGeolocationToPixelLine(GDALDatasetH hDS, double dfGeolocX,
+                                         double dfGeolocY, double *pdfPixel,
+                                         double *pdfLine,
+                                         CSLConstList papszTransformerOptions)
+{
+    VALIDATE_POINTER1(hDS, "GDALDatasetGeolocationToPixelLine", CE_Failure);
+
+    GDALDataset *poDS = GDALDataset::FromHandle(hDS);
+    return poDS->GeolocationToPixelLine(dfGeolocX, dfGeolocY, pdfPixel, pdfLine,
+                                        papszTransformerOptions);
+}

gcore/gdalrasterband.cpp

@@ -9668,8 +9668,8 @@ std::shared_ptr<GDALMDArray> GDALRasterBand::AsMDArray() const
 /************************************************************************/
 
 /**
- * \brief Interpolates the value between pixels using
- * a resampling algorithm
+ * \brief Interpolates the value between pixels using a resampling algorithm,
+ * taking pixel/line coordinates as input.
  *
  * @param dfPixel pixel coordinate as a double, where interpolation should be done.
  * @param dfLine line coordinate as a double, where interpolation should be done.
@@ -9732,3 +9732,75 @@ CPLErr GDALRasterInterpolateAtPoint(GDALRasterBandH hBand, double dfPixel,
     return poBand->InterpolateAtPoint(dfPixel, dfLine, eInterpolation,
                                       pdfRealValue, pdfImagValue);
 }
+
+/************************************************************************/
+/*                    InterpolateAtGeolocation()                        */
+/************************************************************************/
+
+/**
+ * \brief Interpolates the value between pixels using a resampling algorithm,
+ * taking georeferenced coordinates as input.
+ *
+ * Those georeferenced coordinates must be in the SRS of the dataset,
+ * that is the one returned by GDALDataset::GetSpatialRef() if there is a
+ * geotransform, GDALDataset::GetGCPSpatialRef() if there are GCPs, WGS 84 if
+ * there are RPC coefficients, or the SRS of the geolocation array (generally
+ * WGS 84) if there is a geolocation array.
+ *
+ * The GDALDataset::GeolocationToPixelLine() will be used to transform from
+ * (X,Y) georeferenced coordinates to (pixel, line). Refer to it for details
+ * on how that transformation is done.
+ *
+ * @param dfGeolocX X coordinate of the position (longitude or easting), where interpolation should be done.
+ * @param dfGeolocY Y coordinate of the position (latitude or northing), where interpolation should be done.
+ * @param eInterpolation interpolation type. Only near, bilinear, cubic and cubicspline are allowed.
+ * @param pdfRealValue pointer to real part of interpolated value
+ * @param pdfImagValue pointer to imaginary part of interpolated value (may be null if not needed)
+ * @param papszTransformerOptions Options accepted by GDALDataset::GeolocationToPixelLine() (GDALCreateGenImgProjTransformer2()), or nullptr.
+ *
+ * @return CE_None on success, or an error code on failure.
+ * @since GDAL 3.11
+ */
+
+CPLErr GDALRasterBand::InterpolateAtGeolocation(
+    double dfGeolocX, double dfGeolocY, GDALRIOResampleAlg eInterpolation,
+    double *pdfRealValue, double *pdfImagValue,
+    CSLConstList papszTransformerOptions) const
+{
+    double dfPixel;
+    double dfLine;
+    if (poDS->GeolocationToPixelLine(dfGeolocX, dfGeolocY, &dfPixel, &dfLine,
+                                     papszTransformerOptions) != CE_None)
+    {
+        return CE_Failure;
+    }
+    return InterpolateAtPoint(dfPixel, dfLine, eInterpolation, pdfRealValue,
+                              pdfImagValue);
+}
+
+/************************************************************************/
+/*                  GDALRasterInterpolateAtGeolocation()                */
+/************************************************************************/
+
+/**
+ * \brief Interpolates the value between pixels using a resampling algorithm,
+ * taking georeferenced coordinates as input.
+ *
+ * @see GDALRasterBand::InterpolateAtGeolocation()
+ * @since GDAL 3.11
+ */
+
+CPLErr GDALRasterInterpolateAtGeolocation(GDALRasterBandH hBand,
+                                          double dfGeolocX, double dfGeolocY,
+                                          GDALRIOResampleAlg eInterpolation,
+                                          double *pdfRealValue,
+                                          double *pdfImagValue,
+                                          CSLConstList papszTransformerOptions)
+{
+    VALIDATE_POINTER1(hBand, "GDALRasterInterpolateAtGeolocation", CE_Failure);
+
+    GDALRasterBand *poBand = GDALRasterBand::FromHandle(hBand);
+    return poBand->InterpolateAtGeolocation(
+        dfGeolocX, dfGeolocY, eInterpolation, pdfRealValue, pdfImagValue,
+        papszTransformerOptions);
+}

swig/include/Band.i

@@ -668,6 +668,26 @@ CPLErr AdviseRead(  int xoff, int yoff, int xsize, int ysize,
 %clear (CPLErr);
 #endif
 
+
+%apply (double *OUTPUT){double *pdfRealValue, double *pdfImagValue};
+#if !defined(SWIGPYTHON)
+%apply (IF_ERROR_RETURN_NONE) { (CPLErr) };
+#endif
+%apply (char **options) { char ** transformerOptions };
+  CPLErr InterpolateAtGeolocation( double geolocX, double geolocY,
+                             GDALRIOResampleAlg interpolation,
+                             double *pdfRealValue,
+                             double *pdfImagValue, char** transformerOptions = NULL ) {
+    if (pdfRealValue) *pdfRealValue = 0;
+    if (pdfImagValue) *pdfImagValue = 0;
+    return GDALRasterInterpolateAtGeolocation( self, geolocX, geolocY, interpolation, pdfRealValue, pdfImagValue, transformerOptions );
+  }
+#if !defined(SWIGPYTHON)
+%clear (CPLErr);
+%clear char ** transformerOptions;
+#endif
+
+
 %apply (double *OUTPUT){double *pdfMin, double *pdfMax};
 %apply (int *OUTPUT){int *pnMinX, int *pnMinY};
 %apply (int *OUTPUT){int *pnMaxX, int *pnMaxY};

swig/include/python/gdal_python.i

@@ -5024,6 +5024,36 @@ def InterpolateAtPoint(self, *args, **kwargs):
         return ret[1]
 %}
 
+%feature("shadow") InterpolateAtGeolocation %{
+def InterpolateAtGeolocation(self, *args, **kwargs):
+    """Return the interpolated value at georeferenced coordinates.
+       See :cpp:func:`GDALRasterBand::InterpolateAtGeolocation`.
+
+       Parameters
+       ----------
+       geolocX : float
+           Longitude/easting
+       geolocY : float
+           Latitude/northing
+       interpolation : GRIOResampleAlg (nearest, bilinear, cubic, cubicspline)
+
+       Returns
+       -------
+       float:
+           Interpolated value, or ``None`` if it has any error.
+    """
+
+    ret = $action(self, *args, **kwargs)
+    if ret[0] != CE_None:
+        return None
+
+    from . import gdal
+    if gdal.DataTypeIsComplex(self.DataType):
+        return complex(ret[1], ret[2])
+    else:
+        return ret[1]
+%}
+
 %feature("shadow") ComputeMinMaxLocation %{
 def ComputeMinMaxLocation(self, *args, **kwargs):
     """Compute the min/max values for a band, and their location.