Skip to content

Commit

Permalink
createOperations: make it work when transforming from/to a CompoundCR…
Browse files Browse the repository at this point in the history 
…S with a DerivedVerticalCRS with ellipsoidal height

Such CompoundCRS is strictly speaking not OGC Topic 2 conformant, but
can make sense from a practical point of view.

Fixes OSGeo#4175
  • Loading branch information
@rouault
rouault committed Jun 18, 2024
1 parent 44e45b9 commit ae77e39
Show file tree
Hide file tree
Showing 2 changed files with 270 additions and 7 deletions.
94 changes: 87 additions & 7 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 @@ -4739,8 +4747,33 @@ void CoordinateOperationFactory::Private::createOperationsDerivedTo(
res.emplace_back(opFirst);
return;
}

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

// Optimization to remove a no-op "Transformation 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 transf =
dynamic_cast<const Transformation *>(opsSecond.front().get());
if (transf &&
transf->nameStr() ==
buildTransfName(targetCRS->nameStr(), targetCRS->nameStr()) &&
transf->method()->getEPSGCode() ==
EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT &&
transf->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 @@ -5286,7 +5319,7 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeog(

void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark(
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,12 +5353,38 @@ 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(buildTransfName(!bIsSameDatum || factor != 1.0
? sourceCRS->nameStr()
: targetCRS->nameStr(),
targetCRS->nameStr()));

if (!bIsSameDatum) {
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)
Expand All @@ -5334,7 +5393,9 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark(
auto conv = Transformation::createChangeVerticalUnit(
map, sourceCRS, targetCRS,
common::Scale(heightDepthReversal ? -factor : factor), {});
conv->setHasBallparkTransformation(true);
if (!bIsSameDatum) {
conv->setHasBallparkTransformation(true);
}
res.push_back(conv);
}

Expand Down Expand Up @@ -5860,6 +5921,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
183 changes: 183 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5060,6 +5060,88 @@ TEST(operation, compoundCRS_with_vert_bound_to_bound_geog3D) {

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

TEST(operation,
compoundCRS_with_derivedVerticalCRS_ellipsoidal_height_to_geog) {
// Test scenario of https://github.com/OSGeo/PROJ/issues/4175
// Note that the CompoundCRS with a DerivedVerticalCRS with a "Ellipsoid"
// datum is an extension, not OGC Topic 2 compliant.

auto wkt =
"COMPOUNDCRS[\"UTM30 with vertical offset\",\n"
" PROJCRS[\"WGS 84 / UTM zone 30N\",\n"
" BASEGEOGCRS[\"WGS 84\",\n"
" ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n"
" MEMBER[\"World Geodetic System 1984 (Transit)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G730)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G873)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G1150)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G1674)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G1762)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G2139)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G2296)\"],\n"
" ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ENSEMBLEACCURACY[2.0]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" ID[\"EPSG\",4326]],\n"
" CONVERSION[\"UTM zone 30N\",\n"
" METHOD[\"Transverse Mercator\",\n"
" ID[\"EPSG\",9807]],\n"
" PARAMETER[\"Latitude of natural origin\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8801]],\n"
" PARAMETER[\"Longitude of natural origin\",-3,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8802]],\n"
" PARAMETER[\"Scale factor at natural origin\",0.9996,\n"
" SCALEUNIT[\"unity\",1],\n"
" ID[\"EPSG\",8805]],\n"
" PARAMETER[\"False easting\",500000,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8806]],\n"
" PARAMETER[\"False northing\",0,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8807]]],\n"
" CS[Cartesian,2],\n"
" AXIS[\"(E)\",east,\n"
" ORDER[1],\n"
" LENGTHUNIT[\"metre\",1]],\n"
" AXIS[\"(N)\",north,\n"
" ORDER[2],\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ID[\"EPSG\",32630]],\n"
" VERTCRS[\"Derived verticalCRS\",\n"
" BASEVERTCRS[\"Ellipsoid (metre)\",\n"
" VDATUM[\"Ellipsoid\"]],\n"
" DERIVINGCONVERSION[\"Conv Vertical Offset\",\n"
" METHOD[\"Vertical Offset\",\n"
" ID[\"EPSG\",9616]],\n"
" PARAMETER[\"Vertical Offset\",42.3,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8603]]],\n"
" CS[vertical,1],\n"
" AXIS[\"gravity-related height (H)\",up,\n"
" LENGTHUNIT[\"metre\",1,\n"
" ID[\"EPSG\",9001]]]]]";
auto objDst = WKTParser().createFromWKT(wkt);
auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
ASSERT_TRUE(dst != nullptr);

auto op = CoordinateOperationFactory::create()->createOperation(
GeographicCRS::EPSG_4979, NN_NO_CHECK(dst));
ASSERT_TRUE(op != nullptr);
EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +proj=geogoffset +dh=42.3 "
"+step +proj=utm +zone=30 +ellps=WGS84");
EXPECT_STREQ(op->nameStr().c_str(), "Conv Vertical Offset + UTM zone 30N");
}

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

TEST(operation, geocent_to_compoundCRS) {
auto objSrc = PROJStringParser().createFromPROJString(
"+proj=geocent +datum=WGS84 +units=m +type=crs");
Expand Down Expand Up @@ -6325,6 +6407,105 @@ TEST(operation,

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

TEST(operation,
compoundCRS_to_compoundCRS_with_derivedVerticalCRS_ellipsoidal_height) {
// Test scenario of https://github.com/OSGeo/PROJ/issues/4175
// Note that the CompoundCRS with a DerivedVerticalCRS with a "Ellipsoid"
// datum is an extension, not OGC Topic 2 compliant.

auto dbContext = DatabaseContext::create();
auto authFactory = AuthorityFactory::create(dbContext, std::string());
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setGridAvailabilityUse(
CoordinateOperationContext::GridAvailabilityUse::
IGNORE_GRID_AVAILABILITY);

// WGS84 + EGM96
auto objSrc = createFromUserInput("EPSG:4326+5773", dbContext);
auto src = nn_dynamic_pointer_cast<CompoundCRS>(objSrc);
ASSERT_TRUE(src != nullptr);

//
auto wkt =
"COMPOUNDCRS[\"UTM30 with vertical offset\",\n"
" PROJCRS[\"WGS 84 / UTM zone 30N\",\n"
" BASEGEOGCRS[\"WGS 84\",\n"
" ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n"
" MEMBER[\"World Geodetic System 1984 (Transit)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G730)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G873)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G1150)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G1674)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G1762)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G2139)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G2296)\"],\n"
" ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ENSEMBLEACCURACY[2.0]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" ID[\"EPSG\",4326]],\n"
" CONVERSION[\"UTM zone 30N\",\n"
" METHOD[\"Transverse Mercator\",\n"
" ID[\"EPSG\",9807]],\n"
" PARAMETER[\"Latitude of natural origin\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8801]],\n"
" PARAMETER[\"Longitude of natural origin\",-3,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8802]],\n"
" PARAMETER[\"Scale factor at natural origin\",0.9996,\n"
" SCALEUNIT[\"unity\",1],\n"
" ID[\"EPSG\",8805]],\n"
" PARAMETER[\"False easting\",500000,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8806]],\n"
" PARAMETER[\"False northing\",0,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8807]]],\n"
" CS[Cartesian,2],\n"
" AXIS[\"(E)\",east,\n"
" ORDER[1],\n"
" LENGTHUNIT[\"metre\",1]],\n"
" AXIS[\"(N)\",north,\n"
" ORDER[2],\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ID[\"EPSG\",32630]],\n"
" VERTCRS[\"Derived verticalCRS\",\n"
" BASEVERTCRS[\"Ellipsoid (metre)\",\n"
" VDATUM[\"Ellipsoid\"]],\n"
" DERIVINGCONVERSION[\"Conv Vertical Offset\",\n"
" METHOD[\"Vertical Offset\",\n"
" ID[\"EPSG\",9616]],\n"
" PARAMETER[\"Vertical Offset\",42.3,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8603]]],\n"
" CS[vertical,1],\n"
" AXIS[\"gravity-related height (H)\",up,\n"
" LENGTHUNIT[\"metre\",1,\n"
" ID[\"EPSG\",9001]]]]]";
auto objDst =
WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
ASSERT_TRUE(dst != nullptr);

auto list = CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt);
ASSERT_GE(list.size(), 1U);
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1 "
"+step +proj=geogoffset +dh=42.3 "
"+step +proj=utm +zone=30 +ellps=WGS84");
EXPECT_STREQ(list[0]->nameStr().c_str(),
"Inverse of WGS 84 to EGM96 height (1) + Conv Vertical Offset "
"+ UTM zone 30N");
}

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

TEST(operation, vertCRS_to_vertCRS) {

auto vertcrs_m_obj = PROJStringParser().createFromPROJString("+vunits=m");
Expand Down Expand Up @@ -9985,6 +10166,8 @@ TEST(operation, createOperation_derived_projected_crs) {
"+step +proj=axisswap +order=2,1");
}

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

TEST(operation,
geogCRS_to_compoundCRS_with_boundVerticalCRS_and_derivedProjected) {
auto wkt =
Expand Down

0 comments on commit ae77e39

Please sign in to comment.