Skip to content

Commit

Permalink
Merge pull request #2365 from rouault/bound_crs_to_wgs84_vertcrs
Browse files Browse the repository at this point in the history 
proj_crs_create_bound_crs_to_WGS84(): make it work on verticalCRS/compoundCRS such as EPSG:4326+5773 and EPSG:4326+3855
  • Loading branch information
@rouault
rouault authored Oct 5, 2020
2 parents 390d1ee + a227423 commit f9c1739
Show file tree
Hide file tree
Showing 3 changed files with 148 additions and 36 deletions.
30 changes: 21 additions & 9 deletions src/iso19111/coordinateoperation.cpp
View file
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
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
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

0 comments on commit f9c1739

Please sign in to comment.