From f6232db745af1acd2473f51f82d006372c04fc55 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Wed, 20 Feb 2019 10:36:14 +0100
Subject: Add VERTCON grid name alternatives in database, and handle filename
 substitution for VERTCON method

---
 src/iso19111/coordinateoperation.cpp | 40 ++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 2128124b..90266307 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -8001,6 +8001,46 @@ TransformationNNPtr Transformation::substitutePROJAlternativeGridNames(
         }
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_VERTCON) {
+        auto fileParameter =
+            parameterValue(EPSG_NAME_PARAMETER_VERTICAL_OFFSET_FILE,
+                           EPSG_CODE_PARAMETER_VERTICAL_OFFSET_FILE);
+        if (fileParameter &&
+            fileParameter->type() == ParameterValue::Type::FILENAME) {
+
+            auto filename = fileParameter->valueFile();
+            if (databaseContext->lookForGridAlternative(
+                    filename, projFilename, projGridFormat, inverseDirection)) {
+
+                if (filename == projFilename) {
+                    assert(!inverseDirection);
+                    return self;
+                }
+
+                auto parameters = std::vector<OperationParameterNNPtr>{
+                    createOpParamNameEPSGCode(
+                        EPSG_CODE_PARAMETER_VERTICAL_OFFSET_FILE)};
+                if (inverseDirection) {
+                    return create(createPropertiesForInverse(
+                                      self.as_nullable().get(), true, false),
+                                  targetCRS(), sourceCRS(), nullptr,
+                                  createSimilarPropertiesMethod(method()),
+                                  parameters, {ParameterValue::createFilename(
+                                                  projFilename)},
+                                  coordinateOperationAccuracies())
+                        ->inverseAsTransformation();
+                } else {
+                    return create(
+                        createSimilarPropertiesTransformation(self),
+                        sourceCRS(), targetCRS(), nullptr,
+                        createSimilarPropertiesMethod(method()), parameters,
+                        {ParameterValue::createFilename(projFilename)},
+                        coordinateOperationAccuracies());
+                }
+            }
+        }
+    }
+
     return self;
 }
 // ---------------------------------------------------------------------------
-- 
cgit v1.2.3


From 3664cb546811146c588cab6db41b1ccef6fcee7a Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Wed, 20 Feb 2019 11:12:38 +0100
Subject: compoundCRS to compoundCRS: avoid emitting dummy 'Null geographic
 offset from X to X' in transformation name

---
 src/iso19111/coordinateoperation.cpp | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 90266307..ed98832f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -10748,8 +10748,13 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
             interpolationGeogCRS);
 
     bool dummy = false;
-    auto ops = std::vector<CoordinateOperationNNPtr>{
-        opSrcCRSToGeogCRS, verticalTransform, opGeogCRStoDstCRS};
+    auto ops = opSrcCRSToGeogCRS->sourceCRS()->_isEquivalentTo(
+                   opSrcCRSToGeogCRS->targetCRS().get())
+                   ? std::vector<CoordinateOperationNNPtr>{verticalTransform,
+                                                           opGeogCRStoDstCRS}
+                   : std::vector<CoordinateOperationNNPtr>{opSrcCRSToGeogCRS,
+                                                           verticalTransform,
+                                                           opGeogCRStoDstCRS};
     auto extent = getExtent(ops, true, dummy);
     auto properties = util::PropertyMap();
     properties.set(common::IdentifiedObject::NAME_KEY,
-- 
cgit v1.2.3


From 70bc293a43def169fa34ed8e97a5cb06b336f247 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Wed, 20 Feb 2019 12:18:16 +0100
Subject: Vertical CRS transformation: synthetize a vertical unit change
 transformation when needed, and also sort Null geographic offset
 transformation in last

---
 src/iso19111/coordinateoperation.cpp | 153 ++++++++++++++++++++++++++++++-----
 1 file changed, 131 insertions(+), 22 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index ed98832f..d3446460 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -106,6 +106,7 @@ constexpr double UTM_SOUTH_FALSE_NORTHING = 10000000.0;
 static const std::string INVERSE_OF = "Inverse of ";
 static const char *NULL_GEOCENTRIC_TRANSLATION = "Null geocentric translation";
 static const char *NULL_GEOGRAPHIC_OFFSET = "Null geographic offset";
+static const char *APPROXIMATE_TRANSFORMATION = " (approximate transformation)";
 //! @endcond
 
 //! @cond Doxygen_Suppress
@@ -7054,6 +7055,39 @@ TransformationNNPtr Transformation::createVerticalOffset(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Instantiate a transformation based on the Change of Vertical Unit
+ * method.
+ *
+ * This method is defined as [EPSG:1069]
+ * (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::1069)
+ *
+ * @param properties See \ref general_properties of the conversion. If the name
+ * is not provided, it is automatically set.
+ * @param sourceCRSIn Source CRS.
+ * @param targetCRSIn Target CRS.
+ * @param factor Conversion factor
+ * @param accuracies Vector of positional accuracy (might be empty).
+ * @return a new Transformation.
+ */
+TransformationNNPtr Transformation::createChangeVerticalUnit(
+    const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
+    const crs::CRSNNPtr &targetCRSIn, const common::Scale &factor,
+    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
+    return create(
+        properties, sourceCRSIn, targetCRSIn, nullptr,
+        createMethodMapNameEPSGCode(EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT),
+        VectorOfParameters{
+            createOpParamNameEPSGCode(
+                EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR),
+        },
+        VectorOfValues{
+            factor,
+        },
+        accuracies);
+}
+
+// ---------------------------------------------------------------------------
+
 static const char *getCRSQualifierStr(const crs::CRSPtr &crs) {
     auto geod = dynamic_cast<crs::GeodeticCRS *>(crs.get());
     if (geod) {
@@ -7481,6 +7515,17 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
                                  coordinateOperationAccuracies()));
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT) {
+        const double convFactor = parameterValueNumericAsSI(
+            EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR);
+        return d->registerInv(
+            shared_from_this(),
+            createChangeVerticalUnit(
+                createPropertiesForInverse(this, false, false), l_targetCRS,
+                l_sourceCRS, common::Scale(1.0 / convFactor),
+                coordinateOperationAccuracies()));
+    }
+
     return InverseTransformation::create(NN_NO_CHECK(
         util::nn_dynamic_pointer_cast<Transformation>(shared_from_this())));
 }
@@ -8835,6 +8880,31 @@ bool SingleOperation::exportToPROJStringGeneric(
                                       "conversion");
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT) {
+        double convFactor = parameterValueNumericAsSI(
+            EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR);
+        auto uom = common::UnitOfMeasure(std::string(), convFactor,
+                                         common::UnitOfMeasure::Type::LINEAR)
+                       .exportToPROJString();
+        auto reverse_uom =
+            common::UnitOfMeasure(std::string(), 1.0 / convFactor,
+                                  common::UnitOfMeasure::Type::LINEAR)
+                .exportToPROJString();
+        if (!uom.empty()) {
+            formatter->addStep("unitconvert");
+            formatter->addParam("z_in", uom);
+            formatter->addParam("z_out", "m");
+        } else if (!reverse_uom.empty()) {
+            formatter->addStep("unitconvert");
+            formatter->addParam("z_in", "m");
+            formatter->addParam("z_out", reverse_uom);
+        } else {
+            formatter->addStep("affine");
+            formatter->addParam("s33", convFactor);
+        }
+        return true;
+    }
+
     return false;
 }
 
@@ -9667,14 +9737,18 @@ struct PrecomputedOpCharacteristics {
     bool gridsAvailable_ = false;
     bool gridsKnown_ = false;
     size_t stepCount_ = 0;
+    bool isApprox_ = false;
+    bool isNullTransformation_ = false;
 
     PrecomputedOpCharacteristics() = default;
     PrecomputedOpCharacteristics(double area, double accuracy, bool hasGrids,
                                  bool gridsAvailable, bool gridsKnown,
-                                 size_t stepCount)
+                                 size_t stepCount, bool isApprox,
+                                 bool isNullTransformation)
         : area_(area), accuracy_(accuracy), hasGrids_(hasGrids),
           gridsAvailable_(gridsAvailable), gridsKnown_(gridsKnown),
-          stepCount_(stepCount) {}
+          stepCount_(stepCount), isApprox_(isApprox),
+          isNullTransformation_(isNullTransformation) {}
 };
 
 // ---------------------------------------------------------------------------
@@ -9701,6 +9775,22 @@ struct SortFunction {
         // CAUTION: the order of the comparisons is extremely important
         // to get the intended result.
 
+        if (!iterA->second.isApprox_ && iterB->second.isApprox_) {
+            return true;
+        }
+        if (iterA->second.isApprox_ && !iterB->second.isApprox_) {
+            return false;
+        }
+
+        if (!iterA->second.isNullTransformation_ &&
+            iterB->second.isNullTransformation_) {
+            return true;
+        }
+        if (iterA->second.isNullTransformation_ &&
+            !iterB->second.isNullTransformation_) {
+            return false;
+        }
+
         if (iterA->second.hasGrids_ && iterB->second.hasGrids_) {
             // Operations where grids are all available go before other
             if (iterA->second.gridsAvailable_ &&
@@ -9926,6 +10016,8 @@ struct FilterResults {
                     if (name.find(NULL_GEOGRAPHIC_OFFSET) ==
                             std::string::npos &&
                         name.find(NULL_GEOCENTRIC_TRANSLATION) ==
+                            std::string::npos &&
+                        name.find(APPROXIMATE_TRANSFORMATION) ==
                             std::string::npos) {
                         hasOpThatContainsAreaOfInterest = true;
                     }
@@ -9961,6 +10053,8 @@ struct FilterResults {
                     if (name.find(NULL_GEOGRAPHIC_OFFSET) ==
                             std::string::npos &&
                         name.find(NULL_GEOCENTRIC_TRANSLATION) ==
+                            std::string::npos &&
+                        name.find(APPROXIMATE_TRANSFORMATION) ==
                             std::string::npos) {
                         hasOpThatContainsAreaOfInterest = true;
                     }
@@ -10067,9 +10161,17 @@ struct FilterResults {
 
             const auto stepCount = getStepCount(op);
 
+            const bool isApprox =
+                op->nameStr().find(APPROXIMATE_TRANSFORMATION) !=
+                std::string::npos;
+            const bool isNullTransformation =
+                op->nameStr().find(NULL_GEOGRAPHIC_OFFSET) !=
+                    std::string::npos ||
+                op->nameStr().find(NULL_GEOCENTRIC_TRANSLATION) !=
+                    std::string::npos;
             map[op.get()] = PrecomputedOpCharacteristics(
                 area, getAccuracy(op), hasGrids, gridsAvailable, gridsKnown,
-                stepCount);
+                stepCount, isApprox, isNullTransformation);
         }
 
         // Sort !
@@ -10082,7 +10184,8 @@ struct FilterResults {
 
         // If we have more than one result, and than the last result is the
         // default "Null geographic offset" or "Null geocentric translation"
-        // operations we have synthetized, remove it as
+        // operations we have synthetized, and that at least one operation
+        // has the desired area of interest, remove it as
         // all previous results are necessarily better
         if (hasOpThatContainsAreaOfInterest && res.size() > 1) {
             const std::string &name = res.back()->nameStr();
@@ -11571,29 +11674,35 @@ CoordinateOperationFactory::Private::createOperations(
     if (vertSrc && vertDst) {
         const auto &srcDatum = vertSrc->datum();
         const auto &dstDatum = vertDst->datum();
-        if (srcDatum && dstDatum &&
-            srcDatum->_isEquivalentTo(
-                dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) {
-            const double convSrc = vertSrc->coordinateSystem()
-                                       ->axisList()[0]
-                                       ->unit()
-                                       .conversionToSI();
-            const double convDst = vertDst->coordinateSystem()
-                                       ->axisList()[0]
-                                       ->unit()
-                                       .conversionToSI();
-            if (convSrc != convDst) {
-                const double factor = convSrc / convDst;
+        const bool equivalentVDatum =
+            (srcDatum && dstDatum &&
+             srcDatum->_isEquivalentTo(
+                 dstDatum.get(), util::IComparable::Criterion::EQUIVALENT));
+
+        const double convSrc =
+            vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI();
+        const double convDst =
+            vertDst->coordinateSystem()->axisList()[0]->unit().conversionToSI();
+        if (convSrc != convDst) {
+            const double factor = convSrc / convDst;
+            auto name =
+                buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
+            if (!equivalentVDatum) {
+                name += APPROXIMATE_TRANSFORMATION;
+                auto conv = Transformation::createChangeVerticalUnit(
+                    util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                            name),
+                    sourceCRS, targetCRS, common::Scale(factor), {});
+                res.push_back(conv);
+            } else {
                 auto conv = Conversion::createChangeVerticalUnit(
-                    util::PropertyMap().set(
-                        common::IdentifiedObject::NAME_KEY,
-                        buildTransfName(sourceCRS->nameStr(),
-                                        targetCRS->nameStr())),
+                    util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                            name),
                     common::Scale(factor));
                 conv->setCRSs(sourceCRS, targetCRS, nullptr);
                 res.push_back(conv);
-                return res;
             }
+            return res;
         }
     }
 
-- 
cgit v1.2.3


From 545936e04bda914a9caf8c1e4e0eefb4c5c80397 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Wed, 20 Feb 2019 12:57:26 +0100
Subject: Operation sorting: tweak

---
 src/iso19111/coordinateoperation.cpp | 28 ++++++++++++----------------
 1 file changed, 12 insertions(+), 16 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index d3446460..d66afd0d 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -9821,6 +9821,17 @@ struct SortFunction {
             return false;
         }
 
+        if (accuracyA < 0 && accuracyB < 0) {
+            // unknown accuracy ? then prefer operations with grids, which
+            // are likely to have best practical accuracy
+            if (iterA->second.hasGrids_ && !iterB->second.hasGrids_) {
+                return true;
+            }
+            if (!iterA->second.hasGrids_ && iterB->second.hasGrids_) {
+                return false;
+            }
+        }
+
         // Operations with larger non-zero area of use go before those with
         // lower one
         const double areaA = iterA->second.area_;
@@ -9852,15 +9863,6 @@ struct SortFunction {
             if (iterA->second.hasGrids_ && !iterB->second.hasGrids_) {
                 return false;
             }
-        } else if (accuracyA < 0 && accuracyB < 0) {
-            // unknown accuracy ? then prefer operations with grids, which
-            // are likely to have best practical accuracy
-            if (iterA->second.hasGrids_ && !iterB->second.hasGrids_) {
-                return true;
-            }
-            if (!iterA->second.hasGrids_ && iterB->second.hasGrids_) {
-                return false;
-            }
         }
 
         // The less intermediate steps, the better
@@ -10136,13 +10138,7 @@ struct FilterResults {
             bool hasGrids = false;
             bool gridsAvailable = true;
             bool gridsKnown = true;
-            if (context->getAuthorityFactory() &&
-                (gridAvailabilityUse ==
-                     CoordinateOperationContext::GridAvailabilityUse::
-                         USE_FOR_SORTING ||
-                 gridAvailabilityUse ==
-                     CoordinateOperationContext::GridAvailabilityUse::
-                         IGNORE_GRID_AVAILABILITY)) {
+            if (context->getAuthorityFactory()) {
                 const auto gridsNeeded = op->gridsNeeded(
                     context->getAuthorityFactory()->databaseContext());
                 for (const auto &gridDesc : gridsNeeded) {
-- 
cgit v1.2.3


From ca8f21ecbcc404b9e9c648784216846c048a3d69 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Wed, 20 Feb 2019 15:29:38 +0100
Subject: PROJStringFormatting: change order of emission of push/pop w.r.t axis
 swap/unitconvert to avoid useless simplification rules

---
 src/iso19111/coordinateoperation.cpp | 77 ++++++++++++++++--------------------
 1 file changed, 33 insertions(+), 44 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index d66afd0d..fdcb6af8 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -8101,7 +8101,7 @@ static void ThrowExceptionNotGeodeticGeographic(const char *trfrm_name) {
 // ---------------------------------------------------------------------------
 
 static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter,
-                                       const crs::CRSNNPtr &crs,
+                                       const crs::CRSNNPtr &crs, bool addPushV3,
                                        const char *trfrm_name) {
     auto sourceCRSGeog = dynamic_cast<const crs::GeographicCRS *>(crs.get());
     if (sourceCRSGeog) {
@@ -8109,6 +8109,11 @@ static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter,
         sourceCRSGeog->_exportToPROJString(formatter);
         formatter->stopInversion();
 
+        if (addPushV3) {
+            formatter->addStep("push");
+            formatter->addParam("v_3");
+        }
+
         formatter->addStep("cart");
         sourceCRSGeog->ellipsoid()->_exportToPROJString(formatter);
     } else {
@@ -8124,7 +8129,7 @@ static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter,
 // ---------------------------------------------------------------------------
 
 static void setupPROJGeodeticTargetCRS(io::PROJStringFormatter *formatter,
-                                       const crs::CRSNNPtr &crs,
+                                       const crs::CRSNNPtr &crs, bool addPopV3,
                                        const char *trfrm_name) {
     auto targetCRSGeog = dynamic_cast<const crs::GeographicCRS *>(crs.get());
     if (targetCRSGeog) {
@@ -8132,6 +8137,11 @@ static void setupPROJGeodeticTargetCRS(io::PROJStringFormatter *formatter,
         formatter->setCurrentStepInverted(true);
         targetCRSGeog->ellipsoid()->_exportToPROJString(formatter);
 
+        if (addPopV3) {
+            formatter->addStep("pop");
+            formatter->addParam("v_3");
+        }
+
         targetCRSGeog->_exportToPROJString(formatter);
     } else {
         auto targetCRSGeod = dynamic_cast<const crs::GeodeticCRS *>(crs.get());
@@ -8223,19 +8233,19 @@ void Transformation::_exportToPROJString(
         double z =
             parameterValueNumericAsSI(EPSG_CODE_PARAMETER_Z_AXIS_TRANSLATION);
 
-        if (methodEPSGCode == EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_2D ||
-            methodEPSGCode == EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D) {
-            formatter->addStep("push");
-            formatter->addParam("v_3");
-        }
+        bool addPushPopV3 =
+            (methodEPSGCode ==
+                 EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D ||
+             methodEPSGCode ==
+                 EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_2D ||
+             methodEPSGCode == EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
+             methodEPSGCode ==
+                 EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_2D ||
+             methodEPSGCode ==
+                 EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D);
 
-        setupPROJGeodeticSourceCRS(formatter, sourceCRS(), "Helmert");
+        setupPROJGeodeticSourceCRS(formatter, sourceCRS(), addPushPopV3,
+                                   "Helmert");
 
         formatter->addStep("helmert");
         formatter->addParam("x", x);
@@ -8299,19 +8309,8 @@ void Transformation::_exportToPROJString(
             }
         }
 
-        setupPROJGeodeticTargetCRS(formatter, targetCRS(), "Helmert");
-
-        if (methodEPSGCode == EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_2D ||
-            methodEPSGCode == EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D) {
-            formatter->addStep("pop");
-            formatter->addParam("v_3");
-        }
+        setupPROJGeodeticTargetCRS(formatter, targetCRS(), addPushPopV3,
+                                   "Helmert");
 
         return;
     }
@@ -8359,15 +8358,13 @@ void Transformation::_exportToPROJString(
         double pz = parameterValueNumericAsSI(
             EPSG_CODE_PARAMETER_ORDINATE_3_EVAL_POINT);
 
-        if (methodEPSGCode ==
-                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_PV_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_CF_GEOGRAPHIC_2D) {
-            formatter->addStep("push");
-            formatter->addParam("v_3");
-        }
+        bool addPushPopV3 =
+            (methodEPSGCode ==
+                 EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_PV_GEOGRAPHIC_2D ||
+             methodEPSGCode ==
+                 EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_CF_GEOGRAPHIC_2D);
 
-        setupPROJGeodeticSourceCRS(formatter, sourceCRS(),
+        setupPROJGeodeticSourceCRS(formatter, sourceCRS(), addPushPopV3,
                                    "Molodensky-Badekas");
 
         formatter->addStep("molobadekas");
@@ -8387,17 +8384,9 @@ void Transformation::_exportToPROJString(
             formatter->addParam("convention", "coordinate_frame");
         }
 
-        setupPROJGeodeticTargetCRS(formatter, targetCRS(),
+        setupPROJGeodeticTargetCRS(formatter, targetCRS(), addPushPopV3,
                                    "Molodensky-Badekas");
 
-        if (methodEPSGCode ==
-                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_PV_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_CF_GEOGRAPHIC_2D) {
-            formatter->addStep("pop");
-            formatter->addParam("v_3");
-        }
-
         return;
     }
 
-- 
cgit v1.2.3


From 374cc258510428fa3bfb9d8ca61ad7ac83a00db1 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Wed, 20 Feb 2019 16:13:24 +0100
Subject: CompoundCRS to Geog3DCRS: in synthetised transformation, document in
 the name we are lacking an ellipsoid height to vertCRS height correction

---
 src/iso19111/coordinateoperation.cpp | 78 +++++++++++++++++++-----------------
 1 file changed, 42 insertions(+), 36 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index fdcb6af8..ede6b72e 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -106,7 +106,12 @@ constexpr double UTM_SOUTH_FALSE_NORTHING = 10000000.0;
 static const std::string INVERSE_OF = "Inverse of ";
 static const char *NULL_GEOCENTRIC_TRANSLATION = "Null geocentric translation";
 static const char *NULL_GEOGRAPHIC_OFFSET = "Null geographic offset";
+static const char *APPROXIMATE_TRANSFORMATION_PREFIX =
+    " (approximate transformation";
 static const char *APPROXIMATE_TRANSFORMATION = " (approximate transformation)";
+static const char *APPROXIMATE_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT =
+    " (approximate transformation, without ellipsoid height to vertical height "
+    "correction)";
 //! @endcond
 
 //! @cond Doxygen_Suppress
@@ -5645,7 +5650,9 @@ void Conversion::_exportToPROJString(
             common::UnitOfMeasure(std::string(), 1.0 / convFactor,
                                   common::UnitOfMeasure::Type::LINEAR)
                 .exportToPROJString();
-        if (!uom.empty()) {
+        if (uom == "m") {
+            // do nothing
+        } else if (!uom.empty()) {
             formatter->addStep("unitconvert");
             formatter->addParam("z_in", uom);
             formatter->addParam("z_out", "m");
@@ -8879,7 +8886,9 @@ bool SingleOperation::exportToPROJStringGeneric(
             common::UnitOfMeasure(std::string(), 1.0 / convFactor,
                                   common::UnitOfMeasure::Type::LINEAR)
                 .exportToPROJString();
-        if (!uom.empty()) {
+        if (uom == "m") {
+            // do nothing
+        } else if (!uom.empty()) {
             formatter->addStep("unitconvert");
             formatter->addParam("z_in", uom);
             formatter->addParam("z_out", "m");
@@ -10008,7 +10017,7 @@ struct FilterResults {
                             std::string::npos &&
                         name.find(NULL_GEOCENTRIC_TRANSLATION) ==
                             std::string::npos &&
-                        name.find(APPROXIMATE_TRANSFORMATION) ==
+                        name.find(APPROXIMATE_TRANSFORMATION_PREFIX) ==
                             std::string::npos) {
                         hasOpThatContainsAreaOfInterest = true;
                     }
@@ -10045,7 +10054,7 @@ struct FilterResults {
                             std::string::npos &&
                         name.find(NULL_GEOCENTRIC_TRANSLATION) ==
                             std::string::npos &&
-                        name.find(APPROXIMATE_TRANSFORMATION) ==
+                        name.find(APPROXIMATE_TRANSFORMATION_PREFIX) ==
                             std::string::npos) {
                         hasOpThatContainsAreaOfInterest = true;
                     }
@@ -10147,7 +10156,7 @@ struct FilterResults {
             const auto stepCount = getStepCount(op);
 
             const bool isApprox =
-                op->nameStr().find(APPROXIMATE_TRANSFORMATION) !=
+                op->nameStr().find(APPROXIMATE_TRANSFORMATION_PREFIX) !=
                 std::string::npos;
             const bool isNullTransformation =
                 op->nameStr().find(NULL_GEOGRAPHIC_OFFSET) !=
@@ -11668,27 +11677,25 @@ CoordinateOperationFactory::Private::createOperations(
             vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI();
         const double convDst =
             vertDst->coordinateSystem()->axisList()[0]->unit().conversionToSI();
-        if (convSrc != convDst) {
-            const double factor = convSrc / convDst;
-            auto name =
-                buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
-            if (!equivalentVDatum) {
-                name += APPROXIMATE_TRANSFORMATION;
-                auto conv = Transformation::createChangeVerticalUnit(
-                    util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
-                                            name),
-                    sourceCRS, targetCRS, common::Scale(factor), {});
-                res.push_back(conv);
-            } else {
-                auto conv = Conversion::createChangeVerticalUnit(
-                    util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
-                                            name),
-                    common::Scale(factor));
-                conv->setCRSs(sourceCRS, targetCRS, nullptr);
-                res.push_back(conv);
-            }
-            return res;
+
+        const double factor = convSrc / convDst;
+        auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
+        if (!equivalentVDatum) {
+            name += APPROXIMATE_TRANSFORMATION;
+            auto conv = Transformation::createChangeVerticalUnit(
+                util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                        name),
+                sourceCRS, targetCRS, common::Scale(factor), {});
+            res.push_back(conv);
+        } else if (convSrc != convDst) {
+            auto conv = Conversion::createChangeVerticalUnit(
+                util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                        name),
+                common::Scale(factor));
+            conv->setCRSs(sourceCRS, targetCRS, nullptr);
+            res.push_back(conv);
         }
+        return res;
     }
 
     // A bit odd case as we are comparing apples to oranges, but in case
@@ -11701,17 +11708,16 @@ CoordinateOperationFactory::Private::createOperations(
         if (geogAxis.size() == 3) {
             convDst = geogAxis[2]->unit().conversionToSI();
         }
-        if (convSrc != convDst) {
-            const double factor = convSrc / convDst;
-            auto conv = Conversion::createChangeVerticalUnit(
-                util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
-                                        buildTransfName(sourceCRS->nameStr(),
-                                                        targetCRS->nameStr())),
-                common::Scale(factor));
-            conv->setCRSs(sourceCRS, targetCRS, nullptr);
-            res.push_back(conv);
-            return res;
-        }
+
+        const double factor = convSrc / convDst;
+        auto conv = Transformation::createChangeVerticalUnit(
+            util::PropertyMap().set(
+                common::IdentifiedObject::NAME_KEY,
+                buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) +
+                    APPROXIMATE_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT),
+            sourceCRS, targetCRS, common::Scale(factor), {});
+        res.push_back(conv);
+        return res;
     }
 
     // reverse of previous case
-- 
cgit v1.2.3


From 94578ea8ff38f4bc6b1f6f52b80ecf7359f5dfc2 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Wed, 20 Feb 2019 18:04:35 +0100
Subject: CoordinateOperation: add a hasBallparkTransformation() method that
 can be used to know if it includes a very approximative transformation term

---
 src/iso19111/coordinateoperation.cpp | 174 ++++++++++++++++++++++-------------
 1 file changed, 108 insertions(+), 66 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index ede6b72e..65588923 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -104,14 +104,16 @@ constexpr double UTM_NORTH_FALSE_NORTHING = 0.0;
 constexpr double UTM_SOUTH_FALSE_NORTHING = 10000000.0;
 
 static const std::string INVERSE_OF = "Inverse of ";
-static const char *NULL_GEOCENTRIC_TRANSLATION = "Null geocentric translation";
-static const char *NULL_GEOGRAPHIC_OFFSET = "Null geographic offset";
-static const char *APPROXIMATE_TRANSFORMATION_PREFIX =
-    " (approximate transformation";
-static const char *APPROXIMATE_TRANSFORMATION = " (approximate transformation)";
-static const char *APPROXIMATE_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT =
-    " (approximate transformation, without ellipsoid height to vertical height "
-    "correction)";
+static const char *BALLPARK_GEOCENTRIC_TRANSLATION =
+    "Ballpark geocentric translation";
+static const char *BALLPARK_GEOGRAPHIC_OFFSET = "Ballpark geographic offset";
+static const char *BALLPARK_VERTICAL_TRANSFORMATION_PREFIX =
+    " (ballpark vertical transformation";
+static const char *BALLPARK_VERTICAL_TRANSFORMATION =
+    " (ballpark vertical transformation)";
+static const char *BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT =
+    " (ballpark vertical transformation, without ellipsoid height to vertical "
+    "height correction)";
 //! @endcond
 
 //! @cond Doxygen_Suppress
@@ -519,6 +521,7 @@ struct CoordinateOperation::Private {
     crs::CRSPtr interpolationCRS_{};
     util::optional<common::DataEpoch> sourceCoordinateEpoch_{};
     util::optional<common::DataEpoch> targetCoordinateEpoch_{};
+    bool hasBallparkTransformation_ = false;
 
     // do not set this for a ProjectedCRS.definingConversion
     struct CRSStrongRef {
@@ -729,6 +732,27 @@ bool CoordinateOperation::isPROJInstanciable(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Return whether a coordinate operation has a "ballpark"
+ * transformation,
+ * that is a very approximate one, due to lack of more accurate transformations.
+ *
+ * Typically a null geographic offset between two horizontal datum, or a
+ * null vertical offset (or limited to unit changes) between two vertical
+ * datum. Errors of several tens to one hundred meters might be expected,
+ * compared to more accurate transformations.
+ */
+bool CoordinateOperation::hasBallparkTransformation() const {
+    return d->hasBallparkTransformation_;
+}
+
+// ---------------------------------------------------------------------------
+
+void CoordinateOperation::setHasBallparkTransformation(bool b) {
+    d->hasBallparkTransformation_ = b;
+}
+
+// ---------------------------------------------------------------------------
+
 //! @cond Doxygen_Suppress
 struct OperationMethod::Private {
     util::optional<std::string> formula_{};
@@ -1529,11 +1553,12 @@ static SingleOperationNNPtr createPROJBased(
     const util::PropertyMap &properties,
     const io::IPROJStringExportableNNPtr &projExportable,
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
-    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies =
-        std::vector<metadata::PositionalAccuracyNNPtr>()) {
+    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies,
+    bool hasBallparkTransformation) {
     return util::nn_static_pointer_cast<SingleOperation>(
         PROJBasedOperation::create(properties, projExportable, false, sourceCRS,
-                                   targetCRS, accuracies));
+                                   targetCRS, accuracies,
+                                   hasBallparkTransformation));
 }
 //! @endcond
 
@@ -6149,6 +6174,11 @@ TransformationNNPtr Transformation::create(
         accuracies);
     conv->assignSelf(conv);
     conv->setProperties(properties);
+    std::string name;
+    if (properties.getStringValue(common::IdentifiedObject::NAME_KEY, name) &&
+        ci_find(name, "ballpark") != std::string::npos) {
+        conv->setHasBallparkTransformation(true);
+    }
     return conv;
 }
 
@@ -7155,10 +7185,10 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
     // Forge a name for the inverse, either from the forward name, or
     // from the source and target CRS names
     const char *opType;
-    if (starts_with(forwardName, NULL_GEOCENTRIC_TRANSLATION)) {
-        opType = NULL_GEOCENTRIC_TRANSLATION;
-    } else if (starts_with(forwardName, NULL_GEOGRAPHIC_OFFSET)) {
-        opType = NULL_GEOGRAPHIC_OFFSET;
+    if (starts_with(forwardName, BALLPARK_GEOCENTRIC_TRANSLATION)) {
+        opType = BALLPARK_GEOCENTRIC_TRANSLATION;
+    } else if (starts_with(forwardName, BALLPARK_GEOGRAPHIC_OFFSET)) {
+        opType = BALLPARK_GEOGRAPHIC_OFFSET;
     } else if (dynamic_cast<const Transformation *>(op) ||
                starts_with(forwardName, "Transformation from ")) {
         opType = "Transformation";
@@ -9188,7 +9218,9 @@ CoordinateOperationNNPtr ConcatenatedOperation::createComputeMetadata(
     }
 
     std::vector<CoordinateOperationNNPtr> flattenOps;
+    bool hasBallparkTransformation = false;
     for (const auto &subOp : operationsIn) {
+        hasBallparkTransformation |= subOp->hasBallparkTransformation();
         auto subOpConcat =
             dynamic_cast<const ConcatenatedOperation *>(subOp.get());
         if (subOpConcat) {
@@ -9228,6 +9260,7 @@ CoordinateOperationNNPtr ConcatenatedOperation::createComputeMetadata(
     }
 
     auto op = create(properties, flattenOps, accuracies);
+    op->setHasBallparkTransformation(hasBallparkTransformation);
     op->d->computedName_ = true;
     return op;
 }
@@ -10012,13 +10045,7 @@ struct FilterResults {
                 bool extentContains =
                     extent->contains(NN_NO_CHECK(areaOfInterest));
                 if (extentContains) {
-                    const auto &name = op->nameStr();
-                    if (name.find(NULL_GEOGRAPHIC_OFFSET) ==
-                            std::string::npos &&
-                        name.find(NULL_GEOCENTRIC_TRANSLATION) ==
-                            std::string::npos &&
-                        name.find(APPROXIMATE_TRANSFORMATION_PREFIX) ==
-                            std::string::npos) {
+                    if (!op->hasBallparkTransformation()) {
                         hasOpThatContainsAreaOfInterest = true;
                     }
                 }
@@ -10049,13 +10076,7 @@ struct FilterResults {
                     !targetCRSExtent ||
                     extent->contains(NN_NO_CHECK(targetCRSExtent));
                 if (extentContainsSource && extentContainsTarget) {
-                    const auto &name = op->nameStr();
-                    if (name.find(NULL_GEOGRAPHIC_OFFSET) ==
-                            std::string::npos &&
-                        name.find(NULL_GEOCENTRIC_TRANSLATION) ==
-                            std::string::npos &&
-                        name.find(APPROXIMATE_TRANSFORMATION_PREFIX) ==
-                            std::string::npos) {
+                    if (!op->hasBallparkTransformation()) {
                         hasOpThatContainsAreaOfInterest = true;
                     }
                 }
@@ -10156,12 +10177,12 @@ struct FilterResults {
             const auto stepCount = getStepCount(op);
 
             const bool isApprox =
-                op->nameStr().find(APPROXIMATE_TRANSFORMATION_PREFIX) !=
+                op->nameStr().find(BALLPARK_VERTICAL_TRANSFORMATION_PREFIX) !=
                 std::string::npos;
             const bool isNullTransformation =
-                op->nameStr().find(NULL_GEOGRAPHIC_OFFSET) !=
+                op->nameStr().find(BALLPARK_GEOGRAPHIC_OFFSET) !=
                     std::string::npos ||
-                op->nameStr().find(NULL_GEOCENTRIC_TRANSLATION) !=
+                op->nameStr().find(BALLPARK_GEOCENTRIC_TRANSLATION) !=
                     std::string::npos;
             map[op.get()] = PrecomputedOpCharacteristics(
                 area, getAccuracy(op), hasGrids, gridsAvailable, gridsKnown,
@@ -10177,14 +10198,16 @@ struct FilterResults {
     void removeSyntheticNullTransforms() {
 
         // If we have more than one result, and than the last result is the
-        // default "Null geographic offset" or "Null geocentric translation"
+        // default "Ballpark geographic offset" or "Ballpark geocentric
+        // translation"
         // operations we have synthetized, and that at least one operation
         // has the desired area of interest, remove it as
         // all previous results are necessarily better
         if (hasOpThatContainsAreaOfInterest && res.size() > 1) {
             const std::string &name = res.back()->nameStr();
-            if (name.find(NULL_GEOGRAPHIC_OFFSET) != std::string::npos ||
-                name.find(NULL_GEOCENTRIC_TRANSLATION) != std::string::npos) {
+            if (name.find(BALLPARK_GEOGRAPHIC_OFFSET) != std::string::npos ||
+                name.find(BALLPARK_GEOCENTRIC_TRANSLATION) !=
+                    std::string::npos) {
                 std::vector<CoordinateOperationNNPtr> resTemp;
                 for (size_t i = 0; i < res.size() - 1; i++) {
                     resTemp.emplace_back(res[i]);
@@ -10597,9 +10620,9 @@ static std::vector<CoordinateOperationNNPtr> findsOpsInRegistryWithIntermediate(
 
 //! @cond Doxygen_Suppress
 static TransformationNNPtr
-createNullGeographicOffset(const crs::CRSNNPtr &sourceCRS,
-                           const crs::CRSNNPtr &targetCRS) {
-    std::string name(NULL_GEOGRAPHIC_OFFSET);
+createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS,
+                               const crs::CRSNNPtr &targetCRS) {
+    std::string name(BALLPARK_GEOGRAPHIC_OFFSET);
     name += " from ";
     name += sourceCRS->nameStr();
     name += " to ";
@@ -10790,7 +10813,7 @@ createGeodToGeodPROJBased(const crs::CRSNNPtr &geodSrc,
     auto properties = util::PropertyMap().set(
         common::IdentifiedObject::NAME_KEY,
         buildTransfName(geodSrc->nameStr(), geodDst->nameStr()));
-    return createPROJBased(properties, exportable, geodSrc, geodDst);
+    return createPROJBased(properties, exportable, geodSrc, geodDst, {}, false);
 }
 
 // ---------------------------------------------------------------------------
@@ -10827,7 +10850,9 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
     }
 
     return createPROJBased(properties, exportable, sourceCRS, targetCRS,
-                           accuracies);
+                           accuracies,
+                           horizTransform->hasBallparkTransformation() ||
+                               verticalTransform->hasBallparkTransformation());
 }
 
 // ---------------------------------------------------------------------------
@@ -10852,6 +10877,10 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
                    : std::vector<CoordinateOperationNNPtr>{opSrcCRSToGeogCRS,
                                                            verticalTransform,
                                                            opGeogCRStoDstCRS};
+    bool hasBallparkTransformation = false;
+    for (const auto &op : ops) {
+        hasBallparkTransformation |= op->hasBallparkTransformation();
+    }
     auto extent = getExtent(ops, true, dummy);
     auto properties = util::PropertyMap();
     properties.set(common::IdentifiedObject::NAME_KEY,
@@ -10870,7 +10899,7 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
     }
 
     return createPROJBased(properties, exportable, sourceCRS, targetCRS,
-                           accuracies);
+                           accuracies, hasBallparkTransformation);
 }
 
 //! @endcond
@@ -10926,6 +10955,11 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
 
     std::string name(buildTransfName(geogSrc->nameStr(), geogDst->nameStr()));
 
+    const bool sameDatum =
+        geogSrc->datum() != nullptr && geogDst->datum() != nullptr &&
+        geogSrc->datum()->_isEquivalentTo(
+            geogDst->datum().get(), util::IComparable::Criterion::EQUIVALENT);
+
     // Do they differ by vertical units ?
     if (vconvSrc != vconvDst &&
         geogSrc->ellipsoid()->_isEquivalentTo(
@@ -10941,18 +10975,19 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
                                         name),
                 common::Scale(factor));
             conv->setCRSs(sourceCRS, targetCRS, nullptr);
+            conv->setHasBallparkTransformation(!sameDatum);
             res.push_back(conv);
             return res;
         } else {
-            res.emplace_back(createGeodToGeodPROJBased(sourceCRS, targetCRS));
+            auto op = createGeodToGeodPROJBased(sourceCRS, targetCRS);
+            op->setHasBallparkTransformation(!sameDatum);
+            res.emplace_back(op);
             return res;
         }
     }
 
     // Do the CRS differ only by their axis order ?
-    if (geogSrc->datum() != nullptr && geogDst->datum() != nullptr &&
-        geogSrc->datum()->_isEquivalentTo(
-            geogDst->datum().get(), util::IComparable::Criterion::EQUIVALENT) &&
+    if (sameDatum &&
         !srcCS->_isEquivalentTo(dstCS.get(),
                                 util::IComparable::Criterion::EQUIVALENT)) {
         auto srcOrder = srcCS->axisOrder();
@@ -11013,7 +11048,8 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
                          metadata::Extent::WORLD),
                 datum, dstCS));
 
-        steps.emplace_back(createNullGeographicOffset(sourceCRS, interm_crs));
+        steps.emplace_back(
+            createBallparkGeographicOffset(sourceCRS, interm_crs));
 
         steps.emplace_back(Transformation::createLongitudeRotation(
             util::PropertyMap()
@@ -11048,15 +11084,17 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
                          metadata::Extent::WORLD),
                 sourceCRS, interm_crs, offset_pm));
             steps.emplace_back(
-                createNullGeographicOffset(interm_crs, targetCRS));
+                createBallparkGeographicOffset(interm_crs, targetCRS));
         } else {
             steps.emplace_back(
-                createNullGeographicOffset(sourceCRS, targetCRS));
+                createBallparkGeographicOffset(sourceCRS, targetCRS));
         }
     }
 
-    res.emplace_back(ConcatenatedOperation::createComputeMetadata(
-        steps, !allowEmptyIntersection));
+    auto op = ConcatenatedOperation::createComputeMetadata(
+        steps, !allowEmptyIntersection);
+    op->setHasBallparkTransformation(!sameDatum);
+    res.emplace_back(op);
     return res;
 }
 
@@ -11108,8 +11146,8 @@ findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
 
 static bool isNullTransformation(const std::string &name) {
 
-    return starts_with(name, NULL_GEOCENTRIC_TRANSLATION) ||
-           starts_with(name, NULL_GEOGRAPHIC_OFFSET);
+    return starts_with(name, BALLPARK_GEOCENTRIC_TRANSLATION) ||
+           starts_with(name, BALLPARK_GEOGRAPHIC_OFFSET);
 }
 
 // ---------------------------------------------------------------------------
@@ -11227,9 +11265,9 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 // ---------------------------------------------------------------------------
 
 static CoordinateOperationNNPtr
-createNullGeocentricTranslation(const crs::CRSNNPtr &sourceCRS,
-                                const crs::CRSNNPtr &targetCRS) {
-    std::string name(NULL_GEOCENTRIC_TRANSLATION);
+createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS,
+                                    const crs::CRSNNPtr &targetCRS) {
+    std::string name(BALLPARK_GEOCENTRIC_TRANSLATION);
     name += " from ";
     name += sourceCRS->nameStr();
     name += " to ";
@@ -11451,7 +11489,7 @@ CoordinateOperationFactory::Private::createOperations(
                             util::nn_dynamic_pointer_cast<cs::CartesianCS>(
                                 geodSrc->coordinateSystem()))));
                 auto opFirst =
-                    createNullGeocentricTranslation(sourceCRS, interm_crs);
+                    createBallparkGeocentricTranslation(sourceCRS, interm_crs);
                 auto opSecond =
                     createGeographicGeocentric(interm_crs, targetCRS);
                 res.emplace_back(ConcatenatedOperation::createComputeMetadata(
@@ -11466,7 +11504,7 @@ CoordinateOperationFactory::Private::createOperations(
 
         if (isSrcGeocentric && isTargetGeocentric) {
             res.emplace_back(
-                createNullGeocentricTranslation(sourceCRS, targetCRS));
+                createBallparkGeocentricTranslation(sourceCRS, targetCRS));
             return res;
         }
 
@@ -11681,11 +11719,12 @@ CoordinateOperationFactory::Private::createOperations(
         const double factor = convSrc / convDst;
         auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
         if (!equivalentVDatum) {
-            name += APPROXIMATE_TRANSFORMATION;
+            name += BALLPARK_VERTICAL_TRANSFORMATION;
             auto conv = Transformation::createChangeVerticalUnit(
                 util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
                                         name),
                 sourceCRS, targetCRS, common::Scale(factor), {});
+            conv->setHasBallparkTransformation(true);
             res.push_back(conv);
         } else if (convSrc != convDst) {
             auto conv = Conversion::createChangeVerticalUnit(
@@ -11714,8 +11753,9 @@ CoordinateOperationFactory::Private::createOperations(
             util::PropertyMap().set(
                 common::IdentifiedObject::NAME_KEY,
                 buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) +
-                    APPROXIMATE_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT),
+                    BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT),
             sourceCRS, targetCRS, common::Scale(factor), {});
+        conv->setHasBallparkTransformation(true);
         res.push_back(conv);
         return res;
     }
@@ -12022,7 +12062,8 @@ PROJBasedOperationNNPtr PROJBasedOperation::create(
     const util::PropertyMap &properties,
     const io::IPROJStringExportableNNPtr &projExportable, bool inverse,
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
-    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
+    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies,
+    bool hasBallparkTransformation) {
 
     auto formatter = io::PROJStringFormatter::create();
     if (inverse) {
@@ -12048,6 +12089,7 @@ PROJBasedOperationNNPtr PROJBasedOperation::create(
     op->setAccuracies(accuracies);
     op->projStringExportable_ = projExportable.as_nullable();
     op->inverse_ = inverse;
+    op->setHasBallparkTransformation(hasBallparkTransformation);
     return op;
 }
 
@@ -12061,7 +12103,7 @@ CoordinateOperationNNPtr PROJBasedOperation::inverse() const {
                 createPropertiesForInverse(this, false, false),
                 NN_NO_CHECK(projStringExportable_), !inverse_,
                 NN_NO_CHECK(targetCRS()), NN_NO_CHECK(sourceCRS()),
-                coordinateOperationAccuracies()));
+                coordinateOperationAccuracies(), hasBallparkTransformation()));
     }
 
     auto formatter = io::PROJStringFormatter::create();
@@ -12074,11 +12116,11 @@ CoordinateOperationNNPtr PROJBasedOperation::inverse() const {
     }
     formatter->stopInversion();
 
-    return util::nn_static_pointer_cast<CoordinateOperation>(
-        PROJBasedOperation::create(
-            createPropertiesForInverse(this, false, false),
-            formatter->toString(), targetCRS(), sourceCRS(),
-            coordinateOperationAccuracies()));
+    auto op = PROJBasedOperation::create(
+        createPropertiesForInverse(this, false, false), formatter->toString(),
+        targetCRS(), sourceCRS(), coordinateOperationAccuracies());
+    op->setHasBallparkTransformation(hasBallparkTransformation());
+    return util::nn_static_pointer_cast<CoordinateOperation>(op);
 }
 
 // ---------------------------------------------------------------------------
-- 
cgit v1.2.3


From 69ef7449f5f26453a8b6cab1ba02cb870055615f Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Wed, 20 Feb 2019 20:42:26 +0100
Subject: typo fixes: s/Explictly/Explicitly/ and s/instanciat/instantiat/

---
 src/iso19111/coordinateoperation.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 65588923..224b19ef 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -715,7 +715,7 @@ void CoordinateOperation::setAccuracies(
  * a PROJ pipeline, checking in particular that referenced grids are
  * available.
  */
-bool CoordinateOperation::isPROJInstanciable(
+bool CoordinateOperation::isPROJInstantiable(
     const io::DatabaseContextPtr &databaseContext) const {
     try {
         exportToPROJString(io::PROJStringFormatter::create().get());
-- 
cgit v1.2.3


From 34168532d3f3773e1afa2e560bdfba9bd369a0aa Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Thu, 21 Feb 2019 12:47:38 +0100
Subject: Geog2D+Height -> Geog3D of same datum: avoid inserting a useless
 ballpark horizontal transformation

---
 src/iso19111/coordinateoperation.cpp | 70 +++++++++++++++++++++++++-----------
 1 file changed, 49 insertions(+), 21 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 224b19ef..6e2d7f4f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -10829,30 +10829,58 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
     auto exportable = util::nn_make_shared<MyPROJStringExportableHorizVertical>(
         horizTransform, verticalTransform, geogDst);
 
-    bool dummy = false;
-    auto ops = std::vector<CoordinateOperationNNPtr>{horizTransform,
-                                                     verticalTransform};
-    auto extent = getExtent(ops, true, dummy);
-    auto properties = util::PropertyMap();
-    properties.set(common::IdentifiedObject::NAME_KEY,
-                   computeConcatenatedName(ops));
-
-    if (extent) {
-        properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
-                       NN_NO_CHECK(extent));
+    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()));
+        }
     }
 
-    std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
-    const double accuracy = getAccuracy(ops);
-    if (accuracy >= 0.0) {
-        accuracies.emplace_back(
-            metadata::PositionalAccuracy::create(toString(accuracy)));
-    }
+    if (horizTransformIsNoOp) {
+        auto properties = util::PropertyMap();
+        properties.set(common::IdentifiedObject::NAME_KEY,
+                       verticalTransform->nameStr());
+        bool dummy = false;
+        auto extent = getExtent(verticalTransform, true, dummy);
+        if (extent) {
+            properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+                           NN_NO_CHECK(extent));
+        }
+        return createPROJBased(
+            properties, exportable, sourceCRS, targetCRS,
+            verticalTransform->coordinateOperationAccuracies(),
+            verticalTransform->hasBallparkTransformation());
+    } else {
+        bool dummy = false;
+        auto ops = std::vector<CoordinateOperationNNPtr>{horizTransform,
+                                                         verticalTransform};
+        auto extent = getExtent(ops, true, dummy);
+        auto properties = util::PropertyMap();
+        properties.set(common::IdentifiedObject::NAME_KEY,
+                       computeConcatenatedName(ops));
 
-    return createPROJBased(properties, exportable, sourceCRS, targetCRS,
-                           accuracies,
-                           horizTransform->hasBallparkTransformation() ||
-                               verticalTransform->hasBallparkTransformation());
+        if (extent) {
+            properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+                           NN_NO_CHECK(extent));
+        }
+
+        std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
+        const double accuracy = getAccuracy(ops);
+        if (accuracy >= 0.0) {
+            accuracies.emplace_back(
+                metadata::PositionalAccuracy::create(toString(accuracy)));
+        }
+
+        return createPROJBased(
+            properties, exportable, sourceCRS, targetCRS, accuracies,
+            horizTransform->hasBallparkTransformation() ||
+                verticalTransform->hasBallparkTransformation());
+    }
 }
 
 // ---------------------------------------------------------------------------
-- 
cgit v1.2.3


From 287230f86d89a26574c777bb5e5b498084a84897 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Thu, 21 Feb 2019 13:48:28 +0100
Subject: Transformation: reintroduce the term of 'Null geographic offset' for
 transformations between geographic CRS of same datum (typically 3D to 2D)

---
 src/iso19111/coordinateoperation.cpp | 33 +++++++++++++++++++++++++++++----
 1 file changed, 29 insertions(+), 4 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 6e2d7f4f..7b0adc6f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -106,6 +106,7 @@ constexpr double UTM_SOUTH_FALSE_NORTHING = 10000000.0;
 static const std::string INVERSE_OF = "Inverse of ";
 static const char *BALLPARK_GEOCENTRIC_TRANSLATION =
     "Ballpark geocentric translation";
+static const char *NULL_GEOGRAPHIC_OFFSET = "Null geographic offset";
 static const char *BALLPARK_GEOGRAPHIC_OFFSET = "Ballpark geographic offset";
 static const char *BALLPARK_VERTICAL_TRANSFORMATION_PREFIX =
     " (ballpark vertical transformation";
@@ -7189,6 +7190,8 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
         opType = BALLPARK_GEOCENTRIC_TRANSLATION;
     } else if (starts_with(forwardName, BALLPARK_GEOGRAPHIC_OFFSET)) {
         opType = BALLPARK_GEOGRAPHIC_OFFSET;
+    } else if (starts_with(forwardName, NULL_GEOGRAPHIC_OFFSET)) {
+        opType = NULL_GEOGRAPHIC_OFFSET;
     } else if (dynamic_cast<const Transformation *>(op) ||
                starts_with(forwardName, "Transformation from ")) {
         opType = "Transformation";
@@ -10182,6 +10185,8 @@ struct FilterResults {
             const bool isNullTransformation =
                 op->nameStr().find(BALLPARK_GEOGRAPHIC_OFFSET) !=
                     std::string::npos ||
+                op->nameStr().find(NULL_GEOGRAPHIC_OFFSET) !=
+                    std::string::npos ||
                 op->nameStr().find(BALLPARK_GEOCENTRIC_TRANSLATION) !=
                     std::string::npos;
             map[op.get()] = PrecomputedOpCharacteristics(
@@ -10206,6 +10211,7 @@ struct FilterResults {
         if (hasOpThatContainsAreaOfInterest && res.size() > 1) {
             const std::string &name = res.back()->nameStr();
             if (name.find(BALLPARK_GEOGRAPHIC_OFFSET) != std::string::npos ||
+                name.find(NULL_GEOGRAPHIC_OFFSET) != std::string::npos ||
                 name.find(BALLPARK_GEOCENTRIC_TRANSLATION) !=
                     std::string::npos) {
                 std::vector<CoordinateOperationNNPtr> resTemp;
@@ -10622,7 +10628,18 @@ static std::vector<CoordinateOperationNNPtr> findsOpsInRegistryWithIntermediate(
 static TransformationNNPtr
 createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS,
                                const crs::CRSNNPtr &targetCRS) {
-    std::string name(BALLPARK_GEOGRAPHIC_OFFSET);
+
+    const crs::GeographicCRS *geogSrc =
+        dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
+    const crs::GeographicCRS *geogDst =
+        dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
+    const bool isSameDatum =
+        geogSrc && geogDst && geogSrc->datum() && geogDst->datum() &&
+        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 ";
@@ -10641,6 +10658,12 @@ createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS,
              sameExtent ? NN_NO_CHECK(sourceCRSExtent)
                         : metadata::Extent::WORLD);
     const common::Angle angle0(0);
+
+    std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
+    if (isSameDatum) {
+        accuracies.emplace_back(metadata::PositionalAccuracy::create("0"));
+    }
+
     if (dynamic_cast<const crs::SingleCRS *>(sourceCRS.get())
                 ->coordinateSystem()
                 ->axisList()
@@ -10650,10 +10673,11 @@ createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS,
                 ->axisList()
                 .size() == 3) {
         return Transformation::createGeographic3DOffsets(
-            map, sourceCRS, targetCRS, angle0, angle0, common::Length(0), {});
+            map, sourceCRS, targetCRS, angle0, angle0, common::Length(0),
+            accuracies);
     } else {
         return Transformation::createGeographic2DOffsets(
-            map, sourceCRS, targetCRS, angle0, angle0, {});
+            map, sourceCRS, targetCRS, angle0, angle0, accuracies);
     }
 }
 //! @endcond
@@ -11175,7 +11199,8 @@ findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
 static bool isNullTransformation(const std::string &name) {
 
     return starts_with(name, BALLPARK_GEOCENTRIC_TRANSLATION) ||
-           starts_with(name, BALLPARK_GEOGRAPHIC_OFFSET);
+           starts_with(name, BALLPARK_GEOGRAPHIC_OFFSET) ||
+           starts_with(name, NULL_GEOGRAPHIC_OFFSET);
 }
 
 // ---------------------------------------------------------------------------
-- 
cgit v1.2.3


From 22a786693f66c815220861f15fd041584d32a1f1 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sun, 24 Feb 2019 12:06:58 +0100
Subject: ParameterValue::_exportToWKT(): fix null pointer dereference with
 -D_GLIBCXX_ASSERTIONS on WKT with PARAMETERFILE (fixes #1290)

---
 src/iso19111/coordinateoperation.cpp | 85 +++++++++++++++++++-----------------
 1 file changed, 44 insertions(+), 41 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 7b0adc6f..fbb67e6b 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -2156,52 +2156,55 @@ void ParameterValue::_exportToWKT(io::WKTFormatter *formatter) const {
     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
 
     const auto &l_type = type();
-    const auto &l_value = value();
-    if (formatter->abridgedTransformation() && l_type == Type::MEASURE) {
-        const auto &unit = l_value.unit();
-        const auto &unitType = unit.type();
-        if (unitType == common::UnitOfMeasure::Type::LINEAR) {
-            formatter->add(l_value.getSIValue());
-        } else if (unitType == common::UnitOfMeasure::Type::ANGULAR) {
-            formatter->add(
-                l_value.convertToUnit(common::UnitOfMeasure::ARC_SECOND));
-        } else if (unit == common::UnitOfMeasure::PARTS_PER_MILLION) {
-            formatter->add(1.0 + l_value.value() * 1e-6);
-        } else {
-            formatter->add(l_value.value());
-        }
-    } else if (l_type == Type::MEASURE) {
-        const auto &unit = l_value.unit();
-        if (isWKT2) {
-            formatter->add(l_value.value());
-        } else {
-            // In WKT1, as we don't output the natural unit, output to the
-            // registered linear / angular unit.
+    if (l_type == Type::MEASURE) {
+        const auto &l_value = value();
+        if (formatter->abridgedTransformation()) {
+            const auto &unit = l_value.unit();
             const auto &unitType = unit.type();
             if (unitType == common::UnitOfMeasure::Type::LINEAR) {
-                const auto &targetUnit = *(formatter->axisLinearUnit());
-                if (targetUnit.conversionToSI() == 0.0) {
-                    throw io::FormattingException(
-                        "cannot convert value to target linear unit");
-                }
-                formatter->add(l_value.convertToUnit(targetUnit));
+                formatter->add(l_value.getSIValue());
             } else if (unitType == common::UnitOfMeasure::Type::ANGULAR) {
-                const auto &targetUnit = *(formatter->axisAngularUnit());
-                if (targetUnit.conversionToSI() == 0.0) {
-                    throw io::FormattingException(
-                        "cannot convert value to target angular unit");
-                }
-                formatter->add(l_value.convertToUnit(targetUnit));
+                formatter->add(
+                    l_value.convertToUnit(common::UnitOfMeasure::ARC_SECOND));
+            } else if (unit == common::UnitOfMeasure::PARTS_PER_MILLION) {
+                formatter->add(1.0 + l_value.value() * 1e-6);
             } else {
-                formatter->add(l_value.getSIValue());
+                formatter->add(l_value.value());
             }
-        }
-        if (isWKT2 && unit != common::UnitOfMeasure::NONE) {
-            if (!formatter->primeMeridianOrParameterUnitOmittedIfSameAsAxis() ||
-                (unit != common::UnitOfMeasure::SCALE_UNITY &&
-                 unit != *(formatter->axisLinearUnit()) &&
-                 unit != *(formatter->axisAngularUnit()))) {
-                unit._exportToWKT(formatter);
+        } else {
+            const auto &unit = l_value.unit();
+            if (isWKT2) {
+                formatter->add(l_value.value());
+            } else {
+                // In WKT1, as we don't output the natural unit, output to the
+                // registered linear / angular unit.
+                const auto &unitType = unit.type();
+                if (unitType == common::UnitOfMeasure::Type::LINEAR) {
+                    const auto &targetUnit = *(formatter->axisLinearUnit());
+                    if (targetUnit.conversionToSI() == 0.0) {
+                        throw io::FormattingException(
+                            "cannot convert value to target linear unit");
+                    }
+                    formatter->add(l_value.convertToUnit(targetUnit));
+                } else if (unitType == common::UnitOfMeasure::Type::ANGULAR) {
+                    const auto &targetUnit = *(formatter->axisAngularUnit());
+                    if (targetUnit.conversionToSI() == 0.0) {
+                        throw io::FormattingException(
+                            "cannot convert value to target angular unit");
+                    }
+                    formatter->add(l_value.convertToUnit(targetUnit));
+                } else {
+                    formatter->add(l_value.getSIValue());
+                }
+            }
+            if (isWKT2 && unit != common::UnitOfMeasure::NONE) {
+                if (!formatter
+                         ->primeMeridianOrParameterUnitOmittedIfSameAsAxis() ||
+                    (unit != common::UnitOfMeasure::SCALE_UNITY &&
+                     unit != *(formatter->axisLinearUnit()) &&
+                     unit != *(formatter->axisAngularUnit()))) {
+                    unit._exportToWKT(formatter);
+                }
             }
         }
     } else if (l_type == Type::STRING || l_type == Type::FILENAME) {
-- 
cgit v1.2.3


From 1a2513badd1406ef061fd044273268a1e0ab3eed Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Fri, 1 Mar 2019 14:53:22 +0100
Subject: Doc: rename to ISO-19111:2019

And publish link to corresponding promoted and public OGC doc:
http://docs.opengeospatial.org/as/18-005r4/18-005r4.html
---
 src/iso19111/coordinateoperation.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index fbb67e6b..d5bcfa84 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -1,7 +1,7 @@
 /******************************************************************************
  *
  * Project:  PROJ
- * Purpose:  ISO19111:2018 implementation
+ * Purpose:  ISO19111:2019 implementation
  * Author:   Even Rouault <even dot rouault at spatialys dot com>
  *
  ******************************************************************************
-- 
cgit v1.2.3


From a1cc9977decb62b4576c6c7f17a0d64cab9bf36a Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sat, 16 Mar 2019 11:03:07 +0100
Subject: Fix doc generation with Breathe 4.12.0

Breathe 4.12.0 (as pulled by MacOSX builds such as
https://travis-ci.com/OSGeo/proj.4/jobs/185395222)
does not seem to like default initialization in documented C++ structs (regression/bug)

/Users/travis/build/OSGeo/proj.4/docs/source/development/reference/cpp/io.rst:6:Parsing of expression failed. Using fallback parser. Error was:
Error in postfix expression, expected primary expression or type.
If primary expression:
  Invalid definition: Expected identifier in nested name. [error at 67]
    std::string osgeo::proj::io::AuthorityFactory::CRSInfo::authName = {}
    -------------------------------------------------------------------^
If type:
  Invalid definition: Expected identifier in nested name. [error at 67]
    std::string osgeo::proj::io::AuthorityFactory::CRSInfo::authName = {}
    -------------------------------------------------------------------^
---
 src/iso19111/coordinateoperation.cpp | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index d5bcfa84..0eef9954 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -552,7 +552,16 @@ struct CoordinateOperation::Private {
 
 // ---------------------------------------------------------------------------
 
-GridDescription::GridDescription() = default;
+GridDescription::GridDescription():
+    shortName{},
+    fullName{},
+    packageName{},
+    url{},
+    directDownload(false),
+    openLicense(false),
+    available(false)
+{}
+
 
 GridDescription::~GridDescription() = default;
 
-- 
cgit v1.2.3


From a9886a5e4a64610aa605e98afeb467d1c44ed925 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sat, 16 Mar 2019 13:31:09 +0100
Subject: Run scripts/reformat_cpp.sh

---
 src/iso19111/coordinateoperation.cpp | 15 ++++-----------
 1 file changed, 4 insertions(+), 11 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 0eef9954..575eae39 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -552,16 +552,9 @@ struct CoordinateOperation::Private {
 
 // ---------------------------------------------------------------------------
 
-GridDescription::GridDescription():
-    shortName{},
-    fullName{},
-    packageName{},
-    url{},
-    directDownload(false),
-    openLicense(false),
-    available(false)
-{}
-
+GridDescription::GridDescription()
+    : shortName{}, fullName{}, packageName{}, url{}, directDownload(false),
+      openLicense(false), available(false) {}
 
 GridDescription::~GridDescription() = default;
 
@@ -10811,7 +10804,7 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
 
 MyPROJStringExportableHorizVerticalHorizPROJBased::
     ~MyPROJStringExportableHorizVerticalHorizPROJBased() = default;
-}
+} // namespace operation
 NS_PROJ_END
 
 #if 0
-- 
cgit v1.2.3


From d95b0cb8b317d0a8142b664c42928083145194ed Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sat, 16 Mar 2019 15:14:55 +0100
Subject: createOperation(): fix geocent <--> nadgrids+geoidgrids case (fixes
 #1323)

---
 src/iso19111/coordinateoperation.cpp | 31 ++++++++++++++++++++++++++++++-
 1 file changed, 30 insertions(+), 1 deletion(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 575eae39..240833df 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -11914,11 +11914,40 @@ CoordinateOperationFactory::Private::createOperations(
                 return horizTransforms;
             }
         }
+    } else if (compoundSrc && geodDst) {
+        auto datum = geodDst->datum();
+        if (datum) {
+            auto cs =
+                cs::EllipsoidalCS::createLatitudeLongitudeEllipsoidalHeight(
+                    common::UnitOfMeasure::DEGREE,
+                    common::UnitOfMeasure::METRE);
+            auto intermGeog3DCRS = util::nn_static_pointer_cast<crs::CRS>(
+                crs::GeographicCRS::create(
+                    util::PropertyMap()
+                        .set(common::IdentifiedObject::NAME_KEY,
+                             geodDst->nameStr())
+                        .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+                             metadata::Extent::WORLD),
+                    NN_NO_CHECK(datum), cs));
+            auto sourceToGeog3DOps =
+                createOperations(sourceCRS, intermGeog3DCRS, context);
+            auto geog3DToTargetOps =
+                createOperations(intermGeog3DCRS, targetCRS, context);
+            if (!geog3DToTargetOps.empty()) {
+                for (const auto &op : sourceToGeog3DOps) {
+                    res.emplace_back(
+                        ConcatenatedOperation::createComputeMetadata(
+                            {op, geog3DToTargetOps.front()},
+                            !allowEmptyIntersection));
+                }
+                return res;
+            }
+        }
     }
 
     // reverse of previous case
     auto compoundDst = dynamic_cast<const crs::CompoundCRS *>(targetCRS.get());
-    if (geogSrc && compoundDst) {
+    if (geodSrc && compoundDst) {
         return applyInverse(createOperations(targetCRS, sourceCRS, context));
     }
 
-- 
cgit v1.2.3


From 70c9d0f75c9ecbb6f1ae9b5c7cefbe3a6fdbd5b9 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sat, 16 Mar 2019 17:26:17 +0100
Subject: createOperations(): fix nadgrids+geoidgrids -> nadgrids+geoidgrids

---
 src/iso19111/coordinateoperation.cpp | 43 ++++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 240833df..c1b3704c 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -11879,6 +11879,34 @@ CoordinateOperationFactory::Private::createOperations(
             }
         }
 
+        auto vertCRSOfBaseOfBoundSrc =
+            boundSrc->baseCRS()->extractVerticalCRS();
+        auto vertCRSOfBaseOfBoundDst =
+            boundDst->baseCRS()->extractVerticalCRS();
+        if (hubSrcGeog && hubDstGeog &&
+            hubSrcGeog->_isEquivalentTo(
+                hubDstGeog, util::IComparable::Criterion::EQUIVALENT) &&
+            vertCRSOfBaseOfBoundSrc && vertCRSOfBaseOfBoundDst) {
+            auto opsFirst = createOperations(sourceCRS, hubSrc, 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, opLast},
+                                    !allowEmptyIntersection));
+                        } catch (const InvalidOperationEmptyIntersection &) {
+                        }
+                    }
+                }
+                if (!res.empty()) {
+                    return res;
+                }
+            }
+        }
+
         return createOperations(boundSrc->baseCRS(), boundDst->baseCRS(),
                                 context);
     }
@@ -11987,6 +12015,21 @@ CoordinateOperationFactory::Private::createOperations(
                                         nn_interpTransformCRS));
                             }
                         }
+                    } else {
+                        auto compSrc0BoundCrs = dynamic_cast<crs::BoundCRS *>(
+                            componentsSrc[0].get());
+                        auto compDst0BoundCrs = dynamic_cast<crs::BoundCRS *>(
+                            componentsDst[0].get());
+                        if (compSrc0BoundCrs && compDst0BoundCrs &&
+                            dynamic_cast<crs::GeographicCRS *>(
+                                compSrc0BoundCrs->hubCRS().get()) &&
+                            compSrc0BoundCrs->hubCRS()->_isEquivalentTo(
+                                compDst0BoundCrs->hubCRS().get())) {
+                            interpolationGeogCRS =
+                                NN_NO_CHECK(util::nn_dynamic_pointer_cast<
+                                            crs::GeographicCRS>(
+                                    compSrc0BoundCrs->hubCRS()));
+                        }
                     }
                     auto opSrcCRSToGeogCRS = createOperations(
                         componentsSrc[0], interpolationGeogCRS, context);
-- 
cgit v1.2.3


From b5fb5f070793c29bb9b60faec9d221b25e7017d4 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sat, 16 Mar 2019 23:11:25 +0100
Subject: createOperations(): fix nadgrids -> nadgrids+geoidgrids

---
 src/iso19111/coordinateoperation.cpp | 83 ++++++++++++++++++++++++++++++++++++
 1 file changed, 83 insertions(+)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index c1b3704c..7ee20758 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -12054,6 +12054,89 @@ CoordinateOperationFactory::Private::createOperations(
         }
     }
 
+    // '+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +type=crs' to
+    // '+proj=longlat +ellps=GRS80 +nadgrids=@bar.gsb +geoidgrids=@bar.gtx
+    // +type=crs'
+    if (boundSrc && compoundDst) {
+        const auto &componentsDst = compoundDst->componentReferenceSystems();
+        if (!componentsDst.empty()) {
+            auto compDst0BoundCrs =
+                dynamic_cast<crs::BoundCRS *>(componentsDst[0].get());
+            if (compDst0BoundCrs) {
+                auto boundSrcHubAsGeogCRS = dynamic_cast<crs::GeographicCRS *>(
+                    boundSrc->hubCRS().get());
+                auto compDst0BoundCrsHubAsGeogCRS =
+                    dynamic_cast<crs::GeographicCRS *>(
+                        compDst0BoundCrs->hubCRS().get());
+                if (boundSrcHubAsGeogCRS && compDst0BoundCrsHubAsGeogCRS) {
+                    const auto &boundSrcHubAsGeogCRSDatum =
+                        boundSrcHubAsGeogCRS->datum();
+                    const auto &compDst0BoundCrsHubAsGeogCRSDatum =
+                        compDst0BoundCrsHubAsGeogCRS->datum();
+                    if (boundSrcHubAsGeogCRSDatum &&
+                        compDst0BoundCrsHubAsGeogCRSDatum &&
+                        boundSrcHubAsGeogCRSDatum->_isEquivalentTo(
+                            compDst0BoundCrsHubAsGeogCRSDatum.get())) {
+                        auto cs = cs::EllipsoidalCS::
+                            createLatitudeLongitudeEllipsoidalHeight(
+                                common::UnitOfMeasure::DEGREE,
+                                common::UnitOfMeasure::METRE);
+                        auto intermGeog3DCRS = util::nn_static_pointer_cast<
+                            crs::CRS>(crs::GeographicCRS::create(
+                            util::PropertyMap()
+                                .set(common::IdentifiedObject::NAME_KEY,
+                                     boundSrcHubAsGeogCRS->nameStr())
+                                .set(
+                                    common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+                                    metadata::Extent::WORLD),
+                            NN_NO_CHECK(boundSrcHubAsGeogCRSDatum), cs));
+                        auto sourceToGeog3DOps = createOperations(
+                            sourceCRS, intermGeog3DCRS, context);
+                        auto geog3DToTargetOps = createOperations(
+                            intermGeog3DCRS, targetCRS, context);
+                        for (const auto &opSrc : sourceToGeog3DOps) {
+                            for (const auto &opDst : geog3DToTargetOps) {
+                                if (opSrc->targetCRS() && opDst->sourceCRS() &&
+                                    !opSrc->targetCRS()->_isEquivalentTo(
+                                        opDst->sourceCRS().get())) {
+                                    // Shouldn't happen normally, but typically
+                                    // one of them can be 2D and the other 3D
+                                    // due to above createOperations() not
+                                    // exactly setting the expected source and
+                                    // target CRS.
+                                    // So create an adapter operation...
+                                    auto intermOps = createOperations(
+                                        NN_NO_CHECK(opSrc->targetCRS()),
+                                        NN_NO_CHECK(opDst->sourceCRS()),
+                                        context);
+                                    if (!intermOps.empty()) {
+                                        res.emplace_back(
+                                            ConcatenatedOperation::
+                                                createComputeMetadata(
+                                                    {opSrc, intermOps.front(),
+                                                     opDst},
+                                                    !allowEmptyIntersection));
+                                    }
+                                } else {
+                                    res.emplace_back(
+                                        ConcatenatedOperation::
+                                            createComputeMetadata(
+                                                {opSrc, opDst},
+                                                !allowEmptyIntersection));
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // reverse of previous case
+    if (boundDst && compoundSrc) {
+        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+    }
+
     return res;
 }
 //! @endcond
-- 
cgit v1.2.3


From 11e3047990c46cd5b705b20caaa70c37713e229f Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sat, 23 Mar 2019 21:29:22 +0100
Subject: Fix GCC 9 warning about useless std::move()

---
 src/iso19111/coordinateoperation.cpp | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 7ee20758..57cf9832 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -4816,7 +4816,7 @@ ConversionPtr Conversion::convertToOtherMethod(int targetEPSGCode) const {
             common::Length(
                 parameterValueMeasure(EPSG_CODE_PARAMETER_FALSE_NORTHING)));
         conv->setCRSs(this, false);
-        return std::move(conv);
+        return conv.as_nullable();
     }
 
     if (current_epsg_code == EPSG_CODE_METHOD_MERCATOR_VARIANT_B &&
@@ -4837,7 +4837,7 @@ ConversionPtr Conversion::convertToOtherMethod(int targetEPSGCode) const {
             common::Length(
                 parameterValueMeasure(EPSG_CODE_PARAMETER_FALSE_NORTHING)));
         conv->setCRSs(this, false);
-        return std::move(conv);
+        return conv.as_nullable();
     }
 
     if (current_epsg_code == EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_1SP &&
@@ -4872,7 +4872,7 @@ ConversionPtr Conversion::convertToOtherMethod(int targetEPSGCode) const {
                 common::Length(
                     parameterValueMeasure(EPSG_CODE_PARAMETER_FALSE_NORTHING)));
             conv->setCRSs(this, false);
-            return std::move(conv);
+            return conv.as_nullable();
         } else {
             const double K = k0 * m0 / std::pow(t0, n);
             const double phi1 =
@@ -4928,7 +4928,7 @@ ConversionPtr Conversion::convertToOtherMethod(int targetEPSGCode) const {
                             EPSG_CODE_PARAMETER_FALSE_EASTING)),
                         common::Length(FN_corrected_rounded));
                     conv->setCRSs(this, false);
-                    return std::move(conv);
+                    return conv.as_nullable();
                 }
             }
 
@@ -4942,7 +4942,7 @@ ConversionPtr Conversion::convertToOtherMethod(int targetEPSGCode) const {
                     parameterValueMeasure(EPSG_CODE_PARAMETER_FALSE_EASTING)),
                 common::Length(FN));
             conv->setCRSs(this, false);
-            return std::move(conv);
+            return conv.as_nullable();
         }
     }
 
@@ -5006,7 +5006,7 @@ ConversionPtr Conversion::convertToOtherMethod(int targetEPSGCode) const {
                     EPSG_CODE_PARAMETER_NORTHING_FALSE_ORIGIN) +
                 (std::fabs(FN_correction) > 1e-8 ? FN_correction : 0)));
         conv->setCRSs(this, false);
-        return std::move(conv);
+        return conv.as_nullable();
     }
 
     return nullptr;
-- 
cgit v1.2.3


From 8763ea7f01bd349df29c5c4ce3b4edd6252eff37 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sun, 24 Mar 2019 20:01:45 +0100
Subject: WKT2 parser: update to OGC 18-010r6

- Allow ID[] in base CRS of Derived CRS
- Allow VERSION[] in non-conversion coordinate operations
- Use VERSION[] to set operationVersion member of CoordinateOperation
- Export operationVersion in WKT2:2018
---
 src/iso19111/coordinateoperation.cpp | 37 +++++++++++++++++++++++++++++++-----
 1 file changed, 32 insertions(+), 5 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 57cf9832..80c1a572 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -756,6 +756,15 @@ void CoordinateOperation::setHasBallparkTransformation(bool b) {
 
 // ---------------------------------------------------------------------------
 
+void CoordinateOperation::setProperties(
+    const util::PropertyMap &properties) // throw(InvalidValueTypeException)
+{
+    ObjectUsage::setProperties(properties);
+    properties.getStringValue(OPERATION_VERSION_KEY, d->operationVersion_);
+}
+
+// ---------------------------------------------------------------------------
+
 //! @cond Doxygen_Suppress
 struct OperationMethod::Private {
     util::optional<std::string> formula_{};
@@ -6175,17 +6184,17 @@ TransformationNNPtr Transformation::create(
         throw InvalidOperation(
             "Inconsistent number of parameters and parameter values");
     }
-    auto conv = Transformation::nn_make_shared<Transformation>(
+    auto transf = Transformation::nn_make_shared<Transformation>(
         sourceCRSIn, targetCRSIn, interpolationCRSIn, methodIn, values,
         accuracies);
-    conv->assignSelf(conv);
-    conv->setProperties(properties);
+    transf->assignSelf(transf);
+    transf->setProperties(properties);
     std::string name;
     if (properties.getStringValue(common::IdentifiedObject::NAME_KEY, name) &&
         ci_find(name, "ballpark") != std::string::npos) {
-        conv->setHasBallparkTransformation(true);
+        transf->setHasBallparkTransformation(true);
     }
-    return conv;
+    return transf;
 }
 
 // ---------------------------------------------------------------------------
@@ -7653,6 +7662,15 @@ void SingleOperation::exportTransformationToWKT(
 
     formatter->addQuotedString(nameStr());
 
+    if (isWKT2 && formatter->use2018Keywords()) {
+        const auto &version = operationVersion();
+        if (version.has_value()) {
+            formatter->startNode(io::WKTConstants::VERSION, false);
+            formatter->addQuotedString(*version);
+            formatter->endNode();
+        }
+    }
+
     if (!formatter->abridgedTransformation()) {
         formatter->startNode(io::WKTConstants::SOURCECRS, false);
         l_sourceCRS->_exportToWKT(formatter);
@@ -9310,6 +9328,15 @@ void ConcatenatedOperation::_exportToWKT(io::WKTFormatter *formatter) const {
                          !identifiers().empty());
     formatter->addQuotedString(nameStr());
 
+    if (isWKT2 && formatter->use2018Keywords()) {
+        const auto &version = operationVersion();
+        if (version.has_value()) {
+            formatter->startNode(io::WKTConstants::VERSION, false);
+            formatter->addQuotedString(*version);
+            formatter->endNode();
+        }
+    }
+
     formatter->startNode(io::WKTConstants::SOURCECRS, false);
     sourceCRS()->_exportToWKT(formatter);
     formatter->endNode();
-- 
cgit v1.2.3


From 09db4826d4a1e5df900cb4b93a4b3eae2c487cb9 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 25 Mar 2019 13:56:43 +0100
Subject: WKT2_2018: always export ID of SOURCECRS/TARGETCRS and STEPs

even if there is one on upper node

This is a particular logic allowed by paragraph 7.3.3 Identifier
of OGC 18-010r6
---
 src/iso19111/coordinateoperation.cpp | 88 ++++++++++++++++++++++++++++--------
 1 file changed, 69 insertions(+), 19 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 80c1a572..d7f138a4 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -7639,6 +7639,53 @@ void Transformation::_exportToWKT(io::WKTFormatter *formatter) const {
 
 // ---------------------------------------------------------------------------
 
+static void exportSourceCRSAndTargetCRSToWKT(const CoordinateOperation *co,
+                                             io::WKTFormatter *formatter) {
+    auto l_sourceCRS = co->sourceCRS();
+    assert(l_sourceCRS);
+    auto l_targetCRS = co->targetCRS();
+    assert(l_targetCRS);
+    const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
+    const bool canExportCRSId =
+        (isWKT2 && formatter->use2018Keywords() &&
+         !(formatter->idOnTopLevelOnly() && formatter->topLevelHasId()));
+
+    const bool hasDomains = !co->domains().empty();
+    if (hasDomains) {
+        formatter->pushDisableUsage();
+    }
+
+    formatter->startNode(io::WKTConstants::SOURCECRS, false);
+    if (canExportCRSId && !l_sourceCRS->identifiers().empty()) {
+        // fake that top node has no id, so that the sourceCRS id is
+        // considered
+        formatter->pushHasId(false);
+        l_sourceCRS->_exportToWKT(formatter);
+        formatter->popHasId();
+    } else {
+        l_sourceCRS->_exportToWKT(formatter);
+    }
+    formatter->endNode();
+
+    formatter->startNode(io::WKTConstants::TARGETCRS, false);
+    if (canExportCRSId && !l_targetCRS->identifiers().empty()) {
+        // fake that top node has no id, so that the targetCRS id is
+        // considered
+        formatter->pushHasId(false);
+        l_targetCRS->_exportToWKT(formatter);
+        formatter->popHasId();
+    } else {
+        l_targetCRS->_exportToWKT(formatter);
+    }
+    formatter->endNode();
+
+    if (hasDomains) {
+        formatter->popDisableUsage();
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 void SingleOperation::exportTransformationToWKT(
     io::WKTFormatter *formatter) const {
     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
@@ -7647,11 +7694,6 @@ void SingleOperation::exportTransformationToWKT(
             "Transformation can only be exported to WKT2");
     }
 
-    auto l_sourceCRS = sourceCRS();
-    assert(l_sourceCRS);
-    auto l_targetCRS = targetCRS();
-    assert(l_targetCRS);
-
     if (formatter->abridgedTransformation()) {
         formatter->startNode(io::WKTConstants::ABRIDGEDTRANSFORMATION,
                              !identifiers().empty());
@@ -7672,13 +7714,7 @@ void SingleOperation::exportTransformationToWKT(
     }
 
     if (!formatter->abridgedTransformation()) {
-        formatter->startNode(io::WKTConstants::SOURCECRS, false);
-        l_sourceCRS->_exportToWKT(formatter);
-        formatter->endNode();
-
-        formatter->startNode(io::WKTConstants::TARGETCRS, false);
-        l_targetCRS->_exportToWKT(formatter);
-        formatter->endNode();
+        exportSourceCRSAndTargetCRSToWKT(this, formatter);
     }
 
     method()->_exportToWKT(formatter);
@@ -9337,20 +9373,34 @@ void ConcatenatedOperation::_exportToWKT(io::WKTFormatter *formatter) const {
         }
     }
 
-    formatter->startNode(io::WKTConstants::SOURCECRS, false);
-    sourceCRS()->_exportToWKT(formatter);
-    formatter->endNode();
+    exportSourceCRSAndTargetCRSToWKT(this, formatter);
 
-    formatter->startNode(io::WKTConstants::TARGETCRS, false);
-    targetCRS()->_exportToWKT(formatter);
-    formatter->endNode();
+    const bool canExportOperationId =
+        !(formatter->idOnTopLevelOnly() && formatter->topLevelHasId());
+
+    const bool hasDomains = !domains().empty();
+    if (hasDomains) {
+        formatter->pushDisableUsage();
+    }
 
     for (const auto &operation : operations()) {
         formatter->startNode(io::WKTConstants::STEP, false);
-        operation->_exportToWKT(formatter);
+        if (canExportOperationId && !operation->identifiers().empty()) {
+            // fake that top node has no id, so that the operation id is
+            // considered
+            formatter->pushHasId(false);
+            operation->_exportToWKT(formatter);
+            formatter->popHasId();
+        } else {
+            operation->_exportToWKT(formatter);
+        }
         formatter->endNode();
     }
 
+    if (hasDomains) {
+        formatter->popDisableUsage();
+    }
+
     ObjectUsage::baseExportToWKT(formatter);
     formatter->endNode();
 }
-- 
cgit v1.2.3


From 399ebef6a66e57d389a7b7761e7373a4cebe86e5 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 25 Mar 2019 20:51:56 +0100
Subject: lookForGridInfo(): correctly return that a grid is present, if
 present on the file system, but not in the database

---
 src/iso19111/coordinateoperation.cpp | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 80c1a572..a8a25a6a 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -8742,13 +8742,25 @@ void Transformation::_exportToPROJString(
         if (fileParameter &&
             fileParameter->type() == ParameterValue::Type::FILENAME) {
             auto filename = fileParameter->valueFile();
-            if (isMethodInverseOf) {
+            bool doInversion = isMethodInverseOf;
+            if (!identifiers().empty() &&
+                *identifiers().front()->codeSpace() ==
+                    metadata::Identifier::EPSG &&
+                method()->nameStr() ==
+                    "Geographic3D to GravityRelatedHeight (US .gtx)" &&
+                ends_with(filename, ".gtx")) {
+                // gtx files, from straight EPSG definition, must be applied in
+                // reverse order for "Geographic3D to GravityRelatedHeight"
+                // method
+                doInversion = !doInversion;
+            }
+            if (doInversion) {
                 formatter->startInversion();
             }
             formatter->addStep("vgridshift");
             formatter->addParam("grids", filename);
             formatter->addParam("multiplier", 1.0);
-            if (isMethodInverseOf) {
+            if (doInversion) {
                 formatter->stopInversion();
             }
             return;
-- 
cgit v1.2.3


From af0e994486346d89cdbe68fff68ce3c8d27b0c94 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Tue, 26 Mar 2019 13:30:52 +0100
Subject: coordinateoperation.cpp: silent false positive about copy paste
 error. Coverity CID 193519

---
 src/iso19111/coordinateoperation.cpp | 1 +
 1 file changed, 1 insertion(+)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 50a6a003..ff104026 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -9163,6 +9163,7 @@ void ConcatenatedOperation::fixStepsDirection(
                     op = op->inverse();
                 }
             } else if (i + 1 < operationsInOut.size()) {
+                /* coverity[copy_paste_error] */
                 l_targetCRS = operationsInOut[i + 1]->sourceCRS();
                 if (l_targetCRS) {
                     op->setCRSs(concatOpSourceCRS, NN_NO_CHECK(l_targetCRS),
-- 
cgit v1.2.3


From 71c4439617a77051ec8049f08c08c61f89946e95 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Tue, 26 Mar 2019 13:32:46 +0100
Subject: coordinateoperation.cpp: remove dead code. Coverity CID 193522

---
 src/iso19111/coordinateoperation.cpp | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index ff104026..1abb74b3 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -7704,7 +7704,7 @@ void SingleOperation::exportTransformationToWKT(
 
     formatter->addQuotedString(nameStr());
 
-    if (isWKT2 && formatter->use2018Keywords()) {
+    if (formatter->use2018Keywords()) {
         const auto &version = operationVersion();
         if (version.has_value()) {
             formatter->startNode(io::WKTConstants::VERSION, false);
@@ -7719,10 +7719,8 @@ void SingleOperation::exportTransformationToWKT(
 
     method()->_exportToWKT(formatter);
 
-    const MethodMapping *mapping =
-        !isWKT2 ? getMapping(method().get()) : nullptr;
     for (const auto &paramValue : parameterValues()) {
-        paramValue->_exportToWKT(formatter, mapping);
+        paramValue->_exportToWKT(formatter, nullptr);
     }
 
     if (!formatter->abridgedTransformation()) {
-- 
cgit v1.2.3


From 5893c7ddef817df0c42c7ee636e6e083b0c65aed Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Thu, 28 Mar 2019 17:33:26 +0100
Subject: createOperations(): improve BoundCRS<-->non-bound-CRS case

Fixes #1388
Typically helps for
projinfo -s "+proj=longlat +ellps=GRS80 +towgs84=1,2,3 +type=crs" -t EPSG:4258
by researching operations from the pivot WGS84 implied by the towgs84 clause
to EPSG:4258.
---
 src/iso19111/coordinateoperation.cpp | 31 ++++++++++++++++++++++++++++++-
 1 file changed, 30 insertions(+), 1 deletion(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 1abb74b3..043d09bc 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -11689,7 +11689,6 @@ CoordinateOperationFactory::Private::createOperations(
         return applyInverse(createOperations(targetCRS, sourceCRS, context));
     }
 
-    // boundCRS to a geogCRS that is the same as the hubCRS
     auto boundSrc = dynamic_cast<const crs::BoundCRS *>(sourceCRS.get());
     auto geogDst = dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
     if (boundSrc && geogDst) {
@@ -11698,6 +11697,7 @@ CoordinateOperationFactory::Private::createOperations(
             dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
         auto geogCRSOfBaseOfBoundSrc =
             boundSrc->baseCRS()->extractGeographicCRS();
+        // Is it: boundCRS to a geogCRS that is the same as the hubCRS ?
         if (hubSrcGeog && geogCRSOfBaseOfBoundSrc &&
             (hubSrcGeog->_isEquivalentTo(
                  geogDst, util::IComparable::Criterion::EQUIVALENT) ||
@@ -11815,6 +11815,35 @@ CoordinateOperationFactory::Private::createOperations(
             return res;
         }
 
+        if (hubSrcGeog && geogCRSOfBaseOfBoundSrc) {
+            // This one should go to the above 'Is it: boundCRS to a geogCRS
+            // that is the same as the hubCRS ?' case
+            auto opsFirst = createOperations(sourceCRS, hubSrc, context);
+            auto opsLast = createOperations(hubSrc, targetCRS, context);
+            if (!opsFirst.empty() && !opsLast.empty()) {
+                for (const auto &opFirst : opsFirst) {
+                    for (const auto &opLast : opsLast) {
+                        // Exclude artificial transformations from the hub
+                        // to the target CRS
+                        if (!opLast->hasBallparkTransformation()) {
+                            try {
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {opFirst, opLast},
+                                            !allowEmptyIntersection));
+                            } catch (
+                                const InvalidOperationEmptyIntersection &) {
+                            }
+                        }
+                    }
+                }
+                if (!res.empty()) {
+                    return res;
+                }
+            }
+        }
+
         return createOperations(boundSrc->baseCRS(), targetCRS, context);
     }
 
-- 
cgit v1.2.3


From 6a7e24dce79f93b73f4919f267df2fdf3ee95713 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Thu, 28 Mar 2019 15:26:00 +0100
Subject: Add proj_normalize_for_visualization()

Fixes #1301

This function takes the output PJ from proj_create_crs_to_crs(),
and add (or undo) the needed axis swap operations so that the
object returned by proj_normalize_for_visualization() has the usual
GIS axis order.
In this implementation, this does something only if the coordinate
system of the source or target CRS, geographic or projected, has
NORTH, EAST ordering.
CompoundCRS wrapping those objects are also handled.
---
 src/iso19111/coordinateoperation.cpp | 40 ++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 1abb74b3..58e97272 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -765,6 +765,46 @@ void CoordinateOperation::setProperties(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Return a variation of the current coordinate operation whose axis
+ * order is the one expected for visualization purposes.
+ */
+CoordinateOperationNNPtr
+CoordinateOperation::normalizeForVisualization() const {
+    auto l_sourceCRS = sourceCRS();
+    auto l_targetCRS = targetCRS();
+    if (!l_sourceCRS || !l_targetCRS) {
+        throw util::UnsupportedOperationException(
+            "Cannot retrieve source or target CRS");
+    }
+    const bool swapSource =
+        l_sourceCRS->mustAxisOrderBeSwitchedForVisualization();
+    const bool swapTarget =
+        l_targetCRS->mustAxisOrderBeSwitchedForVisualization();
+    auto l_this = NN_NO_CHECK(std::dynamic_pointer_cast<CoordinateOperation>(
+        shared_from_this().as_nullable()));
+    if (!swapSource && !swapTarget) {
+        return l_this;
+    }
+    std::vector<CoordinateOperationNNPtr> subOps;
+    if (swapSource) {
+        auto op = Conversion::createAxisOrderReversal(false);
+        op->setCRSs(l_sourceCRS->normalizeForVisualization(),
+                    NN_NO_CHECK(l_sourceCRS), nullptr);
+        subOps.emplace_back(op);
+    }
+    subOps.emplace_back(l_this);
+    if (swapTarget) {
+        auto op = Conversion::createAxisOrderReversal(false);
+        op->setCRSs(NN_NO_CHECK(l_targetCRS),
+                    l_targetCRS->normalizeForVisualization(), nullptr);
+        subOps.emplace_back(op);
+    }
+    return util::nn_static_pointer_cast<CoordinateOperation>(
+        ConcatenatedOperation::createComputeMetadata(subOps, true));
+}
+
+// ---------------------------------------------------------------------------
+
 //! @cond Doxygen_Suppress
 struct OperationMethod::Private {
     util::optional<std::string> formula_{};
-- 
cgit v1.2.3


From 88f2661c9fd9fbf9d714a184243340cafb64d91d Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Fri, 29 Mar 2019 13:30:57 +0100
Subject: createOperations(): improve behaviour with input CRS from WKT that
 lacks intermediate IDs (fixes #1343)

---
 src/iso19111/coordinateoperation.cpp | 147 ++++++++++++++++++++++++++++++-----
 1 file changed, 127 insertions(+), 20 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 24e5f8ab..15532a89 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -11311,16 +11311,37 @@ static bool hasIdentifiers(const CoordinateOperationNNPtr &op) {
 
 static std::vector<crs::CRSNNPtr>
 findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
-                             const datum::GeodeticReferenceFramePtr &datum) {
+                             const datum::GeodeticReferenceFrame *datum) {
     std::vector<crs::CRSNNPtr> candidates;
-    for (const auto &id : datum->identifiers()) {
-        const auto &authName = *(id->codeSpace());
-        const auto &code = id->code();
-        if (!authName.empty()) {
-            auto l_candidates = authFactory->createGeodeticCRSFromDatum(
-                authName, code, std::string());
-            for (const auto &candidate : l_candidates) {
-                candidates.emplace_back(candidate);
+    assert(datum);
+    const auto &ids = datum->identifiers();
+    const auto &datumName = datum->nameStr();
+    if (!ids.empty()) {
+        for (const auto &id : ids) {
+            const auto &authName = *(id->codeSpace());
+            const auto &code = id->code();
+            if (!authName.empty()) {
+                auto l_candidates = authFactory->createGeodeticCRSFromDatum(
+                    authName, code, std::string());
+                for (const auto &candidate : l_candidates) {
+                    candidates.emplace_back(candidate);
+                }
+            }
+        }
+    } else if (datumName != "unknown" && datumName != "unnamed") {
+        auto matches = authFactory->createObjectsFromName(
+            datumName,
+            {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME}, false,
+            2);
+        if (matches.size() == 1) {
+            const auto &match = matches.front();
+            if (datum->_isEquivalentTo(
+                    match.get(), util::IComparable::Criterion::EQUIVALENT) &&
+                !match->identifiers().empty()) {
+                return findCandidateGeodCRSForDatum(
+                    authFactory,
+                    dynamic_cast<const datum::GeodeticReferenceFrame *>(
+                        match.get()));
             }
         }
     }
@@ -11362,9 +11383,9 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 
     const auto &authFactory = context.context->getAuthorityFactory();
     const auto candidatesSrcGeod(
-        findCandidateGeodCRSForDatum(authFactory, geodSrc->datum()));
+        findCandidateGeodCRSForDatum(authFactory, geodSrc->datum().get()));
     const auto candidatesDstGeod(
-        findCandidateGeodCRSForDatum(authFactory, geodDst->datum()));
+        findCandidateGeodCRSForDatum(authFactory, geodDst->datum().get()));
 
     auto createTransformations = [&](const crs::CRSNNPtr &candidateSrcGeod,
                                      const crs::CRSNNPtr &candidateDstGeod,
@@ -11605,12 +11626,8 @@ CoordinateOperationFactory::Private::createOperations(
             // but transformations are only available between their
             // corresponding geocentric CRS.
             const auto &srcDatum = geodSrc->datum();
-            const bool srcHasDatumWithId =
-                srcDatum && !srcDatum->identifiers().empty();
             const auto &dstDatum = geodDst->datum();
-            const bool dstHasDatumWithId =
-                dstDatum && !dstDatum->identifiers().empty();
-            if (srcHasDatumWithId && dstHasDatumWithId &&
+            if (srcDatum != nullptr && dstDatum != nullptr &&
                 !srcDatum->_isEquivalentTo(
                     dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) {
                 createOperationsWithDatumPivot(res, sourceCRS, targetCRS,
@@ -11955,6 +11972,31 @@ CoordinateOperationFactory::Private::createOperations(
     // A bit odd case as we are comparing apples to oranges, but in case
     // the vertical unit differ, do something useful.
     if (vertSrc && geogDst) {
+
+        if (vertSrc->identifiers().empty()) {
+            const auto &authFactory = context.context->getAuthorityFactory();
+            const auto &vertSrcName = vertSrc->nameStr();
+            if (authFactory != nullptr && vertSrcName != "unnamed" &&
+                vertSrcName != "unknown") {
+                auto matches = authFactory->createObjectsFromName(
+                    vertSrcName,
+                    {io::AuthorityFactory::ObjectType::VERTICAL_CRS}, false, 2);
+                if (matches.size() == 1) {
+                    const auto &match = matches.front();
+                    if (vertSrc->_isEquivalentTo(
+                            match.get(),
+                            util::IComparable::Criterion::EQUIVALENT) &&
+                        !match->identifiers().empty()) {
+                        return createOperations(
+                            NN_NO_CHECK(
+                                util::nn_dynamic_pointer_cast<crs::VerticalCRS>(
+                                    match)),
+                            targetCRS, context);
+                    }
+                }
+            }
+        }
+
         const double convSrc =
             vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI();
         double convDst = 1.0;
@@ -12300,6 +12342,67 @@ CoordinateOperationFactory::Private::createOperations(
 
 // ---------------------------------------------------------------------------
 
+static crs::CRSNNPtr
+getResolvedCRS(const crs::CRSNNPtr &crs,
+               const CoordinateOperationContextNNPtr &context) {
+    const auto &authFactory = context->getAuthorityFactory();
+
+    auto projectedCrs = dynamic_cast<crs::ProjectedCRS *>(crs.get());
+    if (projectedCrs && authFactory) {
+        const auto &ids = projectedCrs->identifiers();
+        if (!ids.empty() && projectedCrs->baseCRS()->identifiers().empty()) {
+            const auto tmpAuthFactory = io::AuthorityFactory::create(
+                authFactory->databaseContext(), *ids.front()->codeSpace());
+            try {
+                crs::CRSNNPtr resolvedCrs(
+                    tmpAuthFactory->createProjectedCRS(ids.front()->code()));
+                if (resolvedCrs->_isEquivalentTo(
+                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
+                    return resolvedCrs;
+                }
+            } catch (const std::exception &) {
+            }
+        }
+    }
+
+    auto compoundCrs = dynamic_cast<crs::CompoundCRS *>(crs.get());
+    // If we get a CompoundCRS that has an EPSG code, but whose component CRS
+    // lack one, typically from WKT2, this might be an issue to get proper
+    // results in createOperations(), so import the CompoundCRS from the
+    // registry, and if equivalent to the original one, then use the version
+    // from the registry.
+    if (compoundCrs && authFactory) {
+        const auto &ids = compoundCrs->identifiers();
+        if (!ids.empty()) {
+            const auto &components = compoundCrs->componentReferenceSystems();
+            bool hasMissingId = false;
+            for (const auto &comp : components) {
+                if (comp->identifiers().empty()) {
+                    hasMissingId = true;
+                    break;
+                }
+            }
+            if (hasMissingId) {
+                const auto tmpAuthFactory = io::AuthorityFactory::create(
+                    authFactory->databaseContext(), *ids.front()->codeSpace());
+                try {
+                    crs::CRSNNPtr resolvedCrs(
+                        tmpAuthFactory->createCompoundCRS(ids.front()->code()));
+                    if (resolvedCrs->_isEquivalentTo(
+                            crs.get(),
+                            util::IComparable::Criterion::EQUIVALENT)) {
+                        return resolvedCrs;
+                    }
+                } catch (const std::exception &) {
+                }
+            }
+        }
+    }
+    return crs;
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Find a list of CoordinateOperation from sourceCRS to targetCRS.
  *
  * The operations are sorted with the most relevant ones first: by
@@ -12328,10 +12431,14 @@ CoordinateOperationFactory::createOperations(
     auto l_sourceCRS = srcBoundCRS ? NN_NO_CHECK(srcBoundCRS) : sourceCRS;
     auto l_targetCRS = targetBoundCRS ? NN_NO_CHECK(targetBoundCRS) : targetCRS;
 
-    Private::Context contextPrivate(sourceCRS, targetCRS, context);
-    return filterAndSort(
-        Private::createOperations(l_sourceCRS, l_targetCRS, contextPrivate),
-        context, l_sourceCRS, l_targetCRS);
+    auto l_resolvedSourceCRS = getResolvedCRS(l_sourceCRS, context);
+    auto l_resolvedTargetCRS = getResolvedCRS(l_targetCRS, context);
+    Private::Context contextPrivate(l_resolvedSourceCRS, l_resolvedTargetCRS,
+                                    context);
+    return filterAndSort(Private::createOperations(l_resolvedSourceCRS,
+                                                   l_resolvedTargetCRS,
+                                                   contextPrivate),
+                         context, l_resolvedSourceCRS, l_resolvedTargetCRS);
 }
 
 // ---------------------------------------------------------------------------
-- 
cgit v1.2.3


From d9e2a15f2e17b6710ccffa3e271595e006ceadf2 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Tue, 16 Apr 2019 19:26:16 +0200
Subject: createOperations(): do not attempt using a unrelated datum
 intermediate when doing geog2D<-->geog3D conversions of same datum

Seen when testing transformations between "CR 05" (EPSG:5365) and "CR-SIRGAS" (EPSG:8907)
which require going through their corresponding 3D GeogCRS to find a Helmert
transformation.
---
 src/iso19111/coordinateoperation.cpp | 326 ++++++++++++++++++++++++++++++-----
 1 file changed, 283 insertions(+), 43 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 15532a89..7eab849c 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -56,6 +56,10 @@
 #include <string>
 #include <vector>
 
+#ifdef DEBUG
+#include <iostream>
+#endif
+
 using namespace NS_PROJ::internal;
 
 #if 0
@@ -805,6 +809,14 @@ CoordinateOperation::normalizeForVisualization() const {
 
 // ---------------------------------------------------------------------------
 
+//! @cond Doxygen_Suppress
+CoordinateOperationNNPtr CoordinateOperation::shallowClone() const {
+    return _shallowClone();
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
 //! @cond Doxygen_Suppress
 struct OperationMethod::Private {
     util::optional<std::string> formula_{};
@@ -2342,6 +2354,10 @@ ConversionNNPtr Conversion::shallowClone() const {
     conv->setCRSs(this, false);
     return conv;
 }
+
+CoordinateOperationNNPtr Conversion::_shallowClone() const {
+    return util::nn_static_pointer_cast<CoordinateOperation>(shallowClone());
+}
 //! @endcond
 
 // ---------------------------------------------------------------------------
@@ -4679,6 +4695,16 @@ InverseConversion::create(const ConversionNNPtr &forward) {
     return conv;
 }
 
+// ---------------------------------------------------------------------------
+
+CoordinateOperationNNPtr InverseConversion::_shallowClone() const {
+    auto op = InverseConversion::nn_make_shared<InverseConversion>(
+        inverseAsConversion()->shallowClone());
+    op->assignSelf(op);
+    op->setCRSs(this, false);
+    return util::nn_static_pointer_cast<CoordinateOperation>(op);
+}
+
 //! @endcond
 
 // ---------------------------------------------------------------------------
@@ -5181,6 +5207,7 @@ const char *Conversion::getWKT1GDALMethodName() const {
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+
 void Conversion::_exportToWKT(io::WKTFormatter *formatter) const {
     const auto &l_method = method();
     const auto &methodName = l_method->nameStr();
@@ -5210,6 +5237,17 @@ void Conversion::_exportToWKT(io::WKTFormatter *formatter) const {
         formatter->pushOutputId(false);
     }
 
+#ifdef DEBUG_CONVERSION_ID
+    if (sourceCRS() && targetCRS()) {
+        formatter->startNode("SOURCECRS_ID", false);
+        sourceCRS()->formatID(formatter);
+        formatter->endNode();
+        formatter->startNode("TARGETCRS_ID", false);
+        targetCRS()->formatID(formatter);
+        formatter->endNode();
+    }
+#endif
+
     bool bAlreadyWritten = false;
     if (!isWKT2 && formatter->useESRIDialect()) {
         const ESRIParamMapping *esriParams = nullptr;
@@ -6036,10 +6074,14 @@ const crs::CRSNNPtr &Transformation::targetCRS() PROJ_PURE_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;
+    auto transf = Transformation::nn_make_shared<Transformation>(*this);
+    transf->assignSelf(transf);
+    transf->setCRSs(this, false);
+    return transf;
+}
+
+CoordinateOperationNNPtr Transformation::_shallowClone() const {
+    return util::nn_static_pointer_cast<CoordinateOperation>(shallowClone());
 }
 //! @endcond
 
@@ -7656,6 +7698,13 @@ InverseTransformation::create(const TransformationNNPtr &forward) {
 
 // ---------------------------------------------------------------------------
 
+TransformationNNPtr InverseTransformation::inverseAsTransformation() const {
+    return NN_NO_CHECK(
+        util::nn_dynamic_pointer_cast<Transformation>(forwardOperation_));
+}
+
+// ---------------------------------------------------------------------------
+
 void InverseTransformation::_exportToWKT(io::WKTFormatter *formatter) const {
 
     auto approxInverse = createApproximateInverseIfPossible(
@@ -7667,6 +7716,16 @@ void InverseTransformation::_exportToWKT(io::WKTFormatter *formatter) const {
     }
 }
 
+// ---------------------------------------------------------------------------
+
+CoordinateOperationNNPtr InverseTransformation::_shallowClone() const {
+    auto op = InverseTransformation::nn_make_shared<InverseTransformation>(
+        inverseAsTransformation()->shallowClone());
+    op->assignSelf(op);
+    op->setCRSs(this, false);
+    return util::nn_static_pointer_cast<CoordinateOperation>(op);
+}
+
 //! @endcond
 
 // ---------------------------------------------------------------------------
@@ -9063,6 +9122,7 @@ struct ConcatenatedOperation::Private {
 
     explicit Private(const std::vector<CoordinateOperationNNPtr> &operationsIn)
         : operations_(operationsIn) {}
+    Private(const Private &) = default;
 };
 //! @endcond
 
@@ -9074,6 +9134,14 @@ ConcatenatedOperation::~ConcatenatedOperation() = default;
 
 // ---------------------------------------------------------------------------
 
+//! @cond Doxygen_Suppress
+ConcatenatedOperation::ConcatenatedOperation(const ConcatenatedOperation &other)
+    : CoordinateOperation(other),
+      d(internal::make_unique<Private>(*(other.d))) {}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
 ConcatenatedOperation::ConcatenatedOperation(
     const std::vector<CoordinateOperationNNPtr> &operationsIn)
     : CoordinateOperation(), d(internal::make_unique<Private>(operationsIn)) {}
@@ -9136,6 +9204,22 @@ ConcatenatedOperationNNPtr ConcatenatedOperation::create(
         }
         if (i >= 1) {
             if (!compareStepCRS(l_sourceCRS.get(), lastTargetCRS.get())) {
+#ifdef DEBUG_CONCATENATED_OPERATION
+                {
+                    auto f(io::WKTFormatter::create(
+                        io::WKTFormatter::Convention::WKT2_2018));
+                    std::cerr << "Source CRS of step " << i << ":" << std::endl;
+                    std::cerr << l_sourceCRS->exportToWKT(f.get()) << std::endl;
+                }
+                {
+                    auto f(io::WKTFormatter::create(
+                        io::WKTFormatter::Convention::WKT2_2018));
+                    std::cerr << "Target CRS of step " << i - 1 << ":"
+                              << std::endl;
+                    std::cerr << lastTargetCRS->exportToWKT(f.get())
+                              << std::endl;
+                }
+#endif
                 throw InvalidOperation(
                     "Inconsistent chaining of CRS in operations");
             }
@@ -9149,6 +9233,14 @@ ConcatenatedOperationNNPtr ConcatenatedOperation::create(
     op->setCRSs(NN_NO_CHECK(operationsIn[0]->sourceCRS()),
                 NN_NO_CHECK(operationsIn.back()->targetCRS()), nullptr);
     op->setAccuracies(accuracies);
+#ifdef DEBUG_CONCATENATED_OPERATION
+    {
+        auto f(
+            io::WKTFormatter::create(io::WKTFormatter::Convention::WKT2_2018));
+        std::cerr << "ConcatenatedOperation::create()" << std::endl;
+        std::cerr << op->exportToWKT(f.get()) << std::endl;
+    }
+#endif
     return op;
 }
 
@@ -9459,6 +9551,18 @@ void ConcatenatedOperation::_exportToWKT(io::WKTFormatter *formatter) const {
 
 // ---------------------------------------------------------------------------
 
+//! @cond Doxygen_Suppress
+CoordinateOperationNNPtr ConcatenatedOperation::_shallowClone() const {
+    auto op =
+        ConcatenatedOperation::nn_make_shared<ConcatenatedOperation>(*this);
+    op->assignSelf(op);
+    op->setCRSs(this, false);
+    return util::nn_static_pointer_cast<CoordinateOperation>(op);
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
 void ConcatenatedOperation::_exportToPROJString(
     io::PROJStringFormatter *formatter) const // throw(FormattingException)
 {
@@ -11359,11 +11463,46 @@ static bool isNullTransformation(const std::string &name) {
 
 // ---------------------------------------------------------------------------
 
+#ifdef DEBUG
+
+static int nCallLevel = 0;
+
+struct EnterDebugLevel {
+    EnterDebugLevel() { ++nCallLevel; }
+    ~EnterDebugLevel() { --nCallLevel; }
+};
+
+static void debugTrace(const std::string &str) {
+    for (int i = 1; i < nCallLevel; i++)
+        std::cerr << " ";
+    std::cerr << str << std::endl;
+}
+
+static std::string objectAsStr(const common::IdentifiedObject *obj) {
+    std::string ret(obj->nameStr());
+    const auto &ids = obj->identifiers();
+    if (!ids.empty()) {
+        ret += " (";
+        ret += (*ids[0]->codeSpace()) + ":" + ids[0]->code();
+        ret += ")";
+    }
+    return ret;
+}
+#endif
+
+// ---------------------------------------------------------------------------
+
 void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
     std::vector<CoordinateOperationNNPtr> &res, const crs::CRSNNPtr &sourceCRS,
     const crs::CRSNNPtr &targetCRS, const crs::GeodeticCRS *geodSrc,
     const crs::GeodeticCRS *geodDst, Private::Context &context) {
 
+#ifdef DEBUG
+    EnterDebugLevel enterFunction;
+    debugTrace("createOperationsWithDatumPivot(" +
+               objectAsStr(sourceCRS.get()) + "," +
+               objectAsStr(targetCRS.get()) + ")");
+#endif
     const bool allowEmptyIntersection = true;
 
     struct CreateOperationsWithDatumPivotAntiRecursion {
@@ -11411,20 +11550,73 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
             }
 
             std::vector<CoordinateOperationNNPtr> subOps;
+            const bool isNullThird =
+                isNullTransformation(opsThird[0]->nameStr());
+            CoordinateOperationNNPtr opSecondCloned(
+                (isNullFirst || isNullThird) ? opSecond->shallowClone()
+                                             : opSecond);
+            CoordinateOperation *invCOForward = nullptr;
+            if (isNullFirst || isNullThird) {
+                if (opSecondCloned->identifiers().size() == 1 &&
+                    (*opSecondCloned->identifiers()[0]->codeSpace())
+                            .find("DERIVED_FROM") == std::string::npos) {
+                    {
+                        util::PropertyMap map;
+                        addModifiedIdentifier(map, opSecondCloned.get(), false,
+                                              true);
+                        opSecondCloned->setProperties(map);
+                    }
+                    auto invCO = dynamic_cast<InverseCoordinateOperation *>(
+                        opSecondCloned.get());
+                    if (invCO) {
+                        invCOForward = invCO->forwardOperation().get();
+                        if (invCOForward->identifiers().size() == 1 &&
+                            (*invCOForward->identifiers()[0]->codeSpace())
+                                    .find("DERIVED_FROM") ==
+                                std::string::npos) {
+                            util::PropertyMap map;
+                            addModifiedIdentifier(map, invCOForward, false,
+                                                  true);
+                            invCOForward->setProperties(map);
+                        }
+                    }
+                }
+            }
             if (isNullFirst) {
-                opSecond->setCRSs(
-                    sourceCRS, NN_CHECK_ASSERT(opSecond->targetCRS()), nullptr);
+                auto oldTarget(NN_CHECK_ASSERT(opSecondCloned->targetCRS()));
+                opSecondCloned->setCRSs(sourceCRS, oldTarget, nullptr);
+                if (invCOForward) {
+                    invCOForward->setCRSs(oldTarget, sourceCRS, nullptr);
+                }
             } else {
                 subOps.emplace_back(opFirst);
             }
-            if (isNullTransformation(opsThird[0]->nameStr())) {
-                opSecond->setCRSs(NN_CHECK_ASSERT(opSecond->sourceCRS()),
-                                  targetCRS, nullptr);
-                subOps.emplace_back(opSecond);
+            if (isNullThird) {
+                auto oldSource(NN_CHECK_ASSERT(opSecondCloned->sourceCRS()));
+                opSecondCloned->setCRSs(oldSource, targetCRS, nullptr);
+                if (invCOForward) {
+                    invCOForward->setCRSs(targetCRS, oldSource, nullptr);
+                }
+                subOps.emplace_back(opSecondCloned);
             } else {
-                subOps.emplace_back(opSecond);
+                subOps.emplace_back(opSecondCloned);
                 subOps.emplace_back(opsThird[0]);
             }
+#ifdef DEBUG
+            std::string debugStr;
+            for (const auto &op : subOps) {
+                if (!debugStr.empty()) {
+                    debugStr += " + ";
+                }
+                debugStr += objectAsStr(op.get());
+                debugStr += " (";
+                debugStr += objectAsStr(op->sourceCRS().get());
+                debugStr += "->";
+                debugStr += objectAsStr(op->targetCRS().get());
+                debugStr += ")";
+            }
+            debugTrace("transformation " + debugStr);
+#endif
             res.emplace_back(ConcatenatedOperation::createComputeMetadata(
                 subOps, !allowEmptyIntersection));
         }
@@ -11436,6 +11628,14 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
         if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
             for (const auto &candidateDstGeod : candidatesDstGeod) {
                 if (candidateDstGeod->nameStr() == targetCRS->nameStr()) {
+#ifdef DEBUG
+                    EnterDebugLevel loopLevel;
+                    debugTrace("try " + objectAsStr(sourceCRS.get()) + "->" +
+                               objectAsStr(candidateSrcGeod.get()) + "->" +
+                               objectAsStr(candidateDstGeod.get()) + "->" +
+                               objectAsStr(targetCRS.get()) + ")");
+                    EnterDebugLevel loopLevel2;
+#endif
                     const auto opsFirst =
                         createOperations(sourceCRS, candidateSrcGeod, context);
                     assert(!opsFirst.empty());
@@ -11454,17 +11654,28 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
     }
 
     for (const auto &candidateSrcGeod : candidatesSrcGeod) {
+#ifdef DEBUG
+        EnterDebugLevel loopLevel;
+#endif
         const auto opsFirst =
             createOperations(sourceCRS, candidateSrcGeod, context);
         assert(!opsFirst.empty());
         const bool isNullFirst = isNullTransformation(opsFirst[0]->nameStr());
 
         for (const auto &candidateDstGeod : candidatesDstGeod) {
+#ifdef DEBUG
+            EnterDebugLevel loopLevel2;
+            debugTrace("try " + objectAsStr(sourceCRS.get()) + "->" +
+                       objectAsStr(candidateSrcGeod.get()) + "->" +
+                       objectAsStr(candidateDstGeod.get()) + "->" +
+                       objectAsStr(targetCRS.get()) + ")");
+            EnterDebugLevel loopLevel3;
+#endif
             createTransformations(candidateSrcGeod, candidateDstGeod,
                                   opsFirst[0], isNullFirst);
-        }
-        if (!res.empty()) {
-            return;
+            if (!res.empty()) {
+                return;
+            }
         }
     }
 }
@@ -11519,6 +11730,13 @@ CoordinateOperationFactory::Private::createOperations(
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
     Private::Context &context) {
 
+#ifdef DEBUG
+    EnterDebugLevel enterFunction;
+    auto debugStr("createOperations(" + objectAsStr(sourceCRS.get()) + "," +
+                  objectAsStr(targetCRS.get()) + ")");
+    debugTrace(debugStr);
+#endif
+
     std::vector<CoordinateOperationNNPtr> res;
     const bool allowEmptyIntersection = true;
 
@@ -11598,15 +11816,48 @@ CoordinateOperationFactory::Private::createOperations(
 
             doFilterAndCheckPerfectOp = false;
 
+            bool sameGeodeticDatum = false;
+            if (geodSrc && geodDst) {
+                const auto &srcDatum = geodSrc->datum();
+                const auto &dstDatum = geodDst->datum();
+                sameGeodeticDatum =
+                    srcDatum != nullptr && dstDatum != nullptr &&
+                    srcDatum->_isEquivalentTo(
+                        dstDatum.get(),
+                        util::IComparable::Criterion::EQUIVALENT);
+            }
+
+            if (res.empty() && !sameGeodeticDatum &&
+                !context.inCreateOperationsWithDatumPivotAntiRecursion &&
+                geodSrc && geodDst) {
+                // If we still didn't find a transformation, and that the source
+                // and target are GeodeticCRS, then go through their underlying
+                // datum to find potential transformations between other
+                // GeodeticRSs
+                // that are made of those datum
+                // The typical example is if transforming between two
+                // GeographicCRS,
+                // but transformations are only available between their
+                // corresponding geocentric CRS.
+                const auto &srcDatum = geodSrc->datum();
+                const auto &dstDatum = geodDst->datum();
+                if (srcDatum != nullptr && dstDatum != nullptr) {
+                    createOperationsWithDatumPivot(res, sourceCRS, targetCRS,
+                                                   geodSrc, geodDst, context);
+                    doFilterAndCheckPerfectOp = !res.empty();
+                }
+            }
+
             // NAD27 to NAD83 has tens of results already. No need to look
             // for a pivot
-            if ((res.empty() &&
+            if (!sameGeodeticDatum &&
+                ((res.empty() &&
+                  context.context->getAllowUseIntermediateCRS() ==
+                      CoordinateOperationContext::IntermediateCRSUse::
+                          IF_NO_DIRECT_TRANSFORMATION) ||
                  context.context->getAllowUseIntermediateCRS() ==
-                     CoordinateOperationContext::IntermediateCRSUse::
-                         IF_NO_DIRECT_TRANSFORMATION) ||
-                context.context->getAllowUseIntermediateCRS() ==
-                    CoordinateOperationContext::IntermediateCRSUse::ALWAYS ||
-                getenv("PROJ_FORCE_SEARCH_PIVOT")) {
+                     CoordinateOperationContext::IntermediateCRSUse::ALWAYS ||
+                 getenv("PROJ_FORCE_SEARCH_PIVOT"))) {
                 auto resWithIntermediate = findsOpsInRegistryWithIntermediate(
                     sourceCRS, targetCRS, context.context);
                 res.insert(res.end(), resWithIntermediate.begin(),
@@ -11615,27 +11866,6 @@ CoordinateOperationFactory::Private::createOperations(
             }
         }
 
-        if (res.empty() &&
-            !context.inCreateOperationsWithDatumPivotAntiRecursion && geodSrc &&
-            geodDst) {
-            // If we still didn't find a transformation, and that the source
-            // and target are GeodeticCRS, then go through their underlying
-            // datum to find potential transformations between other GeodeticRSs
-            // that are made of those datum
-            // The typical example is if transforming between two GeographicCRS,
-            // but transformations are only available between their
-            // corresponding geocentric CRS.
-            const auto &srcDatum = geodSrc->datum();
-            const auto &dstDatum = geodDst->datum();
-            if (srcDatum != nullptr && dstDatum != nullptr &&
-                !srcDatum->_isEquivalentTo(
-                    dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) {
-                createOperationsWithDatumPivot(res, sourceCRS, targetCRS,
-                                               geodSrc, geodDst, context);
-                doFilterAndCheckPerfectOp = !res.empty();
-            }
-        }
-
         if (doFilterAndCheckPerfectOp) {
             // If we get at least a result with perfect accuracy, do not bother
             // generating synthetic transforms.
@@ -12459,8 +12689,9 @@ InverseCoordinateOperation::~InverseCoordinateOperation() = default;
 // ---------------------------------------------------------------------------
 
 InverseCoordinateOperation::InverseCoordinateOperation(
-    const CoordinateOperationNNPtr &forwardOperation, bool wktSupportsInversion)
-    : forwardOperation_(forwardOperation),
+    const CoordinateOperationNNPtr &forwardOperationIn,
+    bool wktSupportsInversion)
+    : forwardOperation_(forwardOperationIn),
       wktSupportsInversion_(wktSupportsInversion) {}
 
 // ---------------------------------------------------------------------------
@@ -12651,6 +12882,15 @@ void PROJBasedOperation::_exportToPROJString(
 
 // ---------------------------------------------------------------------------
 
+CoordinateOperationNNPtr PROJBasedOperation::_shallowClone() const {
+    auto op = PROJBasedOperation::nn_make_shared<PROJBasedOperation>(*this);
+    op->assignSelf(op);
+    op->setCRSs(this, false);
+    return util::nn_static_pointer_cast<CoordinateOperation>(op);
+}
+
+// ---------------------------------------------------------------------------
+
 std::set<GridDescription> PROJBasedOperation::gridsNeeded(
     const io::DatabaseContextPtr &databaseContext) const {
     std::set<GridDescription> res;
-- 
cgit v1.2.3


From d8526d7870e4b91238c9a7b652ed03c21b77e884 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 22 Apr 2019 13:15:02 +0200
Subject: ESRI_WKT: preserve Gauss_Kruger in conversion name for round-tripping

---
 src/iso19111/coordinateoperation.cpp | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 7eab849c..0b4e7912 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -5120,10 +5120,11 @@ static void getESRIMethodNameAndParams(const Conversion *conv,
             }
         } else if (esriMapping->epsg_code ==
                    EPSG_CODE_METHOD_TRANSVERSE_MERCATOR) {
-            if (l_targetCRS &&
-                (ci_find(l_targetCRS->nameStr(), "Gauss") !=
-                     std::string::npos ||
-                 ci_find(l_targetCRS->nameStr(), "GK_") != std::string::npos)) {
+            if (ci_find(conv->nameStr(), "Gauss Kruger") != std::string::npos ||
+                (l_targetCRS && (ci_find(l_targetCRS->nameStr(), "Gauss") !=
+                                     std::string::npos ||
+                                 ci_find(l_targetCRS->nameStr(), "GK_") !=
+                                     std::string::npos))) {
                 esriParams = paramsESRI_Gauss_Kruger;
                 esriMethodName = "Gauss_Kruger";
             } else {
-- 
cgit v1.2.3


From c3c1efedbeae11a09947f2873fe97e50afa87b29 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sun, 28 Apr 2019 13:00:18 +0200
Subject: Fix false-positive -Wnull-dereference GCC 8 warning

---
 src/iso19111/coordinateoperation.cpp | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 7eab849c..6a05c285 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -12584,11 +12584,11 @@ getResolvedCRS(const crs::CRSNNPtr &crs,
             const auto tmpAuthFactory = io::AuthorityFactory::create(
                 authFactory->databaseContext(), *ids.front()->codeSpace());
             try {
-                crs::CRSNNPtr resolvedCrs(
+                auto resolvedCrs(
                     tmpAuthFactory->createProjectedCRS(ids.front()->code()));
-                if (resolvedCrs->_isEquivalentTo(
+                if (resolvedCrs->isEquivalentTo(
                         crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
-                    return resolvedCrs;
+                    return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
                 }
             } catch (const std::exception &) {
             }
@@ -12616,12 +12616,13 @@ getResolvedCRS(const crs::CRSNNPtr &crs,
                 const auto tmpAuthFactory = io::AuthorityFactory::create(
                     authFactory->databaseContext(), *ids.front()->codeSpace());
                 try {
-                    crs::CRSNNPtr resolvedCrs(
+                    auto resolvedCrs(
                         tmpAuthFactory->createCompoundCRS(ids.front()->code()));
-                    if (resolvedCrs->_isEquivalentTo(
+                    if (resolvedCrs->isEquivalentTo(
                             crs.get(),
                             util::IComparable::Criterion::EQUIVALENT)) {
-                        return resolvedCrs;
+                        return util::nn_static_pointer_cast<crs::CRS>(
+                            resolvedCrs);
                     }
                 } catch (const std::exception &) {
                 }
-- 
cgit v1.2.3


From cca27b1fae234a90df42ff5341121759846dc39b Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Tue, 30 Apr 2019 14:31:16 +0200
Subject: Propagate ballpark transformation flag to inverse coordinate
 operations

---
 src/iso19111/coordinateoperation.cpp | 5 +++++
 1 file changed, 5 insertions(+)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 6a05c285..2cab05bd 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -7484,6 +7484,8 @@ Transformation::Private::registerInv(util::BaseObjectNNPtr thisIn,
                                      TransformationNNPtr invTransform) {
     invTransform->d->forwardOperation_ =
         util::nn_dynamic_pointer_cast<Transformation>(thisIn);
+    invTransform->setHasBallparkTransformation(
+        invTransform->d->forwardOperation_->hasBallparkTransformation());
     return invTransform;
 }
 //! @endcond
@@ -9490,6 +9492,7 @@ CoordinateOperationNNPtr ConcatenatedOperation::inverse() const {
     auto op =
         create(properties, inversedOperations, coordinateOperationAccuracies());
     op->d->computedName_ = d->computedName_;
+    op->setHasBallparkTransformation(hasBallparkTransformation());
     return op;
 }
 
@@ -12704,6 +12707,8 @@ void InverseCoordinateOperation::setPropertiesFromForward() {
     if (forwardOperation_->sourceCRS() && forwardOperation_->targetCRS()) {
         setCRSs(forwardOperation_.get(), true);
     }
+    setHasBallparkTransformation(
+        forwardOperation_->hasBallparkTransformation());
 }
 
 // ---------------------------------------------------------------------------
-- 
cgit v1.2.3


From b4a9e65cec051ca3cb16b8cccfa012d70ce10570 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Tue, 30 Apr 2019 14:56:32 +0200
Subject: createOperations(): in SourceTargetCRSExtentUse::INTERSECTION mode,
 early return if the intersection of the areas is empty

---
 src/iso19111/coordinateoperation.cpp | 11 +++++++++++
 1 file changed, 11 insertions(+)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 2cab05bd..ca882f49 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -12669,6 +12669,17 @@ CoordinateOperationFactory::createOperations(
     auto l_resolvedTargetCRS = getResolvedCRS(l_targetCRS, context);
     Private::Context contextPrivate(l_resolvedSourceCRS, l_resolvedTargetCRS,
                                     context);
+
+    if (context->getSourceAndTargetCRSExtentUse() ==
+        CoordinateOperationContext::SourceTargetCRSExtentUse::INTERSECTION) {
+        auto sourceCRSExtent(getExtent(l_resolvedSourceCRS));
+        auto targetCRSExtent(getExtent(l_resolvedTargetCRS));
+        if (sourceCRSExtent && targetCRSExtent &&
+            !sourceCRSExtent->intersects(NN_NO_CHECK(targetCRSExtent))) {
+            return std::vector<CoordinateOperationNNPtr>();
+        }
+    }
+
     return filterAndSort(Private::createOperations(l_resolvedSourceCRS,
                                                    l_resolvedTargetCRS,
                                                    contextPrivate),
-- 
cgit v1.2.3


From 09cbfb85c834d99e5a00f5989dc144613e0cfbf2 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Tue, 30 Apr 2019 17:25:38 +0200
Subject: WKT importer: accepts PROJ-based COORDINATEOPERATION

---
 src/iso19111/coordinateoperation.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index ca882f49..9c393aba 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -12802,7 +12802,7 @@ PROJBasedOperationNNPtr PROJBasedOperation::create(
 
     auto method = OperationMethod::create(
         util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
-                                "PROJ-based operation method (approximate) : " +
+                                "PROJ-based operation method (approximate): " +
                                     projString),
         std::vector<GeneralOperationParameterNNPtr>{});
     auto op = PROJBasedOperation::nn_make_shared<PROJBasedOperation>(method);
-- 
cgit v1.2.3


From 5e98fed78205605ccb01ab4310d3cba292de73b4 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 6 May 2019 12:43:04 +0200
Subject: createOperations(): fix case of ETRS89 3D to proj string with
 nadgrids and geoidgrids

Fixes https://lists.osgeo.org/pipermail/proj/2019-May/008519.html
---
 src/iso19111/coordinateoperation.cpp | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 605004b6..f75f7588 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -12135,6 +12135,37 @@ CoordinateOperationFactory::Private::createOperations(
             }
         }
 
+        auto vertCRSOfBaseOfBoundSrc =
+            dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get());
+        if (vertCRSOfBaseOfBoundSrc && hubSrcGeog &&
+            hubSrcGeog->coordinateSystem()->axisList().size() == 3 &&
+            geogDst->coordinateSystem()->axisList().size() == 3) {
+            auto opsFirst = createOperations(sourceCRS, hubSrc, context);
+            auto opsSecond = createOperations(hubSrc, targetCRS, context);
+            if (!opsFirst.empty() && !opsSecond.empty()) {
+                for (const auto &opFirst : opsFirst) {
+                    for (const auto &opLast : opsSecond) {
+                        // Exclude artificial transformations from the hub
+                        // to the target CRS
+                        if (!opLast->hasBallparkTransformation()) {
+                            try {
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {opFirst, opLast},
+                                            !allowEmptyIntersection));
+                            } catch (
+                                const InvalidOperationEmptyIntersection &) {
+                            }
+                        }
+                    }
+                }
+                if (!res.empty()) {
+                    return res;
+                }
+            }
+        }
+
         return createOperations(boundSrc->baseCRS(), targetCRS, context);
     }
 
-- 
cgit v1.2.3


From 61cf8c5b29c82ab7e46b207bd125eaad49c03021 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 6 May 2019 19:45:27 +0200
Subject: createOperations(): for 'Amersfoort / RD New + NAP height'
 (EPSG:7415) to ETRS89 (EPSG:4937), make sure that the vgridshift is applied
 first (ie on Amersfoort datum) before the hgridshift

---
 src/iso19111/coordinateoperation.cpp | 111 ++++++++++++++++++++++++++++-------
 1 file changed, 89 insertions(+), 22 deletions(-)

(limited to 'src/iso19111/coordinateoperation.cpp')

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index f75f7588..348a776a 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -6958,11 +6958,11 @@ TransformationNNPtr Transformation::createNTv2(
 static TransformationNNPtr _createGravityRelatedHeightToGeographic3D(
     const util::PropertyMap &properties, bool inverse,
     const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn,
-    const std::string &filename,
+    const crs::CRSPtr &interpolationCRSIn, const std::string &filename,
     const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
 
     return Transformation::create(
-        properties, sourceCRSIn, targetCRSIn, nullptr,
+        properties, sourceCRSIn, targetCRSIn, interpolationCRSIn,
         util::PropertyMap().set(
             common::IdentifiedObject::NAME_KEY,
             inverse ? INVERSE_OF + PROJ_WKT2_NAME_METHOD_HEIGHT_TO_GEOG3D
@@ -6981,17 +6981,20 @@ static TransformationNNPtr _createGravityRelatedHeightToGeographic3D(
  * At minimum the name should be defined.
  * @param sourceCRSIn Source CRS.
  * @param targetCRSIn Target CRS.
+ * @param interpolationCRSIn Interpolation CRS. (might be null)
  * @param filename GRID filename.
  * @param accuracies Vector of positional accuracy (might be empty).
  * @return new Transformation.
  */
 TransformationNNPtr Transformation::createGravityRelatedHeightToGeographic3D(
     const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
-    const crs::CRSNNPtr &targetCRSIn, const std::string &filename,
+    const crs::CRSNNPtr &targetCRSIn, const crs::CRSPtr &interpolationCRSIn,
+    const std::string &filename,
     const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
 
     return _createGravityRelatedHeightToGeographic3D(
-        properties, false, sourceCRSIn, targetCRSIn, filename, accuracies);
+        properties, false, sourceCRSIn, targetCRSIn, interpolationCRSIn,
+        filename, accuracies);
 }
 
 // ---------------------------------------------------------------------------
@@ -7302,8 +7305,20 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
     auto targetCRS = op->targetCRS();
     std::string name;
     if (!forwardName.empty()) {
-        if (starts_with(forwardName, INVERSE_OF)) {
-            name = forwardName.substr(INVERSE_OF.size());
+        if (starts_with(forwardName, INVERSE_OF) ||
+            forwardName.find(" + ") != std::string::npos) {
+            auto tokens = split(forwardName, " + ");
+            for (size_t i = tokens.size(); i > 0;) {
+                i--;
+                if (!name.empty()) {
+                    name += " + ";
+                }
+                if (starts_with(tokens[i], INVERSE_OF)) {
+                    name += tokens[i].substr(INVERSE_OF.size());
+                } else {
+                    name += INVERSE_OF + tokens[i];
+                }
+            }
         } else if (!sourceCRS || !targetCRS ||
                    forwardName != buildOpName(opType, sourceCRS, targetCRS)) {
             name = INVERSE_OF + forwardName;
@@ -8195,13 +8210,14 @@ TransformationNNPtr Transformation::substitutePROJAlternativeGridNames(
                 return createGravityRelatedHeightToGeographic3D(
                            createPropertiesForInverse(self.as_nullable().get(),
                                                       true, false),
-                           targetCRS(), sourceCRS(), projFilename,
-                           coordinateOperationAccuracies())
+                           targetCRS(), sourceCRS(), interpolationCRS(),
+                           projFilename, coordinateOperationAccuracies())
                     ->inverseAsTransformation();
             } else {
                 return createGravityRelatedHeightToGeographic3D(
                     createSimilarPropertiesTransformation(self), sourceCRS(),
-                    targetCRS(), projFilename, coordinateOperationAccuracies());
+                    targetCRS(), interpolationCRS(), projFilename,
+                    coordinateOperationAccuracies());
             }
         }
     }
@@ -10972,19 +10988,19 @@ struct MyPROJStringExportableHorizVertical final
         // cppcheck-suppress functionStatic
         _exportToPROJString(io::PROJStringFormatter *formatter) const override {
 
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
         horizTransform->_exportToPROJString(formatter);
 
         formatter->startInversion();
         geogDst->addAngularUnitConvertAndAxisSwap(formatter);
         formatter->stopInversion();
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
 
         verticalTransform->_exportToPROJString(formatter);
 
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
         geogDst->addAngularUnitConvertAndAxisSwap(formatter);
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
     }
 };
 
@@ -11016,7 +11032,7 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
         // cppcheck-suppress functionStatic
         _exportToPROJString(io::PROJStringFormatter *formatter) const override {
 
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
 
         opSrcCRSToGeogCRS->_exportToPROJString(formatter);
 
@@ -11024,17 +11040,17 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
         interpolationGeogCRS->addAngularUnitConvertAndAxisSwap(formatter);
         formatter->stopInversion();
 
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
 
         verticalTransform->_exportToPROJString(formatter);
 
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
 
         interpolationGeogCRS->addAngularUnitConvertAndAxisSwap(formatter);
 
         opGeogCRStoDstCRS->_exportToPROJString(formatter);
 
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
     }
 };
 
@@ -12380,7 +12396,8 @@ CoordinateOperationFactory::Private::createOperations(
         const auto &componentsSrc = compoundSrc->componentReferenceSystems();
         if (!componentsSrc.empty()) {
             std::vector<CoordinateOperationNNPtr> horizTransforms;
-            if (componentsSrc[0]->extractGeographicCRS()) {
+            auto srcGeogCRS = componentsSrc[0]->extractGeographicCRS();
+            if (srcGeogCRS) {
                 horizTransforms =
                     createOperations(componentsSrc[0], targetCRS, context);
             }
@@ -12394,11 +12411,61 @@ CoordinateOperationFactory::Private::createOperations(
                 for (const auto &horizTransform : horizTransforms) {
                     for (const auto &verticalTransform : verticalTransforms) {
 
-                        auto op = createHorizVerticalPROJBased(
-                            sourceCRS, targetCRS, horizTransform,
-                            verticalTransform);
+                        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);
+                                }
+                            }
+                        }
+                        bool done = false;
+                        if (interpolationGeogCRS &&
+                            (interpolationGeogCRS->_isEquivalentTo(
+                                 srcGeogCRS.get(),
+                                 util::IComparable::Criterion::EQUIVALENT) ||
+                             interpolationGeogCRS->is2DPartOf3D(
+                                 NN_NO_CHECK(srcGeogCRS.get())))) {
+                            auto srcToInterp = 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));
+                            }
+                        }
+                        if (!done) {
+                            auto op = createHorizVerticalPROJBased(
+                                sourceCRS, targetCRS, horizTransform,
+                                verticalTransform);
 
-                        res.emplace_back(op);
+                            res.emplace_back(op);
+                        }
                     }
                 }
                 return res;
-- 
cgit v1.2.3