Skip to content

Commit

Permalink
createOperations(): fix NAD83_CSRS_1997_MTMxx_HT2_1997 to NAD83_CSRS_…
Browse files Browse the repository at this point in the history 
…2010_UTMyy_CGVD2013_2010
  • Loading branch information
@rouault
rouault committed Sep 5, 2023
1 parent 601a599 commit 223a112
Show file tree
Hide file tree
Showing 2 changed files with 118 additions and 4 deletions.
84 changes: 80 additions & 4 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<const crs::DerivedCRS *>(componentsSrc[0].get());
Expand Down Expand Up @@ -5607,6 +5611,82 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
return;
}

auto geogSrc = dynamic_cast<const crs::GeographicCRS *>(
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<crs::CRSNNPtr> 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<const crs::BoundCRS *>(componentsSrc[0].get());
if (boundSrc) {
Expand Down Expand Up @@ -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
Expand Down
38 changes: 38 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down

0 comments on commit 223a112

Please sign in to comment.