CompoundCRS::create(): reject combinations of components not allowed by ISO 19111

5 files changed, 169 insertions, 39 deletions

diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp
index 3b760099..30134cb4 100644
--- a/include/proj/crs.hpp
+++ b/include/proj/crs.hpp
@@ -834,6 +834,22 @@ class PROJ_GCC_DLL ParametricCRS : virtual public SingleCRS {
 
 // ---------------------------------------------------------------------------
 
+/** \brief Exception thrown when attempting to create an invalid compound CRS
+ */
+class PROJ_GCC_DLL InvalidCompoundCRSException : public util::Exception {
+  public:
+    //! @cond Doxygen_Suppress
+    PROJ_INTERNAL explicit InvalidCompoundCRSException(const char *message);
+    PROJ_INTERNAL explicit InvalidCompoundCRSException(
+        const std::string &message);
+    PROJ_DLL
+    InvalidCompoundCRSException(const InvalidCompoundCRSException &other);
+    PROJ_DLL ~InvalidCompoundCRSException() override;
+    //! @endcond
+};
+
+// ---------------------------------------------------------------------------
+
 class CompoundCRS;
 /** Shared pointer of CompoundCRS */
 using CompoundCRSPtr = std::shared_ptr<CompoundCRS>;
@@ -873,7 +889,8 @@ class PROJ_GCC_DLL CompoundCRS final : public CRS,
 
     PROJ_DLL static CompoundCRSNNPtr
     create(const util::PropertyMap &properties,
-           const std::vector<CRSNNPtr> &components);
+           const std::vector<CRSNNPtr>
+               &components); // throw InvalidCompoundCRSException
 
   protected:
     // relaxed: standard say SingleCRSNNPtr
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index 50be8d8a..674abfed 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -156,6 +156,8 @@ osgeo::proj::crs::GeographicCRS::create(osgeo::proj::util::PropertyMap const&, s
 osgeo::proj::crs::GeographicCRS::demoteTo2D(std::string const&, std::shared_ptr<osgeo::proj::io::DatabaseContext> const&) const
 osgeo::proj::crs::GeographicCRS::~GeographicCRS()
 osgeo::proj::crs::GeographicCRS::is2DPartOf3D(dropbox::oxygen::nn<osgeo::proj::crs::GeographicCRS const*>) const
+osgeo::proj::crs::InvalidCompoundCRSException::~InvalidCompoundCRSException()
+osgeo::proj::crs::InvalidCompoundCRSException::InvalidCompoundCRSException(osgeo::proj::crs::InvalidCompoundCRSException const&)
 osgeo::proj::crs::ParametricCRS::coordinateSystem() const
 osgeo::proj::crs::ParametricCRS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::datum::ParametricDatum> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::ParametricCS> > const&)
 osgeo::proj::crs::ParametricCRS::datum() const
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 67b0bb00..dbb68f0b 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -3915,6 +3915,28 @@ ProjectedCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+InvalidCompoundCRSException::InvalidCompoundCRSException(const char *message)
+    : Exception(message) {}
+
+// ---------------------------------------------------------------------------
+
+InvalidCompoundCRSException::InvalidCompoundCRSException(
+    const std::string &message)
+    : Exception(message) {}
+
+// ---------------------------------------------------------------------------
+
+InvalidCompoundCRSException::~InvalidCompoundCRSException() = default;
+
+// ---------------------------------------------------------------------------
+
+InvalidCompoundCRSException::InvalidCompoundCRSException(
+    const InvalidCompoundCRSException &) = default;
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 struct CompoundCRS::Private {
     std::vector<CRSNNPtr> components_{};
 };
@@ -3965,9 +3987,67 @@ CompoundCRS::componentReferenceSystems() PROJ_PURE_DEFN {
  * At minimum the name should be defined.
  * @param components the component CRS of the CompoundCRS.
  * @return new CompoundCRS.
+ * @throw InvalidCompoundCRSException
  */
 CompoundCRSNNPtr CompoundCRS::create(const util::PropertyMap &properties,
                                      const std::vector<CRSNNPtr> &components) {
+
+    if (components.size() < 2) {
+        throw InvalidCompoundCRSException(
+            "compound CRS should have at least 2 components");
+    }
+
+    auto comp0 = components[0].get();
+    auto comp0Bound = dynamic_cast<const BoundCRS *>(comp0);
+    if (comp0Bound) {
+        comp0 = comp0Bound->baseCRS().get();
+    }
+    auto comp0Geog = dynamic_cast<const GeographicCRS *>(comp0);
+    auto comp0Proj = dynamic_cast<const ProjectedCRS *>(comp0);
+    auto comp0Eng = dynamic_cast<const EngineeringCRS *>(comp0);
+
+    auto comp1 = components[1].get();
+    auto comp1Bound = dynamic_cast<const BoundCRS *>(comp1);
+    if (comp1Bound) {
+        comp1 = comp1Bound->baseCRS().get();
+    }
+    auto comp1Vert = dynamic_cast<const VerticalCRS *>(comp1);
+    auto comp1Eng = dynamic_cast<const EngineeringCRS *>(comp1);
+    // Loose validation based on
+    // http://docs.opengeospatial.org/as/18-005r4/18-005r4.html#34
+    bool ok = false;
+    if ((comp0Geog && comp0Geog->coordinateSystem()->axisList().size() == 2 &&
+         (comp1Vert ||
+          (comp1Eng &&
+           comp1Eng->coordinateSystem()->axisList().size() == 1))) ||
+        (comp0Proj && comp0Proj->coordinateSystem()->axisList().size() == 2 &&
+         (comp1Vert ||
+          (comp1Eng &&
+           comp1Eng->coordinateSystem()->axisList().size() == 1))) ||
+        (comp0Eng && comp0Eng->coordinateSystem()->axisList().size() <= 2 &&
+         comp1Vert)) {
+        // Spatial compound coordinate reference system
+        ok = true;
+    } else {
+        bool isComp0Spatial = comp0Geog || comp0Proj || comp0Eng ||
+                              dynamic_cast<const GeodeticCRS *>(comp0) ||
+                              dynamic_cast<const VerticalCRS *>(comp0);
+        if (isComp0Spatial && dynamic_cast<const TemporalCRS *>(comp1)) {
+            // Spatio-temporal compound coordinate reference system
+            ok = true;
+        } else if (isComp0Spatial &&
+                   dynamic_cast<const ParametricCRS *>(comp1)) {
+            // Spatio-parametric compound coordinate reference system
+            ok = true;
+        }
+    }
+    if (!ok) {
+        throw InvalidCompoundCRSException(
+            "components of the compound CRS do not belong to one of the "
+            "allowed combinations of "
+            "http://docs.opengeospatial.org/as/18-005r4/18-005r4.html#34");
+    }
+
     auto compoundCRS(CompoundCRS::nn_make_shared<CompoundCRS>(components));
     compoundCRS->assignSelf(compoundCRS);
     compoundCRS->setProperties(properties);
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 542d1226..c988080b 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -3331,6 +3331,46 @@ static CompoundCRSNNPtr createCompoundCRS() {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, compoundCRS_valid) {
+    // geographic 2D + vertical
+    CompoundCRS::create(
+        PropertyMap(),
+        std::vector<CRSNNPtr>{GeographicCRS::EPSG_4326, createVerticalCRS()});
+
+    // projected 2D + vertical
+    CompoundCRS::create(
+        PropertyMap(),
+        std::vector<CRSNNPtr>{createProjected(), createVerticalCRS()});
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(crs, compoundCRS_invalid) {
+    EXPECT_THROW(CompoundCRS::create(PropertyMap(), {}),
+                 InvalidCompoundCRSException);
+
+    // Only one component
+    EXPECT_THROW(CompoundCRS::create(PropertyMap(),
+                                     std::vector<CRSNNPtr>{createProjected()}),
+                 InvalidCompoundCRSException);
+
+    // Two geographic
+    EXPECT_THROW(
+        CompoundCRS::create(PropertyMap(),
+                            std::vector<CRSNNPtr>{GeographicCRS::EPSG_4326,
+                                                  GeographicCRS::EPSG_4326}),
+        InvalidCompoundCRSException);
+
+    // geographic 3D + vertical
+    EXPECT_THROW(
+        CompoundCRS::create(PropertyMap(),
+                            std::vector<CRSNNPtr>{GeographicCRS::EPSG_4979,
+                                                  createVerticalCRS()}),
+        InvalidCompoundCRSException);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, compoundCRS_as_WKT2) {
     auto crs = createCompoundCRS();
     auto expected =
@@ -3385,18 +3425,10 @@ TEST(crs, compoundCRS_isEquivalentTo) {
     EXPECT_TRUE(crs->isEquivalentTo(crs.get()));
     EXPECT_TRUE(crs->shallowClone()->isEquivalentTo(crs.get()));
     EXPECT_FALSE(crs->isEquivalentTo(createUnrelatedObject().get()));
-    auto compoundCRSOfProjCRS =
-        CompoundCRS::create(PropertyMap().set(IdentifiedObject::NAME_KEY, ""),
-                            std::vector<CRSNNPtr>{createProjected()});
-    auto emptyCompoundCRS =
-        CompoundCRS::create(PropertyMap().set(IdentifiedObject::NAME_KEY, ""),
-                            std::vector<CRSNNPtr>{});
-    EXPECT_FALSE(compoundCRSOfProjCRS->isEquivalentTo(emptyCompoundCRS.get()));
-    auto compoundCRSOfGeogCRS =
-        CompoundCRS::create(PropertyMap().set(IdentifiedObject::NAME_KEY, ""),
-                            std::vector<CRSNNPtr>{GeographicCRS::EPSG_4326});
-    EXPECT_FALSE(
-        compoundCRSOfProjCRS->isEquivalentTo(compoundCRSOfGeogCRS.get()));
+    auto otherCompoundCRS = CompoundCRS::create(
+        PropertyMap().set(IdentifiedObject::NAME_KEY, ""),
+        std::vector<CRSNNPtr>{GeographicCRS::EPSG_4326, createVerticalCRS()});
+    EXPECT_FALSE(crs->isEquivalentTo(otherCompoundCRS.get()));
 }
 
 // ---------------------------------------------------------------------------
@@ -4154,11 +4186,11 @@ TEST(crs, extractGeographicCRS) {
               GeographicCRS::EPSG_4326);
     EXPECT_EQ(createProjected()->extractGeographicCRS(),
               GeographicCRS::EPSG_4326);
-    EXPECT_EQ(
-        CompoundCRS::create(PropertyMap(),
-                            std::vector<CRSNNPtr>{GeographicCRS::EPSG_4326})
-            ->extractGeographicCRS(),
-        GeographicCRS::EPSG_4326);
+    EXPECT_EQ(CompoundCRS::create(
+                  PropertyMap(), std::vector<CRSNNPtr>{GeographicCRS::EPSG_4326,
+                                                       createVerticalCRS()})
+                  ->extractGeographicCRS(),
+              GeographicCRS::EPSG_4326);
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 63ed97e7..1af41fde 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -1142,16 +1142,32 @@ TEST(operation, transformation_inverse) {
 
 // ---------------------------------------------------------------------------
 
+static VerticalCRSNNPtr createVerticalCRS() {
+    PropertyMap propertiesVDatum;
+    propertiesVDatum.set(Identifier::CODESPACE_KEY, "EPSG")
+        .set(Identifier::CODE_KEY, 5101)
+        .set(IdentifiedObject::NAME_KEY, "Ordnance Datum Newlyn");
+    auto vdatum = VerticalReferenceFrame::create(propertiesVDatum);
+    PropertyMap propertiesCRS;
+    propertiesCRS.set(Identifier::CODESPACE_KEY, "EPSG")
+        .set(Identifier::CODE_KEY, 5701)
+        .set(IdentifiedObject::NAME_KEY, "ODN height");
+    return VerticalCRS::create(
+        propertiesCRS, vdatum,
+        VerticalCS::createGravityRelatedHeight(UnitOfMeasure::METRE));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, transformation_createTOWGS84) {
 
     EXPECT_THROW(Transformation::createTOWGS84(GeographicCRS::EPSG_4326,
                                                std::vector<double>()),
                  InvalidOperation);
 
-    auto crsIn = CompoundCRS::create(PropertyMap(), std::vector<CRSNNPtr>{});
-    EXPECT_THROW(
-        Transformation::createTOWGS84(crsIn, std::vector<double>(7, 0)),
-        InvalidOperation);
+    EXPECT_THROW(Transformation::createTOWGS84(createVerticalCRS(),
+                                               std::vector<double>(7, 0)),
+                 InvalidOperation);
 }
 
 // ---------------------------------------------------------------------------
@@ -6563,23 +6579,6 @@ TEST(operation, WGS84_G1762_to_compoundCRS_with_bound_vertCRS) {
 
 // ---------------------------------------------------------------------------
 
-static VerticalCRSNNPtr createVerticalCRS() {
-    PropertyMap propertiesVDatum;
-    propertiesVDatum.set(Identifier::CODESPACE_KEY, "EPSG")
-        .set(Identifier::CODE_KEY, 5101)
-        .set(IdentifiedObject::NAME_KEY, "Ordnance Datum Newlyn");
-    auto vdatum = VerticalReferenceFrame::create(propertiesVDatum);
-    PropertyMap propertiesCRS;
-    propertiesCRS.set(Identifier::CODESPACE_KEY, "EPSG")
-        .set(Identifier::CODE_KEY, 5701)
-        .set(IdentifiedObject::NAME_KEY, "ODN height");
-    return VerticalCRS::create(
-        propertiesCRS, vdatum,
-        VerticalCS::createGravityRelatedHeight(UnitOfMeasure::METRE));
-}
-
-// ---------------------------------------------------------------------------
-
 TEST(operation, compoundCRS_to_geogCRS) {
 
     auto compound = CompoundCRS::create(