Merge pull request #1182 from rouault/plug_new_code

Remove data/epsg, IGNF and esri.* files / support legacy +init=epsg:XXXX syntax

1 files changed, 315 insertions, 47 deletions

diff --git a/src/coordinateoperation.cpp b/src/coordinateoperation.cpp
index 42cc806e..bbe84698 100644
--- a/src/coordinateoperation.cpp
+++ b/src/coordinateoperation.cpp
@@ -149,6 +149,8 @@ static std::set<std::string> buildSetEquivalentParameters() {
         {"latitude_of_point_2", "Latitude_Of_2nd_Point", nullptr},
         {"longitude_of_point_2", "Longitude_Of_2nd_Point", nullptr},
 
+        {"satellite_height", "height", nullptr},
+
         {EPSG_NAME_PARAMETER_FALSE_EASTING,
          EPSG_NAME_PARAMETER_EASTING_FALSE_ORIGIN,
          EPSG_NAME_PARAMETER_EASTING_PROJECTION_CENTRE, nullptr},
@@ -160,7 +162,8 @@ static std::set<std::string> buildSetEquivalentParameters() {
         {WKT1_LATITUDE_OF_ORIGIN, WKT1_LATITUDE_OF_CENTER,
          EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN,
          EPSG_NAME_PARAMETER_LATITUDE_FALSE_ORIGIN,
-         EPSG_NAME_PARAMETER_LATITUDE_PROJECTION_CENTRE, nullptr},
+         EPSG_NAME_PARAMETER_LATITUDE_PROJECTION_CENTRE, "Central_Parallel",
+         nullptr},
 
         {WKT1_CENTRAL_MERIDIAN, WKT1_LONGITUDE_OF_CENTER,
          EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN,
@@ -1898,6 +1901,58 @@ ConversionNNPtr Conversion::shallowClone() const {
 
 // ---------------------------------------------------------------------------
 
+//! @cond Doxygen_Suppress
+ConversionNNPtr
+Conversion::alterParametersLinearUnit(const common::UnitOfMeasure &unit,
+                                      bool convertToNewUnit) const {
+
+    std::vector<GeneralParameterValueNNPtr> newValues;
+    bool changesDone = false;
+    for (const auto &genOpParamvalue : parameterValues()) {
+        bool updated = false;
+        auto opParamvalue = dynamic_cast<const OperationParameterValue *>(
+            genOpParamvalue.get());
+        if (opParamvalue) {
+            const auto &paramValue = opParamvalue->parameterValue();
+            if (paramValue->type() == ParameterValue::Type::MEASURE) {
+                const auto &measure = paramValue->value();
+                if (measure.unit().type() ==
+                    common::UnitOfMeasure::Type::LINEAR) {
+                    if (!measure.unit()._isEquivalentTo(
+                            unit, util::IComparable::Criterion::EQUIVALENT)) {
+                        const double newValue =
+                            convertToNewUnit ? measure.convertToUnit(unit)
+                                             : measure.value();
+                        newValues.emplace_back(OperationParameterValue::create(
+                            opParamvalue->parameter(),
+                            ParameterValue::create(
+                                common::Measure(newValue, unit))));
+                        updated = true;
+                    }
+                }
+            }
+        }
+        if (updated) {
+            changesDone = true;
+        } else {
+            newValues.emplace_back(genOpParamvalue);
+        }
+    }
+    if (changesDone) {
+        auto conv = create(util::PropertyMap().set(
+                               common::IdentifiedObject::NAME_KEY, "unknown"),
+                           method(), newValues);
+        conv->setCRSs(this, false);
+        return conv;
+    } else {
+        return NN_NO_CHECK(
+            util::nn_dynamic_pointer_cast<Conversion>(shared_from_this()));
+    }
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
 /** \brief Instanciate a Conversion from a vector of GeneralParameterValue.
  *
  * @param properties See \ref general_properties. At minimum the name should be
@@ -5471,6 +5526,12 @@ Transformation::~Transformation() = default;
 
 // ---------------------------------------------------------------------------
 
+Transformation::Transformation(const Transformation &other)
+    : CoordinateOperation(other), SingleOperation(other),
+      d(internal::make_unique<Private>(*other.d)) {}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Return the source crs::CRS of the transformation.
  *
  * @return the source CRS.
@@ -5492,6 +5553,17 @@ const crs::CRSNNPtr &Transformation::targetCRS() PROJ_CONST_DEFN {
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+TransformationNNPtr Transformation::shallowClone() const {
+    auto conv = Transformation::nn_make_shared<Transformation>(*this);
+    conv->assignSelf(conv);
+    conv->setCRSs(this, false);
+    return conv;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 /** \brief Return the TOWGS84 parameters of the transformation.
  *
  * If this transformation uses Coordinate Frame Rotation, Position Vector
@@ -7595,6 +7667,8 @@ void Transformation::_exportToPROJString(
             "Transformation cannot be exported as a PROJ.4 string");
     }
 
+    formatter->setCoordinateOperationOptimizations(true);
+
     bool positionVectorConvention = true;
     bool sevenParamsTransform = false;
     bool threeParamsTransform = false;
@@ -8220,6 +8294,32 @@ bool SingleOperation::exportToPROJStringGeneric(
     if (isAxisOrderReversal(methodEPSGCode)) {
         formatter->addStep("axisswap");
         formatter->addParam("order", "2,1");
+        auto sourceCRSGeog =
+            dynamic_cast<const crs::GeographicCRS *>(sourceCRS().get());
+        auto targetCRSGeog =
+            dynamic_cast<const crs::GeographicCRS *>(targetCRS().get());
+        if (sourceCRSGeog && targetCRSGeog) {
+            const auto &unitSrc =
+                sourceCRSGeog->coordinateSystem()->axisList()[0]->unit();
+            const auto &unitDst =
+                targetCRSGeog->coordinateSystem()->axisList()[0]->unit();
+            if (!unitSrc._isEquivalentTo(
+                    unitDst, util::IComparable::Criterion::EQUIVALENT)) {
+                formatter->addStep("unitconvert");
+                auto projUnit = unitSrc.exportToPROJString();
+                if (projUnit.empty()) {
+                    formatter->addParam("xy_in", unitSrc.conversionToSI());
+                } else {
+                    formatter->addParam("xy_in", projUnit);
+                }
+                projUnit = unitDst.exportToPROJString();
+                if (projUnit.empty()) {
+                    formatter->addParam("xy_out", unitDst.conversionToSI());
+                } else {
+                    formatter->addParam("xy_out", projUnit);
+                }
+            }
+        }
         return true;
     }
 
@@ -9228,6 +9328,7 @@ struct FilterAndSort {
         // Precompute a number of parameters for each operation that will be
         // useful for the sorting.
         std::map<CoordinateOperation *, PrecomputedOpCharacteristics> map;
+        const auto gridAvailabilityUse = context->getGridAvailabilityUse();
         for (const auto &op : res) {
             bool dummy = false;
             auto extentOp = getExtent(op, true, dummy);
@@ -9260,14 +9361,20 @@ struct FilterAndSort {
             bool gridsAvailable = true;
             bool gridsKnown = true;
             if (context->getAuthorityFactory() &&
-                context->getGridAvailabilityUse() ==
-                    CoordinateOperationContext::GridAvailabilityUse::
-                        USE_FOR_SORTING) {
+                (gridAvailabilityUse ==
+                     CoordinateOperationContext::GridAvailabilityUse::
+                         USE_FOR_SORTING ||
+                 gridAvailabilityUse ==
+                     CoordinateOperationContext::GridAvailabilityUse::
+                         IGNORE_GRID_AVAILABILITY)) {
                 const auto gridsNeeded = op->gridsNeeded(
                     context->getAuthorityFactory()->databaseContext());
                 for (const auto &gridDesc : gridsNeeded) {
                     hasGrids = true;
-                    if (!gridDesc.available) {
+                    if (gridAvailabilityUse ==
+                            CoordinateOperationContext::GridAvailabilityUse::
+                                USE_FOR_SORTING &&
+                        !gridDesc.available) {
                         gridsAvailable = false;
                     }
                     if (gridDesc.packageName.empty()) {
@@ -10130,6 +10237,72 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
     const auto candidatesDstGeod(
         findCandidateGeodCRSForDatum(authFactory, geodDst->datum()));
 
+    auto createTransformations = [&](const crs::CRSNNPtr &candidateSrcGeod,
+                                     const crs::CRSNNPtr &candidateDstGeod,
+                                     const CoordinateOperationNNPtr &opFirst,
+                                     bool isNullFirst) {
+        const auto opsSecond =
+            createOperations(candidateSrcGeod, candidateDstGeod, context);
+        const auto opsThird =
+            createOperations(candidateDstGeod, targetCRS, context);
+        assert(!opsThird.empty());
+
+        for (auto &opSecond : opsSecond) {
+            // Check that it is not a transformation synthetized by
+            // ourselves
+            if (!hasIdentifiers(opSecond)) {
+                continue;
+            }
+            // And even if it is a referenced transformation, check that
+            // it is not a trivial one
+            auto so = dynamic_cast<const SingleOperation *>(opSecond.get());
+            if (so && isAxisOrderReversal(so->method()->getEPSGCode())) {
+                continue;
+            }
+
+            std::vector<CoordinateOperationNNPtr> subOps;
+            if (isNullFirst) {
+                opSecond->setCRSs(
+                    sourceCRS, NN_CHECK_ASSERT(opSecond->targetCRS()), nullptr);
+            } else {
+                subOps.emplace_back(opFirst);
+            }
+            if (isNullTransformation(opsThird[0]->nameStr())) {
+                opSecond->setCRSs(NN_CHECK_ASSERT(opSecond->sourceCRS()),
+                                  targetCRS, nullptr);
+                subOps.emplace_back(opSecond);
+            } else {
+                subOps.emplace_back(opSecond);
+                subOps.emplace_back(opsThird[0]);
+            }
+            res.emplace_back(ConcatenatedOperation::createComputeMetadata(
+                subOps, !allowEmptyIntersection));
+        }
+    };
+
+    // Start in priority with candidates that have exactly the same name as
+    // the sourcCRS and targetCRS. Typically for the case of init=IGNF:XXXX
+    for (const auto &candidateSrcGeod : candidatesSrcGeod) {
+        if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
+            for (const auto &candidateDstGeod : candidatesDstGeod) {
+                if (candidateDstGeod->nameStr() == targetCRS->nameStr()) {
+                    const auto opsFirst =
+                        createOperations(sourceCRS, candidateSrcGeod, context);
+                    assert(!opsFirst.empty());
+                    const bool isNullFirst =
+                        isNullTransformation(opsFirst[0]->nameStr());
+                    createTransformations(candidateSrcGeod, candidateDstGeod,
+                                          opsFirst[0], isNullFirst);
+                    if (!res.empty()) {
+                        return;
+                    }
+                    break;
+                }
+            }
+            break;
+        }
+    }
+
     for (const auto &candidateSrcGeod : candidatesSrcGeod) {
         const auto opsFirst =
             createOperations(sourceCRS, candidateSrcGeod, context);
@@ -10137,44 +10310,8 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
         const bool isNullFirst = isNullTransformation(opsFirst[0]->nameStr());
 
         for (const auto &candidateDstGeod : candidatesDstGeod) {
-            const auto opsSecond =
-                createOperations(candidateSrcGeod, candidateDstGeod, context);
-            const auto opsThird =
-                createOperations(candidateDstGeod, targetCRS, context);
-            assert(!opsThird.empty());
-
-            for (auto &opSecond : opsSecond) {
-                // Check that it is not a transformation synthetized by
-                // ourselves
-                if (!hasIdentifiers(opSecond)) {
-                    continue;
-                }
-                // And even if it is a referenced transformation, check that
-                // it is not a trivial one
-                auto so = dynamic_cast<const SingleOperation *>(opSecond.get());
-                if (so && isAxisOrderReversal(so->method()->getEPSGCode())) {
-                    continue;
-                }
-
-                std::vector<CoordinateOperationNNPtr> subOps;
-                if (isNullFirst) {
-                    opSecond->setCRSs(sourceCRS,
-                                      NN_CHECK_ASSERT(opSecond->targetCRS()),
-                                      nullptr);
-                } else {
-                    subOps.emplace_back(opsFirst[0]);
-                }
-                if (isNullTransformation(opsThird[0]->nameStr())) {
-                    opSecond->setCRSs(NN_CHECK_ASSERT(opSecond->sourceCRS()),
-                                      targetCRS, nullptr);
-                    subOps.emplace_back(opSecond);
-                } else {
-                    subOps.emplace_back(opSecond);
-                    subOps.emplace_back(opsThird[0]);
-                }
-                res.emplace_back(ConcatenatedOperation::createComputeMetadata(
-                    subOps, !allowEmptyIntersection));
-            }
+            createTransformations(candidateSrcGeod, candidateDstGeod,
+                                  opsFirst[0], isNullFirst);
         }
         if (!res.empty()) {
             return;
@@ -10235,6 +10372,56 @@ CoordinateOperationFactory::Private::createOperations(
     std::vector<CoordinateOperationNNPtr> res;
     const bool allowEmptyIntersection = true;
 
+    const auto &sourceProj4Ext = sourceCRS->getExtensionProj4();
+    const auto &targetProj4Ext = targetCRS->getExtensionProj4();
+    if (!sourceProj4Ext.empty() || !targetProj4Ext.empty()) {
+
+        auto sourceProjExportable =
+            dynamic_cast<io::IPROJStringExportable *>(sourceCRS.get());
+        auto targetProjExportable =
+            dynamic_cast<io::IPROJStringExportable *>(targetCRS.get());
+        if (!sourceProjExportable) {
+            throw InvalidOperation("Source CRS is not PROJ exportable");
+        }
+        if (!targetProjExportable) {
+            throw InvalidOperation("Target CRS is not PROJ exportable");
+        }
+        auto projFormatter = io::PROJStringFormatter::create(
+            io::PROJStringFormatter::Convention::PROJ_4);
+        projFormatter->startInversion();
+        sourceProjExportable->_exportToPROJString(projFormatter.get());
+
+        auto geogSrc =
+            dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
+        if (geogSrc) {
+            auto proj5Formatter = io::PROJStringFormatter::create(
+                io::PROJStringFormatter::Convention::PROJ_5);
+            geogSrc->addAngularUnitConvertAndAxisSwap(proj5Formatter.get());
+            projFormatter->ingestPROJString(proj5Formatter->toString());
+        }
+
+        projFormatter->stopInversion();
+
+        targetProjExportable->_exportToPROJString(projFormatter.get());
+
+        auto geogDst =
+            dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
+        if (geogDst) {
+            auto proj5Formatter = io::PROJStringFormatter::create(
+                io::PROJStringFormatter::Convention::PROJ_5);
+            geogDst->addAngularUnitConvertAndAxisSwap(proj5Formatter.get());
+            projFormatter->ingestPROJString(proj5Formatter->toString());
+        }
+
+        const auto PROJString = projFormatter->toString();
+        auto properties = util::PropertyMap().set(
+            common::IdentifiedObject::NAME_KEY,
+            buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()));
+        res.emplace_back(SingleOperation::createPROJBased(
+            properties, PROJString, sourceCRS, targetCRS, {}));
+        return res;
+    }
+
     auto geodSrc = dynamic_cast<const crs::GeodeticCRS *>(sourceCRS.get());
     auto geodDst = dynamic_cast<const crs::GeodeticCRS *>(targetCRS.get());
 
@@ -10291,7 +10478,9 @@ CoordinateOperationFactory::Private::createOperations(
             const auto &dstDatum = geodDst->datum();
             const bool dstHasDatumWithId =
                 dstDatum && !dstDatum->identifiers().empty();
-            if (srcHasDatumWithId && dstHasDatumWithId) {
+            if (srcHasDatumWithId && dstHasDatumWithId &&
+                !srcDatum->_isEquivalentTo(
+                    dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) {
                 createOperationsWithDatumPivot(res, sourceCRS, targetCRS,
                                                geodSrc, geodDst, context);
                 doFilterAndCheckPerfectOp = !res.empty();
@@ -10417,11 +10606,10 @@ CoordinateOperationFactory::Private::createOperations(
             dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
         auto geogCRSOfBaseOfBoundSrc =
             boundSrc->baseCRS()->extractGeographicCRS();
-        if (hubSrcGeog &&
+        if (hubSrcGeog && geogCRSOfBaseOfBoundSrc &&
             (hubSrcGeog->_isEquivalentTo(
                  geogDst, util::IComparable::Criterion::EQUIVALENT) ||
-             hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst))) &&
-            geogCRSOfBaseOfBoundSrc) {
+             hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst)))) {
             if (boundSrc->baseCRS() == geogCRSOfBaseOfBoundSrc) {
                 // Optimization to avoid creating a useless concatenated
                 // operation
@@ -10445,6 +10633,86 @@ CoordinateOperationFactory::Private::createOperations(
                     return res;
                 }
             }
+            // If the datum are equivalent, this is also fine
+        } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
+                   geogDst->datum() &&
+                   hubSrcGeog->datum()->_isEquivalentTo(
+                       geogDst->datum().get(),
+                       util::IComparable::Criterion::EQUIVALENT)) {
+            auto opsFirst =
+                createOperations(boundSrc->baseCRS(),
+                                 NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
+            auto opsLast = createOperations(hubSrc, targetCRS, context);
+            if (!opsFirst.empty() && !opsLast.empty()) {
+                for (const auto &opFirst : opsFirst) {
+                    for (const auto &opLast : opsLast) {
+                        try {
+                            res.emplace_back(
+                                ConcatenatedOperation::createComputeMetadata(
+                                    {opFirst, boundSrc->transformation(),
+                                     opLast},
+                                    !allowEmptyIntersection));
+                        } catch (const InvalidOperationEmptyIntersection &) {
+                        }
+                    }
+                }
+                if (!res.empty()) {
+                    return res;
+                }
+            }
+            // Consider WGS 84 and NAD83 as equivalent in that context if the
+            // geogCRSOfBaseOfBoundSrc ellipsoid is Clarke66 (for NAD27)
+            // Case of "+proj=latlong +ellps=clrk66
+            // +nadgrids=ntv1_can.dat,conus"
+            // to "+proj=latlong +datum=NAD83"
+        } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
+                   geogDst->datum() &&
+                   geogCRSOfBaseOfBoundSrc->ellipsoid()->_isEquivalentTo(
+                       datum::Ellipsoid::CLARKE_1866.get(),
+                       util::IComparable::Criterion::EQUIVALENT) &&
+                   hubSrcGeog->datum()->_isEquivalentTo(
+                       datum::GeodeticReferenceFrame::EPSG_6326.get(),
+                       util::IComparable::Criterion::EQUIVALENT) &&
+                   geogDst->datum()->_isEquivalentTo(
+                       datum::GeodeticReferenceFrame::EPSG_6269.get(),
+                       util::IComparable::Criterion::EQUIVALENT)) {
+            auto nnGeogCRSOfBaseOfBoundSrc =
+                NN_NO_CHECK(geogCRSOfBaseOfBoundSrc);
+            if (boundSrc->baseCRS()->_isEquivalentTo(
+                    nnGeogCRSOfBaseOfBoundSrc.get(),
+                    util::IComparable::Criterion::EQUIVALENT)) {
+                auto transf = boundSrc->transformation()->shallowClone();
+                transf->setProperties(util::PropertyMap().set(
+                    common::IdentifiedObject::NAME_KEY,
+                    buildTransfName(boundSrc->baseCRS()->nameStr(),
+                                    targetCRS->nameStr())));
+                transf->setCRSs(boundSrc->baseCRS(), targetCRS, nullptr);
+                res.emplace_back(transf);
+                return res;
+            } else {
+                auto opsFirst = createOperations(
+                    boundSrc->baseCRS(), nnGeogCRSOfBaseOfBoundSrc, context);
+                auto transf = boundSrc->transformation()->shallowClone();
+                transf->setProperties(util::PropertyMap().set(
+                    common::IdentifiedObject::NAME_KEY,
+                    buildTransfName(nnGeogCRSOfBaseOfBoundSrc->nameStr(),
+                                    targetCRS->nameStr())));
+                transf->setCRSs(nnGeogCRSOfBaseOfBoundSrc, targetCRS, nullptr);
+                if (!opsFirst.empty()) {
+                    for (const auto &opFirst : opsFirst) {
+                        try {
+                            res.emplace_back(
+                                ConcatenatedOperation::createComputeMetadata(
+                                    {opFirst, transf},
+                                    !allowEmptyIntersection));
+                        } catch (const InvalidOperationEmptyIntersection &) {
+                        }
+                    }
+                    if (!res.empty()) {
+                        return res;
+                    }
+                }
+            }
         }
 
         if (hubSrcGeog &&