Fix importing '+proj=topocentric ... +type=crs' by using a geocentric…

… CRS as the base CRS Fixes https://lists.osgeo.org/pipermail/proj/2023-October/011153.html
OSGeo · Oct 18, 2023 · 451a2a4 · 451a2a4
1 parent 3920d63
commit 451a2a4
Show file tree

Hide file tree

Showing 4 changed files with 90 additions and 1 deletion.
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
@@ -10297,6 +10297,12 @@ static bool isGeocentricStep(const std::string &name) {
 
 // ---------------------------------------------------------------------------
 
+static bool isTopocentricStep(const std::string &name) {
+    return name == "topocentric";
+}
+
+// ---------------------------------------------------------------------------
+
 static bool isProjectedStep(const std::string &name) {
     if (name == "etmerc" || name == "utm" ||
         !getMappingsFromPROJName(name).empty()) {
@@ -10973,7 +10979,7 @@ GeodeticCRSNNPtr
 PROJStringParser::Private::buildGeocentricCRS(int iStep, int iUnitConvert) {
     auto &step = steps_[iStep];
 
-    assert(isGeocentricStep(step.name));
+    assert(isGeocentricStep(step.name) || isTopocentricStep(step.name));
     assert(iUnitConvert < 0 ||
            ci_equal(steps_[iUnitConvert].name, "unitconvert"));
 
@@ -11625,6 +11631,13 @@ PROJStringParser::Private::buildProjectedCRS(int iStep,
                   ? CartesianCS::create(emptyPropertyMap, axis[0], axis[1])
                   : CartesianCS::create(emptyPropertyMap, axis[0], axis[1],
                                         csGeodCRS->axisList()[2]);
+    if (isTopocentricStep(step.name)) {
+        cs = CartesianCS::create(
+            emptyPropertyMap,
+            createAxis("topocentric East", "U", AxisDirection::EAST, unit),
+            createAxis("topocentric North", "V", AxisDirection::NORTH, unit),
+            createAxis("topocentric Up", "W", AxisDirection::UP, unit));
+    }
 
     auto props = PropertyMap().set(IdentifiedObject::NAME_KEY,
                                    title.empty() ? "unknown" : title);
@@ -11731,6 +11744,10 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
          (d->steps_.size() == 2 && d->steps_[1].name == "unitconvert")) &&
         !d->steps_[0].inverted && isGeocentricStep(d->steps_[0].name);
 
+    const bool isTopocentricCRS =
+        (d->steps_.size() == 1 && isTopocentricStep(d->steps_[0].name) &&
+         d->getParamValue(d->steps_[0], "type") == "crs");
+
     // +init=xxxx:yyyy syntax
     if (d->steps_.size() == 1 && d->steps_[0].isInit &&
         !d->steps_[0].inverted) {
@@ -12123,6 +12140,18 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
         }
     }
 
+    if (isTopocentricCRS) {
+        // First run is dry run to mark all recognized/unrecognized tokens
+        for (int iter = 0; iter < 2; iter++) {
+            auto obj = d->buildBoundOrCompoundCRSIfNeeded(
+                0,
+                d->buildProjectedCRS(0, d->buildGeocentricCRS(0, -1), -1, -1));
+            if (iter == 1) {
+                return nn_static_pointer_cast<BaseObject>(obj);
+            }
+        }
+    }
+
     if (!unexpectedStructure) {
         if (iFirstGeogStep == 0 && !d->steps_[iFirstGeogStep].inverted &&
             iSecondGeogStep < 0 && iProjStep < 0 &&

diff --git a/test/cli/testvarious b/test/cli/testvarious
@@ -1352,6 +1352,13 @@ $EXE -d 8 --t_epoch 2022 "GDA2020" "ITRF2014" -E >>${OUT} 2>&1 <<EOF
 -30 150
 EOF
 
+echo "##############################################################" >> ${OUT}
+echo "Test +proj=topocentric +datum=WGS84 +X_0=-3982059.42 +Y_0=3331314.88 +Z_0=3692463.58 +no_defs +type=crs +to EPSG:4978" >> ${OUT}
+#
+$EXE -d 3 +proj=topocentric +datum=WGS84 +X_0=-3982059.42 +Y_0=3331314.88 +Z_0=3692463.58 +no_defs +type=crs +to EPSG:4978 -E >>${OUT} 2>&1 <<EOF
+0 0 0
+EOF
+
 
 # Done!
 # do 'diff' with distribution results

diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
@@ -646,3 +646,6 @@ Test --s_epoch
 ##############################################################
 Test --t_epoch
 -30 150	-29.99999901	150.00000044 -0.00031879
+##############################################################
+Test +proj=topocentric +datum=WGS84 +X_0=-3982059.42 +Y_0=3331314.88 +Z_0=3692463.58 +no_defs +type=crs +to EPSG:4978
+0 0 0	-3982059.420	3331314.880 3692463.580
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
@@ -11860,6 +11860,56 @@ TEST(io, projparse_geoc) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_topocentric) {
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=topocentric +datum=WGS84 +X_0=-3982059.42 +Y_0=3331314.88 "
+        "+Z_0=3692463.58 +no_defs +type=crs");
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    auto expected =
+        "PROJCRS[\"unknown\",\n"
+        "    BASEGEODCRS[\"unknown\",\n"
+        "        DATUM[\"World Geodetic System 1984\",\n"
+        "            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
+        "                LENGTHUNIT[\"metre\",1]],\n"
+        "            ID[\"EPSG\",6326]],\n"
+        "        PRIMEM[\"Greenwich\",0,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "            ID[\"EPSG\",8901]]],\n"
+        "    CONVERSION[\"unknown\",\n"
+        "        METHOD[\"Geocentric/topocentric conversions\",\n"
+        "            ID[\"EPSG\",9836]],\n"
+        "        PARAMETER[\"Geocentric X of topocentric "
+        "origin\",-3982059.42,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8837]],\n"
+        "        PARAMETER[\"Geocentric Y of topocentric origin\",3331314.88,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8838]],\n"
+        "        PARAMETER[\"Geocentric Z of topocentric origin\",3692463.58,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8839]]],\n"
+        "    CS[Cartesian,3],\n"
+        "        AXIS[\"topocentric East (U)\",east,\n"
+        "            ORDER[1],\n"
+        "            LENGTHUNIT[\"metre\",1,\n"
+        "                ID[\"EPSG\",9001]]],\n"
+        "        AXIS[\"topocentric North (V)\",north,\n"
+        "            ORDER[2],\n"
+        "            LENGTHUNIT[\"metre\",1,\n"
+        "                ID[\"EPSG\",9001]]],\n"
+        "        AXIS[\"topocentric Up (W)\",up,\n"
+        "            ORDER[3],\n"
+        "            LENGTHUNIT[\"metre\",1,\n"
+        "                ID[\"EPSG\",9001]]]]";
+    EXPECT_EQ(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT2_2019).get()),
+        expected);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_projected_wktext) {
     std::string input("+proj=merc +foo +wktext +type=crs");
     auto obj = PROJStringParser().createFromPROJString(input);