diff --git a/autotest/gcore/interpolateatpoint.py b/autotest/gcore/interpolateatpoint.py
index d3ee6add3249..1409ff6a735e 100755
--- a/autotest/gcore/interpolateatpoint.py
+++ b/autotest/gcore/interpolateatpoint.py
@@ -15,7 +15,7 @@
 import gdaltest
 import pytest
 
-from osgeo import gdal
+from osgeo import gdal, osr
 
 ###############################################################################
 # Test some cases of InterpolateAtPoint
@@ -353,3 +353,83 @@ 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, None, gdal.GRIORA_NearestNeighbour
+    )
+    assert got_nearest == pytest.approx(173, 1e-6)
+    got_bilinear = ds.GetRasterBand(1).InterpolateAtGeolocation(
+        X, Y, None, gdal.GRIORA_Bilinear
+    )
+    assert got_bilinear == pytest.approx(139.75, 1e-6)
+    got_cubic = ds.GetRasterBand(1).InterpolateAtGeolocation(
+        X, Y, None, gdal.GRIORA_CubicSpline
+    )
+    assert got_cubic == pytest.approx(138.02, 1e-2)
+    got_cubic = ds.GetRasterBand(1).InterpolateAtGeolocation(
+        X, Y, None, gdal.GRIORA_Cubic
+    )
+    assert got_cubic == pytest.approx(145.57, 1e-2)
+
+    wgs84_srs_auth_compliant = osr.SpatialReference()
+    wgs84_srs_auth_compliant.SetFromUserInput("WGS84")
+    wgs84_srs_auth_compliant.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
+    ct = osr.CoordinateTransformation(ds.GetSpatialRef(), wgs84_srs_auth_compliant)
+    wgs84_lat, wgs84_lon, _ = ct.TransformPoint(X, Y, 0)
+
+    got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation(
+        wgs84_lat, wgs84_lon, wgs84_srs_auth_compliant, gdal.GRIORA_NearestNeighbour
+    )
+    assert got_nearest == pytest.approx(173, 1e-6)
+
+    wgs84_trad_gis_order = osr.SpatialReference()
+    wgs84_trad_gis_order.SetFromUserInput("WGS84")
+    wgs84_trad_gis_order.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation(
+        wgs84_lon, wgs84_lat, wgs84_trad_gis_order, gdal.GRIORA_NearestNeighbour
+    )
+    assert got_nearest == pytest.approx(173, 1e-6)
+
+    wgs84_lon_lat_srs = osr.SpatialReference()
+    wgs84_lon_lat_srs.SetFromUserInput("WGS84")
+    wgs84_lon_lat_srs.SetDataAxisToSRSAxisMapping([2, 1])
+    got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation(
+        wgs84_lon, wgs84_lat, wgs84_lon_lat_srs, gdal.GRIORA_NearestNeighbour
+    )
+    assert got_nearest == pytest.approx(173, 1e-6)
+
+    wgs84_lat_lon_srs = osr.SpatialReference()
+    wgs84_lat_lon_srs.SetFromUserInput("WGS84")
+    wgs84_lat_lon_srs.SetDataAxisToSRSAxisMapping([1, 2])
+    got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation(
+        wgs84_lat, wgs84_lon, wgs84_lat_lon_srs, gdal.GRIORA_NearestNeighbour
+    )
+    assert got_nearest == pytest.approx(173, 1e-6)
+
+    X = gt[0] + -1 * gt[1]
+    Y = gt[3] + -1 * gt[5]
+    assert (
+        ds.GetRasterBand(1).InterpolateAtGeolocation(
+            X, Y, None, 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, None, gdal.GRIORA_NearestNeighbour
+        )
diff --git a/doc/source/spelling_wordlist.txt b/doc/source/spelling_wordlist.txt
index f39a583c6268..5a9138f49ad7 100644
--- a/doc/source/spelling_wordlist.txt
+++ b/doc/source/spelling_wordlist.txt
@@ -1065,6 +1065,8 @@ geoloc
 geolocated
 geolocation
 Geolocation
+geolocX
+geolocY
 Géologiques
 Geomatics
 Geomedia
diff --git a/gcore/gdal.h b/gcore/gdal.h
index 8f83b5d3bc9b..fea2e37f7d10 100644
--- a/gcore/gdal.h
+++ b/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, OGRSpatialReferenceH hSRS,
+    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,12 @@ CPLErr CPL_DLL GDALRasterInterpolateAtPoint(GDALRasterBandH hBand,
                                             double *pdfRealValue,
                                             double *pdfImagValue);
 
+CPLErr CPL_DLL GDALRasterInterpolateAtGeolocation(
+    GDALRasterBandH hBand, double dfGeolocX, double dfGeolocY,
+    OGRSpatialReferenceH hSRS, GDALRIOResampleAlg eInterpolation,
+    double *pdfRealValue, double *pdfImagValue,
+    CSLConstList papszTransformerOptions);
+
 /** Generic pointer for the working structure of VRTProcessedDataset
  * function. */
 typedef void *VRTPDWorkingDataPtr;
diff --git a/gcore/gdal_priv.h b/gcore/gdal_priv.h
index 5ec00f099baf..6405ebe9f6e3 100644
--- a/gcore/gdal_priv.h
+++ b/gcore/gdal_priv.h
@@ -715,6 +715,11 @@ class CPL_DLL GDALDataset : public GDALMajorObject
     virtual CPLErr GetGeoTransform(double *padfTransform);
     virtual CPLErr SetGeoTransform(double *padfTransform);
 
+    CPLErr GeolocationToPixelLine(
+        double dfGeolocX, double dfGeolocY, const OGRSpatialReference *poSRS,
+        double *pdfPixel, double *pdfLine,
+        CSLConstList papszTransformerOptions = nullptr) const;
+
     virtual CPLErr AddBand(GDALDataType eType, char **papszOptions = nullptr);
 
     virtual void *GetInternalHandle(const char *pszHandleName);
@@ -1874,6 +1879,12 @@ class CPL_DLL GDALRasterBand : public GDALMajorObject
 
     std::shared_ptr<GDALMDArray> AsMDArray() const;
 
+    CPLErr InterpolateAtGeolocation(
+        double dfGeolocX, double dfGeolocY, const OGRSpatialReference *poSRS,
+        GDALRIOResampleAlg eInterpolation, double *pdfRealValue,
+        double *pdfImagValue = nullptr,
+        CSLConstList papszTransfomerOptions = nullptr) const;
+
     virtual CPLErr InterpolateAtPoint(double dfPixel, double dfLine,
                                       GDALRIOResampleAlg eInterpolation,
                                       double *pdfRealValue,
diff --git a/gcore/gdaldataset.cpp b/gcore/gdaldataset.cpp
index a6bf53191298..bbb2343e6f0d 100644
--- a/gcore/gdaldataset.cpp
+++ b/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,129 @@ GDALDataset::Clone(int nScopeFlags, [[maybe_unused]] bool bCanShareState) const
 }
 
 //! @endcond
+
+/************************************************************************/
+/*                    GeolocationToPixelLine()                          */
+/************************************************************************/
+
+/** Transform georeferenced coordinates to pixel/line coordinates.
+ *
+ * When poSRS is null, those georeferenced coordinates (dfGeolocX, dfGeolocY)
+ * must be in the "natural" 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.
+ * If that natural SRS is a geographic one, dfGeolocX must be a longitude, and
+ * dfGeolocY a latitude. If that natural SRS is a projected one, dfGeolocX must
+ * be a easting, and dfGeolocY a northing.
+ *
+ * When poSRS is set to a non-null value, (dfGeolocX, dfGeolocY) must be
+ * expressed in that CRS, and that tuple must be conformant with the
+ * data-axis-to-crs-axis setting of poSRS, that is the one returned by
+ * the OGRSpatialReference::GetDataAxisToSRSAxisMapping(). If you want to be sure
+ * of the axis order, then make sure to call poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER)
+ * before calling this method, and in that case, dfGeolocX must be a longitude
+ * or an easting value, and dfGeolocX a latitude or a northing value.
+ *
+ * This method uses GDALCreateGenImgProjTransformer2() underneath.
+ *
+ * @param dfGeolocX X coordinate of the position (longitude or easting if poSRS
+ * is null, otherwise consistent with poSRS data-axis-to-crs-axis mapping),
+ * where interpolation should be done.
+ * @param dfGeolocY Y coordinate of the position (latitude or northing if poSRS
+ * is null, otherwise consistent with poSRS data-axis-to-crs-axis mapping),
+ * where interpolation should be done.
+ * @param poSRS If set, override the natural CRS in which dfGeolocX, dfGeolocY are expressed
+ * @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,
+                                    const OGRSpatialReference *poSRS,
+                                    double *pdfPixel, double *pdfLine,
+                                    CSLConstList papszTransformerOptions) const
+{
+    CPLStringList aosTO(papszTransformerOptions);
+
+    if (poSRS)
+    {
+        const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
+        const std::string osWKT = poSRS->exportToWkt(apszOptions);
+        aosTO.SetNameValue("DST_SRS", osWKT.c_str());
+        const auto eAxisMappingStrategy = poSRS->GetAxisMappingStrategy();
+        if (eAxisMappingStrategy == OAMS_TRADITIONAL_GIS_ORDER)
+            aosTO.SetNameValue("DST_SRS_AXIS_MAPPING_STRATEGY",
+                               "TRADITIONAL_GIS_ORDER");
+        else if (eAxisMappingStrategy == OAMS_AUTHORITY_COMPLIANT)
+            aosTO.SetNameValue("DST_SRS_AXIS_MAPPING_STRATEGY",
+                               "AUTHORITY_COMPLIANT");
+        else
+        {
+            const auto &anValues = poSRS->GetDataAxisToSRSAxisMapping();
+            std::string osVal;
+            for (int v : anValues)
+            {
+                if (!osVal.empty())
+                    osVal += ',';
+                osVal += std::to_string(v);
+            }
+            aosTO.SetNameValue("DST_SRS_DATA_AXIS_TO_SRS_AXIS_MAPPING",
+                               osVal.c_str());
+        }
+    }
+
+    auto hTransformer = GDALCreateGenImgProjTransformer2(
+        GDALDataset::ToHandle(const_cast<GDALDataset *>(this)), nullptr,
+        aosTO.List());
+    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,
+                                         OGRSpatialReferenceH hSRS,
+                                         double *pdfPixel, double *pdfLine,
+                                         CSLConstList papszTransformerOptions)
+{
+    VALIDATE_POINTER1(hDS, "GDALDatasetGeolocationToPixelLine", CE_Failure);
+
+    GDALDataset *poDS = GDALDataset::FromHandle(hDS);
+    return poDS->GeolocationToPixelLine(
+        dfGeolocX, dfGeolocY, OGRSpatialReference::FromHandle(hSRS), pdfPixel,
+        pdfLine, papszTransformerOptions);
+}
diff --git a/gcore/gdalrasterband.cpp b/gcore/gdalrasterband.cpp
index 187f5dddb872..4185b65cab29 100644
--- a/gcore/gdalrasterband.cpp
+++ b/gcore/gdalrasterband.cpp
@@ -9662,8 +9662,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.
@@ -9726,3 +9726,93 @@ 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.
+ *
+ * When poSRS is null, those georeferenced coordinates (dfGeolocX, dfGeolocY)
+ * must be in the "natural" 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.
+ * If that natural SRS is a geographic one, dfGeolocX must be a longitude, and
+ * dfGeolocY a latitude. If that natural SRS is a projected one, dfGeolocX must
+ * be a easting, and dfGeolocY a northing.
+ *
+ * When poSRS is set to a non-null value, (dfGeolocX, dfGeolocY) must be
+ * expressed in that CRS, and that tuple must be conformant with the
+ * data-axis-to-crs-axis setting of poSRS, that is the one returned by
+ * the OGRSpatialReference::GetDataAxisToSRSAxisMapping(). If you want to be sure
+ * of the axis order, then make sure to call poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER)
+ * before calling this method, and in that case, dfGeolocX must be a longitude
+ * or an easting value, and dfGeolocX a latitude or a northing value.
+ *
+ * The GDALDataset::GeolocationToPixelLine() will be used to transform from
+ * (dfGeolocX,dfGeolocY) 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 if poSRS
+ * is null, otherwise consistent with poSRS data-axis-to-crs-axis mapping),
+ * where interpolation should be done.
+ * @param dfGeolocY Y coordinate of the position (latitude or northing if poSRS
+ * is null, otherwise consistent with poSRS data-axis-to-crs-axis mapping),
+ * where interpolation should be done.
+ * @param poSRS If set, override the natural CRS in which dfGeolocX, dfGeolocY are expressed
+ * @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, const OGRSpatialReference *poSRS,
+    GDALRIOResampleAlg eInterpolation, double *pdfRealValue,
+    double *pdfImagValue, CSLConstList papszTransformerOptions) const
+{
+    double dfPixel;
+    double dfLine;
+    if (poDS->GeolocationToPixelLine(dfGeolocX, dfGeolocY, poSRS, &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,
+                                          OGRSpatialReferenceH hSRS,
+                                          GDALRIOResampleAlg eInterpolation,
+                                          double *pdfRealValue,
+                                          double *pdfImagValue,
+                                          CSLConstList papszTransformerOptions)
+{
+    VALIDATE_POINTER1(hBand, "GDALRasterInterpolateAtGeolocation", CE_Failure);
+
+    GDALRasterBand *poBand = GDALRasterBand::FromHandle(hBand);
+    return poBand->InterpolateAtGeolocation(
+        dfGeolocX, dfGeolocY, OGRSpatialReference::FromHandle(hSRS),
+        eInterpolation, pdfRealValue, pdfImagValue, papszTransformerOptions);
+}
diff --git a/swig/include/Band.i b/swig/include/Band.i
index 66d809e166c3..d2eac7f958e7 100644
--- a/swig/include/Band.i
+++ b/swig/include/Band.i
@@ -668,6 +668,29 @@ 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,
+                                   OSRSpatialReferenceShadow* srs,
+                                   GDALRIOResampleAlg interpolation,
+                                   double *pdfRealValue,
+                                   double *pdfImagValue, char** transformerOptions = NULL ) {
+    if (pdfRealValue) *pdfRealValue = 0;
+    if (pdfImagValue) *pdfImagValue = 0;
+    return GDALRasterInterpolateAtGeolocation( self, geolocX, geolocY,
+                (OGRSpatialReferenceH)srs, 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};
diff --git a/swig/include/java/gdal_java.i b/swig/include/java/gdal_java.i
index 9efea1534fb8..70637de4a2ba 100644
--- a/swig/include/java/gdal_java.i
+++ b/swig/include/java/gdal_java.i
@@ -1122,6 +1122,7 @@ import java.awt.Color;
 
 %typemap(javaimports) GDALRasterBandShadow %{
 import org.gdal.gdalconst.gdalconstConstants;
+import org.gdal.osr.SpatialReference;
 %}
 
 
diff --git a/swig/include/python/gdal_python.i b/swig/include/python/gdal_python.i
index 06eb319fa67c..bfd4a429dd90 100644
--- a/swig/include/python/gdal_python.i
+++ b/swig/include/python/gdal_python.i
@@ -5024,6 +5024,75 @@ 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`.
+
+       When srs is None, those georeferenced coordinates (geolocX, geolocY)
+       must be in the "natural" 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.
+       If that natural SRS is a geographic one, geolocX must be a longitude, and
+       geolocY a latitude. If that natural SRS is a projected one, geolocX must
+       be a easting, and geolocY a northing.
+
+       When srs is set to a non-None value, (geolocX, geolocY) must be
+       expressed in that CRS, and that tuple must be conformant with the
+       data-axis-to-crs-axis setting of srs, that is the one returned by
+       the :py:func:`osgeo.osr.SpatialReference.GetDataAxisToSRSAxisMapping().
+       If you want to be sure of the axis order, then make sure to call
+       ``srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)``
+       before calling this method, and in that case, geolocX must be a longitude
+       or an easting value, and geolocX a latitude or a northing value.
+
+       Parameters
+       ----------
+       geolocX : float
+           X coordinate of the position where interpolation should be done.
+           Longitude or easting in "natural" CRS if `srs` is None,
+           otherwise consistent with first axis of `srs`,
+           taking into account the data-axis-to-crs-axis mapping
+       geolocY : float
+           Y coordinate of the position where interpolation should be done.
+           Latitude or northing in "natural" CRS if `srs` is None,
+           otherwise consistent with second axis of `srs`,
+           taking into account the data-axis-to-crs-axis mapping
+       srs : osgeo.osr.SpatialReference
+           If set, override the natural CRS in which geolocX, geolocY are expressed
+       interpolation : GRIOResampleAlg (nearest, bilinear, cubic, cubicspline)
+
+       Returns
+       -------
+       float:
+           Interpolated value, or ``None`` if it has any error.
+
+       Example
+       -------
+
+       >>> from osgeo import gdal, osr
+       >>> with gdal.Open("my.tif") as ds:
+       ...    wgs84_srs = osr.SpatialReference()
+       ...    wgs84_srs.SetFromUserInput("WGS84")
+       ...    wgs84_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+       ...    val = ds.GetRasterBand(1).InterpolateAtGeolocation(longitude_degree,
+                                                                 latitude_degree,
+                                                                 wgs84_srs,
+                                                                 gdal.GRIORA_Bilinear)
+    """
+
+    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.