From 805b187edd4c9246629e9aeb062118f8c2de2dfe Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 29 Oct 2019 11:48:17 +0100 Subject: Little tweaks in implicit 2D/3D geog conversions in compoundCRS to geogCRS transformations --- src/iso19111/coordinateoperation.cpp | 181 ++++++++++++++++++----------------- src/iso19111/crs.cpp | 2 +- test/unit/test_operation.cpp | 36 ++++--- 3 files changed, 120 insertions(+), 99 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 0f0e216c..4e7ebd38 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -10297,7 +10297,7 @@ struct CoordinateOperationFactory::Private { const CoordinateOperationContextNNPtr &context; bool inCreateOperationsWithDatumPivotAntiRecursion = false; bool inCreateOperationsThroughPreferredHub = false; - bool inCreateOperationsGeogToVertWithIntermediate = false; + bool inCreateOperationsGeogToVertWithAlternativeGeog = false; bool skipHorizontalTransformation = false; Context(const crs::CRSNNPtr &sourceCRSIn, @@ -10330,9 +10330,9 @@ struct CoordinateOperationFactory::Private { Context &context); static std::vector - createOperationsGeogToVertWithIntermediate(const crs::CRSNNPtr &sourceCRS, - const crs::CRSNNPtr &targetCRS, - Context &context); + createOperationsGeogToVertWithAlternativeGeog( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Context &context); static bool hasPerfectAccuracyResult(const std::vector &res, @@ -12327,11 +12327,11 @@ void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub( // --------------------------------------------------------------------------- -std::vector -CoordinateOperationFactory::Private::createOperationsGeogToVertWithIntermediate( - const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS - const crs::CRSNNPtr &targetCRS, // vertical CRS - Private::Context &context) { +std::vector CoordinateOperationFactory::Private:: + createOperationsGeogToVertWithAlternativeGeog( + const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS + const crs::CRSNNPtr &targetCRS, // vertical CRS + Private::Context &context) { std::vector res; @@ -12339,12 +12339,12 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertWithIntermediate( Context &context; explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) { - assert(!context.inCreateOperationsGeogToVertWithIntermediate); - context.inCreateOperationsGeogToVertWithIntermediate = true; + assert(!context.inCreateOperationsGeogToVertWithAlternativeGeog); + context.inCreateOperationsGeogToVertWithAlternativeGeog = true; } ~AntiRecursionGuard() { - context.inCreateOperationsGeogToVertWithIntermediate = false; + context.inCreateOperationsGeogToVertWithAlternativeGeog = false; } }; AntiRecursionGuard guard(context); @@ -12594,69 +12594,63 @@ CoordinateOperationFactory::Private::createOperations( // 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. - if (res.empty() && geogSrc && vertDst && - geogSrc->coordinateSystem()->axisList().size() == 3 && - geogSrc->datum()) { - const auto candidatesSrcGeod(findCandidateGeodCRSForDatum( - authFactory, geogSrc->datum().get())); - for (const auto &candidate : candidatesSrcGeod) { - auto geogCandidate = - util::nn_dynamic_pointer_cast( - candidate); - if (geogCandidate && - geogCandidate->coordinateSystem()->axisList().size() == - 2) { - res = - findOpsInRegistryDirect(NN_NO_CHECK(geogCandidate), - targetCRS, context.context); - if (res.empty()) { - res = applyInverse(findOpsInRegistryDirect( - targetCRS, NN_NO_CHECK(geogCandidate), - context.context)); - } - break; - } - } - } else if (res.empty() && geogDst && vertSrc && - geogDst->coordinateSystem()->axisList().size() == 3 && - geogDst->datum()) { - const auto candidatesDstGeod(findCandidateGeodCRSForDatum( - authFactory, geogDst->datum().get())); - for (const auto &candidate : candidatesDstGeod) { - auto geogCandidate = - util::nn_dynamic_pointer_cast( - candidate); - if (geogCandidate && - geogCandidate->coordinateSystem()->axisList().size() == - 2) { - res = findOpsInRegistryDirect( - sourceCRS, NN_NO_CHECK(geogCandidate), - context.context); - if (res.empty()) { - res = applyInverse(findOpsInRegistryDirect( - NN_NO_CHECK(geogCandidate), sourceCRS, - context.context)); + // transformation, try to down cast to 2D. + const auto geog3DToVertTryThroughGeog2D = [&res, &authFactory, + &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 candidatesSrcGeod(findCandidateGeodCRSForDatum( + authFactory, geogSrcIn->datum().get())); + for (const auto &candidate : candidatesSrcGeod) { + auto geogCandidate = + util::nn_dynamic_pointer_cast( + 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; } - 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.inCreateOperationsGeogToVertWithIntermediate && + !context.inCreateOperationsGeogToVertWithAlternativeGeog && geogSrc && vertDst) { - res = createOperationsGeogToVertWithIntermediate( + res = createOperationsGeogToVertWithAlternativeGeog( sourceCRS, targetCRS, context); } else if (res.empty() && - !context.inCreateOperationsGeogToVertWithIntermediate && + !context + .inCreateOperationsGeogToVertWithAlternativeGeog && geogDst && vertSrc) { - res = applyInverse(createOperationsGeogToVertWithIntermediate( - targetCRS, sourceCRS, context)); + res = + applyInverse(createOperationsGeogToVertWithAlternativeGeog( + targetCRS, sourceCRS, context)); } if (res.empty() && !sameGeodeticDatum && @@ -13219,6 +13213,10 @@ CoordinateOperationFactory::Private::createOperations( createOperations(componentsSrc[0], targetCRS, context); } std::vector verticalTransforms; + + const auto dbContext = + authFactory ? authFactory->databaseContext().as_nullable() + : nullptr; if (componentsSrc.size() >= 2 && componentsSrc[1]->extractVerticalCRS()) { @@ -13237,18 +13235,24 @@ CoordinateOperationFactory::Private::createOperations( }; SetSkipHorizontalTransform setSkipHorizontalTransform(context); - verticalTransforms = - createOperations(componentsSrc[1], targetCRS, context); + verticalTransforms = createOperations( + componentsSrc[1], + targetCRS->promoteTo3D(std::string(), dbContext), context); bool foundRegisteredTransformWithAllGridsAvailable = false; + const bool ignoreMissingGrids = + context.context->getGridAvailabilityUse() == + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY; for (const auto &op : verticalTransforms) { - if (!op->identifiers().empty() && authFactory) { + if (!op->identifiers().empty() && dbContext) { bool missingGrid = false; - const auto gridsNeeded = - op->gridsNeeded(authFactory->databaseContext()); - for (const auto &gridDesc : gridsNeeded) { - if (!gridDesc.available) { - missingGrid = true; - break; + if (!ignoreMissingGrids) { + const auto gridsNeeded = op->gridsNeeded(dbContext); + for (const auto &gridDesc : gridsNeeded) { + if (!gridDesc.available) { + missingGrid = true; + break; + } } } if (!missingGrid) { @@ -13264,18 +13268,23 @@ CoordinateOperationFactory::Private::createOperations( geogDst, util::IComparable::Criterion::EQUIVALENT) && !srcGeogCRS->is2DPartOf3D(NN_NO_CHECK(geogDst))) { auto verticalTransformsTmp = createOperations( - componentsSrc[1], NN_NO_CHECK(srcGeogCRS), context); + componentsSrc[1], + NN_NO_CHECK(srcGeogCRS) + ->promoteTo3D(std::string(), dbContext), + context); bool foundRegisteredTransform = false; foundRegisteredTransformWithAllGridsAvailable = false; for (const auto &op : verticalTransformsTmp) { - if (!op->identifiers().empty() && authFactory) { + if (!op->identifiers().empty() && dbContext) { bool missingGrid = false; - const auto gridsNeeded = - op->gridsNeeded(authFactory->databaseContext()); - for (const auto &gridDesc : gridsNeeded) { - if (!gridDesc.available) { - missingGrid = true; - break; + if (!ignoreMissingGrids) { + const auto gridsNeeded = + op->gridsNeeded(dbContext); + for (const auto &gridDesc : gridsNeeded) { + if (!gridDesc.available) { + missingGrid = true; + break; + } } } foundRegisteredTransform = true; @@ -13335,10 +13344,6 @@ CoordinateOperationFactory::Private::createOperations( if (interpolationGeogCRS->coordinateSystem() ->axisList() .size() == 3) { - io::DatabaseContextPtr dbContext; - if (authFactory) { - dbContext = authFactory->databaseContext(); - } // We need to force the interpolation CRS, which // will // frequently be 3D, to 2D to avoid transformations @@ -13367,7 +13372,8 @@ CoordinateOperationFactory::Private::createOperations( componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS), context); interpToTargetOps = createOperations( - NN_NO_CHECK(interpolationGeogCRS), targetCRS, + NN_NO_CHECK(interpolationGeogCRS), + targetCRS->demoteTo2D(std::string(), dbContext), context); cacheHorizToInterpAndInterpToTarget[key] = PairOfTransforms(srcToInterpOps, @@ -13380,9 +13386,10 @@ CoordinateOperationFactory::Private::createOperations( srcToInterpOps = createOperations( componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS), context); - interpToTargetOps = - createOperations(NN_NO_CHECK(interpolationGeogCRS), - targetCRS, context); + interpToTargetOps = createOperations( + NN_NO_CHECK(interpolationGeogCRS), + targetCRS->demoteTo2D(std::string(), dbContext), + context); } for (const auto &srcToInterp : srcToInterpOps) { diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 5d3af369..f16ed7cf 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -2170,7 +2170,7 @@ GeographicCRS::demoteTo2D(const std::string &newName, const auto &axisList = coordinateSystem()->axisList(); if (axisList.size() == 3) { const auto &l_identifiers = identifiers(); - // First check if there is a Geographic 3D CRS in the database + // First check if there is a Geographic 2D CRS in the database // of the same name. // This is the common practice in the EPSG dataset. if (dbContext && l_identifiers.size() == 1) { diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 739d8ec3..79541d88 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -7129,6 +7129,11 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) { 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 + NAVD88 height auto srcObj = createFromUserInput("EPSG:4269+5703", authFactory->databaseContext(), false); @@ -7136,17 +7141,26 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) { ASSERT_TRUE(src != nullptr); auto nnSrc = NN_NO_CHECK(src); auto dst = authFactory->createCoordinateReferenceSystem("4269"); // NAD83 - dst = dst->promoteTo3D(std::string(), authFactory->databaseContext()); - { - auto list = CoordinateOperationFactory::create()->createOperations( - nnSrc, dst, ctxt); - ASSERT_GE(list.size(), 8U); - } - { - auto list = CoordinateOperationFactory::create()->createOperations( - dst, nnSrc, ctxt); - ASSERT_GE(list.size(), 8U); - } + + auto listCompoundToGeog2D = + CoordinateOperationFactory::create()->createOperations(nnSrc, dst, + ctxt); + // The checked value is not that important, but in case this changes, + // likely due to a EPSG upgrade, worth checking + ASSERT_EQ(listCompoundToGeog2D.size(), 469U); + + auto listGeog2DToCompound = + CoordinateOperationFactory::create()->createOperations(dst, nnSrc, + ctxt); + ASSERT_EQ(listGeog2DToCompound.size(), listCompoundToGeog2D.size()); + + auto listCompoundToGeog3D = + CoordinateOperationFactory::create()->createOperations( + nnSrc, + dst->promoteTo3D(std::string(), authFactory->databaseContext()), + ctxt); + // It includes a trailing ballpart transformation + ASSERT_EQ(listCompoundToGeog3D.size(), listCompoundToGeog2D.size() + 1); } // --------------------------------------------------------------------------- -- cgit v1.2.3 From eed030d96b1b8142e1a1c236555054c32a143e93 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 29 Oct 2019 15:34:10 +0100 Subject: createOperations(): split gigantic method into many smaller ones. No functional change expected --- src/iso19111/coordinateoperation.cpp | 2229 ++++++++++++++++++---------------- 1 file changed, 1205 insertions(+), 1024 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 4e7ebd38..c7581642 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -3197,7 +3197,7 @@ ConversionNNPtr Conversion::createBonne(const util::PropertyMap &properties, * (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::9834) * * \warning The PROJ cea computation code would select the ellipsoidal form if - * a non-spherical ellipsoid is used for the base GeographicalCRS. + * a non-spherical ellipsoid is used for the base GeographicCRS. * * @param properties See \ref general_properties of the conversion. If the name * is not provided, it is automatically set. @@ -10312,6 +10312,86 @@ struct CoordinateOperationFactory::Private { const crs::CRSNNPtr &targetCRS, Context &context); private: + static constexpr bool allowEmptyIntersection = true; + + static void createOperationsFromProj4Ext( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + const crs::BoundCRS *boundSrc, const crs::BoundCRS *boundDst, + std::vector &res); + + static bool createOperationsFromDatabase( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::GeodeticCRS *geodSrc, + const crs::GeodeticCRS *geodDst, const crs::GeographicCRS *geogSrc, + const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc, + const crs::VerticalCRS *vertDst, + std::vector &res); + + static void createOperationsGeodToGeod( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::GeodeticCRS *geodSrc, + const crs::GeodeticCRS *geodDst, + std::vector &res); + + static void createOperationsDerivedTo( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::DerivedCRS *derivedSrc, + std::vector &res); + + static void createOperationsBoundToGeog( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::BoundCRS *boundSrc, + const crs::GeographicCRS *geogDst, + std::vector &res); + + static void createOperationsBoundToVert( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::BoundCRS *boundSrc, + const crs::VerticalCRS *vertDst, + std::vector &res); + + static void createOperationsVertToVert( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::VerticalCRS *vertSrc, + const crs::VerticalCRS *vertDst, + std::vector &res); + + static void createOperationsVertToGeog( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::VerticalCRS *vertSrc, + const crs::GeographicCRS *geogDst, + std::vector &res); + + static void createOperationsBoundToBound( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::BoundCRS *boundSrc, + const crs::BoundCRS *boundDst, + std::vector &res); + + static void createOperationsCompoundToGeog( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::CompoundCRS *compoundSrc, + const crs::GeographicCRS *geogDst, + std::vector &res); + + static void createOperationsCompoundToGeod( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::CompoundCRS *compoundSrc, + const crs::GeodeticCRS *geodDst, + std::vector &res); + + static void createOperationsCompoundToCompound( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::CompoundCRS *compoundSrc, + const crs::CompoundCRS *compoundDst, + std::vector &res); + + static void createOperationsBoundToCompound( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::BoundCRS *boundSrc, + const crs::CompoundCRS *compoundDst, + std::vector &res); + static std::vector createOperationsGeogToGeog( std::vector &res, const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, @@ -11544,6 +11624,9 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final MyPROJStringExportableHorizVerticalHorizPROJBased:: ~MyPROJStringExportableHorizVerticalHorizPROJBased() = default; + +//! @endcond + } // namespace operation NS_PROJ_END @@ -11558,6 +11641,8 @@ template<> nnprimeMeridian()->longitude(); const auto &dst_pm = geogDst->primeMeridian()->longitude(); @@ -11885,12 +11969,8 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog( return res; } -//! @endcond - // --------------------------------------------------------------------------- -//! @cond Doxygen_Suppress - static bool hasIdentifiers(const CoordinateOperationNNPtr &op) { if (!op->identifiers().empty()) { return true; @@ -11905,12 +11985,9 @@ static bool hasIdentifiers(const CoordinateOperationNNPtr &op) { } return false; } -//! @endcond // --------------------------------------------------------------------------- -//! @cond Doxygen_Suppress - static std::vector findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory, const datum::GeodeticReferenceFrame *datum) { @@ -12032,7 +12109,6 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( objectAsStr(sourceCRS.get()) + "," + objectAsStr(targetCRS.get()) + ")"); #endif - const bool allowEmptyIntersection = true; struct CreateOperationsWithDatumPivotAntiRecursion { Context &context; @@ -12298,7 +12374,6 @@ void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub( } } - const bool allowEmptyIntersection = true; for (const auto &intermCRS : candidatesIntermCRS) { #ifdef DEBUG EnterDebugLevel loopLevel; @@ -12390,12 +12465,8 @@ createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS, sourceCRS, targetCRS, 0.0, 0.0, 0.0, {})); } -//! @endcond - // --------------------------------------------------------------------------- -//! @cond Doxygen_Suppress - bool CoordinateOperationFactory::Private::hasPerfectAccuracyResult( const std::vector &res, const Context &context) { auto resTmp = FilterResults(res, context.context, context.sourceCRS, @@ -12410,72 +12481,8 @@ bool CoordinateOperationFactory::Private::hasPerfectAccuracyResult( return false; } -//! @endcond - -// --------------------------------------------------------------------------- - -//! @cond Doxygen_Suppress - -static std::vector -getOps(const CoordinateOperationNNPtr &op) { - auto concatenated = dynamic_cast(op.get()); - if (concatenated) - return concatenated->operations(); - return {op}; -} - -static bool useDifferentTransformationsForSameSourceTarget( - const CoordinateOperationNNPtr &opA, const CoordinateOperationNNPtr &opB) { - auto subOpsA = getOps(opA); - auto subOpsB = getOps(opB); - for (const auto &subOpA : subOpsA) { - if (!dynamic_cast(subOpA.get())) - continue; - if (subOpA->sourceCRS()->nameStr() == "unknown" || - subOpA->targetCRS()->nameStr() == "unknown") - continue; - for (const auto &subOpB : subOpsB) { - if (!dynamic_cast(subOpB.get())) - continue; - if (subOpB->sourceCRS()->nameStr() == "unknown" || - subOpB->targetCRS()->nameStr() == "unknown") - continue; - - if (subOpA->sourceCRS()->nameStr() == - subOpB->sourceCRS()->nameStr() && - subOpA->targetCRS()->nameStr() == - subOpB->targetCRS()->nameStr()) { - if (starts_with(subOpA->nameStr(), NULL_GEOGRAPHIC_OFFSET) && - starts_with(subOpB->nameStr(), NULL_GEOGRAPHIC_OFFSET)) { - continue; - } - - if (!subOpA->isEquivalentTo(subOpB.get())) { - return true; - } - } else if (subOpA->sourceCRS()->nameStr() == - subOpB->targetCRS()->nameStr() && - subOpA->targetCRS()->nameStr() == - subOpB->sourceCRS()->nameStr()) { - if (starts_with(subOpA->nameStr(), NULL_GEOGRAPHIC_OFFSET) && - starts_with(subOpB->nameStr(), NULL_GEOGRAPHIC_OFFSET)) { - continue; - } - - if (!subOpA->isEquivalentTo(subOpB->inverse().get())) { - return true; - } - } - } - } - return false; -} - -//! @endcond - // --------------------------------------------------------------------------- -//! @cond Doxygen_Suppress std::vector CoordinateOperationFactory::Private::createOperations( const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, @@ -12489,7 +12496,6 @@ CoordinateOperationFactory::Private::createOperations( #endif std::vector res; - const bool allowEmptyIntersection = true; auto boundSrc = dynamic_cast(sourceCRS.get()); auto boundDst = dynamic_cast(targetCRS.get()); @@ -12501,49 +12507,8 @@ CoordinateOperationFactory::Private::createOperations( ? boundDst->baseCRS()->getExtensionProj4() : targetCRS->getExtensionProj4(); if (!sourceProj4Ext.empty() || !targetProj4Ext.empty()) { - - auto sourceProjExportable = - dynamic_cast( - boundSrc ? boundSrc : sourceCRS.get()); - auto targetProjExportable = - dynamic_cast( - boundDst ? boundDst : targetCRS.get()); - if (!sourceProjExportable) { - throw InvalidOperation("Source CRS is not PROJ exportable"); - } - if (!targetProjExportable) { - throw InvalidOperation("Target CRS is not PROJ exportable"); - } - auto projFormatter = io::PROJStringFormatter::create(); - projFormatter->setCRSExport(true); - projFormatter->setLegacyCRSToCRSContext(true); - projFormatter->startInversion(); - sourceProjExportable->_exportToPROJString(projFormatter.get()); - auto geogSrc = - dynamic_cast(sourceCRS.get()); - if (geogSrc) { - auto tmpFormatter = io::PROJStringFormatter::create(); - geogSrc->addAngularUnitConvertAndAxisSwap(tmpFormatter.get()); - projFormatter->ingestPROJString(tmpFormatter->toString()); - } - - projFormatter->stopInversion(); - - targetProjExportable->_exportToPROJString(projFormatter.get()); - auto geogDst = - dynamic_cast(targetCRS.get()); - if (geogDst) { - auto tmpFormatter = io::PROJStringFormatter::create(); - geogDst->addAngularUnitConvertAndAxisSwap(tmpFormatter.get()); - projFormatter->ingestPROJString(tmpFormatter->toString()); - } - - const auto PROJString = projFormatter->toString(); - auto properties = util::PropertyMap().set( - common::IdentifiedObject::NAME_KEY, - buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr())); - res.emplace_back(SingleOperation::createPROJBased( - properties, PROJString, sourceCRS, targetCRS, {})); + createOperationsFromProj4Ext(sourceCRS, targetCRS, boundSrc, boundDst, + res); return res; } @@ -12566,216 +12531,17 @@ CoordinateOperationFactory::Private::createOperations( !derivedDst->baseCRS()->_isEquivalentTo( sourceCRS.get(), util::IComparable::Criterion::EQUIVALENT))) { - bool doFilterAndCheckPerfectOp = true; - res = findOpsInRegistryDirect(sourceCRS, targetCRS, context.context); - if (!sourceCRS->_isEquivalentTo(targetCRS.get())) { - auto resFromInverse = applyInverse( - findOpsInRegistryDirect(targetCRS, sourceCRS, context.context)); - res.insert(res.end(), resFromInverse.begin(), resFromInverse.end()); - - // If we get at least a result with perfect accuracy, do not - // bother generating synthetic transforms. - if (hasPerfectAccuracyResult(res, context)) { - return res; - } - - doFilterAndCheckPerfectOp = false; - - bool sameGeodeticDatum = false; - 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, &authFactory, - &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 candidatesSrcGeod(findCandidateGeodCRSForDatum( - authFactory, geogSrcIn->datum().get())); - for (const auto &candidate : candidatesSrcGeod) { - auto geogCandidate = - util::nn_dynamic_pointer_cast( - 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) { - createOperationsWithDatumPivot(res, sourceCRS, targetCRS, - geodSrc, geodDst, context); - doFilterAndCheckPerfectOp = !res.empty(); - } - } - - // NAD27 to NAD83 has tens of results already. No need to look - // for a pivot - if (!sameGeodeticDatum && - ((res.empty() && - context.context->getAllowUseIntermediateCRS() == - CoordinateOperationContext::IntermediateCRSUse:: - IF_NO_DIRECT_TRANSFORMATION) || - context.context->getAllowUseIntermediateCRS() == - CoordinateOperationContext::IntermediateCRSUse::ALWAYS || - getenv("PROJ_FORCE_SEARCH_PIVOT"))) { - auto resWithIntermediate = findsOpsInRegistryWithIntermediate( - sourceCRS, targetCRS, context.context); - res.insert(res.end(), resWithIntermediate.begin(), - resWithIntermediate.end()); - doFilterAndCheckPerfectOp = !res.empty(); - - // If transforming from/to WGS84 (Gxxxx), try through 'neutral' - // WGS84 - if (res.empty() && geodSrc && geodDst && - !context.inCreateOperationsThroughPreferredHub && - !context.inCreateOperationsWithDatumPivotAntiRecursion) { - createOperationsThroughPreferredHub( - res, sourceCRS, targetCRS, geodSrc, geodDst, context); - doFilterAndCheckPerfectOp = !res.empty(); - } - } - } - - if (doFilterAndCheckPerfectOp) { - // If we get at least a result with perfect accuracy, do not bother - // generating synthetic transforms. - if (hasPerfectAccuracyResult(res, context)) { - return res; - } + if (createOperationsFromDatabase(sourceCRS, targetCRS, context, geodSrc, + geodDst, geogSrc, geogDst, vertSrc, + vertDst, res)) { + return res; } } // Special case if both CRS are geodetic if (geodSrc && geodDst && !derivedSrc && !derivedDst) { - - if (geodSrc->ellipsoid()->celestialBody() != - geodDst->ellipsoid()->celestialBody()) { - throw util::UnsupportedOperationException( - "Source and target ellipsoid do not belong to the same " - "celestial body"); - } - - if (geogSrc && geogDst) { - return createOperationsGeogToGeog(res, sourceCRS, targetCRS, - geogSrc, geogDst); - } - - const bool isSrcGeocentric = geodSrc->isGeocentric(); - const bool isSrcGeographic = geogSrc != nullptr; - const bool isTargetGeocentric = geodDst->isGeocentric(); - const bool isTargetGeographic = geogDst != nullptr; - if (((isSrcGeocentric && isTargetGeographic) || - (isSrcGeographic && isTargetGeocentric)) && - geodSrc->datum() != nullptr && geodDst->datum() != nullptr) { - - // Same datum ? - if (geodSrc->datum()->_isEquivalentTo( - geodDst->datum().get(), - util::IComparable::Criterion::EQUIVALENT)) { - res.emplace_back(Conversion::createGeographicGeocentric( - sourceCRS, targetCRS)); - } else if (isSrcGeocentric) { - std::string interm_crs_name(geogDst->nameStr()); - interm_crs_name += " (geocentric)"; - auto interm_crs = util::nn_static_pointer_cast( - crs::GeodeticCRS::create( - addDomains(util::PropertyMap().set( - common::IdentifiedObject::NAME_KEY, - interm_crs_name), - geogDst), - NN_NO_CHECK(geogDst->datum()), - NN_CHECK_ASSERT( - util::nn_dynamic_pointer_cast( - geodSrc->coordinateSystem())))); - auto opFirst = - createBallparkGeocentricTranslation(sourceCRS, interm_crs); - auto opSecond = Conversion::createGeographicGeocentric( - interm_crs, targetCRS); - res.emplace_back(ConcatenatedOperation::createComputeMetadata( - {opFirst, opSecond}, !allowEmptyIntersection)); - } else { - return applyInverse( - createOperations(targetCRS, sourceCRS, context)); - } - - return res; - } - - if (isSrcGeocentric && isTargetGeocentric) { - res.emplace_back( - createBallparkGeocentricTranslation(sourceCRS, targetCRS)); - return res; - } - - // Transformation between two geodetic systems of unknown type - // This should normally not be triggered with "standard" CRS - res.emplace_back(createGeodToGeodPROJBased(sourceCRS, targetCRS)); + createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc, + geodDst, res); return res; } @@ -12783,23 +12549,8 @@ CoordinateOperationFactory::Private::createOperations( // deriving conversion, with transforms from its baseCRS to the // targetCRS if (derivedSrc) { - auto opFirst = derivedSrc->derivingConversion()->inverse(); - // Small optimization if the targetCRS is the baseCRS of the source - // derivedCRS. - if (derivedSrc->baseCRS()->_isEquivalentTo( - targetCRS.get(), util::IComparable::Criterion::EQUIVALENT)) { - res.emplace_back(opFirst); - return res; - } - auto opsSecond = - createOperations(derivedSrc->baseCRS(), targetCRS, context); - for (const auto &opSecond : opsSecond) { - try { - res.emplace_back(ConcatenatedOperation::createComputeMetadata( - {opFirst, opSecond}, !allowEmptyIntersection)); - } catch (const InvalidOperationEmptyIntersection &) { - } - } + createOperationsDerivedTo(sourceCRS, targetCRS, context, derivedSrc, + res); return res; } @@ -12809,196 +12560,9 @@ CoordinateOperationFactory::Private::createOperations( } if (boundSrc && geogDst) { - const auto &hubSrc = boundSrc->hubCRS(); - auto hubSrcGeog = - dynamic_cast(hubSrc.get()); - auto geogCRSOfBaseOfBoundSrc = - boundSrc->baseCRS()->extractGeographicCRS(); - bool triedBoundCrsToGeogCRSSameAsHubCRS = false; - // Is it: boundCRS to a geogCRS that is the same as the hubCRS ? - if (hubSrcGeog && geogCRSOfBaseOfBoundSrc && - (hubSrcGeog->_isEquivalentTo( - geogDst, util::IComparable::Criterion::EQUIVALENT) || - hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst)))) { - triedBoundCrsToGeogCRSSameAsHubCRS = true; - if (boundSrc->baseCRS() == geogCRSOfBaseOfBoundSrc) { - // Optimization to avoid creating a useless concatenated - // operation - res.emplace_back(boundSrc->transformation()); - return res; - } - auto opsFirst = - createOperations(boundSrc->baseCRS(), - NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context); - if (!opsFirst.empty()) { - for (const auto &opFirst : opsFirst) { - try { - res.emplace_back( - ConcatenatedOperation::createComputeMetadata( - {opFirst, boundSrc->transformation()}, - !allowEmptyIntersection)); - } catch (const InvalidOperationEmptyIntersection &) { - } - } - if (!res.empty()) { - return res; - } - } - // If the datum are equivalent, this is also fine - } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() && - geogDst->datum() && - hubSrcGeog->datum()->_isEquivalentTo( - geogDst->datum().get(), - util::IComparable::Criterion::EQUIVALENT)) { - auto opsFirst = - createOperations(boundSrc->baseCRS(), - NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context); - auto opsLast = createOperations(hubSrc, targetCRS, context); - if (!opsFirst.empty() && !opsLast.empty()) { - for (const auto &opFirst : opsFirst) { - for (const auto &opLast : opsLast) { - try { - res.emplace_back( - ConcatenatedOperation::createComputeMetadata( - {opFirst, boundSrc->transformation(), - opLast}, - !allowEmptyIntersection)); - } catch (const InvalidOperationEmptyIntersection &) { - } - } - } - if (!res.empty()) { - return res; - } - } - // Consider WGS 84 and NAD83 as equivalent in that context if the - // geogCRSOfBaseOfBoundSrc ellipsoid is Clarke66 (for NAD27) - // Case of "+proj=latlong +ellps=clrk66 - // +nadgrids=ntv1_can.dat,conus" - // to "+proj=latlong +datum=NAD83" - } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() && - geogDst->datum() && - geogCRSOfBaseOfBoundSrc->ellipsoid()->_isEquivalentTo( - datum::Ellipsoid::CLARKE_1866.get(), - util::IComparable::Criterion::EQUIVALENT) && - hubSrcGeog->datum()->_isEquivalentTo( - datum::GeodeticReferenceFrame::EPSG_6326.get(), - util::IComparable::Criterion::EQUIVALENT) && - geogDst->datum()->_isEquivalentTo( - datum::GeodeticReferenceFrame::EPSG_6269.get(), - util::IComparable::Criterion::EQUIVALENT)) { - auto nnGeogCRSOfBaseOfBoundSrc = - NN_NO_CHECK(geogCRSOfBaseOfBoundSrc); - if (boundSrc->baseCRS()->_isEquivalentTo( - nnGeogCRSOfBaseOfBoundSrc.get(), - util::IComparable::Criterion::EQUIVALENT)) { - auto transf = boundSrc->transformation()->shallowClone(); - transf->setProperties(util::PropertyMap().set( - common::IdentifiedObject::NAME_KEY, - buildTransfName(boundSrc->baseCRS()->nameStr(), - targetCRS->nameStr()))); - transf->setCRSs(boundSrc->baseCRS(), targetCRS, nullptr); - res.emplace_back(transf); - return res; - } else { - auto opsFirst = createOperations( - boundSrc->baseCRS(), nnGeogCRSOfBaseOfBoundSrc, context); - auto transf = boundSrc->transformation()->shallowClone(); - transf->setProperties(util::PropertyMap().set( - common::IdentifiedObject::NAME_KEY, - buildTransfName(nnGeogCRSOfBaseOfBoundSrc->nameStr(), - targetCRS->nameStr()))); - transf->setCRSs(nnGeogCRSOfBaseOfBoundSrc, targetCRS, nullptr); - if (!opsFirst.empty()) { - for (const auto &opFirst : opsFirst) { - try { - res.emplace_back( - ConcatenatedOperation::createComputeMetadata( - {opFirst, transf}, - !allowEmptyIntersection)); - } catch (const InvalidOperationEmptyIntersection &) { - } - } - if (!res.empty()) { - return res; - } - } - } - } - - if (hubSrcGeog && - hubSrcGeog->_isEquivalentTo( - geogDst, util::IComparable::Criterion::EQUIVALENT) && - dynamic_cast(boundSrc->baseCRS().get())) { - res.emplace_back(boundSrc->transformation()); - return res; - } - - if (!triedBoundCrsToGeogCRSSameAsHubCRS && hubSrcGeog && - geogCRSOfBaseOfBoundSrc) { - // This one should go to the above 'Is it: boundCRS to a geogCRS - // that is the same as the hubCRS ?' case - auto opsFirst = createOperations(sourceCRS, hubSrc, context); - auto opsLast = createOperations(hubSrc, targetCRS, context); - if (!opsFirst.empty() && !opsLast.empty()) { - for (const auto &opFirst : opsFirst) { - for (const auto &opLast : opsLast) { - // Exclude artificial transformations from the hub - // to the target CRS - if (!opLast->hasBallparkTransformation()) { - try { - res.emplace_back( - ConcatenatedOperation:: - createComputeMetadata( - {opFirst, opLast}, - !allowEmptyIntersection)); - } catch ( - const InvalidOperationEmptyIntersection &) { - } - } - } - } - if (!res.empty()) { - return res; - } - } - } - - auto vertCRSOfBaseOfBoundSrc = - dynamic_cast(boundSrc->baseCRS().get()); - if (vertCRSOfBaseOfBoundSrc && hubSrcGeog) { - auto opsFirst = createOperations(sourceCRS, hubSrc, context); - if (context.skipHorizontalTransformation) { - if (!opsFirst.empty()) - return opsFirst; - } else { - auto opsSecond = createOperations(hubSrc, targetCRS, context); - if (!opsFirst.empty() && !opsSecond.empty()) { - for (const auto &opFirst : opsFirst) { - for (const auto &opLast : opsSecond) { - // Exclude artificial transformations from the hub - // to the target CRS - if (!opLast->hasBallparkTransformation()) { - try { - res.emplace_back( - ConcatenatedOperation:: - createComputeMetadata( - {opFirst, opLast}, - !allowEmptyIntersection)); - } catch ( - const InvalidOperationEmptyIntersection &) { - } - } - } - } - if (!res.empty()) { - return res; - } - } - } - } - - return createOperations(boundSrc->baseCRS(), targetCRS, context); + createOperationsBoundToGeog(sourceCRS, targetCRS, context, boundSrc, + geogDst, res); + return res; } // reverse of previous case @@ -13009,18 +12573,9 @@ CoordinateOperationFactory::Private::createOperations( // vertCRS (as boundCRS with transformation to target vertCRS) to // vertCRS if (boundSrc && vertDst) { - auto baseSrcVert = - dynamic_cast(boundSrc->baseCRS().get()); - const auto &hubSrc = boundSrc->hubCRS(); - auto hubSrcVert = dynamic_cast(hubSrc.get()); - if (baseSrcVert && hubSrcVert && - vertDst->_isEquivalentTo( - hubSrcVert, util::IComparable::Criterion::EQUIVALENT)) { - res.emplace_back(boundSrc->transformation()); - return res; - } - - return createOperations(boundSrc->baseCRS(), targetCRS, context); + createOperationsBoundToVert(sourceCRS, targetCRS, context, boundSrc, + vertDst, res); + return res; } // reverse of previous case @@ -13029,160 +12584,509 @@ CoordinateOperationFactory::Private::createOperations( } if (vertSrc && vertDst) { - const auto &srcDatum = vertSrc->datum(); - const auto &dstDatum = vertDst->datum(); - const bool equivalentVDatum = - (srcDatum && dstDatum && - srcDatum->_isEquivalentTo( - dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)); - - const double convSrc = - vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI(); - const double convDst = - vertDst->coordinateSystem()->axisList()[0]->unit().conversionToSI(); - - const double factor = convSrc / convDst; - auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()); - if (!equivalentVDatum) { - name += BALLPARK_VERTICAL_TRANSFORMATION; - auto conv = Transformation::createChangeVerticalUnit( - util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, - name), - sourceCRS, targetCRS, common::Scale(factor), {}); - conv->setHasBallparkTransformation(true); - res.push_back(conv); - } else if (convSrc != convDst) { - auto conv = Conversion::createChangeVerticalUnit( - util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, - name), - common::Scale(factor)); - conv->setCRSs(sourceCRS, targetCRS, nullptr); - res.push_back(conv); - } + createOperationsVertToVert(sourceCRS, targetCRS, context, vertSrc, + vertDst, res); return res; } // A bit odd case as we are comparing apples to oranges, but in case // the vertical unit differ, do something useful. if (vertSrc && geogDst) { + createOperationsVertToGeog(sourceCRS, targetCRS, context, vertSrc, + geogDst, res); + return res; + } - if (vertSrc->identifiers().empty()) { - const auto &vertSrcName = vertSrc->nameStr(); - if (authFactory != nullptr && vertSrcName != "unnamed" && - vertSrcName != "unknown") { - auto matches = authFactory->createObjectsFromName( - vertSrcName, - {io::AuthorityFactory::ObjectType::VERTICAL_CRS}, false, 2); - if (matches.size() == 1) { - const auto &match = matches.front(); - if (vertSrc->_isEquivalentTo( - match.get(), - util::IComparable::Criterion::EQUIVALENT) && - !match->identifiers().empty()) { - return createOperations( - NN_NO_CHECK( - util::nn_dynamic_pointer_cast( - match)), - targetCRS, context); + // reverse of previous case + if (vertDst && geogSrc) { + return applyInverse(createOperations(targetCRS, sourceCRS, context)); + } + + // boundCRS to boundCRS using the same geographic hubCRS + if (boundSrc && boundDst) { + createOperationsBoundToBound(sourceCRS, targetCRS, context, boundSrc, + boundDst, res); + return res; + } + + auto compoundSrc = dynamic_cast(sourceCRS.get()); + // Order of comparison between the geogDst vs geodDst is impotant + if (compoundSrc && geogDst) { + createOperationsCompoundToGeog(sourceCRS, targetCRS, context, + compoundSrc, geogDst, res); + return res; + } else if (compoundSrc && geodDst) { + createOperationsCompoundToGeod(sourceCRS, targetCRS, context, + compoundSrc, geodDst, res); + return res; + } + + // reverse of previous cases + auto compoundDst = dynamic_cast(targetCRS.get()); + if (geodSrc && compoundDst) { + return applyInverse(createOperations(targetCRS, sourceCRS, context)); + } + + if (compoundSrc && compoundDst) { + createOperationsCompoundToCompound(sourceCRS, targetCRS, context, + compoundSrc, compoundDst, res); + return res; + } + + // '+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +type=crs' to + // '+proj=longlat +ellps=GRS80 +nadgrids=@bar.gsb +geoidgrids=@bar.gtx + // +type=crs' + if (boundSrc && compoundDst) { + createOperationsBoundToCompound(sourceCRS, targetCRS, context, boundSrc, + compoundDst, res); + return res; + } + + // reverse of previous case + if (boundDst && compoundSrc) { + return applyInverse(createOperations(targetCRS, sourceCRS, context)); + } + + return res; +} + +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private::createOperationsFromProj4Ext( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + const crs::BoundCRS *boundSrc, const crs::BoundCRS *boundDst, + std::vector &res) { + + auto sourceProjExportable = dynamic_cast( + boundSrc ? boundSrc : sourceCRS.get()); + auto targetProjExportable = dynamic_cast( + boundDst ? boundDst : targetCRS.get()); + if (!sourceProjExportable) { + throw InvalidOperation("Source CRS is not PROJ exportable"); + } + if (!targetProjExportable) { + throw InvalidOperation("Target CRS is not PROJ exportable"); + } + auto projFormatter = io::PROJStringFormatter::create(); + projFormatter->setCRSExport(true); + projFormatter->setLegacyCRSToCRSContext(true); + projFormatter->startInversion(); + sourceProjExportable->_exportToPROJString(projFormatter.get()); + auto geogSrc = dynamic_cast(sourceCRS.get()); + if (geogSrc) { + auto tmpFormatter = io::PROJStringFormatter::create(); + geogSrc->addAngularUnitConvertAndAxisSwap(tmpFormatter.get()); + projFormatter->ingestPROJString(tmpFormatter->toString()); + } + + projFormatter->stopInversion(); + + targetProjExportable->_exportToPROJString(projFormatter.get()); + auto geogDst = dynamic_cast(targetCRS.get()); + if (geogDst) { + auto tmpFormatter = io::PROJStringFormatter::create(); + geogDst->addAngularUnitConvertAndAxisSwap(tmpFormatter.get()); + projFormatter->ingestPROJString(tmpFormatter->toString()); + } + + const auto PROJString = projFormatter->toString(); + auto properties = util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr())); + res.emplace_back(SingleOperation::createPROJBased( + properties, PROJString, sourceCRS, targetCRS, {})); +} + +// --------------------------------------------------------------------------- + +bool CoordinateOperationFactory::Private::createOperationsFromDatabase( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::GeodeticCRS *geodSrc, + const crs::GeodeticCRS *geodDst, const crs::GeographicCRS *geogSrc, + const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc, + const crs::VerticalCRS *vertDst, + std::vector &res) { + + bool doFilterAndCheckPerfectOp = true; + res = findOpsInRegistryDirect(sourceCRS, targetCRS, context.context); + if (!sourceCRS->_isEquivalentTo(targetCRS.get())) { + auto resFromInverse = applyInverse( + findOpsInRegistryDirect(targetCRS, sourceCRS, context.context)); + res.insert(res.end(), resFromInverse.begin(), resFromInverse.end()); + + // If we get at least a result with perfect accuracy, do not + // bother generating synthetic transforms. + if (hasPerfectAccuracyResult(res, context)) { + return true; + } + + doFilterAndCheckPerfectOp = false; + + bool sameGeodeticDatum = false; + 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( + 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) { + createOperationsWithDatumPivot(res, sourceCRS, targetCRS, + geodSrc, geodDst, context); + doFilterAndCheckPerfectOp = !res.empty(); } } - const double convSrc = - vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI(); - double convDst = 1.0; - const auto &geogAxis = geogDst->coordinateSystem()->axisList(); - if (geogAxis.size() == 3) { - convDst = geogAxis[2]->unit().conversionToSI(); + // NAD27 to NAD83 has tens of results already. No need to look + // for a pivot + if (!sameGeodeticDatum && + ((res.empty() && + context.context->getAllowUseIntermediateCRS() == + CoordinateOperationContext::IntermediateCRSUse:: + IF_NO_DIRECT_TRANSFORMATION) || + context.context->getAllowUseIntermediateCRS() == + CoordinateOperationContext::IntermediateCRSUse::ALWAYS || + getenv("PROJ_FORCE_SEARCH_PIVOT"))) { + auto resWithIntermediate = findsOpsInRegistryWithIntermediate( + sourceCRS, targetCRS, context.context); + res.insert(res.end(), resWithIntermediate.begin(), + resWithIntermediate.end()); + doFilterAndCheckPerfectOp = !res.empty(); + + // If transforming from/to WGS84 (Gxxxx), try through 'neutral' + // WGS84 + if (res.empty() && geodSrc && geodDst && + !context.inCreateOperationsThroughPreferredHub && + !context.inCreateOperationsWithDatumPivotAntiRecursion) { + createOperationsThroughPreferredHub(res, sourceCRS, targetCRS, + geodSrc, geodDst, context); + doFilterAndCheckPerfectOp = !res.empty(); + } } + } - const double factor = convSrc / convDst; - auto conv = Transformation::createChangeVerticalUnit( - util::PropertyMap().set( - common::IdentifiedObject::NAME_KEY, - buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) + - BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT), - sourceCRS, targetCRS, common::Scale(factor), {}); - conv->setHasBallparkTransformation(true); - res.push_back(conv); - return res; + if (doFilterAndCheckPerfectOp) { + // If we get at least a result with perfect accuracy, do not bother + // generating synthetic transforms. + if (hasPerfectAccuracyResult(res, context)) { + return true; + } } + return false; +} - // reverse of previous case - if (vertDst && geogSrc) { - return applyInverse(createOperations(targetCRS, sourceCRS, context)); +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private::createOperationsGeodToGeod( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::GeodeticCRS *geodSrc, + const crs::GeodeticCRS *geodDst, + std::vector &res) { + + if (geodSrc->ellipsoid()->celestialBody() != + geodDst->ellipsoid()->celestialBody()) { + throw util::UnsupportedOperationException( + "Source and target ellipsoid do not belong to the same " + "celestial body"); } - // boundCRS to boundCRS using the same geographic hubCRS - if (boundSrc && boundDst) { - const auto &hubSrc = boundSrc->hubCRS(); - auto hubSrcGeog = - dynamic_cast(hubSrc.get()); - const auto &hubDst = boundDst->hubCRS(); - auto hubDstGeog = - dynamic_cast(hubDst.get()); - auto geogCRSOfBaseOfBoundSrc = - boundSrc->baseCRS()->extractGeographicCRS(); - auto geogCRSOfBaseOfBoundDst = - boundDst->baseCRS()->extractGeographicCRS(); - if (hubSrcGeog && hubDstGeog && - hubSrcGeog->_isEquivalentTo( - hubDstGeog, util::IComparable::Criterion::EQUIVALENT) && - geogCRSOfBaseOfBoundSrc && geogCRSOfBaseOfBoundDst) { - const bool firstIsNoOp = geogCRSOfBaseOfBoundSrc->_isEquivalentTo( - boundSrc->baseCRS().get(), - util::IComparable::Criterion::EQUIVALENT); - const bool lastIsNoOp = geogCRSOfBaseOfBoundDst->_isEquivalentTo( - boundDst->baseCRS().get(), - util::IComparable::Criterion::EQUIVALENT); - auto opsFirst = - createOperations(boundSrc->baseCRS(), - NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context); - auto opsLast = - createOperations(NN_NO_CHECK(geogCRSOfBaseOfBoundDst), - boundDst->baseCRS(), context); - if (!opsFirst.empty() && !opsLast.empty()) { - const auto &opSecond = boundSrc->transformation(); - auto opThird = boundDst->transformation()->inverse(); + auto geogSrc = dynamic_cast(geodSrc); + auto geogDst = dynamic_cast(geodDst); + + if (geogSrc && geogDst) { + createOperationsGeogToGeog(res, sourceCRS, targetCRS, geogSrc, geogDst); + return; + } + + const bool isSrcGeocentric = geodSrc->isGeocentric(); + const bool isSrcGeographic = geogSrc != nullptr; + const bool isTargetGeocentric = geodDst->isGeocentric(); + const bool isTargetGeographic = geogDst != nullptr; + if (((isSrcGeocentric && isTargetGeographic) || + (isSrcGeographic && isTargetGeocentric)) && + geodSrc->datum() != nullptr && geodDst->datum() != nullptr) { + + // Same datum ? + if (geodSrc->datum()->_isEquivalentTo( + geodDst->datum().get(), + util::IComparable::Criterion::EQUIVALENT)) { + res.emplace_back( + Conversion::createGeographicGeocentric(sourceCRS, targetCRS)); + } else if (isSrcGeocentric) { + std::string interm_crs_name(geogDst->nameStr()); + interm_crs_name += " (geocentric)"; + auto interm_crs = + util::nn_static_pointer_cast(crs::GeodeticCRS::create( + addDomains(util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + interm_crs_name), + geogDst), + NN_NO_CHECK(geogDst->datum()), + NN_CHECK_ASSERT( + util::nn_dynamic_pointer_cast( + geodSrc->coordinateSystem())))); + auto opFirst = + createBallparkGeocentricTranslation(sourceCRS, interm_crs); + auto opSecond = + Conversion::createGeographicGeocentric(interm_crs, targetCRS); + res.emplace_back(ConcatenatedOperation::createComputeMetadata( + {opFirst, opSecond}, !allowEmptyIntersection)); + } else { + res = applyInverse(createOperations(targetCRS, sourceCRS, context)); + } + + return; + } + + if (isSrcGeocentric && isTargetGeocentric) { + res.emplace_back( + createBallparkGeocentricTranslation(sourceCRS, targetCRS)); + return; + } + + // Transformation between two geodetic systems of unknown type + // This should normally not be triggered with "standard" CRS + res.emplace_back(createGeodToGeodPROJBased(sourceCRS, targetCRS)); +} + +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private::createOperationsDerivedTo( + const crs::CRSNNPtr & /*sourceCRS*/, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::DerivedCRS *derivedSrc, + std::vector &res) { + + auto opFirst = derivedSrc->derivingConversion()->inverse(); + // Small optimization if the targetCRS is the baseCRS of the source + // derivedCRS. + if (derivedSrc->baseCRS()->_isEquivalentTo( + targetCRS.get(), util::IComparable::Criterion::EQUIVALENT)) { + res.emplace_back(opFirst); + return; + } + auto opsSecond = + createOperations(derivedSrc->baseCRS(), targetCRS, context); + for (const auto &opSecond : opsSecond) { + try { + res.emplace_back(ConcatenatedOperation::createComputeMetadata( + {opFirst, opSecond}, !allowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } + } +} + +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private::createOperationsBoundToGeog( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::BoundCRS *boundSrc, + const crs::GeographicCRS *geogDst, + std::vector &res) { + + const auto &hubSrc = boundSrc->hubCRS(); + auto hubSrcGeog = dynamic_cast(hubSrc.get()); + auto geogCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractGeographicCRS(); + bool triedBoundCrsToGeogCRSSameAsHubCRS = false; + // Is it: boundCRS to a geogCRS that is the same as the hubCRS ? + if (hubSrcGeog && geogCRSOfBaseOfBoundSrc && + (hubSrcGeog->_isEquivalentTo( + geogDst, util::IComparable::Criterion::EQUIVALENT) || + hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst)))) { + triedBoundCrsToGeogCRSSameAsHubCRS = true; + if (boundSrc->baseCRS() == geogCRSOfBaseOfBoundSrc) { + // Optimization to avoid creating a useless concatenated + // operation + res.emplace_back(boundSrc->transformation()); + return; + } + auto opsFirst = createOperations( + boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context); + if (!opsFirst.empty()) { + for (const auto &opFirst : opsFirst) { + try { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + {opFirst, boundSrc->transformation()}, + !allowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } + } + if (!res.empty()) { + return; + } + } + // If the datum are equivalent, this is also fine + } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() && + geogDst->datum() && + hubSrcGeog->datum()->_isEquivalentTo( + geogDst->datum().get(), + util::IComparable::Criterion::EQUIVALENT)) { + auto opsFirst = createOperations( + boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context); + auto opsLast = createOperations(hubSrc, targetCRS, context); + if (!opsFirst.empty() && !opsLast.empty()) { + for (const auto &opFirst : opsFirst) { + for (const auto &opLast : opsLast) { + try { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + {opFirst, boundSrc->transformation(), opLast}, + !allowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } + } + } + if (!res.empty()) { + return; + } + } + // Consider WGS 84 and NAD83 as equivalent in that context if the + // geogCRSOfBaseOfBoundSrc ellipsoid is Clarke66 (for NAD27) + // Case of "+proj=latlong +ellps=clrk66 + // +nadgrids=ntv1_can.dat,conus" + // to "+proj=latlong +datum=NAD83" + } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() && + geogDst->datum() && + geogCRSOfBaseOfBoundSrc->ellipsoid()->_isEquivalentTo( + datum::Ellipsoid::CLARKE_1866.get(), + util::IComparable::Criterion::EQUIVALENT) && + hubSrcGeog->datum()->_isEquivalentTo( + datum::GeodeticReferenceFrame::EPSG_6326.get(), + util::IComparable::Criterion::EQUIVALENT) && + geogDst->datum()->_isEquivalentTo( + datum::GeodeticReferenceFrame::EPSG_6269.get(), + util::IComparable::Criterion::EQUIVALENT)) { + auto nnGeogCRSOfBaseOfBoundSrc = NN_NO_CHECK(geogCRSOfBaseOfBoundSrc); + if (boundSrc->baseCRS()->_isEquivalentTo( + nnGeogCRSOfBaseOfBoundSrc.get(), + util::IComparable::Criterion::EQUIVALENT)) { + auto transf = boundSrc->transformation()->shallowClone(); + transf->setProperties(util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + buildTransfName(boundSrc->baseCRS()->nameStr(), + targetCRS->nameStr()))); + transf->setCRSs(boundSrc->baseCRS(), targetCRS, nullptr); + res.emplace_back(transf); + return; + } else { + auto opsFirst = createOperations( + boundSrc->baseCRS(), nnGeogCRSOfBaseOfBoundSrc, context); + auto transf = boundSrc->transformation()->shallowClone(); + transf->setProperties(util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + buildTransfName(nnGeogCRSOfBaseOfBoundSrc->nameStr(), + targetCRS->nameStr()))); + transf->setCRSs(nnGeogCRSOfBaseOfBoundSrc, targetCRS, nullptr); + if (!opsFirst.empty()) { for (const auto &opFirst : opsFirst) { - for (const auto &opLast : opsLast) { - try { - std::vector ops; - if (!firstIsNoOp) { - ops.push_back(opFirst); - } - ops.push_back(opSecond); - ops.push_back(opThird); - if (!lastIsNoOp) { - ops.push_back(opLast); - } - res.emplace_back( - ConcatenatedOperation::createComputeMetadata( - ops, !allowEmptyIntersection)); - } catch (const InvalidOperationEmptyIntersection &) { - } + try { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + {opFirst, transf}, !allowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { } } if (!res.empty()) { - return res; + return; } } } + } - auto vertCRSOfBaseOfBoundSrc = - boundSrc->baseCRS()->extractVerticalCRS(); - auto vertCRSOfBaseOfBoundDst = - boundDst->baseCRS()->extractVerticalCRS(); - if (hubSrcGeog && hubDstGeog && - hubSrcGeog->_isEquivalentTo( - hubDstGeog, util::IComparable::Criterion::EQUIVALENT) && - vertCRSOfBaseOfBoundSrc && vertCRSOfBaseOfBoundDst) { - auto opsFirst = createOperations(sourceCRS, hubSrc, context); - auto opsLast = createOperations(hubSrc, targetCRS, context); - if (!opsFirst.empty() && !opsLast.empty()) { - for (const auto &opFirst : opsFirst) { - for (const auto &opLast : opsLast) { + if (hubSrcGeog && + hubSrcGeog->_isEquivalentTo(geogDst, + util::IComparable::Criterion::EQUIVALENT) && + dynamic_cast(boundSrc->baseCRS().get())) { + res.emplace_back(boundSrc->transformation()); + return; + } + + if (!triedBoundCrsToGeogCRSSameAsHubCRS && hubSrcGeog && + geogCRSOfBaseOfBoundSrc) { + // This one should go to the above 'Is it: boundCRS to a geogCRS + // that is the same as the hubCRS ?' case + auto opsFirst = createOperations(sourceCRS, hubSrc, context); + auto opsLast = createOperations(hubSrc, targetCRS, context); + if (!opsFirst.empty() && !opsLast.empty()) { + for (const auto &opFirst : opsFirst) { + for (const auto &opLast : opsLast) { + // Exclude artificial transformations from the hub + // to the target CRS + if (!opLast->hasBallparkTransformation()) { try { res.emplace_back( ConcatenatedOperation::createComputeMetadata( @@ -13192,58 +13096,379 @@ CoordinateOperationFactory::Private::createOperations( } } } + } + if (!res.empty()) { + return; + } + } + } + + auto vertCRSOfBaseOfBoundSrc = + dynamic_cast(boundSrc->baseCRS().get()); + if (vertCRSOfBaseOfBoundSrc && hubSrcGeog) { + auto opsFirst = createOperations(sourceCRS, hubSrc, context); + if (context.skipHorizontalTransformation) { + if (!opsFirst.empty()) + res = opsFirst; + return; + } else { + auto opsSecond = createOperations(hubSrc, targetCRS, context); + if (!opsFirst.empty() && !opsSecond.empty()) { + for (const auto &opFirst : opsFirst) { + for (const auto &opLast : opsSecond) { + // Exclude artificial transformations from the hub + // to the target CRS + if (!opLast->hasBallparkTransformation()) { + try { + res.emplace_back( + ConcatenatedOperation:: + createComputeMetadata( + {opFirst, opLast}, + !allowEmptyIntersection)); + } catch ( + const InvalidOperationEmptyIntersection &) { + } + } + } + } if (!res.empty()) { - return res; + return; } } } + } + + res = createOperations(boundSrc->baseCRS(), targetCRS, context); +} - return createOperations(boundSrc->baseCRS(), boundDst->baseCRS(), - context); +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private::createOperationsBoundToVert( + const crs::CRSNNPtr & /*sourceCRS*/, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::BoundCRS *boundSrc, + const crs::VerticalCRS *vertDst, + std::vector &res) { + + auto baseSrcVert = + dynamic_cast(boundSrc->baseCRS().get()); + const auto &hubSrc = boundSrc->hubCRS(); + auto hubSrcVert = dynamic_cast(hubSrc.get()); + if (baseSrcVert && hubSrcVert && + vertDst->_isEquivalentTo(hubSrcVert, + util::IComparable::Criterion::EQUIVALENT)) { + res.emplace_back(boundSrc->transformation()); + return; } - auto compoundSrc = dynamic_cast(sourceCRS.get()); - if (compoundSrc && geogDst) { - const auto &componentsSrc = compoundSrc->componentReferenceSystems(); - if (!componentsSrc.empty()) { - std::vector horizTransforms; - auto srcGeogCRS = componentsSrc[0]->extractGeographicCRS(); - if (srcGeogCRS) { - horizTransforms = - createOperations(componentsSrc[0], targetCRS, context); + res = createOperations(boundSrc->baseCRS(), targetCRS, context); +} + +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private::createOperationsVertToVert( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context & /*context*/, const crs::VerticalCRS *vertSrc, + const crs::VerticalCRS *vertDst, + std::vector &res) { + + const auto &srcDatum = vertSrc->datum(); + const auto &dstDatum = vertDst->datum(); + const bool equivalentVDatum = + (srcDatum && dstDatum && + srcDatum->_isEquivalentTo(dstDatum.get(), + util::IComparable::Criterion::EQUIVALENT)); + + const double convSrc = + vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI(); + const double convDst = + vertDst->coordinateSystem()->axisList()[0]->unit().conversionToSI(); + + const double factor = convSrc / convDst; + auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()); + if (!equivalentVDatum) { + name += BALLPARK_VERTICAL_TRANSFORMATION; + auto conv = Transformation::createChangeVerticalUnit( + util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name), + sourceCRS, targetCRS, common::Scale(factor), {}); + conv->setHasBallparkTransformation(true); + res.push_back(conv); + } else if (convSrc != convDst) { + auto conv = Conversion::createChangeVerticalUnit( + util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name), + common::Scale(factor)); + conv->setCRSs(sourceCRS, targetCRS, nullptr); + res.push_back(conv); + } +} + +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private::createOperationsVertToGeog( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::VerticalCRS *vertSrc, + const crs::GeographicCRS *geogDst, + std::vector &res) { + + if (vertSrc->identifiers().empty()) { + const auto &vertSrcName = vertSrc->nameStr(); + const auto &authFactory = context.context->getAuthorityFactory(); + if (authFactory != nullptr && vertSrcName != "unnamed" && + vertSrcName != "unknown") { + auto matches = authFactory->createObjectsFromName( + vertSrcName, {io::AuthorityFactory::ObjectType::VERTICAL_CRS}, + false, 2); + if (matches.size() == 1) { + const auto &match = matches.front(); + if (vertSrc->_isEquivalentTo( + match.get(), + util::IComparable::Criterion::EQUIVALENT) && + !match->identifiers().empty()) { + res = createOperations( + NN_NO_CHECK( + util::nn_dynamic_pointer_cast( + match)), + targetCRS, context); + return; + } } - std::vector verticalTransforms; + } + } - const auto dbContext = - authFactory ? authFactory->databaseContext().as_nullable() - : nullptr; - if (componentsSrc.size() >= 2 && - componentsSrc[1]->extractVerticalCRS()) { + const double convSrc = + vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI(); + double convDst = 1.0; + const auto &geogAxis = geogDst->coordinateSystem()->axisList(); + if (geogAxis.size() == 3) { + convDst = geogAxis[2]->unit().conversionToSI(); + } + + const double factor = convSrc / convDst; + auto conv = Transformation::createChangeVerticalUnit( + util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) + + BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT), + sourceCRS, targetCRS, common::Scale(factor), {}); + conv->setHasBallparkTransformation(true); + res.push_back(conv); +} - struct SetSkipHorizontalTransform { - Context &context; +// --------------------------------------------------------------------------- - explicit SetSkipHorizontalTransform(Context &contextIn) - : context(contextIn) { - assert(!context.skipHorizontalTransformation); - context.skipHorizontalTransformation = true; +void CoordinateOperationFactory::Private::createOperationsBoundToBound( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::BoundCRS *boundSrc, + const crs::BoundCRS *boundDst, std::vector &res) { + + const auto &hubSrc = boundSrc->hubCRS(); + auto hubSrcGeog = dynamic_cast(hubSrc.get()); + const auto &hubDst = boundDst->hubCRS(); + auto hubDstGeog = dynamic_cast(hubDst.get()); + auto geogCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractGeographicCRS(); + auto geogCRSOfBaseOfBoundDst = boundDst->baseCRS()->extractGeographicCRS(); + if (hubSrcGeog && hubDstGeog && + hubSrcGeog->_isEquivalentTo(hubDstGeog, + util::IComparable::Criterion::EQUIVALENT) && + geogCRSOfBaseOfBoundSrc && geogCRSOfBaseOfBoundDst) { + const bool firstIsNoOp = geogCRSOfBaseOfBoundSrc->_isEquivalentTo( + boundSrc->baseCRS().get(), + util::IComparable::Criterion::EQUIVALENT); + const bool lastIsNoOp = geogCRSOfBaseOfBoundDst->_isEquivalentTo( + boundDst->baseCRS().get(), + util::IComparable::Criterion::EQUIVALENT); + auto opsFirst = createOperations( + boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context); + auto opsLast = createOperations(NN_NO_CHECK(geogCRSOfBaseOfBoundDst), + boundDst->baseCRS(), context); + if (!opsFirst.empty() && !opsLast.empty()) { + const auto &opSecond = boundSrc->transformation(); + auto opThird = boundDst->transformation()->inverse(); + for (const auto &opFirst : opsFirst) { + for (const auto &opLast : opsLast) { + try { + std::vector ops; + if (!firstIsNoOp) { + ops.push_back(opFirst); + } + ops.push_back(opSecond); + ops.push_back(opThird); + if (!lastIsNoOp) { + ops.push_back(opLast); + } + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + ops, !allowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { } + } + } + if (!res.empty()) { + return; + } + } + } - ~SetSkipHorizontalTransform() { - context.skipHorizontalTransformation = false; + auto vertCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractVerticalCRS(); + auto vertCRSOfBaseOfBoundDst = boundDst->baseCRS()->extractVerticalCRS(); + if (hubSrcGeog && hubDstGeog && + hubSrcGeog->_isEquivalentTo(hubDstGeog, + util::IComparable::Criterion::EQUIVALENT) && + vertCRSOfBaseOfBoundSrc && vertCRSOfBaseOfBoundDst) { + auto opsFirst = createOperations(sourceCRS, hubSrc, context); + auto opsLast = createOperations(hubSrc, targetCRS, context); + if (!opsFirst.empty() && !opsLast.empty()) { + for (const auto &opFirst : opsFirst) { + for (const auto &opLast : opsLast) { + try { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + {opFirst, opLast}, !allowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { } - }; - SetSkipHorizontalTransform setSkipHorizontalTransform(context); + } + } + if (!res.empty()) { + return; + } + } + } - verticalTransforms = createOperations( + res = createOperations(boundSrc->baseCRS(), boundDst->baseCRS(), context); +} + +// --------------------------------------------------------------------------- + +static std::vector +getOps(const CoordinateOperationNNPtr &op) { + auto concatenated = dynamic_cast(op.get()); + if (concatenated) + return concatenated->operations(); + return {op}; +} + +static bool useDifferentTransformationsForSameSourceTarget( + const CoordinateOperationNNPtr &opA, const CoordinateOperationNNPtr &opB) { + auto subOpsA = getOps(opA); + auto subOpsB = getOps(opB); + for (const auto &subOpA : subOpsA) { + if (!dynamic_cast(subOpA.get())) + continue; + if (subOpA->sourceCRS()->nameStr() == "unknown" || + subOpA->targetCRS()->nameStr() == "unknown") + continue; + for (const auto &subOpB : subOpsB) { + if (!dynamic_cast(subOpB.get())) + continue; + if (subOpB->sourceCRS()->nameStr() == "unknown" || + subOpB->targetCRS()->nameStr() == "unknown") + continue; + + if (subOpA->sourceCRS()->nameStr() == + subOpB->sourceCRS()->nameStr() && + subOpA->targetCRS()->nameStr() == + subOpB->targetCRS()->nameStr()) { + if (starts_with(subOpA->nameStr(), NULL_GEOGRAPHIC_OFFSET) && + starts_with(subOpB->nameStr(), NULL_GEOGRAPHIC_OFFSET)) { + continue; + } + + if (!subOpA->isEquivalentTo(subOpB.get())) { + return true; + } + } else if (subOpA->sourceCRS()->nameStr() == + subOpB->targetCRS()->nameStr() && + subOpA->targetCRS()->nameStr() == + subOpB->sourceCRS()->nameStr()) { + if (starts_with(subOpA->nameStr(), NULL_GEOGRAPHIC_OFFSET) && + starts_with(subOpB->nameStr(), NULL_GEOGRAPHIC_OFFSET)) { + continue; + } + + if (!subOpA->isEquivalentTo(subOpB->inverse().get())) { + return true; + } + } + } + } + return false; +} + +void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::CompoundCRS *compoundSrc, + const crs::GeographicCRS *geogDst, + std::vector &res) { + + const auto &authFactory = context.context->getAuthorityFactory(); + const auto &componentsSrc = compoundSrc->componentReferenceSystems(); + if (!componentsSrc.empty()) { + std::vector horizTransforms; + auto srcGeogCRS = componentsSrc[0]->extractGeographicCRS(); + if (srcGeogCRS) { + horizTransforms = + createOperations(componentsSrc[0], targetCRS, context); + } + std::vector verticalTransforms; + + const auto dbContext = + authFactory ? authFactory->databaseContext().as_nullable() + : nullptr; + if (componentsSrc.size() >= 2 && + componentsSrc[1]->extractVerticalCRS()) { + + struct SetSkipHorizontalTransform { + Context &context; + + explicit SetSkipHorizontalTransform(Context &contextIn) + : context(contextIn) { + assert(!context.skipHorizontalTransformation); + context.skipHorizontalTransformation = true; + } + + ~SetSkipHorizontalTransform() { + context.skipHorizontalTransformation = false; + } + }; + SetSkipHorizontalTransform setSkipHorizontalTransform(context); + + verticalTransforms = createOperations( + componentsSrc[1], + targetCRS->promoteTo3D(std::string(), dbContext), context); + bool foundRegisteredTransformWithAllGridsAvailable = false; + const bool ignoreMissingGrids = + context.context->getGridAvailabilityUse() == + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY; + for (const auto &op : verticalTransforms) { + if (!op->identifiers().empty() && dbContext) { + bool missingGrid = false; + if (!ignoreMissingGrids) { + const auto gridsNeeded = op->gridsNeeded(dbContext); + for (const auto &gridDesc : gridsNeeded) { + if (!gridDesc.available) { + missingGrid = true; + break; + } + } + } + if (!missingGrid) { + foundRegisteredTransformWithAllGridsAvailable = true; + break; + } + } + } + if (!foundRegisteredTransformWithAllGridsAvailable && srcGeogCRS && + !srcGeogCRS->_isEquivalentTo( + geogDst, util::IComparable::Criterion::EQUIVALENT) && + !srcGeogCRS->is2DPartOf3D(NN_NO_CHECK(geogDst))) { + auto verticalTransformsTmp = createOperations( componentsSrc[1], - targetCRS->promoteTo3D(std::string(), dbContext), context); - bool foundRegisteredTransformWithAllGridsAvailable = false; - const bool ignoreMissingGrids = - context.context->getGridAvailabilityUse() == - CoordinateOperationContext::GridAvailabilityUse:: - IGNORE_GRID_AVAILABILITY; - for (const auto &op : verticalTransforms) { + NN_NO_CHECK(srcGeogCRS) + ->promoteTo3D(std::string(), dbContext), + context); + bool foundRegisteredTransform = false; + foundRegisteredTransformWithAllGridsAvailable = false; + for (const auto &op : verticalTransformsTmp) { if (!op->identifiers().empty() && dbContext) { bool missingGrid = false; if (!ignoreMissingGrids) { @@ -13255,6 +13480,7 @@ CoordinateOperationFactory::Private::createOperations( } } } + foundRegisteredTransform = true; if (!missingGrid) { foundRegisteredTransformWithAllGridsAvailable = true; @@ -13262,127 +13488,77 @@ CoordinateOperationFactory::Private::createOperations( } } } - if (!foundRegisteredTransformWithAllGridsAvailable && - srcGeogCRS && - !srcGeogCRS->_isEquivalentTo( - geogDst, util::IComparable::Criterion::EQUIVALENT) && - !srcGeogCRS->is2DPartOf3D(NN_NO_CHECK(geogDst))) { - auto verticalTransformsTmp = createOperations( - componentsSrc[1], - NN_NO_CHECK(srcGeogCRS) - ->promoteTo3D(std::string(), dbContext), - context); - bool foundRegisteredTransform = false; - foundRegisteredTransformWithAllGridsAvailable = false; - for (const auto &op : verticalTransformsTmp) { - if (!op->identifiers().empty() && dbContext) { - bool missingGrid = false; - if (!ignoreMissingGrids) { - const auto gridsNeeded = - op->gridsNeeded(dbContext); - for (const auto &gridDesc : gridsNeeded) { - if (!gridDesc.available) { - missingGrid = true; - break; - } - } - } - foundRegisteredTransform = true; - if (!missingGrid) { - foundRegisteredTransformWithAllGridsAvailable = - true; - break; - } - } - } - if (foundRegisteredTransformWithAllGridsAvailable) { - verticalTransforms = verticalTransformsTmp; - } else if (foundRegisteredTransform) { - verticalTransforms.insert(verticalTransforms.end(), - verticalTransformsTmp.begin(), - verticalTransformsTmp.end()); - } + if (foundRegisteredTransformWithAllGridsAvailable) { + verticalTransforms = verticalTransformsTmp; + } else if (foundRegisteredTransform) { + verticalTransforms.insert(verticalTransforms.end(), + verticalTransformsTmp.begin(), + verticalTransformsTmp.end()); } } + } - if (horizTransforms.empty() || verticalTransforms.empty()) { - return horizTransforms; - } - - typedef std::pair, - std::vector> - PairOfTransforms; - std::map - cacheHorizToInterpAndInterpToTarget; + if (horizTransforms.empty() || verticalTransforms.empty()) { + res = horizTransforms; + return; + } - for (const auto &verticalTransform : verticalTransforms) { - crs::GeographicCRSPtr interpolationGeogCRS; - auto transformationVerticalTransform = - dynamic_cast( - verticalTransform.get()); - if (transformationVerticalTransform && - !transformationVerticalTransform - ->hasBallparkTransformation()) { - auto interpTransformCRS = - transformationVerticalTransform->interpolationCRS(); - if (interpTransformCRS) { - interpolationGeogCRS = - std::dynamic_pointer_cast( - 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( - transformationVerticalTransform->targetCRS() - .as_nullable()); - } + typedef std::pair, + std::vector> + PairOfTransforms; + std::map + cacheHorizToInterpAndInterpToTarget; + + for (const auto &verticalTransform : verticalTransforms) { + crs::GeographicCRSPtr interpolationGeogCRS; + auto transformationVerticalTransform = + dynamic_cast(verticalTransform.get()); + if (transformationVerticalTransform && + !transformationVerticalTransform->hasBallparkTransformation()) { + auto interpTransformCRS = + transformationVerticalTransform->interpolationCRS(); + if (interpTransformCRS) { + interpolationGeogCRS = + std::dynamic_pointer_cast( + 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( + 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(); - } + 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 srcToInterpOps; - std::vector interpToTargetOps; + std::vector srcToInterpOps; + std::vector interpToTargetOps; - std::string key; - const auto &ids = interpolationGeogCRS->identifiers(); - if (!ids.empty()) { - key = (*ids.front()->codeSpace()) + ':' + - ids.front()->code(); - } - if (!key.empty()) { - auto iter = - cacheHorizToInterpAndInterpToTarget.find(key); - if (iter == cacheHorizToInterpAndInterpToTarget.end()) { - srcToInterpOps = createOperations( - componentsSrc[0], - NN_NO_CHECK(interpolationGeogCRS), context); - interpToTargetOps = createOperations( - NN_NO_CHECK(interpolationGeogCRS), - targetCRS->demoteTo2D(std::string(), dbContext), - context); - cacheHorizToInterpAndInterpToTarget[key] = - PairOfTransforms(srcToInterpOps, - interpToTargetOps); - } else { - srcToInterpOps = iter->second.first; - interpToTargetOps = iter->second.second; - } - } else { + std::string key; + const auto &ids = interpolationGeogCRS->identifiers(); + if (!ids.empty()) { + key = + (*ids.front()->codeSpace()) + ':' + ids.front()->code(); + } + if (!key.empty()) { + auto iter = cacheHorizToInterpAndInterpToTarget.find(key); + if (iter == cacheHorizToInterpAndInterpToTarget.end()) { srcToInterpOps = createOperations( componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS), context); @@ -13390,233 +13566,245 @@ CoordinateOperationFactory::Private::createOperations( NN_NO_CHECK(interpolationGeogCRS), targetCRS->demoteTo2D(std::string(), dbContext), context); + cacheHorizToInterpAndInterpToTarget[key] = + PairOfTransforms(srcToInterpOps, interpToTargetOps); + } else { + srcToInterpOps = iter->second.first; + interpToTargetOps = iter->second.second; } + } else { + srcToInterpOps = createOperations( + componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS), + context); + interpToTargetOps = createOperations( + NN_NO_CHECK(interpolationGeogCRS), + targetCRS->demoteTo2D(std::string(), dbContext), + context); + } - for (const auto &srcToInterp : srcToInterpOps) { - for (const auto &interpToTarget : interpToTargetOps) { - - if (useDifferentTransformationsForSameSourceTarget( - srcToInterp, interpToTarget)) { - continue; - } + for (const auto &srcToInterp : srcToInterpOps) { + for (const auto &interpToTarget : interpToTargetOps) { - try { - auto op = createHorizVerticalHorizPROJBased( - sourceCRS, targetCRS, srcToInterp, - verticalTransform, interpToTarget, - interpolationGeogCRS, true); - res.emplace_back(op); - } catch (const std::exception &) { - } + if (useDifferentTransformationsForSameSourceTarget( + srcToInterp, interpToTarget)) { + continue; } - } - } else { - // This case is probably only correct if - // verticalTransform and horizTransform are independant - // and in particular that verticalTransform does not - // involve a grid, because of the rather arbitrary order - // horizontal then vertical applied - for (const auto &horizTransform : horizTransforms) { + try { - auto op = createHorizVerticalPROJBased( - sourceCRS, targetCRS, horizTransform, - verticalTransform); + auto op = createHorizVerticalHorizPROJBased( + sourceCRS, targetCRS, srcToInterp, + verticalTransform, interpToTarget, + interpolationGeogCRS, true); res.emplace_back(op); } catch (const std::exception &) { } } } - } - - return res; - } - } else if (compoundSrc && geodDst) { - auto datum = geodDst->datum(); - if (datum) { - auto cs = - cs::EllipsoidalCS::createLatitudeLongitudeEllipsoidalHeight( - common::UnitOfMeasure::DEGREE, - common::UnitOfMeasure::METRE); - auto intermGeog3DCRS = util::nn_static_pointer_cast( - crs::GeographicCRS::create( - util::PropertyMap() - .set(common::IdentifiedObject::NAME_KEY, - geodDst->nameStr()) - .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY, - metadata::Extent::WORLD), - NN_NO_CHECK(datum), cs)); - auto sourceToGeog3DOps = - createOperations(sourceCRS, intermGeog3DCRS, context); - auto geog3DToTargetOps = - createOperations(intermGeog3DCRS, targetCRS, context); - if (!geog3DToTargetOps.empty()) { - for (const auto &op : sourceToGeog3DOps) { - res.emplace_back( - ConcatenatedOperation::createComputeMetadata( - {op, geog3DToTargetOps.front()}, - !allowEmptyIntersection)); + } else { + // This case is probably only correct if + // verticalTransform and horizTransform are independant + // and in particular that verticalTransform does not + // involve a grid, because of the rather arbitrary order + // horizontal then vertical applied + for (const auto &horizTransform : horizTransforms) { + try { + auto op = createHorizVerticalPROJBased( + sourceCRS, targetCRS, horizTransform, + verticalTransform); + res.emplace_back(op); + } catch (const std::exception &) { + } } - return res; } } } +} - // reverse of previous case - auto compoundDst = dynamic_cast(targetCRS.get()); - if (geodSrc && compoundDst) { - return applyInverse(createOperations(targetCRS, sourceCRS, context)); +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private::createOperationsCompoundToGeod( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::CompoundCRS * /*compoundSrc*/, + const crs::GeodeticCRS *geodDst, + std::vector &res) { + auto datum = geodDst->datum(); + if (datum) { + auto cs = cs::EllipsoidalCS::createLatitudeLongitudeEllipsoidalHeight( + common::UnitOfMeasure::DEGREE, common::UnitOfMeasure::METRE); + auto intermGeog3DCRS = + util::nn_static_pointer_cast(crs::GeographicCRS::create( + util::PropertyMap() + .set(common::IdentifiedObject::NAME_KEY, geodDst->nameStr()) + .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY, + metadata::Extent::WORLD), + NN_NO_CHECK(datum), cs)); + auto sourceToGeog3DOps = + createOperations(sourceCRS, intermGeog3DCRS, context); + auto geog3DToTargetOps = + createOperations(intermGeog3DCRS, targetCRS, context); + if (!geog3DToTargetOps.empty()) { + for (const auto &op : sourceToGeog3DOps) { + res.emplace_back(ConcatenatedOperation::createComputeMetadata( + {op, geog3DToTargetOps.front()}, !allowEmptyIntersection)); + } + } } +} - if (compoundSrc && compoundDst) { - const auto &componentsSrc = compoundSrc->componentReferenceSystems(); - const auto &componentsDst = compoundDst->componentReferenceSystems(); - if (!componentsSrc.empty() && - componentsSrc.size() == componentsDst.size()) { - if (componentsSrc[0]->extractGeographicCRS() && - componentsDst[0]->extractGeographicCRS()) { - - std::vector verticalTransforms; - if (componentsSrc.size() >= 2 && - componentsSrc[1]->extractVerticalCRS() && - componentsDst[1]->extractVerticalCRS()) { - verticalTransforms = createOperations( - componentsSrc[1], componentsDst[1], context); - } +// --------------------------------------------------------------------------- - for (const auto &verticalTransform : verticalTransforms) { - auto interpolationGeogCRS = - NN_NO_CHECK(componentsSrc[0]->extractGeographicCRS()); - auto transformationVerticalTransform = - dynamic_cast( - verticalTransform.get()); - if (transformationVerticalTransform) { - auto interpTransformCRS = - transformationVerticalTransform->interpolationCRS(); - if (interpTransformCRS) { - auto nn_interpTransformCRS = - NN_NO_CHECK(interpTransformCRS); - if (dynamic_cast( - nn_interpTransformCRS.get())) { - interpolationGeogCRS = - NN_NO_CHECK(util::nn_dynamic_pointer_cast< - crs::GeographicCRS>( - nn_interpTransformCRS)); - } - } - } else { - auto compSrc0BoundCrs = dynamic_cast( - componentsSrc[0].get()); - auto compDst0BoundCrs = dynamic_cast( - componentsDst[0].get()); - if (compSrc0BoundCrs && compDst0BoundCrs && - dynamic_cast( - compSrc0BoundCrs->hubCRS().get()) && - compSrc0BoundCrs->hubCRS()->_isEquivalentTo( - compDst0BoundCrs->hubCRS().get())) { - interpolationGeogCRS = - NN_NO_CHECK(util::nn_dynamic_pointer_cast< - crs::GeographicCRS>( - compSrc0BoundCrs->hubCRS())); +void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::CompoundCRS *compoundSrc, + const crs::CompoundCRS *compoundDst, + std::vector &res) { + + const auto &componentsSrc = compoundSrc->componentReferenceSystems(); + const auto &componentsDst = compoundDst->componentReferenceSystems(); + if (!componentsSrc.empty() && + componentsSrc.size() == componentsDst.size()) { + if (componentsSrc[0]->extractGeographicCRS() && + componentsDst[0]->extractGeographicCRS()) { + + std::vector verticalTransforms; + if (componentsSrc.size() >= 2 && + componentsSrc[1]->extractVerticalCRS() && + componentsDst[1]->extractVerticalCRS()) { + verticalTransforms = createOperations( + componentsSrc[1], componentsDst[1], context); + } + + for (const auto &verticalTransform : verticalTransforms) { + auto interpolationGeogCRS = + NN_NO_CHECK(componentsSrc[0]->extractGeographicCRS()); + auto transformationVerticalTransform = + dynamic_cast( + verticalTransform.get()); + if (transformationVerticalTransform) { + auto interpTransformCRS = + transformationVerticalTransform->interpolationCRS(); + if (interpTransformCRS) { + auto nn_interpTransformCRS = + NN_NO_CHECK(interpTransformCRS); + if (dynamic_cast( + nn_interpTransformCRS.get())) { + interpolationGeogCRS = NN_NO_CHECK( + util::nn_dynamic_pointer_cast< + crs::GeographicCRS>(nn_interpTransformCRS)); } } - auto opSrcCRSToGeogCRS = createOperations( - componentsSrc[0], interpolationGeogCRS, context); - auto opGeogCRStoDstCRS = createOperations( - interpolationGeogCRS, componentsDst[0], context); - for (const auto &opSrc : opSrcCRSToGeogCRS) { - for (const auto &opDst : opGeogCRStoDstCRS) { + } else { + auto compSrc0BoundCrs = + dynamic_cast(componentsSrc[0].get()); + auto compDst0BoundCrs = + dynamic_cast(componentsDst[0].get()); + if (compSrc0BoundCrs && compDst0BoundCrs && + dynamic_cast( + compSrc0BoundCrs->hubCRS().get()) && + compSrc0BoundCrs->hubCRS()->_isEquivalentTo( + compDst0BoundCrs->hubCRS().get())) { + interpolationGeogCRS = NN_NO_CHECK( + util::nn_dynamic_pointer_cast( + compSrc0BoundCrs->hubCRS())); + } + } + auto opSrcCRSToGeogCRS = createOperations( + componentsSrc[0], interpolationGeogCRS, context); + auto opGeogCRStoDstCRS = createOperations( + interpolationGeogCRS, componentsDst[0], context); + for (const auto &opSrc : opSrcCRSToGeogCRS) { + for (const auto &opDst : opGeogCRStoDstCRS) { - try { - auto op = createHorizVerticalHorizPROJBased( - sourceCRS, targetCRS, opSrc, - verticalTransform, opDst, - interpolationGeogCRS, true); - res.emplace_back(op); - } catch ( - const InvalidOperationEmptyIntersection &) { - } + try { + auto op = createHorizVerticalHorizPROJBased( + sourceCRS, targetCRS, opSrc, verticalTransform, + opDst, interpolationGeogCRS, true); + res.emplace_back(op); + } catch (const InvalidOperationEmptyIntersection &) { } } } + } - if (verticalTransforms.empty()) { - return createOperations(componentsSrc[0], componentsDst[0], - context); - } + if (verticalTransforms.empty()) { + res = createOperations(componentsSrc[0], componentsDst[0], + context); } } } +} - // '+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +type=crs' to - // '+proj=longlat +ellps=GRS80 +nadgrids=@bar.gsb +geoidgrids=@bar.gtx - // +type=crs' - if (boundSrc && compoundDst) { - const auto &componentsDst = compoundDst->componentReferenceSystems(); - if (!componentsDst.empty()) { - auto compDst0BoundCrs = - dynamic_cast(componentsDst[0].get()); - if (compDst0BoundCrs) { - auto boundSrcHubAsGeogCRS = dynamic_cast( - boundSrc->hubCRS().get()); - auto compDst0BoundCrsHubAsGeogCRS = - dynamic_cast( - compDst0BoundCrs->hubCRS().get()); - if (boundSrcHubAsGeogCRS && compDst0BoundCrsHubAsGeogCRS) { - const auto &boundSrcHubAsGeogCRSDatum = - boundSrcHubAsGeogCRS->datum(); - const auto &compDst0BoundCrsHubAsGeogCRSDatum = - compDst0BoundCrsHubAsGeogCRS->datum(); - if (boundSrcHubAsGeogCRSDatum && - compDst0BoundCrsHubAsGeogCRSDatum && - boundSrcHubAsGeogCRSDatum->_isEquivalentTo( - compDst0BoundCrsHubAsGeogCRSDatum.get())) { - auto cs = cs::EllipsoidalCS:: - createLatitudeLongitudeEllipsoidalHeight( - common::UnitOfMeasure::DEGREE, - common::UnitOfMeasure::METRE); - auto intermGeog3DCRS = util::nn_static_pointer_cast< - crs::CRS>(crs::GeographicCRS::create( - util::PropertyMap() - .set(common::IdentifiedObject::NAME_KEY, - boundSrcHubAsGeogCRS->nameStr()) - .set( - common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY, - metadata::Extent::WORLD), - NN_NO_CHECK(boundSrcHubAsGeogCRSDatum), cs)); - auto sourceToGeog3DOps = createOperations( - sourceCRS, intermGeog3DCRS, context); - auto geog3DToTargetOps = createOperations( - intermGeog3DCRS, targetCRS, context); - for (const auto &opSrc : sourceToGeog3DOps) { - for (const auto &opDst : geog3DToTargetOps) { - if (opSrc->targetCRS() && opDst->sourceCRS() && - !opSrc->targetCRS()->_isEquivalentTo( - opDst->sourceCRS().get())) { - // Shouldn't happen normally, but typically - // one of them can be 2D and the other 3D - // due to above createOperations() not - // exactly setting the expected source and - // target CRS. - // So create an adapter operation... - auto intermOps = createOperations( - NN_NO_CHECK(opSrc->targetCRS()), - NN_NO_CHECK(opDst->sourceCRS()), - context); - if (!intermOps.empty()) { - res.emplace_back( - ConcatenatedOperation:: - createComputeMetadata( - {opSrc, intermOps.front(), - opDst}, - !allowEmptyIntersection)); - } - } else { +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private::createOperationsBoundToCompound( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::BoundCRS *boundSrc, + const crs::CompoundCRS *compoundDst, + std::vector &res) { + + const auto &componentsDst = compoundDst->componentReferenceSystems(); + if (!componentsDst.empty()) { + auto compDst0BoundCrs = + dynamic_cast(componentsDst[0].get()); + if (compDst0BoundCrs) { + auto boundSrcHubAsGeogCRS = + dynamic_cast(boundSrc->hubCRS().get()); + auto compDst0BoundCrsHubAsGeogCRS = + dynamic_cast( + compDst0BoundCrs->hubCRS().get()); + if (boundSrcHubAsGeogCRS && compDst0BoundCrsHubAsGeogCRS) { + const auto &boundSrcHubAsGeogCRSDatum = + boundSrcHubAsGeogCRS->datum(); + const auto &compDst0BoundCrsHubAsGeogCRSDatum = + compDst0BoundCrsHubAsGeogCRS->datum(); + if (boundSrcHubAsGeogCRSDatum && + compDst0BoundCrsHubAsGeogCRSDatum && + boundSrcHubAsGeogCRSDatum->_isEquivalentTo( + compDst0BoundCrsHubAsGeogCRSDatum.get())) { + auto cs = cs::EllipsoidalCS:: + createLatitudeLongitudeEllipsoidalHeight( + common::UnitOfMeasure::DEGREE, + common::UnitOfMeasure::METRE); + auto intermGeog3DCRS = util::nn_static_pointer_cast< + crs::CRS>(crs::GeographicCRS::create( + util::PropertyMap() + .set(common::IdentifiedObject::NAME_KEY, + boundSrcHubAsGeogCRS->nameStr()) + .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY, + metadata::Extent::WORLD), + NN_NO_CHECK(boundSrcHubAsGeogCRSDatum), cs)); + auto sourceToGeog3DOps = + createOperations(sourceCRS, intermGeog3DCRS, context); + auto geog3DToTargetOps = + createOperations(intermGeog3DCRS, targetCRS, context); + for (const auto &opSrc : sourceToGeog3DOps) { + for (const auto &opDst : geog3DToTargetOps) { + if (opSrc->targetCRS() && opDst->sourceCRS() && + !opSrc->targetCRS()->_isEquivalentTo( + opDst->sourceCRS().get())) { + // Shouldn't happen normally, but typically + // one of them can be 2D and the other 3D + // due to above createOperations() not + // exactly setting the expected source and + // target CRS. + // So create an adapter operation... + auto intermOps = createOperations( + NN_NO_CHECK(opSrc->targetCRS()), + NN_NO_CHECK(opDst->sourceCRS()), context); + if (!intermOps.empty()) { res.emplace_back( ConcatenatedOperation:: createComputeMetadata( - {opSrc, opDst}, + {opSrc, intermOps.front(), + opDst}, !allowEmptyIntersection)); } + } else { + res.emplace_back( + ConcatenatedOperation:: + createComputeMetadata( + {opSrc, opDst}, + !allowEmptyIntersection)); } } } @@ -13624,13 +13812,6 @@ CoordinateOperationFactory::Private::createOperations( } } } - - // reverse of previous case - if (boundDst && compoundSrc) { - return applyInverse(createOperations(targetCRS, sourceCRS, context)); - } - - return res; } //! @endcond -- cgit v1.2.3 From ffc865a41aa540673eaedb2552565cf9f8d18679 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 29 Oct 2019 22:20:24 +0100 Subject: Vertical transformations: improve situations similar to transforming from 'NAVD88 (ftUS)' to X, where we now consider the available transformations from 'NAVD88' to X that might exist in the database --- include/proj/io.hpp | 4 + src/iso19111/coordinateoperation.cpp | 479 +++++++++++++++++++++++------------ src/iso19111/factory.cpp | 24 ++ 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 + createVerticalCRSFromDatum(const std::string &datum_auth_name, + const std::string &datum_code) const; + PROJ_INTERNAL std::list 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 &res); + static std::vector + createOperationsGeogToVertWithIntermediateVert( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + const crs::VerticalCRS *vertDst, Context &context); + + static std::vector + 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 &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 - createOperationsGeogToVertWithAlternativeGeog( - const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, - Context &context); - static bool hasPerfectAccuracyResult(const std::vector &res, const Context &context); @@ -12402,51 +12415,6 @@ void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub( // --------------------------------------------------------------------------- -std::vector CoordinateOperationFactory::Private:: - createOperationsGeogToVertWithAlternativeGeog( - const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS - const crs::CRSNNPtr &targetCRS, // vertical CRS - Private::Context &context) { - - std::vector 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(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( - 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 +findCandidateVertCRSForDatum(const io::AuthorityFactoryPtr &authFactory, + const datum::VerticalReferenceFrame *datum) { + std::vector 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( + match.get())); + } + } + } + return candidates; +} + +// --------------------------------------------------------------------------- + +std::vector CoordinateOperationFactory::Private:: + createOperationsGeogToVertWithIntermediateVert( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + const crs::VerticalCRS *vertDst, Private::Context &context) { + + std::vector 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 CoordinateOperationFactory::Private:: + createOperationsGeogToVertWithAlternativeGeog( + const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS + const crs::CRSNNPtr &targetCRS, // vertical CRS + Private::Context &context) { + + std::vector 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(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 &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( + 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(verticalTransform.get()); + if (transformationVerticalTransform == nullptr) { + const auto concat = dynamic_cast( + verticalTransform.get()); + if (concat) { + const auto &steps = concat->operations(); + // Is this change of unit + transformation ? + if (steps.size() == 2 && + dynamic_cast(steps[0].get())) { + transformationVerticalTransform = + dynamic_cast(steps[1].get()); + } + } + } + if (transformationVerticalTransform && + !transformationVerticalTransform->hasBallparkTransformation()) { + auto interpTransformCRS = + transformationVerticalTransform->interpolationCRS(); + if (interpTransformCRS) { + interpolationGeogCRS = + std::dynamic_pointer_cast( + 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( + 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(verticalTransform.get()); - if (transformationVerticalTransform && - !transformationVerticalTransform->hasBallparkTransformation()) { - auto interpTransformCRS = - transformationVerticalTransform->interpolationCRS(); - if (interpTransformCRS) { - interpolationGeogCRS = - std::dynamic_pointer_cast( - 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( - 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 srcToInterpOps; std::vector 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 @@ -4602,6 +4602,30 @@ std::list AuthorityFactory::createGeodeticCRSFromDatum( // --------------------------------------------------------------------------- +//! @cond Doxygen_Suppress +std::list 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 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 AuthorityFactory::createGeodeticCRSFromEllipsoid( 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(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(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(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(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"); -- cgit v1.2.3 From fc769bbd9a4fb61e96e500788d24d1d12019a4d0 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 29 Oct 2019 23:34:59 +0100 Subject: news.rst: fix hyperlink --- docs/source/news.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/news.rst b/docs/source/news.rst index 2a6e8e47..1fba43b9 100644 --- a/docs/source/news.rst +++ b/docs/source/news.rst @@ -10,7 +10,7 @@ News Updates ------- - * Introduced PROJJSON, a JSON encoding of WKT2 (`#1547 `_) + * Introduced PROJJSON, a JSON encoding of WKT2 (`#1547 `_) * Support CRS instantiation of OGC URN's (`#1505 `_) * Expose scope and remarks of database objects (`#1537 `_) -- cgit v1.2.3 From 05a7bcfa56a03437b2ba73616a6bc21c9347d2a7 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 30 Oct 2019 11:28:43 +0100 Subject: Rework importing of Vertical unit change from EPSG db, add support for Height Depth Reversal and use it in createOperations() --- data/sql/conversion.sql | 21 --- data/sql/other_transformation.sql | 21 +++ include/proj/coordinateoperation.hpp | 3 + .../internal/coordinateoperation_constants.hpp | 4 + scripts/build_db.py | 33 +++- scripts/reference_exported_symbols.txt | 1 + src/iso19111/coordinateoperation.cpp | 96 ++++++++++-- src/iso19111/factory.cpp | 31 +++- src/proj_constants.h | 5 + test/unit/test_factory.cpp | 13 +- test/unit/test_operation.cpp | 169 ++++++++++++++++++++- 11 files changed, 350 insertions(+), 47 deletions(-) diff --git a/data/sql/conversion.sql b/data/sql/conversion.sql index 02e96c7d..f42f0e79 100644 --- a/data/sql/conversion.sql +++ b/data/sql/conversion.sql @@ -739,15 +739,6 @@ INSERT INTO "conversion" VALUES('EPSG','7823','CS63 zone X6',NULL,NULL,'EPSG','4 INSERT INTO "conversion" VALUES('EPSG','7824','CS63 zone X7',NULL,NULL,'EPSG','4434','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',0.05,'EPSG','9110','EPSG','8802','Longitude of natural origin',41.3,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.0,'EPSG','9201','EPSG','8806','False easting',7300000.0,'EPSG','9001','EPSG','8807','False northing',0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','7875','St. Helena Local Grid 1971',NULL,NULL,'EPSG','3183','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',-15.58,'EPSG','9110','EPSG','8802','Longitude of natural origin',-5.43,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.0,'EPSG','9201','EPSG','8806','False easting',300000.0,'EPSG','9001','EPSG','8807','False northing',2000000.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','7876','St. Helena Local Grid (Tritan)',NULL,NULL,'EPSG','3183','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',-15.58,'EPSG','9110','EPSG','8802','Longitude of natural origin',-5.43,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.0,'EPSG','9201','EPSG','8806','False easting',299483.737,'EPSG','9001','EPSG','8807','False northing',2000527.879,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','7963','Poolbeg height (ft(Br36)) to Poolbeg height (m)',NULL,NULL,'EPSG','1305','EPSG','1069','Change of Vertical Unit','EPSG','1051','Unit conversion scalar',0.3048007491,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','7972','NGVD29 height (ftUS) to NGVD29 height (m)',NULL,NULL,'EPSG','1323','EPSG','1069','Change of Vertical Unit','EPSG','1051','Unit conversion scalar',0.304800609601219,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','7978','NGVD29 height (ftUS) to NGVD29 depth (ftUS)',NULL,NULL,'EPSG','1323','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','7982','HKPD height to HKPD depth',NULL,NULL,'EPSG','3334','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','7984','KOC WD height to KOC WD depth',NULL,NULL,'EPSG','3267','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','7985','KOC WD depth to KOC WD depth (ft)',NULL,NULL,'EPSG','3267','EPSG','1069','Change of Vertical Unit','EPSG','1051','Unit conversion scalar',3.28083989501312,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','7988','NAVD88 height to NAVD88 height (ftUS)',NULL,NULL,'EPSG','3664','EPSG','1069','Change of Vertical Unit','EPSG','1051','Unit conversion scalar',3.28083333333333,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','7989','NAVD88 height to NAVD88 depth',NULL,NULL,'EPSG','4161','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','7990','NAVD88 height (ftUS) to NAVD88 depth (ftUS)',NULL,NULL,'EPSG','3664','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','7993','Albany Grid 2020',NULL,NULL,'EPSG','4439','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',0.0,'EPSG','9110','EPSG','8802','Longitude of natural origin',117.53,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.0000044,'EPSG','9201','EPSG','8806','False easting',50000.0,'EPSG','9001','EPSG','8807','False northing',4100000.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','7994','Barrow Island Grid 2020',NULL,NULL,'EPSG','4438','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',0.0,'EPSG','9110','EPSG','8802','Longitude of natural origin',115.15,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.0000022,'EPSG','9201','EPSG','8806','False easting',60000.0,'EPSG','9001','EPSG','8807','False northing',2700000.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','7995','Broome Grid 2020',NULL,NULL,'EPSG','4441','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',0.0,'EPSG','9110','EPSG','8802','Longitude of natural origin',122.2,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.00000298,'EPSG','9201','EPSG','8806','False easting',50000.0,'EPSG','9001','EPSG','8807','False northing',2300000.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); @@ -770,15 +761,8 @@ INSERT INTO "conversion" VALUES('EPSG','8011','Perth Coastal Grid',NULL,NULL,'EP INSERT INTO "conversion" VALUES('EPSG','8012','Port Hedland Grid 2020',NULL,NULL,'EPSG','4466','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',0.0,'EPSG','9110','EPSG','8802','Longitude of natural origin',118.36,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.00000135,'EPSG','9201','EPSG','8806','False easting',50000.0,'EPSG','9001','EPSG','8807','False northing',2500000.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8033','TM Zone 20N (US survey feet)',NULL,NULL,'EPSG','4467','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',0.0,'EPSG','9102','EPSG','8802','Longitude of natural origin',-63.0,'EPSG','9102','EPSG','8805','Scale factor at natural origin',0.9996,'EPSG','9201','EPSG','8806','False easting',1640416.667,'EPSG','9003','EPSG','8807','False northing',0.0,'EPSG','9003',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8034','TM Zone 21N (US survey feet)',NULL,NULL,'EPSG','4468','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',0.0,'EPSG','9102','EPSG','8802','Longitude of natural origin',-57.0,'EPSG','9102','EPSG','8805','Scale factor at natural origin',0.9996,'EPSG','9201','EPSG','8806','False easting',1640416.667,'EPSG','9003','EPSG','8807','False northing',0.0,'EPSG','9003',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8038','Instantaneous Water Level height to Instantaneous Water Level depth',NULL,NULL,'EPSG','1262','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8039','MSL height to MSL depth',NULL,NULL,'EPSG','1262','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8040','Gusterberg Grid',NULL,NULL,'EPSG','4455','EPSG','9806','Cassini-Soldner','EPSG','8801','Latitude of natural origin',48.021847,'EPSG','9110','EPSG','8802','Longitude of natural origin',31.481505,'EPSG','9110','EPSG','8806','False easting',0.0,'EPSG','9001','EPSG','8807','False northing',0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8041','St. Stephen Grid',NULL,NULL,'EPSG','4456','EPSG','9806','Cassini-Soldner','EPSG','8801','Latitude of natural origin',48.123154,'EPSG','9110','EPSG','8802','Longitude of natural origin',34.022732,'EPSG','9110','EPSG','8806','False easting',0.0,'EPSG','9001','EPSG','8807','False northing',0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8054','MSL height to MSL height (ft)',NULL,NULL,'EPSG','1262','EPSG','1069','Change of Vertical Unit','EPSG','1051','Unit conversion scalar',3.28083989501312,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8055','MSL height to MSL height (ftUS)',NULL,NULL,'EPSG','1245','EPSG','1069','Change of Vertical Unit','EPSG','1051','Unit conversion scalar',3.28083333333333,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8056','MSL depth to MSL depth (ft)',NULL,NULL,'EPSG','1262','EPSG','1069','Change of Vertical Unit','EPSG','1051','Unit conversion scalar',3.28083989501312,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8057','MSL depth to MSL depth (ftUS)',NULL,NULL,'EPSG','1245','EPSG','1069','Change of Vertical Unit','EPSG','1051','Unit conversion scalar',3.28083333333333,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8060','Baltic 1977 height to Baltic 1977 depth',NULL,NULL,'EPSG','2423','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8061','Pima County zone 1 East (ft)',NULL,NULL,'EPSG','4472','EPSG','9815','Hotine Oblique Mercator (variant B)','EPSG','8811','Latitude of projection centre',32.15,'EPSG','9110','EPSG','8812','Longitude of projection centre',-111.24,'EPSG','9110','EPSG','8813','Azimuth of initial line',45.0,'EPSG','9102','EPSG','8814','Angle from Rectified to Skew Grid',45.0,'EPSG','9102','EPSG','8815','Scale factor on initial line',1.00011,'EPSG','9201','EPSG','8816','Easting at projection centre',160000.0,'EPSG','9002','EPSG','8817','Northing at projection centre',800000.0,'EPSG','9002',0); INSERT INTO "conversion" VALUES('EPSG','8062','Pima County zone 2 Central (ft)',NULL,NULL,'EPSG','4460','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',31.15,'EPSG','9110','EPSG','8802','Longitude of natural origin',-112.1,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.00009,'EPSG','9201','EPSG','8806','False easting',1800000.0,'EPSG','9002','EPSG','8807','False northing',1000000.0,'EPSG','9002',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8063','Pima County zone 3 West (ft)',NULL,NULL,'EPSG','4450','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',31.3,'EPSG','9110','EPSG','8802','Longitude of natural origin',-113.1,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.000055,'EPSG','9201','EPSG','8806','False easting',600000.0,'EPSG','9002','EPSG','8807','False northing',0.0,'EPSG','9002',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); @@ -786,7 +770,6 @@ INSERT INTO "conversion" VALUES('EPSG','8064','Pima County zone 4 Mt. Lemmon (ft INSERT INTO "conversion" VALUES('EPSG','8080','MTM Nova Scotia zone 4 v2',NULL,NULL,'EPSG','1534','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',0.0,'EPSG','9110','EPSG','8802','Longitude of natural origin',-61.3,'EPSG','9110','EPSG','8805','Scale factor at natural origin',0.9999,'EPSG','9201','EPSG','8806','False easting',24500000.0,'EPSG','9001','EPSG','8807','False northing',0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8081','MTM Nova Scotia zone 5 v2',NULL,NULL,'EPSG','1535','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',0.0,'EPSG','9110','EPSG','8802','Longitude of natural origin',-64.3,'EPSG','9110','EPSG','8805','Scale factor at natural origin',0.9999,'EPSG','9201','EPSG','8806','False easting',25500000.0,'EPSG','9001','EPSG','8807','False northing',0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8087','Iceland Lambert 2016',NULL,NULL,'EPSG','1120','EPSG','9802','Lambert Conic Conformal (2SP)','EPSG','8821','Latitude of false origin',65.0,'EPSG','9110','EPSG','8822','Longitude of false origin',-19.0,'EPSG','9110','EPSG','8823','Latitude of 1st standard parallel',64.15,'EPSG','9110','EPSG','8824','Latitude of 2nd standard parallel',65.45,'EPSG','9110','EPSG','8826','Easting at false origin',2700000.0,'EPSG','9001','EPSG','8827','Northing at false origin',300000.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8229','NAVD88 height to NAVD88 height (ft)',NULL,NULL,'EPSG','4464','EPSG','1069','Change of Vertical Unit','EPSG','1051','Unit conversion scalar',3.28083989501312,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8273','Oregon Burns-Harper zone (meters)',NULL,NULL,'EPSG','4459','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',43.3,'EPSG','9110','EPSG','8802','Longitude of natural origin',-117.4,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.00014,'EPSG','9201','EPSG','8806','False easting',90000.0,'EPSG','9001','EPSG','8807','False northing',0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8274','Oregon Burns-Harper zone (International feet)',NULL,NULL,'EPSG','4459','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',43.3,'EPSG','9110','EPSG','8802','Longitude of natural origin',-117.4,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.00014,'EPSG','9201','EPSG','8806','False easting',295275.5906,'EPSG','9002','EPSG','8807','False northing',0.0,'EPSG','9002',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8275','Oregon Canyon City-Burns zone (meters)',NULL,NULL,'EPSG','4465','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',43.3,'EPSG','9110','EPSG','8802','Longitude of natural origin',-119.0,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.00022,'EPSG','9201','EPSG','8806','False easting',20000.0,'EPSG','9001','EPSG','8807','False northing',0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); @@ -825,10 +808,6 @@ INSERT INTO "conversion" VALUES('EPSG','8307','Oregon Warner Highway zone (meter INSERT INTO "conversion" VALUES('EPSG','8308','Oregon Warner Highway zone (International feet)',NULL,NULL,'EPSG','4486','EPSG','9801','Lambert Conic Conformal (1SP)','EPSG','8801','Latitude of natural origin',42.3,'EPSG','9110','EPSG','8802','Longitude of natural origin',-120.0,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.000245,'EPSG','9201','EPSG','8806','False easting',131233.5958,'EPSG','9002','EPSG','8807','False northing',196850.3937,'EPSG','9002',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8309','Oregon Willamette Pass zone (meters)',NULL,NULL,'EPSG','4488','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',43.0,'EPSG','9110','EPSG','8802','Longitude of natural origin',-122.0,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.000223,'EPSG','9201','EPSG','8806','False easting',20000.0,'EPSG','9001','EPSG','8807','False northing',0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8310','Oregon Willamette Pass zone (International feet)',NULL,NULL,'EPSG','4488','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',43.0,'EPSG','9110','EPSG','8802','Longitude of natural origin',-122.0,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.000223,'EPSG','9201','EPSG','8806','False easting',65616.7979,'EPSG','9002','EPSG','8807','False northing',0.0,'EPSG','9002',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8354','Black Sea height to Black Sea depth',NULL,NULL,'EPSG','3251','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8355','AIOC95 height to AIOC95 depth',NULL,NULL,'EPSG','2592','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8356','Caspian height to Caspian depth',NULL,NULL,'EPSG','1291','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "conversion" VALUES('EPSG','8359','Baltic 1957 height to Baltic 1957 depth',NULL,NULL,'EPSG','1306','EPSG','1068','Height Depth Reversal',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8373','NCRS Las Vegas zone (m)',NULL,NULL,'EPSG','4485','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',36.15,'EPSG','9110','EPSG','8802','Longitude of natural origin',-114.58,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.0001,'EPSG','9201','EPSG','8806','False easting',100000.0,'EPSG','9001','EPSG','8807','False northing',200000.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8374','NCRS Las Vegas zone (ftUS)',NULL,NULL,'EPSG','4485','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',36.15,'EPSG','9110','EPSG','8802','Longitude of natural origin',-114.58,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.0001,'EPSG','9201','EPSG','8806','False easting',328083.3333,'EPSG','9003','EPSG','8807','False northing',656166.6667,'EPSG','9003',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); INSERT INTO "conversion" VALUES('EPSG','8375','NCRS Las Vegas high elevation zone (m)',NULL,NULL,'EPSG','4487','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',36.15,'EPSG','9110','EPSG','8802','Longitude of natural origin',-114.58,'EPSG','9110','EPSG','8805','Scale factor at natural origin',1.000135,'EPSG','9201','EPSG','8806','False easting',300000.0,'EPSG','9001','EPSG','8807','False northing',400000.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); diff --git a/data/sql/other_transformation.sql b/data/sql/other_transformation.sql index 152e5f82..b4ad598c 100644 --- a/data/sql/other_transformation.sql +++ b/data/sql/other_transformation.sql @@ -192,11 +192,32 @@ INSERT INTO "other_transformation" VALUES('EPSG','7653','EGM96 height to Kumul 3 INSERT INTO "other_transformation" VALUES('EPSG','7654','EGM2008 height to Kiunga height (1)','','Derivation of heights referenced to Kiunga vertical CRS (CRS code 7652).','EPSG','9616','Vertical Offset','EPSG','3855','EPSG','7652','EPSG','4383',0.0,'EPSG','8603','Vertical Offset',-3.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'QC-Png Kiunga',0); INSERT INTO "other_transformation" VALUES('EPSG','7873','EGM96 height to POM96 height (1)','','Derivation of POM96 heights.','EPSG','9616','Vertical Offset','EPSG','5773','EPSG','7832','EPSG','4425',0.0,'EPSG','8603','Vertical Offset',-1.58,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'QC-Png Gulf-Cen',0); INSERT INTO "other_transformation" VALUES('EPSG','7874','EGM2008 height to POM08 height (1)','','Derivation of POM08 heights.','EPSG','9616','Vertical Offset','EPSG','3855','EPSG','7841','EPSG','4425',0.0,'EPSG','8603','Vertical Offset',-0.93,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'QC-Png Gulf-Cen',0); +INSERT INTO "other_transformation" VALUES('EPSG','7963','Poolbeg height (ft(Br36)) to Poolbeg height (m)','Change of unit from British foot (1936) [ft(BR36)] to metre.','Change of unit to facilitate transformation of heights or depths through concatenated operations.','EPSG','1069','Change of Vertical Unit','EPSG','5754','EPSG','7962','EPSG','1305',NULL,'EPSG','1051','Unit conversion scalar',0.3048007491,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); INSERT INTO "other_transformation" VALUES('EPSG','7964','Poolbeg height (m) to Malin Head height (1)','Poolbeg vertical reference surface is 2.7m below Malin Head surface. Transformations are subject to localised anomalies.','Change of height to a different vertical reference surface for topographic mapping and engineering survey. Accuracy 0.1m.','EPSG','9616','Vertical Offset','EPSG','7962','EPSG','5731','EPSG','1305',0.1,'EPSG','8603','Vertical Offset',-2.7,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'OSI-Ire',0); INSERT INTO "other_transformation" VALUES('EPSG','7966','Poolbeg height (m) to Belfast height (1)','Poolbeg vertical reference surface is 2.7m below Belfast surface. Transformations are subject to localised anomalies.','Change of height to a different vertical reference surface for topographic mapping and engineering survey. Accuracy 0.1m.','EPSG','9616','Vertical Offset','EPSG','7962','EPSG','5732','EPSG','1305',0.1,'EPSG','8603','Vertical Offset',-2.7,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'OSI-Ire',0); +INSERT INTO "other_transformation" VALUES('EPSG','7972','NGVD29 height (ftUS) to NGVD29 height (m)','Change of unit from US survey foot (ftUS) to metre. 1 ftUS = (12/39.37)m ≈ 0.304800609601219m.','Change of unit to facilitate transformation of heights or depths through concatenated operations.','EPSG','1069','Change of Vertical Unit','EPSG','5702','EPSG','7968','EPSG','1323',NULL,'EPSG','1051','Unit conversion scalar',0.304800609601219,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); INSERT INTO "other_transformation" VALUES('EPSG','7977','HKPD depth to HKCD depth (1)','The HKPD vertical reference surface is 0.146m above the HKCD surface.','Change of depth to a different vertical reference surface.','EPSG','9616','Vertical Offset','EPSG','7976','EPSG','5739','EPSG','3335',0.0,'EPSG','8603','Vertical Offset',-0.146,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'SMO-HK',0); +INSERT INTO "other_transformation" VALUES('EPSG','7978','NGVD29 height (ftUS) to NGVD29 depth (ftUS)','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','5702','EPSG','6359','EPSG','1323',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); INSERT INTO "other_transformation" VALUES('EPSG','7980','KOC CD height to KOC WD height (1)','The KOC CD vertical reference surface is 4.74m (15.55ft) below KOC WD surface.','Change of height to a different vertical reference surface.','EPSG','9616','Vertical Offset','EPSG','5790','EPSG','7979','EPSG','3267',0.1,'EPSG','8603','Vertical Offset',-4.74,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'KOC-Kwt',0); INSERT INTO "other_transformation" VALUES('EPSG','7981','Kuwait PWD height to KOC WD height (1)','The KOC WD vertical reference surface is 4.25m above the Kuwait PWD surface.','Change of height to a different vertical reference surface.','EPSG','9616','Vertical Offset','EPSG','5788','EPSG','7979','EPSG','3267',0.1,'EPSG','8603','Vertical Offset',-4.25,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'KOC-Kwt',0); +INSERT INTO "other_transformation" VALUES('EPSG','7982','HKPD height to HKPD depth','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','5738','EPSG','7976','EPSG','3334',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','7984','KOC WD height to KOC WD depth','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','7979','EPSG','5789','EPSG','3267',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','7985','KOC WD depth to KOC WD depth (ft)','Change of unit from metre to International foot (ft). 1ft = 0.3048m.','Change of unit to facilitate transformation of heights or depths through concatenated operations.','EPSG','1069','Change of Vertical Unit','EPSG','5789','EPSG','5614','EPSG','3267',NULL,'EPSG','1051','Unit conversion scalar',3.28083989501312,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','7988','NAVD88 height to NAVD88 height (ftUS)','Change of unit from metre to US survey foot. 1 ftUS = (12/39.37)m ≈ 0.304800609601219m.','Change of unit to facilitate transformation of heights or depths through concatenated operations.','EPSG','1069','Change of Vertical Unit','EPSG','5703','EPSG','6360','EPSG','3664',NULL,'EPSG','1051','Unit conversion scalar',3.28083333333333,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','7989','NAVD88 height to NAVD88 depth','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','5703','EPSG','6357','EPSG','4161',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','7990','NAVD88 height (ftUS) to NAVD88 depth (ftUS)','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','6360','EPSG','6358','EPSG','3664',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8038','Instantaneous Water Level height to Instantaneous Water Level depth','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','5829','EPSG','5831','EPSG','1262',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8039','MSL height to MSL depth','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','5714','EPSG','5715','EPSG','1262',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8054','MSL height to MSL height (ft)','Change of unit from metre to International foot (ft). 1ft = 0.3048m.','Change of unit to facilitate transformation of heights or depths through concatenated operations.','EPSG','1069','Change of Vertical Unit','EPSG','5714','EPSG','8050','EPSG','1262',NULL,'EPSG','1051','Unit conversion scalar',3.28083989501312,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8055','MSL height to MSL height (ftUS)','Change of unit from metre to US survey foot (ftUS). 1 ftUS = (12/39.37)m ≈ 0.304800609601219m.','Change of unit to facilitate transformation of heights or depths through concatenated operations.','EPSG','1069','Change of Vertical Unit','EPSG','5714','EPSG','8052','EPSG','1245',NULL,'EPSG','1051','Unit conversion scalar',3.28083333333333,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8056','MSL depth to MSL depth (ft)','Change of unit from metre to International foot (ft). 1ft = 0.3048m.','Change of unit to facilitate transformation of heights or depths through concatenated operations.','EPSG','1069','Change of Vertical Unit','EPSG','5715','EPSG','8051','EPSG','1262',NULL,'EPSG','1051','Unit conversion scalar',3.28083989501312,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8057','MSL depth to MSL depth (ftUS)','Change of unit from metre to US survey foot (ftUS). 1 ftUS = (12/39.37)m ≈ 0.304800609601219m.','Change of unit to facilitate transformation of heights or depths through concatenated operations.','EPSG','1069','Change of Vertical Unit','EPSG','5715','EPSG','8053','EPSG','1245',NULL,'EPSG','1051','Unit conversion scalar',3.28083333333333,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8060','Baltic 1977 height to Baltic 1977 depth','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','5705','EPSG','5612','EPSG','2423',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8229','NAVD88 height to NAVD88 height (ft)','Change of unit from metre to International foot (ft). 1ft = 0.3048m.','Change of unit to facilitate transformation of heights or depths through concatenated operations for States of the US which have adopted International feet for their State Plane coordinate systems.','EPSG','1069','Change of Vertical Unit','EPSG','5703','EPSG','8228','EPSG','4464',NULL,'EPSG','1051','Unit conversion scalar',3.28083989501312,'EPSG','9201',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8354','Black Sea height to Black Sea depth','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','5735','EPSG','5336','EPSG','3251',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8355','AIOC95 height to AIOC95 depth','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','5797','EPSG','5734','EPSG','2592',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8356','Caspian height to Caspian depth','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','5611','EPSG','5706','EPSG','1291',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); +INSERT INTO "other_transformation" VALUES('EPSG','8359','Baltic 1957 height to Baltic 1957 depth','Change of axis positive direction from up to down.','Change of axis positive direction to facilitate transformation of heights or depths through concatenated operations.','EPSG','1068','Height Depth Reversal','EPSG','8357','EPSG','8358','EPSG','1306',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0); INSERT INTO "other_transformation" VALUES('EPSG','10087','Jamaica 1875 / Jamaica (Old Grid) to JAD69 / Jamaica National Grid (1)','Derived by least squares fit at primary triangulation stations. Accuracy will be less outside of this network due to extrapolation.','Topographic mapping.','EPSG','9624','Affine parametric transformation','EPSG','24100','EPSG','24200','EPSG','3342',1.5,'EPSG','8623','A0',82357.457,'EPSG','9001','EPSG','8624','A1',0.304794369,'EPSG','9203','EPSG','8625','A2',1.5417425e-05,'EPSG','9203','EPSG','8639','B0',28091.324,'EPSG','9001','EPSG','8640','B1',-1.5417425e-05,'EPSG','9203','EPSG','8641','B2',0.304794369,'EPSG','9203',NULL,NULL,NULL,NULL,NULL,NULL,'SD-Jam',0); INSERT INTO "other_transformation" VALUES('EPSG','10088','JAD69 / Jamaica National Grid to Jamaica 1875 / Jamaica (Old Grid) (1)','Derived by least squares fit at primary triangulation stations. Accuracy will be less outside of this network due to extrapolation.','Topographic mapping.','EPSG','9624','Affine parametric transformation','EPSG','24200','EPSG','24100','EPSG','3342',1.5,'EPSG','8623','A0',-270201.96,'EPSG','9005','EPSG','8624','A1',0.00016595792,'EPSG','9203','EPSG','8625','A2',3.2809005,'EPSG','9203','EPSG','8639','B0',-92178.51,'EPSG','9005','EPSG','8640','B1',3.2809005,'EPSG','9203','EPSG','8641','B2',-0.00016595792,'EPSG','9203',NULL,NULL,NULL,NULL,NULL,NULL,'SD-Jam',1); INSERT INTO "other_transformation" VALUES('EPSG','10095','Mauritania 1999 / UTM zone 28N to WGS 84 / UTM zone 28N (1)','Parameter values consistent with the OGP Affine parametric transformation method derived by OGP from the published Helmert 2D parameter values.','Minerals management.','EPSG','9624','Affine parametric transformation','EPSG','3103','EPSG','32628','EPSG','2971',40.0,'EPSG','8623','A0',NULL,'EPSG',NULL,'EPSG','8624','A1',NULL,'EPSG',NULL,'EPSG','8625','A2',NULL,'EPSG',NULL,'EPSG','8639','B0',NULL,'EPSG',NULL,'EPSG','8640','B1',NULL,'EPSG',NULL,'EPSG','8641','B2',NULL,'EPSG',NULL,NULL,NULL,NULL,NULL,NULL,NULL,'MMI-Mau W',1); diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp index 6c4c25c2..1ced5333 100644 --- a/include/proj/coordinateoperation.hpp +++ b/include/proj/coordinateoperation.hpp @@ -1323,6 +1323,9 @@ class PROJ_GCC_DLL Conversion : public SingleOperation { createChangeVerticalUnit(const util::PropertyMap &properties, const common::Scale &factor); + PROJ_DLL static ConversionNNPtr + createHeightDepthReversal(const util::PropertyMap &properties); + PROJ_DLL static ConversionNNPtr createAxisOrderReversal(bool is3D); PROJ_DLL static ConversionNNPtr diff --git a/include/proj/internal/coordinateoperation_constants.hpp b/include/proj/internal/coordinateoperation_constants.hpp index f1925c9b..eb0bb8c5 100644 --- a/include/proj/internal/coordinateoperation_constants.hpp +++ b/include/proj/internal/coordinateoperation_constants.hpp @@ -826,6 +826,7 @@ static const struct MethodNameCode { METHOD_NAME_CODE(VERTICAL_PERSPECTIVE), // Other conversions METHOD_NAME_CODE(CHANGE_VERTICAL_UNIT), + METHOD_NAME_CODE(HEIGHT_DEPTH_REVERSAL), METHOD_NAME_CODE(AXIS_ORDER_REVERSAL_2D), METHOD_NAME_CODE(AXIS_ORDER_REVERSAL_3D), METHOD_NAME_CODE(GEOGRAPHIC_GEOCENTRIC), @@ -1153,6 +1154,9 @@ static const MethodMapping otherMethodMappings[] = { {EPSG_NAME_METHOD_CHANGE_VERTICAL_UNIT, EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT, nullptr, nullptr, nullptr, paramsChangeVerticalUnit}, + {EPSG_NAME_METHOD_HEIGHT_DEPTH_REVERSAL, + EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL, nullptr, nullptr, nullptr, + paramsChangeVerticalUnit}, {EPSG_NAME_METHOD_AXIS_ORDER_REVERSAL_2D, EPSG_CODE_METHOD_AXIS_ORDER_REVERSAL_2D, nullptr, nullptr, nullptr, nullptr}, diff --git a/scripts/build_db.py b/scripts/build_db.py index f52883cf..ba1502ae 100755 --- a/scripts/build_db.py +++ b/scripts/build_db.py @@ -137,7 +137,12 @@ BEFORE INSERT ON conversion BEGIN """ - proj_db_cursor.execute("SELECT coord_op_code, coord_op_name, area_of_use_code, coord_op_method_code, coord_op_method_name, epsg_coordoperation.deprecated FROM epsg.epsg_coordoperation LEFT JOIN epsg.epsg_coordoperationmethod USING (coord_op_method_code) WHERE coord_op_type = 'conversion' AND coord_op_name NOT LIKE '%to DMSH'") + # 1068 and 1069 are Height Depth Reversal and Change of Vertical Unit + # In EPSG, there is one generic instance of those as 7812 and 7813 that + # don't refer to particular CRS, and instances pointing to CRS names + # The later are imported in the other_transformation table since we recover + # the source/target CRS names from the transformation name. + proj_db_cursor.execute("SELECT coord_op_code, coord_op_name, area_of_use_code, coord_op_method_code, coord_op_method_name, epsg_coordoperation.deprecated FROM epsg.epsg_coordoperation LEFT JOIN epsg.epsg_coordoperationmethod USING (coord_op_method_code) WHERE coord_op_type = 'conversion' AND coord_op_name NOT LIKE '%to DMSH' AND (coord_op_method_code NOT IN (1068, 1069) OR coord_op_code IN (7812,7813))") for (code, name, area_of_use_code, method_code, method_name, deprecated) in proj_db_cursor.fetchall(): expected_order = 1 max_n_params = 7 @@ -423,8 +428,31 @@ def fill_other_transformation(proj_db_cursor): # 9619: Geographic2D offsets # 9624: Affine Parametric Transformation # 9660: Geographic3D offsets - proj_db_cursor.execute("SELECT coord_op_code, coord_op_name, coord_op_method_code, coord_op_method_name, source_crs_code, target_crs_code, area_of_use_code, coord_op_accuracy, coord_tfm_version, epsg_coordoperation.deprecated, coord_op_scope, epsg_coordoperation.remarks FROM epsg.epsg_coordoperation LEFT JOIN epsg.epsg_coordoperationmethod USING (coord_op_method_code) WHERE coord_op_type = 'transformation' AND coord_op_method_code IN (9601, 9616, 9618, 9619, 9624, 9660)") + # 1068: Height Depth Reversal + # 1069: Change of Vertical Unit + proj_db_cursor.execute("SELECT coord_op_code, coord_op_name, coord_op_method_code, coord_op_method_name, source_crs_code, target_crs_code, area_of_use_code, coord_op_accuracy, coord_tfm_version, epsg_coordoperation.deprecated, coord_op_scope, epsg_coordoperation.remarks FROM epsg.epsg_coordoperation LEFT JOIN epsg.epsg_coordoperationmethod USING (coord_op_method_code) WHERE coord_op_method_code IN (9601, 9616, 9618, 9619, 9624, 9660, 1068, 1069)") for (code, name, method_code, method_name, source_crs_code, target_crs_code, area_of_use_code, coord_op_accuracy, coord_tfm_version, deprecated, scope, remarks) in proj_db_cursor.fetchall(): + + # 1068 and 1069 are Height Depth Reversal and Change of Vertical Unit + # In EPSG, there is one generic instance of those as 7812 and 7813 that + # don't refer to particular CRS, and instances pointing to CRS names + # The later are imported in the other_transformation table since we recover + # the source/target CRS names from the transformation name. + if method_code in (1068, 1069) and source_crs_code is None and target_crs_code is None: + parts = name.split(" to ") + if len(parts) != 2: + continue + + proj_db_cursor.execute("SELECT coord_ref_sys_code FROM epsg_coordinatereferencesystem WHERE coord_ref_sys_name = ?", (parts[0],)) + source_codes = proj_db_cursor.fetchall() + proj_db_cursor.execute("SELECT coord_ref_sys_code FROM epsg_coordinatereferencesystem WHERE coord_ref_sys_name = ?", (parts[1],)) + target_codes = proj_db_cursor.fetchall() + if len(source_codes) != 1 and len(target_codes) != 1: + continue + + source_crs_code = source_codes[0][0] + target_crs_code = target_codes[0][0] + expected_order = 1 max_n_params = 7 param_auth_name = [None for i in range(max_n_params)] @@ -471,6 +499,7 @@ def fill_other_transformation(proj_db_cursor): deprecated) #proj_db_cursor.execute("INSERT INTO coordinate_operation VALUES (?,?,'other_transformation')", (EPSG_AUTHORITY, code)) + #print(arg) proj_db_cursor.execute('INSERT INTO other_transformation VALUES (' + '?,?,?, ?,?, ?,?,?, ?,?, ?,?, ?,?, ?, ?,?,?,?,?,?, ?,?,?,?,?,?, ?,?,?,?,?,?, ' + '?,?,?,?,?,?, ?,?,?,?,?,?, ?,?,?,?,?,?, ?,?,?,?,?,?, ?,?)', arg) diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt index 73a36a70..34e23b22 100644 --- a/scripts/reference_exported_symbols.txt +++ b/scripts/reference_exported_symbols.txt @@ -495,6 +495,7 @@ osgeo::proj::operation::Conversion::createGeostationarySatelliteSweepY(osgeo::pr osgeo::proj::operation::Conversion::createGnomonic(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&) osgeo::proj::operation::Conversion::createGoodeHomolosine(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&) osgeo::proj::operation::Conversion::createGuamProjection(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&) +osgeo::proj::operation::Conversion::createHeightDepthReversal(osgeo::proj::util::PropertyMap const&) osgeo::proj::operation::Conversion::createHotineObliqueMercatorTwoPointNaturalOrigin(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Scale const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&) osgeo::proj::operation::Conversion::createHotineObliqueMercatorVariantA(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Scale const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&) osgeo::proj::operation::Conversion::createHotineObliqueMercatorVariantB(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Scale const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 27ac712c..5b38a6a4 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -4661,6 +4661,26 @@ Conversion::createChangeVerticalUnit(const util::PropertyMap &properties, // --------------------------------------------------------------------------- +/** \brief Instantiate a conversion based on the Height Depth Reversal + * method. + * + * This method is defined as [EPSG:1068] + * (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::1068) + * + * @param properties See \ref general_properties of the conversion. If the name + * is not provided, it is automatically set. + * @return a new Conversion. + * @since 7.0 + */ +ConversionNNPtr +Conversion::createHeightDepthReversal(const util::PropertyMap &properties) { + return create(properties, createMethodMapNameEPSGCode( + EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL), + {}, {}); +} + +// --------------------------------------------------------------------------- + /** \brief Instantiate a conversion based on the Axis order reversal method * * This swaps the longitude, latitude axis. @@ -4927,6 +4947,14 @@ CoordinateOperationNNPtr Conversion::inverse() const { return conv; } + if (methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL) { + + auto conv = createHeightDepthReversal( + createPropertiesForInverse(this, false, false)); + conv->setCRSs(this, true); + return conv; + } + return InverseConversion::create(NN_NO_CHECK( util::nn_dynamic_pointer_cast(shared_from_this()))); } @@ -5808,9 +5836,12 @@ void Conversion::_exportToPROJString( methodEPSGCode == EPSG_CODE_METHOD_AFFINE_PARAMETRIC_TRANSFORMATION; const bool isGeographicGeocentric = methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_GEOCENTRIC; + const bool isHeightDepthReversal = + methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL; const bool applySourceCRSModifiers = !isZUnitConversion && !isAffineParametric && - !isAxisOrderReversal(methodEPSGCode) && !isGeographicGeocentric; + !isAxisOrderReversal(methodEPSGCode) && !isGeographicGeocentric && + !isHeightDepthReversal; bool applyTargetCRSModifiers = applySourceCRSModifiers; auto l_sourceCRS = sourceCRS(); @@ -7876,6 +7907,17 @@ TransformationNNPtr Transformation::inverseAsTransformation() const { coordinateOperationAccuracies())); } +#ifdef notdef + // We dont need that currently, but we might... + if (methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL) { + return d->registerInv( + shared_from_this(), + createHeightDepthReversal( + createPropertiesForInverse(this, false, false), l_targetCRS, + l_sourceCRS, coordinateOperationAccuracies())); + } +#endif + return InverseTransformation::create(NN_NO_CHECK( util::nn_dynamic_pointer_cast(shared_from_this()))); } @@ -9389,6 +9431,12 @@ bool SingleOperation::exportToPROJStringGeneric( return true; } + if (methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL) { + formatter->addStep("axisswap"); + formatter->addParam("order", "1,2,-3"); + return true; + } + return false; } @@ -13306,10 +13354,16 @@ void CoordinateOperationFactory::Private::createOperationsVertToVert( srcDatum->_isEquivalentTo(dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)); - const double convSrc = - vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI(); - const double convDst = - vertDst->coordinateSystem()->axisList()[0]->unit().conversionToSI(); + const auto &srcAxis = vertSrc->coordinateSystem()->axisList()[0]; + const double convSrc = srcAxis->unit().conversionToSI(); + const auto &dstAxis = vertDst->coordinateSystem()->axisList()[0]; + const double convDst = dstAxis->unit().conversionToSI(); + const bool srcIsUp = srcAxis->direction() == cs::AxisDirection::UP; + const bool srcIsDown = srcAxis->direction() == cs::AxisDirection::DOWN; + const bool dstIsUp = dstAxis->direction() == cs::AxisDirection::UP; + const bool dstIsDown = dstAxis->direction() == cs::AxisDirection::DOWN; + const bool heightDepthReversal = + ((srcIsUp && dstIsDown) || (srcIsDown && dstIsUp)); const double factor = convSrc / convDst; auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()); @@ -13317,13 +13371,23 @@ void CoordinateOperationFactory::Private::createOperationsVertToVert( name += BALLPARK_VERTICAL_TRANSFORMATION; auto conv = Transformation::createChangeVerticalUnit( util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name), - sourceCRS, targetCRS, common::Scale(factor), {}); + sourceCRS, targetCRS, + // In case of a height depth reversal, we should probably have + // 2 steps instead of putting a negative factor... + common::Scale(heightDepthReversal ? -factor : factor), {}); conv->setHasBallparkTransformation(true); res.push_back(conv); } else if (convSrc != convDst) { auto conv = Conversion::createChangeVerticalUnit( util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name), - common::Scale(factor)); + // In case of a height depth reversal, we should probably have + // 2 steps instead of putting a negative factor... + common::Scale(heightDepthReversal ? -factor : factor)); + conv->setCRSs(sourceCRS, targetCRS, nullptr); + res.push_back(conv); + } else if (heightDepthReversal) { + auto conv = Conversion::createHeightDepthReversal( + util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name)); conv->setCRSs(sourceCRS, targetCRS, nullptr); res.push_back(conv); } @@ -13536,11 +13600,19 @@ getInterpolationGeogCRS(const CoordinateOperationNNPtr &verticalTransform, verticalTransform.get()); if (concat) { const auto &steps = concat->operations(); - // Is this change of unit + transformation ? - if (steps.size() == 2 && - dynamic_cast(steps[0].get())) { - transformationVerticalTransform = - dynamic_cast(steps[1].get()); + // Is this change of unit and/or height depth reversal + + // transformation ? + for (const auto &step : steps) { + const auto transf = + dynamic_cast(step.get()); + if (transf) { + // Only support a single Transformation in the steps + if (transformationVerticalTransform != nullptr) { + transformationVerticalTransform = nullptr; + break; + } + transformationVerticalTransform = transf; + } } } } diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index d3ea8f26..44e44dbe 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -2295,6 +2295,20 @@ AuthorityFactory::createConversion(const std::string &code) const { auto res = d->runWithCodeParam(sql, code); if (res.empty()) { + try { + // Conversions using methods Change of Vertical Unit or + // Height Depth Reveral are stored in other_transformation + auto op = createCoordinateOperation( + code, false /* allowConcatenated */, + false /* usePROJAlternativeGridNames */, + "other_transformation"); + auto conv = + util::nn_dynamic_pointer_cast(op); + if (conv) { + return NN_NO_CHECK(conv); + } + } catch (const std::exception &) { + } throw NoSuchAuthorityCodeException("conversion not found", d->authority(), code); } @@ -3118,13 +3132,16 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation( .set(metadata::Identifier::CODE_KEY, method_code) .set(common::IdentifiedObject::NAME_KEY, method_name); - if (method_auth_name == metadata::Identifier::EPSG && - operation::isAxisOrderReversal( - std::atoi(method_code.c_str()))) { - auto op = operation::Conversion::create(props, propsMethod, - parameters, values); - op->setCRSs(sourceCRS, targetCRS, nullptr); - return op; + if (method_auth_name == metadata::Identifier::EPSG) { + int method_code_int = std::atoi(method_code.c_str()); + if (operation::isAxisOrderReversal(method_code_int) || + method_code_int == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT || + method_code_int == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL) { + auto op = operation::Conversion::create(props, propsMethod, + parameters, values); + op->setCRSs(sourceCRS, targetCRS, nullptr); + return op; + } } return operation::Transformation::create( props, sourceCRS, targetCRS, nullptr, propsMethod, parameters, diff --git a/src/proj_constants.h b/src/proj_constants.h index 619d9d91..62cf94be 100644 --- a/src/proj_constants.h +++ b/src/proj_constants.h @@ -592,4 +592,9 @@ #define EPSG_NAME_METHOD_AXIS_ORDER_REVERSAL_3D \ "Axis Order Reversal (Geographic3D horizontal)" +/* ------------------------------------------------------------------------ */ + +#define EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL 1068 +#define EPSG_NAME_METHOD_HEIGHT_DEPTH_REVERSAL "Height Depth Reversal" + #endif /* PROJ_CONSTANTS_INCLUDED */ diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp index 47cee060..94010135 100644 --- a/test/unit/test_factory.cpp +++ b/test/unit/test_factory.cpp @@ -542,6 +542,15 @@ TEST(factory, AuthorityFactory_createConversion) { // --------------------------------------------------------------------------- +TEST(factory, AuthorityFactory_createConversion_from_other_transformation) { + auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto op = factory->createCoordinateOperation("7984", false); + auto conversion = nn_dynamic_pointer_cast(op); + ASSERT_TRUE(conversion != nullptr); +} + +// --------------------------------------------------------------------------- + TEST(factory, AuthorityFactory_createProjectedCRS) { auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); EXPECT_THROW(factory->createProjectedCRS("-1"), @@ -1157,9 +1166,7 @@ TEST(factory, AuthorityFactory_build_all_concatenated) { AuthorityFactory::ObjectType::CONCATENATED_OPERATION, false); EXPECT_LT(setConcatenatedNoDeprecated.size(), setConcatenated.size()); for (const auto &code : setConcatenated) { - if (in(code, {"8422", "8481", "8482", "8565", "8566", "8572", - // the issue with 7987 is the chaining of two conversions - "7987"})) { + if (in(code, {"8422", "8481", "8482", "8565", "8566", "8572"})) { EXPECT_THROW(factory->createCoordinateOperation(code, false), FactoryException) << code; diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 7fbbf70e..f7a86eb4 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -6872,6 +6872,43 @@ TEST(operation, vertCRS_to_vertCRS) { EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), "+proj=affine +s33=0.999998"); } + + auto vertCRSMetreUp = + nn_dynamic_pointer_cast(WKTParser().createFromWKT( + "VERTCRS[\"my height\",VDATUM[\"my datum\"],CS[vertical,1]," + "AXIS[\"gravity-related height (H)\",up," + "LENGTHUNIT[\"metre\",1]]]")); + ASSERT_TRUE(vertCRSMetreUp != nullptr); + + auto vertCRSMetreDown = + nn_dynamic_pointer_cast(WKTParser().createFromWKT( + "VERTCRS[\"my depth\",VDATUM[\"my datum\"],CS[vertical,1]," + "AXIS[\"depth (D)\",down,LENGTHUNIT[\"metre\",1]]]")); + ASSERT_TRUE(vertCRSMetreDown != nullptr); + + auto vertCRSMetreDownFtUS = + nn_dynamic_pointer_cast(WKTParser().createFromWKT( + "VERTCRS[\"my depth (ftUS)\",VDATUM[\"my datum\"],CS[vertical,1]," + "AXIS[\"depth (D)\",down,LENGTHUNIT[\"US survey " + "foot\",0.304800609601219]]]")); + ASSERT_TRUE(vertCRSMetreDownFtUS != nullptr); + + { + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(vertCRSMetreUp), NN_CHECK_ASSERT(vertCRSMetreDown)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=axisswap +order=1,2,-3"); + } + + { + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(vertCRSMetreUp), + NN_CHECK_ASSERT(vertCRSMetreDownFtUS)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=affine +s33=-3.28083333333333"); + } } // --------------------------------------------------------------------------- @@ -7202,7 +7239,7 @@ TEST(operation, compoundCRS_to_geogCRS_with_vertical_unit_change) { ASSERT_EQ(listCompoundToGeog.size(), listCompoundMetreToGeog.size()); EXPECT_EQ(listCompoundToGeog[0]->nameStr(), - "Transformation from NAVD88 height (ftUS) to NAVD88 height + " + + "Inverse of NAVD88 height to NAVD88 height (ftUS) + " + listCompoundMetreToGeog[0]->nameStr()); EXPECT_EQ( listCompoundToGeog[0]->exportToPROJString( @@ -7267,7 +7304,7 @@ TEST( ASSERT_GE(listCompoundToGeog.size(), 1U); EXPECT_EQ(listCompoundToGeog[0]->nameStr(), - "Transformation from NAVD88 height (ftUS) to NAVD88 height + " + + "Inverse of NAVD88 height to NAVD88 height (ftUS) + " + listCompoundMetreToGeog[0]->nameStr()); EXPECT_EQ( listCompoundToGeog[0]->exportToPROJString( @@ -7286,6 +7323,134 @@ TEST( // --------------------------------------------------------------------------- +TEST(operation, compoundCRS_to_geogCRS_with_height_depth_reversal) { + 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 depth + auto srcObj = createFromUserInput("EPSG:6318+6357", + authFactory->databaseContext(), false); + auto src = nn_dynamic_pointer_cast(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(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 depth variant + ASSERT_EQ(listCompoundToGeog.size(), listCompoundMetreToGeog.size()); + + EXPECT_EQ(listCompoundToGeog[0]->nameStr(), + "Inverse of NAVD88 height to NAVD88 depth + " + + 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 +xy_out=rad " + "+step +proj=axisswap +order=1,2,-3")); + + // 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_height_depth_reversal) { + 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 depth (ftUS) + auto srcObj = createFromUserInput("EPSG:6318+6358", + authFactory->databaseContext(), false); + auto src = nn_dynamic_pointer_cast(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(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 depth (ftUS) variant + ASSERT_EQ(listCompoundToGeog.size(), listCompoundMetreToGeog.size()); + + EXPECT_EQ(listCompoundToGeog[0]->nameStr(), + "Inverse of NAVD88 height (ftUS) to NAVD88 depth (ftUS) + " + "Inverse of NAVD88 height to NAVD88 height (ftUS) + " + + 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 +xy_out=rad " + "+step +proj=axisswap +order=1,2,-3 " + "+step +proj=unitconvert +z_in=us-ft +z_out=m")); + + // Check reverse path + auto listGeogToCompound = + CoordinateOperationFactory::create()->createOperations(dst, nnSrc, + ctxt); + EXPECT_EQ(listGeogToCompound.size(), listCompoundToGeog.size()); +} + +// --------------------------------------------------------------------------- + TEST(operation, compoundCRS_from_WKT2_to_geogCRS_3D_context) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); -- cgit v1.2.3 From 6ce8ef30389480b3fabc3991bdf2b476d9435b60 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 30 Oct 2019 16:22:17 +0100 Subject: Document PROJJSON More could probably be written, but at least this can serve as a landing/reference page for other documents/specifications to point to. --- docs/source/news.rst | 2 +- docs/source/usage/index.rst | 2 + docs/source/usage/projjson.rst | 261 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 264 insertions(+), 1 deletion(-) create mode 100644 docs/source/usage/projjson.rst diff --git a/docs/source/news.rst b/docs/source/news.rst index 1fba43b9..bb8f63f9 100644 --- a/docs/source/news.rst +++ b/docs/source/news.rst @@ -10,7 +10,7 @@ News Updates ------- - * Introduced PROJJSON, a JSON encoding of WKT2 (`#1547 `_) + * Introduced :ref:`PROJJSON`, a JSON encoding of WKT2 (`#1547 `_) * Support CRS instantiation of OGC URN's (`#1505 `_) * Expose scope and remarks of database objects (`#1537 `_) diff --git a/docs/source/usage/index.rst b/docs/source/usage/index.rst index 823e8fe7..c31c6dce 100644 --- a/docs/source/usage/index.rst +++ b/docs/source/usage/index.rst @@ -17,3 +17,5 @@ command line applications or the C API that is a part of the software package. transformation environmentvars differences + projjson + diff --git a/docs/source/usage/projjson.rst b/docs/source/usage/projjson.rst new file mode 100644 index 00000000..09123711 --- /dev/null +++ b/docs/source/usage/projjson.rst @@ -0,0 +1,261 @@ +.. _projjson: + +================================================================================ +PROJJSON +================================================================================ + +PROJJSON is a JSON encoding of +`WKT2:2019 / ISO-19162:2019 `_, +which itself implements the model of +`OGC Topic 2: Referencing by coordinates `_. +Apart from the difference of encodings, the semantics is intented to be exactly +the same as WKT2:2019. + +PROJJSON is available as input and output of PROJ since PROJ 6.2. + +The current version is 0.1. + +Schema +------ + +A JSON schema of its grammar is available at +https://proj.org/schemas/v0.1/projjson.schema.json + +Content +------- + +The high level objects are: + +* Coordinate Reference Systems (CRS): + + - Common ones: + + + ``GeographicCRS`` + + ``GeodeticCRS`` + + ``ProjectedCRS`` + + ``CompoundCRS`` + + ``BoundCRS`` + + - More esoteric ones: + + + ``VerticalCRS`` + + ``EngineeringCRS`` + + ``TemporalCRS`` + + ``ParametricCRS`` + + ``DerivedGeographicCRS`` + + ``DerivedGeodeticCRS`` + + ``DerivedProjectedCRS`` + + ``DerivedVerticalCRS`` + + ``DerivedEngineeringCRS`` + + ``DerivedTemporalCRS`` + + ``DerivedParametricCRS`` + +* Coordinate operations: + + - ``Transformation`` + - ``Conversion`` + - ``ConcatenatedOperation`` + +* Others: + + - ``PrimeMeridian`` + - ``Ellipsoid`` + - ``Datum`` + - ``DatumEnsemble`` + +Examples +-------- + +GeographicCRS ++++++++++++++ + +The following invokation + +:: + + projinfo EPSG:4326 -o PROJJSON -q + +will output: + +.. code-block:: json + + { + "$schema": "https://proj.org/schemas/v0.1/projjson.schema.json", + "type": "GeographicCRS", + "name": "WGS 84", + "datum": { + "type": "GeodeticReferenceFrame", + "name": "World Geodetic System 1984", + "ellipsoid": { + "name": "WGS 84", + "semi_major_axis": 6378137, + "inverse_flattening": 298.257223563 + } + }, + "coordinate_system": { + "subtype": "ellipsoidal", + "axis": [ + { + "name": "Geodetic latitude", + "abbreviation": "Lat", + "direction": "north", + "unit": "degree" + }, + { + "name": "Geodetic longitude", + "abbreviation": "Lon", + "direction": "east", + "unit": "degree" + } + ] + }, + "area": "World", + "bbox": { + "south_latitude": -90, + "west_longitude": -180, + "north_latitude": 90, + "east_longitude": 180 + }, + "id": { + "authority": "EPSG", + "code": 4326 + } + } + + +ProjectedCRS +++++++++++++ + +The following invokation + +:: + + projinfo EPSG:32631 -o PROJJSON -q + +will output: + +.. code-block:: json + + { + "$schema": "https://proj.org/schemas/v0.1/projjson.schema.json", + "type": "ProjectedCRS", + "name": "WGS 84 / UTM zone 31N", + "base_crs": { + "name": "WGS 84", + "datum": { + "type": "GeodeticReferenceFrame", + "name": "World Geodetic System 1984", + "ellipsoid": { + "name": "WGS 84", + "semi_major_axis": 6378137, + "inverse_flattening": 298.257223563 + } + }, + "coordinate_system": { + "subtype": "ellipsoidal", + "axis": [ + { + "name": "Geodetic latitude", + "abbreviation": "Lat", + "direction": "north", + "unit": "degree" + }, + { + "name": "Geodetic longitude", + "abbreviation": "Lon", + "direction": "east", + "unit": "degree" + } + ] + }, + "id": { + "authority": "EPSG", + "code": 4326 + } + }, + "conversion": { + "name": "UTM zone 31N", + "method": { + "name": "Transverse Mercator", + "id": { + "authority": "EPSG", + "code": 9807 + } + }, + "parameters": [ + { + "name": "Latitude of natural origin", + "value": 0, + "unit": "degree", + "id": { + "authority": "EPSG", + "code": 8801 + } + }, + { + "name": "Longitude of natural origin", + "value": 3, + "unit": "degree", + "id": { + "authority": "EPSG", + "code": 8802 + } + }, + { + "name": "Scale factor at natural origin", + "value": 0.9996, + "unit": "unity", + "id": { + "authority": "EPSG", + "code": 8805 + } + }, + { + "name": "False easting", + "value": 500000, + "unit": "metre", + "id": { + "authority": "EPSG", + "code": 8806 + } + }, + { + "name": "False northing", + "value": 0, + "unit": "metre", + "id": { + "authority": "EPSG", + "code": 8807 + } + } + ] + }, + "coordinate_system": { + "subtype": "Cartesian", + "axis": [ + { + "name": "Easting", + "abbreviation": "E", + "direction": "east", + "unit": "metre" + }, + { + "name": "Northing", + "abbreviation": "N", + "direction": "north", + "unit": "metre" + } + ] + }, + "area": "World - N hemisphere - 0°E to 6°E - by country", + "bbox": { + "south_latitude": 0, + "west_longitude": 0, + "north_latitude": 84, + "east_longitude": 6 + }, + "id": { + "authority": "EPSG", + "code": 32631 + } + } -- cgit v1.2.3 From 1920c5cee79f258a69f1000124fd6fbde20f1490 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 30 Oct 2019 17:41:10 +0100 Subject: Doc: add a redirects extension, and make a projjson.html redirect at the top to its current location --- docs/source/_extensions/redirects.py | 39 ++++++++++++++++++++++++++++++++++++ docs/source/conf.py | 5 +++++ 2 files changed, 44 insertions(+) create mode 100644 docs/source/_extensions/redirects.py diff --git a/docs/source/_extensions/redirects.py b/docs/source/_extensions/redirects.py new file mode 100644 index 00000000..6a59b622 --- /dev/null +++ b/docs/source/_extensions/redirects.py @@ -0,0 +1,39 @@ +import os + +# https://tech.signavio.com/2017/managing-sphinx-redirects + + +template=""" + + + + +""" + + +def gather_redirects(): + output = {} + + output.update({ 'projjson.html' : 'usage/projjson.html' }) + return output + +from shutil import copyfile +# copy legacy redirects +def copy_legacy_redirects(app, docname): # Sphinx expects two arguments + if app.builder.name == 'html': + for key in app.config.redirect_files: + src = key + tgt = app.config.redirect_files[key] + html = template % (tgt, tgt) + with open(os.path.join(app.outdir, src), 'wb') as f: + f.write(html.encode('utf-8')) + f.close() + + + +def setup(app): + app.add_config_value('redirect_files', {}, 'html') + app.connect('build-finished', copy_legacy_redirects) + return { 'parallel_read_safe': False, 'parallel_write_safe': True } diff --git a/docs/source/conf.py b/docs/source/conf.py index 243c93b7..cd004777 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -20,6 +20,7 @@ import datetime # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. sys.path.insert(0, os.path.abspath('.')) +sys.path.insert(0, os.path.abspath('_extensions')) import bibstyle @@ -35,6 +36,7 @@ extensions = [ 'sphinx.ext.mathjax', 'sphinxcontrib.bibtex', 'breathe', + 'redirects', ] # Add any paths that contain templates here, relative to this directory. @@ -370,3 +372,6 @@ texinfo_documents = [ breathe_projects = { "cpp_stuff":"../build/xml/", } + +import redirects +redirect_files = redirects.gather_redirects() -- cgit v1.2.3 From 9095cf6fa351b5e6208cec811b86eb3d958c6f06 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 30 Oct 2019 16:56:41 +0100 Subject: createFromWkt(): be tolerant to missing scale_factor parameter (fixes #1700) This is invalid WKT, but GDAL 2.4 used to accept it and make a reasonable use of it... Currently we default it to 0 which is non sensical. Better use 1 as GDAL 2.4 did, and emit a warning. Other fix: proj_create_from_wkt() was documented to operate by default in non-strict validation mode, but it was actually in strict mode. So do as documented. --- src/iso19111/c_api.cpp | 20 ++++++++++--- src/iso19111/coordinateoperation.cpp | 33 +++++++++++++--------- src/iso19111/io.cpp | 55 +++++++++++++++++++++++++++++++++++- test/unit/test_c_api.cpp | 6 ++-- test/unit/test_io.cpp | 34 ++++++++++++++++++++++ 5 files changed, 128 insertions(+), 20 deletions(-) diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index f5f7ba55..337de313 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -498,7 +498,7 @@ template static PROJ_STRING_LIST to_string_list(T &&set) { * proj_string_list_destroy(). * @param out_grammar_errors Pointer to a PROJ_STRING_LIST object, or NULL. * If provided, *out_grammar_errors will contain a list of errors regarding the - * WKT grammaer. It must be freed with proj_string_list_destroy(). + * WKT grammar. It must be freed with proj_string_list_destroy(). * @return Object that must be unreferenced with proj_destroy(), or NULL in * case of error. */ @@ -522,6 +522,7 @@ PJ *proj_create_from_wkt(PJ_CONTEXT *ctx, const char *wkt, if (dbContext) { parser.attachDatabaseContext(NN_NO_CHECK(dbContext)); } + parser.setStrict(false); for (auto iter = options; iter && iter[0]; ++iter) { const char *value; if ((value = getOptionValue(*iter, "STRICT="))) { @@ -536,10 +537,19 @@ PJ *proj_create_from_wkt(PJ_CONTEXT *ctx, const char *wkt, auto obj = nn_dynamic_pointer_cast( parser.createFromWKT(wkt)); + std::vector warningsFromParsing; if (out_grammar_errors) { - auto warnings = parser.warningList(); - if (!warnings.empty()) { - *out_grammar_errors = to_string_list(warnings); + auto rawWarnings = parser.warningList(); + std::vector grammarWarnings; + for (const auto &msg : rawWarnings) { + if (msg.find("Default it to") != std::string::npos) { + warningsFromParsing.push_back(msg); + } else { + grammarWarnings.push_back(msg); + } + } + if (!grammarWarnings.empty()) { + *out_grammar_errors = to_string_list(grammarWarnings); } } @@ -548,6 +558,8 @@ PJ *proj_create_from_wkt(PJ_CONTEXT *ctx, const char *wkt, if (derivedCRS) { auto warnings = derivedCRS->derivingConversionRef()->validateParameters(); + warnings.insert(warnings.end(), warningsFromParsing.begin(), + warningsFromParsing.end()); if (!warnings.empty()) { *out_warnings = to_string_list(warnings); } diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index c7581642..5ac81aa1 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -6056,24 +6056,31 @@ void Conversion::_exportToPROJString( if (!param->proj_name) { continue; } - auto value = + const auto value = parameterValueMeasure(param->wkt2_name, param->epsg_code); + double valueConverted = 0; + if (value == nullMeasure) { + // Deal with missing values. In an ideal world, this would + // not happen + if (param->epsg_code == + EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN) { + valueConverted = 1.0; + } + } else if (param->unit_type == + common::UnitOfMeasure::Type::ANGULAR) { + valueConverted = + value.convertToUnit(common::UnitOfMeasure::DEGREE); + } else { + valueConverted = value.getSIValue(); + } + if (mapping->epsg_code == EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_1SP && strcmp(param->proj_name, "lat_1") == 0) { - formatter->addParam( - param->proj_name, - value.convertToUnit(common::UnitOfMeasure::DEGREE)); - formatter->addParam( - "lat_0", - value.convertToUnit(common::UnitOfMeasure::DEGREE)); - } else if (param->unit_type == - common::UnitOfMeasure::Type::ANGULAR) { - formatter->addParam( - param->proj_name, - value.convertToUnit(common::UnitOfMeasure::DEGREE)); + formatter->addParam(param->proj_name, valueConverted); + formatter->addParam("lat_0", valueConverted); } else { - formatter->addParam(param->proj_name, value.getSIValue()); + formatter->addParam(param->proj_name, valueConverted); } } diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index fc66b6c9..50ad5a87 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -3476,6 +3476,14 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( } propertiesMethod.set(IdentifiedObject::NAME_KEY, projectionName); + std::vector foundParameters; + if (mapping) { + size_t countParams = 0; + while (mapping->params[countParams] != nullptr) { + ++countParams; + } + foundParameters.resize(countParams); + } for (const auto &childNode : projCRSNode->GP()->children()) { if (ci_equal(childNode->GP()->value(), WKTConstants::PARAMETER)) { const auto &childNodeChildren = childNode->GP()->children(); @@ -3496,11 +3504,18 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( continue; } } - const auto *paramMapping = + auto *paramMapping = mapping ? getMappingFromWKT1(mapping, parameterName) : nullptr; if (mapping && mapping->epsg_code == EPSG_CODE_METHOD_MERCATOR_VARIANT_B && ci_equal(parameterName, "latitude_of_origin")) { + for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) { + if (mapping->params[idx]->epsg_code == + EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN) { + foundParameters[idx] = true; + break; + } + } parameterName = EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN; propertiesParameter.set( Identifier::CODE_KEY, @@ -3508,6 +3523,12 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( propertiesParameter.set(Identifier::CODESPACE_KEY, Identifier::EPSG); } else if (paramMapping) { + for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) { + if (mapping->params[idx] == paramMapping) { + foundParameters[idx] = true; + break; + } + } parameterName = paramMapping->wkt2_name; if (paramMapping->epsg_code != 0) { propertiesParameter.set(Identifier::CODE_KEY, @@ -3531,6 +3552,38 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( } } + // Add back important parameters that should normally be present, but + // are sometimes missing. Currently we only deal with Scale factor at + // natural origin. This is to avoid a default value of 0 to slip in later. + // But such WKT should be considered invalid. + if (mapping) { + for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) { + if (!foundParameters[idx] && + mapping->params[idx]->epsg_code == + EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN) { + + emitRecoverableWarning( + "The WKT string lacks a value " + "for " EPSG_NAME_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN + ". Default it to 1."); + + PropertyMap propertiesParameter; + propertiesParameter.set( + Identifier::CODE_KEY, + EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN); + propertiesParameter.set(Identifier::CODESPACE_KEY, + Identifier::EPSG); + propertiesParameter.set( + IdentifiedObject::NAME_KEY, + EPSG_NAME_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN); + parameters.push_back( + OperationParameter::create(propertiesParameter)); + values.push_back(ParameterValue::create( + Measure(1.0, UnitOfMeasure::SCALE_UNITY))); + } + } + } + return Conversion::create( PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"), propertiesMethod, parameters, values) diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index 8ac5a511..b8310ce5 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -231,7 +231,8 @@ TEST_F(CApi, proj_create_from_wkt) { " PRIMEM[\"Greenwich\",0],\n" " UNIT[\"degree\",0.0174532925199433]]", nullptr, nullptr, nullptr); - EXPECT_EQ(obj, nullptr); + ObjectKeeper keeper(obj); + EXPECT_NE(obj, nullptr); } { PROJ_STRING_LIST warningList = nullptr; @@ -244,7 +245,8 @@ TEST_F(CApi, proj_create_from_wkt) { " PRIMEM[\"Greenwich\",0],\n" " UNIT[\"degree\",0.0174532925199433]]", nullptr, &warningList, &errorList); - EXPECT_EQ(obj, nullptr); + ObjectKeeper keeper(obj); + EXPECT_NE(obj, nullptr); EXPECT_EQ(warningList, nullptr); proj_string_list_destroy(warningList); EXPECT_NE(errorList, nullptr); diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 074a59b9..8aff0908 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -1174,6 +1174,40 @@ TEST(wkt_parse, wkt1_Mercator_1SP_with_latitude_origin_0) { // --------------------------------------------------------------------------- +TEST(wkt_parse, wkt1_Mercator_1SP_without_scale_factor) { + // See https://github.com/OSGeo/PROJ/issues/1700 + auto wkt = "PROJCS[\"unnamed\",\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"unknown\",\n" + " SPHEROID[\"WGS84\",6378137,298.257223563]],\n" + " PRIMEM[\"Greenwich\",0],\n" + " UNIT[\"degree\",0.0174532925199433]],\n" + " PROJECTION[\"Mercator_1SP\"],\n" + " PARAMETER[\"central_meridian\",0],\n" + " PARAMETER[\"false_easting\",0],\n" + " PARAMETER[\"false_northing\",0],\n" + " UNIT[\"Meter\",1],\n" + " AXIS[\"Easting\",EAST],\n" + " AXIS[\"Northing\",NORTH]]"; + WKTParser parser; + parser.setStrict(false).attachDatabaseContext(DatabaseContext::create()); + auto obj = parser.createFromWKT(wkt); + EXPECT_TRUE(!parser.warningList().empty()); + + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + auto got_wkt = crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()); + EXPECT_TRUE(got_wkt.find("PARAMETER[\"scale_factor\",1]") != + std::string::npos) + << got_wkt; + EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +units=m " + "+no_defs +type=crs"); +} + +// --------------------------------------------------------------------------- + TEST(wkt_parse, wkt1_krovak_south_west) { auto wkt = "PROJCS[\"S-JTSK / Krovak\"," -- cgit v1.2.3 From 81e06f42c7552494bcb3f466b0b1317341187679 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 30 Oct 2019 18:20:17 +0100 Subject: Add a test to check we can use a VerticalCRS from its name only without the EPSG code --- test/unit/test_operation.cpp | 53 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index f7a86eb4..d61c74c0 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -7202,6 +7202,59 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) { // --------------------------------------------------------------------------- +TEST(operation, compoundCRS_from_wkt_without_id_to_geogCRS) { + 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); + auto wkt = + "COMPOUNDCRS[\"NAD83(2011) + NAVD88 height\",\n" + " GEOGCRS[\"NAD83(2011)\",\n" + " DATUM[\"NAD83 (National Spatial Reference System 2011)\",\n" + " ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n" + " LENGTHUNIT[\"metre\",1]]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " CS[ellipsoidal,2],\n" + " AXIS[\"geodetic latitude (Lat)\",north,\n" + " ORDER[1],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " AXIS[\"geodetic longitude (Lon)\",east,\n" + " ORDER[2],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]]],\n" + " VERTCRS[\"NAVD88 height\",\n" + " VDATUM[\"North American Vertical Datum 1988\"],\n" + " CS[vertical,1],\n" + " AXIS[\"gravity-related height (H)\",up,\n" + " LENGTHUNIT[\"metre\",1]]]]"; + auto srcObj = + createFromUserInput(wkt, authFactory->databaseContext(), false); + auto src = nn_dynamic_pointer_cast(srcObj); + ASSERT_TRUE(src != nullptr); + auto dst = + authFactory->createCoordinateReferenceSystem("6319"); // NAD83(2011) + + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(src), dst, ctxt); + // NAD83(2011) + NAVD88 height + auto srcRefObj = createFromUserInput("EPSG:6318+5703", + authFactory->databaseContext(), false); + auto srcRef = nn_dynamic_pointer_cast(srcRefObj); + ASSERT_TRUE(srcRef != nullptr); + ASSERT_TRUE( + src->isEquivalentTo(srcRef.get(), IComparable::Criterion::EQUIVALENT)); + auto listRef = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcRef), dst, ctxt); + + EXPECT_EQ(list.size(), listRef.size()); +} + +// --------------------------------------------------------------------------- + TEST(operation, compoundCRS_to_geogCRS_with_vertical_unit_change) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); -- cgit v1.2.3