From fa4c932d2a14ab861b1ab15ad7ea6cdaf7d79bf0 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 16:51:12 +0200 Subject: [PATCH 1/2] ogr2ogr: implement GetInverse() in CompositeCT and AxisMappingCoordinateTransformation --- apps/ogr2ogr_lib.cpp | 34 +++++++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/apps/ogr2ogr_lib.cpp b/apps/ogr2ogr_lib.cpp index 16b3539ddd42..42dea3bce9dc 100644 --- a/apps/ogr2ogr_lib.cpp +++ b/apps/ogr2ogr_lib.cpp @@ -43,6 +43,7 @@ #include #include #include +#include #include #include #include @@ -1205,6 +1206,14 @@ class GCPCoordTransformation : public OGRCoordinateTransformation virtual OGRCoordinateTransformation *GetInverse() const override { + static std::once_flag flag; + std::call_once(flag, + []() + { + CPLDebug("OGR2OGR", + "GCPCoordTransformation::GetInverse() " + "called, but not implemented"); + }); return nullptr; } }; @@ -1292,7 +1301,21 @@ class CompositeCT : public OGRCoordinateTransformation virtual OGRCoordinateTransformation *GetInverse() const override { - return nullptr; + if (!poCT1 && !poCT2) + return nullptr; + if (!poCT2) + return poCT1->GetInverse(); + if (!poCT1) + return poCT2->GetInverse(); + auto poInvCT1 = + std::unique_ptr(poCT1->GetInverse()); + auto poInvCT2 = + std::unique_ptr(poCT2->GetInverse()); + if (!poInvCT1 || !poInvCT2) + return nullptr; + return std::make_unique(poInvCT2.release(), true, + poInvCT1.release(), true) + .release(); } }; @@ -1302,9 +1325,14 @@ class CompositeCT : public OGRCoordinateTransformation class AxisMappingCoordinateTransformation : public OGRCoordinateTransformation { - public: bool bSwapXY = false; + explicit AxisMappingCoordinateTransformation(bool bSwapXYIn) + : bSwapXY(bSwapXYIn) + { + } + + public: AxisMappingCoordinateTransformation(const std::vector &mappingIn, const std::vector &mappingOut) { @@ -1360,7 +1388,7 @@ class AxisMappingCoordinateTransformation : public OGRCoordinateTransformation virtual OGRCoordinateTransformation *GetInverse() const override { - return nullptr; + return new AxisMappingCoordinateTransformation(bSwapXY); } }; From ffb6004c25503d80d00575524c44c4adb8bccadd Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 15:40:08 +0200 Subject: [PATCH 2/2] OGRGeometryFactory::transformWithOptions(): deal with polar or anti-meridian discontinuities when going from projected to (any) geographic CRS The logic already existed but was restricted to WGS 84 without any good reason. --- autotest/cpp/test_ogr.cpp | 69 ++++++++++++++++++++ ogr/ogrgeometryfactory.cpp | 125 ++++++++++++++++++++++++------------- 2 files changed, 150 insertions(+), 44 deletions(-) diff --git a/autotest/cpp/test_ogr.cpp b/autotest/cpp/test_ogr.cpp index e75d0b6253b9..8fe4e7469977 100644 --- a/autotest/cpp/test_ogr.cpp +++ b/autotest/cpp/test_ogr.cpp @@ -4201,4 +4201,73 @@ TEST_F(test_ogr, OGRCurve_reversePoints) } } +// Test OGRGeometryFactory::transformWithOptions() +TEST_F(test_ogr, transformWithOptions) +{ + // Projected CRS to national geographic CRS (not including poles or antimeridian) + { + OGRGeometry *poGeom = nullptr; + OGRGeometryFactory::createFromWkt( + "LINESTRING(700000 6600000, 700001 6600001)", nullptr, &poGeom); + ASSERT_NE(poGeom, nullptr); + + OGRSpatialReference oEPSG_2154; + oEPSG_2154.importFromEPSG(2154); // "RGF93 v1 / Lambert-93" + OGRSpatialReference oEPSG_4171; + oEPSG_4171.importFromEPSG(4171); // "RGF93 v1" + oEPSG_4171.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + auto poCT = std::unique_ptr( + OGRCreateCoordinateTransformation(&oEPSG_2154, &oEPSG_4171)); + OGRGeometryFactory::TransformWithOptionsCache oCache; + poGeom = OGRGeometryFactory::transformWithOptions(poGeom, poCT.get(), + nullptr, oCache); + EXPECT_NEAR(poGeom->toLineString()->getX(0), 3, 1e-8); + EXPECT_NEAR(poGeom->toLineString()->getY(0), 46.5, 1e-8); + + delete poGeom; + } + +#ifdef HAVE_GEOS + // Projected CRS to national geographic CRS including antimeridian + { + OGRGeometry *poGeom = nullptr; + OGRGeometryFactory::createFromWkt( + "LINESTRING(657630.64 4984896.17,815261.43 4990738.26)", nullptr, + &poGeom); + ASSERT_NE(poGeom, nullptr); + + OGRSpatialReference oEPSG_6329; + oEPSG_6329.importFromEPSG(6329); // "NAD83(2011) / UTM zone 60N" + OGRSpatialReference oEPSG_6318; + oEPSG_6318.importFromEPSG(6318); // "NAD83(2011)" + oEPSG_6318.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + auto poCT = std::unique_ptr( + OGRCreateCoordinateTransformation(&oEPSG_6329, &oEPSG_6318)); + OGRGeometryFactory::TransformWithOptionsCache oCache; + poGeom = OGRGeometryFactory::transformWithOptions(poGeom, poCT.get(), + nullptr, oCache); + EXPECT_EQ(poGeom->getGeometryType(), wkbMultiLineString); + if (poGeom->getGeometryType() == wkbMultiLineString) + { + const auto poMLS = poGeom->toMultiLineString(); + EXPECT_EQ(poMLS->getNumGeometries(), 2); + if (poMLS->getNumGeometries() == 2) + { + const auto poLS = poMLS->getGeometryRef(0); + EXPECT_EQ(poLS->getNumPoints(), 2); + if (poLS->getNumPoints() == 2) + { + EXPECT_NEAR(poLS->getX(0), 179, 1e-6); + EXPECT_NEAR(poLS->getY(0), 45, 1e-6); + EXPECT_NEAR(poLS->getX(1), 180, 1e-6); + EXPECT_NEAR(poLS->getY(1), 45.004384301691303, 1e-6); + } + } + } + + delete poGeom; + } +#endif +} + } // namespace diff --git a/ogr/ogrgeometryfactory.cpp b/ogr/ogrgeometryfactory.cpp index 90583391997e..532028a25507 100644 --- a/ogr/ogrgeometryfactory.cpp +++ b/ogr/ogrgeometryfactory.cpp @@ -3302,15 +3302,15 @@ static void AlterPole(OGRGeometry *poGeom, OGRPoint *poPole, } /************************************************************************/ -/* IsPolarToWGS84() */ +/* IsPolarToGeographic() */ /* */ /* Returns true if poCT transforms from a projection that includes one */ /* of the pole in a continuous way. */ /************************************************************************/ -static bool IsPolarToWGS84(OGRCoordinateTransformation *poCT, - OGRCoordinateTransformation *poRevCT, - bool &bIsNorthPolarOut) +static bool IsPolarToGeographic(OGRCoordinateTransformation *poCT, + OGRCoordinateTransformation *poRevCT, + bool &bIsNorthPolarOut) { bool bIsNorthPolar = false; bool bIsSouthPolar = false; @@ -3366,14 +3366,14 @@ static bool IsPolarToWGS84(OGRCoordinateTransformation *poCT, } /************************************************************************/ -/* TransformBeforePolarToWGS84() */ +/* TransformBeforePolarToGeographic() */ /* */ /* Transform the geometry (by intersection), so as to cut each geometry */ /* that crosses the pole, in 2 parts. Do also tricks for geometries */ /* that just touch the pole. */ /************************************************************************/ -static std::unique_ptr TransformBeforePolarToWGS84( +static std::unique_ptr TransformBeforePolarToGeographic( OGRCoordinateTransformation *poRevCT, bool bIsNorthPolar, std::unique_ptr poDstGeom, bool &bNeedPostCorrectionOut) { @@ -3449,15 +3449,15 @@ static std::unique_ptr TransformBeforePolarToWGS84( } /************************************************************************/ -/* IsAntimeridianProjToWGS84() */ +/* IsAntimeridianProjToGeographic() */ /* */ /* Returns true if poCT transforms from a projection that includes the */ /* antimeridian in a continuous way. */ /************************************************************************/ -static bool IsAntimeridianProjToWGS84(OGRCoordinateTransformation *poCT, - OGRCoordinateTransformation *poRevCT, - OGRGeometry *poDstGeometry) +static bool IsAntimeridianProjToGeographic(OGRCoordinateTransformation *poCT, + OGRCoordinateTransformation *poRevCT, + OGRGeometry *poDstGeometry) { const bool bBackupEmitErrors = poCT->GetEmitErrors(); poRevCT->SetEmitErrors(false); @@ -3633,13 +3633,13 @@ struct SortPointsByAscendingY }; /************************************************************************/ -/* TransformBeforeAntimeridianToWGS84() */ +/* TransformBeforeAntimeridianToGeographic() */ /* */ /* Transform the geometry (by intersection), so as to cut each geometry */ /* that crosses the antimeridian, in 2 parts. */ /************************************************************************/ -static std::unique_ptr TransformBeforeAntimeridianToWGS84( +static std::unique_ptr TransformBeforeAntimeridianToGeographic( OGRCoordinateTransformation *poCT, OGRCoordinateTransformation *poRevCT, std::unique_ptr poDstGeom, bool &bNeedPostCorrectionOut) { @@ -3820,9 +3820,22 @@ static void SnapCoordsCloseToLatLongBounds(OGRGeometry *poGeom) struct OGRGeometryFactory::TransformWithOptionsCache::Private { + const OGRSpatialReference *poSourceCRS = nullptr; + const OGRSpatialReference *poTargetCRS = nullptr; + const OGRCoordinateTransformation *poCT = nullptr; std::unique_ptr poRevCT{}; bool bIsPolar = false; bool bIsNorthPolar = false; + + void clear() + { + poSourceCRS = nullptr; + poTargetCRS = nullptr; + poCT = nullptr; + poRevCT.reset(); + bIsPolar = false; + bIsNorthPolar = false; + } }; /************************************************************************/ @@ -3858,52 +3871,76 @@ OGRGeometry *OGRGeometryFactory::transformWithOptions( char **papszOptions, CPL_UNUSED const TransformWithOptionsCache &cache) { auto poDstGeom = std::unique_ptr(poSrcGeom->clone()); - if (poCT != nullptr) + if (poCT) { #ifdef HAVE_GEOS bool bNeedPostCorrection = false; - - auto poSourceCRS = poCT->GetSourceCS(); - auto poTargetCRS = poCT->GetTargetCS(); - if (poSourceCRS != nullptr && poTargetCRS != nullptr && - poSourceCRS->IsProjected() && poTargetCRS->IsGeographic()) + const auto poSourceCRS = poCT->GetSourceCS(); + const auto poTargetCRS = poCT->GetTargetCS(); + const auto eSrcGeomType = wkbFlatten(poSrcGeom->getGeometryType()); + // Check if we are transforming from projected coordinates to + // geographic coordinates, with a chance that there might be polar or + // anti-meridian discontinuities. If so, create the inverse transform. + if (eSrcGeomType != wkbPoint && eSrcGeomType != wkbMultiPoint && + (poSourceCRS != cache.d->poSourceCRS || + poTargetCRS != cache.d->poTargetCRS || poCT != cache.d->poCT)) { - OGRSpatialReference oSRSWGS84; - oSRSWGS84.SetWellKnownGeogCS("WGS84"); - oSRSWGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); - if (poTargetCRS->IsSame(&oSRSWGS84)) + cache.d->clear(); + cache.d->poSourceCRS = poSourceCRS; + cache.d->poTargetCRS = poTargetCRS; + cache.d->poCT = poCT; + if (poSourceCRS && poTargetCRS && poSourceCRS->IsProjected() && + poTargetCRS->IsGeographic() && + poTargetCRS->GetAxisMappingStrategy() == + OAMS_TRADITIONAL_GIS_ORDER && + // check that angular units is degree + std::fabs(poTargetCRS->GetAngularUnits(nullptr) - + CPLAtof(SRS_UA_DEGREE_CONV)) <= + 1e-8 * CPLAtof(SRS_UA_DEGREE_CONV)) { - if (cache.d->poRevCT == nullptr || - !cache.d->poRevCT->GetTargetCS()->IsSame(poSourceCRS)) + double dfWestLong = 0.0; + double dfSouthLat = 0.0; + double dfEastLong = 0.0; + double dfNorthLat = 0.0; + if (poSourceCRS->GetAreaOfUse(&dfWestLong, &dfSouthLat, + &dfEastLong, &dfNorthLat, + nullptr) && + !(dfSouthLat == -90.0 || dfNorthLat == 90.0 || + dfWestLong == -180.0 || dfEastLong == 180.0 || + dfWestLong > dfEastLong)) + { + // Not a global geographic CRS + } + else { cache.d->poRevCT.reset(OGRCreateCoordinateTransformation( - &oSRSWGS84, poSourceCRS)); + poTargetCRS, poSourceCRS)); cache.d->bIsNorthPolar = false; cache.d->bIsPolar = false; + cache.d->poRevCT.reset(poCT->GetInverse()); if (cache.d->poRevCT && - IsPolarToWGS84(poCT, cache.d->poRevCT.get(), - cache.d->bIsNorthPolar)) + IsPolarToGeographic(poCT, cache.d->poRevCT.get(), + cache.d->bIsNorthPolar)) { cache.d->bIsPolar = true; } } - auto poRevCT = cache.d->poRevCT.get(); - if (poRevCT != nullptr) - { - if (cache.d->bIsPolar) - { - poDstGeom = TransformBeforePolarToWGS84( - poRevCT, cache.d->bIsNorthPolar, - std::move(poDstGeom), bNeedPostCorrection); - } - else if (IsAntimeridianProjToWGS84(poCT, poRevCT, - poDstGeom.get())) - { - poDstGeom = TransformBeforeAntimeridianToWGS84( - poCT, poRevCT, std::move(poDstGeom), - bNeedPostCorrection); - } - } + } + } + + if (auto poRevCT = cache.d->poRevCT.get()) + { + if (cache.d->bIsPolar) + { + poDstGeom = TransformBeforePolarToGeographic( + poRevCT, cache.d->bIsNorthPolar, std::move(poDstGeom), + bNeedPostCorrection); + } + else if (IsAntimeridianProjToGeographic(poCT, poRevCT, + poDstGeom.get())) + { + poDstGeom = TransformBeforeAntimeridianToGeographic( + poCT, poRevCT, std::move(poDstGeom), bNeedPostCorrection); } } #endif