Merge pull request #4124 from rouault/fix_gdal_9732

DerivedProjectedCRS: deal with lack of explicit CS in BASEPROJCRS in WKT
OSGeo · May 7, 2024 · a715623 · a715623
2 parents b59c4b5 + c873334
commit a715623
Show file tree

Hide file tree

Showing 7 changed files with 1,402 additions and 1,049 deletions.
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
@@ -6706,6 +6706,37 @@ void DerivedProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
         formatter->endNode();
 
         l_baseProjCRS->derivingConversionRef()->_exportToWKT(formatter);
+
+        const auto &baseCSAxisList =
+            l_baseProjCRS->coordinateSystem()->axisList();
+        // Current WKT grammar (as of WKT2 18-010r11) does not allow a
+        // BASEPROJCRS.CS node, but in situations where this is ambiguous, emit
+        // one. Cf WKTParser::Private::buildProjectedCRS() for more details
+        if (!baseCSAxisList.empty() &&
+            baseCSAxisList[0]->unit() != common::UnitOfMeasure::METRE &&
+            l_baseProjCRS->identifiers().empty()) {
+            bool knownBaseCRS = false;
+            auto &dbContext = formatter->databaseContext();
+            if (dbContext) {
+                auto authFactory = io::AuthorityFactory::create(
+                    NN_NO_CHECK(dbContext), std::string());
+                auto res = authFactory->createObjectsFromName(
+                    l_baseProjCRS->nameStr(),
+                    {io::AuthorityFactory::ObjectType::PROJECTED_CRS}, false,
+                    2);
+                if (res.size() == 1) {
+                    knownBaseCRS = true;
+                }
+            }
+            if (!knownBaseCRS) {
+                l_baseProjCRS->coordinateSystem()->_exportToWKT(formatter);
+            }
+        }
+
+        if (identifiers().empty() && !l_baseProjCRS->identifiers().empty()) {
+            l_baseProjCRS->formatID(formatter);
+        }
+
         formatter->endNode();
     }
 

diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
@@ -4490,10 +4490,80 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
         !ci_equal(nodeValue, WKTConstants::BASEPROJCRS)) {
         ThrowMissing(WKTConstants::CS_);
     }
-    auto cs = buildCS(csNode, node, UnitOfMeasure::NONE);
-    auto cartesianCS = nn_dynamic_pointer_cast<CartesianCS>(cs);
 
     const std::string projCRSName = stripQuotes(nodeP->children()[0]);
+
+    auto cs = [this, &projCRSName, &nodeP, &csNode, &node, &nodeValue,
+               &conversionNode]() -> CoordinateSystemNNPtr {
+        if (isNull(csNode) && ci_equal(nodeValue, WKTConstants::BASEPROJCRS) &&
+            !isNull(conversionNode)) {
+            // A BASEPROJCRS (as of WKT2 18-010r11) normally lacks an explicit
+            // CS[] which cause issues to properly instanciate it. So we first
+            // start by trying to identify the BASEPROJCRS by its id or name.
+            // And fallback to exploring the conversion parameters to infer the
+            // CS AXIS unit from the linear parameter unit... Not fully bullet
+            // proof.
+            if (dbContext_) {
+                // Get official name from database if ID is present
+                auto &idNode = nodeP->lookForChild(WKTConstants::ID);
+                if (!isNull(idNode)) {
+                    try {
+                        auto id = buildId(idNode, false, false);
+                        auto authFactory = AuthorityFactory::create(
+                            NN_NO_CHECK(dbContext_), *id->codeSpace());
+                        auto projCRS =
+                            authFactory->createProjectedCRS(id->code());
+                        return projCRS->coordinateSystem();
+                    } catch (const std::exception &) {
+                    }
+                }
+
+                auto authFactory = AuthorityFactory::create(
+                    NN_NO_CHECK(dbContext_), std::string());
+                auto res = authFactory->createObjectsFromName(
+                    projCRSName, {AuthorityFactory::ObjectType::PROJECTED_CRS},
+                    false, 2);
+                if (res.size() == 1) {
+                    auto projCRS =
+                        dynamic_cast<const ProjectedCRS *>(res.front().get());
+                    if (projCRS) {
+                        return projCRS->coordinateSystem();
+                    }
+                }
+            }
+
+            auto conv = buildConversion(conversionNode, UnitOfMeasure::METRE,
+                                        UnitOfMeasure::DEGREE);
+            UnitOfMeasure linearUOM = UnitOfMeasure::NONE;
+            for (const auto &genOpParamvalue : conv->parameterValues()) {
+                auto opParamvalue =
+                    dynamic_cast<const operation::OperationParameterValue *>(
+                        genOpParamvalue.get());
+                if (opParamvalue) {
+                    const auto &parameterValue = opParamvalue->parameterValue();
+                    if (parameterValue->type() ==
+                        operation::ParameterValue::Type::MEASURE) {
+                        const auto &measure = parameterValue->value();
+                        const auto &unit = measure.unit();
+                        if (unit.type() == UnitOfMeasure::Type::LINEAR) {
+                            if (linearUOM == UnitOfMeasure::NONE) {
+                                linearUOM = unit;
+                            } else if (linearUOM != unit) {
+                                linearUOM = UnitOfMeasure::NONE;
+                                break;
+                            }
+                        }
+                    }
+                }
+            }
+            if (linearUOM != UnitOfMeasure::NONE) {
+                return CartesianCS::createEastingNorthing(linearUOM);
+            }
+        }
+        return buildCS(csNode, node, UnitOfMeasure::NONE);
+    }();
+    auto cartesianCS = nn_dynamic_pointer_cast<CartesianCS>(cs);
+
     if (esriStyle_ && dbContext_) {
         if (cartesianCS) {
             std::string outTableName;

diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
@@ -344,7 +344,7 @@ add_custom_target(check_wkt2_grammar_md5 ALL
                   COMMAND ${CMAKE_COMMAND}
                       "-DIN_FILE=wkt2_grammar.y"
                       "-DTARGET=generate_wkt2_parser"
-                      "-DEXPECTED_MD5SUM=289572eebe9dab3c7225bd48c445c287"
+                      "-DEXPECTED_MD5SUM=94337ca471caf74688db939009652ddb"
                       -P "${CMAKE_CURRENT_SOURCE_DIR}/check_md5sum.cmake"
                   WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
                   DEPENDS "${CMAKE_CURRENT_SOURCE_DIR}/wkt2_grammar.y"