diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index b34f8468a3..b8673d056b 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) @@ -699,7 +707,7 @@ struct CoordinateOperationFactory::Private { const crs::GeographicCRS *geogDst, std::vector &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, @@ -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(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{horizTransform, @@ -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); } } @@ -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(derivedSrc->baseCRS().get()) && + !dynamic_cast(derivedSrc->baseCRS().get()) && + dynamic_cast(targetCRS.get()) && + !dynamic_cast(targetCRS.get())) { + auto conv = dynamic_cast(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( @@ -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 &res) { @@ -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); + } } // --------------------------------------------------------------------------- @@ -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(), targetCRS->promoteTo3D(std::string(), dbContext), diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 42c5a2190b..5c43b39af7 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -5060,6 +5060,91 @@ 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"); + EXPECT_FALSE(op->hasBallparkTransformation()); + ASSERT_EQ(op->coordinateOperationAccuracies().size(), 1U); + EXPECT_EQ(op->coordinateOperationAccuracies()[0]->value(), "0"); +} + +// --------------------------------------------------------------------------- + TEST(operation, geocent_to_compoundCRS) { auto objSrc = PROJStringParser().createFromPROJString( "+proj=geocent +datum=WGS84 +units=m +type=crs"); @@ -6325,6 +6410,106 @@ 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"); + EXPECT_FALSE(list[0]->hasBallparkTransformation()); +} + +// --------------------------------------------------------------------------- + TEST(operation, vertCRS_to_vertCRS) { auto vertcrs_m_obj = PROJStringParser().createFromPROJString("+vunits=m"); @@ -9985,6 +10170,8 @@ TEST(operation, createOperation_derived_projected_crs) { "+step +proj=axisswap +order=2,1"); } +// --------------------------------------------------------------------------- + TEST(operation, geogCRS_to_compoundCRS_with_boundVerticalCRS_and_derivedProjected) { auto wkt =