From 940186f7afce353befd67cb8303999c8fce66f5c Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 19 Mar 2022 12:20:21 +0100 Subject: Merge pull request #3123 from rouault/fix_compound_crs_nad83_csrs createOperations(): fix transformation involving CompoundCRS, ToWGS84 and PROJ4_GRIDS --- .../operation/coordinateoperationfactory.cpp | 22 +++++++++ test/unit/test_operationfactory.cpp | 55 ++++++++++++++++++++++ 2 files changed, 77 insertions(+) diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 20042f22..e8b1a371 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -5084,6 +5084,28 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS), context); + // e.g when doing COMPOUND_CRS[ + // NAD83(CRS)+TOWGS84[0,0,0], + // CGVD28 height + EXTENSION["PROJ4_GRIDS","HT2_0.gtx"] + // to NAD83(CRS) 3D + const auto boundSrc = + dynamic_cast(componentsSrc[0].get()); + if (boundSrc && + boundSrc->baseCRS()->isEquivalentTo( + targetCRS->demoteTo2D(std::string(), dbContext) + .get(), + util::IComparable::Criterion::EQUIVALENT) && + boundSrc->hubCRS()->isEquivalentTo( + interpolationGeogCRS + ->demoteTo2D(std::string(), dbContext) + .get(), + util::IComparable::Criterion::EQUIVALENT)) { + // Make sure to use the same horizontal transformation + // (likely a null shift) + interpToTargetOps = applyInverse(srcToInterpOps); + return; + } + // But do the interpolation CRS to targetCRS in 3D // to have proper ellipsoid height transformation. // We need to force the vertical axis of this 3D'ified diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 0c3ecae4..3476aa9a 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -5970,6 +5970,61 @@ TEST(operation, compoundCRS_of_vertCRS_with_geoid_model_to_geogCRS) { // --------------------------------------------------------------------------- +TEST(operation, + compoundCRS_of_horizCRS_with_TOWGS84_vertCRS_with_geoid_model_to_geogCRS) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + auto wkt = "COMPD_CS[\"NAD83(CSRS) + CGVD28 height - HT2_0\",\n" + " GEOGCS[\"NAD83(CSRS)\",\n" + " DATUM[\"NAD83_Canadian_Spatial_Reference_System\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6140\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4617\"]],\n" + " VERT_CS[\"CGVD28 height - HT2_0\",\n" + " VERT_DATUM[\"Canadian Geodetic Vertical Datum of " + "1928\",2005,\n" + " EXTENSION[\"PROJ4_GRIDS\",\"HT2_0.gtx\"],\n" + " AUTHORITY[\"EPSG\",\"5114\"]],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Gravity-related height\",UP],\n" + " AUTHORITY[\"EPSG\",\"5713\"]]]"; + auto srcObj = + createFromUserInput(wkt, authFactory->databaseContext(), false); + auto src = nn_dynamic_pointer_cast(srcObj); + ASSERT_TRUE(src != nullptr); + // NAD83(CSRS) 3D + auto dst = authFactory->createCoordinateReferenceSystem("4955"); + + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(src), dst, ctxt); + ASSERT_EQ(list.size(), 1U); + auto op_proj = + list[0]->exportToPROJString(PROJStringFormatter::create().get()); + EXPECT_EQ(op_proj, "+proj=pipeline " + "+step +proj=push +v_1 +v_2 " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=HT2_0.gtx +multiplier=1 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1 " + "+step +proj=pop +v_1 +v_2"); +} + +// --------------------------------------------------------------------------- + TEST(operation, compoundCRS_from_WKT2_to_geogCRS_3D_context) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); -- cgit v1.2.3