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

Jump to bottom

proj_crs_create_bound_crs_to_WGS84(): make it work on verticalCRS/compoundCRS such as EPSG:4326+5773 and EPSG:4326+3855 #2365

Merged
merged 1 commit into from
Oct 5, 2020
Merged

proj_crs_create_bound_crs_to_WGS84(): make it work on verticalCRS/compoundCRS such as EPSG:4326+5773 and EPSG:4326+3855 #2365

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
30 changes: 21 additions & 9 deletions src/iso19111/coordinateoperation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8701,15 +8701,6 @@ _getHeightToGeographic3DFilename(const Transformation *op, bool allowInverse) {

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

//! @cond Doxygen_Suppress
const std::string &Transformation::getHeightToGeographic3DFilename() const {

return _getHeightToGeographic3DFilename(this, false);
}
//! @endcond

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

//! @cond Doxygen_Suppress
static bool
isGeographic3DToGravityRelatedHeight(const OperationMethodNNPtr &method,
Expand Down Expand Up @@ -8764,6 +8755,27 @@ isGeographic3DToGravityRelatedHeight(const OperationMethodNNPtr &method,

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

//! @cond Doxygen_Suppress
const std::string &Transformation::getHeightToGeographic3DFilename() const {

const std::string &ret = _getHeightToGeographic3DFilename(this, false);
if (!ret.empty())
return ret;
if (isGeographic3DToGravityRelatedHeight(method(), false)) {
const auto &fileParameter =
parameterValue(EPSG_NAME_PARAMETER_GEOID_CORRECTION_FILENAME,
EPSG_CODE_PARAMETER_GEOID_CORRECTION_FILENAME);
if (fileParameter &&
fileParameter->type() == ParameterValue::Type::FILENAME) {
return fileParameter->valueFile();
}
}
return nullString;
}
//! @endcond

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

//! @cond Doxygen_Suppress
static util::PropertyMap
createSimilarPropertiesMethod(common::IdentifiedObjectNNPtr obj) {
Expand Down
104 changes: 89 additions & 15 deletions src/iso19111/crs.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -401,23 +401,20 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
}
}

auto geodCRS = util::nn_dynamic_pointer_cast<GeodeticCRS>(thisAsCRS);
auto geogCRS = extractGeographicCRS();
auto hubCRS = util::nn_static_pointer_cast<CRS>(GeographicCRS::EPSG_4326);
if (geodCRS && !geogCRS) {
if (geodCRS->_isEquivalentTo(GeographicCRS::EPSG_4978.get(),
util::IComparable::Criterion::EQUIVALENT,
dbContext)) {
return thisAsCRS;
auto compoundCRS = dynamic_cast<const CompoundCRS *>(this);
if (compoundCRS) {
const auto &comps = compoundCRS->componentReferenceSystems();
if (comps.size() == 2) {
auto horiz = comps[0]->createBoundCRSToWGS84IfPossible(
dbContext, allowIntermediateCRSUse);
auto vert = comps[1]->createBoundCRSToWGS84IfPossible(
dbContext, allowIntermediateCRSUse);
if (horiz.get() != comps[0].get() || vert.get() != comps[1].get()) {
return CompoundCRS::create(createPropertyMap(this),
{horiz, vert});
}
}
hubCRS = util::nn_static_pointer_cast<CRS>(GeodeticCRS::EPSG_4978);
} else if (!geogCRS ||
geogCRS->_isEquivalentTo(
GeographicCRS::EPSG_4326.get(),
util::IComparable::Criterion::EQUIVALENT, dbContext)) {
return thisAsCRS;
} else {
geodCRS = geogCRS;
}

if (!dbContext) {
Expand All @@ -443,6 +440,83 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
if (authorities.empty()) {
authorities.emplace_back();
}

// Vertical CRS ?
auto vertCRS = dynamic_cast<const VerticalCRS *>(this);
if (vertCRS) {
auto hubCRS =
util::nn_static_pointer_cast<CRS>(GeographicCRS::EPSG_4979);
for (const auto &authority : authorities) {
try {

auto authFactory = io::AuthorityFactory::create(
NN_NO_CHECK(dbContext),
authority == "any" ? std::string() : authority);
auto ctxt = operation::CoordinateOperationContext::create(
authFactory, extent, 0.0);
ctxt->setAllowUseIntermediateCRS(allowIntermediateCRSUse);
// ctxt->setSpatialCriterion(
// operation::CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
auto list = operation::CoordinateOperationFactory::create()
->createOperations(hubCRS, thisAsCRS, ctxt);
CRSPtr candidateBoundCRS;
for (const auto &op : list) {
auto transf = util::nn_dynamic_pointer_cast<
operation::Transformation>(op);
// Only keep transformations that use a known grid
if (transf && !transf->hasBallparkTransformation()) {
auto gridsNeeded = transf->gridsNeeded(dbContext, true);
bool gridsKnown = !gridsNeeded.empty();
for (const auto &gridDesc : gridsNeeded) {
if (gridDesc.packageName.empty() &&
!(!gridDesc.url.empty() &&
gridDesc.openLicense) &&
!gridDesc.available) {
gridsKnown = false;
break;
}
}
if (gridsKnown) {
if (candidateBoundCRS) {
candidateBoundCRS = nullptr;
break;
}
candidateBoundCRS =
BoundCRS::create(thisAsCRS, hubCRS,
NN_NO_CHECK(transf))
.as_nullable();
}
}
}
if (candidateBoundCRS) {
return NN_NO_CHECK(candidateBoundCRS);
}
} catch (const std::exception &) {
}
}
return thisAsCRS;
}

// Geodetic/geographic CRS ?
auto geodCRS = util::nn_dynamic_pointer_cast<GeodeticCRS>(thisAsCRS);
auto geogCRS = extractGeographicCRS();
auto hubCRS = util::nn_static_pointer_cast<CRS>(GeographicCRS::EPSG_4326);
if (geodCRS && !geogCRS) {
if (geodCRS->_isEquivalentTo(GeographicCRS::EPSG_4978.get(),
util::IComparable::Criterion::EQUIVALENT,
dbContext)) {
return thisAsCRS;
}
hubCRS = util::nn_static_pointer_cast<CRS>(GeodeticCRS::EPSG_4978);
} else if (!geogCRS ||
geogCRS->_isEquivalentTo(
GeographicCRS::EPSG_4326.get(),
util::IComparable::Criterion::EQUIVALENT, dbContext)) {
return thisAsCRS;
} else {
geodCRS = geogCRS;
}

for (const auto &authority : authorities) {
try {

Expand Down
50 changes: 38 additions & 12 deletions test/unit/test_crs.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5452,24 +5452,43 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
"+proj=longlat +ellps=clrk80ign +pm=paris "
"+towgs84=-168,-60,320,0,0,0,0 +no_defs +type=crs");
}
{
// WGS 84 + EGM2008 height
auto obj = createFromUserInput("EPSG:4326+3855", dbContext);
auto crs = nn_dynamic_pointer_cast<CRS>(obj);
ASSERT_TRUE(crs != nullptr);
auto res = crs->createBoundCRSToWGS84IfPossible(
dbContext, CoordinateOperationContext::IntermediateCRSUse::NEVER);
EXPECT_NE(res, crs);
EXPECT_EQ(res->createBoundCRSToWGS84IfPossible(
dbContext,
CoordinateOperationContext::IntermediateCRSUse::NEVER),
res);
auto compoundCRS = nn_dynamic_pointer_cast<CompoundCRS>(res);
ASSERT_TRUE(compoundCRS != nullptr);
EXPECT_EQ(compoundCRS->exportToPROJString(
PROJStringFormatter::create().get()),
"+proj=longlat +datum=WGS84 +geoidgrids=us_nga_egm08_25.tif "
"+vunits=m +no_defs +type=crs");
}
{
// NTF (Paris) / Lambert zone II + NGF-IGN69 height
auto crs_7421 = factory->createCoordinateReferenceSystem("7421");
auto bound = crs_7421->createBoundCRSToWGS84IfPossible(
auto res = crs_7421->createBoundCRSToWGS84IfPossible(
dbContext, CoordinateOperationContext::IntermediateCRSUse::NEVER);
EXPECT_NE(bound, crs_7421);
EXPECT_EQ(bound->createBoundCRSToWGS84IfPossible(
EXPECT_NE(res, crs_7421);
EXPECT_EQ(res->createBoundCRSToWGS84IfPossible(
dbContext,
CoordinateOperationContext::IntermediateCRSUse::NEVER),
bound);
auto boundCRS = nn_dynamic_pointer_cast<BoundCRS>(bound);
ASSERT_TRUE(boundCRS != nullptr);
EXPECT_EQ(
boundCRS->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 "
"+x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris "
"+towgs84=-168,-60,320,0,0,0,0 +units=m "
"+vunits=m +no_defs +type=crs");
res);
auto compoundCRS = nn_dynamic_pointer_cast<CompoundCRS>(res);
ASSERT_TRUE(compoundCRS != nullptr);
EXPECT_EQ(compoundCRS->exportToPROJString(
PROJStringFormatter::create().get()),
"+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 "
"+x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris "
"+towgs84=-168,-60,320,0,0,0,0 +units=m "
"+geoidgrids=fr_ign_RAF18.tif +vunits=m +no_defs +type=crs");
}
{
auto crs = createVerticalCRS();
Expand All @@ -5478,6 +5497,13 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
CoordinateOperationContext::IntermediateCRSUse::NEVER),
crs);
}
{
auto crs = createCompoundCRS();
EXPECT_EQ(crs->createBoundCRSToWGS84IfPossible(
dbContext,
CoordinateOperationContext::IntermediateCRSUse::NEVER),
crs);
}
{
auto factoryIGNF =
AuthorityFactory::create(DatabaseContext::create(), "IGNF");
Expand Down