aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--include/proj/crs.hpp4
-rw-r--r--scripts/reference_exported_symbols.txt1
-rw-r--r--src/iso19111/coordinateoperation.cpp320
-rw-r--r--src/iso19111/crs.cpp52
-rw-r--r--test/unit/test_operation.cpp185
5 files changed, 443 insertions, 119 deletions
diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp
index 65459058..c07904c8 100644
--- a/include/proj/crs.hpp
+++ b/include/proj/crs.hpp
@@ -345,6 +345,10 @@ class PROJ_GCC_DLL GeographicCRS : public GeodeticCRS {
const datum::DatumEnsemblePtr &datumEnsemble,
const cs::EllipsoidalCSNNPtr &cs);
+ PROJ_DLL GeographicCRSNNPtr
+ demoteTo2D(const std::string &newName,
+ const io::DatabaseContextPtr &dbContext) const;
+
PROJ_DLL static const GeographicCRSNNPtr EPSG_4267; // NAD27
PROJ_DLL static const GeographicCRSNNPtr EPSG_4269; // NAD83
PROJ_DLL static const GeographicCRSNNPtr EPSG_4326; // WGS 84 2D
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index b70d3b9d..c9f44ead 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -151,6 +151,7 @@ osgeo::proj::crs::GeodeticCRS::velocityModel() const
osgeo::proj::crs::GeographicCRS::coordinateSystem() const
osgeo::proj::crs::GeographicCRS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::datum::GeodeticReferenceFrame> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::EllipsoidalCS> > const&)
osgeo::proj::crs::GeographicCRS::create(osgeo::proj::util::PropertyMap const&, std::shared_ptr<osgeo::proj::datum::GeodeticReferenceFrame> const&, std::shared_ptr<osgeo::proj::datum::DatumEnsemble> const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::EllipsoidalCS> > const&)
+osgeo::proj::crs::GeographicCRS::demoteTo2D(std::string const&, std::shared_ptr<osgeo::proj::io::DatabaseContext> const&) const
osgeo::proj::crs::GeographicCRS::~GeographicCRS()
osgeo::proj::crs::GeographicCRS::is2DPartOf3D(dropbox::oxygen::nn<osgeo::proj::crs::GeographicCRS const*>) const
osgeo::proj::crs::ParametricCRS::coordinateSystem() const
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index bf3b5146..d4f2903f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -58,7 +58,9 @@
// #define DEBUG
// #define DEBUG_SORT
-#if defined(DEBUG) || defined(DEBUG_SORT)
+// #define DEBUG_CONCATENATED_OPERATION
+#if defined(DEBUG) || defined(DEBUG_SORT) || \
+ defined(DEBUG_CONCATENATED_OPERATION)
#include <iostream>
#endif
@@ -121,6 +123,8 @@ static const char *BALLPARK_VERTICAL_TRANSFORMATION =
static const char *BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT =
" (ballpark vertical transformation, without ellipsoid height to vertical "
"height correction)";
+
+static const std::string AXIS_ORDER_CHANGE_2D_NAME = "axis order change (2D)";
//! @endcond
//! @cond Doxygen_Suppress
@@ -400,7 +404,7 @@ static double getAccuracy(const CoordinateOperationNNPtr &op) {
// ---------------------------------------------------------------------------
-// Returns the accuracy of a set of concantenated operations, or -1 if unknown
+// Returns the accuracy of a set of concatenated operations, or -1 if unknown
static double getAccuracy(const std::vector<CoordinateOperationNNPtr> &ops) {
double accuracy = -1.0;
for (const auto &subop : ops) {
@@ -4631,7 +4635,7 @@ ConversionNNPtr Conversion::createAxisOrderReversal(bool is3D) {
EPSG_CODE_METHOD_AXIS_ORDER_REVERSAL_3D),
{}, {});
} else {
- return create(createMapNameEPSGCode("axis order change (2D)", 15498),
+ return create(createMapNameEPSGCode(AXIS_ORDER_CHANGE_2D_NAME, 15498),
createMethodMapNameEPSGCode(
EPSG_CODE_METHOD_AXIS_ORDER_REVERSAL_2D),
{}, {});
@@ -7421,6 +7425,8 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
}
if (starts_with(tokens[i], INVERSE_OF)) {
name += tokens[i].substr(INVERSE_OF.size());
+ } else if (tokens[i] == AXIS_ORDER_CHANGE_2D_NAME) {
+ name += tokens[i];
} else {
name += INVERSE_OF + tokens[i];
}
@@ -9405,6 +9411,10 @@ ConcatenatedOperationNNPtr ConcatenatedOperation::create(
if (i >= 1) {
if (!compareStepCRS(l_sourceCRS.get(), lastTargetCRS.get())) {
#ifdef DEBUG_CONCATENATED_OPERATION
+ std::cerr << "Step " << i - 1 << ": "
+ << operationsIn[i - 1]->nameStr() << std::endl;
+ std::cerr << "Step " << i << ": " << operationsIn[i]->nameStr()
+ << std::endl;
{
auto f(io::WKTFormatter::create(
io::WKTFormatter::Convention::WKT2_2019));
@@ -9560,7 +9570,7 @@ void ConcatenatedOperation::fixStepsDirection(
!compareStepCRS(l_sourceCRS.get(), concatOpSourceCRS.get())) {
throw InvalidOperation("The source CRS of the first step of "
"concatenated operation is not the same "
- "as the source CRS of the concantenated "
+ "as the source CRS of the concatenated "
"operation itself");
}
@@ -9569,7 +9579,7 @@ void ConcatenatedOperation::fixStepsDirection(
!compareStepCRS(l_targetCRS.get(), concatOpTargetCRS.get())) {
throw InvalidOperation("The target CRS of the last step of "
"concatenated operation is not the same "
- "as the target CRS of the concantenated "
+ "as the target CRS of the concatenated "
"operation itself");
}
}
@@ -9648,7 +9658,7 @@ CoordinateOperationNNPtr ConcatenatedOperation::createComputeMetadata(
auto extent = getExtent(flattenOps, false, emptyIntersection);
if (checkExtent && emptyIntersection) {
std::string msg(
- "empty intersection of area of validity of concantenated "
+ "empty intersection of area of validity of concatenated "
"operations");
throw InvalidOperationEmptyIntersection(msg);
}
@@ -11265,12 +11275,9 @@ createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS,
geogSrc->datum()->_isEquivalentTo(
geogDst->datum().get(), util::IComparable::Criterion::EQUIVALENT);
- std::string name(isSameDatum ? NULL_GEOGRAPHIC_OFFSET
- : BALLPARK_GEOGRAPHIC_OFFSET);
- name += " from ";
- name += sourceCRS->nameStr();
- name += " to ";
- name += targetCRS->nameStr();
+ auto name = buildOpName(isSameDatum ? NULL_GEOGRAPHIC_OFFSET
+ : BALLPARK_GEOGRAPHIC_OFFSET,
+ sourceCRS, targetCRS);
const auto &sourceCRSExtent = getExtent(sourceCRS);
const auto &targetCRSExtent = getExtent(targetCRS);
@@ -11480,18 +11487,8 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
auto exportable = util::nn_make_shared<MyPROJStringExportableHorizVertical>(
horizTransform, verticalTransform, geogDst);
- bool horizTransformIsNoOp = horizTransform->sourceCRS()->_isEquivalentTo(
- horizTransform->targetCRS().get());
- if (!horizTransformIsNoOp) {
- const crs::GeographicCRS *geogSrc =
- dynamic_cast<const crs::GeographicCRS *>(
- horizTransform->sourceCRS().get());
- if (geogSrc) {
- horizTransformIsNoOp =
- geogSrc->is2DPartOf3D(NN_NO_CHECK(geogDst.get()));
- }
- }
-
+ const bool horizTransformIsNoOp =
+ starts_with(horizTransform->nameStr(), NULL_GEOGRAPHIC_OFFSET);
if (horizTransformIsNoOp) {
auto properties = util::PropertyMap();
properties.set(common::IdentifiedObject::NAME_KEY,
@@ -11541,26 +11538,34 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
const operation::CoordinateOperationNNPtr &opSrcCRSToGeogCRS,
const operation::CoordinateOperationNNPtr &verticalTransform,
const operation::CoordinateOperationNNPtr &opGeogCRStoDstCRS,
- const crs::GeographicCRSPtr &interpolationGeogCRS) {
+ const crs::GeographicCRSPtr &interpolationGeogCRS, bool checkExtent) {
auto exportable =
util::nn_make_shared<MyPROJStringExportableHorizVerticalHorizPROJBased>(
opSrcCRSToGeogCRS, verticalTransform, opGeogCRStoDstCRS,
interpolationGeogCRS);
- bool dummy = false;
- auto ops = opSrcCRSToGeogCRS->sourceCRS()->_isEquivalentTo(
- opSrcCRSToGeogCRS->targetCRS().get())
- ? std::vector<CoordinateOperationNNPtr>{verticalTransform,
- opGeogCRStoDstCRS}
- : std::vector<CoordinateOperationNNPtr>{opSrcCRSToGeogCRS,
- verticalTransform,
- opGeogCRStoDstCRS};
+ std::vector<CoordinateOperationNNPtr> ops;
+ if (!starts_with(opSrcCRSToGeogCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
+ ops.emplace_back(opSrcCRSToGeogCRS);
+ }
+ ops.emplace_back(verticalTransform);
+ if (!starts_with(opGeogCRStoDstCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
+ ops.emplace_back(opGeogCRStoDstCRS);
+ }
+
bool hasBallparkTransformation = false;
for (const auto &op : ops) {
hasBallparkTransformation |= op->hasBallparkTransformation();
}
- auto extent = getExtent(ops, true, dummy);
+ bool emptyIntersection = false;
+ auto extent = getExtent(ops, false, emptyIntersection);
+ if (checkExtent && emptyIntersection) {
+ std::string msg(
+ "empty intersection of area of validity of concatenated "
+ "operations");
+ throw InvalidOperationEmptyIntersection(msg);
+ }
auto properties = util::PropertyMap();
properties.set(common::IdentifiedObject::NAME_KEY,
computeConcatenatedName(ops));
@@ -12311,6 +12316,67 @@ bool CoordinateOperationFactory::Private::hasPerfectAccuracyResult(
// ---------------------------------------------------------------------------
//! @cond Doxygen_Suppress
+
+static std::vector<CoordinateOperationNNPtr>
+getOps(const CoordinateOperationNNPtr &op) {
+ auto concatenated = dynamic_cast<const ConcatenatedOperation *>(op.get());
+ if (concatenated)
+ return concatenated->operations();
+ return {op};
+}
+
+static bool useDifferentTransformationsForSameSourceTarget(
+ const CoordinateOperationNNPtr &opA, const CoordinateOperationNNPtr &opB) {
+ auto subOpsA = getOps(opA);
+ auto subOpsB = getOps(opB);
+ for (const auto &subOpA : subOpsA) {
+ if (!dynamic_cast<const Transformation *>(subOpA.get()))
+ continue;
+ if (subOpA->sourceCRS()->nameStr() == "unknown" ||
+ subOpA->targetCRS()->nameStr() == "unknown")
+ continue;
+ for (const auto &subOpB : subOpsB) {
+ if (!dynamic_cast<const Transformation *>(subOpB.get()))
+ continue;
+ if (subOpB->sourceCRS()->nameStr() == "unknown" ||
+ subOpB->targetCRS()->nameStr() == "unknown")
+ continue;
+
+ if (subOpA->sourceCRS()->nameStr() ==
+ subOpB->sourceCRS()->nameStr() &&
+ subOpA->targetCRS()->nameStr() ==
+ subOpB->targetCRS()->nameStr()) {
+ if (starts_with(subOpA->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
+ starts_with(subOpB->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
+ continue;
+ }
+
+ if (!subOpA->isEquivalentTo(subOpB.get())) {
+ return true;
+ }
+ } else if (subOpA->sourceCRS()->nameStr() ==
+ subOpB->targetCRS()->nameStr() &&
+ subOpA->targetCRS()->nameStr() ==
+ subOpB->sourceCRS()->nameStr()) {
+ if (starts_with(subOpA->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
+ starts_with(subOpB->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
+ continue;
+ }
+
+ if (!subOpA->isEquivalentTo(subOpB->inverse().get())) {
+ return true;
+ }
+ }
+ }
+ }
+ return false;
+}
+
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
std::vector<CoordinateOperationNNPtr>
CoordinateOperationFactory::Private::createOperations(
const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
@@ -13129,71 +13195,132 @@ CoordinateOperationFactory::Private::createOperations(
}
}
}
- if (!horizTransforms.empty() && !verticalTransforms.empty()) {
- for (const auto &horizTransform : horizTransforms) {
- for (const auto &verticalTransform : verticalTransforms) {
-
- crs::GeographicCRSPtr interpolationGeogCRS;
- auto transformationVerticalTransform =
- dynamic_cast<const Transformation *>(
- verticalTransform.get());
- if (transformationVerticalTransform) {
- auto interpTransformCRS =
- transformationVerticalTransform
- ->interpolationCRS();
- if (interpTransformCRS) {
- auto nn_interpTransformCRS =
- NN_NO_CHECK(interpTransformCRS);
- if (dynamic_cast<const crs::GeographicCRS *>(
- nn_interpTransformCRS.get())) {
- interpolationGeogCRS =
- util::nn_dynamic_pointer_cast<
- crs::GeographicCRS>(
- nn_interpTransformCRS);
- }
- }
+
+ if (horizTransforms.empty() || verticalTransforms.empty()) {
+ return horizTransforms;
+ }
+
+ typedef std::pair<std::vector<CoordinateOperationNNPtr>,
+ std::vector<CoordinateOperationNNPtr>>
+ PairOfTransforms;
+ std::map<std::string, PairOfTransforms>
+ cacheHorizToInterpAndInterpToTarget;
+
+ for (const auto &verticalTransform : verticalTransforms) {
+ crs::GeographicCRSPtr interpolationGeogCRS;
+ auto transformationVerticalTransform =
+ dynamic_cast<const Transformation *>(
+ verticalTransform.get());
+ if (transformationVerticalTransform &&
+ !transformationVerticalTransform
+ ->hasBallparkTransformation()) {
+ auto interpTransformCRS =
+ transformationVerticalTransform->interpolationCRS();
+ if (interpTransformCRS) {
+ interpolationGeogCRS =
+ std::dynamic_pointer_cast<crs::GeographicCRS>(
+ interpTransformCRS);
+ } else {
+ // If no explicit interpolation CRS, then
+ // this will be the geographic CRS of the
+ // vertical to geog transformation
+ interpolationGeogCRS =
+ std::dynamic_pointer_cast<crs::GeographicCRS>(
+ transformationVerticalTransform->targetCRS()
+ .as_nullable());
+ }
+ }
+
+ if (interpolationGeogCRS) {
+ if (interpolationGeogCRS->coordinateSystem()
+ ->axisList()
+ .size() == 3) {
+ io::DatabaseContextPtr dbContext;
+ if (authFactory) {
+ dbContext = authFactory->databaseContext();
}
- bool done = false;
- if (interpolationGeogCRS &&
- (interpolationGeogCRS->_isEquivalentTo(
- srcGeogCRS.get(),
- util::IComparable::Criterion::EQUIVALENT) ||
- interpolationGeogCRS->is2DPartOf3D(
- NN_NO_CHECK(srcGeogCRS.get())))) {
- auto srcToInterp = createOperations(
+ // We need to force the interpolation CRS, which
+ // will
+ // frequently be 3D, to 2D to avoid transformations
+ // between source CRS and interpolation CRS to have
+ // 3D terms.
+ interpolationGeogCRS =
+ interpolationGeogCRS
+ ->demoteTo2D(std::string(), dbContext)
+ .as_nullable();
+ }
+
+ std::vector<CoordinateOperationNNPtr> srcToInterpOps;
+ std::vector<CoordinateOperationNNPtr> interpToTargetOps;
+
+ std::string key;
+ const auto &ids = interpolationGeogCRS->identifiers();
+ if (!ids.empty()) {
+ key = (*ids.front()->codeSpace()) + ':' +
+ ids.front()->code();
+ }
+ if (!key.empty()) {
+ auto iter =
+ cacheHorizToInterpAndInterpToTarget.find(key);
+ if (iter == cacheHorizToInterpAndInterpToTarget.end()) {
+ srcToInterpOps = createOperations(
componentsSrc[0],
NN_NO_CHECK(interpolationGeogCRS), context);
- auto interpToCompoundHoriz = createOperations(
- NN_NO_CHECK(interpolationGeogCRS),
- componentsSrc[0], context);
- if (!srcToInterp.empty() &&
- !interpToCompoundHoriz.empty()) {
- auto op = createHorizVerticalHorizPROJBased(
- sourceCRS, componentsSrc[0],
- srcToInterp.front(), verticalTransform,
- interpToCompoundHoriz.front(),
- interpolationGeogCRS);
- done = true;
- res.emplace_back(
- ConcatenatedOperation::
- createComputeMetadata(
- {op, horizTransform},
- !allowEmptyIntersection));
- }
+ interpToTargetOps = createOperations(
+ NN_NO_CHECK(interpolationGeogCRS), targetCRS,
+ context);
+ cacheHorizToInterpAndInterpToTarget[key] =
+ PairOfTransforms(srcToInterpOps,
+ interpToTargetOps);
+ } else {
+ srcToInterpOps = iter->second.first;
+ interpToTargetOps = iter->second.second;
}
- if (!done) {
- auto op = createHorizVerticalPROJBased(
- sourceCRS, targetCRS, horizTransform,
- verticalTransform);
+ } else {
+ srcToInterpOps = createOperations(
+ componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS),
+ context);
+ interpToTargetOps =
+ createOperations(NN_NO_CHECK(interpolationGeogCRS),
+ targetCRS, context);
+ }
+
+ for (const auto &srcToInterp : srcToInterpOps) {
+ for (const auto &interpToTarget : interpToTargetOps) {
+
+ if (useDifferentTransformationsForSameSourceTarget(
+ srcToInterp, interpToTarget)) {
+ continue;
+ }
- res.emplace_back(op);
+ try {
+ auto op = createHorizVerticalHorizPROJBased(
+ sourceCRS, targetCRS, srcToInterp,
+ verticalTransform, interpToTarget,
+ interpolationGeogCRS, true);
+ res.emplace_back(op);
+ } catch (
+ const InvalidOperationEmptyIntersection &) {
+ continue;
+ }
}
}
+ } else {
+ // This case is probably only correct if
+ // verticalTransform and horizTransform are independant
+ // and in particular that verticalTransform does not
+ // involve a grid, because of the rather arbitrary order
+ // horizontal then vertical applied
+ for (const auto &horizTransform : horizTransforms) {
+ auto op = createHorizVerticalPROJBased(
+ sourceCRS, targetCRS, horizTransform,
+ verticalTransform);
+ res.emplace_back(op);
+ }
}
- return res;
- } else {
- return horizTransforms;
}
+
+ return res;
}
} else if (compoundSrc && geodDst) {
auto datum = geodDst->datum();
@@ -13291,10 +13418,15 @@ CoordinateOperationFactory::Private::createOperations(
for (const auto &opSrc : opSrcCRSToGeogCRS) {
for (const auto &opDst : opGeogCRStoDstCRS) {
- auto op = createHorizVerticalHorizPROJBased(
- sourceCRS, targetCRS, opSrc, verticalTransform,
- opDst, interpolationGeogCRS);
- res.emplace_back(op);
+ try {
+ auto op = createHorizVerticalHorizPROJBased(
+ sourceCRS, targetCRS, opSrc,
+ verticalTransform, opDst,
+ interpolationGeogCRS, true);
+ res.emplace_back(op);
+ } catch (
+ const InvalidOperationEmptyIntersection &) {
+ }
}
}
}
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 1b7078c3..e4e05166 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -2111,6 +2111,58 @@ GeographicCRSNNPtr GeographicCRS::createEPSG_4807() {
// ---------------------------------------------------------------------------
+/** \brief Return a variant of this CRS "demoted" to a 2D one, if not already
+ * the case.
+ *
+ *
+ * @param newName Name of the new CRS. If empty, nameStr() will be used.
+ * @param dbContext Database context to look for potentially already registered
+ * 2D CRS. May be nullptr.
+ * @return a new CRS demoted to 2D, or the current one if already 2D or not
+ * applicable.
+ * @since 7.0
+ */
+GeographicCRSNNPtr
+GeographicCRS::demoteTo2D(const std::string &newName,
+ const io::DatabaseContextPtr &dbContext) const {
+
+ const auto &axisList = coordinateSystem()->axisList();
+ if (axisList.size() == 3) {
+ const auto &l_identifiers = identifiers();
+ // First check if there is a Geographic 3D CRS in the database
+ // of the same name.
+ // This is the common practice in the EPSG dataset.
+ if (dbContext && l_identifiers.size() == 1) {
+ auto authFactory = io::AuthorityFactory::create(
+ NN_NO_CHECK(dbContext), *(l_identifiers[0]->codeSpace()));
+ auto res = authFactory->createObjectsFromName(
+ nameStr(),
+ {io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS}, false);
+ if (!res.empty()) {
+ const auto &firstRes = res.front();
+ auto firstResAsGeogCRS =
+ util::nn_dynamic_pointer_cast<GeographicCRS>(firstRes);
+ if (firstResAsGeogCRS &&
+ firstResAsGeogCRS->is2DPartOf3D(NN_NO_CHECK(this))) {
+ return NN_NO_CHECK(firstResAsGeogCRS);
+ }
+ }
+ }
+
+ auto cs = cs::EllipsoidalCS::create(util::PropertyMap(), axisList[0],
+ axisList[1]);
+ return GeographicCRS::create(
+ util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+ !newName.empty() ? newName : nameStr()),
+ datum(), datumEnsemble(), cs);
+ }
+
+ return NN_NO_CHECK(std::dynamic_pointer_cast<GeographicCRS>(
+ shared_from_this().as_nullable()));
+}
+
+// ---------------------------------------------------------------------------
+
//! @cond Doxygen_Suppress
void GeographicCRS::addAngularUnitConvertAndAxisSwap(
io::PROJStringFormatter *formatter) const {
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index c3ed371e..ddbb1cc9 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -6115,10 +6115,10 @@ TEST(operation, WGS84_G1762_to_compoundCRS_with_bound_vertCRS) {
src, NN_NO_CHECK(dst), ctxt);
ASSERT_GE(list.size(), 53U);
EXPECT_EQ(list[0]->nameStr(),
- "Inverse of unknown to WGS84 ellipsoidal height + "
"Inverse of WGS 84 to WGS 84 (G1762) + "
+ "Inverse of unknown to WGS84 ellipsoidal height + "
"Inverse of NAD83 to WGS 84 (1) + "
- "Inverse of axis order change (2D)");
+ "axis order change (2D)");
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
@@ -6465,15 +6465,19 @@ TEST(operation, compoundCRS_with_boundGeogCRS_and_boundVerticalCRS_to_geogCRS) {
// Not completely sure the order of horizontal and vertical operations
// makes sense
EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
- "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
- "+proj=unitconvert +xy_in=grad +xy_out=rad +step +inv "
- "+proj=longlat +ellps=clrk80ign +pm=paris +step +proj=push +v_3 "
- "+step +proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 "
- "+z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step "
- "+inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
- "+proj=vgridshift +grids=egm08_25.gtx +multiplier=1 +step "
- "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
- "+order=2,1");
+ "+proj=pipeline "
+ "+step +proj=axisswap +order=2,1 "
+ "+step +proj=unitconvert +xy_in=grad +xy_out=rad "
+ "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris "
+ "+step +proj=push +v_3 "
+ "+step +proj=cart +ellps=clrk80ign "
+ "+step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
+ "+convention=position_vector "
+ "+step +inv +proj=cart +ellps=WGS84 "
+ "+step +proj=pop +v_3 "
+ "+step +proj=vgridshift +grids=egm08_25.gtx +multiplier=1 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+ "+step +proj=axisswap +order=2,1");
auto grids = op->gridsNeeded(DatabaseContext::create());
EXPECT_EQ(grids.size(), 1U);
@@ -6503,13 +6507,17 @@ TEST(operation, compoundCRS_with_boundProjCRS_and_boundVerticalCRS_to_geogCRS) {
// Not completely sure the order of horizontal and vertical operations
// makes sense
EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
- "+proj=pipeline +step +inv +proj=utm +zone=31 +ellps=clrk80ign "
- "+pm=paris +step +proj=push +v_3 +step +proj=cart "
- "+ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
- "+rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
- "+ellps=WGS84 +step +proj=pop +v_3 +step +proj=vgridshift "
- "+grids=egm08_25.gtx +multiplier=1 +step +proj=unitconvert "
- "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
+ "+proj=pipeline "
+ "+step +inv +proj=utm +zone=31 +ellps=clrk80ign +pm=paris "
+ "+step +proj=push +v_3 "
+ "+step +proj=cart +ellps=clrk80ign "
+ "+step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
+ "+convention=position_vector "
+ "+step +inv +proj=cart +ellps=WGS84 "
+ "+step +proj=pop +v_3 "
+ "+step +proj=vgridshift +grids=egm08_25.gtx +multiplier=1 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+ "+step +proj=axisswap +order=2,1");
auto opInverse = CoordinateOperationFactory::create()->createOperation(
GeographicCRS::EPSG_4979, compound);
@@ -6890,8 +6898,9 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
"+proj=axisswap +order=2,1");
}
- // CompoundCRS to Geog3DCRS, with same vertical unit, but without
- // ellipsoid height <--> vertical height correction
+ // CompoundCRS to Geog3DCRS, with same vertical unit, and with
+ // direct ellipsoid height <--> vertical height correction and
+ // direct horizontal transform (no-op here)
{
auto ctxt =
CoordinateOperationContext::create(authFactory, nullptr, 0.0);
@@ -6905,16 +6914,142 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
ctxt);
ASSERT_GE(list.size(), 1U);
EXPECT_EQ(list[0]->nameStr(),
- "NAD83(NSRS2007) to WGS 84 (1) + "
- "Inverse of NAD83(2011) to NAVD88 height (1)");
+ "Inverse of NAD83(NSRS2007) to NAVD88 height (1) + "
+ "NAD83(NSRS2007) to WGS 84 (1)");
EXPECT_EQ(list[0]->exportToPROJString(
PROJStringFormatter::create(
PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
- "+proj=pipeline +step +proj=axisswap +order=2,1 "
+ "+proj=pipeline "
+ "+step +proj=axisswap +order=2,1 "
+ "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+ // Inv here since the grid is not known
+ "+step +inv +proj=vgridshift +grids=geoid09_conus.bin "
+ "+multiplier=1 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+ "+step +proj=axisswap +order=2,1");
+ }
+
+ // CompoundCRS to Geog3DCRS, with same vertical unit, and with
+ // ellipsoid height <--> vertical height correction that requires a
+ // horizontal adjustment before and after (which is empty in practice here
+ // as NAD83 to NAD83(2011) is no-op)
+ {
+ auto ctxt =
+ CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+ ctxt->setSpatialCriterion(
+ CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+ ctxt->setGridAvailabilityUse(
+ CoordinateOperationContext::GridAvailabilityUse::
+ IGNORE_GRID_AVAILABILITY);
+ // NAD83 + NAVD88 height
+ auto srcObj = createFromUserInput(
+ "EPSG:4269+5703", authFactory->databaseContext(), false);
+ auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+ ASSERT_TRUE(src != nullptr);
+ auto nnSrc = NN_NO_CHECK(src);
+
+ auto list = CoordinateOperationFactory::create()->createOperations(
+ nnSrc,
+ authFactory->createCoordinateReferenceSystem("4979"), // WGS 84
+ ctxt);
+ ASSERT_GE(list.size(), 2U);
+
+ EXPECT_EQ(list[0]->nameStr(),
+ "NAD83 to NAD83(2011) (1) + "
+ "Inverse of NAD83(2011) to NAVD88 height (1) + "
+ "Inverse of NAD83 to NAD83(2011) (1) + "
+ "NAD83 to WGS 84 (1)");
+ EXPECT_EQ(list[0]->exportToPROJString(
+ PROJStringFormatter::create(
+ PROJStringFormatter::Convention::PROJ_5,
+ authFactory->databaseContext())
+ .get()),
+ "+proj=pipeline "
+ "+step +proj=axisswap +order=2,1 "
+ "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+ "+step +proj=vgridshift +grids=g2012bu0.gtx "
+ "+multiplier=1 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+ "+step +proj=axisswap +order=2,1");
+
+ // Shows vertical step, and then horizontal step
+ EXPECT_EQ(list[1]->nameStr(),
+ "NAD83 to NAD83(2011) (1) + "
+ "Inverse of NAD83(2011) to NAVD88 height (1) + "
+ "Inverse of NAD83 to NAD83(2011) (1) + "
+ "NAD83 to WGS 84 (18)");
+ EXPECT_EQ(list[1]->exportToPROJString(
+ PROJStringFormatter::create(
+ PROJStringFormatter::Convention::PROJ_5,
+ authFactory->databaseContext())
+ .get()),
+ "+proj=pipeline "
+ "+step +proj=axisswap +order=2,1 "
+ "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+ "+step +proj=vgridshift +grids=g2012bu0.gtx "
+ "+multiplier=1 "
+ "+step +proj=hgridshift +grids=FL "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+ "+step +proj=axisswap +order=2,1");
+ }
+
+ // Another variation, but post horizontal adjustment is in two steps
+ {
+ auto ctxt =
+ CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+ ctxt->setSpatialCriterion(
+ CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+ ctxt->setGridAvailabilityUse(
+ CoordinateOperationContext::GridAvailabilityUse::
+ IGNORE_GRID_AVAILABILITY);
+ // NAD83(2011) + NAVD88 height
+ auto srcObj = createFromUserInput(
+ "EPSG:6318+5703", authFactory->databaseContext(), false);
+ auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+ ASSERT_TRUE(src != nullptr);
+ auto nnSrc = NN_NO_CHECK(src);
+
+ auto list = CoordinateOperationFactory::create()->createOperations(
+ nnSrc,
+ authFactory->createCoordinateReferenceSystem("4979"), // WGS 84
+ ctxt);
+ ASSERT_GE(list.size(), 2U);
+
+ EXPECT_EQ(list[0]->nameStr(),
+ "Inverse of NAD83(2011) to NAVD88 height (1) + "
+ "Inverse of NAD83 to NAD83(2011) (1) + "
+ "NAD83 to WGS 84 (1)");
+ EXPECT_EQ(list[0]->exportToPROJString(
+ PROJStringFormatter::create(
+ PROJStringFormatter::Convention::PROJ_5,
+ authFactory->databaseContext())
+ .get()),
+ "+proj=pipeline "
+ "+step +proj=axisswap +order=2,1 "
+ "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+ "+step +proj=vgridshift +grids=g2012bu0.gtx "
+ "+multiplier=1 "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+ "+step +proj=axisswap +order=2,1");
+
+ // Shows vertical step, and then horizontal step
+ EXPECT_EQ(list[1]->nameStr(),
+ "Inverse of NAD83(2011) to NAVD88 height (1) + "
+ "Inverse of NAD83 to NAD83(2011) (1) + "
+ "NAD83 to WGS 84 (18)");
+ EXPECT_EQ(list[1]->exportToPROJString(
+ PROJStringFormatter::create(
+ PROJStringFormatter::Convention::PROJ_5,
+ authFactory->databaseContext())
+ .get()),
+ "+proj=pipeline "
+ "+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
- "+step +proj=vgridshift +grids=g2012bu0.gtx +multiplier=1 "
+ "+step +proj=vgridshift +grids=g2012bu0.gtx "
+ "+multiplier=1 "
+ "+step +proj=hgridshift +grids=FL "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}
@@ -7002,7 +7137,7 @@ TEST(operation, compoundCRS_from_WKT2_no_id_to_geogCRS_3D_context) {
"+ellps=bessel "
"+step +proj=vgridshift +grids=naptrans2008.gtx +multiplier=1 "
"+step +proj=hgridshift +grids=rdtrans2008.gsb "
- "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
+ "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-16 04:54:08 +0000