From 451a2a48d3ec461b56094cce623f9f91c5c457bf Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 18 Oct 2023 18:40:18 +0200 Subject: [PATCH] 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 --- src/iso19111/io.cpp | 31 ++++++++++++++++++++++++++- test/cli/testvarious | 7 ++++++ test/cli/tv_out.dist | 3 +++ test/unit/test_io.cpp | 50 +++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 90 insertions(+), 1 deletion(-) diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index be2d10f036..31c770487e 100644 --- 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(obj); + } + } + } + if (!unexpectedStructure) { if (iFirstGeogStep == 0 && !d->steps_[iFirstGeogStep].inverted && iSecondGeogStep < 0 && iProjStep < 0 && diff --git a/test/cli/testvarious b/test/cli/testvarious index 8134e54dca..a9589f500e 100755 --- a/test/cli/testvarious +++ b/test/cli/testvarious @@ -1352,6 +1352,13 @@ $EXE -d 8 --t_epoch 2022 "GDA2020" "ITRF2014" -E >>${OUT} 2>&1 <> ${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 <(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);