Skip to content

Commit

Permalink
createOperations(): tune for 'NAD83_CSRS_1997_xxxx_HT2_1997 to NAD83_…
Browse files Browse the repository at this point in the history 
…CSRS_1997_yyyy_CGVD2013_1997' or 'NAD83_CSRS_2002_xxxx_HT2_2002 to NAD83_CSRS_2002_yyyy_CGVD2013_2002'
  • Loading branch information
@rouault
rouault committed Sep 5, 2023
1 parent 223a112 commit be5501d
Show file tree
Hide file tree
Showing 2 changed files with 219 additions and 4 deletions.
79 changes: 75 additions & 4 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6298,8 +6298,40 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
const auto vertDst = componentsDst[1]->extractVerticalCRS();
if (vertSrc && vertDst &&
!componentsSrc[1]->_isEquivalentTo(componentsDst[1].get())) {
if (!vertSrc->geoidModel().empty() ||
!vertDst->geoidModel().empty()) {
if ((!vertSrc->geoidModel().empty() ||
!vertDst->geoidModel().empty()) &&
// To be able to use "CGVD28 height to
// CGVD2013a(1997|2002|2010)" single grid
!(vertSrc->nameStr() == "CGVD28 height" &&
!vertSrc->geoidModel().empty() &&
vertSrc->geoidModel().front()->nameStr() == "HT2_1997" &&
vertDst->nameStr() == "CGVD2013a(1997) height" &&
vertDst->geoidModel().empty()) &&
!(vertSrc->nameStr() == "CGVD28 height" &&
!vertSrc->geoidModel().empty() &&
vertSrc->geoidModel().front()->nameStr() == "HT2_2002" &&
vertDst->nameStr() == "CGVD2013a(2002) height" &&
vertDst->geoidModel().empty()) &&
!(vertSrc->nameStr() == "CGVD28 height" &&
!vertSrc->geoidModel().empty() &&
vertSrc->geoidModel().front()->nameStr() == "HT2_2010" &&
vertDst->nameStr() == "CGVD2013a(2010) height" &&
vertDst->geoidModel().empty()) &&
!(vertDst->nameStr() == "CGVD28 height" &&
!vertDst->geoidModel().empty() &&
vertDst->geoidModel().front()->nameStr() == "HT2_1997" &&
vertSrc->nameStr() == "CGVD2013a(1997) height" &&
vertSrc->geoidModel().empty()) &&
!(vertDst->nameStr() == "CGVD28 height" &&
!vertDst->geoidModel().empty() &&
vertDst->geoidModel().front()->nameStr() == "HT2_2002" &&
vertSrc->nameStr() == "CGVD2013a(2002) height" &&
vertSrc->geoidModel().empty()) &&
!(vertDst->nameStr() == "CGVD28 height" &&
!vertDst->geoidModel().empty() &&
vertDst->geoidModel().front()->nameStr() == "HT2_2010" &&
vertSrc->nameStr() == "CGVD2013a(2010) height" &&
vertSrc->geoidModel().empty())) {
// If we have a geoid model, force using through it
bTryThroughIntermediateGeogCRS = true;
} else {
Expand Down Expand Up @@ -6607,10 +6639,49 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
}
}

auto opSrcCRSToGeogCRS = createOperations(
// Hack for
// NAD83_CSRS_1997_xxxx_HT2_1997 to NAD83_CSRS_1997_yyyy_CGVD2013_1997
// NAD83_CSRS_2002_xxxx_HT2_2002 to NAD83_CSRS_2002_yyyy_CGVD2013_2002
if (sourceEpoch.has_value() && targetEpoch.has_value() &&
sourceEpoch->coordinateEpoch()._isEquivalentTo(
targetEpoch->coordinateEpoch()) &&
srcGeog->_isEquivalentTo(
dstGeog.get(), util::IComparable::Criterion::EQUIVALENT) &&
srcGeog->nameStr() == "NAD83(CSRS)v7") {
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");
auto nad83CSRSv7 = authFactoryEPSG->createGeographicCRS("8255");
if (srcGeog->_isEquivalentTo(nad83CSRSv7.get(),
util::IComparable::Criterion::
EQUIVALENT) && // NAD83(CSRS)v7
// 2D
((is1997 &&
verticalTransform->nameStr().find(
"CGVD28 height to CGVD2013a(1997) height (1)") !=
std::string::npos) ||
(is2002 &&
verticalTransform->nameStr().find(
"CGVD28 height to CGVD2013a(2002) height (1)") !=
std::string::npos))) {
interpolationGeogCRS = nad83CSRSv7;
}
} catch (const std::exception &) {
}
}

const auto opSrcCRSToGeogCRS = createOperations(
componentsSrc[0], util::optional<common::DataEpoch>(),
interpolationGeogCRS, util::optional<common::DataEpoch>(), context);
auto opGeogCRStoDstCRS = createOperations(
const auto opGeogCRStoDstCRS = createOperations(
interpolationGeogCRS, util::optional<common::DataEpoch>(),
componentsDst[0], util::optional<common::DataEpoch>(), context);
for (const auto &opSrc : opSrcCRSToGeogCRS) {
Expand Down
144 changes: 144 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9608,6 +9608,150 @@ TEST(operation,

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

TEST(operation,
createOperation_compound_to_compound_with_epoch_HT2_1997_CGG2013a) {
auto dbContext = DatabaseContext::create();
auto factoryNRCAN = AuthorityFactory::create(dbContext, "NRCAN");
auto sourceCM =
factoryNRCAN->createCoordinateMetadata("NAD83_CSRS_1997_MTM7_HT2_1997");
auto targetCM = factoryNRCAN->createCoordinateMetadata(
"NAD83_CSRS_1997_UTM19_CGVD2013_1997");
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=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 "
"+x_0=304800 +y_0=0 +ellps=GRS80 "
"+step +proj=vgridshift +grids=ca_nrc_HT2_1997_CGG2013a.tif "
"+multiplier=1 "
"+step +proj=utm +zone=19 +ellps=GRS80");
}
{
auto list = CoordinateOperationFactory::create()->createOperations(
targetCM, sourceCM, 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=19 +ellps=GRS80 "
"+step +inv +proj=vgridshift "
"+grids=ca_nrc_HT2_1997_CGG2013a.tif "
"+multiplier=1 "
"+step +proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 "
"+x_0=304800 +y_0=0 +ellps=GRS80");
}
}

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

TEST(operation,
createOperation_compound_to_compound_with_epoch_HT2_2002_CGG2013a) {
auto dbContext = DatabaseContext::create();
auto factoryNRCAN = AuthorityFactory::create(dbContext, "NRCAN");
auto sourceCM =
factoryNRCAN->createCoordinateMetadata("NAD83_CSRS_2002_MTM7_HT2_2002");
auto targetCM = factoryNRCAN->createCoordinateMetadata(
"NAD83_CSRS_2002_UTM19_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=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 "
"+x_0=304800 +y_0=0 +ellps=GRS80 "
"+step +proj=vgridshift +grids=ca_nrc_HT2_2002v70_CGG2013a.tif "
"+multiplier=1 "
"+step +proj=utm +zone=19 +ellps=GRS80");
}
{
auto list = CoordinateOperationFactory::create()->createOperations(
targetCM, sourceCM, 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=19 +ellps=GRS80 "
"+step +inv +proj=vgridshift "
"+grids=ca_nrc_HT2_2002v70_CGG2013a.tif "
"+multiplier=1 "
"+step +proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 "
"+x_0=304800 +y_0=0 +ellps=GRS80");
}
}

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

TEST(operation,
createOperation_compound_to_compound_with_epoch_HT2_2010_CGG2013a) {
auto dbContext = DatabaseContext::create();
auto factoryNRCAN = AuthorityFactory::create(dbContext, "NRCAN");
auto sourceCM =
factoryNRCAN->createCoordinateMetadata("NAD83_CSRS_2010_MTM7_HT2_2010");
auto targetCM = factoryNRCAN->createCoordinateMetadata(
"NAD83_CSRS_2010_UTM19_CGVD2013_2010");
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=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 "
"+x_0=304800 +y_0=0 +ellps=GRS80 "
"+step +proj=vgridshift +grids=ca_nrc_HT2_2010v70_CGG2013a.tif "
"+multiplier=1 "
"+step +proj=utm +zone=19 +ellps=GRS80");
}
{
auto list = CoordinateOperationFactory::create()->createOperations(
targetCM, sourceCM, 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=19 +ellps=GRS80 "
"+step +inv +proj=vgridshift "
"+grids=ca_nrc_HT2_2010v70_CGG2013a.tif "
"+multiplier=1 "
"+step +proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 "
"+x_0=304800 +y_0=0 +ellps=GRS80");
}
}

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

TEST(
operation,
createOperation_compound_to_compound_with_Geographic3D_Offset_by_velocity_grid) {
Expand Down

0 comments on commit be5501d

Please sign in to comment.