Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

createOperations: make it work when transforming from/to a CompoundCRS with a DerivedVerticalCRS with ellipsoidal height #4176

Merged
merged 1 commit into from
Jun 19, 2024
Merged

createOperations: make it work when transforming from/to a CompoundCRS with a DerivedVerticalCRS with ellipsoidal height #4176

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
133 changes: 112 additions & 21 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -565,6 +565,14 @@ struct CoordinateOperationFactory::Private {
std::list<std::pair<std::string, std::string>>>
cacheNameToCRS{};

// Normally there should be at most one element in this stack
// This is set when computing a CompoundCRS to GeogCRS operation to
// relate the VerticalCRS of the CompoundCRS to the geographicCRS of
// the horizontal component of the CompoundCRS. This is especially
// used if the VerticalCRS is a DerivedVerticalCRS using a
// (non-standard) "Ellipsoid" vertical datum
std::vector<crs::GeographicCRSNNPtr> geogCRSOfVertCRSStack{};

Context(const metadata::ExtentPtr &extent1In,
const metadata::ExtentPtr &extent2In,
const CoordinateOperationContextNNPtr &contextIn)
Expand Down Expand Up @@ -699,7 +707,7 @@ struct CoordinateOperationFactory::Private {
const crs::GeographicCRS *geogDst,
std::vector<CoordinateOperationNNPtr> &res);

static void createOperationsVertToGeogBallpark(
static void createOperationsVertToGeogSynthetized(
const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
Private::Context &context, const crs::VerticalCRS *vertSrc,
const crs::GeographicCRS *geogDst,
Expand Down Expand Up @@ -2584,10 +2592,14 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
if (!remarks.empty()) {
properties.set(common::IdentifiedObject::REMARKS_KEY, remarks);
}
return createPROJBased(
properties, exportable, sourceCRS, targetCRS, nullptr,
verticalTransform->coordinateOperationAccuracies(),
verticalTransform->hasBallparkTransformation());
auto accuracies = verticalTransform->coordinateOperationAccuracies();
if (accuracies.empty() &&
dynamic_cast<const Conversion *>(verticalTransform.get())) {
accuracies.emplace_back(metadata::PositionalAccuracy::create("0"));
}
return createPROJBased(properties, exportable, sourceCRS, targetCRS,
nullptr, accuracies,
verticalTransform->hasBallparkTransformation());
} else {
bool emptyIntersection = false;
auto ops = std::vector<CoordinateOperationNNPtr>{horizTransform,
Expand Down Expand Up @@ -3735,8 +3747,8 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
res = applyInverse(createOperationsGeogToVertFromGeoid(
targetCRS, sourceCRS, vertSrc, context));
if (!res.empty()) {
createOperationsVertToGeogBallpark(sourceCRS, targetCRS, context,
vertSrc, geogDst, res);
createOperationsVertToGeogSynthetized(sourceCRS, targetCRS, context,
vertSrc, geogDst, res);
}
}

Expand Down Expand Up @@ -4739,8 +4751,32 @@ void CoordinateOperationFactory::Private::createOperationsDerivedTo(
res.emplace_back(opFirst);
return;
}

auto opsSecond = createOperations(derivedSrc->baseCRS(), sourceEpoch,
targetCRS, targetEpoch, context);

// Optimization to remove a no-op "Conversion from WGS 84 to WGS 84"
// when transforming from EPSG:4979 to a CompoundCRS whose vertical CRS
// is a DerivedVerticalCRS of a datum with ellipsoid height.
if (opsSecond.size() == 1 &&
!opsSecond.front()->hasBallparkTransformation() &&
dynamic_cast<const crs::VerticalCRS *>(derivedSrc->baseCRS().get()) &&
!dynamic_cast<const crs::DerivedCRS *>(derivedSrc->baseCRS().get()) &&
dynamic_cast<const crs::GeographicCRS *>(targetCRS.get()) &&
!dynamic_cast<const crs::DerivedCRS *>(targetCRS.get())) {
auto conv = dynamic_cast<const Conversion *>(opsSecond.front().get());
if (conv &&
conv->nameStr() ==
buildConvName(targetCRS->nameStr(), targetCRS->nameStr()) &&
conv->method()->getEPSGCode() ==
EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT &&
conv->parameterValueNumericAsSI(
EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR) == 1.0) {
res.emplace_back(opFirst);
return;
}
}

for (const auto &opSecond : opsSecond) {
try {
res.emplace_back(ConcatenatedOperation::createComputeMetadata(
Expand Down Expand Up @@ -5278,15 +5314,15 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeog(
}
}

createOperationsVertToGeogBallpark(sourceCRS, targetCRS, context, vertSrc,
geogDst, res);
createOperationsVertToGeogSynthetized(sourceCRS, targetCRS, context,
vertSrc, geogDst, res);
}

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

void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark(
void CoordinateOperationFactory::Private::createOperationsVertToGeogSynthetized(
const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
Private::Context &, const crs::VerticalCRS *vertSrc,
Private::Context &context, const crs::VerticalCRS *vertSrc,
const crs::GeographicCRS *geogDst,
std::vector<CoordinateOperationNNPtr> &res) {

Expand Down Expand Up @@ -5320,22 +5356,58 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark(
sourceCRSExtent->_isEquivalentTo(
targetCRSExtent.get(), util::IComparable::Criterion::EQUIVALENT);

const auto &authFactory = context.context->getAuthorityFactory();
const auto dbContext =
authFactory ? authFactory->databaseContext().as_nullable() : nullptr;

const auto vertDatumName = vertSrc->datumNonNull(dbContext)->nameStr();
const auto geogDstDatumName = geogDst->datumNonNull(dbContext)->nameStr();
// We accept a vertical CRS whose datum name is the same datum name as the
// source geographic CRS, or whose datum name is "Ellipsoid" if it is part
// of a CompoundCRS whose horizontal CRS has a geodetic datum of the same
// datum name as the source geographic CRS, to mean an ellipsoidal height.
// This is against OGC Topic 2, and an extension needed for use case of
// https://github.com/OSGeo/PROJ/issues/4175
const bool bIsSameDatum = vertDatumName != "unknown" &&
(vertDatumName == geogDstDatumName ||
(vertDatumName == "Ellipsoid" &&
!context.geogCRSOfVertCRSStack.empty() &&
context.geogCRSOfVertCRSStack.back()
->datumNonNull(dbContext)
->nameStr() == geogDstDatumName));

std::string transfName;

if (bIsSameDatum) {
transfName = buildConvName(factor != 1.0 ? sourceCRS->nameStr()
: targetCRS->nameStr(),
targetCRS->nameStr());
} else {
transfName =
buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
transfName += " (";
transfName += BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT;
transfName += ')';
}

util::PropertyMap map;
std::string transfName(
buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()));
transfName += " (";
transfName += BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT;
transfName += ')';
map.set(common::IdentifiedObject::NAME_KEY, transfName)
.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
sameExtent ? NN_NO_CHECK(sourceCRSExtent)
: metadata::Extent::WORLD);

auto conv = Transformation::createChangeVerticalUnit(
map, sourceCRS, targetCRS,
common::Scale(heightDepthReversal ? -factor : factor), {});
conv->setHasBallparkTransformation(true);
res.push_back(conv);
if (bIsSameDatum) {
auto conv = Conversion::createChangeVerticalUnit(
map, common::Scale(heightDepthReversal ? -factor : factor));
conv->setCRSs(sourceCRS, targetCRS, nullptr);
res.push_back(conv);
} else {
auto transf = Transformation::createChangeVerticalUnit(
map, sourceCRS, targetCRS,
common::Scale(heightDepthReversal ? -factor : factor), {});
transf->setHasBallparkTransformation(true);
res.push_back(transf);
}
}

// ---------------------------------------------------------------------------
Expand Down Expand Up @@ -5860,6 +5932,25 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
};
SetSkipHorizontalTransform setSkipHorizontalTransform(context);

struct SetGeogCRSOfVertCRS {
Context &context;
const bool hasPushed;

explicit SetGeogCRSOfVertCRS(
Context &contextIn, const crs::GeographicCRSPtr &geogCRS)
: context(contextIn), hasPushed(geogCRS != nullptr) {
if (geogCRS)
context.geogCRSOfVertCRSStack.push_back(
NN_NO_CHECK(geogCRS));
}

~SetGeogCRSOfVertCRS() {
if (hasPushed)
context.geogCRSOfVertCRSStack.pop_back();
}
};
SetGeogCRSOfVertCRS setGeogCRSOfVertCRS(context, srcGeogCRS);

verticalTransforms = createOperations(
componentsSrc[1], util::optional<common::DataEpoch>(),
targetCRS->promoteTo3D(std::string(), dbContext),
Expand Down
Loading