aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--include/proj/coordinateoperation.hpp11
-rw-r--r--include/proj/crs.hpp2
-rw-r--r--include/proj/internal/internal.hpp2
-rw-r--r--include/proj/io.hpp1
-rw-r--r--src/c_api.cpp147
-rw-r--r--src/coordinateoperation.cpp264
-rw-r--r--src/crs.cpp8
-rw-r--r--src/cs2cs.cpp327
-rw-r--r--src/io.cpp441
-rw-r--r--src/pj_apply_gridshift.c2
-rw-r--r--src/pj_fwd.c9
-rw-r--r--src/proj.h13
-rw-r--r--src/proj_4D_api.c7
-rw-r--r--test/gie/more_builtins.gie10
-rw-r--r--test/old/proj_outIGNF.dist42
-rwxr-xr-xtest/old/testIGNF12
-rwxr-xr-xtest/old/testdatumfile4
-rwxr-xr-xtest/old/testvarious24
-rw-r--r--test/old/tv_out.dist12
-rw-r--r--test/unit/test_c_api.cpp81
-rw-r--r--test/unit/test_io.cpp103
-rw-r--r--test/unit/test_operation.cpp128
22 files changed, 1278 insertions, 372 deletions
diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp
index 293ba480..4005176a 100644
--- a/include/proj/coordinateoperation.hpp
+++ b/include/proj/coordinateoperation.hpp
@@ -818,10 +818,6 @@ class PROJ_GCC_DLL Conversion : public SingleOperation {
PROJ_DLL ~Conversion() override;
//! @endcond
- //! @cond Doxygen_Suppress
- PROJ_INTERNAL ConversionNNPtr shallowClone() const;
- //! @endcond
-
PROJ_DLL CoordinateOperationNNPtr inverse() const override;
//! @cond Doxygen_Suppress
@@ -1287,6 +1283,9 @@ class PROJ_GCC_DLL Conversion : public SingleOperation {
PROJ_INTERNAL const char *getESRIMethodName() const;
PROJ_INTERNAL const char *getWKT1GDALMethodName() const;
+
+ PROJ_INTERNAL ConversionNNPtr shallowClone() const;
+
//! @endcond
protected:
@@ -1490,6 +1489,8 @@ class PROJ_GCC_DLL Transformation : public SingleOperation {
PROJ_INTERNAL void _exportToWKT(io::WKTFormatter *formatter)
const override; // throw(io::FormattingException)
+ PROJ_INTERNAL TransformationNNPtr shallowClone() const;
+
//! @endcond
protected:
@@ -1499,6 +1500,7 @@ class PROJ_GCC_DLL Transformation : public SingleOperation {
const OperationMethodNNPtr &methodIn,
const std::vector<GeneralParameterValueNNPtr> &values,
const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies);
+ PROJ_INTERNAL Transformation(const Transformation &other);
INLINED_MAKE_SHARED
PROJ_INTERNAL void _exportToPROJString(io::PROJStringFormatter *formatter)
@@ -1508,7 +1510,6 @@ class PROJ_GCC_DLL Transformation : public SingleOperation {
private:
PROJ_OPAQUE_PRIVATE_DATA
- Transformation(const Transformation &) = delete;
};
// ---------------------------------------------------------------------------
diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp
index 82b2bd49..10e0a639 100644
--- a/include/proj/crs.hpp
+++ b/include/proj/crs.hpp
@@ -108,6 +108,8 @@ class PROJ_GCC_DLL CRS : public common::ObjectUsage {
PROJ_FOR_TEST CRSNNPtr shallowClone() const;
+ PROJ_INTERNAL const std::string &getExtensionProj4() const noexcept;
+
//! @endcond
protected:
diff --git a/include/proj/internal/internal.hpp b/include/proj/internal/internal.hpp
index 85dd5ac3..b16e12bc 100644
--- a/include/proj/internal/internal.hpp
+++ b/include/proj/internal/internal.hpp
@@ -104,7 +104,7 @@ PROJ_FOR_TEST std::string replaceAll(const std::string &str,
const std::string &before,
const std::string &after);
-size_t ci_find(const std::string &osStr, const char *needle) noexcept;
+PROJ_DLL size_t ci_find(const std::string &osStr, const char *needle) noexcept;
size_t ci_find(const std::string &osStr, const std::string &needle,
size_t startPos = 0) noexcept;
diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index 3b5019c1..653edfbc 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -365,6 +365,7 @@ class PROJ_GCC_DLL PROJStringFormatter {
PROJ_DLL void stopInversion();
PROJ_INTERNAL bool isInverted() const;
PROJ_INTERNAL bool getUseETMercForTMerc(bool &settingSetOut) const;
+ PROJ_INTERNAL void setCoordinateOperationOptimizations(bool enable);
PROJ_DLL void
ingestPROJString(const std::string &str); // throw ParsingException
diff --git a/src/c_api.cpp b/src/c_api.cpp
index ba1b9534..3b51d905 100644
--- a/src/c_api.cpp
+++ b/src/c_api.cpp
@@ -39,6 +39,7 @@
#include "proj/common.hpp"
#include "proj/coordinateoperation.hpp"
+#include "proj/coordinatesystem.hpp"
#include "proj/crs.hpp"
#include "proj/datum.hpp"
#include "proj/io.hpp"
@@ -56,6 +57,7 @@
using namespace NS_PROJ::common;
using namespace NS_PROJ::crs;
+using namespace NS_PROJ::cs;
using namespace NS_PROJ::datum;
using namespace NS_PROJ::io;
using namespace NS_PROJ::internal;
@@ -1354,8 +1356,8 @@ int proj_obj_prime_meridian_get_parameters(PJ_OBJ *prime_meridian,
// ---------------------------------------------------------------------------
-/** \brief Return the base CRS of a BoundCRS or the source CRS of a
- * CoordinateOperation.
+/** \brief Return the base CRS of a BoundCRS or a DerivedCRS/ProjectedCRS, or
+ * the source CRS of a CoordinateOperation.
*
* The returned object must be unreferenced with proj_obj_unref() after
* use.
@@ -1372,6 +1374,10 @@ PJ_OBJ *proj_obj_get_source_crs(PJ_OBJ *obj) {
if (boundCRS) {
return PJ_OBJ::create(obj->ctx, boundCRS->baseCRS());
}
+ auto derivedCRS = dynamic_cast<const DerivedCRS *>(ptr);
+ if (derivedCRS) {
+ return PJ_OBJ::create(obj->ctx, derivedCRS->baseCRS());
+ }
auto co = dynamic_cast<const CoordinateOperation *>(ptr);
if (co) {
auto sourceCRS = co->sourceCRS();
@@ -4071,3 +4077,140 @@ double proj_coordoperation_get_accuracy(PJ_OBJ *coordoperation) {
}
return -1;
}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns the coordinate system of a SingleCRS.
+ *
+ * The returned object must be unreferenced with proj_obj_unref() after
+ * use.
+ * It should be used by at most one thread at a time.
+ *
+ * @param crs Objet of type SingleCRS (must not be NULL)
+ * @return Object that must be unreferenced with proj_obj_unref(), or NULL
+ * in case of error.
+ */
+PJ_OBJ *proj_obj_crs_get_coordinate_system(PJ_OBJ *crs) {
+ assert(crs);
+ auto l_crs = dynamic_cast<const SingleCRS *>(crs->obj.get());
+ if (!l_crs) {
+ proj_log_error(crs->ctx, __FUNCTION__, "Object is not a SingleCRS");
+ return nullptr;
+ }
+ return PJ_OBJ::create(crs->ctx, l_crs->coordinateSystem());
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns the type of the coordinate system.
+ *
+ * @param cs Objet of type CoordinateSystem (must not be NULL)
+ * @return type, or NULL in case of error.
+ */
+const char *proj_obj_cs_get_type(PJ_OBJ *cs) {
+ assert(cs);
+ auto l_cs = dynamic_cast<const CoordinateSystem *>(cs->obj.get());
+ if (!l_cs) {
+ proj_log_error(cs->ctx, __FUNCTION__,
+ "Object is not a CoordinateSystem");
+ return nullptr;
+ }
+ if (dynamic_cast<const CartesianCS *>(l_cs)) {
+ return "Cartesian";
+ }
+ if (dynamic_cast<const EllipsoidalCS *>(l_cs)) {
+ return "Ellipsoidal";
+ }
+ if (dynamic_cast<const VerticalCS *>(l_cs)) {
+ return "Vertical";
+ }
+ if (dynamic_cast<const SphericalCS *>(l_cs)) {
+ return "Spherical";
+ }
+ if (dynamic_cast<const OrdinalCS *>(l_cs)) {
+ return "Ordinal";
+ }
+ if (dynamic_cast<const ParametricCS *>(l_cs)) {
+ return "Parametric";
+ }
+ if (dynamic_cast<const DateTimeTemporalCS *>(l_cs)) {
+ return "DateTimeTemporal";
+ }
+ if (dynamic_cast<const TemporalCountCS *>(l_cs)) {
+ return "TemporalCount";
+ }
+ if (dynamic_cast<const TemporalMeasureCS *>(l_cs)) {
+ return "TemporalMeasure";
+ }
+ return "unknown";
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns the number of axis of the coordinate system.
+ *
+ * @param cs Objet of type CoordinateSystem (must not be NULL)
+ * @return number of axis, or -1 in case of error.
+ */
+int proj_obj_cs_get_axis_count(PJ_OBJ *cs) {
+ assert(cs);
+ auto l_cs = dynamic_cast<const CoordinateSystem *>(cs->obj.get());
+ if (!l_cs) {
+ proj_log_error(cs->ctx, __FUNCTION__,
+ "Object is not a CoordinateSystem");
+ return -1;
+ }
+ return static_cast<int>(l_cs->axisList().size());
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns information on an axis
+ *
+ * @param cs Objet of type CoordinateSystem (must not be NULL)
+ * @param index Index of the coordinate system (between 0 and
+ * proj_obj_cs_get_axis_count() - 1)
+ * @param pName Pointer to a string value to store the axis name. or NULL
+ * @param pAbbrev Pointer to a string value to store the axis abbreviation. or
+ * NULL
+ * @param pDirection Pointer to a string value to store the axis direction. or
+ * NULL
+ * @param pUnitConvFactor Pointer to a double value to store the axis
+ * unit conversion factor. or NULL
+ * @param pUnitName Pointer to a string value to store the axis
+ * unit name. or NULL
+ * @return TRUE in case of success
+ */
+int proj_obj_cs_get_axis_info(PJ_OBJ *cs, int index, const char **pName,
+ const char **pAbbrev, const char **pDirection,
+ double *pUnitConvFactor, const char **pUnitName) {
+ assert(cs);
+ auto l_cs = dynamic_cast<const CoordinateSystem *>(cs->obj.get());
+ if (!l_cs) {
+ proj_log_error(cs->ctx, __FUNCTION__,
+ "Object is not a CoordinateSystem");
+ return false;
+ }
+ const auto &axisList = l_cs->axisList();
+ if (index < 0 || static_cast<size_t>(index) >= axisList.size()) {
+ proj_log_error(cs->ctx, __FUNCTION__, "Invalid index");
+ return false;
+ }
+ const auto &axis = axisList[index];
+ if (pName) {
+ *pName = axis->nameStr().c_str();
+ }
+ if (pAbbrev) {
+ *pAbbrev = axis->abbreviation().c_str();
+ }
+ if (pDirection) {
+ *pDirection = axis->direction().toString().c_str();
+ }
+ if (pUnitConvFactor) {
+ *pUnitConvFactor = axis->unit().conversionToSI();
+ }
+ if (pUnitName) {
+ *pUnitName = axis->unit().name().c_str();
+ }
+ return true;
+}
diff --git a/src/coordinateoperation.cpp b/src/coordinateoperation.cpp
index e0e02931..893b52d3 100644
--- a/src/coordinateoperation.cpp
+++ b/src/coordinateoperation.cpp
@@ -5471,6 +5471,12 @@ Transformation::~Transformation() = default;
// ---------------------------------------------------------------------------
+Transformation::Transformation(const Transformation &other)
+ : CoordinateOperation(other), SingleOperation(other),
+ d(internal::make_unique<Private>(*other.d)) {}
+
+// ---------------------------------------------------------------------------
+
/** \brief Return the source crs::CRS of the transformation.
*
* @return the source CRS.
@@ -5492,6 +5498,17 @@ const crs::CRSNNPtr &Transformation::targetCRS() PROJ_CONST_DEFN {
// ---------------------------------------------------------------------------
//! @cond Doxygen_Suppress
+TransformationNNPtr Transformation::shallowClone() const {
+ auto conv = Transformation::nn_make_shared<Transformation>(*this);
+ conv->assignSelf(conv);
+ conv->setCRSs(this, false);
+ return conv;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
/** \brief Return the TOWGS84 parameters of the transformation.
*
* If this transformation uses Coordinate Frame Rotation, Position Vector
@@ -7595,6 +7612,8 @@ void Transformation::_exportToPROJString(
"Transformation cannot be exported as a PROJ.4 string");
}
+ formatter->setCoordinateOperationOptimizations(true);
+
bool positionVectorConvention = true;
bool sevenParamsTransform = false;
bool threeParamsTransform = false;
@@ -10156,6 +10175,72 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
const auto candidatesDstGeod(
findCandidateGeodCRSForDatum(authFactory, geodDst->datum()));
+ auto createTransformations = [&](const crs::CRSNNPtr &candidateSrcGeod,
+ const crs::CRSNNPtr &candidateDstGeod,
+ const CoordinateOperationNNPtr &opFirst,
+ bool isNullFirst) {
+ const auto opsSecond =
+ createOperations(candidateSrcGeod, candidateDstGeod, context);
+ const auto opsThird =
+ createOperations(candidateDstGeod, targetCRS, context);
+ assert(!opsThird.empty());
+
+ for (auto &opSecond : opsSecond) {
+ // Check that it is not a transformation synthetized by
+ // ourselves
+ if (!hasIdentifiers(opSecond)) {
+ continue;
+ }
+ // And even if it is a referenced transformation, check that
+ // it is not a trivial one
+ auto so = dynamic_cast<const SingleOperation *>(opSecond.get());
+ if (so && isAxisOrderReversal(so->method()->getEPSGCode())) {
+ continue;
+ }
+
+ std::vector<CoordinateOperationNNPtr> subOps;
+ if (isNullFirst) {
+ opSecond->setCRSs(
+ sourceCRS, NN_CHECK_ASSERT(opSecond->targetCRS()), nullptr);
+ } else {
+ subOps.emplace_back(opFirst);
+ }
+ if (isNullTransformation(opsThird[0]->nameStr())) {
+ opSecond->setCRSs(NN_CHECK_ASSERT(opSecond->sourceCRS()),
+ targetCRS, nullptr);
+ subOps.emplace_back(opSecond);
+ } else {
+ subOps.emplace_back(opSecond);
+ subOps.emplace_back(opsThird[0]);
+ }
+ res.emplace_back(ConcatenatedOperation::createComputeMetadata(
+ subOps, !allowEmptyIntersection));
+ }
+ };
+
+ // Start in priority with candidates that have exactly the same name as
+ // the sourcCRS and targetCRS. Typically for the case of init=IGNF:XXXX
+ for (const auto &candidateSrcGeod : candidatesSrcGeod) {
+ if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
+ for (const auto &candidateDstGeod : candidatesDstGeod) {
+ if (candidateDstGeod->nameStr() == targetCRS->nameStr()) {
+ const auto opsFirst =
+ createOperations(sourceCRS, candidateSrcGeod, context);
+ assert(!opsFirst.empty());
+ const bool isNullFirst =
+ isNullTransformation(opsFirst[0]->nameStr());
+ createTransformations(candidateSrcGeod, candidateDstGeod,
+ opsFirst[0], isNullFirst);
+ if (!res.empty()) {
+ return;
+ }
+ break;
+ }
+ }
+ break;
+ }
+ }
+
for (const auto &candidateSrcGeod : candidatesSrcGeod) {
const auto opsFirst =
createOperations(sourceCRS, candidateSrcGeod, context);
@@ -10163,44 +10248,8 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
const bool isNullFirst = isNullTransformation(opsFirst[0]->nameStr());
for (const auto &candidateDstGeod : candidatesDstGeod) {
- const auto opsSecond =
- createOperations(candidateSrcGeod, candidateDstGeod, context);
- const auto opsThird =
- createOperations(candidateDstGeod, targetCRS, context);
- assert(!opsThird.empty());
-
- for (auto &opSecond : opsSecond) {
- // Check that it is not a transformation synthetized by
- // ourselves
- if (!hasIdentifiers(opSecond)) {
- continue;
- }
- // And even if it is a referenced transformation, check that
- // it is not a trivial one
- auto so = dynamic_cast<const SingleOperation *>(opSecond.get());
- if (so && isAxisOrderReversal(so->method()->getEPSGCode())) {
- continue;
- }
-
- std::vector<CoordinateOperationNNPtr> subOps;
- if (isNullFirst) {
- opSecond->setCRSs(sourceCRS,
- NN_CHECK_ASSERT(opSecond->targetCRS()),
- nullptr);
- } else {
- subOps.emplace_back(opsFirst[0]);
- }
- if (isNullTransformation(opsThird[0]->nameStr())) {
- opSecond->setCRSs(NN_CHECK_ASSERT(opSecond->sourceCRS()),
- targetCRS, nullptr);
- subOps.emplace_back(opSecond);
- } else {
- subOps.emplace_back(opSecond);
- subOps.emplace_back(opsThird[0]);
- }
- res.emplace_back(ConcatenatedOperation::createComputeMetadata(
- subOps, !allowEmptyIntersection));
- }
+ createTransformations(candidateSrcGeod, candidateDstGeod,
+ opsFirst[0], isNullFirst);
}
if (!res.empty()) {
return;
@@ -10261,6 +10310,56 @@ CoordinateOperationFactory::Private::createOperations(
std::vector<CoordinateOperationNNPtr> res;
const bool allowEmptyIntersection = true;
+ const auto &sourceProj4Ext = sourceCRS->getExtensionProj4();
+ const auto &targetProj4Ext = targetCRS->getExtensionProj4();
+ if (!sourceProj4Ext.empty() || !targetProj4Ext.empty()) {
+
+ auto sourceProjExportable =
+ dynamic_cast<io::IPROJStringExportable *>(sourceCRS.get());
+ auto targetProjExportable =
+ dynamic_cast<io::IPROJStringExportable *>(targetCRS.get());
+ if (!sourceProjExportable) {
+ throw InvalidOperation("Source CRS is not PROJ exportable");
+ }
+ if (!targetProjExportable) {
+ throw InvalidOperation("Target CRS is not PROJ exportable");
+ }
+ auto projFormatter = io::PROJStringFormatter::create(
+ io::PROJStringFormatter::Convention::PROJ_4);
+ projFormatter->startInversion();
+ sourceProjExportable->_exportToPROJString(projFormatter.get());
+
+ auto geogSrc =
+ dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
+ if (geogSrc) {
+ auto proj5Formatter = io::PROJStringFormatter::create(
+ io::PROJStringFormatter::Convention::PROJ_5);
+ geogSrc->addAngularUnitConvertAndAxisSwap(proj5Formatter.get());
+ projFormatter->ingestPROJString(proj5Formatter->toString());
+ }
+
+ projFormatter->stopInversion();
+
+ targetProjExportable->_exportToPROJString(projFormatter.get());
+
+ auto geogDst =
+ dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
+ if (geogDst) {
+ auto proj5Formatter = io::PROJStringFormatter::create(
+ io::PROJStringFormatter::Convention::PROJ_5);
+ geogDst->addAngularUnitConvertAndAxisSwap(proj5Formatter.get());
+ projFormatter->ingestPROJString(proj5Formatter->toString());
+ }
+
+ const auto PROJString = projFormatter->toString();
+ auto properties = util::PropertyMap().set(
+ common::IdentifiedObject::NAME_KEY,
+ buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()));
+ res.emplace_back(SingleOperation::createPROJBased(
+ properties, PROJString, sourceCRS, targetCRS, {}));
+ return res;
+ }
+
auto geodSrc = dynamic_cast<const crs::GeodeticCRS *>(sourceCRS.get());
auto geodDst = dynamic_cast<const crs::GeodeticCRS *>(targetCRS.get());
@@ -10317,7 +10416,9 @@ CoordinateOperationFactory::Private::createOperations(
const auto &dstDatum = geodDst->datum();
const bool dstHasDatumWithId =
dstDatum && !dstDatum->identifiers().empty();
- if (srcHasDatumWithId && dstHasDatumWithId) {
+ if (srcHasDatumWithId && dstHasDatumWithId &&
+ !srcDatum->_isEquivalentTo(
+ dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) {
createOperationsWithDatumPivot(res, sourceCRS, targetCRS,
geodSrc, geodDst, context);
doFilterAndCheckPerfectOp = !res.empty();
@@ -10443,11 +10544,10 @@ CoordinateOperationFactory::Private::createOperations(
dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
auto geogCRSOfBaseOfBoundSrc =
boundSrc->baseCRS()->extractGeographicCRS();
- if (hubSrcGeog &&
+ if (hubSrcGeog && geogCRSOfBaseOfBoundSrc &&
(hubSrcGeog->_isEquivalentTo(
geogDst, util::IComparable::Criterion::EQUIVALENT) ||
- hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst))) &&
- geogCRSOfBaseOfBoundSrc) {
+ hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst)))) {
if (boundSrc->baseCRS() == geogCRSOfBaseOfBoundSrc) {
// Optimization to avoid creating a useless concatenated
// operation
@@ -10471,6 +10571,86 @@ CoordinateOperationFactory::Private::createOperations(
return res;
}
}
+ // If the datum are equivalent, this is also fine
+ } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
+ geogDst->datum() &&
+ hubSrcGeog->datum()->_isEquivalentTo(
+ geogDst->datum().get(),
+ util::IComparable::Criterion::EQUIVALENT)) {
+ auto opsFirst =
+ createOperations(boundSrc->baseCRS(),
+ NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
+ auto opsLast = createOperations(hubSrc, targetCRS, context);
+ if (!opsFirst.empty() && !opsLast.empty()) {
+ for (const auto &opFirst : opsFirst) {
+ for (const auto &opLast : opsLast) {
+ try {
+ res.emplace_back(
+ ConcatenatedOperation::createComputeMetadata(
+ {opFirst, boundSrc->transformation(),
+ opLast},
+ !allowEmptyIntersection));
+ } catch (const InvalidOperationEmptyIntersection &) {
+ }
+ }
+ }
+ if (!res.empty()) {
+ return res;
+ }
+ }
+ // Consider WGS 84 and NAD83 as equivalent in that context if the
+ // geogCRSOfBaseOfBoundSrc ellipsoid is Clarke66 (for NAD27)
+ // Case of "+proj=latlong +ellps=clrk66
+ // +nadgrids=ntv1_can.dat,conus"
+ // to "+proj=latlong +datum=NAD83"
+ } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
+ geogDst->datum() &&
+ geogCRSOfBaseOfBoundSrc->ellipsoid()->_isEquivalentTo(
+ datum::Ellipsoid::CLARKE_1866.get(),
+ util::IComparable::Criterion::EQUIVALENT) &&
+ hubSrcGeog->datum()->_isEquivalentTo(
+ datum::GeodeticReferenceFrame::EPSG_6326.get(),
+ util::IComparable::Criterion::EQUIVALENT) &&
+ geogDst->datum()->_isEquivalentTo(
+ datum::GeodeticReferenceFrame::EPSG_6269.get(),
+ util::IComparable::Criterion::EQUIVALENT)) {
+ auto nnGeogCRSOfBaseOfBoundSrc =
+ NN_NO_CHECK(geogCRSOfBaseOfBoundSrc);
+ if (boundSrc->baseCRS()->_isEquivalentTo(
+ nnGeogCRSOfBaseOfBoundSrc.get(),
+ util::IComparable::Criterion::EQUIVALENT)) {
+ auto transf = boundSrc->transformation()->shallowClone();
+ transf->setProperties(util::PropertyMap().set(
+ common::IdentifiedObject::NAME_KEY,
+ buildTransfName(boundSrc->baseCRS()->nameStr(),
+ targetCRS->nameStr())));
+ transf->setCRSs(boundSrc->baseCRS(), targetCRS, nullptr);
+ res.emplace_back(transf);
+ return res;
+ } else {
+ auto opsFirst = createOperations(
+ boundSrc->baseCRS(), nnGeogCRSOfBaseOfBoundSrc, context);
+ auto transf = boundSrc->transformation()->shallowClone();
+ transf->setProperties(util::PropertyMap().set(
+ common::IdentifiedObject::NAME_KEY,
+ buildTransfName(nnGeogCRSOfBaseOfBoundSrc->nameStr(),
+ targetCRS->nameStr())));
+ transf->setCRSs(nnGeogCRSOfBaseOfBoundSrc, targetCRS, nullptr);
+ if (!opsFirst.empty()) {
+ for (const auto &opFirst : opsFirst) {
+ try {
+ res.emplace_back(
+ ConcatenatedOperation::createComputeMetadata(
+ {opFirst, transf},
+ !allowEmptyIntersection));
+ } catch (const InvalidOperationEmptyIntersection &) {
+ }
+ }
+ if (!res.empty()) {
+ return res;
+ }
+ }
+ }
}
if (hubSrcGeog &&
diff --git a/src/crs.cpp b/src/crs.cpp
index a204a037..eec7a926 100644
--- a/src/crs.cpp
+++ b/src/crs.cpp
@@ -170,6 +170,14 @@ const GeodeticCRS *CRS::extractGeodeticCRSRaw() const {
// ---------------------------------------------------------------------------
+//! @cond Doxygen_Suppress
+const std::string &CRS::getExtensionProj4() const noexcept {
+ return d->extensionProj4_;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
/** \brief Return the GeographicCRS of the CRS.
*
* Returns the GeographicCRS contained in a CRS. This works currently with
diff --git a/src/cs2cs.cpp b/src/cs2cs.cpp
index 6f302ae3..6f4c4a55 100644
--- a/src/cs2cs.cpp
+++ b/src/cs2cs.cpp
@@ -26,6 +26,8 @@
* DEALINGS IN THE SOFTWARE.
*****************************************************************************/
+#define FROM_PROJ_CPP
+
#include <ctype.h>
#include <locale.h>
#include <math.h>
@@ -33,6 +35,11 @@
#include <stdlib.h>
#include <string.h>
+#include <cassert>
+#include <string>
+
+#include <proj/internal/internal.hpp>
+
// PROJ include order is sensitive
// clang-format off
#include "proj.h"
@@ -41,9 +48,15 @@
// clang-format on
#define MAX_LINE 1000
-#define MAX_PARGS 100
-static projPJ fromProj, toProj;
+static PJ *transformation = nullptr;
+
+static bool srcIsGeog = false;
+static double srcToRadians = 0.0;
+
+static bool destIsGeog = false;
+static double destToRadians = 0.0;
+static bool destIsLatLong = false;
static int reversein = 0, /* != 0 reverse input arguments */
reverseout = 0, /* != 0 reverse output arguments */
@@ -116,19 +129,45 @@ static void process(FILE *fid)
}
if (data.u != HUGE_VAL) {
- if (pj_transform(fromProj, toProj, 1, 0, &(data.u), &(data.v),
- &z) != 0) {
- data.u = HUGE_VAL;
- data.v = HUGE_VAL;
- emess(-3, "pj_transform(): %s", pj_strerrno(pj_errno));
+
+ if (srcIsGeog) {
+ /* dmstor gives values to radians. Convert now to the SRS unit
+ */
+ data.u /= srcToRadians;
+ data.v /= srcToRadians;
}
+
+ PJ_COORD coord;
+ coord.xyzt.x = data.u;
+ coord.xyzt.y = data.v;
+ coord.xyzt.z = z;
+ coord.xyzt.t = HUGE_VAL;
+ coord = proj_trans(transformation, PJ_FWD, coord);
+ data.u = coord.xyz.x;
+ data.v = coord.xyz.y;
+ z = coord.xyz.z;
}
if (data.u == HUGE_VAL) /* error output */
fputs(oterr, stdout);
- else if (pj_is_latlong(toProj) && !oform) { /*ascii DMS output */
- if (reverseout) {
+ else if (destIsGeog && !oform) { /*ascii DMS output */
+
+ // rtodms() expect radians: convert from the output SRS unit
+ data.u *= destToRadians;
+ data.v *= destToRadians;
+
+ if (destIsLatLong) {
+ if (reverseout) {
+ fputs(rtodms(pline, data.v, 'E', 'W'), stdout);
+ putchar('\t');
+ fputs(rtodms(pline, data.u, 'N', 'S'), stdout);
+ } else {
+ fputs(rtodms(pline, data.u, 'N', 'S'), stdout);
+ putchar('\t');
+ fputs(rtodms(pline, data.v, 'E', 'W'), stdout);
+ }
+ } else if (reverseout) {
fputs(rtodms(pline, data.v, 'N', 'S'), stdout);
putchar('\t');
fputs(rtodms(pline, data.u, 'E', 'W'), stdout);
@@ -139,9 +178,9 @@ static void process(FILE *fid)
}
} else { /* x-y or decimal degree ascii output */
- if (proj_angular_output(toProj, PJ_FWD)) {
- data.v *= RAD_TO_DEG;
- data.u *= RAD_TO_DEG;
+ if (destIsGeog) {
+ data.v *= destToRadians * RAD_TO_DEG;
+ data.u *= destToRadians * RAD_TO_DEG;
}
if (reverseout) {
printf(oform, data.v);
@@ -167,17 +206,116 @@ static void process(FILE *fid)
}
/************************************************************************/
+/* instanciate_crs() */
+/************************************************************************/
+
+static PJ_OBJ *instanciate_crs(const std::string &definition,
+ const char *const *optionsImportCRS,
+ bool &isGeog, double &toRadians,
+ bool &isLatFirst) {
+ PJ_OBJ *crs = proj_obj_create_from_user_input(nullptr, definition.c_str(),
+ optionsImportCRS);
+ if (!crs) {
+ return nullptr;
+ }
+
+ isGeog = false;
+ toRadians = 0.0;
+ isLatFirst = false;
+
+ auto type = proj_obj_get_type(crs);
+ if (type == PJ_OBJ_TYPE_BOUND_CRS) {
+ auto base = proj_obj_get_source_crs(crs);
+ proj_obj_unref(crs);
+ crs = base;
+ type = proj_obj_get_type(crs);
+ }
+ if (type == PJ_OBJ_TYPE_GEOGRAPHIC_2D_CRS ||
+ type == PJ_OBJ_TYPE_GEOGRAPHIC_3D_CRS) {
+ auto cs = proj_obj_crs_get_coordinate_system(crs);
+ assert(cs);
+
+ isGeog = true;
+ const char *axisName = "";
+ proj_obj_cs_get_axis_info(cs, 0,
+ &axisName, // name,
+ nullptr, // abbrev
+ nullptr, // direction
+ &toRadians,
+ nullptr // unit name
+ );
+ isLatFirst =
+ NS_PROJ::internal::ci_find(std::string(axisName), "latitude") !=
+ std::string::npos;
+
+ proj_obj_unref(cs);
+ }
+
+ return crs;
+}
+
+/************************************************************************/
+/* get_geog_crs_proj_string_from_proj_crs() */
+/************************************************************************/
+
+static std::string get_geog_crs_proj_string_from_proj_crs(PJ_OBJ *src,
+ double &toRadians,
+ bool &isLatFirst) {
+ auto srcType = proj_obj_get_type(src);
+ if (srcType == PJ_OBJ_TYPE_BOUND_CRS) {
+ auto base = proj_obj_get_source_crs(src);
+ assert(base);
+ proj_obj_unref(src);
+ src = base;
+ srcType = proj_obj_get_type(src);
+ }
+ if (srcType != PJ_OBJ_TYPE_PROJECTED_CRS) {
+ return std::string();
+ }
+
+ auto base = proj_obj_get_source_crs(src);
+ assert(base);
+ auto baseType = proj_obj_get_type(base);
+ if (baseType != PJ_OBJ_TYPE_GEOGRAPHIC_2D_CRS &&
+ baseType != PJ_OBJ_TYPE_GEOGRAPHIC_3D_CRS) {
+ proj_obj_unref(base);
+ return std::string();
+ }
+
+ auto cs = proj_obj_crs_get_coordinate_system(base);
+ assert(cs);
+
+ const char *axisName = "";
+ proj_obj_cs_get_axis_info(cs, 0,
+ &axisName, // name,
+ nullptr, // abbrev
+ nullptr, // direction
+ &toRadians,
+ nullptr // unit name
+ );
+ isLatFirst = NS_PROJ::internal::ci_find(std::string(axisName),
+ "latitude") != std::string::npos;
+
+ proj_obj_unref(cs);
+
+ auto retCStr = proj_obj_as_proj_string(base, PJ_PROJ_5, nullptr);
+ std::string ret(retCStr ? retCStr : "");
+ proj_obj_unref(base);
+ return ret;
+}
+
+/************************************************************************/
/* main() */
/************************************************************************/
int main(int argc, char **argv) {
char *arg;
char **eargv = argv;
- char *from_argv[MAX_PARGS];
- char *to_argv[MAX_PARGS];
+ std::string fromStr;
+ std::string toStr;
FILE *fid;
- int from_argc = 0, to_argc = 0, eargc = 0, mon = 0;
- int have_to_flag = 0, inverse = 0, i;
+ int eargc = 0, mon = 0;
+ int have_to_flag = 0, inverse = 0;
int use_env_locale = 0;
/* This is just to check that pj_init() is locale-safe */
@@ -185,6 +323,11 @@ int main(int argc, char **argv) {
if (getenv("PROJ_USE_ENV_LOCALE") != nullptr)
use_env_locale = 1;
+ /* Enable compatibility mode for init=epsg:XXXX by default */
+ if (getenv("PROJ_USE_PROJ4_INIT_RULES") == nullptr) {
+ proj_context_use_proj4_init_rules(nullptr, true);
+ }
+
if ((emess_dat.Prog_name = strrchr(*argv, DIR_CHAR)) != nullptr)
++emess_dat.Prog_name;
else
@@ -194,9 +337,27 @@ int main(int argc, char **argv) {
(void)fprintf(stderr, usage, pj_get_release(), emess_dat.Prog_name);
exit(0);
}
+
+ // First pass to check if we have "cs2cs [-bla]* <SRC> <DEST>" syntax
+ int countNonOptionArg = 0;
+ for (int i = 1; i < argc; i++) {
+ if (argv[i][0] == '-') {
+ if (argv[i][1] == '\0') {
+ countNonOptionArg++;
+ }
+ } else {
+ if (strcmp(argv[i], "+to") == 0) {
+ countNonOptionArg = -1;
+ break;
+ }
+ countNonOptionArg++;
+ }
+ }
+ const bool isSrcDestSyntax = (countNonOptionArg == 2);
+
/* process run line arguments */
while (--argc > 0) { /* collect run line arguments */
- if (**++argv == '-')
+ if (**++argv == '-') {
for (arg = *argv;;) {
switch (*++arg) {
case '\0': /* position of "stdin" */
@@ -325,21 +486,28 @@ int main(int argc, char **argv) {
}
break;
}
- else if (strcmp(*argv, "+to") == 0) {
+ } else if (isSrcDestSyntax) {
+ if (fromStr.empty())
+ fromStr = *argv;
+ else
+ toStr = *argv;
+ } else if (strcmp(*argv, "+to") == 0) {
have_to_flag = 1;
} else if (**argv == '+') { /* + argument */
if (have_to_flag) {
- if (to_argc < MAX_PARGS)
- to_argv[to_argc++] = *argv + 1;
- else
- emess(1, "overflowed + argument table");
+ if (!toStr.empty())
+ toStr += ' ';
+ toStr += *argv;
} else {
- if (from_argc < MAX_PARGS)
- from_argv[from_argc++] = *argv + 1;
- else
- emess(1, "overflowed + argument table");
+ if (!fromStr.empty())
+ fromStr += ' ';
+ fromStr += *argv;
}
+ } else if (!have_to_flag) {
+ fromStr = *argv;
+ } else if (toStr.empty()) {
+ toStr = *argv;
} else /* assumed to be input file name(s) */
eargv[eargc++] = *argv;
}
@@ -351,17 +519,7 @@ int main(int argc, char **argv) {
* coordinate systems.
*/
if (inverse) {
- int argcount;
-
- for (i = 0; i < MAX_PARGS; i++) {
- arg = from_argv[i];
- from_argv[i] = to_argv[i];
- to_argv[i] = arg;
- }
-
- argcount = from_argc;
- from_argc = to_argc;
- to_argc = argcount;
+ std::swap(fromStr, toStr);
}
if (use_env_locale) {
@@ -369,48 +527,64 @@ int main(int argc, char **argv) {
setlocale(LC_ALL, "");
}
- if (from_argc == 0 && to_argc != 0) {
- /* we will generate the from proj as the latlong of the +to in a bit */
- } else if (!(fromProj = pj_init(from_argc, from_argv))) {
- printf("Using from definition: ");
- for (i = 0; i < from_argc; i++)
- printf("%s ", from_argv[i]);
- printf("\n");
-
- emess(3, "projection initialization failure\ncause: %s",
- pj_strerrno(pj_errno));
+ if (fromStr.empty() && toStr.empty()) {
+ emess(3, "missing source and target coordinate systems");
}
- if (to_argc == 0) {
- if (!(toProj = pj_latlong_from_proj(fromProj))) {
- printf("Using to definition: ");
- for (i = 0; i < to_argc; i++)
- printf("%s ", to_argv[i]);
- printf("\n");
+ const char *const optionsProj4Mode[] = {"USE_PROJ4_INIT_RULES=YES",
+ nullptr};
+ const char *const *optionsImportCRS =
+ proj_context_get_use_proj4_init_rules(nullptr) ? optionsProj4Mode
+ : nullptr;
+
+ PJ_OBJ *src = nullptr;
+ if (!fromStr.empty()) {
+ bool ignored;
+ src = instanciate_crs(fromStr, optionsImportCRS, srcIsGeog,
+ srcToRadians, ignored);
+ if (!src) {
+ emess(3, "cannot instanciate source coordinate system");
+ }
+ }
- emess(3, "projection initialization failure\ncause: %s",
- pj_strerrno(pj_errno));
+ PJ_OBJ *dst = nullptr;
+ if (!toStr.empty()) {
+ dst = instanciate_crs(toStr, optionsImportCRS, destIsGeog,
+ destToRadians, destIsLatLong);
+ if (!dst) {
+ emess(3, "cannot instanciate target coordinate system");
}
- } else if (!(toProj = pj_init(to_argc, to_argv))) {
- printf("Using to definition: ");
- for (i = 0; i < to_argc; i++)
- printf("%s ", to_argv[i]);
- printf("\n");
+ }
- emess(3, "projection initialization failure\ncause: %s",
- pj_strerrno(pj_errno));
+ if (toStr.empty()) {
+ assert(src);
+ toStr = get_geog_crs_proj_string_from_proj_crs(src, destToRadians,
+ destIsLatLong);
+ if (toStr.empty()) {
+ emess(3,
+ "missing target CRS and source CRS is not a projected CRS");
+ }
+ destIsGeog = true;
+ } else if (fromStr.empty()) {
+ assert(dst);
+ bool ignored;
+ fromStr =
+ get_geog_crs_proj_string_from_proj_crs(dst, srcToRadians, ignored);
+ if (fromStr.empty()) {
+ emess(3,
+ "missing source CRS and target CRS is not a projected CRS");
+ }
+ srcIsGeog = true;
}
- if (from_argc == 0 && toProj != nullptr) {
- if (!(fromProj = pj_latlong_from_proj(toProj))) {
- printf("Using to definition: ");
- for (i = 0; i < to_argc; i++)
- printf("%s ", to_argv[i]);
- printf("\n");
+ proj_obj_unref(src);
+ proj_obj_unref(dst);
- emess(3, "projection initialization failure\ncause: %s",
- pj_strerrno(pj_errno));
- }
+ transformation = proj_create_crs_to_crs(nullptr, fromStr.c_str(),
+ toStr.c_str(), nullptr);
+ if (!transformation) {
+ emess(3, "cannot initialize transformation\ncause: %s",
+ pj_strerrno(pj_errno));
}
if (use_env_locale) {
@@ -420,19 +594,19 @@ int main(int argc, char **argv) {
if (mon) {
printf("%c ---- From Coordinate System ----\n", tag);
- pj_pr_list(fromProj);
+ printf("%s\n", fromStr.c_str());
printf("%c ---- To Coordinate System ----\n", tag);
- pj_pr_list(toProj);
+ printf("%s\n", toStr.c_str());
}
/* set input formatting control */
- if (!fromProj->is_latlong)
+ if (!srcIsGeog)
informat = strtod;
else {
informat = dmstor;
}
- if (!toProj->is_latlong && !oform)
+ if (!destIsGeog && !oform)
oform = "%.2f";
/* process input file list */
@@ -454,8 +628,7 @@ int main(int argc, char **argv) {
emess_dat.File_name = nullptr;
}
- pj_free(fromProj);
- pj_free(toProj);
+ proj_destroy(transformation);
pj_deallocate_grids();
diff --git a/src/io.cpp b/src/io.cpp
index 0d220e13..15df312b 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -4422,6 +4422,7 @@ struct PROJStringFormatter::Private {
bool useETMercForTMerc_ = false;
bool useETMercForTMercSet_ = false;
bool addNoDefs_ = true;
+ bool coordOperationOptimizations_ = false;
std::string result_{};
@@ -4504,6 +4505,12 @@ const std::string &PROJStringFormatter::toString() const {
step.paramValues[6].equals("s", "0") &&
step.paramValues[7].keyEquals("convention")))) {
iter = d->steps_.erase(iter);
+ } else if (d->coordOperationOptimizations_ &&
+ step.name == "unitconvert" && paramCount == 2 &&
+ step.paramValues[0].keyEquals("xy_in") &&
+ step.paramValues[1].keyEquals("xy_out") &&
+ step.paramValues[0].value == step.paramValues[1].value) {
+ iter = d->steps_.erase(iter);
} else {
++iter;
}
@@ -4925,6 +4932,12 @@ bool PROJStringFormatter::getUseETMercForTMerc(bool &settingSetOut) const {
// ---------------------------------------------------------------------------
+void PROJStringFormatter::setCoordinateOperationOptimizations(bool enable) {
+ d->coordOperationOptimizations_ = enable;
+}
+
+// ---------------------------------------------------------------------------
+
void PROJStringFormatter::Private::appendToResult(const char *str) {
if (!result_.empty()) {
result_ += " ";
@@ -4984,11 +4997,13 @@ PROJStringSyntaxParser(const std::string &projString, std::vector<Step> &steps,
}
prevWasStep = false;
} else if (starts_with(word, "proj=")) {
+ auto stepName = word.substr(strlen("proj="));
if (prevWasInit) {
- throw ParsingException("+init= found at unexpected place");
+ steps.back() = Step();
+ prevWasInit = false;
+ } else {
+ steps.push_back(Step());
}
- auto stepName = word.substr(strlen("proj="));
- steps.push_back(Step());
steps.back().name = stepName;
steps.back().inverted = inverted;
prevWasStep = false;
@@ -5041,7 +5056,7 @@ void PROJStringFormatter::ingestPROJString(
std::string vto_meter;
PROJStringSyntaxParser(str, steps, d->globalParamValues_, title, vunits,
vto_meter);
- d->steps_.insert(d->steps_.begin(), steps.begin(), steps.end());
+ d->steps_.insert(d->steps_.end(), steps.begin(), steps.end());
}
// ---------------------------------------------------------------------------
@@ -5633,11 +5648,12 @@ PROJStringParser::Private::buildPrimeMeridian(const Step &step) {
PrimeMeridianNNPtr pm = PrimeMeridian::GREENWICH;
const auto &pmStr = getParamValue(step, "pm");
if (!pmStr.empty()) {
- try {
- double pmValue = c_locale_stod(pmStr);
+ char *end;
+ double pmValue = dmstor(pmStr.c_str(), &end) * RAD_TO_DEG;
+ if (pmValue != HUGE_VAL && *end == '\0') {
pm = PrimeMeridian::create(createMapWithUnknownName(),
Angle(pmValue));
- } catch (const std::invalid_argument &) {
+ } else {
bool found = false;
if (pmStr == "paris") {
found = true;
@@ -5650,9 +5666,8 @@ PROJStringParser::Private::buildPrimeMeridian(const Step &step) {
found = true;
std::string name = static_cast<char>(::toupper(pmStr[0])) +
pmStr.substr(1);
- double pmValue =
- dmstor(proj_prime_meridians[i].defn, nullptr) *
- RAD_TO_DEG;
+ pmValue = dmstor(proj_prime_meridians[i].defn, nullptr) *
+ RAD_TO_DEG;
pm = PrimeMeridian::create(
PropertyMap().set(IdentifiedObject::NAME_KEY, name),
Angle(pmValue));
@@ -5681,12 +5696,20 @@ PROJStringParser::Private::buildDatum(const Step &step,
const auto &ellpsStr = getParamValue(step, "ellps");
const auto &datumStr = getParamValue(step, "datum");
+ const auto &RStr = getParamValue(step, "R");
const auto &aStr = getParamValue(step, "a");
const auto &bStr = getParamValue(step, "b");
const auto &rfStr = getParamValue(step, "rf");
const auto &fStr = getParamValue(step, "f");
- const auto &RStr = getParamValue(step, "R");
+ const auto &esStr = getParamValue(step, "es");
+ const auto &eStr = getParamValue(step, "e");
+ double a = -1.0;
+ double b = -1.0;
+ double rf = -1.0;
const util::optional<std::string> optionalEmptyString{};
+ const bool numericParamPresent =
+ !RStr.empty() || !aStr.empty() || !bStr.empty() || !rfStr.empty() ||
+ !fStr.empty() || !esStr.empty() || !eStr.empty();
PrimeMeridianNNPtr pm(buildPrimeMeridian(step));
PropertyMap grfMap;
@@ -5709,104 +5732,151 @@ PROJStringParser::Private::buildDatum(const Step &step,
}
};
+ // R take precedence
+ if (!RStr.empty()) {
+ double R;
+ try {
+ R = c_locale_stod(RStr);
+ } catch (const std::invalid_argument &) {
+ throw ParsingException("Invalid R value");
+ }
+ auto ellipsoid = Ellipsoid::createSphere(createMapWithUnknownName(),
+ Length(R), guessBodyName(R));
+ return GeodeticReferenceFrame::create(
+ grfMap.set(IdentifiedObject::NAME_KEY,
+ title.empty() ? "unknown" : title.c_str()),
+ ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
+ }
+
if (!datumStr.empty()) {
- if (datumStr == "WGS84") {
- return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6326);
- } else if (datumStr == "NAD83") {
- return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6269);
- } else if (datumStr == "NAD27") {
- return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6267);
- } else {
+ auto l_datum = [&datumStr, &overridePmIfNeeded, &grfMap,
+ &optionalEmptyString, &pm]() {
+ if (datumStr == "WGS84") {
+ return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6326);
+ } else if (datumStr == "NAD83") {
+ return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6269);
+ } else if (datumStr == "NAD27") {
+ return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6267);
+ } else {
- for (const auto &datumDesc : datumDescs) {
- if (datumStr == datumDesc.projName) {
- (void)datumDesc.gcsName; // to please cppcheck
- (void)datumDesc.gcsCode; // to please cppcheck
- auto ellipsoid = Ellipsoid::createFlattenedSphere(
- grfMap
- .set(IdentifiedObject::NAME_KEY,
- datumDesc.ellipsoidName)
- .set(Identifier::CODESPACE_KEY, Identifier::EPSG)
- .set(Identifier::CODE_KEY, datumDesc.ellipsoidCode),
- Length(datumDesc.a), Scale(datumDesc.rf));
- return GeodeticReferenceFrame::create(
- grfMap
- .set(IdentifiedObject::NAME_KEY,
- datumDesc.datumName)
- .set(Identifier::CODESPACE_KEY, Identifier::EPSG)
- .set(Identifier::CODE_KEY, datumDesc.datumCode),
- ellipsoid, optionalEmptyString, pm);
+ for (const auto &datumDesc : datumDescs) {
+ if (datumStr == datumDesc.projName) {
+ (void)datumDesc.gcsName; // to please cppcheck
+ (void)datumDesc.gcsCode; // to please cppcheck
+ auto ellipsoid = Ellipsoid::createFlattenedSphere(
+ grfMap
+ .set(IdentifiedObject::NAME_KEY,
+ datumDesc.ellipsoidName)
+ .set(Identifier::CODESPACE_KEY,
+ Identifier::EPSG)
+ .set(Identifier::CODE_KEY,
+ datumDesc.ellipsoidCode),
+ Length(datumDesc.a), Scale(datumDesc.rf));
+ return GeodeticReferenceFrame::create(
+ grfMap
+ .set(IdentifiedObject::NAME_KEY,
+ datumDesc.datumName)
+ .set(Identifier::CODESPACE_KEY,
+ Identifier::EPSG)
+ .set(Identifier::CODE_KEY, datumDesc.datumCode),
+ ellipsoid, optionalEmptyString, pm);
+ }
}
}
+ throw ParsingException("unknown datum " + datumStr);
+ }();
+ if (!numericParamPresent) {
+ return l_datum;
}
- throw ParsingException("unknown datum " + datumStr);
+ a = l_datum->ellipsoid()->semiMajorAxis().getSIValue();
+ rf = l_datum->ellipsoid()->computedInverseFlattening();
}
else if (!ellpsStr.empty()) {
- if (ellpsStr == "WGS84") {
- return GeodeticReferenceFrame::create(
- grfMap.set(IdentifiedObject::NAME_KEY,
- title.empty() ? "Unknown based on WGS84 ellipsoid"
- : title.c_str()),
- Ellipsoid::WGS84, optionalEmptyString, pm);
- } else if (ellpsStr == "GRS80") {
- return GeodeticReferenceFrame::create(
- grfMap.set(IdentifiedObject::NAME_KEY,
- title.empty() ? "Unknown based on GRS80 ellipsoid"
- : title.c_str()),
- Ellipsoid::GRS1980, optionalEmptyString, pm);
- } else {
- auto proj_ellps = proj_list_ellps();
- for (int i = 0; proj_ellps[i].id != nullptr; i++) {
- if (ellpsStr == proj_ellps[i].id) {
- assert(strncmp(proj_ellps[i].major, "a=", 2) == 0);
- const double a_iter =
- c_locale_stod(proj_ellps[i].major + 2);
- EllipsoidPtr ellipsoid;
- PropertyMap ellpsMap;
- if (strncmp(proj_ellps[i].ell, "b=", 2) == 0) {
- const double b_iter =
- c_locale_stod(proj_ellps[i].ell + 2);
- ellipsoid = Ellipsoid::createTwoAxis(
- ellpsMap.set(IdentifiedObject::NAME_KEY,
- proj_ellps[i].name),
- Length(a_iter), Length(b_iter))
- .as_nullable();
- } else {
- assert(strncmp(proj_ellps[i].ell, "rf=", 3) == 0);
- const double rf_iter =
- c_locale_stod(proj_ellps[i].ell + 3);
- ellipsoid = Ellipsoid::createFlattenedSphere(
- ellpsMap.set(IdentifiedObject::NAME_KEY,
- proj_ellps[i].name),
- Length(a_iter), Scale(rf_iter))
- .as_nullable();
+ auto l_datum = [&ellpsStr, &title, &grfMap, &optionalEmptyString,
+ &pm]() {
+ if (ellpsStr == "WGS84") {
+ return GeodeticReferenceFrame::create(
+ grfMap.set(IdentifiedObject::NAME_KEY,
+ title.empty()
+ ? "Unknown based on WGS84 ellipsoid"
+ : title.c_str()),
+ Ellipsoid::WGS84, optionalEmptyString, pm);
+ } else if (ellpsStr == "GRS80") {
+ return GeodeticReferenceFrame::create(
+ grfMap.set(IdentifiedObject::NAME_KEY,
+ title.empty()
+ ? "Unknown based on GRS80 ellipsoid"
+ : title.c_str()),
+ Ellipsoid::GRS1980, optionalEmptyString, pm);
+ } else {
+ auto proj_ellps = proj_list_ellps();
+ for (int i = 0; proj_ellps[i].id != nullptr; i++) {
+ if (ellpsStr == proj_ellps[i].id) {
+ assert(strncmp(proj_ellps[i].major, "a=", 2) == 0);
+ const double a_iter =
+ c_locale_stod(proj_ellps[i].major + 2);
+ EllipsoidPtr ellipsoid;
+ PropertyMap ellpsMap;
+ if (strncmp(proj_ellps[i].ell, "b=", 2) == 0) {
+ const double b_iter =
+ c_locale_stod(proj_ellps[i].ell + 2);
+ ellipsoid =
+ Ellipsoid::createTwoAxis(
+ ellpsMap.set(IdentifiedObject::NAME_KEY,
+ proj_ellps[i].name),
+ Length(a_iter), Length(b_iter))
+ .as_nullable();
+ } else {
+ assert(strncmp(proj_ellps[i].ell, "rf=", 3) == 0);
+ const double rf_iter =
+ c_locale_stod(proj_ellps[i].ell + 3);
+ ellipsoid =
+ Ellipsoid::createFlattenedSphere(
+ ellpsMap.set(IdentifiedObject::NAME_KEY,
+ proj_ellps[i].name),
+ Length(a_iter), Scale(rf_iter))
+ .as_nullable();
+ }
+ return GeodeticReferenceFrame::create(
+ grfMap.set(IdentifiedObject::NAME_KEY,
+ title.empty()
+ ? std::string("Unknown based on ") +
+ proj_ellps[i].name +
+ " ellipsoid"
+ : title),
+ NN_NO_CHECK(ellipsoid), optionalEmptyString, pm);
}
- return GeodeticReferenceFrame::create(
- grfMap.set(IdentifiedObject::NAME_KEY,
- title.empty()
- ? std::string("Unknown based on ") +
- proj_ellps[i].name + " ellipsoid"
- : title),
- NN_NO_CHECK(ellipsoid), optionalEmptyString, pm);
}
+ throw ParsingException("unknown ellipsoid " + ellpsStr);
}
- throw ParsingException("unknown ellipsoid " + ellpsStr);
+ }();
+ if (!numericParamPresent) {
+ return l_datum;
+ }
+ a = l_datum->ellipsoid()->semiMajorAxis().getSIValue();
+ if (l_datum->ellipsoid()->semiMinorAxis().has_value()) {
+ b = l_datum->ellipsoid()->semiMinorAxis()->getSIValue();
+ } else {
+ rf = l_datum->ellipsoid()->computedInverseFlattening();
}
}
- else if (!aStr.empty() && !bStr.empty()) {
- double a;
+ if (!aStr.empty()) {
try {
a = c_locale_stod(aStr);
} catch (const std::invalid_argument &) {
throw ParsingException("Invalid a value");
}
- double b;
- try {
- b = c_locale_stod(bStr);
- } catch (const std::invalid_argument &) {
- throw ParsingException("Invalid b value");
+ }
+
+ if (a > 0 && (b > 0 || !bStr.empty())) {
+ if (!bStr.empty()) {
+ try {
+ b = c_locale_stod(bStr);
+ } catch (const std::invalid_argument &) {
+ throw ParsingException("Invalid b value");
+ }
}
auto ellipsoid =
Ellipsoid::createTwoAxis(createMapWithUnknownName(), Length(a),
@@ -5818,18 +5888,13 @@ PROJStringParser::Private::buildDatum(const Step &step,
ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
}
- else if (!aStr.empty() && !rfStr.empty()) {
- double a;
- try {
- a = c_locale_stod(aStr);
- } catch (const std::invalid_argument &) {
- throw ParsingException("Invalid a value");
- }
- double rf;
- try {
- rf = c_locale_stod(rfStr);
- } catch (const std::invalid_argument &) {
- throw ParsingException("Invalid rf value");
+ else if (a > 0 && (rf >= 0 || !rfStr.empty())) {
+ if (!rfStr.empty()) {
+ try {
+ rf = c_locale_stod(rfStr);
+ } catch (const std::invalid_argument &) {
+ throw ParsingException("Invalid rf value");
+ }
}
auto ellipsoid = Ellipsoid::createFlattenedSphere(
createMapWithUnknownName(), Length(a), Scale(rf),
@@ -5841,13 +5906,7 @@ PROJStringParser::Private::buildDatum(const Step &step,
ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
}
- else if (!aStr.empty() && !fStr.empty()) {
- double a;
- try {
- a = c_locale_stod(aStr);
- } catch (const std::invalid_argument &) {
- throw ParsingException("Invalid a value");
- }
+ else if (a > 0 && !fStr.empty()) {
double f;
try {
f = c_locale_stod(fStr);
@@ -5864,23 +5923,52 @@ PROJStringParser::Private::buildDatum(const Step &step,
ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
}
- else if (!RStr.empty()) {
- double R;
+ else if (a > 0 && !eStr.empty()) {
+ double e;
try {
- R = c_locale_stod(RStr);
+ e = c_locale_stod(eStr);
} catch (const std::invalid_argument &) {
- throw ParsingException("Invalid R value");
+ throw ParsingException("Invalid e value");
}
- auto ellipsoid = Ellipsoid::createSphere(createMapWithUnknownName(),
- Length(R), guessBodyName(R));
+ double alpha = asin(e); /* angular eccentricity */
+ double f = 1 - cos(alpha); /* = 1 - sqrt (1 - es); */
+ auto ellipsoid = Ellipsoid::createFlattenedSphere(
+ createMapWithUnknownName(), Length(a),
+ Scale(f != 0.0 ? 1.0 / f : 0.0), guessBodyName(a))
+ ->identify();
return GeodeticReferenceFrame::create(
grfMap.set(IdentifiedObject::NAME_KEY,
title.empty() ? "unknown" : title.c_str()),
ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
}
- if (!aStr.empty() && bStr.empty() && rfStr.empty()) {
- throw ParsingException("a found, but b, f or rf missing");
+ else if (a > 0 && !esStr.empty()) {
+ double es;
+ try {
+ es = c_locale_stod(esStr);
+ } catch (const std::invalid_argument &) {
+ throw ParsingException("Invalid es value");
+ }
+ double f = 1 - sqrt(1 - es);
+ auto ellipsoid = Ellipsoid::createFlattenedSphere(
+ createMapWithUnknownName(), Length(a),
+ Scale(f != 0.0 ? 1.0 / f : 0.0), guessBodyName(a))
+ ->identify();
+ return GeodeticReferenceFrame::create(
+ grfMap.set(IdentifiedObject::NAME_KEY,
+ title.empty() ? "unknown" : title.c_str()),
+ ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
+ }
+
+ // If only a is specified, create a sphere
+ if (a > 0 && bStr.empty() && rfStr.empty() && eStr.empty() &&
+ esStr.empty()) {
+ auto ellipsoid = Ellipsoid::createSphere(createMapWithUnknownName(),
+ Length(a), guessBodyName(a));
+ return GeodeticReferenceFrame::create(
+ grfMap.set(IdentifiedObject::NAME_KEY,
+ title.empty() ? "unknown" : title.c_str()),
+ ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
}
if (!bStr.empty() && aStr.empty()) {
@@ -5895,6 +5983,14 @@ PROJStringParser::Private::buildDatum(const Step &step,
throw ParsingException("f found, but a missing");
}
+ if (!eStr.empty() && aStr.empty()) {
+ throw ParsingException("e found, but a missing");
+ }
+
+ if (!esStr.empty() && aStr.empty()) {
+ throw ParsingException("es found, but a missing");
+ }
+
return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6326);
}
@@ -6064,6 +6160,22 @@ PROJStringParser::Private::buildEllipsoidalCS(int iStep, int iUnitConvert,
// ---------------------------------------------------------------------------
+static double getNumericValue(const std::string &paramValue,
+ bool *pHasError = nullptr) {
+ try {
+ double value = c_locale_stod(paramValue);
+ if (pHasError)
+ *pHasError = false;
+ return value;
+ } catch (const std::invalid_argument &) {
+ if (pHasError)
+ *pHasError = true;
+ return 0.0;
+ }
+}
+
+// ---------------------------------------------------------------------------
+
GeographicCRSNNPtr
PROJStringParser::Private::buildGeographicCRS(int iStep, int iUnitConvert,
int iAxisSwap, bool ignoreVUnits,
@@ -6077,7 +6189,10 @@ PROJStringParser::Private::buildGeographicCRS(int iStep, int iUnitConvert,
auto props = PropertyMap().set(IdentifiedObject::NAME_KEY,
title.empty() ? "unknown" : title);
- if (l_isGeographicStep && hasParamValue(step, "wktext")) {
+ if (l_isGeographicStep &&
+ (hasParamValue(step, "wktext") ||
+ hasParamValue(step, "lon_wrap") | hasParamValue(step, "geoc") ||
+ getNumericValue(getParamValue(step, "lon_0")) != 0.0)) {
props.set("EXTENSION_PROJ4", projString_);
}
@@ -6214,22 +6329,6 @@ static double getAngularValue(const std::string &paramValue,
// ---------------------------------------------------------------------------
-static double getNumericValue(const std::string &paramValue,
- bool *pHasError = nullptr) {
- try {
- double value = c_locale_stod(paramValue);
- if (pHasError)
- *pHasError = false;
- return value;
- } catch (const std::invalid_argument &) {
- if (pHasError)
- *pHasError = true;
- return 0.0;
- }
-}
-
-// ---------------------------------------------------------------------------
-
CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
int iStep, GeographicCRSNNPtr geogCRS, int iUnitConvert, int iAxisSwap) {
auto &step = steps_[iStep];
@@ -6285,7 +6384,10 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
mapping =
getMapping(PROJ_WKT2_NAME_METHOD_GEOSTATIONARY_SATELLITE_SWEEP_Y);
} else if (step.name == "omerc") {
- if (hasParamValue(step, "no_uoff") || hasParamValue(step, "no_off")) {
+ if (hasParamValue(step, "no_rot")) {
+ mapping = nullptr;
+ } else if (hasParamValue(step, "no_uoff") ||
+ hasParamValue(step, "no_off")) {
mapping =
getMapping(EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_A);
} else if (hasParamValue(step, "lat_1") &&
@@ -6839,6 +6941,21 @@ PROJStringParser::Private::buildMolodenskyTransformation(
// ---------------------------------------------------------------------------
+//! @cond Doxygen_Suppress
+static const metadata::ExtentPtr nullExtent{};
+
+static const metadata::ExtentPtr &getExtent(const crs::CRS *crs) {
+ const auto &domains = crs->domains();
+ if (!domains.empty()) {
+ return domains[0]->domainOfValidity();
+ }
+ return nullExtent;
+}
+
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
/** \brief Instanciate a sub-class of BaseObject from a PROJ string.
* @throw ParsingException
*/
@@ -6909,30 +7026,36 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
}
auto obj =
createFromUserInput(d->steps_[0].name, d->dbContext_, true);
- auto geogCRS = dynamic_cast<GeographicCRS *>(obj.get());
- if (geogCRS) {
- // Override with longitude latitude in radian
- return GeographicCRS::create(
- PropertyMap().set(IdentifiedObject::NAME_KEY,
- geogCRS->nameStr()),
- geogCRS->datum(), geogCRS->datumEnsemble(),
- EllipsoidalCS::createLongitudeLatitude(
- UnitOfMeasure::RADIAN));
- }
- auto projCRS = dynamic_cast<ProjectedCRS *>(obj.get());
- if (projCRS) {
- // Override with easting northing orer
- auto conv = projCRS->derivingConversionRef();
- if (conv->method()->getEPSGCode() !=
- EPSG_CODE_METHOD_TRANSVERSE_MERCATOR_SOUTH_ORIENTATED) {
- return ProjectedCRS::create(
- PropertyMap().set(IdentifiedObject::NAME_KEY,
- projCRS->nameStr()),
- projCRS->baseCRS(), conv,
- CartesianCS::createEastingNorthing(
- projCRS->coordinateSystem()
- ->axisList()[0]
- ->unit()));
+ auto crs = dynamic_cast<CRS *>(obj.get());
+ if (crs) {
+ PropertyMap properties;
+ properties.set(IdentifiedObject::NAME_KEY, crs->nameStr());
+ const auto &extent = getExtent(crs);
+ if (extent) {
+ properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+ NN_NO_CHECK(extent));
+ }
+ auto geogCRS = dynamic_cast<GeographicCRS *>(crs);
+ if (geogCRS) {
+ // Override with longitude latitude in radian
+ return GeographicCRS::create(
+ properties, geogCRS->datum(), geogCRS->datumEnsemble(),
+ EllipsoidalCS::createLongitudeLatitude(
+ UnitOfMeasure::RADIAN));
+ }
+ auto projCRS = dynamic_cast<ProjectedCRS *>(crs);
+ if (projCRS) {
+ // Override with easting northing order
+ const auto &conv = projCRS->derivingConversionRef();
+ if (conv->method()->getEPSGCode() !=
+ EPSG_CODE_METHOD_TRANSVERSE_MERCATOR_SOUTH_ORIENTATED) {
+ return ProjectedCRS::create(
+ properties, projCRS->baseCRS(), conv,
+ CartesianCS::createEastingNorthing(
+ projCRS->coordinateSystem()
+ ->axisList()[0]
+ ->unit()));
+ }
}
}
return obj;
diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c
index 31e0124e..45ce5c8e 100644
--- a/src/pj_apply_gridshift.c
+++ b/src/pj_apply_gridshift.c
@@ -349,7 +349,7 @@ LP proj_hgrid_apply(PJ *P, LP lp, PJ_DIRECTION direction) {
out = nad_cvt(lp, inverse, ct);
if (out.lam == HUGE_VAL || out.phi == HUGE_VAL)
- pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+ pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
return out;
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 1a970374..38443f07 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -103,7 +103,6 @@ static PJ_COORD fwd_prepare (PJ *P, PJ_COORD coo) {
}
-
static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
switch (OUTPUT_UNITS) {
@@ -138,6 +137,14 @@ static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
case PJ_IO_UNITS_ANGULAR:
coo.lpz.z = P->vfr_meter * (coo.lpz.z + P->z0);
+
+ if( P->is_long_wrap_set ) {
+ if( coo.lpz.lam != HUGE_VAL ) {
+ coo.lpz.lam = P->long_wrap_center +
+ adjlon(coo.lpz.lam - P->long_wrap_center);
+ }
+ }
+
break;
}
diff --git a/src/proj.h b/src/proj.h
index 4b599eba..31cd730c 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -1381,6 +1381,19 @@ PJ_OBJ PROJ_DLL *proj_obj_crs_get_sub_crs(PJ_OBJ *crs, int index);
PJ_OBJ PROJ_DLL *proj_obj_crs_create_bound_crs_to_WGS84(PJ_OBJ *crs);
+PJ_OBJ PROJ_DLL *proj_obj_crs_get_coordinate_system(PJ_OBJ *crs);
+
+const char PROJ_DLL *proj_obj_cs_get_type(PJ_OBJ* cs);
+
+int PROJ_DLL proj_obj_cs_get_axis_count(PJ_OBJ *cs);
+
+int PROJ_DLL proj_obj_cs_get_axis_info(PJ_OBJ *cs, int index,
+ const char **pName,
+ const char **pAbbrev,
+ const char **pDirection,
+ double *pUnitConvFactor,
+ const char **pUnitName);
+
PJ_OBJ PROJ_DLL *proj_obj_get_ellipsoid(PJ_OBJ *obj);
int PROJ_DLL proj_obj_ellipsoid_get_parameters(PJ_OBJ *ellipsoid,
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 72e1a2d6..b7b500a7 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -794,7 +794,12 @@ PJ *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *srid_from, const char
return NULL;
}
- P = proj_create(ctx, proj_string);
+ if( proj_string[0] == '\0' ) {
+ /* Null transform ? */
+ P = proj_create(ctx, "proj=affine");
+ } else {
+ P = proj_create(ctx, proj_string);
+ }
proj_obj_unref(op);
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 272e503c..00c9d4a0 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -709,4 +709,14 @@ direction reverse
accept 0 0 0 0
expect failure
+-------------------------------------------------------------------------------
+# Test lon_wrap
+operation +proj=longlat +ellps=WGS84 +lon_wrap=180
+-------------------------------------------------------------------------------
+direction forward
+accept -1 10 0
+expect 359 10 0
+
+-------------------------------------------------------------------------------
+
</gie>
diff --git a/test/old/proj_outIGNF.dist b/test/old/proj_outIGNF.dist
index 45112f60..611144b8 100644
--- a/test/old/proj_outIGNF.dist
+++ b/test/old/proj_outIGNF.dist
@@ -1,24 +1,24 @@
-+init=./IGNF:NTFG +to +init=./IGNF:RGF93G
++init=IGNF:NTFG +to +init=IGNF:RGF93G
3.300866856 43.4477976569 0.0000 3d18'0.915"E 43d26'52.077"N 0.000
-+init=./IGNF:LAMBE +to +init=./IGNF:LAMB93
- 600000.0000 2600545.4523 0.0000 652760.737 7033791.243 0.000
- 135638.3592 2418760.4094 0.0000 187194.062 6855928.882 0.000
- 998137.3947 2413822.2844 0.0000 1049052.258 6843776.562 0.000
++init=IGNF:LAMBE +to +init=IGNF:LAMB93
+ 600000.0000 2600545.4523 0.0000 652759.036 7033588.609 0.000
+ 135638.3592 2418760.4094 0.0000 187444.148 6856142.911 0.000
+ 998137.3947 2413822.2844 0.0000 1048843.997 6843923.913 0.000
600000.0000 2200000.0000 0.0000 649398.872 6633524.191 0.000
- 311552.5340 1906457.4840 0.0000 358799.172 6342652.486 0.000
- 960488.4138 1910172.8812 0.0000 1007068.686 6340907.237 0.000
- 600000.0000 1699510.8340 0.0000 645204.279 6133556.746 0.000
-1203792.5981 626873.17210 0.0000 1238875.764 5057405.016 0.000
-+init=./IGNF:LAMBE +to +init=./IGNF:GEOPORTALFXX
- 600000.0000 2600545.4523 0.0000 179040.148 5610495.275 0.000
- 135638.3592 2418760.4094 0.0000 -303729.363 5410118.356 0.000
- 998137.3947 2413822.2844 0.0000 592842.792 5410120.554 0.000
+ 311552.5340 1906457.4840 0.0000 358593.374 6342647.465 0.000
+ 960488.4138 1910172.8812 0.0000 1007324.119 6340956.093 0.000
+ 600000.0000 1699510.8340 0.0000 645201.753 6133255.515 0.000
+1203792.5981 626873.17210 0.0000 * * inf
++init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX
+ 600000.0000 2600545.4523 0.0000 179040.150 5610292.766 0.000
+ 135638.3592 2418760.4094 0.0000 -303490.059 5410353.890 0.000
+ 998137.3947 2413822.2844 0.0000 592635.926 5410280.335 0.000
600000.0000 2200000.0000 0.0000 179041.670 5209746.080 0.000
- 311552.5340 1906457.4840 0.0000 -96825.465 4909184.136 0.000
- 960488.4138 1910172.8812 0.0000 523880.019 4909191.141 0.000
- 600000.0000 1699510.8340 0.0000 179047.633 4708817.007 0.000
-1203792.5981 626873.17210 0.0000 658287.395 3623739.237 0.000
-+init=./IGNF:RGF93G +to +init=./IGNF:GEOPORTALFXX
+ 311552.5340 1906457.4840 0.0000 -97021.878 4909167.981 0.000
+ 960488.4138 1910172.8812 0.0000 524126.466 4909227.598 0.000
+ 600000.0000 1699510.8340 0.0000 179047.637 4708515.623 0.000
+1203792.5981 626873.17210 0.0000 * * inf
++init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX
2d20'11.4239243" 50d23'59.7718445" 0.0 179040.151 5610495.281 0.000
-3d57'49.4051448" 48d35'59.7121716" 0.0 -303729.365 5410118.352 0.000
7d44'12.1439796" 48d35'59.7832558" 0.0 592842.794 5410120.550 0.000
@@ -27,7 +27,7 @@
6d50'12.2276489" 44d06'00.0517019" 0.0 523880.022 4909191.143 0.000
2d20'11.7754730" 42d18'00.0824436" 0.0 179047.634 4708817.010 0.000
9d32'12.6680218" 41d24'00.3542556" 0.0 730783.054 4608637.873 0.000
-+init=./IGNF:RGF93G +to +init=./IGNF:MILLER
++init=IGNF:RGF93G +to +init=IGNF:MILLER
2d20'11.4239243" 50d23'59.7718445" 0.0 260098.730 6140682.441 0.000
-3d57'49.4051448" 48d35'59.7121716" 0.0 -441239.699 5880610.004 0.000
7d44'12.1439796" 48d35'59.7832558" 0.0 861246.246 5880612.827 0.000
@@ -36,5 +36,5 @@
6d50'12.2276489" 44d06'00.0517019" 0.0 761061.291 5252498.745 0.000
2d20'11.7754730" 42d18'00.0824436" 0.0 260109.601 5009175.714 0.000
9d32'12.6680218" 41d24'00.3542556" 0.0 1061637.534 4889066.592 0.000
-+init=./IGNF:RGR92 +to +init=./IGNF:REUN47
-3356123.5400 1303218.3090 5247430.6050 3353421.833 1304074.314 5248935.607
++init=IGNF:RGR92 +to +init=IGNF:REUN47
+3356123.5400 1303218.3090 5247430.6050 3353420.949 1304075.021 5248935.144
diff --git a/test/old/testIGNF b/test/old/testIGNF
index f785641e..abbaed94 100755
--- a/test/old/testIGNF
+++ b/test/old/testIGNF
@@ -41,13 +41,7 @@ echo "============================================"
OUT=proj_outIGNF
-REMOVE_IGNF=NO
-if test ! -f IGNF; then
- cp ${PROJ_LIB}/IGNF .
- REMOVE_IGNF=YES
-fi
-
-INIT_FILE=./IGNF
+INIT_FILE=IGNF
RES="-f %.3f"
#
echo "doing tests into file ${OUT}, please wait"
@@ -162,10 +156,6 @@ $EXE +init=${INIT_FILE}:RGR92 +to +init=${INIT_FILE}:REUN47 -E $RES >>${OUT} <<E
3356123.5400 1303218.3090 5247430.6050
EOF
-if test "${REMOVE_IGNF}" = "YES"; then
- rm ./IGNF
-fi
-
#
# do 'diff' with distribution results
echo "diff ${OUT} with ${OUT}.dist"
diff --git a/test/old/testdatumfile b/test/old/testdatumfile
index 6cd25b76..6520453d 100755
--- a/test/old/testdatumfile
+++ b/test/old/testdatumfile
@@ -77,7 +77,7 @@ echo "edge or even a wee bit outside (#141)." >> ${OUT}
# Our test points are (1) right on mesh corner, (2) outside but within
# epsilon (3) inside a bit (4) outside by more than epsilon
#
-$EXE +proj=latlong +ellps=WGS84 +nadgrids=ntf_r93.gsb \
+$EXE +proj=latlong +ellps=WGS84 +nadgrids=ntf_r93.gsb,null \
+to +proj=latlong +datum=WGS84 \
-E -f "%.12f" >>${OUT} <<EOF
-5.5 52.0
@@ -87,7 +87,7 @@ $EXE +proj=latlong +ellps=WGS84 +nadgrids=ntf_r93.gsb \
EOF
#
$EXE +proj=latlong +datum=WGS84 \
- +to +proj=latlong +ellps=WGS84 +nadgrids=ntf_r93.gsb \
+ +to +proj=latlong +ellps=WGS84 +nadgrids=ntf_r93.gsb,null \
-E -f "%.12f" >>${OUT} <<EOF
-5.5 52.0
-5.5000000000001 52.0000000000001
diff --git a/test/old/testvarious b/test/old/testvarious
index dbaeb0ce..b3ba4040 100755
--- a/test/old/testvarious
+++ b/test/old/testvarious
@@ -523,7 +523,7 @@ EOF
#
echo "##############################################################" >> ${OUT}
echo "Test the Natural Earth Projection" >> ${OUT}
-$EXE +proj=latlong +a=1 +lon_0=0 \
+$EXE +proj=latlong +a=6371008.7714 +b=6371008.7714 \
+to +proj=natearth +a=6371008.7714 +b=6371008.7714 -f '%.'7'f' \
-E >>${OUT} <<EOF
0.0 0.0 0 0.0 0.0
@@ -554,7 +554,7 @@ $EXE +proj=latlong +a=1 +lon_0=0 \
EOF
echo "##############################################################" >> ${OUT}
echo "Test the Natural Earth II Projection" >> ${OUT}
-$EXE +proj=latlong +a=1 +lon_0=0 \
+$EXE +proj=latlong +a=6371008.7714 +b=6371008.7714 \
+to +proj=natearth2 +a=6371008.7714 +b=6371008.7714 -f '%.'7'f' \
-E >>${OUT} <<EOF
0.0 0.0 0 0.00000000 0.00000000
@@ -585,7 +585,7 @@ $EXE +proj=latlong +a=1 +lon_0=0 \
EOF
echo "##############################################################" >> ${OUT}
echo "Test the Compact Miller projection" >> ${OUT}
-$EXE +proj=latlong +a=1 +lon_0=0 \
+$EXE +proj=latlong +a=6371008.7714 +b=6371008.7714 \
+to +proj=comill +a=6371008.7714 +b=6371008.7714 -f '%.'7'f' \
-E >>${OUT} <<EOF
0.0 0.0 0 0.0 0.0
@@ -656,14 +656,14 @@ EOF
echo "##############################################################" >> ${OUT}
echo "Check inverse error handling with ob_tran (#225)" >> ${OUT}
$EXE +proj=ob_tran \
- +o_proj=moll +a=6378137 +es=0 +o_lon_p=LON_POLE +o_lat_p=LAT_POLE +lon_0=180 \
+ +o_proj=moll +a=6378137 +es=0 +o_lon_p=0 +o_lat_p=0 +lon_0=180 \
-E >>${OUT} <<EOF
300000 400000
20000000 30000000
EOF
echo "Test inverse handling" >> ${OUT}
$EXE -I +proj=ob_tran \
- +o_proj=moll +a=6378137 +es=0 +o_lon_p=LON_POLE +o_lat_p=LAT_POLE +lon_0=180 \
+ +o_proj=moll +a=6378137 +es=0 +o_lon_p=0 +o_lat_p=0 +lon_0=180 \
-E >>${OUT} <<EOF
10 20
EOF
@@ -939,6 +939,20 @@ $EXE -f '%0.3f' \
0 0 1000
EOF
+echo "##############################################################" >> ${OUT}
+echo "Test EPSG:4326 to EPSG:32631" >> ${OUT}
+# Input is latitude, longitude order
+$EXE EPSG:4326 +to EPSG:32631 -E >> ${OUT} <<EOF
+49 3 0
+EOF
+
+echo "##############################################################" >> ${OUT}
+echo "Test EPSG:32631 to EPSG:4326" >> ${OUT}
+# Input is latitude, longitude order
+$EXE EPSG:32631 EPSG:4326 -E >> ${OUT} <<EOF
+400000 5000000 0
+EOF
+
# Done!
# do 'diff' with distribution results
diff --git a/test/old/tv_out.dist b/test/old/tv_out.dist
index 222c0c5f..148d413d 100644
--- a/test/old/tv_out.dist
+++ b/test/old/tv_out.dist
@@ -139,7 +139,7 @@ Test healpix forward projection on ellipsoid
0 60.0 0.00000 3.35128 0.00000
0 -60.0 0.00000 -3.35128 0.00000
Test healpix inverse projection on ellipsoid
-0 0.7853981633974483 * * 0.00000
+0 0.7853981633974483 * * inf
-1.5707963267948966 0 -90.10072 0.00000 0.00000
0.0 0.0 0.00000 0.00000 0.00000
0.0 2.0547874222147415 0.00000 39.58811 0.00000
@@ -347,7 +347,7 @@ Test inverse calcofi projection
##############################################################
Check inverse error handling with ob_tran (#225)
300000 400000 42d45'22.377"W 85d35'28.083"N 0.000
-20000000 30000000 * * 0.000
+20000000 30000000 * * inf
Test inverse handling
10 20 -1384841.19 7581707.88 0.00
##############################################################
@@ -447,10 +447,16 @@ Test patterson inverse projection
20015114.352186374 -11409566.822831295 180.000 -90.000 0.000
##############################################################
Test Web Mercator to avoid issue #834 in the future
-487147.594520173 4934316.46263998 -10370728.796 5552839.742 -0.000
+487147.594520173 4934316.46263998 -10370728.796 5552839.742 0.000
##############################################################
Test vto_meter
0 0 1 0.000 0.000 1000.000
0 0 1 0.000 0.000 1000.000
0 0 1000 0.000 0.000 1.000
0 0 1000 0.000 0.000 1.000
+##############################################################
+Test EPSG:4326 to EPSG:32631
+49 3 0 500000.00 5427455.78 0.00
+##############################################################
+Test EPSG:32631 to EPSG:4326
+400000 5000000 0 45d8'47.014"N 1d43'40.681"E 0.000
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 8c9f114b..6c22cade 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -809,6 +809,22 @@ TEST_F(CApi, proj_obj_get_source_target_crs_transformation) {
// ---------------------------------------------------------------------------
+TEST_F(CApi, proj_obj_get_source_crs_of_projected_crs) {
+ auto crs = proj_obj_create_from_wkt(
+ m_ctxt,
+ createProjectedCRS()->exportToWKT(WKTFormatter::create().get()).c_str(),
+ nullptr);
+ ASSERT_NE(crs, nullptr);
+ ObjectKeeper keeper(crs);
+
+ auto sourceCRS = proj_obj_get_source_crs(crs);
+ ASSERT_NE(sourceCRS, nullptr);
+ ObjectKeeper keeper_sourceCRS(sourceCRS);
+ EXPECT_EQ(std::string(proj_obj_get_name(sourceCRS)), "WGS 84");
+}
+
+// ---------------------------------------------------------------------------
+
TEST_F(CApi, proj_obj_get_source_target_crs_conversion_without_crs) {
auto obj = proj_obj_create_from_database(
m_ctxt, "EPSG", "16031", PJ_OBJ_CATEGORY_COORDINATE_OPERATION, false,
@@ -1955,4 +1971,69 @@ TEST_F(CApi, proj_obj_create_projections) {
/* END: Generated by scripts/create_c_api_projections.py*/
}
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_obj_cs_get_axis_info) {
+ {
+ auto crs = proj_obj_create_from_database(
+ m_ctxt, "EPSG", "4326", PJ_OBJ_CATEGORY_CRS, false, nullptr);
+ ASSERT_NE(crs, nullptr);
+ ObjectKeeper keeper(crs);
+
+ auto cs = proj_obj_crs_get_coordinate_system(crs);
+ ASSERT_NE(cs, nullptr);
+ ObjectKeeper keeperCs(cs);
+
+ const char *type = proj_obj_cs_get_type(cs);
+ ASSERT_NE(type, nullptr);
+ EXPECT_EQ(std::string(type), "Ellipsoidal");
+
+ EXPECT_EQ(proj_obj_cs_get_axis_count(cs), 2);
+
+ EXPECT_FALSE(proj_obj_cs_get_axis_info(cs, -1, nullptr, nullptr,
+ nullptr, nullptr, nullptr));
+
+ EXPECT_FALSE(proj_obj_cs_get_axis_info(cs, 2, nullptr, nullptr, nullptr,
+ nullptr, nullptr));
+
+ EXPECT_TRUE(proj_obj_cs_get_axis_info(cs, 0, nullptr, nullptr, nullptr,
+ nullptr, nullptr));
+
+ const char *name = nullptr;
+ const char *abbrev = nullptr;
+ const char *direction = nullptr;
+ double unitConvFactor = 0.0;
+ const char *unitName = nullptr;
+
+ EXPECT_TRUE(proj_obj_cs_get_axis_info(cs, 0, &name, &abbrev, &direction,
+ &unitConvFactor, &unitName));
+ ASSERT_NE(name, nullptr);
+ ASSERT_NE(abbrev, nullptr);
+ ASSERT_NE(direction, nullptr);
+ ASSERT_NE(unitName, nullptr);
+ EXPECT_EQ(std::string(name), "Geodetic latitude");
+ EXPECT_EQ(std::string(abbrev), "Lat");
+ EXPECT_EQ(std::string(direction), "north");
+ EXPECT_EQ(unitConvFactor, 0.017453292519943295) << unitConvFactor;
+ EXPECT_EQ(std::string(unitName), "degree");
+ }
+
+ // Non CRS object
+ {
+ auto obj = proj_obj_create_from_database(
+ m_ctxt, "EPSG", "1170", PJ_OBJ_CATEGORY_COORDINATE_OPERATION, false,
+ nullptr);
+ ASSERT_NE(obj, nullptr);
+ ObjectKeeper keeper(obj);
+ EXPECT_EQ(proj_obj_crs_get_coordinate_system(obj), nullptr);
+
+ EXPECT_EQ(proj_obj_cs_get_type(obj), nullptr);
+
+ EXPECT_EQ(proj_obj_cs_get_axis_count(obj), -1);
+
+ EXPECT_FALSE(proj_obj_cs_get_axis_info(obj, 0, nullptr, nullptr,
+ nullptr, nullptr, nullptr));
+ }
+}
+
} // namespace
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 99f58739..a4974234 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -5686,23 +5686,33 @@ TEST(io, projparse_longlat_a_f_non_zero) {
PROJStringParser().createFromPROJString("+proj=longlat +a=2 +f=0.5");
auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
ASSERT_TRUE(crs != nullptr);
- WKTFormatterNNPtr f(WKTFormatter::create());
- f->simulCurNodeHasId();
- crs->exportToWKT(f.get());
- auto expected = "GEODCRS[\"unknown\",\n"
- " DATUM[\"unknown\",\n"
- " ELLIPSOID[\"unknown\",2,2,\n"
- " LENGTHUNIT[\"metre\",1]]],\n"
- " PRIMEM[\"Reference meridian\",0,\n"
- " ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
- " CS[ellipsoidal,2],\n"
- " AXIS[\"longitude\",east,\n"
- " ORDER[1],\n"
- " ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
- " AXIS[\"latitude\",north,\n"
- " ORDER[2],\n"
- " ANGLEUNIT[\"degree\",0.0174532925199433]]]";
- EXPECT_EQ(f->toString(), expected);
+ EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+ auto rf = crs->ellipsoid()->computedInverseFlattening();
+ EXPECT_EQ(rf, 2) << rf;
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_longlat_a_e) {
+ auto obj =
+ PROJStringParser().createFromPROJString("+proj=longlat +a=2 +e=0.5");
+ auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+ ASSERT_TRUE(crs != nullptr);
+ EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+ auto rf = crs->ellipsoid()->computedInverseFlattening();
+ EXPECT_NEAR(rf, 7.46410161513775, 1e-14) << rf;
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_longlat_a_es) {
+ auto obj =
+ PROJStringParser().createFromPROJString("+proj=longlat +a=2 +es=0.5");
+ auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+ ASSERT_TRUE(crs != nullptr);
+ EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+ auto rf = crs->ellipsoid()->computedInverseFlattening();
+ EXPECT_EQ(rf, 3.4142135623730958) << rf;
}
// ---------------------------------------------------------------------------
@@ -5711,23 +5721,31 @@ TEST(io, projparse_longlat_R) {
auto obj = PROJStringParser().createFromPROJString("+proj=longlat +R=2");
auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
ASSERT_TRUE(crs != nullptr);
- WKTFormatterNNPtr f(WKTFormatter::create());
- f->simulCurNodeHasId();
- crs->exportToWKT(f.get());
- auto expected = "GEODCRS[\"unknown\",\n"
- " DATUM[\"unknown\",\n"
- " ELLIPSOID[\"unknown\",2,0,\n"
- " LENGTHUNIT[\"metre\",1]]],\n"
- " PRIMEM[\"Reference meridian\",0,\n"
- " ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
- " CS[ellipsoidal,2],\n"
- " AXIS[\"longitude\",east,\n"
- " ORDER[1],\n"
- " ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
- " AXIS[\"latitude\",north,\n"
- " ORDER[2],\n"
- " ANGLEUNIT[\"degree\",0.0174532925199433]]]";
- EXPECT_EQ(f->toString(), expected);
+ EXPECT_TRUE(crs->ellipsoid()->isSphere());
+ EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_longlat_a) {
+ auto obj = PROJStringParser().createFromPROJString("+proj=longlat +a=2");
+ auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+ ASSERT_TRUE(crs != nullptr);
+ EXPECT_TRUE(crs->ellipsoid()->isSphere());
+ EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_longlat_a_override_ellps) {
+ auto obj = PROJStringParser().createFromPROJString(
+ "+proj=longlat +a=2 +ellps=WGS84");
+ auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+ ASSERT_TRUE(crs != nullptr);
+ EXPECT_TRUE(!crs->ellipsoid()->isSphere());
+ EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+ EXPECT_EQ(crs->ellipsoid()->computedInverseFlattening(), 298.25722356300003)
+ << crs->ellipsoid()->computedInverseFlattening();
}
// ---------------------------------------------------------------------------
@@ -7505,6 +7523,16 @@ TEST(io, projparse_init) {
EXPECT_EQ(co->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline +step +init=epsg:4326 +step +proj=longlat");
}
+
+ {
+ auto obj = PROJStringParser().createFromPROJString(
+ "init=epsg:4326 proj=longlat ellps=GRS80");
+ auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+ ASSERT_TRUE(crs != nullptr);
+ EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
+ "+proj=pipeline +step +proj=longlat +ellps=GRS80 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+ }
}
// ---------------------------------------------------------------------------
@@ -7536,10 +7564,6 @@ TEST(io, projparse_errors) {
EXPECT_THROW(PROJStringParser().createFromPROJString(
"proj=pipeline step init=epsg:4326 init=epsg:4326"),
ParsingException);
-
- EXPECT_THROW(PROJStringParser().createFromPROJString(
- "proj=pipeline step init=epsg:4326 proj=longlat"),
- ParsingException);
}
// ---------------------------------------------------------------------------
@@ -7577,9 +7601,6 @@ TEST(io, projparse_longlat_errors) {
PROJStringParser().createFromPROJString("+proj=longlat +R=invalid"),
ParsingException);
- EXPECT_THROW(PROJStringParser().createFromPROJString("+proj=longlat +a=1"),
- ParsingException);
-
EXPECT_THROW(PROJStringParser().createFromPROJString("+proj=longlat +b=1"),
ParsingException);
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 528aa6b8..818c10ec 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4477,6 +4477,38 @@ TEST(operation, geogCRS_to_geogCRS_CH1903_to_CH1903plus_context) {
// ---------------------------------------------------------------------------
+TEST(operation, geogCRS_to_geogCRS_init_IGNF_to_init_IGNF_context) {
+
+ auto dbContext = DatabaseContext::create();
+
+ auto sourceCRS_obj = PROJStringParser()
+ .attachDatabaseContext(dbContext)
+ .setUsePROJ4InitRules(true)
+ .createFromPROJString("+init=IGNF:NTFG");
+ auto sourceCRS = nn_dynamic_pointer_cast<CRS>(sourceCRS_obj);
+ ASSERT_TRUE(sourceCRS != nullptr);
+
+ auto targetCRS_obj = PROJStringParser()
+ .attachDatabaseContext(dbContext)
+ .setUsePROJ4InitRules(true)
+ .createFromPROJString("+init=IGNF:RGF93G");
+ auto targetCRS = nn_dynamic_pointer_cast<CRS>(targetCRS_obj);
+ ASSERT_TRUE(targetCRS != nullptr);
+
+ auto authFactory = AuthorityFactory::create(dbContext, std::string());
+ auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+ auto list = CoordinateOperationFactory::create()->createOperations(
+ NN_CHECK_ASSERT(sourceCRS), NN_CHECK_ASSERT(targetCRS), ctxt);
+ ASSERT_EQ(list.size(), 1);
+
+ EXPECT_EQ(list[0]->nameStr(),
+ "NOUVELLE TRIANGULATION DE LA FRANCE (NTF) vers RGF93 (ETRS89)");
+ EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+ "+proj=hgridshift +grids=ntf_r93.gsb");
+}
+
+// ---------------------------------------------------------------------------
+
TEST(operation, geogCRS_to_geogCRS_3D) {
auto geogcrs_m_obj =
@@ -5110,6 +5142,83 @@ TEST(operation, boundCRS_of_geogCRS_to_unrelated_geogCRS) {
// ---------------------------------------------------------------------------
+TEST(operation, createOperation_boundCRS_identified_by_datum) {
+ auto objSrc =
+ PROJStringParser().createFromPROJString("+proj=longlat +datum=WGS84");
+ auto src = nn_dynamic_pointer_cast<GeographicCRS>(objSrc);
+ ASSERT_TRUE(src != nullptr);
+
+ auto objDest = PROJStringParser().createFromPROJString(
+ "+proj=utm +zone=32 +a=6378249.2 +b=6356515 "
+ "+towgs84=-263.0,6.0,431.0 +no_defs");
+ 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 +proj=cart +ellps=WGS84 "
+ "+step +proj=helmert +x=263 +y=-6 +z=-431 "
+ "+step +inv +proj=cart +ellps=clrk80ign "
+ "+step +proj=utm +zone=32 +ellps=clrk80ign");
+
+ auto authFactory =
+ AuthorityFactory::create(DatabaseContext::create(), std::string());
+ auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+ auto list = CoordinateOperationFactory::create()->createOperations(
+ NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest), ctxt);
+ ASSERT_EQ(list.size(), 1U);
+ EXPECT_TRUE(list[0]->isEquivalentTo(op.get()));
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, boundCRS_of_clrk_66_geogCRS_to_nad83_geogCRS) {
+ auto objSrc = PROJStringParser().createFromPROJString(
+ "+proj=latlong +ellps=clrk66 +nadgrids=ntv1_can.dat,conus");
+ auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+ ASSERT_TRUE(src != nullptr);
+
+ auto objDest =
+ PROJStringParser().createFromPROJString("+proj=latlong +datum=NAD83");
+ 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 +proj=hgridshift +grids=ntv1_can.dat,conus "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, boundCRS_of_clrk_66_projCRS_to_nad83_geogCRS) {
+ auto objSrc = PROJStringParser().createFromPROJString(
+ "+proj=utm +zone=17 +ellps=clrk66 +nadgrids=ntv1_can.dat,conus");
+ auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+ ASSERT_TRUE(src != nullptr);
+
+ auto objDest =
+ PROJStringParser().createFromPROJString("+proj=latlong +datum=NAD83");
+ 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 +inv +proj=utm +zone=17 +ellps=clrk66 "
+ "+step +proj=hgridshift +grids=ntv1_can.dat,conus "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
TEST(operation, boundCRS_of_projCRS_to_geogCRS) {
auto utm31 = ProjectedCRS::create(
PropertyMap(), GeographicCRS::EPSG_4807,
@@ -5864,6 +5973,25 @@ TEST(operation, createOperation_on_crs_with_canonical_bound_crs) {
// ---------------------------------------------------------------------------
+TEST(operation, createOperation_fallback_to_proj4_strings) {
+ auto objDest = PROJStringParser().createFromPROJString(
+ "+proj=longlat +geoc +ellps=WGS84");
+ auto dest = nn_dynamic_pointer_cast<GeographicCRS>(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 +inv +proj=longlat +datum=WGS84 "
+ "+step +proj=longlat +geoc +ellps=WGS84 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
TEST(operation, mercator_variant_A_to_variant_B) {
auto projCRS = ProjectedCRS::create(
PropertyMap(), GeographicCRS::EPSG_4326,
generated by cgit v1.2.3 (git 2.25.1) at 2026-01-06 19:11:09 +0000