createFromUserInput(): support OGC URN to create projectedCRS, for example to instanciate a projected 3D CRS

4 files changed, 212 insertions, 2 deletions

diff --git a/data/sql/customizations.sql b/data/sql/customizations.sql
index b9fcac31..5e038de4 100644
--- a/data/sql/customizations.sql
+++ b/data/sql/customizations.sql
@@ -59,3 +59,9 @@ INSERT INTO "ellipsoid" VALUES('PROJ','WGS60','WGS 60',NULL,'PROJ','EARTH',63781
 -- Extra ellipsoids from IAU2000 dictionary (see https://github.com/USGS-Astrogeology/GDAL_scripts/blob/master/OGC_IAU2000_WKT_v2/naifcodes_radii_m_wAsteroids_IAU2000.csv)
 
 INSERT INTO "ellipsoid" VALUES('PROJ','EARTH2000','Earth2000',NULL,'PROJ','EARTH',6378140.0,'EPSG','9001',NULL,6356750.0,0);
+
+-- Coordinate system ENh for ProjectedCRS 3D. Should be removed once EPSG has such a coordinate system
+INSERT INTO "coordinate_system" VALUES('PROJ','ENh','Cartesian',3);
+INSERT INTO "axis" VALUES('PROJ','1','Easting','E','east','PROJ','ENh',1,'EPSG','9001');
+INSERT INTO "axis" VALUES('PROJ','2','Northing','N','north','PROJ','ENh',2,'EPSG','9001');
+INSERT INTO "axis" VALUES('PROJ','3','Ellipsoidal height','h','up','PROJ','ENh',2,'EPSG','9001');
diff --git a/docs/source/apps/projinfo.rst b/docs/source/apps/projinfo.rst
index 105b1e1b..dabcc0fe 100644
--- a/docs/source/apps/projinfo.rst
+++ b/docs/source/apps/projinfo.rst
@@ -38,6 +38,10 @@ Synopsis
     - a OGC URN combining references for compound coordinate reference systems
       (e.g "urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717" or custom abbreviated
       syntax "EPSG:2393+5717"),
+    - a OGC URN combining references for references for projected or derived CRSs
+      e.g. for Projected 3D CRS "UTM zone 31N / WGS 84 (3D)":
+      "urn:ogc:def:crs,crs:EPSG::4979,cs:PROJ::ENh,coordinateOperation:EPSG::16031"
+      (*added in 6.2*)
     - a OGC URN combining references for concatenated operations
       (e.g. "urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618")
 
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 7329758a..7b8e270b 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -4449,13 +4449,118 @@ static BaseObjectNNPtr createFromUserInput(const std::string &text,
         }
     }
 
-    // OGC 07-092r2: para 7.5.2
-    // URN combined references for compound coordinate reference systems
     if (starts_with(text, "urn:ogc:def:crs,")) {
         if (!dbContext) {
             throw ParsingException("no database context specified");
         }
         auto tokensComma = split(text, ',');
+        if (tokensComma.size() == 4 && starts_with(tokensComma[1], "crs:") &&
+            starts_with(tokensComma[2], "cs:") &&
+            starts_with(tokensComma[3], "coordinateOperation:")) {
+            // OGC 07-092r2: para 7.5.4
+            // URN combined references for projected or derived CRSs
+            const auto &crsPart = tokensComma[1];
+            const auto tokensCRS = split(crsPart, ':');
+            if (tokensCRS.size() != 4) {
+                throw ParsingException(
+                    concat("invalid crs component: ", crsPart));
+            }
+            auto factoryCRS =
+                AuthorityFactory::create(NN_NO_CHECK(dbContext), tokensCRS[1]);
+            auto baseCRS =
+                factoryCRS->createCoordinateReferenceSystem(tokensCRS[3], true);
+
+            const auto &csPart = tokensComma[2];
+            auto tokensCS = split(csPart, ':');
+            if (tokensCS.size() != 4) {
+                throw ParsingException(
+                    concat("invalid cs component: ", csPart));
+            }
+            auto factoryCS =
+                AuthorityFactory::create(NN_NO_CHECK(dbContext), tokensCS[1]);
+            auto cs = factoryCS->createCoordinateSystem(tokensCS[3]);
+
+            const auto &opPart = tokensComma[3];
+            auto tokensOp = split(opPart, ':');
+            if (tokensOp.size() != 4) {
+                throw ParsingException(
+                    concat("invalid coordinateOperation component: ", opPart));
+            }
+            auto factoryOp =
+                AuthorityFactory::create(NN_NO_CHECK(dbContext), tokensOp[1]);
+            auto op = factoryOp->createCoordinateOperation(tokensOp[3], true);
+
+            if (dynamic_cast<GeographicCRS *>(baseCRS.get()) &&
+                dynamic_cast<Conversion *>(op.get()) &&
+                dynamic_cast<CartesianCS *>(cs.get())) {
+                auto geogCRS = NN_NO_CHECK(
+                    util::nn_dynamic_pointer_cast<GeographicCRS>(baseCRS));
+                auto name = op->nameStr() + " / " + baseCRS->nameStr();
+                if (geogCRS->coordinateSystem()->axisList().size() == 3 &&
+                    baseCRS->nameStr().find("3D") == std::string::npos) {
+                    name += " (3D)";
+                }
+                return ProjectedCRS::create(
+                    util::PropertyMap().set(IdentifiedObject::NAME_KEY, name),
+                    geogCRS,
+                    NN_NO_CHECK(util::nn_dynamic_pointer_cast<Conversion>(op)),
+                    NN_NO_CHECK(
+                        util::nn_dynamic_pointer_cast<CartesianCS>(cs)));
+            } else if (dynamic_cast<GeodeticCRS *>(baseCRS.get()) &&
+                       !dynamic_cast<GeographicCRS *>(baseCRS.get()) &&
+                       dynamic_cast<Conversion *>(op.get()) &&
+                       dynamic_cast<CartesianCS *>(cs.get())) {
+                return DerivedGeodeticCRS::create(
+                    util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                            op->nameStr() + " / " +
+                                                baseCRS->nameStr()),
+                    NN_NO_CHECK(
+                        util::nn_dynamic_pointer_cast<GeodeticCRS>(baseCRS)),
+                    NN_NO_CHECK(util::nn_dynamic_pointer_cast<Conversion>(op)),
+                    NN_NO_CHECK(
+                        util::nn_dynamic_pointer_cast<CartesianCS>(cs)));
+            } else if (dynamic_cast<GeographicCRS *>(baseCRS.get()) &&
+                       dynamic_cast<Conversion *>(op.get()) &&
+                       dynamic_cast<EllipsoidalCS *>(cs.get())) {
+                return DerivedGeographicCRS::create(
+                    util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                            op->nameStr() + " / " +
+                                                baseCRS->nameStr()),
+                    NN_NO_CHECK(
+                        util::nn_dynamic_pointer_cast<GeodeticCRS>(baseCRS)),
+                    NN_NO_CHECK(util::nn_dynamic_pointer_cast<Conversion>(op)),
+                    NN_NO_CHECK(
+                        util::nn_dynamic_pointer_cast<EllipsoidalCS>(cs)));
+            } else if (dynamic_cast<ProjectedCRS *>(baseCRS.get()) &&
+                       dynamic_cast<Conversion *>(op.get())) {
+                return DerivedProjectedCRS::create(
+                    util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                            op->nameStr() + " / " +
+                                                baseCRS->nameStr()),
+                    NN_NO_CHECK(
+                        util::nn_dynamic_pointer_cast<ProjectedCRS>(baseCRS)),
+                    NN_NO_CHECK(util::nn_dynamic_pointer_cast<Conversion>(op)),
+                    cs);
+            } else if (dynamic_cast<VerticalCRS *>(baseCRS.get()) &&
+                       dynamic_cast<Conversion *>(op.get()) &&
+                       dynamic_cast<VerticalCS *>(cs.get())) {
+                return DerivedVerticalCRS::create(
+                    util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                            op->nameStr() + " / " +
+                                                baseCRS->nameStr()),
+                    NN_NO_CHECK(
+                        util::nn_dynamic_pointer_cast<VerticalCRS>(baseCRS)),
+                    NN_NO_CHECK(util::nn_dynamic_pointer_cast<Conversion>(op)),
+                    NN_NO_CHECK(util::nn_dynamic_pointer_cast<VerticalCS>(cs)));
+            } else {
+                throw ParsingException("unsupported combination of baseCRS, CS "
+                                       "and coordinateOperation for a "
+                                       "DerivedCRS");
+            }
+        }
+
+        // OGC 07-092r2: para 7.5.2
+        // URN combined references for compound coordinate reference systems
         std::vector<CRSNNPtr> components;
         std::string name;
         for (size_t i = 1; i < tokensComma.size(); i++) {
@@ -4622,6 +4727,11 @@ static BaseObjectNNPtr createFromUserInput(const std::string &text,
  *      e.g. "urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717"
  *      We also accept a custom abbreviated syntax EPSG:2393+5717
  * </li>
+ * <li> OGC URN combining references for references for projected or derived
+ * CRSs
+ *      e.g. for Projected 3D CRS "UTM zone 31N / WGS 84 (3D)"
+ *      "urn:ogc:def:crs,crs:EPSG::4979,cs:PROJ::ENh,coordinateOperation:EPSG::16031"
+ * </li>
  * <li> OGC URN combining references for concatenated operations
  *      e.g.
  * "urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618"</li>
@@ -4665,6 +4775,11 @@ BaseObjectNNPtr createFromUserInput(const std::string &text,
  *      e.g. "urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717"
  *      We also accept a custom abbreviated syntax EPSG:2393+5717
  * </li>
+ * <li> OGC URN combining references for references for projected or derived
+ * CRSs
+ *      e.g. for Projected 3D CRS "UTM zone 31N / WGS 84 (3D)"
+ *      "urn:ogc:def:crs,crs:EPSG::4979,cs:PROJ::ENh,coordinateOperation:EPSG::16031"
+ * </li>
  * <li> OGC URN combining references for concatenated operations
  *      e.g.
  * "urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618"</li>
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 294e765d..cea59460 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -8874,6 +8874,91 @@ TEST(io, createFromUserInput) {
         EXPECT_EQ(crs->nameStr(),
                   "KKJ / Finland Uniform Coordinate System + N60 height");
     }
+
+    {
+        auto obj = createFromUserInput("urn:ogc:def:crs,crs:EPSG::4979,"
+                                       "cs:PROJ::ENh,"
+                                       "coordinateOperation:EPSG::16031",
+                                       dbContext);
+        auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        EXPECT_EQ(crs->nameStr(), "UTM zone 31N / WGS 84 (3D)");
+        EXPECT_EQ(crs->baseCRS()->getEPSGCode(), 4979);
+        EXPECT_EQ(crs->coordinateSystem()->axisList().size(), 3U);
+        EXPECT_EQ(crs->derivingConversion()->getEPSGCode(), 16031);
+    }
+
+    EXPECT_THROW(createFromUserInput(
+                     "urn:ogc:def:crs,crs:EPSG::4979,"
+                     "cs:PROJ::ENh,"
+                     "coordinateOperation:EPSG::1024", // not a conversion
+                     dbContext),
+                 ParsingException);
+    EXPECT_THROW(createFromUserInput("urn:ogc:def:crs,crs:,"
+                                     "cs:PROJ::ENh,"
+                                     "coordinateOperation:EPSG::16031",
+                                     dbContext),
+                 ParsingException);
+    EXPECT_THROW(createFromUserInput("urn:ogc:def:crs,crs:EPSG::4979,"
+                                     "cs:,"
+                                     "coordinateOperation:EPSG::16031",
+                                     dbContext),
+                 ParsingException);
+    EXPECT_THROW(createFromUserInput("urn:ogc:def:crs,crs:EPSG::4979,"
+                                     "cs:PROJ::ENh,"
+                                     "coordinateOperation:",
+                                     dbContext),
+                 ParsingException);
+
+    {
+        // Completely non-sensical from a geodesic point of view...
+        auto obj = createFromUserInput("urn:ogc:def:crs,crs:EPSG::4978,"
+                                       "cs:EPSG::6500,"
+                                       "coordinateOperation:EPSG::16031",
+                                       dbContext);
+        auto crs = nn_dynamic_pointer_cast<DerivedGeodeticCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        EXPECT_EQ(crs->baseCRS()->getEPSGCode(), 4978);
+        EXPECT_EQ(crs->coordinateSystem()->getEPSGCode(), 6500);
+        EXPECT_EQ(crs->derivingConversion()->getEPSGCode(), 16031);
+    }
+    {
+        // Completely non-sensical from a geodesic point of view...
+        auto obj = createFromUserInput("urn:ogc:def:crs,crs:EPSG::4979,"
+                                       "cs:EPSG::6423,"
+                                       "coordinateOperation:EPSG::16031",
+                                       dbContext);
+        auto crs = nn_dynamic_pointer_cast<DerivedGeographicCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        EXPECT_EQ(crs->baseCRS()->getEPSGCode(), 4979);
+        EXPECT_EQ(crs->coordinateSystem()->getEPSGCode(), 6423);
+        EXPECT_EQ(crs->derivingConversion()->getEPSGCode(), 16031);
+    }
+    {
+        // Completely non-sensical from a geodesic point of view...
+        auto obj = createFromUserInput("urn:ogc:def:crs,crs:EPSG::32631,"
+                                       "cs:EPSG::4400,"
+                                       "coordinateOperation:EPSG::16031",
+                                       dbContext);
+        auto crs = nn_dynamic_pointer_cast<DerivedProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        EXPECT_EQ(crs->baseCRS()->getEPSGCode(), 32631);
+        EXPECT_EQ(crs->coordinateSystem()->getEPSGCode(), 4400);
+        EXPECT_EQ(crs->derivingConversion()->getEPSGCode(), 16031);
+    }
+    {
+        // Completely non-sensical from a geodesic point of view...
+        auto obj = createFromUserInput("urn:ogc:def:crs,crs:EPSG::3855,"
+                                       "cs:EPSG::6499,"
+                                       "coordinateOperation:EPSG::16031",
+                                       dbContext);
+        auto crs = nn_dynamic_pointer_cast<DerivedVerticalCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        EXPECT_EQ(crs->baseCRS()->getEPSGCode(), 3855);
+        EXPECT_EQ(crs->coordinateSystem()->getEPSGCode(), 6499);
+        EXPECT_EQ(crs->derivingConversion()->getEPSGCode(), 16031);
+    }
+
     {
         auto obj = createFromUserInput("urn:ogc:def:coordinateOperation,"
                                        "coordinateOperation:EPSG::3895,"