Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

createOperations(): fix wrong pipeline generation with CRS that has +nadgrids= and +pm= #1998

Merged
merged 1 commit into from
Feb 29, 2020
Merged

createOperations(): fix wrong pipeline generation with CRS that has +nadgrids= and +pm= #1998

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 35 additions & 4 deletions src/iso19111/coordinateoperation.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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