From ae77e39de5b34cb16a40698ff2a40a738787d2a8 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 18 Jun 2024 23:05:09 +0200 Subject: [PATCH] createOperations: make it work when transforming from/to a CompoundCRS 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 #4175 --- .../operation/coordinateoperationfactory.cpp | 94 ++++++++- test/unit/test_operationfactory.cpp | 183 ++++++++++++++++++ 2 files changed, 270 insertions(+), 7 deletions(-) diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index b34f8468a3..03f45f73c9 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -565,6 +565,14 @@ struct CoordinateOperationFactory::Private { std::list>> 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 geogCRSOfVertCRSStack{}; + Context(const metadata::ExtentPtr &extent1In, const metadata::ExtentPtr &extent2In, const CoordinateOperationContextNNPtr &contextIn) @@ -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(derivedSrc->baseCRS().get()) && + !dynamic_cast(derivedSrc->baseCRS().get()) && + dynamic_cast(targetCRS.get()) && + !dynamic_cast(targetCRS.get())) { + auto transf = + dynamic_cast(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( @@ -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 &res) { @@ -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) @@ -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); } @@ -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(), targetCRS->promoteTo3D(std::string(), dbContext), diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 42c5a2190b..398a1698f5 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -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(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"); @@ -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(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(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"); @@ -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 =