diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/iso19111/operation/coordinateoperationfactory.cpp | 110 |
1 files changed, 96 insertions, 14 deletions
diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index b8a9bcdf..51262ff6 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -3285,21 +3285,62 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid( ENTER_FUNCTION(); - const auto useTransf = [&targetCRS, &context, + const auto useTransf = [&sourceCRS, &targetCRS, &context, vertDst](const CoordinateOperationNNPtr &op) { + + // If the source geographic CRS has a non-metre vertical unit, we need + // to create an intermediate and operation to do the vertical unit + // conversion from that vertical unit to the one of the geographic CRS + // of the source of the operation + const auto geogCRS = + dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get()); + assert(geogCRS); + const auto &srcAxisList = geogCRS->coordinateSystem()->axisList(); + CoordinateOperationPtr opPtr; + const auto opSourceCRSGeog = + dynamic_cast<const crs::GeographicCRS *>(op->sourceCRS().get()); + // I assume opSourceCRSGeog should always be null in practice... + if (opSourceCRSGeog && srcAxisList.size() == 3 && + srcAxisList[2]->unit().conversionToSI() != 1) { + const auto &authFactory = context.context->getAuthorityFactory(); + const auto dbContext = + authFactory ? authFactory->databaseContext().as_nullable() + : nullptr; + auto tmpCRSWithSrcZ = + opSourceCRSGeog->demoteTo2D(std::string(), dbContext) + ->promoteTo3D(std::string(), dbContext, srcAxisList[2]); + + std::vector<CoordinateOperationNNPtr> opsUnitConvert; + createOperationsGeogToGeog( + opsUnitConvert, tmpCRSWithSrcZ, NN_NO_CHECK(op->sourceCRS()), + context, + dynamic_cast<const crs::GeographicCRS *>(tmpCRSWithSrcZ.get()), + opSourceCRSGeog); + assert(opsUnitConvert.size() == 1); + opPtr = opsUnitConvert.front().as_nullable(); + } + + std::vector<CoordinateOperationNNPtr> ops; + if (opPtr) + ops.emplace_back(NN_NO_CHECK(opPtr)); + ops.emplace_back(op); + const auto targetOp = dynamic_cast<const crs::VerticalCRS *>(op->targetCRS().get()); assert(targetOp); if (targetOp->_isEquivalentTo( vertDst, util::IComparable::Criterion::EQUIVALENT)) { - return op; + auto ret = ConcatenatedOperation::createComputeMetadata( + ops, disallowEmptyIntersection); + return ret; } std::vector<CoordinateOperationNNPtr> tmp; createOperationsVertToVert(NN_NO_CHECK(op->targetCRS()), targetCRS, context, targetOp, vertDst, tmp); assert(!tmp.empty()); + ops.emplace_back(tmp.front()); auto ret = ConcatenatedOperation::createComputeMetadata( - {op, tmp.front()}, disallowEmptyIntersection); + ops, disallowEmptyIntersection); return ret; }; @@ -3331,10 +3372,16 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid( }; const auto &axis = vertDst->coordinateSystem()->axisList()[0]; + const auto &authFactory = context.context->getAuthorityFactory(); + const auto dbContext = + authFactory ? authFactory->databaseContext().as_nullable() + : nullptr; + const auto geogSrcCRS = dynamic_cast<crs::GeographicCRS *>(model->interpolationCRS().get()) ? NN_NO_CHECK(model->interpolationCRS()) - : sourceCRS; + : sourceCRS->demoteTo2D(std::string(), dbContext) + ->promoteTo3D(std::string(), dbContext); const auto vertCRSMetre = axis->unit() == common::UnitOfMeasure::METRE && axis->direction() == cs::AxisDirection::UP @@ -3356,7 +3403,6 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid( std::vector<metadata::PositionalAccuracyNNPtr> accuracies; const auto &modelAccuracies = model->coordinateOperationAccuracies(); if (modelAccuracies.empty()) { - const auto &authFactory = context.context->getAuthorityFactory(); if (authFactory) { const auto transformationsForGrid = io::DatabaseContext::getTransformationsForGridName( @@ -3476,8 +3522,8 @@ std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private:: std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private:: createOperationsGeogToVertWithAlternativeGeog( - const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS - const crs::CRSNNPtr &targetCRS, // vertical CRS + const crs::CRSNNPtr &sourceCRS, // geographic CRS + const crs::CRSNNPtr &targetCRS, // vertical CRS Private::Context &context) { ENTER_FUNCTION(); @@ -3501,10 +3547,39 @@ std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private:: // Generally EPSG has operations from GeogCrs to VertCRS auto ops = findOpsInRegistryDirectTo(targetCRS, context); + const auto geogCRS = + dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get()); + assert(geogCRS); + const auto &srcAxisList = geogCRS->coordinateSystem()->axisList(); for (const auto &op : ops) { - const auto tmpCRS = op->sourceCRS(); - if (tmpCRS && dynamic_cast<const crs::GeographicCRS *>(tmpCRS.get())) { - res.emplace_back(op); + const auto tmpCRS = + dynamic_cast<const crs::GeographicCRS *>(op->sourceCRS().get()); + if (tmpCRS) { + if (srcAxisList.size() == 3 && + srcAxisList[2]->unit().conversionToSI() != 1) { + + const auto &authFactory = + context.context->getAuthorityFactory(); + const auto dbContext = + authFactory->databaseContext().as_nullable(); + auto tmpCRSWithSrcZ = + tmpCRS->demoteTo2D(std::string(), dbContext) + ->promoteTo3D(std::string(), dbContext, srcAxisList[2]); + + std::vector<CoordinateOperationNNPtr> opsUnitConvert; + createOperationsGeogToGeog( + opsUnitConvert, tmpCRSWithSrcZ, + NN_NO_CHECK(op->sourceCRS()), context, + dynamic_cast<const crs::GeographicCRS *>( + tmpCRSWithSrcZ.get()), + tmpCRS); + assert(opsUnitConvert.size() == 1); + auto concat = ConcatenatedOperation::createComputeMetadata( + {opsUnitConvert.front(), op}, disallowEmptyIntersection); + res.emplace_back(concat); + } else { + res.emplace_back(op); + } } } @@ -4564,11 +4639,18 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( !srcGeogCRS->_isEquivalentTo( geogDst, util::IComparable::Criterion::EQUIVALENT) && !srcGeogCRS->is2DPartOf3D(NN_NO_CHECK(geogDst), dbContext)) { - auto verticalTransformsTmp = createOperations( - componentsSrc[1], + auto geogCRSTmp = NN_NO_CHECK(srcGeogCRS) - ->promoteTo3D(std::string(), dbContext), - context); + ->demoteTo2D(std::string(), dbContext) + ->promoteTo3D( + std::string(), dbContext, + geogDst->coordinateSystem()->axisList().size() == 3 + ? geogDst->coordinateSystem()->axisList()[2] + : cs::VerticalCS::createGravityRelatedHeight( + common::UnitOfMeasure::METRE) + ->axisList()[0]); + auto verticalTransformsTmp = + createOperations(componentsSrc[1], geogCRSTmp, context); bool foundRegisteredTransform = false; foundRegisteredTransformWithAllGridsAvailable = false; for (const auto &op : verticalTransformsTmp) { |