Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[pull] master from OSGeo:master #155

Merged
merged 10 commits into from
Jan 15, 2025
Merged

[pull] master from OSGeo:master #155

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
006a750
GDALCreateGenImgProjTransformer2(): recognize [SRC_SRS_|DST_SRS_][_DA…
rouault Jan 10, 2025
419d118
Add convenience GDALDataset::GeolocationToPixelLine() and GDALRasterB…
rouault Jan 9, 2025
e734b40
PDF: remove write support for GEO_ENCODING=OGC_BP
rouault Oct 19, 2024
eb94875
CreateGeometryFromGML(): recognize aixm:ElevatedPoint (ignoring the e…
rouault Jan 13, 2025
64de599
ogr2ogr: fix -clipsrc/-clipdst when a input geometry is within the en…
rouault Jan 14, 2025
4f6cdf8
Docker: add muparser dependency [ci skip]
rouault Jan 15, 2025
8a81b58
Merge pull request #11655 from rouault/fix_11652
rouault Jan 15, 2025
8096d0b
Merge pull request #11048 from rouault/pdf_remove_write_ogc_bp
rouault Jan 15, 2025
db84620
Merge pull request #11644 from rouault/fix_11582
rouault Jan 15, 2025
4f6fde6
Merge pull request #11622 from rouault/InterpolateAtGeolocation
rouault Jan 15, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 42 additions & 2 deletions alg/gdaltransformer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1951,30 +1951,70 @@ void *GDALCreateGenImgProjTransformer2(GDALDatasetH hSrcDS, GDALDatasetH hDstDS,

const char *pszCO = CSLFetchNameValue(papszOptions, "COORDINATE_OPERATION");

const auto SetAxisMapping =
[papszOptions](OGRSpatialReference &oSRS, const char *pszPrefix)
{
const char *pszMapping = CSLFetchNameValue(
papszOptions, std::string(pszPrefix)
.append("_DATA_AXIS_TO_SRS_AXIS_MAPPING")
.c_str());
if (pszMapping)
{
CPLStringList aosTokens(CSLTokenizeString2(pszMapping, ",", 0));
std::vector<int> anMapping;
for (int i = 0; i < aosTokens.size(); ++i)
anMapping.push_back(atoi(aosTokens[i]));
oSRS.SetDataAxisToSRSAxisMapping(anMapping);
}
else
{
const char *pszStrategy = CSLFetchNameValueDef(
papszOptions,
std::string(pszPrefix).append("_AXIS_MAPPING_STRATEGY").c_str(),
"TRADITIONAL_GIS_ORDER");
if (EQUAL(pszStrategy, "TRADITIONAL_GIS_ORDER"))
oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
else if (EQUAL(pszStrategy, "AUTHORITY_COMPLIANT"))
oSRS.SetAxisMappingStrategy(OAMS_AUTHORITY_COMPLIANT);
else
{
CPLError(CE_Warning, CPLE_AppDefined,
"Unrecognized value '%s' for %s", pszStrategy,
std::string(pszPrefix)
.append("_AXIS_MAPPING_STRATEGY")
.c_str());
return false;
}
}
return true;
};

OGRSpatialReference oSrcSRS;
if (pszSrcSRS)
{
oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
if (pszSrcSRS[0] != '\0' &&
oSrcSRS.SetFromUserInput(pszSrcSRS) != OGRERR_NONE)
{
CPLError(CE_Failure, CPLE_AppDefined,
"Failed to import coordinate system `%s'.", pszSrcSRS);
return nullptr;
}
if (!SetAxisMapping(oSrcSRS, "SRC_SRS"))
return nullptr;
}

OGRSpatialReference oDstSRS;
if (pszDstSRS)
{
oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
if (pszDstSRS[0] != '\0' &&
oDstSRS.SetFromUserInput(pszDstSRS) != OGRERR_NONE)
{
CPLError(CE_Failure, CPLE_AppDefined,
"Failed to import coordinate system `%s'.", pszDstSRS);
return nullptr;
}
if (!SetAxisMapping(oDstSRS, "DST_SRS"))
return nullptr;
}

const char *pszSrcGeolocArray =
Expand Down
75 changes: 48 additions & 27 deletions apps/ogr2ogr_lib.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -587,12 +587,14 @@ class LayerTranslator
std::unique_ptr<OGRGeometry> m_poClipSrcReprojectedToSrcSRS{};
const OGRSpatialReference *m_poClipSrcReprojectedToSrcSRS_SRS = nullptr;
OGREnvelope m_oClipSrcEnv{};
bool m_bClipSrcIsRectangle = false;

OGRGeometry *m_poClipDstOri = nullptr;
bool m_bWarnedClipDstSRS = false;
std::unique_ptr<OGRGeometry> m_poClipDstReprojectedToDstSRS{};
const OGRSpatialReference *m_poClipDstReprojectedToDstSRS_SRS = nullptr;
OGREnvelope m_oClipDstEnv{};
bool m_bClipDstIsRectangle = false;

bool m_bExplodeCollections = false;
bool m_bNativeData = false;
Expand All @@ -606,10 +608,15 @@ class LayerTranslator
const GDALVectorTranslateOptions *psOptions);

private:
std::pair<const OGRGeometry *, const OGREnvelope *>
GetDstClipGeom(const OGRSpatialReference *poGeomSRS);
std::pair<const OGRGeometry *, const OGREnvelope *>
GetSrcClipGeom(const OGRSpatialReference *poGeomSRS);
struct ClipGeomDesc
{
const OGRGeometry *poGeom = nullptr;
const OGREnvelope *poEnv = nullptr;
bool bGeomIsRectangle = false;
};

ClipGeomDesc GetDstClipGeom(const OGRSpatialReference *poGeomSRS);
ClipGeomDesc GetSrcClipGeom(const OGRSpatialReference *poGeomSRS);
};

static OGRLayer *GetLayerAndOverwriteIfNecessary(GDALDataset *poDstDS,
Expand Down Expand Up @@ -6493,16 +6500,17 @@ bool LayerTranslator::Translate(
if (poStolenGeometry->IsEmpty())
goto end_loop;

const auto [poClipGeom, poClipGeomEnvelope] =
const auto clipGeomDesc =
GetSrcClipGeom(poStolenGeometry->getSpatialReference());

if (poClipGeom && poClipGeomEnvelope)
if (clipGeomDesc.poGeom && clipGeomDesc.poEnv)
{
OGREnvelope oEnv;
poStolenGeometry->getEnvelope(&oEnv);
if (!poClipGeomEnvelope->Contains(oEnv) &&
!(poClipGeomEnvelope->Intersects(oEnv) &&
poClipGeom->Intersects(poStolenGeometry.get())))
if (!clipGeomDesc.poEnv->Contains(oEnv) &&
!(clipGeomDesc.poEnv->Intersects(oEnv) &&
clipGeomDesc.poGeom->Intersects(
poStolenGeometry.get())))
{
goto end_loop;
}
Expand Down Expand Up @@ -6744,22 +6752,23 @@ bool LayerTranslator::Translate(
if (poDstGeometry->IsEmpty())
goto end_loop;

const auto [poClipGeom, poClipGeomEnvelope] =
const auto clipGeomDesc =
GetSrcClipGeom(poDstGeometry->getSpatialReference());

if (!(poClipGeom && poClipGeomEnvelope))
if (!(clipGeomDesc.poGeom && clipGeomDesc.poEnv))
goto end_loop;

OGREnvelope oDstEnv;
poDstGeometry->getEnvelope(&oDstEnv);

if (!poClipGeomEnvelope->Contains(oDstEnv))
if (!(clipGeomDesc.bGeomIsRectangle &&
clipGeomDesc.poEnv->Contains(oDstEnv)))
{
std::unique_ptr<OGRGeometry> poClipped;
if (poClipGeomEnvelope->Intersects(oDstEnv))
if (clipGeomDesc.poEnv->Intersects(oDstEnv))
{
poClipped.reset(
poClipGeom->Intersection(poDstGeometry.get()));
poClipped.reset(clipGeomDesc.poGeom->Intersection(
poDstGeometry.get()));
}
if (poClipped == nullptr || poClipped->IsEmpty())
{
Expand Down Expand Up @@ -6926,23 +6935,25 @@ bool LayerTranslator::Translate(
if (poDstGeometry->IsEmpty())
goto end_loop;

auto [poClipGeom, poClipGeomEnvelope] = GetDstClipGeom(
const auto clipGeomDesc = GetDstClipGeom(
poDstGeometry->getSpatialReference());
if (!poClipGeom || !poClipGeomEnvelope)
if (!clipGeomDesc.poGeom || !clipGeomDesc.poEnv)
{
goto end_loop;
}

OGREnvelope oDstEnv;
poDstGeometry->getEnvelope(&oDstEnv);

if (!poClipGeomEnvelope->Contains(oDstEnv))
if (!(clipGeomDesc.bGeomIsRectangle &&
clipGeomDesc.poEnv->Contains(oDstEnv)))
{
std::unique_ptr<OGRGeometry> poClipped;
if (poClipGeomEnvelope->Intersects(oDstEnv))
if (clipGeomDesc.poEnv->Intersects(oDstEnv))
{
poClipped.reset(poClipGeom->Intersection(
poDstGeometry.get()));
poClipped.reset(
clipGeomDesc.poGeom->Intersection(
poDstGeometry.get()));
}

if (poClipped == nullptr || poClipped->IsEmpty())
Expand Down Expand Up @@ -7151,7 +7162,7 @@ bool LayerTranslator::Translate(
* expressed.
* @return the destination clip geometry and its envelope, or (nullptr, nullptr)
*/
std::pair<const OGRGeometry *, const OGREnvelope *>
LayerTranslator::ClipGeomDesc
LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS)
{
if (m_poClipDstReprojectedToDstSRS_SRS != poGeomSRS)
Expand All @@ -7164,7 +7175,7 @@ LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS)
if (m_poClipDstReprojectedToDstSRS->transformTo(poGeomSRS) !=
OGRERR_NONE)
{
return std::make_pair(nullptr, nullptr);
return ClipGeomDesc();
}
m_poClipDstReprojectedToDstSRS_SRS = poGeomSRS;
}
Expand All @@ -7190,8 +7201,13 @@ LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS)
if (poGeom && !m_oClipDstEnv.IsInit())
{
poGeom->getEnvelope(&m_oClipDstEnv);
m_bClipDstIsRectangle = poGeom->IsRectangle();
}
return std::make_pair(poGeom, poGeom ? &m_oClipDstEnv : nullptr);
ClipGeomDesc ret;
ret.poGeom = poGeom;
ret.poEnv = poGeom ? &m_oClipDstEnv : nullptr;
ret.bGeomIsRectangle = m_bClipDstIsRectangle;
return ret;
}

/************************************************************************/
Expand All @@ -7204,7 +7220,7 @@ LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS)
* expressed.
* @return the source clip geometry and its envelope, or (nullptr, nullptr)
*/
std::pair<const OGRGeometry *, const OGREnvelope *>
LayerTranslator::ClipGeomDesc
LayerTranslator::GetSrcClipGeom(const OGRSpatialReference *poGeomSRS)
{
if (m_poClipSrcReprojectedToSrcSRS_SRS != poGeomSRS)
Expand All @@ -7217,7 +7233,7 @@ LayerTranslator::GetSrcClipGeom(const OGRSpatialReference *poGeomSRS)
if (m_poClipSrcReprojectedToSrcSRS->transformTo(poGeomSRS) !=
OGRERR_NONE)
{
return std::make_pair(nullptr, nullptr);
return ClipGeomDesc();
}
m_poClipSrcReprojectedToSrcSRS_SRS = poGeomSRS;
}
Expand All @@ -7242,8 +7258,13 @@ LayerTranslator::GetSrcClipGeom(const OGRSpatialReference *poGeomSRS)
if (poGeom && !m_oClipSrcEnv.IsInit())
{
poGeom->getEnvelope(&m_oClipSrcEnv);
m_bClipSrcIsRectangle = poGeom->IsRectangle();
}
return std::make_pair(poGeom, poGeom ? &m_oClipSrcEnv : nullptr);
ClipGeomDesc ret;
ret.poGeom = poGeom;
ret.poEnv = poGeom ? &m_oClipSrcEnv : nullptr;
ret.bGeomIsRectangle = m_bClipDstIsRectangle;
return ret;
}

/************************************************************************/
Expand Down
82 changes: 81 additions & 1 deletion autotest/gcore/interpolateatpoint.py
View file
Open in desktop
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,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
)
Binary file modified autotest/gdrivers/data/pdf/test_pdf_composition.pdf
View file
Open in desktop
Binary file not shown.
Binary file modified autotest/gdrivers/data/pdf/test_pdf_composition_libpng_1_6_40.pdf
View file
Open in desktop
Binary file not shown.
Loading
Loading