Import from WKT: add tweaks for Lidar WKT1 VERT_CS that embeds geoid model in CRS name

2 files changed, 123 insertions, 0 deletions

diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 116a8b78..d4e6132f 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -3821,6 +3821,64 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
     }
 
     auto &props = buildProperties(node);
+
+    // Deal with Lidar WKT1 VertCRS that embeds geoid model in CRS name,
+    // following conventions from
+    // https://pubs.usgs.gov/tm/11b4/pdf/tm11-B4.pdf
+    // page 9
+    if (ci_equal(nodeValue, WKTConstants::VERT_CS)) {
+        std::string name;
+        if (props.getStringValue(IdentifiedObject::NAME_KEY, name)) {
+            std::string geoidName;
+            for (const char *prefix :
+                 {"NAVD88 - ", "NAVD88 via ", "NAVD88 height - ",
+                  "NAVD88 height (ftUS) - "}) {
+                if (starts_with(name, prefix)) {
+                    geoidName = name.substr(strlen(prefix));
+                    auto pos = geoidName.find_first_of(" (");
+                    if (pos != std::string::npos) {
+                        geoidName.resize(pos);
+                    }
+                    break;
+                }
+            }
+            if (!geoidName.empty()) {
+                const auto &axis = verticalCS->axisList()[0];
+                const auto &dir = axis->direction();
+                if (dir == cs::AxisDirection::UP) {
+                    if (axis->unit() == common::UnitOfMeasure::METRE) {
+                        props.set(IdentifiedObject::NAME_KEY, "NAVD88 height");
+                        props.set(Identifier::CODE_KEY, 5703);
+                        props.set(Identifier::CODESPACE_KEY, Identifier::EPSG);
+                    } else if (axis->unit().name() == "US survey foot") {
+                        props.set(IdentifiedObject::NAME_KEY,
+                                  "NAVD88 height (ftUS)");
+                        props.set(Identifier::CODE_KEY, 6360);
+                        props.set(Identifier::CODESPACE_KEY, Identifier::EPSG);
+                    }
+                }
+                PropertyMap propsModel;
+                propsModel.set(IdentifiedObject::NAME_KEY, toupper(geoidName));
+                PropertyMap propsDatum;
+                propsDatum.set(IdentifiedObject::NAME_KEY,
+                               "North American Vertical Datum 1988");
+                propsDatum.set(Identifier::CODE_KEY, 5103);
+                propsDatum.set(Identifier::CODESPACE_KEY, Identifier::EPSG);
+                datum =
+                    VerticalReferenceFrame::create(propsDatum).as_nullable();
+                const auto dummyCRS =
+                    VerticalCRS::create(PropertyMap(), datum, datumEnsemble,
+                                        NN_NO_CHECK(verticalCS));
+                const auto model(Transformation::create(
+                    propsModel, dummyCRS, dummyCRS, nullptr,
+                    OperationMethod::create(
+                        PropertyMap(), std::vector<OperationParameterNNPtr>()),
+                    {}, {}));
+                props.set("GEOID_MODEL", model);
+            }
+        }
+    }
+
     auto &geoidModelNode = nodeP->lookForChild(WKTConstants::GEOIDMODEL);
     if (!isNull(geoidModelNode)) {
         auto &propsModel = buildProperties(geoidModelNode);
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 6bb93e97..79851695 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -2031,6 +2031,71 @@ TEST(wkt_parse, vertcrs_WKT1_GDAL_minimum) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, vertcrs_WKT1_LAS_ftUS) {
+    auto wkt = "VERT_CS[\"NAVD88 - Geoid03 (Feet)\","
+               "    VERT_DATUM[\"unknown\",2005],"
+               "    UNIT[\"US survey foot\",0.3048006096012192,"
+               "        AUTHORITY[\"EPSG\",\"9003\"]],"
+               "    AXIS[\"Up\",UP]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<VerticalCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->nameStr(), "NAVD88 height (ftUS)");
+    ASSERT_EQ(crs->identifiers().size(), 1U);
+    EXPECT_EQ(crs->identifiers()[0]->code(), "6360");
+    EXPECT_EQ(*(crs->identifiers()[0]->codeSpace()), "EPSG");
+
+    const auto &geoidModel = crs->geoidModel();
+    ASSERT_TRUE(!geoidModel.empty());
+    EXPECT_EQ(geoidModel[0]->nameStr(), "GEOID03");
+
+    auto datum = crs->datum();
+    EXPECT_EQ(datum->nameStr(), "North American Vertical Datum 1988");
+    ASSERT_EQ(datum->identifiers().size(), 1U);
+    EXPECT_EQ(datum->identifiers()[0]->code(), "5103");
+    EXPECT_EQ(*(datum->identifiers()[0]->codeSpace()), "EPSG");
+
+    const auto &axis = crs->coordinateSystem()->axisList()[0];
+    EXPECT_EQ(axis->direction(), AxisDirection::UP);
+    EXPECT_EQ(axis->unit().name(), "US survey foot");
+    EXPECT_NEAR(axis->unit().conversionToSI(), 0.3048006096012192, 1e-16);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, vertcrs_WKT1_LAS_metre) {
+    auto wkt = "VERT_CS[\"NAVD88 via Geoid09\","
+               "    VERT_DATUM[\"unknown\",2005],"
+               "    UNIT[\"metre\",1.0,"
+               "        AUTHORITY[\"EPSG\",\"9001\"]],"
+               "    AXIS[\"Up\",UP]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<VerticalCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->nameStr(), "NAVD88 height");
+    ASSERT_EQ(crs->identifiers().size(), 1U);
+    EXPECT_EQ(crs->identifiers()[0]->code(), "5703");
+    EXPECT_EQ(*(crs->identifiers()[0]->codeSpace()), "EPSG");
+
+    const auto &geoidModel = crs->geoidModel();
+    ASSERT_TRUE(!geoidModel.empty());
+    EXPECT_EQ(geoidModel[0]->nameStr(), "GEOID09");
+
+    auto datum = crs->datum();
+    EXPECT_EQ(datum->nameStr(), "North American Vertical Datum 1988");
+    ASSERT_EQ(datum->identifiers().size(), 1U);
+    EXPECT_EQ(datum->identifiers()[0]->code(), "5103");
+    EXPECT_EQ(*(datum->identifiers()[0]->codeSpace()), "EPSG");
+
+    const auto &axis = crs->coordinateSystem()->axisList()[0];
+    EXPECT_EQ(axis->direction(), AxisDirection::UP);
+    EXPECT_EQ(axis->unit(), UnitOfMeasure::METRE);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, dynamic_vertical_reference_frame) {
     auto obj = WKTParser().createFromWKT(
         "VERTCRS[\"RH2000\","