aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/iso19111/coordinateoperation.cpp30
-rw-r--r--src/iso19111/crs.cpp104
-rw-r--r--test/unit/test_crs.cpp50
3 files changed, 148 insertions, 36 deletions
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 6593d3e9..6fcf4d30 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -8702,15 +8702,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,
bool allowInverse) {
@@ -8765,6 +8756,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) {
util::PropertyMap map;
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index e96b3cc9..ecbd39e1 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -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) {
@@ -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 {
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 0e340560..e8758932 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -5453,23 +5453,42 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
"+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();
@@ -5479,6 +5498,13 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
crs);
}
{
+ auto crs = createCompoundCRS();
+ EXPECT_EQ(crs->createBoundCRSToWGS84IfPossible(
+ dbContext,
+ CoordinateOperationContext::IntermediateCRSUse::NEVER),
+ crs);
+ }
+ {
auto factoryIGNF =
AuthorityFactory::create(DatabaseContext::create(), "IGNF");
auto crs = factoryIGNF->createCoordinateReferenceSystem("TERA50STEREO");
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 22:09:36 +0000