Merge pull request #2389 from rouault/wkt_parse_implict_compound_crs_esri

WKT1_ESRI: fix import and export of CompoundCRS

3 files changed, 197 insertions, 32 deletions

diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 52b10119..edc8a71f 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -1578,6 +1578,49 @@ GeodeticCRS::create(const util::PropertyMap &properties,
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+
+// Try to format a Geographic/ProjectedCRS 3D CRS as a
+// GEOGCS[]/PROJCS[],VERTCS[...,DATUM[],...] if we find corresponding objects
+static bool exportAsESRIWktCompoundCRSWithEllipsoidalHeight(
+    const CRS *self, const GeodeticCRS *geodCRS, io::WKTFormatter *formatter) {
+    const auto &dbContext = formatter->databaseContext();
+    if (!dbContext) {
+        return false;
+    }
+    const auto l_datum = geodCRS->datumNonNull(formatter->databaseContext());
+    auto l_alias = dbContext->getAliasFromOfficialName(
+        l_datum->nameStr(), "geodetic_datum", "ESRI");
+    if (l_alias.empty()) {
+        return false;
+    }
+    auto authFactory =
+        io::AuthorityFactory::create(NN_NO_CHECK(dbContext), std::string());
+    auto list = authFactory->createObjectsFromName(
+        l_alias, {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME},
+        false /* approximate=false*/);
+    if (list.empty()) {
+        return false;
+    }
+    auto gdatum = util::nn_dynamic_pointer_cast<datum::Datum>(list.front());
+    if (gdatum == nullptr || gdatum->identifiers().empty()) {
+        return false;
+    }
+    const auto &gdatum_ids = gdatum->identifiers();
+    auto vertCRSList = authFactory->createVerticalCRSFromDatum(
+        "ESRI", "from_geogdatum_" + *gdatum_ids[0]->codeSpace() + '_' +
+                    gdatum_ids[0]->code());
+    if (vertCRSList.size() != 1) {
+        return false;
+    }
+    self->demoteTo2D(std::string(), dbContext)->_exportToWKT(formatter);
+    vertCRSList.front()->_exportToWKT(formatter);
+    return true;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
     const bool isGeographic =
@@ -1589,11 +1632,21 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
     auto l_name = nameStr();
     const auto &dbContext = formatter->databaseContext();
 
-    if (formatter->useESRIDialect()) {
-        if (axisList.size() != 2) {
+    if (!isWKT2 && formatter->useESRIDialect() && axisList.size() == 3) {
+        if (!isGeographic) {
             io::FormattingException::Throw(
-                "Only export of Geographic 2D CRS is supported in WKT1_ESRI");
+                "Geocentric CRS not supported in WKT1_ESRI");
+        }
+        // Try to format the Geographic 3D CRS as a GEOGCS[],VERTCS[...,DATUM[]]
+        // if we find corresponding objects
+        if (dbContext) {
+            if (exportAsESRIWktCompoundCRSWithEllipsoidalHeight(this, this,
+                                                                formatter)) {
+                return;
+            }
         }
+        io::FormattingException::Throw(
+            "Cannot export this Geographic 3D CRS in WKT1_ESRI");
     }
 
     if (!isWKT2 && formatter->isStrict() && isGeographic &&
@@ -2853,9 +2906,9 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
                          !identifiers().empty());
 
     auto l_name = nameStr();
+    const auto &dbContext = formatter->databaseContext();
     if (formatter->useESRIDialect()) {
         bool aliasFound = false;
-        const auto &dbContext = formatter->databaseContext();
         if (dbContext) {
             auto l_alias = dbContext->getAliasFromOfficialName(
                 l_name, "vertical_crs", "ESRI");
@@ -2870,7 +2923,34 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
     }
 
     formatter->addQuotedString(l_name);
-    exportDatumOrDatumEnsembleToWkt(formatter);
+
+    const auto l_datum = datum();
+    if (formatter->useESRIDialect() && l_datum &&
+        l_datum->getWKT1DatumType() == "2002") {
+        bool foundMatch = false;
+        if (dbContext) {
+            auto authFactory = io::AuthorityFactory::create(
+                NN_NO_CHECK(dbContext), std::string());
+            auto list = authFactory->createObjectsFromName(
+                l_datum->nameStr(),
+                {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME},
+                false /* approximate=false*/);
+            if (!list.empty()) {
+                auto gdatum =
+                    util::nn_dynamic_pointer_cast<datum::Datum>(list.front());
+                if (gdatum) {
+                    gdatum->_exportToWKT(formatter);
+                    foundMatch = true;
+                }
+            }
+        }
+        if (!foundMatch) {
+            // We should export a geodetic datum, but we cannot really do better
+            l_datum->_exportToWKT(formatter);
+        }
+    } else {
+        exportDatumOrDatumEnsembleToWkt(formatter);
+    }
     const auto &cs = SingleCRS::getPrivate()->coordinateSystem;
     const auto &axisList = cs->axisList();
 
@@ -3497,6 +3577,16 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
             }
         }
 
+        if (formatter->useESRIDialect() && dbContext) {
+            // Try to format the ProjecteD 3D CRS as a
+            // PROJCS[],VERTCS[...,DATUM[]]
+            // if we find corresponding objects
+            if (exportAsESRIWktCompoundCRSWithEllipsoidalHeight(
+                    this, baseCRS().as_nullable().get(), formatter)) {
+                return;
+            }
+        }
+
         if (!formatter->useESRIDialect() &&
             CRS::getPrivate()->allowNonConformantWKT1Export_) {
             formatter->startNode(io::WKTConstants::COMPD_CS, false);
@@ -4528,15 +4618,21 @@ CRSNNPtr CompoundCRS::createLax(const util::PropertyMap &properties,
 //! @cond Doxygen_Suppress
 void CompoundCRS::_exportToWKT(io::WKTFormatter *formatter) const {
     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
-    formatter->startNode(isWKT2 ? io::WKTConstants::COMPOUNDCRS
-                                : io::WKTConstants::COMPD_CS,
-                         !identifiers().empty());
-    formatter->addQuotedString(nameStr());
-    for (const auto &crs : componentReferenceSystems()) {
-        crs->_exportToWKT(formatter);
+    const auto &l_components = componentReferenceSystems();
+    if (!isWKT2 && formatter->useESRIDialect() && l_components.size() == 2) {
+        l_components[0]->_exportToWKT(formatter);
+        l_components[1]->_exportToWKT(formatter);
+    } else {
+        formatter->startNode(isWKT2 ? io::WKTConstants::COMPOUNDCRS
+                                    : io::WKTConstants::COMPD_CS,
+                             !identifiers().empty());
+        formatter->addQuotedString(nameStr());
+        for (const auto &crs : l_components) {
+            crs->_exportToWKT(formatter);
+        }
+        ObjectUsage::baseExportToWKT(formatter);
+        formatter->endNode();
     }
-    ObjectUsage::baseExportToWKT(formatter);
-    formatter->endNode();
 }
 //! @endcond
 
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index ef5c6e02..211eb586 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -2100,6 +2100,9 @@ AuthorityFactory::createVerticalDatum(const std::string &code) const {
         if (!publication_date.empty()) {
             props.set("PUBLICATION_DATE", publication_date);
         }
+        if (d->authority() == "ESRI" && starts_with(code, "from_geogdatum_")) {
+            props.set("VERT_DATUM_TYPE", "2002");
+        }
         auto anchor = util::optional<std::string>();
         if (frame_reference_epoch.empty()) {
             return datum::VerticalReferenceFrame::create(props, anchor);
@@ -5719,9 +5722,12 @@ AuthorityFactory::createObjectsFromNameEx(
         "SELECT table_name, auth_name, code, name, deprecated, is_alias "
         "FROM (");
 
-    const auto getTableAndTypeConstraints = [&allowedObjectTypes]() {
+    const auto getTableAndTypeConstraints = [&allowedObjectTypes,
+                                             &searchedName]() {
         typedef std::pair<std::string, std::string> TableType;
         std::list<TableType> res;
+        // Hide ESRI D_ vertical datums
+        const bool startsWithDUnderscore = starts_with(searchedName, "D_");
         if (allowedObjectTypes.empty()) {
             for (const auto &tableName :
                  {"prime_meridian", "ellipsoid", "geodetic_datum",
@@ -5729,7 +5735,10 @@ AuthorityFactory::createObjectsFromNameEx(
                   "vertical_crs", "compound_crs", "conversion",
                   "helmert_transformation", "grid_transformation",
                   "other_transformation", "concatenated_operation"}) {
-                res.emplace_back(TableType(tableName, std::string()));
+                if (!(startsWithDUnderscore &&
+                      strcmp(tableName, "vertical_datum") == 0)) {
+                    res.emplace_back(TableType(tableName, std::string()));
+                }
             }
         } else {
             for (const auto type : allowedObjectTypes) {
@@ -5744,8 +5753,10 @@ AuthorityFactory::createObjectsFromNameEx(
                 case ObjectType::DATUM:
                     res.emplace_back(
                         TableType("geodetic_datum", std::string()));
-                    res.emplace_back(
-                        TableType("vertical_datum", std::string()));
+                    if (!startsWithDUnderscore) {
+                        res.emplace_back(
+                            TableType("vertical_datum", std::string()));
+                    }
                     break;
                 case ObjectType::GEODETIC_REFERENCE_FRAME:
                     res.emplace_back(
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 0f4ffba0..b8e835d7 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -4145,11 +4145,19 @@ createBoundCRSSourceTransformationCRS(const crs::CRSPtr &sourceCRS,
 
 CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
     const auto *nodeP = node->GP();
-    auto &datumNode =
+    const auto &nodeValue = nodeP->value();
+    auto &vdatumNode =
         nodeP->lookForChild(WKTConstants::VDATUM, WKTConstants::VERT_DATUM,
                             WKTConstants::VERTICALDATUM, WKTConstants::VRF);
     auto &ensembleNode = nodeP->lookForChild(WKTConstants::ENSEMBLE);
-    if (isNull(datumNode) && isNull(ensembleNode)) {
+    // like in ESRI  VERTCS["WGS_1984",DATUM["D_WGS_1984",
+    //               SPHEROID["WGS_1984",6378137.0,298.257223563]],
+    //               PARAMETER["Vertical_Shift",0.0],
+    //               PARAMETER["Direction",1.0],UNIT["Meter",1.0]
+    auto &geogDatumNode = ci_equal(nodeValue, WKTConstants::VERTCS)
+                              ? nodeP->lookForChild(WKTConstants::DATUM)
+                              : null_node;
+    if (isNull(vdatumNode) && isNull(geogDatumNode) && isNull(ensembleNode)) {
         throw ParsingException("Missing VDATUM or ENSEMBLE node");
     }
 
@@ -4164,28 +4172,50 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
     }
 
     auto &dynamicNode = nodeP->lookForChild(WKTConstants::DYNAMIC);
-    auto datum =
-        !isNull(datumNode)
-            ? buildVerticalReferenceFrame(datumNode, dynamicNode).as_nullable()
-            : nullptr;
+    auto vdatum =
+        !isNull(geogDatumNode)
+            ? VerticalReferenceFrame::create(
+                  PropertyMap()
+                      .set(IdentifiedObject::NAME_KEY,
+                           buildGeodeticReferenceFrame(geogDatumNode,
+                                                       PrimeMeridian::GREENWICH,
+                                                       null_node)
+                               ->nameStr())
+                      .set("VERT_DATUM_TYPE", "2002"))
+                  .as_nullable()
+            : !isNull(vdatumNode)
+                  ? buildVerticalReferenceFrame(vdatumNode, dynamicNode)
+                        .as_nullable()
+                  : nullptr;
     auto datumEnsemble =
         !isNull(ensembleNode)
             ? buildDatumEnsemble(ensembleNode, nullptr, false).as_nullable()
             : nullptr;
 
     auto &csNode = nodeP->lookForChild(WKTConstants::CS_);
-    const auto &nodeValue = nodeP->value();
     if (isNull(csNode) && !ci_equal(nodeValue, WKTConstants::VERT_CS) &&
         !ci_equal(nodeValue, WKTConstants::VERTCS) &&
         !ci_equal(nodeValue, WKTConstants::BASEVERTCRS)) {
         ThrowMissing(WKTConstants::CS_);
     }
-    auto cs = buildCS(csNode, node, UnitOfMeasure::NONE);
-    auto verticalCS = nn_dynamic_pointer_cast<VerticalCS>(cs);
+    auto verticalCS = nn_dynamic_pointer_cast<VerticalCS>(
+        buildCS(csNode, node, UnitOfMeasure::NONE));
     if (!verticalCS) {
         ThrowNotExpectedCSType("vertical");
     }
 
+    if (vdatum && vdatum->getWKT1DatumType() == "2002" &&
+        &(verticalCS->axisList()[0]->direction()) == &(AxisDirection::UP)) {
+        verticalCS =
+            VerticalCS::create(
+                util::PropertyMap(),
+                CoordinateSystemAxis::create(
+                    util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                            "ellipsoidal height"),
+                    "h", AxisDirection::UP, verticalCS->axisList()[0]->unit()))
+                .as_nullable();
+    }
+
     auto &props = buildProperties(node);
 
     if (esriStyle_ && dbContext_) {
@@ -4246,10 +4276,10 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
                                "North American Vertical Datum 1988");
                 propsDatum.set(Identifier::CODE_KEY, 5103);
                 propsDatum.set(Identifier::CODESPACE_KEY, Identifier::EPSG);
-                datum =
+                vdatum =
                     VerticalReferenceFrame::create(propsDatum).as_nullable();
                 const auto dummyCRS =
-                    VerticalCRS::create(PropertyMap(), datum, datumEnsemble,
+                    VerticalCRS::create(PropertyMap(), vdatum, datumEnsemble,
                                         NN_NO_CHECK(verticalCS));
                 const auto model(Transformation::create(
                     propsModel, dummyCRS, dummyCRS, nullptr,
@@ -4265,7 +4295,7 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
     if (!isNull(geoidModelNode)) {
         auto &propsModel = buildProperties(geoidModelNode);
         const auto dummyCRS = VerticalCRS::create(
-            PropertyMap(), datum, datumEnsemble, NN_NO_CHECK(verticalCS));
+            PropertyMap(), vdatum, datumEnsemble, NN_NO_CHECK(verticalCS));
         const auto model(Transformation::create(
             propsModel, dummyCRS, dummyCRS, nullptr,
             OperationMethod::create(PropertyMap(),
@@ -4275,10 +4305,10 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
     }
 
     auto crs = nn_static_pointer_cast<CRS>(VerticalCRS::create(
-        props, datum, datumEnsemble, NN_NO_CHECK(verticalCS)));
+        props, vdatum, datumEnsemble, NN_NO_CHECK(verticalCS)));
 
-    if (!isNull(datumNode)) {
-        auto &extensionNode = datumNode->lookForChild(WKTConstants::EXTENSION);
+    if (!isNull(vdatumNode)) {
+        auto &extensionNode = vdatumNode->lookForChild(WKTConstants::EXTENSION);
         const auto &extensionChildren = extensionNode->GP()->children();
         if (extensionChildren.size() == 2) {
             if (ci_equal(stripQuotes(extensionChildren[0]), "PROJ4_GRIDS")) {
@@ -6507,6 +6537,34 @@ BaseObjectNNPtr WKTParser::createFromWKT(const std::string &wkt) {
             }
             return d->buildGeodeticReferenceFrame(root, primeMeridian,
                                                   null_node);
+        } else if (ci_equal(name, WKTConstants::GEOGCS) ||
+                   ci_equal(name, WKTConstants::PROJCS)) {
+            // Parse implicit compoundCRS from ESRI that is
+            // "PROJCS[...],VERTCS[...]" or "GEOGCS[...],VERTCS[...]"
+            if (indexEnd < wkt.size()) {
+                indexEnd = skipSpace(wkt, indexEnd);
+                if (indexEnd < wkt.size() && wkt[indexEnd] == ',') {
+                    ++indexEnd;
+                    indexEnd = skipSpace(wkt, indexEnd);
+                    if (indexEnd < wkt.size() &&
+                        ci_starts_with(wkt.c_str() + indexEnd,
+                                       WKTConstants::VERTCS.c_str())) {
+                        auto horizCRS = d->buildCRS(root);
+                        if (horizCRS) {
+                            auto vertCRS =
+                                d->buildVerticalCRS(WKTNode::createFrom(
+                                    wkt, indexEnd, 0, indexEnd));
+                            return CompoundCRS::createLax(
+                                util::PropertyMap().set(
+                                    IdentifiedObject::NAME_KEY,
+                                    horizCRS->nameStr() + " + " +
+                                        vertCRS->nameStr()),
+                                {NN_NO_CHECK(horizCRS), vertCRS},
+                                d->dbContext_);
+                        }
+                    }
+                }
+            }
         }
         return d->build(root);
     };