aboutsummaryrefslogtreecommitdiff
path: root/src/iso19111
diff options
context:
space:
mode:
Diffstat (limited to 'src/iso19111')
-rw-r--r--src/iso19111/coordinatesystem.cpp21
-rw-r--r--src/iso19111/crs.cpp176
-rw-r--r--src/iso19111/io.cpp185
-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
7 files changed, 543 insertions, 149 deletions
diff --git a/src/iso19111/coordinatesystem.cpp b/src/iso19111/coordinatesystem.cpp
index 498e3035..f9db5406 100644
--- a/src/iso19111/coordinatesystem.cpp
+++ b/src/iso19111/coordinatesystem.cpp
@@ -667,6 +667,27 @@ SphericalCSNNPtr SphericalCS::create(const util::PropertyMap &properties,
// ---------------------------------------------------------------------------
+/** \brief Instantiate a SphericalCS with 2 axis.
+ *
+ * This is an extension to ISO19111 to support (planet)-ocentric CS with
+ * geocentric latitude.
+ *
+ * @param properties See \ref general_properties.
+ * @param axis1 The first axis.
+ * @param axis2 The second axis.
+ * @return a new SphericalCS.
+ */
+SphericalCSNNPtr SphericalCS::create(const util::PropertyMap &properties,
+ const CoordinateSystemAxisNNPtr &axis1,
+ const CoordinateSystemAxisNNPtr &axis2) {
+ std::vector<CoordinateSystemAxisNNPtr> axis{axis1, axis2};
+ auto cs(SphericalCS::nn_make_shared<SphericalCS>(axis));
+ cs->setProperties(properties);
+ return cs;
+}
+
+// ---------------------------------------------------------------------------
+
//! @cond Doxygen_Suppress
EllipsoidalCS::~EllipsoidalCS() = default;
//! @endcond
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 7c8fcd81..be593878 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -1608,7 +1608,7 @@ GeodeticCRS::velocityModel() PROJ_PURE_DEFN {
// ---------------------------------------------------------------------------
-/** \brief Return whether the CRS is a geocentric one.
+/** \brief Return whether the CRS is a Cartesian geocentric one.
*
* A geocentric CRS is a geodetic CRS that has a Cartesian coordinate system
* with three axis, whose direction is respectively
@@ -1629,6 +1629,31 @@ bool GeodeticCRS::isGeocentric() PROJ_PURE_DEFN {
// ---------------------------------------------------------------------------
+/** \brief Return whether the CRS is a Spherical planetocentric one.
+ *
+ * A Spherical planetocentric CRS is a geodetic CRS that has a spherical
+ * (angular) coordinate system with 2 axis, which represent geocentric latitude/
+ * longitude or longitude/geocentric latitude.
+ *
+ * Such CRS are typically used in use case that apply to non-Earth bodies.
+ *
+ * @return true if the CRS is a Spherical planetocentric CRS.
+ *
+ * @since 8.2
+ */
+bool GeodeticCRS::isSphericalPlanetocentric() PROJ_PURE_DEFN {
+ const auto &cs = coordinateSystem();
+ const auto &axisList = cs->axisList();
+ return axisList.size() == 2 &&
+ dynamic_cast<cs::SphericalCS *>(cs.get()) != nullptr &&
+ ((ci_equal(axisList[0]->nameStr(), "planetocentric latitude") &&
+ ci_equal(axisList[1]->nameStr(), "planetocentric longitude")) ||
+ (ci_equal(axisList[0]->nameStr(), "planetocentric longitude") &&
+ ci_equal(axisList[1]->nameStr(), "planetocentric latitude")));
+}
+
+// ---------------------------------------------------------------------------
+
/** \brief Instantiate a GeodeticCRS from a datum::GeodeticReferenceFrame and a
* cs::SphericalCS.
*
@@ -1971,6 +1996,61 @@ void GeodeticCRS::addGeocentricUnitConversionIntoPROJString(
// ---------------------------------------------------------------------------
//! @cond Doxygen_Suppress
+void GeodeticCRS::addAngularUnitConvertAndAxisSwap(
+ io::PROJStringFormatter *formatter) const {
+ const auto &axisList = coordinateSystem()->axisList();
+
+ formatter->addStep("unitconvert");
+ formatter->addParam("xy_in", "rad");
+ if (axisList.size() == 3 && !formatter->omitZUnitConversion()) {
+ formatter->addParam("z_in", "m");
+ }
+ {
+ const auto &unitHoriz = axisList[0]->unit();
+ const auto projUnit = unitHoriz.exportToPROJString();
+ if (projUnit.empty()) {
+ formatter->addParam("xy_out", unitHoriz.conversionToSI());
+ } else {
+ formatter->addParam("xy_out", projUnit);
+ }
+ }
+ if (axisList.size() == 3 && !formatter->omitZUnitConversion()) {
+ const auto &unitZ = axisList[2]->unit();
+ auto projVUnit = unitZ.exportToPROJString();
+ if (projVUnit.empty()) {
+ formatter->addParam("z_out", unitZ.conversionToSI());
+ } else {
+ formatter->addParam("z_out", projVUnit);
+ }
+ }
+
+ const char *order[2] = {nullptr, nullptr};
+ const char *one = "1";
+ const char *two = "2";
+ for (int i = 0; i < 2; i++) {
+ const auto &dir = axisList[i]->direction();
+ if (&dir == &cs::AxisDirection::WEST) {
+ order[i] = "-1";
+ } else if (&dir == &cs::AxisDirection::EAST) {
+ order[i] = one;
+ } else if (&dir == &cs::AxisDirection::SOUTH) {
+ order[i] = "-2";
+ } else if (&dir == &cs::AxisDirection::NORTH) {
+ order[i] = two;
+ }
+ }
+ if (order[0] && order[1] && (order[0] != one || order[1] != two)) {
+ formatter->addStep("axisswap");
+ char orderStr[10];
+ sprintf(orderStr, "%.2s,%.2s", order[0], order[1]);
+ formatter->addParam("order", orderStr);
+ }
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
void GeodeticCRS::_exportToPROJString(
io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
{
@@ -1982,19 +2062,38 @@ void GeodeticCRS::_exportToPROJString(
return;
}
- if (!isGeocentric()) {
- io::FormattingException::Throw(
- "GeodeticCRS::exportToPROJString() only "
- "supports geocentric coordinate systems");
- }
+ if (isGeocentric()) {
+ if (!formatter->getCRSExport()) {
+ formatter->addStep("cart");
+ } else {
+ formatter->addStep("geocent");
+ }
- if (!formatter->getCRSExport()) {
- formatter->addStep("cart");
+ addDatumInfoToPROJString(formatter);
+ addGeocentricUnitConversionIntoPROJString(formatter);
+ } else if (isSphericalPlanetocentric()) {
+ if (!formatter->getCRSExport()) {
+ formatter->addStep("geoc");
+ addDatumInfoToPROJString(formatter);
+ addAngularUnitConvertAndAxisSwap(formatter);
+ } else {
+ io::FormattingException::Throw(
+ "GeodeticCRS::exportToPROJString() not supported on spherical "
+ "planetocentric coordinate systems");
+ // The below code now works as input to PROJ, but I'm not sure we
+ // want to propagate this, given that we got cs2cs doing conversion
+ // in the wrong direction in past versions.
+ /*formatter->addStep("longlat");
+ formatter->addParam("geoc");
+
+ addDatumInfoToPROJString(formatter);*/
+ }
} else {
- formatter->addStep("geocent");
+ io::FormattingException::Throw(
+ "GeodeticCRS::exportToPROJString() only "
+ "supports geocentric or spherical planetocentric "
+ "coordinate systems");
}
- addDatumInfoToPROJString(formatter);
- addGeocentricUnitConversionIntoPROJString(formatter);
}
//! @endcond
@@ -2859,61 +2958,6 @@ GeographicCRS::demoteTo2D(const std::string &newName,
// ---------------------------------------------------------------------------
//! @cond Doxygen_Suppress
-void GeographicCRS::addAngularUnitConvertAndAxisSwap(
- io::PROJStringFormatter *formatter) const {
- const auto &axisList = coordinateSystem()->axisList();
-
- formatter->addStep("unitconvert");
- formatter->addParam("xy_in", "rad");
- if (axisList.size() == 3 && !formatter->omitZUnitConversion()) {
- formatter->addParam("z_in", "m");
- }
- {
- const auto &unitHoriz = axisList[0]->unit();
- const auto projUnit = unitHoriz.exportToPROJString();
- if (projUnit.empty()) {
- formatter->addParam("xy_out", unitHoriz.conversionToSI());
- } else {
- formatter->addParam("xy_out", projUnit);
- }
- }
- if (axisList.size() == 3 && !formatter->omitZUnitConversion()) {
- const auto &unitZ = axisList[2]->unit();
- auto projVUnit = unitZ.exportToPROJString();
- if (projVUnit.empty()) {
- formatter->addParam("z_out", unitZ.conversionToSI());
- } else {
- formatter->addParam("z_out", projVUnit);
- }
- }
-
- const char *order[2] = {nullptr, nullptr};
- const char *one = "1";
- const char *two = "2";
- for (int i = 0; i < 2; i++) {
- const auto &dir = axisList[i]->direction();
- if (&dir == &cs::AxisDirection::WEST) {
- order[i] = "-1";
- } else if (&dir == &cs::AxisDirection::EAST) {
- order[i] = one;
- } else if (&dir == &cs::AxisDirection::SOUTH) {
- order[i] = "-2";
- } else if (&dir == &cs::AxisDirection::NORTH) {
- order[i] = two;
- }
- }
- if (order[0] && order[1] && (order[0] != one || order[1] != two)) {
- formatter->addStep("axisswap");
- char orderStr[10];
- sprintf(orderStr, "%.2s,%.2s", order[0], order[1]);
- formatter->addParam("order", orderStr);
- }
-}
-//! @endcond
-
-// ---------------------------------------------------------------------------
-
-//! @cond Doxygen_Suppress
void GeographicCRS::_exportToPROJString(
io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
{
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index a1b06ae7..5d7e951e 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -2798,7 +2798,11 @@ WKTParser::Private::buildCS(const WKTNodeNNPtr &node, /* maybe null */
return VerticalCS::create(csMap, axisList[0]);
}
} else if (ci_equal(csType, "spherical")) {
- if (axisCount == 3) {
+ if (axisCount == 2) {
+ // Extension to ISO19111 to support (planet)-ocentric CS with
+ // geocentric latitude
+ return SphericalCS::create(csMap, axisList[0], axisList[1]);
+ } else if (axisCount == 3) {
return SphericalCS::create(csMap, axisList[0], axisList[1],
axisList[2]);
}
@@ -5945,62 +5949,67 @@ CoordinateSystemNNPtr JSONParser::buildCS(const json &j) {
axisList.emplace_back(buildAxis(axis));
}
const PropertyMap &csMap = emptyPropertyMap;
+ const auto axisCount = axisList.size();
if (subtype == "ellipsoidal") {
- if (axisList.size() == 2) {
+ if (axisCount == 2) {
return EllipsoidalCS::create(csMap, axisList[0], axisList[1]);
}
- if (axisList.size() == 3) {
+ if (axisCount == 3) {
return EllipsoidalCS::create(csMap, axisList[0], axisList[1],
axisList[2]);
}
throw ParsingException("Expected 2 or 3 axis");
}
if (subtype == "Cartesian") {
- if (axisList.size() == 2) {
+ if (axisCount == 2) {
return CartesianCS::create(csMap, axisList[0], axisList[1]);
}
- if (axisList.size() == 3) {
+ if (axisCount == 3) {
return CartesianCS::create(csMap, axisList[0], axisList[1],
axisList[2]);
}
throw ParsingException("Expected 2 or 3 axis");
}
if (subtype == "vertical") {
- if (axisList.size() == 1) {
+ if (axisCount == 1) {
return VerticalCS::create(csMap, axisList[0]);
}
throw ParsingException("Expected 1 axis");
}
if (subtype == "spherical") {
- if (axisList.size() == 3) {
+ if (axisCount == 2) {
+ // Extension to ISO19111 to support (planet)-ocentric CS with
+ // geocentric latitude
+ return SphericalCS::create(csMap, axisList[0], axisList[1]);
+ } else if (axisCount == 3) {
return SphericalCS::create(csMap, axisList[0], axisList[1],
axisList[2]);
}
- throw ParsingException("Expected 3 axis");
+ throw ParsingException("Expected 2 or 3 axis");
}
if (subtype == "ordinal") {
return OrdinalCS::create(csMap, axisList);
}
if (subtype == "parametric") {
- if (axisList.size() == 1) {
+ if (axisCount == 1) {
return ParametricCS::create(csMap, axisList[0]);
}
throw ParsingException("Expected 1 axis");
}
if (subtype == "TemporalDateTime") {
- if (axisList.size() == 1) {
+ if (axisCount == 1) {
return DateTimeTemporalCS::create(csMap, axisList[0]);
}
throw ParsingException("Expected 1 axis");
}
if (subtype == "TemporalCount") {
- if (axisList.size() == 1) {
+ if (axisCount == 1) {
return TemporalCountCS::create(csMap, axisList[0]);
}
throw ParsingException("Expected 1 axis");
}
if (subtype == "TemporalMeasure") {
- if (axisList.size() == 1) {
+ if (axisCount == 1) {
return TemporalMeasureCS::create(csMap, axisList[0]);
}
throw ParsingException("Expected 1 axis");
@@ -8630,10 +8639,10 @@ struct PROJStringParser::Private {
PrimeMeridianNNPtr buildPrimeMeridian(Step &step);
GeodeticReferenceFrameNNPtr buildDatum(Step &step,
const std::string &title);
- GeographicCRSNNPtr buildGeographicCRS(int iStep, int iUnitConvert,
- int iAxisSwap, bool ignorePROJAxis);
+ GeodeticCRSNNPtr buildGeodeticCRS(int iStep, int iUnitConvert,
+ int iAxisSwap, bool ignorePROJAxis);
GeodeticCRSNNPtr buildGeocentricCRS(int iStep, int iUnitConvert);
- CRSNNPtr buildProjectedCRS(int iStep, GeographicCRSNNPtr geogCRS,
+ CRSNNPtr buildProjectedCRS(int iStep, const GeodeticCRSNNPtr &geogCRS,
int iUnitConvert, int iAxisSwap);
CRSNNPtr buildBoundOrCompoundCRSIfNeeded(int iStep, CRSNNPtr crs);
UnitOfMeasure buildUnit(Step &step, const std::string &unitsParamName,
@@ -8647,6 +8656,9 @@ struct PROJStringParser::Private {
EllipsoidalCSNNPtr buildEllipsoidalCS(int iStep, int iUnitConvert,
int iAxisSwap, bool ignorePROJAxis);
+
+ SphericalCSNNPtr buildSphericalCS(int iStep, int iUnitConvert,
+ int iAxisSwap, bool ignorePROJAxis);
};
// ---------------------------------------------------------------------------
@@ -9273,10 +9285,13 @@ PROJStringParser::Private::processAxisSwap(Step &step,
assert(iAxisSwap < 0 || ci_equal(steps_[iAxisSwap].name, "axisswap"));
const bool isGeographic = unit.type() == UnitOfMeasure::Type::ANGULAR;
+ const bool isSpherical = isGeographic && hasParamValue(step, "geoc");
const auto &eastName =
- isGeographic ? AxisName::Longitude : AxisName::Easting;
- const auto &eastAbbev =
- isGeographic ? AxisAbbreviation::lon : AxisAbbreviation::E;
+ isSpherical ? "Planetocentric longitude"
+ : isGeographic ? AxisName::Longitude : AxisName::Easting;
+ const auto &eastAbbev = isSpherical ? "V"
+ : isGeographic ? AxisAbbreviation::lon
+ : AxisAbbreviation::E;
const auto &eastDir = isGeographic
? AxisDirection::EAST
: (axisType == AxisType::NORTH_POLE)
@@ -9292,9 +9307,11 @@ PROJStringParser::Private::processAxisSwap(Step &step,
: nullMeridian);
const auto &northName =
- isGeographic ? AxisName::Latitude : AxisName::Northing;
- const auto &northAbbev =
- isGeographic ? AxisAbbreviation::lat : AxisAbbreviation::N;
+ isSpherical ? "Planetocentric latitude"
+ : isGeographic ? AxisName::Latitude : AxisName::Northing;
+ const auto &northAbbev = isSpherical ? "U"
+ : isGeographic ? AxisAbbreviation::lat
+ : AxisAbbreviation::N;
const auto &northDir = isGeographic
? AxisDirection::NORTH
: (axisType == AxisType::NORTH_POLE)
@@ -9311,15 +9328,19 @@ PROJStringParser::Private::processAxisSwap(Step &step,
.as_nullable()
: nullMeridian);
- CoordinateSystemAxisNNPtr west =
- createAxis(isGeographic ? AxisName::Longitude : AxisName::Westing,
- isGeographic ? AxisAbbreviation::lon : std::string(),
- AxisDirection::WEST, unit);
+ CoordinateSystemAxisNNPtr west = createAxis(
+ isSpherical ? "Planetocentric longitude"
+ : isGeographic ? AxisName::Longitude : AxisName::Westing,
+ isSpherical ? "V"
+ : isGeographic ? AxisAbbreviation::lon : std::string(),
+ AxisDirection::WEST, unit);
- CoordinateSystemAxisNNPtr south =
- createAxis(isGeographic ? AxisName::Latitude : AxisName::Southing,
- isGeographic ? AxisAbbreviation::lat : std::string(),
- AxisDirection::SOUTH, unit);
+ CoordinateSystemAxisNNPtr south = createAxis(
+ isSpherical ? "Planetocentric latitude"
+ : isGeographic ? AxisName::Latitude : AxisName::Southing,
+ isSpherical ? "U"
+ : isGeographic ? AxisAbbreviation::lat : std::string(),
+ AxisDirection::SOUTH, unit);
std::vector<CoordinateSystemAxisNNPtr> axis{east, north};
@@ -9419,6 +9440,42 @@ EllipsoidalCSNNPtr PROJStringParser::Private::buildEllipsoidalCS(
// ---------------------------------------------------------------------------
+SphericalCSNNPtr PROJStringParser::Private::buildSphericalCS(
+ int iStep, int iUnitConvert, int iAxisSwap, bool ignorePROJAxis) {
+ auto &step = steps_[iStep];
+ assert(iUnitConvert < 0 ||
+ ci_equal(steps_[iUnitConvert].name, "unitconvert"));
+
+ UnitOfMeasure angularUnit = UnitOfMeasure::DEGREE;
+ if (iUnitConvert >= 0) {
+ auto &stepUnitConvert = steps_[iUnitConvert];
+ const std::string *xy_in = &getParamValue(stepUnitConvert, "xy_in");
+ const std::string *xy_out = &getParamValue(stepUnitConvert, "xy_out");
+ if (stepUnitConvert.inverted) {
+ std::swap(xy_in, xy_out);
+ }
+ if (iUnitConvert < iStep) {
+ std::swap(xy_in, xy_out);
+ }
+ if (xy_in->empty() || xy_out->empty() || *xy_in != "rad" ||
+ (*xy_out != "rad" && *xy_out != "deg" && *xy_out != "grad")) {
+ throw ParsingException("unhandled values for xy_in and/or xy_out");
+ }
+ if (*xy_out == "rad") {
+ angularUnit = UnitOfMeasure::RADIAN;
+ } else if (*xy_out == "grad") {
+ angularUnit = UnitOfMeasure::GRAD;
+ }
+ }
+
+ std::vector<CoordinateSystemAxisNNPtr> axis = processAxisSwap(
+ step, angularUnit, iAxisSwap, AxisType::REGULAR, ignorePROJAxis);
+
+ return SphericalCS::create(emptyPropertyMap, axis[0], axis[1]);
+}
+
+// ---------------------------------------------------------------------------
+
static double getNumericValue(const std::string &paramValue,
bool *pHasError = nullptr) {
try {
@@ -9438,7 +9495,7 @@ namespace {
template <class T> inline void ignoreRetVal(T) {}
} // namespace
-GeographicCRSNNPtr PROJStringParser::Private::buildGeographicCRS(
+GeodeticCRSNNPtr PROJStringParser::Private::buildGeodeticCRS(
int iStep, int iUnitConvert, int iAxisSwap, bool ignorePROJAxis) {
auto &step = steps_[iStep];
@@ -9453,8 +9510,6 @@ GeographicCRSNNPtr PROJStringParser::Private::buildGeographicCRS(
auto props = PropertyMap().set(IdentifiedObject::NAME_KEY,
title.empty() ? "unknown" : title);
- auto cs =
- buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, ignorePROJAxis);
if (l_isGeographicStep &&
(hasUnusedParameters(step) ||
@@ -9463,7 +9518,17 @@ GeographicCRSNNPtr PROJStringParser::Private::buildGeographicCRS(
}
props.set("IMPLICIT_CS", true);
- return GeographicCRS::create(props, datum, cs);
+ if (!hasParamValue(step, "geoc")) {
+ auto cs =
+ buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, ignorePROJAxis);
+
+ return GeographicCRS::create(props, datum, cs);
+ } else {
+ auto cs =
+ buildSphericalCS(iStep, iUnitConvert, iAxisSwap, ignorePROJAxis);
+
+ return GeodeticCRS::create(props, datum, cs);
+ }
}
// ---------------------------------------------------------------------------
@@ -9614,8 +9679,10 @@ static bool is_in_stringlist(const std::string &str, const char *stringlist) {
// ---------------------------------------------------------------------------
-CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
- int iStep, GeographicCRSNNPtr geogCRS, int iUnitConvert, int iAxisSwap) {
+CRSNNPtr
+PROJStringParser::Private::buildProjectedCRS(int iStep,
+ const GeodeticCRSNNPtr &geodCRS,
+ int iUnitConvert, int iAxisSwap) {
auto &step = steps_[iStep];
const auto mappings = getMappingsFromPROJName(step.name);
const MethodMapping *mapping = mappings.empty() ? nullptr : mappings[0];
@@ -9641,7 +9708,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
}
if (mapping) {
- mapping = selectSphericalOrEllipsoidal(mapping, geogCRS);
+ mapping = selectSphericalOrEllipsoidal(mapping, geodCRS);
}
assert(isProjectedStep(step.name));
@@ -9651,7 +9718,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
const auto &title = title_;
if (!buildPrimeMeridian(step)->longitude()._isEquivalentTo(
- geogCRS->primeMeridian()->longitude(),
+ geodCRS->primeMeridian()->longitude(),
util::IComparable::Criterion::EQUIVALENT)) {
throw ParsingException("inconsistent pm values between projectedCRS "
"and its base geographicalCRS");
@@ -9917,7 +9984,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
}
if (k >= 0 && k <= 1) {
const double es =
- geogCRS->ellipsoid()->squaredEccentricity();
+ geodCRS->ellipsoid()->squaredEccentricity();
if (es < 0 || es == 1) {
throw ParsingException("Invalid flattening");
}
@@ -10028,10 +10095,16 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
{"PROJ ob_tran o_proj=longlat", "PROJ ob_tran o_proj=lonlat",
"PROJ ob_tran o_proj=latlon", "PROJ ob_tran o_proj=latlong"}) {
if (starts_with(methodName, substr)) {
- return DerivedGeographicCRS::create(
- PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"),
- geogCRS, NN_NO_CHECK(conv),
- buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, false));
+ auto geogCRS =
+ util::nn_dynamic_pointer_cast<GeographicCRS>(geodCRS);
+ if (geogCRS) {
+ return DerivedGeographicCRS::create(
+ PropertyMap().set(IdentifiedObject::NAME_KEY,
+ "unnamed"),
+ NN_NO_CHECK(geogCRS), NN_NO_CHECK(conv),
+ buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap,
+ false));
+ }
}
}
}
@@ -10039,11 +10112,11 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
std::vector<CoordinateSystemAxisNNPtr> axis =
processAxisSwap(step, unit, iAxisSwap, axisType, false);
- auto csGeogCRS = geogCRS->coordinateSystem();
- auto cs = csGeogCRS->axisList().size() == 2
+ auto csGeodCRS = geodCRS->coordinateSystem();
+ auto cs = csGeodCRS->axisList().size() == 2
? CartesianCS::create(emptyPropertyMap, axis[0], axis[1])
: CartesianCS::create(emptyPropertyMap, axis[0], axis[1],
- csGeogCRS->axisList()[2]);
+ csGeodCRS->axisList()[2]);
auto props = PropertyMap().set(IdentifiedObject::NAME_KEY,
title.empty() ? "unknown" : title);
@@ -10057,7 +10130,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
bWebMercator
? createPseudoMercator(
props.set(IdentifiedObject::NAME_KEY, webMercatorName), cs)
- : ProjectedCRS::create(props, geogCRS, NN_NO_CHECK(conv), cs);
+ : ProjectedCRS::create(props, geodCRS, NN_NO_CHECK(conv), cs);
return crs;
}
@@ -10528,8 +10601,8 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
// First run is dry run to mark all recognized/unrecognized tokens
for (int iter = 0; iter < 2; iter++) {
auto obj = d->buildBoundOrCompoundCRSIfNeeded(
- 0, d->buildGeographicCRS(iFirstGeogStep, iFirstUnitConvert,
- iFirstAxisSwap, false));
+ 0, d->buildGeodeticCRS(iFirstGeogStep, iFirstUnitConvert,
+ iFirstAxisSwap, false));
if (iter == 1) {
return nn_static_pointer_cast<BaseObject>(obj);
}
@@ -10546,14 +10619,14 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
iProjStep,
d->buildProjectedCRS(
iProjStep,
- d->buildGeographicCRS(iFirstGeogStep,
- iFirstUnitConvert < iFirstGeogStep
- ? iFirstUnitConvert
- : -1,
- iFirstAxisSwap < iFirstGeogStep
- ? iFirstAxisSwap
- : -1,
- true),
+ d->buildGeodeticCRS(iFirstGeogStep,
+ iFirstUnitConvert < iFirstGeogStep
+ ? iFirstUnitConvert
+ : -1,
+ iFirstAxisSwap < iFirstGeogStep
+ ? iFirstAxisSwap
+ : -1,
+ true),
iFirstUnitConvert < iFirstGeogStep ? iSecondUnitConvert
: iFirstUnitConvert,
iFirstAxisSwap < iFirstGeogStep ? iSecondAxisSwap
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:53:19 +0000