diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp index 09f14a57a7..96429b3dbe 100644 --- a/include/proj/coordinateoperation.hpp +++ b/include/proj/coordinateoperation.hpp @@ -177,12 +177,15 @@ class PROJ_GCC_DLL CoordinateOperation : public common::ObjectUsage, PROJ_FRIEND(io::AuthorityFactory); PROJ_FRIEND(CoordinateOperationFactory); PROJ_FRIEND(ConcatenatedOperation); + PROJ_FRIEND(io::WKTParser); PROJ_INTERNAL void setWeakSourceTargetCRS(std::weak_ptr sourceCRSIn, std::weak_ptr targetCRSIn); PROJ_INTERNAL void setCRSs(const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn, const crs::CRSPtr &interpolationCRSIn); + PROJ_INTERNAL void + setInterpolationCRS(const crs::CRSPtr &interpolationCRSIn); PROJ_INTERNAL void setCRSs(const CoordinateOperation *in, bool inverseSourceTarget); PROJ_INTERNAL diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index e3e10d6c2b..cf4119a66e 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -3334,7 +3334,12 @@ WKTParser::Private::buildConversion(const WKTNodeNNPtr &node, Conversion::create(invConvProps, invMethodProps, parameters, values) ->inverse())); } - return Conversion::create(convProps, methodProps, parameters, values); + auto conv = Conversion::create(convProps, methodProps, parameters, values); + auto interpolationCRS = + dealWithEPSGCodeForInterpolationCRSParameter(parameters, values); + if (interpolationCRS) + conv->setInterpolationCRS(interpolationCRS); + return conv; } // --------------------------------------------------------------------------- @@ -3346,8 +3351,13 @@ CRSPtr WKTParser::Private::dealWithEPSGCodeForInterpolationCRSParameter( // crs_epsg_code] into proper interpolation CRS if (dbContext_ != nullptr) { for (size_t i = 0; i < parameters.size(); ++i) { - if (parameters[i]->nameStr() == - EPSG_NAME_PARAMETER_EPSG_CODE_FOR_INTERPOLATION_CRS) { + const auto &l_name = parameters[i]->nameStr(); + const auto epsgCode = parameters[i]->getEPSGCode(); + if (l_name == EPSG_NAME_PARAMETER_EPSG_CODE_FOR_INTERPOLATION_CRS || + epsgCode == + EPSG_CODE_PARAMETER_EPSG_CODE_FOR_INTERPOLATION_CRS || + l_name == EPSG_NAME_PARAMETER_EPSG_CODE_FOR_HORIZONTAL_CRS || + epsgCode == EPSG_CODE_PARAMETER_EPSG_CODE_FOR_HORIZONTAL_CRS) { const int code = static_cast(values[i]->value().getSIValue()); try { diff --git a/src/iso19111/operation/parammappings.cpp b/src/iso19111/operation/parammappings.cpp index 7a9f59b273..b19881ed8a 100644 --- a/src/iso19111/operation/parammappings.cpp +++ b/src/iso19111/operation/parammappings.cpp @@ -1053,6 +1053,8 @@ const struct ParamNameCode paramNameCodes[] = { PARAM_NAME_CODE(GEOCENTRIC_TRANSLATION_FILE), PARAM_NAME_CODE(INCLINATION_IN_LATITUDE), PARAM_NAME_CODE(INCLINATION_IN_LONGITUDE), + PARAM_NAME_CODE(EPSG_CODE_FOR_HORIZONTAL_CRS), + PARAM_NAME_CODE(EPSG_CODE_FOR_INTERPOLATION_CRS), }; const ParamNameCode *getParamNameCodes(size_t &nElts) { diff --git a/src/iso19111/operation/singleoperation.cpp b/src/iso19111/operation/singleoperation.cpp index 3e079953b4..60e90f2f6b 100644 --- a/src/iso19111/operation/singleoperation.cpp +++ b/src/iso19111/operation/singleoperation.cpp @@ -233,6 +233,14 @@ void CoordinateOperation::setCRSs(const crs::CRSNNPtr &sourceCRSIn, d->targetCRSWeak_ = targetCRSIn.as_nullable(); d->interpolationCRS_ = interpolationCRSIn; } + +// --------------------------------------------------------------------------- + +void CoordinateOperation::setInterpolationCRS( + const crs::CRSPtr &interpolationCRSIn) { + d->interpolationCRS_ = interpolationCRSIn; +} + // --------------------------------------------------------------------------- void CoordinateOperation::setCRSs(const CoordinateOperation *in, diff --git a/src/proj_constants.h b/src/proj_constants.h index cbdccd8f57..432152349c 100644 --- a/src/proj_constants.h +++ b/src/proj_constants.h @@ -658,6 +658,9 @@ #define EPSG_NAME_PARAMETER_INCLINATION_IN_LONGITUDE "Inclination in longitude" #define EPSG_CODE_PARAMETER_INCLINATION_IN_LONGITUDE 8731 +#define EPSG_NAME_PARAMETER_EPSG_CODE_FOR_HORIZONTAL_CRS "EPSG code for Horizontal CRS" +#define EPSG_CODE_PARAMETER_EPSG_CODE_FOR_HORIZONTAL_CRS 1037 + /* ------------------------------------------------------------------------ */ #define EPSG_CODE_METHOD_AFFINE_PARAMETRIC_TRANSFORMATION 9624 diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 5018ba8bd9..3688c245b3 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -7381,6 +7381,95 @@ TEST(operation, compoundCRS_to_proj_string_with_non_metre_height) { // --------------------------------------------------------------------------- +TEST(operation, compoundCRS_with_derived_vertical_CRS) { + + auto dbContext = DatabaseContext::create(); + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + + // ETRS89 + EGM2008 height + auto objSrc = createFromUserInput("EPSG:4258+3855", dbContext, false); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto wkt = + "COMPOUNDCRS[\"WGS 84 + Custom Vertical\",\n" + " GEOGCRS[\"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" + " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n" + " LENGTHUNIT[\"metre\",1]],\n" + " ENSEMBLEACCURACY[2.0]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " CS[ellipsoidal,2],\n" + " AXIS[\"geodetic latitude (Lat)\",north,\n" + " ORDER[1],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " AXIS[\"geodetic longitude (Lon)\",east,\n" + " ORDER[2],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]]],\n" + " VERTCRS[\"Custom Vertical\",\n" + " BASEVERTCRS[\"EGM2008 height\",\n" + " VDATUM[\"EGM2008 geoid\"]],\n" + " DERIVINGCONVERSION[\"vertical offs. and slope\",\n" + " METHOD[\"Vertical Offset and Slope\",\n" + " ID[\"EPSG\",1046]],\n" + " PARAMETER[\"Ordinate 1 of evaluation point\",47,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433],\n" + " ID[\"EPSG\",8617]],\n" + " PARAMETER[\"Ordinate 2 of evaluation point\",8,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433],\n" + " ID[\"EPSG\",8618]],\n" + " PARAMETER[\"Vertical Offset\",-0.245,\n" + " LENGTHUNIT[\"metre\",1],\n" + " ID[\"EPSG\",8603]],\n" + " PARAMETER[\"Inclination in latitude\",-0.21,\n" + " ANGLEUNIT[\"arc-second\",4.84813681109536E-06],\n" + " ID[\"EPSG\",8730]],\n" + " PARAMETER[\"Inclination in longitude\",-0.032,\n" + " ANGLEUNIT[\"arc-second\",4.84813681109536E-06],\n" + " ID[\"EPSG\",8731]],\n" + " PARAMETER[\"EPSG code for Horizontal CRS\",4326,\n" + " UNIT[\"none\",1],\n" + " ID[\"EPSG\",1037]]],\n" + " CS[vertical,1],\n" + " AXIS[\"gravity-related height (H)\",up,\n" + " LENGTHUNIT[\"metre\",1]],\n" + " USAGE[\n" + " SCOPE[\"unknown\"],\n" + " AREA[\"World\"],\n" + " BBOX[-90,-180,90,180]]]]"; + auto objDst = createFromUserInput(wkt, dbContext, false); + auto dst = nn_dynamic_pointer_cast(objDst); + ASSERT_TRUE(dst != nullptr); + + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->nameStr(), + "ETRS89 to WGS 84 (1) + vertical offs. and slope"); + 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=vertoffset +lat_0=47 +lon_0=8 +dh=-0.245 +slope_lat=-0.21 " + "+slope_lon=-0.032 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + TEST(operation, compoundCRS_to_PROJJSON_with_non_metre_height) { auto srcPROJJSON = "{\n"