aboutsummaryrefslogtreecommitdiff
path: root/src/iso19111/operation
diff options
context:
space:
mode:
Diffstat (limited to 'src/iso19111/operation')
-rw-r--r--src/iso19111/operation/conversion.cpp116
-rw-r--r--src/iso19111/operation/coordinateoperationfactory.cpp143
-rw-r--r--src/iso19111/operation/parammappings.cpp2
-rw-r--r--src/iso19111/operation/transformation.cpp49
4 files changed, 283 insertions, 27 deletions
diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp
index 3d09c9fa..f7bb8ec8 100644
--- a/src/iso19111/operation/conversion.cpp
+++ b/src/iso19111/operation/conversion.cpp
@@ -79,6 +79,7 @@ constexpr double UTM_SCALE_FACTOR = 0.9996;
constexpr double UTM_FALSE_EASTING = 500000.0;
constexpr double UTM_NORTH_FALSE_NORTHING = 0.0;
constexpr double UTM_SOUTH_FALSE_NORTHING = 10000000.0;
+
//! @endcond
// ---------------------------------------------------------------------------
@@ -264,7 +265,8 @@ createConversion(const util::PropertyMap &properties,
const std::vector<ParameterValueNNPtr> &values) {
std::vector<OperationParameterNNPtr> parameters;
- for (int i = 0; mapping->params[i] != nullptr; i++) {
+ for (int i = 0; mapping->params != nullptr && mapping->params[i] != nullptr;
+ i++) {
const auto *param = mapping->params[i];
auto paramProperties = util::PropertyMap().set(
common::IdentifiedObject::NAME_KEY, param->wkt2_name);
@@ -2402,6 +2404,28 @@ Conversion::createGeographicGeocentric(const crs::CRSNNPtr &sourceCRS,
// ---------------------------------------------------------------------------
+/** \brief Instantiate a conversion between a GeographicCRS and a spherical
+ * planetocentric GeodeticCRS
+ *
+ * This method peforms conversion between geodetic latitude and geocentric
+ * latitude
+ *
+ * @return a new Conversion.
+ */
+ConversionNNPtr
+Conversion::createGeographicGeocentricLatitude(const crs::CRSNNPtr &sourceCRS,
+ const crs::CRSNNPtr &targetCRS) {
+ auto properties = util::PropertyMap().set(
+ common::IdentifiedObject::NAME_KEY,
+ buildOpName("Conversion", sourceCRS, targetCRS));
+ auto conv = create(
+ properties, PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE, {});
+ conv->setCRSs(sourceCRS, targetCRS, nullptr);
+ return conv;
+}
+
+// ---------------------------------------------------------------------------
+
InverseConversion::InverseConversion(const ConversionNNPtr &forward)
: Conversion(
OperationMethod::create(createPropertiesForInverse(forward->method()),
@@ -2508,6 +2532,15 @@ CoordinateOperationNNPtr Conversion::inverse() const {
return conv;
}
+ if (method()->nameStr() ==
+ PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE) {
+ auto conv =
+ create(createPropertiesForInverse(this, false, false),
+ PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE, {});
+ conv->setCRSs(this, true);
+ return conv;
+ }
+
return InverseConversion::create(NN_NO_CHECK(
util::nn_dynamic_pointer_cast<Conversion>(shared_from_this())));
}
@@ -3446,6 +3479,77 @@ void Conversion::_exportToPROJString(
auto l_sourceCRS = sourceCRS();
auto l_targetCRS = targetCRS();
+ if (methodName == PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE) {
+
+ const auto extractGeodeticCRSIfGeodeticCRSOrEquivalent =
+ [](const crs::CRSPtr &crs) {
+ auto geodCRS = std::dynamic_pointer_cast<crs::GeodeticCRS>(crs);
+ if (!geodCRS) {
+ auto compoundCRS =
+ std::dynamic_pointer_cast<crs::CompoundCRS>(crs);
+ if (compoundCRS) {
+ const auto &components =
+ compoundCRS->componentReferenceSystems();
+ if (!components.empty()) {
+ geodCRS =
+ util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
+ components[0]);
+ if (!geodCRS) {
+ auto boundCRS = util::nn_dynamic_pointer_cast<
+ crs::BoundCRS>(components[0]);
+ if (boundCRS) {
+ geodCRS = util::nn_dynamic_pointer_cast<
+ crs::GeodeticCRS>(boundCRS->baseCRS());
+ }
+ }
+ }
+ } else {
+ auto boundCRS =
+ std::dynamic_pointer_cast<crs::BoundCRS>(crs);
+ if (boundCRS) {
+ geodCRS =
+ util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
+ boundCRS->baseCRS());
+ }
+ }
+ }
+ return geodCRS;
+ };
+
+ auto sourceCRSGeod = dynamic_cast<const crs::GeodeticCRS *>(
+ extractGeodeticCRSIfGeodeticCRSOrEquivalent(l_sourceCRS).get());
+ auto targetCRSGeod = dynamic_cast<const crs::GeodeticCRS *>(
+ extractGeodeticCRSIfGeodeticCRSOrEquivalent(l_targetCRS).get());
+ if (sourceCRSGeod && targetCRSGeod) {
+ auto sourceCRSGeog =
+ dynamic_cast<const crs::GeographicCRS *>(sourceCRSGeod);
+ auto targetCRSGeog =
+ dynamic_cast<const crs::GeographicCRS *>(targetCRSGeod);
+ bool isSrcGeocentricLat =
+ sourceCRSGeod->isSphericalPlanetocentric();
+ bool isSrcGeographic = sourceCRSGeog != nullptr;
+ bool isTargetGeocentricLat =
+ targetCRSGeod->isSphericalPlanetocentric();
+ bool isTargetGeographic = targetCRSGeog != nullptr;
+ if ((isSrcGeocentricLat && isTargetGeographic) ||
+ (isSrcGeographic && isTargetGeocentricLat)) {
+
+ formatter->startInversion();
+ sourceCRSGeod->_exportToPROJString(formatter);
+ formatter->stopInversion();
+
+ targetCRSGeod->_exportToPROJString(formatter);
+
+ return;
+ }
+ }
+
+ throw io::FormattingException("Invalid nature of source and/or "
+ "targetCRS for Geographic latitude / "
+ "Geocentric latitude"
+ "conversion");
+ }
+
crs::GeographicCRSPtr srcGeogCRS;
if (!formatter->getCRSExport() && l_sourceCRS && applySourceCRSModifiers) {
@@ -3459,11 +3563,15 @@ void Conversion::_exportToPROJString(
}
}
- srcGeogCRS = std::dynamic_pointer_cast<crs::GeographicCRS>(horiz);
- if (srcGeogCRS) {
+ auto srcGeodCRS = dynamic_cast<const crs::GeodeticCRS *>(horiz.get());
+ if (srcGeodCRS) {
+ srcGeogCRS = std::dynamic_pointer_cast<crs::GeographicCRS>(horiz);
+ }
+ if (srcGeodCRS &&
+ (srcGeogCRS || srcGeodCRS->isSphericalPlanetocentric())) {
formatter->setOmitProjLongLatIfPossible(true);
formatter->startInversion();
- srcGeogCRS->_exportToPROJString(formatter);
+ srcGeodCRS->_exportToPROJString(formatter);
formatter->stopInversion();
formatter->setOmitProjLongLatIfPossible(false);
}
diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
index fc64bd2e..1b1cae9b 100644
--- a/src/iso19111/operation/coordinateoperationfactory.cpp
+++ b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -558,6 +558,17 @@ struct CoordinateOperationFactory::Private {
const crs::GeodeticCRS *geodDst,
std::vector<CoordinateOperationNNPtr> &res);
+ static void createOperationsFromSphericalPlanetocentric(
+ const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+ Private::Context &context, const crs::GeodeticCRS *geodSrc,
+ std::vector<CoordinateOperationNNPtr> &res);
+
+ static void createOperationsFromBoundOfSphericalPlanetocentric(
+ const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+ Private::Context &context, const crs::BoundCRS *boundSrc,
+ const crs::GeodeticCRSNNPtr &geodSrcBase,
+ std::vector<CoordinateOperationNNPtr> &res);
+
static void createOperationsDerivedTo(
const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
Private::Context &context, const crs::DerivedCRS *derivedSrc,
@@ -2982,6 +2993,14 @@ CoordinateOperationFactory::Private::createOperations(
}
}
+ if (geodSrc && geodSrc->isSphericalPlanetocentric()) {
+ createOperationsFromSphericalPlanetocentric(sourceCRS, targetCRS,
+ context, geodSrc, res);
+ return res;
+ } else if (geodDst && geodDst->isSphericalPlanetocentric()) {
+ return applyInverse(createOperations(targetCRS, sourceCRS, context));
+ }
+
// Special case if both CRS are geodetic
if (geodSrc && geodDst && !derivedSrc && !derivedDst) {
createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc,
@@ -2989,6 +3008,24 @@ CoordinateOperationFactory::Private::createOperations(
return res;
}
+ if (boundSrc) {
+ auto geodSrcBase = util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
+ boundSrc->baseCRS());
+ if (geodSrcBase && geodSrcBase->isSphericalPlanetocentric()) {
+ createOperationsFromBoundOfSphericalPlanetocentric(
+ sourceCRS, targetCRS, context, boundSrc,
+ NN_NO_CHECK(geodSrcBase), res);
+ return res;
+ }
+ } else if (boundDst) {
+ auto geodDstBase = util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
+ boundDst->baseCRS());
+ if (geodDstBase && geodDstBase->isSphericalPlanetocentric()) {
+ return applyInverse(
+ createOperations(targetCRS, sourceCRS, context));
+ }
+ }
+
// If the source is a derived CRS, then chain the inverse of its
// deriving conversion, with transforms from its baseCRS to the
// targetCRS
@@ -3895,6 +3932,112 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod(
// ---------------------------------------------------------------------------
+void CoordinateOperationFactory::Private::
+ createOperationsFromSphericalPlanetocentric(
+ const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+ Private::Context &context, const crs::GeodeticCRS *geodSrc,
+ std::vector<CoordinateOperationNNPtr> &res) {
+
+ ENTER_FUNCTION();
+
+ const auto IsSameDatum = [&context,
+ &geodSrc](const crs::GeodeticCRS *geodDst) {
+ const auto &authFactory = context.context->getAuthorityFactory();
+ const auto dbContext =
+ authFactory ? authFactory->databaseContext().as_nullable()
+ : nullptr;
+
+ return geodSrc->datumNonNull(dbContext)->_isEquivalentTo(
+ geodDst->datumNonNull(dbContext).get(),
+ util::IComparable::Criterion::EQUIVALENT);
+ };
+ auto geogDst = dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
+ if (geogDst && IsSameDatum(geogDst)) {
+ res.emplace_back(Conversion::createGeographicGeocentricLatitude(
+ sourceCRS, targetCRS));
+ return;
+ }
+
+ // Create an intermediate geographic CRS with the same datum as the
+ // source spherical planetocentric one
+ std::string interm_crs_name(geodSrc->nameStr());
+ interm_crs_name += " (geographic)";
+ auto interm_crs =
+ util::nn_static_pointer_cast<crs::CRS>(crs::GeographicCRS::create(
+ addDomains(util::PropertyMap().set(
+ common::IdentifiedObject::NAME_KEY, interm_crs_name),
+ geodSrc),
+ geodSrc->datum(), geodSrc->datumEnsemble(),
+ cs::EllipsoidalCS::createLatitudeLongitude(
+ common::UnitOfMeasure::DEGREE)));
+
+ auto opFirst =
+ Conversion::createGeographicGeocentricLatitude(sourceCRS, interm_crs);
+ auto opsSecond = createOperations(interm_crs, targetCRS, context);
+ for (const auto &opSecond : opsSecond) {
+ try {
+ res.emplace_back(ConcatenatedOperation::createComputeMetadata(
+ {opFirst, opSecond}, disallowEmptyIntersection));
+ } catch (const InvalidOperationEmptyIntersection &) {
+ }
+ }
+}
+
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::
+ createOperationsFromBoundOfSphericalPlanetocentric(
+ const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+ Private::Context &context, const crs::BoundCRS *boundSrc,
+ const crs::GeodeticCRSNNPtr &geodSrcBase,
+ std::vector<CoordinateOperationNNPtr> &res) {
+
+ ENTER_FUNCTION();
+
+ // Create an intermediate geographic CRS with the same datum as the
+ // source spherical planetocentric one
+ std::string interm_crs_name(geodSrcBase->nameStr());
+ interm_crs_name += " (geographic)";
+ auto intermGeog =
+ util::nn_static_pointer_cast<crs::CRS>(crs::GeographicCRS::create(
+ addDomains(util::PropertyMap().set(
+ common::IdentifiedObject::NAME_KEY, interm_crs_name),
+ geodSrcBase.get()),
+ geodSrcBase->datum(), geodSrcBase->datumEnsemble(),
+ cs::EllipsoidalCS::createLatitudeLongitude(
+ common::UnitOfMeasure::DEGREE)));
+
+ // Create an intermediate boundCRS wrapping the above intermediate
+ // geographic CRS
+ auto transf = boundSrc->transformation()->shallowClone();
+ // keep a reference to the target before patching it with itself
+ // (this is due to our abuse of passing shared_ptr by reference
+ auto transfTarget = transf->targetCRS();
+ setCRSs(transf.get(), intermGeog, transfTarget);
+
+ auto intermBoundCRS =
+ crs::BoundCRS::create(intermGeog, boundSrc->hubCRS(), transf);
+
+ auto opFirst =
+ Conversion::createGeographicGeocentricLatitude(geodSrcBase, intermGeog);
+ setCRSs(opFirst.get(), sourceCRS, intermBoundCRS);
+ auto opsSecond = createOperations(intermBoundCRS, targetCRS, context);
+ for (const auto &opSecond : opsSecond) {
+ try {
+ auto opSecondClone = opSecond->shallowClone();
+ // In theory, we should not need that setCRSs() forcing, but due
+ // how BoundCRS transformations are implemented currently, we
+ // need it in practice.
+ setCRSs(opSecondClone.get(), intermBoundCRS, targetCRS);
+ res.emplace_back(ConcatenatedOperation::createComputeMetadata(
+ {opFirst, opSecondClone}, disallowEmptyIntersection));
+ } catch (const InvalidOperationEmptyIntersection &) {
+ }
+ }
+}
+
+// ---------------------------------------------------------------------------
+
void CoordinateOperationFactory::Private::createOperationsDerivedTo(
const crs::CRSNNPtr & /*sourceCRS*/, const crs::CRSNNPtr &targetCRS,
Private::Context &context, const crs::DerivedCRS *derivedSrc,
diff --git a/src/iso19111/operation/parammappings.cpp b/src/iso19111/operation/parammappings.cpp
index 9acc9e93..df5ae7c5 100644
--- a/src/iso19111/operation/parammappings.cpp
+++ b/src/iso19111/operation/parammappings.cpp
@@ -1340,6 +1340,8 @@ static const MethodMapping otherMethodMappings[] = {
{EPSG_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC,
EPSG_CODE_METHOD_GEOGRAPHIC_GEOCENTRIC, nullptr, nullptr, nullptr,
nullptr},
+ {PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE, 0, nullptr, nullptr,
+ nullptr, nullptr},
{EPSG_NAME_METHOD_LONGITUDE_ROTATION, EPSG_CODE_METHOD_LONGITUDE_ROTATION,
nullptr, nullptr, nullptr, paramsLongitudeRotation},
{EPSG_NAME_METHOD_AFFINE_PARAMETRIC_TRANSFORMATION,
diff --git a/src/iso19111/operation/transformation.cpp b/src/iso19111/operation/transformation.cpp
index 273b636f..a30e1248 100644
--- a/src/iso19111/operation/transformation.cpp
+++ b/src/iso19111/operation/transformation.cpp
@@ -474,13 +474,16 @@ static void getTransformationType(const crs::CRSNNPtr &sourceCRSIn,
dynamic_cast<const crs::GeographicCRS *>(sourceCRSIn.get());
auto targetCRSGeog =
dynamic_cast<const crs::GeographicCRS *>(targetCRSIn.get());
- if (!sourceCRSGeog || !targetCRSGeog) {
+ if (!(sourceCRSGeog ||
+ (sourceCRSGeod && sourceCRSGeod->isSphericalPlanetocentric())) ||
+ !(targetCRSGeog ||
+ (targetCRSGeod && targetCRSGeod->isSphericalPlanetocentric()))) {
throw InvalidOperation("Inconsistent CRS type");
}
const auto nSrcAxisCount =
- sourceCRSGeog->coordinateSystem()->axisList().size();
+ sourceCRSGeod->coordinateSystem()->axisList().size();
const auto nTargetAxisCount =
- targetCRSGeog->coordinateSystem()->axisList().size();
+ targetCRSGeod->coordinateSystem()->axisList().size();
isGeog2D = nSrcAxisCount == 2 && nTargetAxisCount == 2;
isGeog3D = !isGeog2D && nSrcAxisCount >= 2 && nTargetAxisCount >= 2;
}
@@ -1003,36 +1006,36 @@ TransformationNNPtr Transformation::createTOWGS84(
"Invalid number of elements in TOWGS84Parameters");
}
- crs::CRSPtr transformSourceCRS = sourceCRSIn->extractGeodeticCRS();
- if (!transformSourceCRS) {
+ auto transformSourceGeodCRS = sourceCRSIn->extractGeodeticCRS();
+ if (!transformSourceGeodCRS) {
throw InvalidOperation(
"Cannot find GeodeticCRS in sourceCRS of TOWGS84 transformation");
}
util::PropertyMap properties;
properties.set(common::IdentifiedObject::NAME_KEY,
- concat("Transformation from ", transformSourceCRS->nameStr(),
- " to WGS84"));
-
- auto targetCRS =
- dynamic_cast<const crs::GeographicCRS *>(transformSourceCRS.get())
- ? util::nn_static_pointer_cast<crs::CRS>(
- crs::GeographicCRS::EPSG_4326)
- : util::nn_static_pointer_cast<crs::CRS>(
- crs::GeodeticCRS::EPSG_4978);
-
+ concat("Transformation from ",
+ transformSourceGeodCRS->nameStr(), " to WGS84"));
+
+ auto targetCRS = dynamic_cast<const crs::GeographicCRS *>(
+ transformSourceGeodCRS.get()) ||
+ transformSourceGeodCRS->isSphericalPlanetocentric()
+ ? util::nn_static_pointer_cast<crs::CRS>(
+ crs::GeographicCRS::EPSG_4326)
+ : util::nn_static_pointer_cast<crs::CRS>(
+ crs::GeodeticCRS::EPSG_4978);
+
+ crs::CRSNNPtr transformSourceCRS = NN_NO_CHECK(transformSourceGeodCRS);
if (TOWGS84Parameters.size() == 3) {
return createGeocentricTranslations(
- properties, NN_NO_CHECK(transformSourceCRS), targetCRS,
- TOWGS84Parameters[0], TOWGS84Parameters[1], TOWGS84Parameters[2],
- {});
+ properties, transformSourceCRS, targetCRS, TOWGS84Parameters[0],
+ TOWGS84Parameters[1], TOWGS84Parameters[2], {});
}
- return createPositionVector(properties, NN_NO_CHECK(transformSourceCRS),
- targetCRS, TOWGS84Parameters[0],
- TOWGS84Parameters[1], TOWGS84Parameters[2],
- TOWGS84Parameters[3], TOWGS84Parameters[4],
- TOWGS84Parameters[5], TOWGS84Parameters[6], {});
+ return createPositionVector(
+ properties, transformSourceCRS, targetCRS, TOWGS84Parameters[0],
+ TOWGS84Parameters[1], TOWGS84Parameters[2], TOWGS84Parameters[3],
+ TOWGS84Parameters[4], TOWGS84Parameters[5], TOWGS84Parameters[6], {});
}
// ---------------------------------------------------------------------------
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 20:52:20 +0000