import/export PROJ strings from ISO19111 code: require/output +type=crs for CRS objects (refs #1214)

6 files changed, 123 insertions, 68 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index ba292262..710b59e3 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -27,6 +27,8 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
+#define FROM_PROJ_CPP
+
 #include <errno.h>
 #include <stddef.h>
 #include <stdio.h>
@@ -44,7 +46,9 @@
 
 #include "proj/common.hpp"
 #include "proj/coordinateoperation.hpp"
+#include "proj/internal/internal.hpp"
 
+using namespace NS_PROJ::internal;
 
 /* Initialize PJ_COORD struct */
 PJ_COORD proj_coord (double x, double y, double z, double t) {
@@ -726,6 +730,20 @@ int proj_context_get_use_proj4_init_rules(PJ_CONTEXT *ctx, int from_legacy_code_
     return from_legacy_code_path;
 }
 
+/** Adds a " +type=crs" suffix to a PROJ string (if it is a PROJ string) */
+std::string pj_add_type_crs_if_needed(const std::string& str)
+{
+    std::string ret(str);
+    if( (starts_with(str, "proj=") ||
+            starts_with(str, "+proj=") ||
+            starts_with(str, "+init=") ||
+            starts_with(str, "+title=")) &&
+        str.find("type=crs") == std::string::npos )
+    {
+        ret += " +type=crs";
+    }
+    return ret;
+}
 
 /*****************************************************************************/
 PJ  *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *source_crs, const char *target_crs, PJ_AREA *area) {
@@ -763,79 +781,89 @@ PJ  *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *source_crs, const char
     const char* const* optionsImportCRS =
         proj_context_get_use_proj4_init_rules(ctx, FALSE) ? optionsProj4Mode : nullptr;
 
-    auto src = proj_create_from_user_input(ctx, source_crs, optionsImportCRS);
-    if( !src ) {
-        proj_context_log_debug(ctx, "Cannot instantiate source_crs");
-        return nullptr;
-    }
+    try
+    {
+        std::string source_crs_modified(pj_add_type_crs_if_needed(source_crs));
+        std::string target_crs_modified(pj_add_type_crs_if_needed(target_crs));
 
-    auto dst = proj_create_from_user_input(ctx, target_crs, optionsImportCRS);
-    if( !dst ) {
-        proj_context_log_debug(ctx, "Cannot instantiate target_crs");
-        proj_destroy(src);
-        return nullptr;
-    }
+        auto src = proj_create_from_user_input(ctx, source_crs_modified.c_str(), optionsImportCRS);
+        if( !src ) {
+            proj_context_log_debug(ctx, "Cannot instantiate source_crs");
+            return nullptr;
+        }
 
-    auto operation_ctx = proj_create_operation_factory_context(ctx, nullptr);
-    if( !operation_ctx ) {
-        proj_destroy(src);
-        proj_destroy(dst);
-        return nullptr;
-    }
+        auto dst = proj_create_from_user_input(ctx, target_crs_modified.c_str(), optionsImportCRS);
+        if( !dst ) {
+            proj_context_log_debug(ctx, "Cannot instantiate target_crs");
+            proj_destroy(src);
+            return nullptr;
+        }
 
-    if( area && area->bbox_set ) {
-        proj_operation_factory_context_set_area_of_interest(
-                                            ctx,
-                                            operation_ctx,
-                                            area->west_lon_degree,
-                                            area->south_lat_degree,
-                                            area->east_lon_degree,
-                                            area->north_lat_degree);
-    }
+        auto operation_ctx = proj_create_operation_factory_context(ctx, nullptr);
+        if( !operation_ctx ) {
+            proj_destroy(src);
+            proj_destroy(dst);
+            return nullptr;
+        }
 
-    proj_operation_factory_context_set_grid_availability_use(
-        ctx, operation_ctx, PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
+        if( area && area->bbox_set ) {
+            proj_operation_factory_context_set_area_of_interest(
+                                                ctx,
+                                                operation_ctx,
+                                                area->west_lon_degree,
+                                                area->south_lat_degree,
+                                                area->east_lon_degree,
+                                                area->north_lat_degree);
+        }
 
-    auto op_list = proj_create_operations(ctx, src, dst, operation_ctx);
+        proj_operation_factory_context_set_grid_availability_use(
+            ctx, operation_ctx, PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
 
-    proj_operation_factory_context_destroy(operation_ctx);
-    proj_destroy(src);
-    proj_destroy(dst);
+        auto op_list = proj_create_operations(ctx, src, dst, operation_ctx);
 
-    if( !op_list ) {
-        return nullptr;
-    }
+        proj_operation_factory_context_destroy(operation_ctx);
+        proj_destroy(src);
+        proj_destroy(dst);
 
-    if( proj_list_get_count(op_list) == 0 ) {
+        if( !op_list ) {
+            return nullptr;
+        }
+
+        if( proj_list_get_count(op_list) == 0 ) {
+            proj_list_destroy(op_list);
+            proj_context_log_debug(ctx, "No operation found matching criteria");
+            return nullptr;
+        }
+
+        auto op = proj_list_get(ctx, op_list, 0);
         proj_list_destroy(op_list);
-        proj_context_log_debug(ctx, "No operation found matching criteria");
-        return nullptr;
-    }
+        if( !op ) {
+            return nullptr;
+        }
 
-    auto op = proj_list_get(ctx, op_list, 0);
-    proj_list_destroy(op_list);
-    if( !op ) {
-        return nullptr;
-    }
+        proj_string = proj_as_proj_string(ctx, op, PJ_PROJ_5, nullptr);
+        if( !proj_string) {
+            proj_destroy(op);
+            proj_context_log_debug(ctx, "Cannot export operation as a PROJ string");
+            return nullptr;
+        }
+
+        PJ* P;
+        if( proj_string[0] == '\0' ) {
+            /* Null transform ? */
+            P = proj_create(ctx, "proj=affine");
+        } else {
+            P = proj_create(ctx, proj_string);
+        }
 
-    proj_string = proj_as_proj_string(ctx, op, PJ_PROJ_5, nullptr);
-    if( !proj_string) {
         proj_destroy(op);
-        proj_context_log_debug(ctx, "Cannot export operation as a PROJ string");
-        return nullptr;
-    }
 
-    PJ* P;
-    if( proj_string[0] == '\0' ) {
-        /* Null transform ? */
-        P = proj_create(ctx, "proj=affine");
-    } else {
-        P = proj_create(ctx, proj_string);
+        return P;
+    }
+    catch( const std::exception& )
+    {
+        return nullptr;
     }
-
-    proj_destroy(op);
-
-    return P;
 }
 
 PJ *proj_destroy (PJ *P) {
diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp
index af53a051..bde4f813 100644
--- a/src/apps/cs2cs.cpp
+++ b/src/apps/cs2cs.cpp
@@ -213,8 +213,9 @@ static PJ *instanciate_crs(const std::string &definition,
                                const char *const *optionsImportCRS,
                                bool &isGeog, double &toRadians,
                                bool &isLatFirst) {
-    PJ *crs = proj_create_from_user_input(nullptr, definition.c_str(),
-                                                  optionsImportCRS);
+    PJ *crs = proj_create_from_user_input(nullptr,
+                                          pj_add_type_crs_if_needed(definition).c_str(),
+                                          optionsImportCRS);
     if (!crs) {
         return nullptr;
     }
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 81c70b7a..a003cd9c 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -1176,7 +1176,8 @@ void GeodeticCRS::_exportToPROJString(
 {
     const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_;
     if (!extensionProj4.empty()) {
-        formatter->ingestPROJString(extensionProj4);
+        formatter->ingestPROJString(
+            replaceAll(extensionProj4, " +type=crs", ""));
         formatter->addNoDefs(false);
         return;
     }
@@ -1907,7 +1908,8 @@ void GeographicCRS::_exportToPROJString(
 {
     const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_;
     if (!extensionProj4.empty()) {
-        formatter->ingestPROJString(extensionProj4);
+        formatter->ingestPROJString(
+            replaceAll(extensionProj4, " +type=crs", ""));
         formatter->addNoDefs(false);
         return;
     }
@@ -2719,7 +2721,8 @@ void ProjectedCRS::_exportToPROJString(
 {
     const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_;
     if (!extensionProj4.empty()) {
-        formatter->ingestPROJString(extensionProj4);
+        formatter->ingestPROJString(
+            replaceAll(extensionProj4, " +type=crs", ""));
         formatter->addNoDefs(false);
         return;
     }
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 3b6563f3..565f43e3 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -2086,7 +2086,8 @@ AuthorityFactory::createGeodeticCRS(const std::string &code,
 
         if (!text_definition.empty()) {
             DatabaseContext::Private::RecursionDetector detector(d->context());
-            auto obj = createFromUserInput(text_definition, d->context());
+            auto obj = createFromUserInput(
+                pj_add_type_crs_if_needed(text_definition), d->context());
             auto geodCRS = util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(obj);
             if (geodCRS) {
                 return cloneWithProps(NN_NO_CHECK(geodCRS), props);
@@ -2333,7 +2334,8 @@ AuthorityFactory::createProjectedCRS(const std::string &code) const {
 
         if (!text_definition.empty()) {
             DatabaseContext::Private::RecursionDetector detector(d->context());
-            auto obj = createFromUserInput(text_definition, d->context());
+            auto obj = createFromUserInput(
+                pj_add_type_crs_if_needed(text_definition), d->context());
             auto projCRS = dynamic_cast<const crs::ProjectedCRS *>(obj.get());
             if (projCRS) {
                 const auto &conv = projCRS->derivingConversionRef();
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 5be02ffe..d403495b 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -3344,6 +3344,9 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard(
                ci_equal(stripQuotes(extensionChildren[0]), "PROJ4")) {
         std::string projString = stripQuotes(extensionChildren[1]);
         if (starts_with(projString, "+proj=")) {
+            if (projString.find(" +type=crs") == std::string::npos) {
+                projString += " +type=crs";
+            }
             try {
                 auto projObj =
                     PROJStringParser().createFromPROJString(projString);
@@ -4619,6 +4622,9 @@ std::string IPROJStringExportable::exportToPROJString(
         }
     }
     if (bIsCRS) {
+        if (!formatter->hasParam("type")) {
+            formatter->addParam("type", "crs");
+        }
         formatter->setCRSExport(false);
     }
     return formatter->toString();
@@ -5303,6 +5309,9 @@ PROJStringSyntaxParser(const std::string &projString, std::vector<Step> &steps,
             vunits = word.substr(strlen("vunits="));
         } else if (starts_with(word, "vto_meter=")) {
             vto_meter = word.substr(strlen("vto_meter="));
+        } else if (word == "type=crs" &&
+                   (!vunits.empty() || !vto_meter.empty())) {
+            // ok
         } else if (starts_with(word, "title=")) {
             title = word.substr(strlen("title="));
             prevWasTitle = true;
@@ -6902,7 +6911,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                 param.key == "towgs84" || param.key == "nadgrids" ||
                 param.key == "geoidgrids" || param.key == "units" ||
                 param.key == "to_meter" || param.key == "vunits" ||
-                param.key == "vto_meter") {
+                param.key == "vto_meter" || param.key == "type") {
                 continue;
             }
             if (param.value.empty()) {
@@ -7252,6 +7261,10 @@ static const metadata::ExtentPtr &getExtent(const crs::CRS *crs) {
 // ---------------------------------------------------------------------------
 
 /** \brief Instantiate a sub-class of BaseObject from a PROJ string.
+ *
+ * The projString must contain +type=crs for the object to be detected as a
+ * CRS instead of a CoordinateOperation.
+ *
  * @throw ParsingException
  */
 BaseObjectNNPtr
@@ -7300,7 +7313,8 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
 
     // +init=xxxx:yyyy syntax
     if (d->steps_.size() == 1 && d->steps_[0].isInit &&
-        d->steps_[0].paramValues.size() == 0) {
+        (d->steps_[0].paramValues.size() == 0 ||
+         d->getParamValue(d->steps_[0], "type") == "crs")) {
 
         // Those used to come from a text init file
         // We only support them in compatibility mode
@@ -7483,6 +7497,11 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
         }
     }
 
+    if (d->steps_.size() == 1 && iHelmert < 0 && iMolodensky < 0 &&
+        d->getParamValue(d->steps_[0], "type") != "crs") {
+        unexpectedStructure = true;
+    }
+
     if (!unexpectedStructure) {
         if (iFirstGeogStep == 0 && iSecondGeogStep < 0 && iProjStep < 0 &&
             iHelmert < 0 && iFirstCart < 0 && iMolodensky < 0 &&
diff --git a/src/proj_internal.h b/src/proj_internal.h
index bbadb376..0cc99c81 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -899,6 +899,8 @@ void   pj_freeup_plain (PJ *P);
 
 PJ* pj_init_ctx_with_allow_init_epsg( projCtx_t *ctx, int argc, char **argv, int allow_init_epsg );
 
+std::string PROJ_DLL pj_add_type_crs_if_needed(const std::string& str);
+
 /* classic public API */
 #include "proj_api.h"