wip: support custom non-crs proj transforms again

closes: #4306
mapnik · Apr 18, 2022 · 10801d4 · 10801d4
1 parent 1ba1278
commit 10801d4
Show file tree

Hide file tree

Showing 2 changed files with 100 additions and 23 deletions.
diff --git a/src/proj_transform.cpp b/src/proj_transform.cpp
@@ -35,6 +35,7 @@
 #ifdef MAPNIK_USE_PROJ
 // proj
 #include <proj.h>
+#include <proj/io.hpp>
 #endif
 
 // stl
@@ -92,6 +93,29 @@ auto envelope_points(box2d<T> const& env, std::size_t num_points) -> geometry::m
 
 } // namespace
 
+// from src/4D_api.cpp
+/** 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 ((str.rfind("proj=", 0) == 0 || str.rfind("+proj=", 0) == 0 ||
+         str.rfind("+init=", 0) == 0 || str.rfind("+title=", 0) == 0) &&
+        str.find("type=crs") == std::string::npos) {
+        ret += " +type=crs";
+    }
+    return ret;
+}
+
+void throw_proj_exception(PJ_CONTEXT *ctx, std::string msg)
+{
+    throw std::runtime_error(msg
+#if MAPNIK_PROJ_VERSION >= 80000
+                             + " because of "
+                             + std::string(proj_context_errno_string(ctx, proj_context_errno(ctx))));
+#else
+                             );
+#endif
+}
+
 proj_transform::proj_transform(projection const& source, projection const& dest)
     : is_source_longlat_(false)
     , is_dest_longlat_(false)
@@ -125,30 +149,80 @@ proj_transform::proj_transform(projection const& source, projection const& dest)
 #ifdef MAPNIK_USE_PROJ
             ctx_ = proj_context_create();
             proj_log_level(ctx_, PJ_LOG_ERROR);
-            transform_ = proj_create_crs_to_crs(ctx_, source.params().c_str(), dest.params().c_str(), nullptr);
-            if (transform_ == nullptr)
-            {
-                throw std::runtime_error(std::string("Cannot initialize proj_transform  (crs_to_crs) for given projections: '") +
-                                         source.params() + "'->'" + dest.params() +
-#if MAPNIK_PROJ_VERSION >= 80000
-                                         "' because of " + std::string(proj_context_errno_string(ctx_, proj_context_errno(ctx_))));
-#else
-                                         "'");
-#endif
+            // we replicate what proj_create_crs_to_crs() does and then call
+            // proj_create_crs_to_crs_from_pj because we need a PJ object for
+            // proj_is_crs and we want to avoid running proj_create twice as
+            // as much (once ourselves and once in proj_create_crs_to_crs)
+            std::string crs_source = pj_add_type_crs_if_needed(source.params());
+            std::string crs_dest = pj_add_type_crs_if_needed(dest.params());
+            PJ* src;
+            PJ* dst;
+            try {
+                src = proj_create(ctx_, crs_source.c_str());
+            } catch( const std::exception& ) {
+                src = nullptr;
             }
-            PJ* transform_gis = proj_normalize_for_visualization(ctx_, transform_);
-            if (transform_gis == nullptr)
-            {
-                throw std::runtime_error(std::string("Cannot initialize proj_transform (normalize) for given projections: '") +
-                                         source.params() + "'->'" + dest.params() +
-#if MAPNIK_PROJ_VERSION >= 80000
-                                         "' because of " + std::string(proj_context_errno_string(ctx_, proj_context_errno(ctx_))));
-#else
-                                         "'");
-#endif
+            if (!src) {
+                throw_proj_exception(ctx_,
+                    std::string("Cannot instantiate source projection for '") + source.params() + "'");
+            }
+            try {
+                dst = proj_create(ctx_, crs_dest.c_str());
+            } catch( const std::exception& ) {
+                dst = nullptr;
+            }
+            if (!dst) {
+                throw_proj_exception(ctx_,
+                    std::string("Cannot instantiate dest projection for '") + dest.params() + "'");
+            }
+            if (proj_is_crs(src) && proj_is_crs(dst)) {
+                // if both source and destination are a CRS, we can use
+                // proj_create_crs_to_crs
+                transform_ = proj_create_crs_to_crs_from_pj(ctx_, src, dst, nullptr, nullptr);
+                if (transform_ == nullptr)
+                {
+                    throw_proj_exception(ctx_,
+                        std::string("Cannot initialize proj_transform (crs_to_crs) for given projections: '") +
+                        source.params() + "'->'" + dest.params() + "'");
+                }
+                PJ* transform_gis = proj_normalize_for_visualization(ctx_, transform_);
+                if (transform_gis == nullptr)
+                {
+                    throw_proj_exception(ctx_,
+                        std::string("Cannot initialize proj_transform (normalize) for given projections: '") +
+                                             source.params() + "'->'" + dest.params() + "'");
+                }
+                proj_destroy(transform_);
+                transform_ = transform_gis;
+            } else {
+                // if one of the projections is not a CRS, then this was
+                // probably a user-supplied transform and we have to create
+                // the transform ourselves
+                auto formatter = osgeo::proj::io::PROJStringFormatter::create();
+                formatter->startInversion();
+                try {
+                    formatter->ingestPROJString(source.params());
+                } catch (const osgeo::proj::io::ParsingException &e) {
+                    throw_proj_exception(ctx_, std::string("Failed to ingest source: '") + source.params() + "'");
+                }
+                formatter->stopInversion();
+                try {
+                    formatter->ingestPROJString(dest.params());
+                } catch (const osgeo::proj::io::ParsingException &e) {
+                    throw_proj_exception(ctx_, std::string("Failed to ingest dest: '") + dest.params() + "'");
+                }
+                try {
+                    transform_ = proj_create(ctx_, formatter->toString().c_str());
+                } catch( const std::exception& ) {
+                    transform_ = nullptr;
+                }
+                if (!transform_) {
+                    throw_proj_exception(ctx_,
+                        std::string("Cannot instantiate custom projection pipeline: '") + formatter->toString() + "'");
+                }
             }
-            proj_destroy(transform_);
-            transform_ = transform_gis;
+            proj_destroy(src);
+            proj_destroy(dst);
 #else
             throw std::runtime_error(
               std::string(

diff --git a/test/unit/projection/proj_transform.cpp b/test/unit/projection/proj_transform.cpp
@@ -88,13 +88,16 @@ TEST_CASE("projection transform")
         mapnik::projection proj_4937("epsg:4937");
         //"Webmercator" WGS 84 / Pseudo-Mercator -- Spherical Mercator, Google Maps, OpenStreetMap, Bing, ArcGIS, ESRI
         mapnik::projection proj_3857("epsg:3857");
+        //axisswap, which is not a CRS
+        mapnik::projection proj_axisswap("+proj=axisswap +order=2,1");
 
         mapnik::proj_transform tr1(proj_4236, proj_27700);
         mapnik::proj_transform tr2(proj_4236, proj_8857);
         mapnik::proj_transform tr3(proj_4236, proj_4236);
         mapnik::proj_transform tr4(proj_4236, proj_4937);
         mapnik::proj_transform tr5(proj_4236, proj_3857);
-        std::initializer_list<std::reference_wrapper<mapnik::proj_transform>> transforms = {tr1, tr2, tr3, tr4, tr5};
+        mapnik::proj_transform tr6(proj_4236, proj_axisswap);
+        std::initializer_list<std::reference_wrapper<mapnik::proj_transform>> transforms = {tr1, tr2, tr3, tr4, tr5, tr6};
 
         std::initializer_list<std::tuple<double, double>> coords = {
           {-4.0278869, 57.8796955}, // Dórnach, Highland