WKT import: use 'EPSG code for Horizontal CRS' parameter to derive in…

…terpolation CRS

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<crs::CRS> sourceCRSIn,
                            std::weak_ptr<crs::CRS> 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

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<int>(values[i]->value().getSIValue());
                 try {

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) {

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,

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

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<CRS>(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<CRS>(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"