Skip to content

Commit

Permalink
Merge pull request OSGeo#10716 from rouault/transformWithOptions
Browse files Browse the repository at this point in the history 
OGRGeometryFactory::transformWithOptions(): deal with polar or anti-meridian discontinuities when going from projected to (any) geographic CRS
  • Loading branch information
@rouault
rouault authored Sep 16, 2024
2 parents 6ccf77f + ffb6004 commit 2742d2b
Show file tree
Hide file tree
Showing 3 changed files with 181 additions and 47 deletions.
34 changes: 31 additions & 3 deletions apps/ogr2ogr_lib.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@
#include <limits>
#include <map>
#include <memory>
#include <mutex>
#include <set>
#include <unordered_set>
#include <string>
Expand Down Expand Up @@ -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;
}
};
Expand Down Expand Up @@ -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<OGRCoordinateTransformation>(poCT1->GetInverse());
auto poInvCT2 =
std::unique_ptr<OGRCoordinateTransformation>(poCT2->GetInverse());
if (!poInvCT1 || !poInvCT2)
return nullptr;
return std::make_unique<CompositeCT>(poInvCT2.release(), true,
poInvCT1.release(), true)
.release();
}
};

Expand All @@ -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<int> &mappingIn,
const std::vector<int> &mappingOut)
{
Expand Down Expand Up @@ -1360,7 +1388,7 @@ class AxisMappingCoordinateTransformation : public OGRCoordinateTransformation

virtual OGRCoordinateTransformation *GetInverse() const override
{
return nullptr;
return new AxisMappingCoordinateTransformation(bSwapXY);
}
};

Expand Down
69 changes: 69 additions & 0 deletions autotest/cpp/test_ogr.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<OGRCoordinateTransformation>(
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<OGRCoordinateTransformation>(
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
125 changes: 81 additions & 44 deletions ogr/ogrgeometryfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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<OGRGeometry> TransformBeforePolarToWGS84(
static std::unique_ptr<OGRGeometry> TransformBeforePolarToGeographic(
OGRCoordinateTransformation *poRevCT, bool bIsNorthPolar,
std::unique_ptr<OGRGeometry> poDstGeom, bool &bNeedPostCorrectionOut)
{
Expand Down Expand Up @@ -3449,15 +3449,15 @@ static std::unique_ptr<OGRGeometry> 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);
Expand Down Expand Up @@ -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<OGRGeometry> TransformBeforeAntimeridianToWGS84(
static std::unique_ptr<OGRGeometry> TransformBeforeAntimeridianToGeographic(
OGRCoordinateTransformation *poCT, OGRCoordinateTransformation *poRevCT,
std::unique_ptr<OGRGeometry> poDstGeom, bool &bNeedPostCorrectionOut)
{
Expand Down Expand Up @@ -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<OGRCoordinateTransformation> poRevCT{};
bool bIsPolar = false;
bool bIsNorthPolar = false;

void clear()
{
poSourceCRS = nullptr;
poTargetCRS = nullptr;
poCT = nullptr;
poRevCT.reset();
bIsPolar = false;
bIsNorthPolar = false;
}
};

/************************************************************************/
Expand Down Expand Up @@ -3858,52 +3871,76 @@ OGRGeometry *OGRGeometryFactory::transformWithOptions(
char **papszOptions, CPL_UNUSED const TransformWithOptionsCache &cache)
{
auto poDstGeom = std::unique_ptr<OGRGeometry>(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
Expand Down

0 comments on commit 2742d2b

Please sign in to comment.