Skip to content

Commit

Permalink
Merge pull request #4176 from rouault/fix_4175
Browse files Browse the repository at this point in the history 
createOperations: make it work when transforming from/to a CompoundCRS with a DerivedVerticalCRS with ellipsoidal height
  • Loading branch information
@rouault
rouault authored Jun 19, 2024
2 parents 90e148d + 53be9fb commit 36815b4
Show file tree
Hide file tree
Showing 2 changed files with 299 additions and 21 deletions.
133 changes: 112 additions & 21 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
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

0 comments on commit 36815b4

Please sign in to comment.