Skip to content

Commit

Permalink
createOperations(): fix bad PROJ pipeline when converting between Geo…
Browse files Browse the repository at this point in the history 
…g3D with non-metre height to CompoundCRS (fixes #3938)
  • Loading branch information
@rouault
rouault committed Nov 2, 2023
1 parent f6765ed commit 0a9417f
Show file tree
Hide file tree
Showing 2 changed files with 207 additions and 4 deletions.
4 changes: 0 additions & 4 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2259,16 +2259,12 @@ struct MyPROJStringExportableHorizVertical final
// cppcheck-suppress functionStatic
_exportToPROJString(io::PROJStringFormatter *formatter) const override {

formatter->pushOmitZUnitConversion();

horizTransform->_exportToPROJString(formatter);

formatter->startInversion();
geogDst->addAngularUnitConvertAndAxisSwap(formatter);
formatter->stopInversion();

formatter->popOmitZUnitConversion();

formatter->pushOmitHorizontalConversionInVertTransformation();
verticalTransform->_exportToPROJString(formatter);
formatter->popOmitHorizontalConversionInVertTransformation();
Expand Down
207 changes: 207 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7134,6 +7134,213 @@ TEST(operation, compoundCRS_to_geogCRS_with_vertical_unit_change) {

// ---------------------------------------------------------------------------

// Use case of https://github.com/OSGeo/PROJ/issues/3938
TEST(operation, compoundCRS_ftUS_to_geogCRS_ft) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
ctxt->setGridAvailabilityUse(
CoordinateOperationContext::GridAvailabilityUse::
IGNORE_GRID_AVAILABILITY);
// NAD83(2011) + NAVD88 height (ftUS)
auto srcObj = createFromUserInput("EPSG:6318+6360",
authFactory->databaseContext(), false);
auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
ASSERT_TRUE(src != nullptr);
auto nnSrc = NN_NO_CHECK(src);
auto dst =
authFactory->createCoordinateReferenceSystem("6319")->alterCSLinearUnit(
UnitOfMeasure::FOOT); // NAD83(2011) with foot

auto res = CoordinateOperationFactory::create()->createOperations(
nnSrc, dst, ctxt);
ASSERT_TRUE(!res.empty());

EXPECT_EQ(
res[0]->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad +z_out=m "
"+step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=ft "
"+step +proj=axisswap +order=2,1");

EXPECT_EQ(
res.back()->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=deg +z_out=ft");

auto resInv = CoordinateOperationFactory::create()->createOperations(
dst, nnSrc, ctxt);
ASSERT_TRUE(!resInv.empty());

EXPECT_EQ(
resInv[0]->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +z_in=ft +xy_out=rad +z_out=m "
"+step +inv +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=us-ft "
"+step +proj=axisswap +order=2,1");

EXPECT_EQ(
resInv.back()->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=unitconvert +xy_in=deg +z_in=ft +xy_out=deg +z_out=us-ft");
}

// ---------------------------------------------------------------------------

// Use case of https://github.com/OSGeo/PROJ/issues/3938
TEST(operation, compoundCRS_ft_to_geogCRS_ft) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
ctxt->setGridAvailabilityUse(
CoordinateOperationContext::GridAvailabilityUse::
IGNORE_GRID_AVAILABILITY);
// NAD83(2011) + NAVD88 height (ft)
auto srcObj = createFromUserInput("EPSG:6318+8228",
authFactory->databaseContext(), false);
auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
ASSERT_TRUE(src != nullptr);
auto nnSrc = NN_NO_CHECK(src);
auto dst =
authFactory->createCoordinateReferenceSystem("6319")->alterCSLinearUnit(
UnitOfMeasure::FOOT); // NAD83(2011) with foot

auto res = CoordinateOperationFactory::create()->createOperations(
nnSrc, dst, ctxt);
ASSERT_TRUE(!res.empty());

EXPECT_EQ(
res[0]->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +z_in=ft +xy_out=rad +z_out=m "
"+step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=ft "
"+step +proj=axisswap +order=2,1");

EXPECT_EQ(
res.back()->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=noop");

auto resInv = CoordinateOperationFactory::create()->createOperations(
dst, nnSrc, ctxt);
ASSERT_TRUE(!resInv.empty());

EXPECT_EQ(
resInv[0]->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +z_in=ft +xy_out=rad +z_out=m "
"+step +inv +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=ft "
"+step +proj=axisswap +order=2,1");

EXPECT_EQ(
resInv.back()->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=noop");
}

// ---------------------------------------------------------------------------

// Use case of https://github.com/OSGeo/PROJ/issues/3938
TEST(operation, compoundCRS_m_to_geogCRS_ft) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
ctxt->setGridAvailabilityUse(
CoordinateOperationContext::GridAvailabilityUse::
IGNORE_GRID_AVAILABILITY);
// NAD83(2011) + NAVD88 height
auto srcObj = createFromUserInput("EPSG:6318+5703",
authFactory->databaseContext(), false);
auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
ASSERT_TRUE(src != nullptr);
auto nnSrc = NN_NO_CHECK(src);
auto dst =
authFactory->createCoordinateReferenceSystem("6319")->alterCSLinearUnit(
UnitOfMeasure::FOOT); // NAD83(2011) with foot

auto res = CoordinateOperationFactory::create()->createOperations(
nnSrc, dst, ctxt);
ASSERT_TRUE(!res.empty());

EXPECT_EQ(
res[0]->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=ft "
"+step +proj=axisswap +order=2,1");

EXPECT_EQ(
res.back()->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=unitconvert +z_in=m +z_out=ft");

auto resInv = CoordinateOperationFactory::create()->createOperations(
dst, nnSrc, ctxt);
ASSERT_TRUE(!resInv.empty());

EXPECT_EQ(
resInv[0]->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +z_in=ft +xy_out=rad +z_out=m "
"+step +inv +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");

EXPECT_EQ(
resInv.back()->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=unitconvert +z_in=ft +z_out=m");
}

// ---------------------------------------------------------------------------

TEST(
operation,
compoundCRS_to_geogCRS_with_vertical_unit_change_and_complex_horizontal_change) {
Expand Down

0 comments on commit 0a9417f

Please sign in to comment.