aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--include/proj/coordinateoperation.hpp4
-rw-r--r--src/iso19111/operation/conversion.cpp103
-rw-r--r--src/iso19111/operation/coordinateoperationfactory.cpp38
-rw-r--r--src/iso19111/operation/parammappings.cpp2
-rw-r--r--src/proj_constants.h4
-rw-r--r--test/unit/test_operationfactory.cpp23
6 files changed, 165 insertions, 9 deletions
diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp
index 3d2ab800..fc1d5b8a 100644
--- a/include/proj/coordinateoperation.hpp
+++ b/include/proj/coordinateoperation.hpp
@@ -1385,6 +1385,10 @@ class PROJ_GCC_DLL Conversion : public SingleOperation {
createGeographicGeocentric(const crs::CRSNNPtr &sourceCRS,
const crs::CRSNNPtr &targetCRS);
+ PROJ_INTERNAL static ConversionNNPtr
+ createGeographicGeocentricLatitude(const crs::CRSNNPtr &sourceCRS,
+ const crs::CRSNNPtr &targetCRS);
+
//! @endcond
protected:
diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp
index c1e5cb44..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
// ---------------------------------------------------------------------------
@@ -2403,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()),
@@ -2509,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())));
}
@@ -3447,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) {
diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
index 60365713..1b1cae9b 100644
--- a/src/iso19111/operation/coordinateoperationfactory.cpp
+++ b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -2993,13 +2993,6 @@ CoordinateOperationFactory::Private::createOperations(
}
}
- // Special case if both CRS are geodetic
- if (geodSrc && geodDst && !derivedSrc && !derivedDst) {
- createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc,
- geodDst, res);
- return res;
- }
-
if (geodSrc && geodSrc->isSphericalPlanetocentric()) {
createOperationsFromSphericalPlanetocentric(sourceCRS, targetCRS,
context, geodSrc, res);
@@ -3008,6 +3001,13 @@ CoordinateOperationFactory::Private::createOperations(
return applyInverse(createOperations(targetCRS, sourceCRS, context));
}
+ // Special case if both CRS are geodetic
+ if (geodSrc && geodDst && !derivedSrc && !derivedDst) {
+ createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc,
+ geodDst, res);
+ return res;
+ }
+
if (boundSrc) {
auto geodSrcBase = util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
boundSrc->baseCRS());
@@ -3940,6 +3940,24 @@ void CoordinateOperationFactory::Private::
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());
@@ -3953,7 +3971,8 @@ void CoordinateOperationFactory::Private::
cs::EllipsoidalCS::createLatitudeLongitude(
common::UnitOfMeasure::DEGREE)));
- auto opFirst = createGeodToGeodPROJBased(sourceCRS, interm_crs);
+ auto opFirst =
+ Conversion::createGeographicGeocentricLatitude(sourceCRS, interm_crs);
auto opsSecond = createOperations(interm_crs, targetCRS, context);
for (const auto &opSecond : opsSecond) {
try {
@@ -3999,7 +4018,8 @@ void CoordinateOperationFactory::Private::
auto intermBoundCRS =
crs::BoundCRS::create(intermGeog, boundSrc->hubCRS(), transf);
- auto opFirst = createGeodToGeodPROJBased(geodSrcBase, intermGeog);
+ auto opFirst =
+ Conversion::createGeographicGeocentricLatitude(geodSrcBase, intermGeog);
setCRSs(opFirst.get(), sourceCRS, intermBoundCRS);
auto opsSecond = createOperations(intermBoundCRS, targetCRS, context);
for (const auto &opSecond : opsSecond) {
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/proj_constants.h b/src/proj_constants.h
index f67a0583..5985f1b1 100644
--- a/src/proj_constants.h
+++ b/src/proj_constants.h
@@ -693,5 +693,9 @@
#define EPSG_NAME_METHOD_GEOGRAPHIC_TOPOCENTRIC "Geographic/topocentric conversions"
#define EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC 9837
+/* ------------------------------------------------------------------------ */
+
+#define PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE \
+ "Geographic latitude / Geocentric latitude"
#endif /* PROJ_CONSTANTS_INCLUDED */
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
index ec7d8700..b50542c6 100644
--- a/test/unit/test_operationfactory.cpp
+++ b/test/unit/test_operationfactory.cpp
@@ -6873,6 +6873,29 @@ TEST(operation, createOperation_geographic_to_spherical_ocentric) {
// ---------------------------------------------------------------------------
+TEST(operation, createOperation_spherical_ocentric_to_geocentric) {
+ auto objSrc = PROJStringParser().createFromPROJString(
+ "+proj=longlat +geoc +datum=WGS84 +type=crs");
+ auto src = nn_dynamic_pointer_cast<GeodeticCRS>(objSrc);
+ ASSERT_TRUE(src != nullptr);
+
+ auto objDest = PROJStringParser().createFromPROJString(
+ "+proj=geocent +datum=WGS84 +type=crs");
+ auto dest = nn_dynamic_pointer_cast<GeodeticCRS>(objDest);
+ ASSERT_TRUE(dest != nullptr);
+
+ auto op = CoordinateOperationFactory::create()->createOperation(
+ NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest));
+ ASSERT_TRUE(op != nullptr);
+ EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+ "+proj=pipeline "
+ "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+ "+step +inv +proj=geoc +ellps=WGS84 "
+ "+step +proj=cart +ellps=WGS84");
+}
+
+// ---------------------------------------------------------------------------
+
TEST(operation, createOperation_bound_of_spherical_ocentric_to_same_type) {
auto objSrc = PROJStringParser().createFromPROJString(
"+proj=longlat +geoc +ellps=GRS80 +towgs84=1,2,3 +type=crs");
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 11:47:22 +0000