Skip to content

Commit

Permalink
OGRGeometryFactory::transformWithOptions(): deal with polar or anti-m…
Browse files Browse the repository at this point in the history 
…eridian discontinuities when going from projected to (any) geographic CRS

The logic already existed but was restricted to WGS 84 without any good
reason.
  • Loading branch information
@rouault
rouault committed Sep 2, 2024
1 parent a0ae60f commit 1e1d2fe
Show file tree
Hide file tree
Showing 2 changed files with 148 additions and 45 deletions.
67 changes: 67 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,71 @@ 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;
}

// 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;
}
}

} // namespace
126 changes: 81 additions & 45 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,17 +3366,17 @@ 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 OGRGeometry *
TransformBeforePolarToWGS84(OGRCoordinateTransformation *poRevCT,
bool bIsNorthPolar, OGRGeometry *poDstGeom,
bool &bNeedPostCorrectionOut)
TransformBeforePolarToGeographic(OGRCoordinateTransformation *poRevCT,
bool bIsNorthPolar, OGRGeometry *poDstGeom,
bool &bNeedPostCorrectionOut)
{
const int nSign = (bIsNorthPolar) ? 1 : -1;

Expand Down Expand Up @@ -3450,15 +3450,15 @@ TransformBeforePolarToWGS84(OGRCoordinateTransformation *poRevCT,
}

/************************************************************************/
/* 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 @@ -3634,13 +3634,13 @@ struct SortPointsByAscendingY
};

/************************************************************************/
/* TransformBeforeAntimeridianToWGS84() */
/* TransformBeforeAntimeridianToGeographic() */
/* */
/* Transform the geometry (by intersection), so as to cut each geometry */
/* that crosses the antimeridian, in 2 parts. */
/************************************************************************/

static OGRGeometry *TransformBeforeAntimeridianToWGS84(
static OGRGeometry *TransformBeforeAntimeridianToGeographic(
OGRCoordinateTransformation *poCT, OGRCoordinateTransformation *poRevCT,
OGRGeometry *poDstGeom, bool &bNeedPostCorrectionOut)
{
Expand Down Expand Up @@ -3821,9 +3821,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 @@ -3859,51 +3872,74 @@ OGRGeometry *OGRGeometryFactory::transformWithOptions(
char **papszOptions, CPL_UNUSED const TransformWithOptionsCache &cache)
{
OGRGeometry *poDstGeom = 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 (poTargetCRS->GetAreaOfUse(&dfWestLong, &dfSouthLat,
&dfEastLong, &dfNorthLat,
nullptr) &&
!(dfSouthLat == -90.0 || dfNorthLat == 90.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, poDstGeom,
bNeedPostCorrection);
}
else if (IsAntimeridianProjToWGS84(poCT, poRevCT,
poDstGeom))
{
poDstGeom = TransformBeforeAntimeridianToWGS84(
poCT, poRevCT, poDstGeom, bNeedPostCorrection);
}
}
}
}

if (auto poRevCT = cache.d->poRevCT.get())
{
if (cache.d->bIsPolar)
{
poDstGeom = TransformBeforePolarToGeographic(
poRevCT, cache.d->bIsNorthPolar, poDstGeom,
bNeedPostCorrection);
}
else if (IsAntimeridianProjToGeographic(poCT, poRevCT, poDstGeom))
{
poDstGeom = TransformBeforeAntimeridianToGeographic(
poCT, poRevCT, poDstGeom, bNeedPostCorrection);
}
}
#endif
Expand Down

0 comments on commit 1e1d2fe

Please sign in to comment.