From 223a1122127bdab77cf0c85ef1d19cb7afbc3370 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 5 Sep 2023 19:25:45 +0200 Subject: [PATCH] createOperations(): fix NAD83_CSRS_1997_MTMxx_HT2_1997 to NAD83_CSRS_2010_UTMyy_CGVD2013_2010 --- .../operation/coordinateoperationfactory.cpp | 84 ++++++++++++++++++- test/unit/test_operationfactory.cpp | 38 +++++++++ 2 files changed, 118 insertions(+), 4 deletions(-) diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 7d0147cf9b..a4237d87d9 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -5576,6 +5576,10 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( const auto &componentsSrc = compoundSrc->componentReferenceSystems(); if (!componentsSrc.empty()) { + const auto dbContext = + authFactory ? authFactory->databaseContext().as_nullable() + : nullptr; + if (componentsSrc.size() == 2) { auto derivedHSrc = dynamic_cast(componentsSrc[0].get()); @@ -5607,6 +5611,82 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( return; } + auto geogSrc = dynamic_cast( + componentsSrc[0].get()); + // Sorry for this hack... aimed at transforming + // "NAD83(CSRS)v7 + CGVD2013a(1997) height @ 1997" to "NAD83(CSRS)v7 + // @ 1997" to "NAD83(CSRS)v3 + CGVD2013a(1997) height" to + // "NAD83(CSRS)v3" OR "NAD83(CSRS)v7 + CGVD2013a(2002) height @ + // 2002" to "NAD83(CSRS)v7 @ 2002" to "NAD83(CSRS)v4 + + // CGVD2013a(2002) height" to "NAD83(CSRS)v4" + if (dbContext && geogSrc && geogSrc->nameStr() == "NAD83(CSRS)v7" && + sourceEpoch.has_value() && + geogDst->coordinateSystem()->axisList().size() == 3U && + geogDst->nameStr() == geogSrc->nameStr() && + targetEpoch.has_value() && + sourceEpoch->coordinateEpoch()._isEquivalentTo( + targetEpoch->coordinateEpoch())) { + const bool is1997 = + std::abs(sourceEpoch->coordinateEpoch().convertToUnit( + common::UnitOfMeasure::YEAR) - + 1997) < 1e-10; + const bool is2002 = + std::abs(sourceEpoch->coordinateEpoch().convertToUnit( + common::UnitOfMeasure::YEAR) - + 2002) < 1e-10; + try { + auto authFactoryEPSG = io::AuthorityFactory::create( + authFactory->databaseContext(), "EPSG"); + if (geogSrc->_isEquivalentTo( + authFactoryEPSG + ->createCoordinateReferenceSystem("8255") + .get(), + util::IComparable::Criterion:: + EQUIVALENT) && // NAD83(CSRS)v7 + // 2D + geogDst->_isEquivalentTo( + authFactoryEPSG + ->createCoordinateReferenceSystem("8254") + .get(), + util::IComparable::Criterion:: + EQUIVALENT) && // NAD83(CSRS)v7 + // 3D + ((is1997 && componentsSrc[1]->nameStr() == + "CGVD2013a(1997) height") || + (is2002 && componentsSrc[1]->nameStr() == + "CGVD2013a(2002) height"))) { + const auto newGeogCRS_2D( + authFactoryEPSG->createCoordinateReferenceSystem( + is1997 ? "8240" : // NAD83(CSRS)v3 2D + "8246" // NAD83(CSRS)v4 2D + )); + const auto newGeogCRS_3D( + authFactoryEPSG->createCoordinateReferenceSystem( + is1997 ? "8239" : // NAD83(CSRS)v3 3D + "8244" // NAD83(CSRS)v4 3D + )); + std::vector intermComponents{ + newGeogCRS_2D, componentsSrc[1]}; + auto properties = util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + intermComponents[0]->nameStr() + " + " + + intermComponents[1]->nameStr()); + auto newCompound = crs::CompoundCRS::create( + properties, intermComponents); + auto ops = createOperations(newCompound, sourceEpoch, + newGeogCRS_3D, sourceEpoch, + context); + for (const auto &op : ops) { + auto opClone = op->shallowClone(); + setCRSs(opClone.get(), sourceCRS, targetCRS); + res.emplace_back(opClone); + } + return; + } + } catch (const std::exception &) { + } + } + auto boundSrc = dynamic_cast(componentsSrc[0].get()); if (boundSrc) { @@ -5645,10 +5725,6 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( } } - const auto dbContext = - authFactory ? authFactory->databaseContext().as_nullable() - : nullptr; - // Deal with "+proj=something +geoidgrids +nadgrids/+towgs84" to // another CRS whose datum is not the same as the horizontal datum // of the source diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index ba5f5ff08f..818c3a2816 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -9696,6 +9696,44 @@ TEST( // --------------------------------------------------------------------------- +TEST( + operation, + createOperation_compound_to_compound_with_point_motion_operation_special_case_CGVD2013a) { + auto dbContext = DatabaseContext::create(); + auto factoryNRCAN = AuthorityFactory::create(dbContext, "NRCAN"); + auto sourceCM = factoryNRCAN->createCoordinateMetadata( + "NAD83_CSRS_1997_UTM17_CGVD2013_1997"); + auto targetCM = factoryNRCAN->createCoordinateMetadata( + "NAD83_CSRS_2002_UTM17_CGVD2013_2002"); + auto ctxt = CoordinateOperationContext::create( + AuthorityFactory::create(dbContext, std::string()), nullptr, 0); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + auto list = CoordinateOperationFactory::create()->createOperations( + sourceCM, targetCM, ctxt); + ASSERT_GE(list.size(), 1U); + EXPECT_FALSE(list[0]->hasBallparkTransformation()); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +inv +proj=utm +zone=17 +ellps=GRS80 " + "+step +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif " + "+multiplier=1 " + "+step +proj=cart +ellps=GRS80 " + "+step +proj=set +v_4=1997 +omit_fwd " + "+step +proj=deformation +dt=5 +grids=ca_nrc_NAD83v70VG.tif " + "+ellps=GRS80 " + "+step +proj=set +v_4=2002 +omit_inv " + "+step +inv +proj=cart +ellps=GRS80 " + "+step +inv +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif " + "+multiplier=1 " + "+step +proj=utm +zone=17 +ellps=GRS80"); +} + +// --------------------------------------------------------------------------- + TEST(operation, createOperation_Geographic3D_Offset_by_velocity_grid) { auto dbContext = DatabaseContext::create(); auto factoryEPSG = AuthorityFactory::create(dbContext, "EPSG");