Skip to content

Commit

Permalink
Fix importing '+proj=topocentric ... +type=crs' by using a geocentric…
Browse files Browse the repository at this point in the history 
… CRS as the base CRS

Fixes https://lists.osgeo.org/pipermail/proj/2023-October/011153.html
  • Loading branch information
@rouault
rouault committed Nov 19, 2023
1 parent c0a118f commit 72de8f7
Show file tree
Hide file tree
Showing 4 changed files with 90 additions and 1 deletion.
31 changes: 30 additions & 1 deletion src/iso19111/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()) {
Expand Down Expand Up @@ -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"));

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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<BaseObject>(obj);
}
}
}

if (!unexpectedStructure) {
if (iFirstGeogStep == 0 && !d->steps_[iFirstGeogStep].inverted &&
iSecondGeogStep < 0 && iProjStep < 0 &&
Expand Down
7 changes: 7 additions & 0 deletions test/cli/testvarious
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1338,6 +1338,13 @@ $EXE -d 3 EPSG:3912+EPSG:5779 EPSG:3794+EPSG:8690 -E >>${OUT} 2>&1 <<EOF
477134.28 95134.21 5
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

echo "##############################################################" >> ${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}
Expand Down
3 changes: 3 additions & 0 deletions test/cli/tv_out.dist
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
50 changes: 50 additions & 0 deletions test/unit/test_io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<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);
Expand Down

0 comments on commit 72de8f7

Please sign in to comment.