Skip to content

Commit

Permalink
Add convenience GDALDataset::GeolocationToPixelLine() and GDALRasterB…
Browse files Browse the repository at this point in the history 
…and::InterpolateAtGeolocation()

and corresponding C functions GDALDatasetGeolocationToPixelLine() and
GDALRasterInterpolateAtGeolocation()

and Python Band.InterpolateAtGeolocation()
  • Loading branch information
@rouault
rouault committed Jan 10, 2025
1 parent 006a750 commit ce2769c
Show file tree
Hide file tree
Showing 8 changed files with 413 additions and 3 deletions.
81 changes: 80 additions & 1 deletion autotest/gcore/interpolateatpoint.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
import gdaltest
import pytest

from osgeo import gdal
from osgeo import gdal, osr

###############################################################################
# Test some cases of InterpolateAtPoint
Expand Down Expand Up @@ -353,3 +353,82 @@ 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 = osr.SpatialReference()
wgs84_srs.SetFromUserInput("WGS84")
ct = osr.CoordinateTransformation(ds.GetSpatialRef(), wgs84_srs)
wgs84_lat, wgs84_lon, _ = ct.TransformPoint(X, Y, 0)

got_nearest = ds.GetRasterBand(1).InterpolateAtGeolocation(
wgs84_lat, wgs84_lon, wgs84_srs, 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.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
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_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
)
10 changes: 10 additions & 0 deletions gcore/gdal.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand Down
11 changes: 11 additions & 0 deletions gcore/gdal_priv.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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,
Expand Down
127 changes: 127 additions & 0 deletions gcore/gdaldataset.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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);
}
94 changes: 92 additions & 2 deletions gcore/gdalrasterband.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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);
}
Loading

0 comments on commit ce2769c

Please sign in to comment.