diff options
| -rw-r--r-- | include/proj/io.hpp | 4 | ||||
| -rw-r--r-- | src/iso19111/coordinateoperation.cpp | 479 | ||||
| -rw-r--r-- | src/iso19111/factory.cpp | 24 | ||||
| -rw-r--r-- | test/unit/test_operation.cpp | 121 |
4 files changed, 465 insertions, 163 deletions
diff --git a/include/proj/io.hpp b/include/proj/io.hpp index 12b3b111..37b94c1e 100644 --- a/include/proj/io.hpp +++ b/include/proj/io.hpp @@ -1089,6 +1089,10 @@ class PROJ_GCC_DLL AuthorityFactory { const std::string &datum_code, const std::string &geodetic_crs_type) const; + PROJ_INTERNAL std::list<crs::VerticalCRSNNPtr> + createVerticalCRSFromDatum(const std::string &datum_auth_name, + const std::string &datum_code) const; + PROJ_INTERNAL std::list<crs::GeodeticCRSNNPtr> createGeodeticCRSFromEllipsoid(const std::string &ellipsoid_auth_name, const std::string &ellipsoid_code, diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index c7581642..27ac712c 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -10298,6 +10298,7 @@ struct CoordinateOperationFactory::Private { bool inCreateOperationsWithDatumPivotAntiRecursion = false; bool inCreateOperationsThroughPreferredHub = false; bool inCreateOperationsGeogToVertWithAlternativeGeog = false; + bool inCreateOperationsGeogToVertWithIntermediateVert = false; bool skipHorizontalTransformation = false; Context(const crs::CRSNNPtr &sourceCRSIn, @@ -10327,6 +10328,23 @@ struct CoordinateOperationFactory::Private { const crs::VerticalCRS *vertDst, std::vector<CoordinateOperationNNPtr> &res); + static std::vector<CoordinateOperationNNPtr> + createOperationsGeogToVertWithIntermediateVert( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + const crs::VerticalCRS *vertDst, Context &context); + + static std::vector<CoordinateOperationNNPtr> + createOperationsGeogToVertWithAlternativeGeog( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Context &context); + + static void createOperationsFromDatabaseWithVertCRS( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::GeographicCRS *geogSrc, + const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc, + const crs::VerticalCRS *vertDst, + std::vector<CoordinateOperationNNPtr> &res); + static void createOperationsGeodToGeod( const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, Private::Context &context, const crs::GeodeticCRS *geodSrc, @@ -10409,11 +10427,6 @@ struct CoordinateOperationFactory::Private { const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst, Context &context); - static std::vector<CoordinateOperationNNPtr> - createOperationsGeogToVertWithAlternativeGeog( - const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, - Context &context); - static bool hasPerfectAccuracyResult(const std::vector<CoordinateOperationNNPtr> &res, const Context &context); @@ -12402,51 +12415,6 @@ void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub( // --------------------------------------------------------------------------- -std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private:: - createOperationsGeogToVertWithAlternativeGeog( - const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS - const crs::CRSNNPtr &targetCRS, // vertical CRS - Private::Context &context) { - - std::vector<CoordinateOperationNNPtr> res; - - struct AntiRecursionGuard { - Context &context; - - explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) { - assert(!context.inCreateOperationsGeogToVertWithAlternativeGeog); - context.inCreateOperationsGeogToVertWithAlternativeGeog = true; - } - - ~AntiRecursionGuard() { - context.inCreateOperationsGeogToVertWithAlternativeGeog = false; - } - }; - AntiRecursionGuard guard(context); - - for (int i = 0; i < 2; i++) { - - // Generally EPSG has operations from GeogCrs to VertCRS - auto ops = - i == 0 ? findOpsInRegistryDirectTo(targetCRS, context.context) - : findOpsInRegistryDirectFrom(targetCRS, context.context); - - for (const auto &op : ops) { - const auto tmpCRS = i == 0 ? op->sourceCRS() : op->targetCRS(); - if (tmpCRS && - dynamic_cast<const crs::GeographicCRS *>(tmpCRS.get())) { - res.emplace_back(i == 0 ? op : op->inverse()); - } - } - if (!res.empty()) - break; - } - - return res; -} - -// --------------------------------------------------------------------------- - static CoordinateOperationNNPtr createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS) { @@ -12723,88 +12691,32 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase( doFilterAndCheckPerfectOp = false; bool sameGeodeticDatum = false; - if (geodSrc && geodDst) { + + if (vertSrc || vertDst) { + createOperationsFromDatabaseWithVertCRS(sourceCRS, targetCRS, + context, geogSrc, geogDst, + vertSrc, vertDst, res); + } else if (geodSrc && geodDst) { + const auto &srcDatum = geodSrc->datum(); const auto &dstDatum = geodDst->datum(); sameGeodeticDatum = srcDatum != nullptr && dstDatum != nullptr && srcDatum->_isEquivalentTo( dstDatum.get(), util::IComparable::Criterion::EQUIVALENT); - } - // NAD83 only exists in 2D version in EPSG, so if it has been - // promotted to 3D, when researching a vertical to geog - // transformation, try to down cast to 2D. - const auto geog3DToVertTryThroughGeog2D = [&res, &context]( - const crs::GeographicCRS *geogSrcIn, - const crs::VerticalCRS *vertDstIn, - const crs::CRSNNPtr &targetCRSIn) { - if (res.empty() && geogSrcIn && vertDstIn && - geogSrcIn->coordinateSystem()->axisList().size() == 3 && - geogSrcIn->datum()) { - const auto &authFactory = - context.context->getAuthorityFactory(); - const auto candidatesSrcGeod(findCandidateGeodCRSForDatum( - authFactory, geogSrcIn->datum().get())); - for (const auto &candidate : candidatesSrcGeod) { - auto geogCandidate = - util::nn_dynamic_pointer_cast<crs::GeographicCRS>( - candidate); - if (geogCandidate && - geogCandidate->coordinateSystem()->axisList().size() == - 2) { - res = findOpsInRegistryDirect( - NN_NO_CHECK(geogCandidate), targetCRSIn, - context.context); - if (res.empty()) { - res = applyInverse(findOpsInRegistryDirect( - targetCRSIn, NN_NO_CHECK(geogCandidate), - context.context)); - } - break; - } - } - return true; - } - return false; - }; - - if (geog3DToVertTryThroughGeog2D(geogSrc, vertDst, targetCRS)) { - // do nothing - } else if (geog3DToVertTryThroughGeog2D(geogDst, vertSrc, sourceCRS)) { - res = applyInverse(res); - } - - // There's no direct transformation from NAVD88 height to WGS84, - // so try to research all transformations from NAVD88 to another - // intermediate GeographicCRS. - if (res.empty() && - !context.inCreateOperationsGeogToVertWithAlternativeGeog && - geogSrc && vertDst) { - res = createOperationsGeogToVertWithAlternativeGeog( - sourceCRS, targetCRS, context); - } else if (res.empty() && - !context.inCreateOperationsGeogToVertWithAlternativeGeog && - geogDst && vertSrc) { - res = applyInverse(createOperationsGeogToVertWithAlternativeGeog( - targetCRS, sourceCRS, context)); - } - - if (res.empty() && !sameGeodeticDatum && - !context.inCreateOperationsWithDatumPivotAntiRecursion && geodSrc && - geodDst) { - // If we still didn't find a transformation, and that the source - // and target are GeodeticCRS, then go through their underlying - // datum to find potential transformations between other - // GeodeticRSs - // that are made of those datum - // The typical example is if transforming between two - // GeographicCRS, - // but transformations are only available between their - // corresponding geocentric CRS. - const auto &srcDatum = geodSrc->datum(); - const auto &dstDatum = geodDst->datum(); - if (srcDatum != nullptr && dstDatum != nullptr) { + if (res.empty() && !sameGeodeticDatum && + !context.inCreateOperationsWithDatumPivotAntiRecursion && + srcDatum != nullptr && dstDatum != nullptr) { + // If we still didn't find a transformation, and that the source + // and target are GeodeticCRS, then go through their underlying + // datum to find potential transformations between other + // GeodeticRSs + // that are made of those datum + // The typical example is if transforming between two + // GeographicCRS, + // but transformations are only available between their + // corresponding geocentric CRS. createOperationsWithDatumPivot(res, sourceCRS, targetCRS, geodSrc, geodDst, context); doFilterAndCheckPerfectOp = !res.empty(); @@ -12851,6 +12763,222 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase( // --------------------------------------------------------------------------- +static std::vector<crs::CRSNNPtr> +findCandidateVertCRSForDatum(const io::AuthorityFactoryPtr &authFactory, + const datum::VerticalReferenceFrame *datum) { + std::vector<crs::CRSNNPtr> candidates; + assert(datum); + const auto &ids = datum->identifiers(); + const auto &datumName = datum->nameStr(); + if (!ids.empty()) { + for (const auto &id : ids) { + const auto &authName = *(id->codeSpace()); + const auto &code = id->code(); + if (!authName.empty()) { + auto l_candidates = + authFactory->createVerticalCRSFromDatum(authName, code); + for (const auto &candidate : l_candidates) { + candidates.emplace_back(candidate); + } + } + } + } else if (datumName != "unknown" && datumName != "unnamed") { + auto matches = authFactory->createObjectsFromName( + datumName, + {io::AuthorityFactory::ObjectType::VERTICAL_REFERENCE_FRAME}, false, + 2); + if (matches.size() == 1) { + const auto &match = matches.front(); + if (datum->_isEquivalentTo( + match.get(), util::IComparable::Criterion::EQUIVALENT) && + !match->identifiers().empty()) { + return findCandidateVertCRSForDatum( + authFactory, + dynamic_cast<const datum::VerticalReferenceFrame *>( + match.get())); + } + } + } + return candidates; +} + +// --------------------------------------------------------------------------- + +std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private:: + createOperationsGeogToVertWithIntermediateVert( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + const crs::VerticalCRS *vertDst, Private::Context &context) { + + std::vector<CoordinateOperationNNPtr> res; + + struct AntiRecursionGuard { + Context &context; + + explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) { + assert(!context.inCreateOperationsGeogToVertWithIntermediateVert); + context.inCreateOperationsGeogToVertWithIntermediateVert = true; + } + + ~AntiRecursionGuard() { + context.inCreateOperationsGeogToVertWithIntermediateVert = false; + } + }; + AntiRecursionGuard guard(context); + const auto &authFactory = context.context->getAuthorityFactory(); + auto candidatesVert = + findCandidateVertCRSForDatum(authFactory, vertDst->datum().get()); + for (const auto &candidateVert : candidatesVert) { + auto resTmp = createOperations(sourceCRS, candidateVert, context); + if (!resTmp.empty()) { + const auto opsSecond = + createOperations(candidateVert, targetCRS, context); + if (!opsSecond.empty()) { + // The transformation from candidateVert to targetCRS should + // be just a unit change typically, so take only the first one, + // which is likely/hopefully the only one. + for (const auto &opFirst : resTmp) { + if (hasIdentifiers(opFirst)) { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + {opFirst, opsSecond.front()}, + !allowEmptyIntersection)); + } + } + if (!res.empty()) + break; + } + } + } + + return res; +} + +// --------------------------------------------------------------------------- + +std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private:: + createOperationsGeogToVertWithAlternativeGeog( + const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS + const crs::CRSNNPtr &targetCRS, // vertical CRS + Private::Context &context) { + + std::vector<CoordinateOperationNNPtr> res; + + struct AntiRecursionGuard { + Context &context; + + explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) { + assert(!context.inCreateOperationsGeogToVertWithAlternativeGeog); + context.inCreateOperationsGeogToVertWithAlternativeGeog = true; + } + + ~AntiRecursionGuard() { + context.inCreateOperationsGeogToVertWithAlternativeGeog = false; + } + }; + AntiRecursionGuard guard(context); + + for (int i = 0; i < 2; i++) { + + // Generally EPSG has operations from GeogCrs to VertCRS + auto ops = + i == 0 ? findOpsInRegistryDirectTo(targetCRS, context.context) + : findOpsInRegistryDirectFrom(targetCRS, context.context); + + for (const auto &op : ops) { + const auto tmpCRS = i == 0 ? op->sourceCRS() : op->targetCRS(); + if (tmpCRS && + dynamic_cast<const crs::GeographicCRS *>(tmpCRS.get())) { + res.emplace_back(i == 0 ? op : op->inverse()); + } + } + if (!res.empty()) + break; + } + + return res; +} + +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private:: + createOperationsFromDatabaseWithVertCRS( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::GeographicCRS *geogSrc, + const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc, + const crs::VerticalCRS *vertDst, + std::vector<CoordinateOperationNNPtr> &res) { + + // Typically to transform from "NAVD88 height (ftUS)" to a geog CRS + // by using transformations of "NAVD88 height" (metre) to that geog CRS + if (res.empty() && + !context.inCreateOperationsGeogToVertWithIntermediateVert && geogSrc && + vertDst) { + res = createOperationsGeogToVertWithIntermediateVert( + sourceCRS, targetCRS, vertDst, context); + } else if (res.empty() && + !context.inCreateOperationsGeogToVertWithIntermediateVert && + geogDst && vertSrc) { + res = applyInverse(createOperationsGeogToVertWithIntermediateVert( + targetCRS, sourceCRS, vertSrc, context)); + } + + // NAD83 only exists in 2D version in EPSG, so if it has been + // promotted to 3D, when researching a vertical to geog + // transformation, try to down cast to 2D. + const auto geog3DToVertTryThroughGeog2D = [&res, &context]( + const crs::GeographicCRS *geogSrcIn, const crs::VerticalCRS *vertDstIn, + const crs::CRSNNPtr &targetCRSIn) { + if (res.empty() && geogSrcIn && vertDstIn && + geogSrcIn->coordinateSystem()->axisList().size() == 3 && + geogSrcIn->datum()) { + const auto &authFactory = context.context->getAuthorityFactory(); + const auto candidatesSrcGeod(findCandidateGeodCRSForDatum( + authFactory, geogSrcIn->datum().get())); + for (const auto &candidate : candidatesSrcGeod) { + auto geogCandidate = + util::nn_dynamic_pointer_cast<crs::GeographicCRS>( + candidate); + if (geogCandidate && + geogCandidate->coordinateSystem()->axisList().size() == 2) { + res = findOpsInRegistryDirect(NN_NO_CHECK(geogCandidate), + targetCRSIn, context.context); + if (res.empty()) { + res = applyInverse(findOpsInRegistryDirect( + targetCRSIn, NN_NO_CHECK(geogCandidate), + context.context)); + } + break; + } + } + return true; + } + return false; + }; + + if (geog3DToVertTryThroughGeog2D(geogSrc, vertDst, targetCRS)) { + // do nothing + } else if (geog3DToVertTryThroughGeog2D(geogDst, vertSrc, sourceCRS)) { + res = applyInverse(res); + } + + // There's no direct transformation from NAVD88 height to WGS84, + // so try to research all transformations from NAVD88 to another + // intermediate GeographicCRS. + if (res.empty() && + !context.inCreateOperationsGeogToVertWithAlternativeGeog && geogSrc && + vertDst) { + res = createOperationsGeogToVertWithAlternativeGeog(sourceCRS, + targetCRS, context); + } else if (res.empty() && + !context.inCreateOperationsGeogToVertWithAlternativeGeog && + geogDst && vertSrc) { + res = applyInverse(createOperationsGeogToVertWithAlternativeGeog( + targetCRS, sourceCRS, context)); + } +} + +// --------------------------------------------------------------------------- + void CoordinateOperationFactory::Private::createOperationsGeodToGeod( const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, Private::Context &context, const crs::GeodeticCRS *geodSrc, @@ -13346,6 +13474,8 @@ getOps(const CoordinateOperationNNPtr &op) { return {op}; } +// --------------------------------------------------------------------------- + static bool useDifferentTransformationsForSameSourceTarget( const CoordinateOperationNNPtr &opA, const CoordinateOperationNNPtr &opB) { auto subOpsA = getOps(opA); @@ -13393,6 +13523,63 @@ static bool useDifferentTransformationsForSameSourceTarget( return false; } +// --------------------------------------------------------------------------- + +static crs::GeographicCRSPtr +getInterpolationGeogCRS(const CoordinateOperationNNPtr &verticalTransform, + const io::DatabaseContextPtr &dbContext) { + crs::GeographicCRSPtr interpolationGeogCRS; + auto transformationVerticalTransform = + dynamic_cast<const Transformation *>(verticalTransform.get()); + if (transformationVerticalTransform == nullptr) { + const auto concat = dynamic_cast<const ConcatenatedOperation *>( + verticalTransform.get()); + if (concat) { + const auto &steps = concat->operations(); + // Is this change of unit + transformation ? + if (steps.size() == 2 && + dynamic_cast<const Conversion *>(steps[0].get())) { + transformationVerticalTransform = + dynamic_cast<const Transformation *>(steps[1].get()); + } + } + } + if (transformationVerticalTransform && + !transformationVerticalTransform->hasBallparkTransformation()) { + auto interpTransformCRS = + transformationVerticalTransform->interpolationCRS(); + if (interpTransformCRS) { + interpolationGeogCRS = + std::dynamic_pointer_cast<crs::GeographicCRS>( + interpTransformCRS); + } else { + // If no explicit interpolation CRS, then + // this will be the geographic CRS of the + // vertical to geog transformation + interpolationGeogCRS = + std::dynamic_pointer_cast<crs::GeographicCRS>( + transformationVerticalTransform->targetCRS().as_nullable()); + } + } + + if (interpolationGeogCRS) { + if (interpolationGeogCRS->coordinateSystem()->axisList().size() == 3) { + // We need to force the interpolation CRS, which + // will + // frequently be 3D, to 2D to avoid transformations + // between source CRS and interpolation CRS to have + // 3D terms. + interpolationGeogCRS = + interpolationGeogCRS->demoteTo2D(std::string(), dbContext) + .as_nullable(); + } + } + + return interpolationGeogCRS; +} + +// --------------------------------------------------------------------------- + void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, Private::Context &context, const crs::CompoundCRS *compoundSrc, @@ -13440,7 +13627,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( CoordinateOperationContext::GridAvailabilityUse:: IGNORE_GRID_AVAILABILITY; for (const auto &op : verticalTransforms) { - if (!op->identifiers().empty() && dbContext) { + if (hasIdentifiers(op) && dbContext) { bool missingGrid = false; if (!ignoreMissingGrids) { const auto gridsNeeded = op->gridsNeeded(dbContext); @@ -13469,7 +13656,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( bool foundRegisteredTransform = false; foundRegisteredTransformWithAllGridsAvailable = false; for (const auto &op : verticalTransformsTmp) { - if (!op->identifiers().empty() && dbContext) { + if (hasIdentifiers(op) && dbContext) { bool missingGrid = false; if (!ignoreMissingGrids) { const auto gridsNeeded = op->gridsNeeded(dbContext); @@ -13510,43 +13697,9 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( cacheHorizToInterpAndInterpToTarget; for (const auto &verticalTransform : verticalTransforms) { - crs::GeographicCRSPtr interpolationGeogCRS; - auto transformationVerticalTransform = - dynamic_cast<const Transformation *>(verticalTransform.get()); - if (transformationVerticalTransform && - !transformationVerticalTransform->hasBallparkTransformation()) { - auto interpTransformCRS = - transformationVerticalTransform->interpolationCRS(); - if (interpTransformCRS) { - interpolationGeogCRS = - std::dynamic_pointer_cast<crs::GeographicCRS>( - interpTransformCRS); - } else { - // If no explicit interpolation CRS, then - // this will be the geographic CRS of the - // vertical to geog transformation - interpolationGeogCRS = - std::dynamic_pointer_cast<crs::GeographicCRS>( - transformationVerticalTransform->targetCRS() - .as_nullable()); - } - } - + crs::GeographicCRSPtr interpolationGeogCRS = + getInterpolationGeogCRS(verticalTransform, dbContext); if (interpolationGeogCRS) { - if (interpolationGeogCRS->coordinateSystem() - ->axisList() - .size() == 3) { - // We need to force the interpolation CRS, which - // will - // frequently be 3D, to 2D to avoid transformations - // between source CRS and interpolation CRS to have - // 3D terms. - interpolationGeogCRS = - interpolationGeogCRS - ->demoteTo2D(std::string(), dbContext) - .as_nullable(); - } - std::vector<CoordinateOperationNNPtr> srcToInterpOps; std::vector<CoordinateOperationNNPtr> interpToTargetOps; diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 993b4eeb..d3ea8f26 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -4603,6 +4603,30 @@ std::list<crs::GeodeticCRSNNPtr> AuthorityFactory::createGeodeticCRSFromDatum( // --------------------------------------------------------------------------- //! @cond Doxygen_Suppress +std::list<crs::VerticalCRSNNPtr> AuthorityFactory::createVerticalCRSFromDatum( + const std::string &datum_auth_name, const std::string &datum_code) const { + std::string sql( + "SELECT auth_name, code FROM vertical_crs WHERE " + "datum_auth_name = ? AND datum_code = ? AND deprecated = 0"); + ListOfParams params{datum_auth_name, datum_code}; + if (d->hasAuthorityRestriction()) { + sql += " AND auth_name = ?"; + params.emplace_back(d->authority()); + } + auto sqlRes = d->run(sql, params); + std::list<crs::VerticalCRSNNPtr> res; + for (const auto &row : sqlRes) { + const auto &auth_name = row[0]; + const auto &code = row[1]; + res.emplace_back(d->createFactory(auth_name)->createVerticalCRS(code)); + } + return res; +} +//! @endcond + +// --------------------------------------------------------------------------- + +//! @cond Doxygen_Suppress std::list<crs::GeodeticCRSNNPtr> AuthorityFactory::createGeodeticCRSFromEllipsoid( const std::string &ellipsoid_auth_name, const std::string &ellipsoid_code, diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 79541d88..7fbbf70e 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -7165,6 +7165,127 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) { // --------------------------------------------------------------------------- +TEST(operation, compoundCRS_to_geogCRS_with_vertical_unit_change) { + 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"); // NAD83(2011) 3D + + auto listCompoundToGeog = + CoordinateOperationFactory::create()->createOperations(nnSrc, dst, + ctxt); + ASSERT_TRUE(!listCompoundToGeog.empty()); + + // NAD83(2011) + NAVD88 height + auto srcObjCompoundVMetre = createFromUserInput( + "EPSG:6318+5703", authFactory->databaseContext(), false); + auto srcCompoundVMetre = nn_dynamic_pointer_cast<CRS>(srcObjCompoundVMetre); + ASSERT_TRUE(srcCompoundVMetre != nullptr); + auto listCompoundMetreToGeog = + CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcCompoundVMetre), dst, ctxt); + + // Check that we get the same and similar results whether we start from + // regular NAVD88 height or its ftUs variant + ASSERT_EQ(listCompoundToGeog.size(), listCompoundMetreToGeog.size()); + + EXPECT_EQ(listCompoundToGeog[0]->nameStr(), + "Transformation from NAVD88 height (ftUS) to NAVD88 height + " + + listCompoundMetreToGeog[0]->nameStr()); + EXPECT_EQ( + listCompoundToGeog[0]->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5, + authFactory->databaseContext()) + .get()), + replaceAll(listCompoundMetreToGeog[0]->exportToPROJString( + PROJStringFormatter::create( + PROJStringFormatter::Convention::PROJ_5, + authFactory->databaseContext()) + .get()), + "+step +proj=unitconvert +xy_in=deg +xy_out=rad", + "+step +proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad " + "+z_out=m")); + + // Check reverse path + auto listGeogToCompound = + CoordinateOperationFactory::create()->createOperations(dst, nnSrc, + ctxt); + EXPECT_EQ(listGeogToCompound.size(), listCompoundToGeog.size()); +} + +// --------------------------------------------------------------------------- + +TEST( + operation, + compoundCRS_to_geogCRS_with_vertical_unit_change_and_complex_horizontal_change) { + 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("7665"); // WGS84(G1762) 3D + + auto listCompoundToGeog = + CoordinateOperationFactory::create()->createOperations(nnSrc, dst, + ctxt); + + // NAD83(2011) + NAVD88 height + auto srcObjCompoundVMetre = createFromUserInput( + "EPSG:6318+5703", authFactory->databaseContext(), false); + auto srcCompoundVMetre = nn_dynamic_pointer_cast<CRS>(srcObjCompoundVMetre); + ASSERT_TRUE(srcCompoundVMetre != nullptr); + auto listCompoundMetreToGeog = + CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcCompoundVMetre), dst, ctxt); + + // Check that we get the same and similar results whether we start from + // regular NAVD88 height or its ftUs variant + ASSERT_EQ(listCompoundToGeog.size(), listCompoundMetreToGeog.size()); + + ASSERT_GE(listCompoundToGeog.size(), 1U); + + EXPECT_EQ(listCompoundToGeog[0]->nameStr(), + "Transformation from NAVD88 height (ftUS) to NAVD88 height + " + + listCompoundMetreToGeog[0]->nameStr()); + EXPECT_EQ( + listCompoundToGeog[0]->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5, + authFactory->databaseContext()) + .get()), + replaceAll(listCompoundMetreToGeog[0]->exportToPROJString( + PROJStringFormatter::create( + PROJStringFormatter::Convention::PROJ_5, + authFactory->databaseContext()) + .get()), + "+step +proj=unitconvert +xy_in=deg +xy_out=rad", + "+step +proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad " + "+z_out=m")); +} + +// --------------------------------------------------------------------------- + TEST(operation, compoundCRS_from_WKT2_to_geogCRS_3D_context) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); |