aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--include/proj/coordinateoperation.hpp4
-rw-r--r--include/proj/coordinatesystem.hpp5
-rw-r--r--include/proj/crs.hpp13
-rw-r--r--scripts/reference_exported_symbols.txt2
-rw-r--r--src/apps/cs2cs.cpp33
-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
-rw-r--r--src/proj_constants.h4
-rwxr-xr-xtest/cli/testvarious12
-rw-r--r--test/cli/tv_out.dist4
-rw-r--r--test/unit/test_io.cpp111
-rw-r--r--test/unit/test_operationfactory.cpp149
17 files changed, 844 insertions, 185 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/include/proj/coordinatesystem.hpp b/include/proj/coordinatesystem.hpp
index e1650168..b40b038d 100644
--- a/include/proj/coordinatesystem.hpp
+++ b/include/proj/coordinatesystem.hpp
@@ -314,6 +314,11 @@ class PROJ_GCC_DLL SphericalCS final : public CoordinateSystem {
const CoordinateSystemAxisNNPtr &axis2,
const CoordinateSystemAxisNNPtr &axis3);
+ PROJ_DLL static SphericalCSNNPtr
+ create(const util::PropertyMap &properties,
+ const CoordinateSystemAxisNNPtr &axis1,
+ const CoordinateSystemAxisNNPtr &axis2);
+
protected:
PROJ_INTERNAL explicit SphericalCS(
const std::vector<CoordinateSystemAxisNNPtr> &axisIn);
diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp
index dcab094a..5a8e75ae 100644
--- a/include/proj/crs.hpp
+++ b/include/proj/crs.hpp
@@ -270,6 +270,8 @@ class PROJ_GCC_DLL GeodeticCRS : virtual public SingleCRS,
PROJ_DLL bool isGeocentric() PROJ_PURE_DECL;
+ PROJ_DLL bool isSphericalPlanetocentric() PROJ_PURE_DECL;
+
PROJ_DLL static GeodeticCRSNNPtr
create(const util::PropertyMap &properties,
const datum::GeodeticReferenceFrameNNPtr &datum,
@@ -308,6 +310,9 @@ class PROJ_GCC_DLL GeodeticCRS : virtual public SingleCRS,
PROJ_INTERNAL void addGeocentricUnitConversionIntoPROJString(
io::PROJStringFormatter *formatter) const;
+ PROJ_INTERNAL void
+ addAngularUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter) const;
+
PROJ_INTERNAL void _exportToWKT(io::WKTFormatter *formatter)
const override; // throw(io::FormattingException)
@@ -400,12 +405,10 @@ class PROJ_GCC_DLL GeographicCRS : public GeodeticCRS {
PROJ_PRIVATE :
//! @cond Doxygen_Suppress
- PROJ_INTERNAL void
- addAngularUnitConvertAndAxisSwap(
- io::PROJStringFormatter *formatter) const;
- PROJ_INTERNAL void _exportToPROJString(io::PROJStringFormatter *formatter)
- const override; // throw(FormattingException)
+ PROJ_INTERNAL void
+ _exportToPROJString(io::PROJStringFormatter *formatter)
+ const override; // throw(FormattingException)
PROJ_INTERNAL void _exportToJSON(io::JSONFormatter *formatter)
const override; // throw(FormattingException)
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index d49a59b2..6cb734e2 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -152,6 +152,7 @@ osgeo::proj::crs::GeodeticCRS::ellipsoid() const
osgeo::proj::crs::GeodeticCRS::~GeodeticCRS()
osgeo::proj::crs::GeodeticCRS::identify(std::shared_ptr<osgeo::proj::io::AuthorityFactory> const&) const
osgeo::proj::crs::GeodeticCRS::isGeocentric() const
+osgeo::proj::crs::GeodeticCRS::isSphericalPlanetocentric() const
osgeo::proj::crs::GeodeticCRS::primeMeridian() const
osgeo::proj::crs::GeodeticCRS::velocityModel() const
osgeo::proj::crs::GeographicCRS::coordinateSystem() const
@@ -225,6 +226,7 @@ osgeo::proj::cs::OrdinalCS::create(osgeo::proj::util::PropertyMap const&, std::v
osgeo::proj::cs::OrdinalCS::~OrdinalCS()
osgeo::proj::cs::ParametricCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&)
osgeo::proj::cs::ParametricCS::~ParametricCS()
+osgeo::proj::cs::SphericalCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&)
osgeo::proj::cs::SphericalCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&)
osgeo::proj::cs::SphericalCS::~SphericalCS()
osgeo::proj::cs::TemporalCountCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&)
diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp
index 19f924d4..03fd4dba 100644
--- a/src/apps/cs2cs.cpp
+++ b/src/apps/cs2cs.cpp
@@ -59,10 +59,10 @@
static PJ *transformation = nullptr;
-static bool srcIsGeog = false;
+static bool srcIsLongLat = false;
static double srcToRadians = 0.0;
-static bool destIsGeog = false;
+static bool destIsLongLat = false;
static double destToRadians = 0.0;
static bool destIsLatLong = false;
@@ -158,7 +158,7 @@ static void process(FILE *fid)
if (data.u != HUGE_VAL) {
- if (srcIsGeog) {
+ if (srcIsLongLat) {
/* dmstor gives values to radians. Convert now to the SRS unit
*/
data.u /= srcToRadians;
@@ -179,7 +179,7 @@ static void process(FILE *fid)
if (data.u == HUGE_VAL) /* error output */
fputs(oterr, stdout);
- else if (destIsGeog && !oform) { /*ascii DMS output */
+ else if (destIsLongLat && !oform) { /*ascii DMS output */
// rtodms() expect radians: convert from the output SRS unit
data.u *= destToRadians;
@@ -206,7 +206,7 @@ static void process(FILE *fid)
}
} else { /* x-y or decimal degree ascii output */
- if (destIsGeog) {
+ if (destIsLongLat) {
data.v *= destToRadians * RAD_TO_DEG;
data.u *= destToRadians * RAD_TO_DEG;
}
@@ -239,7 +239,7 @@ static void process(FILE *fid)
/************************************************************************/
static PJ *instantiate_crs(const std::string &definition,
- bool &isGeog, double &toRadians,
+ bool &isLongLatCS, double &toRadians,
bool &isLatFirst) {
PJ *crs = proj_create(nullptr,
pj_add_type_crs_if_needed(definition).c_str());
@@ -247,7 +247,7 @@ static PJ *instantiate_crs(const std::string &definition,
return nullptr;
}
- isGeog = false;
+ isLongLatCS = false;
toRadians = 0.0;
isLatFirst = false;
@@ -259,11 +259,11 @@ static PJ *instantiate_crs(const std::string &definition,
type = proj_get_type(crs);
}
if (type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
- type == PJ_TYPE_GEOGRAPHIC_3D_CRS) {
+ type == PJ_TYPE_GEOGRAPHIC_3D_CRS ||
+ type == PJ_TYPE_GEODETIC_CRS) {
auto cs = proj_crs_get_coordinate_system(nullptr, crs);
assert(cs);
- isGeog = true;
const char *axisName = "";
proj_cs_get_axis_info(nullptr, cs, 0,
&axisName, // name,
@@ -277,6 +277,9 @@ static PJ *instantiate_crs(const std::string &definition,
isLatFirst =
NS_PROJ::internal::ci_find(std::string(axisName), "latitude") !=
std::string::npos;
+ isLongLatCS = isLatFirst ||
+ NS_PROJ::internal::ci_find(std::string(axisName), "longitude") !=
+ std::string::npos;
proj_destroy(cs);
}
@@ -736,7 +739,7 @@ int main(int argc, char **argv) {
PJ *src = nullptr;
if (!fromStr.empty()) {
bool ignored;
- src = instantiate_crs(fromStr, srcIsGeog,
+ src = instantiate_crs(fromStr, srcIsLongLat,
srcToRadians, ignored);
if (!src) {
emess(3, "cannot instantiate source coordinate system");
@@ -745,7 +748,7 @@ int main(int argc, char **argv) {
PJ *dst = nullptr;
if (!toStr.empty()) {
- dst = instantiate_crs(toStr, destIsGeog,
+ dst = instantiate_crs(toStr, destIsLongLat,
destToRadians, destIsLatLong);
if (!dst) {
emess(3, "cannot instantiate target coordinate system");
@@ -760,7 +763,7 @@ int main(int argc, char **argv) {
emess(3,
"missing target CRS and source CRS is not a projected CRS");
}
- destIsGeog = true;
+ destIsLongLat = true;
} else if (fromStr.empty()) {
assert(dst);
bool ignored;
@@ -770,7 +773,7 @@ int main(int argc, char **argv) {
emess(3,
"missing source CRS and target CRS is not a projected CRS");
}
- srcIsGeog = true;
+ srcIsLongLat = true;
}
proj_destroy(src);
proj_destroy(dst);
@@ -835,13 +838,13 @@ int main(int argc, char **argv) {
}
/* set input formatting control */
- if (!srcIsGeog)
+ if (!srcIsLongLat)
informat = strtod;
else {
informat = dmstor;
}
- if (!destIsGeog && !oform)
+ if (!destIsLongLat && !oform)
oform = "%.2f";
/* process input file list */
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], {});
}
// ---------------------------------------------------------------------------
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/cli/testvarious b/test/cli/testvarious
index 6e9cd43c..80340369 100755
--- a/test/cli/testvarious
+++ b/test/cli/testvarious
@@ -144,10 +144,10 @@ $EXE +proj=geocent +datum=WGS84 \
EOF
#
echo "#############################################################" >> ${OUT}
-echo Test conversion from geocentric latlong to geodetic latlong >> ${OUT}
+echo Test conversion from geodetic latlong to geocentric latlong >> ${OUT}
#
-$EXE +proj=latlong +datum=WGS84 +geoc \
- +to +proj=latlong +datum=WGS84 \
+$EXE +proj=latlong +datum=WGS84 \
+ +to +proj=latlong +datum=WGS84 +geoc \
-E >>${OUT} <<EOF
0d00'00.000"W 0d00'00.000"N 0.0
79d00'00.000"W 45d00'00.000"N 0.0
@@ -156,10 +156,10 @@ $EXE +proj=latlong +datum=WGS84 +geoc \
EOF
#
echo "#############################################################" >> ${OUT}
-echo Test conversion from geodetic latlong to geocentric latlong >> ${OUT}
+echo Test conversion from geocentric latlong to geodetic latlong >> ${OUT}
#
-$EXE +proj=latlong +datum=WGS84 \
- +to +proj=latlong +datum=WGS84 +geoc \
+$EXE +proj=latlong +datum=WGS84 +geoc \
+ +to +proj=latlong +datum=WGS84 \
-E >>${OUT} <<EOF
0d00'00.000"W 0d00'00.000"N 0.0
79d00'00.000"W 44d48'27.276"N 0.000
diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
index 9bb85933..06aa1a2d 100644
--- a/test/cli/tv_out.dist
+++ b/test/cli/tv_out.dist
@@ -43,13 +43,13 @@ Test geocentric x/y/z consumption.
861996.98 -4434590.01 4487348.41 79dW 45dN 0.001
0.00 -0.00 6356752.31 0dE 90dN -0.004
#############################################################
-Test conversion from geocentric latlong to geodetic latlong
+Test conversion from geodetic latlong to geocentric latlong
0d00'00.000"W 0d00'00.000"N 0.0 0dE 0dN 0.000
79d00'00.000"W 45d00'00.000"N 0.0 79dW 44d48'27.276"N 0.000
12d00'00.000"W 45d00'00.000"N 0.0 12dW 44d48'27.276"N 0.000
0d00'00.000"W 90d00'00.000"N 0.0 0dE 90dN 0.000
#############################################################
-Test conversion from geodetic latlong to geocentric latlong
+Test conversion from geocentric latlong to geodetic latlong
0d00'00.000"W 0d00'00.000"N 0.0 0dE 0dN 0.000
79d00'00.000"W 44d48'27.276"N 0.000 79dW 45dN 0.000
12d00'00.000"W 44d48'27.276"N 0.0 12dW 45dN 0.000
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 72ac1868..791a8b9e 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -894,6 +894,37 @@ TEST(wkt_parse, wkt2_EPSG_4979) {
// ---------------------------------------------------------------------------
+TEST(wkt_parse, wkt2_spherical_planetocentric) {
+ const auto wkt =
+ "GEODCRS[\"Mercury (2015) / Ocentric\",\n"
+ " DATUM[\"Mercury (2015)\",\n"
+ " ELLIPSOID[\"Mercury (2015)\",2440530,1075.12334801762,\n"
+ " LENGTHUNIT[\"metre\",1]],\n"
+ " ANCHOR[\"Hun Kal: 20.0\"]],\n"
+ " PRIMEM[\"Reference Meridian\",0,\n"
+ " ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+ " CS[spherical,2],\n"
+ " AXIS[\"planetocentric latitude (U)\",north,\n"
+ " ORDER[1],\n"
+ " ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+ " AXIS[\"planetocentric longitude (V)\",east,\n"
+ " ORDER[2],\n"
+ " ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+ " ID[\"IAU\",19902,2015],\n"
+ " REMARK[\"Source of IAU Coordinate systems: "
+ "doi://10.1007/s10569-017-9805-5\"]]";
+ auto obj = WKTParser().createFromWKT(wkt);
+ auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+ ASSERT_TRUE(crs != nullptr);
+ EXPECT_TRUE(crs->isSphericalPlanetocentric());
+ EXPECT_EQ(
+ crs->exportToWKT(
+ WKTFormatter::create(WKTFormatter::Convention::WKT2_2019).get()),
+ wkt);
+}
+
+// ---------------------------------------------------------------------------
+
static void checkGeocentric(GeodeticCRSPtr crs) {
// Explicitly check this is NOT a GeographicCRS
EXPECT_TRUE(!dynamic_cast<GeographicCRS *>(crs.get()));
@@ -10032,6 +10063,29 @@ TEST(io, projparse_geocent_wktext) {
// ---------------------------------------------------------------------------
+TEST(io, projparse_geoc) {
+ std::string input("+proj=longlat +geoc +datum=WGS84 +no_defs +type=crs");
+ auto obj = PROJStringParser().createFromPROJString(input);
+ auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+ ASSERT_TRUE(crs != nullptr);
+ EXPECT_TRUE(crs->isSphericalPlanetocentric());
+#if 1
+ EXPECT_THROW(
+ crs->exportToPROJString(
+ PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+ .get()),
+ FormattingException);
+#else
+ EXPECT_EQ(
+ crs->exportToPROJString(
+ PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+ .get()),
+ input);
+#endif
+}
+
+// ---------------------------------------------------------------------------
+
TEST(io, projparse_projected_wktext) {
std::string input("+proj=merc +foo +wktext +type=crs");
auto obj = PROJStringParser().createFromPROJString(input);
@@ -10178,13 +10232,13 @@ TEST(io, projparse_init) {
}
{
- auto obj = createFromUserInput("+title=mytitle +geoc +init=epsg:4326",
+ auto obj = createFromUserInput("+title=mytitle +over +init=epsg:4326",
dbContext, true);
auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
ASSERT_TRUE(crs != nullptr);
EXPECT_EQ(crs->nameStr(), "mytitle");
EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
- "+proj=longlat +geoc +datum=WGS84 +no_defs +type=crs");
+ "+proj=longlat +over +datum=WGS84 +no_defs +type=crs");
}
{
@@ -11476,6 +11530,59 @@ TEST(json_import, geographic_crs) {
// ---------------------------------------------------------------------------
+TEST(json_import, spherical_planetocentric) {
+ const auto json = "{\n"
+ " \"$schema\": \"foo\",\n"
+ " \"type\": \"GeodeticCRS\",\n"
+ " \"name\": \"Mercury (2015) / Ocentric\",\n"
+ " \"datum\": {\n"
+ " \"type\": \"GeodeticReferenceFrame\",\n"
+ " \"name\": \"Mercury (2015)\",\n"
+ " \"anchor\": \"Hun Kal: 20.0\",\n"
+ " \"ellipsoid\": {\n"
+ " \"name\": \"Mercury (2015)\",\n"
+ " \"semi_major_axis\": 2440530,\n"
+ " \"inverse_flattening\": 1075.12334801762\n"
+ " },\n"
+ " \"prime_meridian\": {\n"
+ " \"name\": \"Reference Meridian\",\n"
+ " \"longitude\": 0\n"
+ " }\n"
+ " },\n"
+ " \"coordinate_system\": {\n"
+ " \"subtype\": \"spherical\",\n"
+ " \"axis\": [\n"
+ " {\n"
+ " \"name\": \"Planetocentric latitude\",\n"
+ " \"abbreviation\": \"U\",\n"
+ " \"direction\": \"north\",\n"
+ " \"unit\": \"degree\"\n"
+ " },\n"
+ " {\n"
+ " \"name\": \"Planetocentric longitude\",\n"
+ " \"abbreviation\": \"V\",\n"
+ " \"direction\": \"east\",\n"
+ " \"unit\": \"degree\"\n"
+ " }\n"
+ " ]\n"
+ " },\n"
+ " \"id\": {\n"
+ " \"authority\": \"IAU\",\n"
+ " \"code\": 19902\n"
+ " },\n"
+ " \"remarks\": \"Source of IAU Coordinate systems: "
+ "doi://10.1007/s10569-017-9805-5\"\n"
+ "}";
+ auto obj = createFromUserInput(json, nullptr);
+ auto gcrs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+ ASSERT_TRUE(gcrs != nullptr);
+ EXPECT_TRUE(gcrs->isSphericalPlanetocentric());
+ EXPECT_EQ(gcrs->exportToJSON(&(JSONFormatter::create()->setSchema("foo"))),
+ json);
+}
+
+// ---------------------------------------------------------------------------
+
TEST(json_import, geographic_crs_errors) {
EXPECT_THROW(
createFromUserInput(
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
index 56d9f743..b50542c6 100644
--- a/test/unit/test_operationfactory.cpp
+++ b/test/unit/test_operationfactory.cpp
@@ -6092,7 +6092,7 @@ TEST(operation, createOperation_on_crs_with_canonical_bound_crs) {
TEST(operation, createOperation_fallback_to_proj4_strings) {
auto objDest = PROJStringParser().createFromPROJString(
- "+proj=longlat +geoc +datum=WGS84 +type=crs");
+ "+proj=longlat +over +datum=WGS84 +type=crs");
auto dest = nn_dynamic_pointer_cast<GeographicCRS>(objDest);
ASSERT_TRUE(dest != nullptr);
@@ -6102,7 +6102,7 @@ TEST(operation, createOperation_fallback_to_proj4_strings) {
EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline +step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
- "+step +proj=longlat +geoc +datum=WGS84 "
+ "+step +proj=longlat +over +datum=WGS84 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg");
}
@@ -6110,12 +6110,12 @@ TEST(operation, createOperation_fallback_to_proj4_strings) {
TEST(operation, createOperation_fallback_to_proj4_strings_bound_of_geog) {
auto objSrc = PROJStringParser().createFromPROJString(
- "+proj=longlat +geoc +ellps=GRS80 +towgs84=0,0,0 +type=crs");
+ "+proj=longlat +over +ellps=GRS80 +towgs84=0,0,0 +type=crs");
auto src = nn_dynamic_pointer_cast<BoundCRS>(objSrc);
ASSERT_TRUE(src != nullptr);
auto objDest = PROJStringParser().createFromPROJString(
- "+proj=longlat +geoc +ellps=clrk66 +towgs84=0,0,0 +type=crs");
+ "+proj=longlat +over +ellps=clrk66 +towgs84=0,0,0 +type=crs");
auto dest = nn_dynamic_pointer_cast<BoundCRS>(objDest);
ASSERT_TRUE(dest != nullptr);
@@ -6125,8 +6125,8 @@ TEST(operation, createOperation_fallback_to_proj4_strings_bound_of_geog) {
EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
- "+step +inv +proj=longlat +geoc +ellps=GRS80 +towgs84=0,0,0 "
- "+step +proj=longlat +geoc +ellps=clrk66 +towgs84=0,0,0 "
+ "+step +inv +proj=longlat +over +ellps=GRS80 +towgs84=0,0,0 "
+ "+step +proj=longlat +over +ellps=clrk66 +towgs84=0,0,0 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg");
}
@@ -6832,3 +6832,140 @@ TEST(operation, derivedGeographicCRS_with_to_wgs84_to_geographicCRS) {
"+proj=pipeline +step +proj=axisswap +order=2,1 " + pipeline);
}
}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, createOperation_spherical_ocentric_to_geographic) {
+ auto objSrc = PROJStringParser().createFromPROJString(
+ "+proj=longlat +geoc +datum=WGS84 +type=crs");
+ auto src = nn_dynamic_pointer_cast<GeodeticCRS>(objSrc);
+ ASSERT_TRUE(src != nullptr);
+
+ auto op = CoordinateOperationFactory::create()->createOperation(
+ NN_CHECK_ASSERT(src), GeographicCRS::EPSG_4326);
+ 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=unitconvert +xy_in=rad +xy_out=deg "
+ "+step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, createOperation_geographic_to_spherical_ocentric) {
+ auto objDest = PROJStringParser().createFromPROJString(
+ "+proj=longlat +geoc +datum=WGS84 +type=crs");
+ auto dest = nn_dynamic_pointer_cast<GeodeticCRS>(objDest);
+ ASSERT_TRUE(dest != nullptr);
+
+ auto op = CoordinateOperationFactory::create()->createOperation(
+ GeographicCRS::EPSG_4326, NN_CHECK_ASSERT(dest));
+ ASSERT_TRUE(op != nullptr);
+ EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+ "+proj=pipeline "
+ "+step +proj=axisswap +order=2,1 "
+ "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+ "+step +proj=geoc +ellps=WGS84 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+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");
+ auto src = nn_dynamic_pointer_cast<BoundCRS>(objSrc);
+ ASSERT_TRUE(src != nullptr);
+
+ auto objDest = PROJStringParser().createFromPROJString(
+ "+proj=longlat +geoc +ellps=clrk66 +towgs84=4,5,6 +type=crs");
+ auto dest = nn_dynamic_pointer_cast<BoundCRS>(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=GRS80 "
+ "+step +proj=push +v_3 "
+ "+step +proj=cart +ellps=GRS80 "
+ "+step +proj=helmert +x=-3 +y=-3 +z=-3 "
+ "+step +inv +proj=cart +ellps=clrk66 "
+ "+step +proj=pop +v_3 "
+ "+step +proj=geoc +ellps=clrk66 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, createOperation_spherical_ocentric_to_projected) {
+ auto objSrc = PROJStringParser().createFromPROJString(
+ "+proj=longlat +geoc +datum=WGS84 +type=crs");
+ auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+ ASSERT_TRUE(src != nullptr);
+
+ auto objDest = PROJStringParser().createFromPROJString(
+ "+proj=utm +zone=11 +datum=WGS84 +type=crs");
+ auto dest = nn_dynamic_pointer_cast<CRS>(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=utm +zone=11 +ellps=WGS84");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation,
+ createOperation_spherical_ocentric_to_projected_of_spherical_ocentric) {
+ auto objSrc = PROJStringParser().createFromPROJString(
+ "+proj=longlat +geoc +datum=WGS84 +type=crs");
+ auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+ ASSERT_TRUE(src != nullptr);
+
+ auto objDest = PROJStringParser().createFromPROJString(
+ "+proj=utm +geoc +zone=11 +datum=WGS84 +type=crs");
+ auto dest = nn_dynamic_pointer_cast<CRS>(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=utm +zone=11 +ellps=WGS84");
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 14:00:02 +0000