From 72de8f757c4f4585c328c3daf3f33ea04af08a36 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 f5266d9920..9b0298f46b 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -10415,6 +10415,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()) { @@ -11091,7 +11097,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")); @@ -11743,6 +11749,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); @@ -11849,6 +11862,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) { @@ -12241,6 +12258,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 38c5a75686..70b027e3b6 100755 --- a/test/cli/testvarious +++ b/test/cli/testvarious @@ -1338,6 +1338,13 @@ $EXE -d 3 EPSG:3912+EPSG:5779 EPSG:3794+EPSG:8690 -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 <> ${OUT} echo "Test cs2cs EPSG:5488 (RGAF09) to EPSG:4559+5757 (RRAF 1991 / UTM zone 20N + Guadeloupe 1988 height)" >> ${OUT} echo "Check that we use the horizontal transformation for Guadeloupe (RRAF 1991 to RGAF09 (2)), and not the one for Martinique" >> ${OUT} diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist index b2418c52f5..56ba877ce0 100644 --- a/test/cli/tv_out.dist +++ b/test/cli/tv_out.dist @@ -641,6 +641,9 @@ Test Similarity Transformation through CompoundCRS 477134.28 95134.21 476763.303 95620.222 0.000 477134.28 95134.21 5 476763.303 95620.222 5.000 ############################################################## +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 +############################################################## Test cs2cs EPSG:5488 (RGAF09) to EPSG:4559+5757 (RRAF 1991 / UTM zone 20N + Guadeloupe 1988 height) Check that we use the horizontal transformation for Guadeloupe (RRAF 1991 to RGAF09 (2)), and not the one for Martinique 661991.318 1796999.201 93.846 diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 7c62dda885..5d8e817429 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -12085,6 +12085,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(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);