Skip to content

Commit

Permalink
createOperations(): fix wrong pipeline generation with CRS that has +…
Browse files Browse the repository at this point in the history 
…nadgrids= and +pm=

Fixes issue reported at
https://lists.osgeo.org/pipermail/gdal-dev/2020-February/051749.html

The generated pipeline assumes that the input coordinates for the grid transformation
were related to the non-Greenwich based datum, so we must compensate for that and
add logic to go back to Greenwich.
  • Loading branch information
@rouault
rouault committed Feb 28, 2020
1 parent e81f938 commit 47cece3
Show file tree
Hide file tree
Showing 4 changed files with 168 additions and 42 deletions.
39 changes: 35 additions & 4 deletions src/iso19111/coordinateoperation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13850,12 +13850,27 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog(
auto opsFirst = createOperations(
boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
if (!opsFirst.empty()) {
CoordinateOperationPtr opIntermediate;
if (!geogCRSOfBaseOfBoundSrc->_isEquivalentTo(
boundSrc->transformation()->sourceCRS().get(),
util::IComparable::Criterion::EQUIVALENT)) {
auto opsIntermediate = createOperations(
NN_NO_CHECK(geogCRSOfBaseOfBoundSrc),
boundSrc->transformation()->sourceCRS(), context);
assert(!opsIntermediate.empty());
opIntermediate = opsIntermediate.front();
}
for (const auto &opFirst : opsFirst) {
try {
std::vector<CoordinateOperationNNPtr> subops;
subops.emplace_back(opFirst);
if (opIntermediate) {
subops.emplace_back(NN_NO_CHECK(opIntermediate));
}
subops.emplace_back(boundSrc->transformation());
res.emplace_back(
ConcatenatedOperation::createComputeMetadata(
{opFirst, boundSrc->transformation()},
!allowEmptyIntersection));
subops, !allowEmptyIntersection));
} catch (const InvalidOperationEmptyIntersection &) {
}
}
Expand All @@ -13873,13 +13888,29 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog(
boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
auto opsLast = createOperations(hubSrc, targetCRS, context);
if (!opsFirst.empty() && !opsLast.empty()) {
CoordinateOperationPtr opIntermediate;
if (!geogCRSOfBaseOfBoundSrc->_isEquivalentTo(
boundSrc->transformation()->sourceCRS().get(),
util::IComparable::Criterion::EQUIVALENT)) {
auto opsIntermediate = createOperations(
NN_NO_CHECK(geogCRSOfBaseOfBoundSrc),
boundSrc->transformation()->sourceCRS(), context);
assert(!opsIntermediate.empty());
opIntermediate = opsIntermediate.front();
}
for (const auto &opFirst : opsFirst) {
for (const auto &opLast : opsLast) {
try {
std::vector<CoordinateOperationNNPtr> subops;
subops.emplace_back(opFirst);
if (opIntermediate) {
subops.emplace_back(NN_NO_CHECK(opIntermediate));
}
subops.emplace_back(boundSrc->transformation());
subops.emplace_back(opLast);
res.emplace_back(
ConcatenatedOperation::createComputeMetadata(
{opFirst, boundSrc->transformation(), opLast},
!allowEmptyIntersection));
subops, !allowEmptyIntersection));
} catch (const InvalidOperationEmptyIntersection &) {
}
}
Expand Down
26 changes: 22 additions & 4 deletions src/iso19111/crs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4525,9 +4525,27 @@ BoundCRS::createFromTOWGS84(const CRSNNPtr &baseCRSIn,
*/
BoundCRSNNPtr BoundCRS::createFromNadgrids(const CRSNNPtr &baseCRSIn,
const std::string &filename) {
const CRSPtr sourceGeographicCRS = baseCRSIn->extractGeographicCRS();
const auto sourceGeographicCRS = baseCRSIn->extractGeographicCRS();
auto transformationSourceCRS =
sourceGeographicCRS ? sourceGeographicCRS : baseCRSIn.as_nullable();
sourceGeographicCRS
? NN_NO_CHECK(std::static_pointer_cast<CRS>(sourceGeographicCRS))
: baseCRSIn;
if (sourceGeographicCRS != nullptr &&
sourceGeographicCRS->datum() != nullptr &&
sourceGeographicCRS->primeMeridian()->longitude().getSIValue() != 0.0) {
transformationSourceCRS = GeographicCRS::create(
util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
sourceGeographicCRS->nameStr() +
" (with Greenwich prime meridian)"),
datum::GeodeticReferenceFrame::create(
util::PropertyMap().set(
common::IdentifiedObject::NAME_KEY,
sourceGeographicCRS->datum()->nameStr() +
" (with Greenwich prime meridian)"),
sourceGeographicCRS->datum()->ellipsoid(),
util::optional<std::string>(), datum::PrimeMeridian::GREENWICH),
sourceGeographicCRS->coordinateSystem());
}
std::string transformationName = transformationSourceCRS->nameStr();
transformationName += " to WGS84";

Expand All @@ -4536,8 +4554,8 @@ BoundCRSNNPtr BoundCRS::createFromNadgrids(const CRSNNPtr &baseCRSIn,
operation::Transformation::createNTv2(
util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
transformationName),
NN_NO_CHECK(transformationSourceCRS), GeographicCRS::EPSG_4326,
filename, std::vector<metadata::PositionalAccuracyNNPtr>()));
transformationSourceCRS, GeographicCRS::EPSG_4326, filename,
std::vector<metadata::PositionalAccuracyNNPtr>()));
}

// ---------------------------------------------------------------------------
Expand Down
86 changes: 53 additions & 33 deletions src/iso19111/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4060,6 +4060,52 @@ WKTParser::Private::buildCompoundCRS(const WKTNodeNNPtr &node) {

// ---------------------------------------------------------------------------

static CRSNNPtr
createBoundCRSSourceTransformationCRS(const crs::CRSPtr &sourceCRS,
const crs::CRSPtr &targetCRS) {
CRSPtr sourceTransformationCRS;
if (dynamic_cast<GeographicCRS *>(targetCRS.get())) {
GeographicCRSPtr sourceGeographicCRS =
sourceCRS->extractGeographicCRS();
sourceTransformationCRS = sourceGeographicCRS;
if (sourceGeographicCRS) {
if (sourceGeographicCRS->datum() != nullptr &&
sourceGeographicCRS->primeMeridian()
->longitude()
.getSIValue() != 0.0) {
sourceTransformationCRS =
GeographicCRS::create(
util::PropertyMap().set(
common::IdentifiedObject::NAME_KEY,
sourceGeographicCRS->nameStr() +
" (with Greenwich prime meridian)"),
datum::GeodeticReferenceFrame::create(
util::PropertyMap().set(
common::IdentifiedObject::NAME_KEY,
sourceGeographicCRS->datum()->nameStr() +
" (with Greenwich prime meridian)"),
sourceGeographicCRS->datum()->ellipsoid(),
util::optional<std::string>(),
datum::PrimeMeridian::GREENWICH),
sourceGeographicCRS->coordinateSystem())
.as_nullable();
}
} else {
sourceTransformationCRS =
std::dynamic_pointer_cast<VerticalCRS>(sourceCRS);
if (!sourceTransformationCRS) {
throw ParsingException(
"Cannot find GeographicCRS or VerticalCRS in sourceCRS");
}
}
} else {
sourceTransformationCRS = sourceCRS;
}
return NN_NO_CHECK(sourceTransformationCRS);
}

// ---------------------------------------------------------------------------

BoundCRSNNPtr WKTParser::Private::buildBoundCRS(const WKTNodeNNPtr &node) {
const auto *nodeP = node->GP();
auto &abridgedNode =
Expand Down Expand Up @@ -4103,23 +4149,10 @@ BoundCRSNNPtr WKTParser::Private::buildBoundCRS(const WKTNodeNNPtr &node) {
consumeParameters(abridgedNode, true, parameters, values, defaultLinearUnit,
defaultAngularUnit);

CRSPtr sourceTransformationCRS;
if (dynamic_cast<GeographicCRS *>(targetCRS.get())) {
sourceTransformationCRS = sourceCRS->extractGeographicCRS();
if (!sourceTransformationCRS) {
sourceTransformationCRS =
std::dynamic_pointer_cast<VerticalCRS>(sourceCRS);
if (!sourceTransformationCRS) {
throw ParsingException(
"Cannot find GeographicCRS or VerticalCRS in sourceCRS");
}
}
} else {
sourceTransformationCRS = sourceCRS;
}

const auto sourceTransformationCRS(
createBoundCRSSourceTransformationCRS(sourceCRS, targetCRS));
auto transformation = Transformation::create(
buildProperties(abridgedNode), NN_NO_CHECK(sourceTransformationCRS),
buildProperties(abridgedNode), sourceTransformationCRS,
NN_NO_CHECK(targetCRS), nullptr, buildProperties(methodNode),
parameters, values, std::vector<PositionalAccuracyNNPtr>());

Expand Down Expand Up @@ -5265,24 +5298,11 @@ BoundCRSNNPtr JSONParser::buildBoundCRS(const json &j) {
values.emplace_back(ParameterValue::create(getMeasure(param)));
}

CRSPtr sourceTransformationCRS;
if (dynamic_cast<GeographicCRS *>(targetCRS.get())) {
sourceTransformationCRS = sourceCRS->extractGeographicCRS();
if (!sourceTransformationCRS) {
sourceTransformationCRS =
std::dynamic_pointer_cast<VerticalCRS>(sourceCRS.as_nullable());
if (!sourceTransformationCRS) {
throw ParsingException(
"Cannot find GeographicCRS or VerticalCRS in sourceCRS");
}
}
} else {
sourceTransformationCRS = sourceCRS;
}

const auto sourceTransformationCRS(
createBoundCRSSourceTransformationCRS(sourceCRS, targetCRS));
auto transformation = Transformation::create(
buildProperties(transformationJ), NN_NO_CHECK(sourceTransformationCRS),
targetCRS, nullptr, buildProperties(methodJ), parameters, values,
buildProperties(transformationJ), sourceTransformationCRS, targetCRS,
nullptr, buildProperties(methodJ), parameters, values,
std::vector<PositionalAccuracyNNPtr>());

return BoundCRS::create(sourceCRS, targetCRS, transformation);
Expand Down
59 changes: 58 additions & 1 deletion test/unit/test_operation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6567,6 +6567,64 @@ TEST(operation, ETRS89_3D_to_proj_string_with_geoidgrids_nadgrids) {

// ---------------------------------------------------------------------------

TEST(operation, nadgrids_with_pm) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
auto objSrc = PROJStringParser().createFromPROJString(
"+proj=tmerc +lat_0=39.66666666666666 +lon_0=1 +k=1 +x_0=200000 "
"+y_0=300000 +ellps=intl +nadgrids=foo.gsb +pm=lisbon "
"+units=m +type=crs");
auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
ASSERT_TRUE(src != nullptr);
auto dst = authFactory->createCoordinateReferenceSystem("4326");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
auto list = CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(src), dst, ctxt);
ASSERT_EQ(list.size(), 1U);
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +inv +proj=tmerc +lat_0=39.6666666666667 +lon_0=1 "
"+k=1 +x_0=200000 +y_0=300000 +ellps=intl +pm=lisbon "
// Check that there is no extra +step +proj=longlat +pm=lisbon
"+step +proj=hgridshift +grids=foo.gsb "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");

// ETRS89
dst = authFactory->createCoordinateReferenceSystem("4258");
list = CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(src), dst, ctxt);
ASSERT_EQ(list.size(), 1U);
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +inv +proj=tmerc +lat_0=39.6666666666667 +lon_0=1 "
"+k=1 +x_0=200000 +y_0=300000 +ellps=intl +pm=lisbon "
// Check that there is no extra +step +proj=longlat +pm=lisbon
"+step +proj=hgridshift +grids=foo.gsb "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");

// From WKT BOUNDCRS
auto formatter = WKTFormatter::create(WKTFormatter::Convention::WKT2_2019);
auto src_wkt = src->exportToWKT(formatter.get());
auto objFromWkt = WKTParser().createFromWKT(src_wkt);
auto crsFromWkt = nn_dynamic_pointer_cast<BoundCRS>(objFromWkt);
ASSERT_TRUE(crsFromWkt);
list = CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(crsFromWkt), dst, ctxt);
ASSERT_EQ(list.size(), 1U);
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +inv +proj=tmerc +lat_0=39.6666666666667 +lon_0=1 "
"+k=1 +x_0=200000 +y_0=300000 +ellps=intl +pm=lisbon "
// Check that there is no extra +step +proj=longlat +pm=lisbon
"+step +proj=hgridshift +grids=foo.gsb "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}

// ---------------------------------------------------------------------------

TEST(operation, WGS84_G1762_to_compoundCRS_with_bound_vertCRS) {
auto authFactoryEPSG =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
Expand Down Expand Up @@ -6767,7 +6825,6 @@ TEST(operation, transformation_BEV_AT_to_PROJ_string) {
"+multiplier=1");
}


// ---------------------------------------------------------------------------

TEST(operation, transformation_longitude_rotation_to_PROJ_string) {
Expand Down

0 comments on commit 47cece3

Please sign in to comment.