Merge remote-tracking branch 'osgeo/master'6.2.0

13 files changed, 232 insertions, 48 deletions

diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst
index e5f48439..c010fd13 100644
--- a/docs/source/development/reference/functions.rst
+++ b/docs/source/development/reference/functions.rst
@@ -696,6 +696,19 @@ Various
     :type `direction`: PJ_DIRECTION
     :returns: :c:type:`int` 1 if output units is expected in radians, otherwise 0
 
+
+Cleanup
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+.. c:function:: void proj_cleanup()
+
+    .. versionadded:: 6.2.0
+
+    This function frees global resources (grids, cache of +init files). It
+    should be called typically before process termination, and *after* having
+    freed PJ and PJ_CONTEXT objects.
+
+
 C API for ISO-19111 functionality
 +++++++++++++++++++++++++++++++++
 
diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index 9c99fc0e..5b1a8588 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -432,8 +432,8 @@ class PROJ_GCC_DLL PROJStringFormatter {
     PROJ_INTERNAL void popOmitZUnitConversion();
     PROJ_INTERNAL bool omitZUnitConversion() const;
 
-    PROJ_INTERNAL void setDropEarlyBindingsTerms(bool drop);
-    PROJ_INTERNAL bool getDropEarlyBindingsTerms() const;
+    PROJ_INTERNAL void setLegacyCRSToCRSContext(bool legacyContext);
+    PROJ_INTERNAL bool getLegacyCRSToCRSContext() const;
 
     PROJ_INTERNAL const DatabaseContextPtr &databaseContext() const;
 
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index 6df84099..71de2d0e 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -773,6 +773,7 @@ proj_as_projjson
 proj_as_proj_string
 proj_assign_context
 proj_as_wkt
+proj_cleanup
 proj_clone
 proj_concatoperation_get_step
 proj_concatoperation_get_step_count
diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp
index fcd9fa01..e6c7b2b1 100644
--- a/src/apply_gridshift.cpp
+++ b/src/apply_gridshift.cpp
@@ -344,7 +344,14 @@ PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction) {
     ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);
 
     if (ct == nullptr || ct->cvs == nullptr) {
-        pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
+        if( P->gridlist_count == 1 &&
+            strcmp(P->gridlist[0]->gridname, "null") == 0) {
+          // TODO: remove this particular case that is put there just to be
+          // able to handle longitudes outside of -180,180
+          out = lp;
+        } else {
+          pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
+        }
         return out;
     }
 
diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp
index 20e5e73c..3099c3e8 100644
--- a/src/apps/cs2cs.cpp
+++ b/src/apps/cs2cs.cpp
@@ -645,7 +645,7 @@ int main(int argc, char **argv) {
 
     proj_destroy(transformation);
 
-    pj_deallocate_grids();
+    proj_cleanup();
 
     exit(0); /* normal completion */
 }
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index dc700445..e9524163 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -12064,7 +12064,7 @@ CoordinateOperationFactory::Private::createOperations(
         }
         auto projFormatter = io::PROJStringFormatter::create();
         projFormatter->setCRSExport(true);
-        projFormatter->setDropEarlyBindingsTerms(true);
+        projFormatter->setLegacyCRSToCRSContext(true);
         projFormatter->startInversion();
         sourceProjExportable->_exportToPROJString(projFormatter.get());
         auto geogSrc =
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index b9694ba2..60b316f1 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -1308,7 +1308,7 @@ void GeodeticCRS::addDatumInfoToPROJString(
     const auto &nadgrids = formatter->getHDatumExtension();
     const auto &l_datum = datum();
     if (formatter->getCRSExport() && l_datum && TOWGS84Params.empty() &&
-        nadgrids.empty() && !formatter->getDropEarlyBindingsTerms()) {
+        nadgrids.empty()) {
         if (l_datum->_isEquivalentTo(
                 datum::GeodeticReferenceFrame::EPSG_6326.get(),
                 util::IComparable::Criterion::EQUIVALENT)) {
@@ -1323,7 +1323,12 @@ void GeodeticCRS::addDatumInfoToPROJString(
                        datum::GeodeticReferenceFrame::EPSG_6269.get(),
                        util::IComparable::Criterion::EQUIVALENT)) {
             datumWritten = true;
-            formatter->addParam("datum", "NAD83");
+            if (formatter->getLegacyCRSToCRSContext()) {
+                // We do not want datum=NAD83 to cause a useless towgs84=0,0,0
+                formatter->addParam("ellps", "GRS80");
+            } else {
+                formatter->addParam("datum", "NAD83");
+            }
         }
     }
     if (!datumWritten) {
@@ -2083,8 +2088,31 @@ void GeographicCRS::_exportToPROJString(
         primeMeridian()->longitude().getSIValue() != 0.0 ||
         !formatter->getTOWGS84Parameters().empty() ||
         !formatter->getHDatumExtension().empty()) {
+
         formatter->addStep("longlat");
-        addDatumInfoToPROJString(formatter);
+        bool done = false;
+        if (formatter->getLegacyCRSToCRSContext() &&
+            formatter->getHDatumExtension().empty() &&
+            formatter->getTOWGS84Parameters().empty()) {
+            const auto &l_datum = datum();
+            if (l_datum &&
+                l_datum->_isEquivalentTo(
+                    datum::GeodeticReferenceFrame::EPSG_6326.get(),
+                    util::IComparable::Criterion::EQUIVALENT)) {
+                done = true;
+                formatter->addParam("ellps", "WGS84");
+            } else if (l_datum &&
+                       l_datum->_isEquivalentTo(
+                           datum::GeodeticReferenceFrame::EPSG_6269.get(),
+                           util::IComparable::Criterion::EQUIVALENT)) {
+                done = true;
+                // We do not want datum=NAD83 to cause a useless towgs84=0,0,0
+                formatter->addParam("ellps", "GRS80");
+            }
+        }
+        if (!done) {
+            addDatumInfoToPROJString(formatter);
+        }
     }
     if (!formatter->getCRSExport()) {
         addAngularUnitConvertAndAxisSwap(formatter);
@@ -3132,7 +3160,8 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter,
                 formatter->addParam("units", projUnit);
             }
         }
-    } else if (formatter->getCRSExport()) {
+    } else if (formatter->getCRSExport() &&
+               !formatter->getLegacyCRSToCRSContext()) {
         formatter->addParam("units", "m");
     }
 
@@ -4143,11 +4172,6 @@ void BoundCRS::_exportToPROJString(
             "baseCRS of BoundCRS cannot be exported as a PROJ string");
     }
 
-    if (formatter->getDropEarlyBindingsTerms()) {
-        crs_exportable->_exportToPROJString(formatter);
-        return;
-    }
-
     auto vdatumProj4GridName = getVDatumPROJ4GRIDS();
     if (!vdatumProj4GridName.empty()) {
         formatter->setVDatumExtension(vdatumProj4GridName);
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index f4d93ef4..20a90946 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -6145,7 +6145,7 @@ struct PROJStringFormatter::Private {
     bool addNoDefs_ = true;
     bool coordOperationOptimizations_ = false;
     bool crsExport_ = false;
-    bool dropEarlyBindingsTerms_ = false;
+    bool legacyCRSToCRSContext_ = false;
 
     std::string result_{};
 
@@ -6723,7 +6723,7 @@ void PROJStringFormatter::Private::appendToResult(const char *str) {
 static void
 PROJStringSyntaxParser(const std::string &projString, std::vector<Step> &steps,
                        std::vector<Step::KeyValue> &globalParamValues,
-                       std::string &title, bool dropEarlyBindingsTerms) {
+                       std::string &title) {
     const char *c_str = projString.c_str();
     std::vector<std::string> tokens;
 
@@ -6811,27 +6811,14 @@ PROJStringSyntaxParser(const std::string &projString, std::vector<Step> &steps,
             } else if (word != "step") {
                 const auto pos = word.find('=');
                 auto key = word.substr(0, pos);
-                if (!(dropEarlyBindingsTerms &&
-                      (key == "towgs84" || key == "nadgrids" ||
-                       key == "geoidgrids"))) {
-                    auto pair = (pos != std::string::npos)
-                                    ? Step::KeyValue(key, word.substr(pos + 1))
-                                    : Step::KeyValue(key);
-                    if (dropEarlyBindingsTerms && key == "datum") {
-                        const auto datums = pj_get_datums_ref();
-                        for (int i = 0; datums[i].id != nullptr; i++) {
-                            if (pair.value == datums[i].id) {
-                                pair.key = "ellps";
-                                pair.value = datums[i].ellipse_id;
-                                break;
-                            }
-                        }
-                    }
-                    if (steps.empty()) {
-                        globalParamValues.push_back(pair);
-                    } else {
-                        steps.back().paramValues.push_back(pair);
-                    }
+
+                auto pair = (pos != std::string::npos)
+                                ? Step::KeyValue(key, word.substr(pos + 1))
+                                : Step::KeyValue(key);
+                if (steps.empty()) {
+                    globalParamValues.push_back(pair);
+                } else {
+                    steps.back().paramValues.push_back(pair);
                 }
             }
         }
@@ -6906,8 +6893,7 @@ void PROJStringFormatter::ingestPROJString(
 {
     std::vector<Step> steps;
     std::string title;
-    PROJStringSyntaxParser(str, steps, d->globalParamValues_, title,
-                           d->dropEarlyBindingsTerms_);
+    PROJStringSyntaxParser(str, steps, d->globalParamValues_, title);
     d->steps_.insert(d->steps_.end(), steps.begin(), steps.end());
 }
 
@@ -7175,14 +7161,14 @@ bool PROJStringFormatter::omitZUnitConversion() const {
 
 // ---------------------------------------------------------------------------
 
-void PROJStringFormatter::setDropEarlyBindingsTerms(bool drop) {
-    d->dropEarlyBindingsTerms_ = drop;
+void PROJStringFormatter::setLegacyCRSToCRSContext(bool legacyContext) {
+    d->legacyCRSToCRSContext_ = legacyContext;
 }
 
 // ---------------------------------------------------------------------------
 
-bool PROJStringFormatter::getDropEarlyBindingsTerms() const {
-    return d->dropEarlyBindingsTerms_;
+bool PROJStringFormatter::getLegacyCRSToCRSContext() const {
+    return d->legacyCRSToCRSContext_;
 }
 
 // ---------------------------------------------------------------------------
@@ -8269,6 +8255,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
     }
 
     auto axisType = AxisType::REGULAR;
+    bool bWebMercator = false;
 
     if (step.name == "tmerc" &&
         ((getParamValue(step, "axis") == "wsu" && iAxisSwap < 0) ||
@@ -8351,8 +8338,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                 }
             }
             if (getNumericValue(getParamValue(step, "a")) == 6378137) {
-                return createPseudoMercator(PropertyMap().set(
-                    IdentifiedObject::NAME_KEY, "WGS 84 / Pseudo-Mercator"));
+                bWebMercator = true;
             }
         } else if (hasParamValue(step, "lat_ts")) {
             mapping = getMapping(EPSG_CODE_METHOD_MERCATOR_VARIANT_B);
@@ -8613,7 +8599,11 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
 
     props.set("IMPLICIT_CS", true);
 
-    CRSNNPtr crs = ProjectedCRS::create(props, geogCRS, NN_NO_CHECK(conv), cs);
+    CRSNNPtr crs =
+        bWebMercator
+            ? createPseudoMercator(props.set(IdentifiedObject::NAME_KEY,
+                                             "WGS 84 / Pseudo-Mercator"))
+            : ProjectedCRS::create(props, geogCRS, NN_NO_CHECK(conv), cs);
 
     if (!hasParamValue(step, "geoidgrids") &&
         (hasParamValue(step, "vunits") || hasParamValue(step, "vto_meter"))) {
@@ -8671,7 +8661,7 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
     d->globalParamValues_.clear();
     d->projString_ = projString;
     PROJStringSyntaxParser(projString, d->steps_, d->globalParamValues_,
-                           d->title_, false);
+                           d->title_);
 
     if (d->steps_.empty()) {
         const auto &vunits = d->getGlobalParamValue("vunits");
diff --git a/src/malloc.cpp b/src/malloc.cpp
index 817f8f20..393437e3 100644
--- a/src/malloc.cpp
+++ b/src/malloc.cpp
@@ -255,3 +255,10 @@ PJ *pj_default_destructor (PJ *P, int errlev) {   /* Destructor */
     delete P;
     return nullptr;
 }
+
+/*****************************************************************************/
+void proj_cleanup() {
+/*****************************************************************************/
+  pj_clear_initcache();
+  pj_deallocate_grids();
+}
diff --git a/src/proj.h b/src/proj.h
index 57fb3940..14d315b5 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -458,6 +458,8 @@ double PROJ_DLL proj_todeg (double angle_in_radians);
 double PROJ_DLL proj_dmstor(const char *is, char **rs);
 char PROJ_DLL * proj_rtodms(char *s, double r, int pos, int neg);
 
+void PROJ_DLL proj_cleanup(void);
+
 /*! @endcond */
 
 /* ------------------------------------------------------------------------- */
diff --git a/test/cli/testvarious b/test/cli/testvarious
index 836158ca..6c00ba13 100755
--- a/test/cli/testvarious
+++ b/test/cli/testvarious
@@ -985,6 +985,18 @@ $EXE -f %.7f -I +proj=longlat +ellps=GRS80 +to +proj=ob_tran +o_proj=longlat +lo
 -122 46
 EOF
 
+echo "##############################################################" >> ${OUT}
+echo "Check +init=epsg:4326 +over +to +init=epsg:3857 +over" >> ${OUT}
+$EXE -f %.7f +init=epsg:4326 +over +to +init=epsg:3857 +over -E >>${OUT} <<EOF
+-181 49
+EOF
+
+echo "##############################################################" >> ${OUT}
+echo "Check +proj=longlat +over +datum=WGS84 +to proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +over" >> ${OUT}
+$EXE -f %.7f +proj=longlat +over +datum=WGS84 +to proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +over -E >>${OUT} <<EOF
+-181 49
+EOF
+
 
 # Done!
 # do 'diff' with distribution results
diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
index 13987332..8720a124 100644
--- a/test/cli/tv_out.dist
+++ b/test/cli/tv_out.dist
@@ -473,3 +473,9 @@ Check ob_tran with o_proj=longlat (#1525)
 -122 46	-132.0000000	46.0000000 0.0000000
 Test inverse handling
 -122 46	-112.0000000	46.0000000 0.0000000
+##############################################################
+Check +init=epsg:4326 +over +to +init=epsg:3857 +over
+-181 49	-20148827.8335825	6274861.3940066 0.0000000
+##############################################################
+Check +proj=longlat +over +datum=WGS84 +to proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +over
+-181 49	-20148827.8335825	6274861.3940066 0.0000000
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 1a9b69d5..89edfca4 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -7059,7 +7059,129 @@ TEST(operation, createOperation_fallback_to_proj4_strings) {
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +proj=axisswap +order=2,1 "
               "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
-              "+step +proj=longlat +geoc +ellps=WGS84 "
+              "+step +proj=longlat +geoc +datum=WGS84 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(
+    operation,
+    createOperation_fallback_to_proj4_strings_regular_with_datum_to_projliteral) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=utm +zone=11 +datum=NAD27 +type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +datum=WGS84 +over +type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +inv +proj=utm +zone=11 +datum=NAD27 "
+              "+step +proj=longlat +datum=WGS84 +over "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation,
+     createOperation_fallback_to_proj4_strings_proj_NAD83_to_projliteral) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=utm +zone=11 +datum=NAD83 +type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +datum=WGS84 +over +type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +inv +proj=utm +zone=11 +ellps=GRS80 "
+              "+step +proj=longlat +datum=WGS84 +over "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation,
+     createOperation_fallback_to_proj4_strings_geog_NAD83_to_projliteral) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=longlat +datum=NAD83 +type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +datum=WGS84 +over +type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=longlat +datum=WGS84 +over "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(
+    operation,
+    createOperation_fallback_to_proj4_strings_regular_with_nadgrids_to_projliteral) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=utm +zone=11 +ellps=clrk66 +nadgrids=@conus +type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +datum=WGS84 +over +type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +inv +proj=utm +zone=11 +ellps=clrk66 +nadgrids=@conus "
+              "+step +proj=longlat +datum=WGS84 +over "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation,
+     createOperation_fallback_to_proj4_strings_projliteral_to_projliteral) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=utm +zone=11 +datum=NAD27 +over +type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +datum=WGS84 +over +type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +inv +proj=utm +zone=11 +datum=NAD27 +over "
+              "+step +proj=longlat +datum=WGS84 +over "
               "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
 }