Merge pull request #4468 from mapnik/bounding-box-split

Split bounding box along (0, 0) before reprojecting  [WIP] [skip ci]

include/mapnik/projection.hpp

@@ -26,7 +26,7 @@
 // mapnik
 #include <mapnik/config.hpp>
 #include <mapnik/well_known_srs.hpp>
-
+#include <mapnik/geometry/box2d.hpp>
 #include <mapnik/warning.hpp>
 MAPNIK_DISABLE_WARNING_PUSH
 #include <mapnik/warning_ignore.hpp>
@@ -79,7 +79,8 @@ class MAPNIK_DECL projection
     void inverse(double& x, double& y) const;
     std::string definition() const;
     std::string description() const;
-    void init_proj() const;
+    void init_proj();
+    std::optional<box2d<double>> area_of_use() const;
 
   private:
     void swap(projection& rhs);
@@ -87,7 +88,8 @@ class MAPNIK_DECL projection
   private:
     std::string params_;
     bool defer_proj_init_;
-    mutable bool is_geographic_;
+    bool is_geographic_;
+    std::optional<box2d<double>> area_of_use_;
     mutable PJ* proj_;
     mutable PJ_CONTEXT* proj_ctx_;
 };

src/proj_transform.cpp

@@ -42,7 +42,6 @@
 #include <stdexcept>
 
 namespace mapnik {
-
 namespace { // (local)
 
 // Returns points in clockwise order. This allows us to do anti-meridian checks.
@@ -90,6 +89,36 @@ auto envelope_points(box2d<T> const& env, std::size_t num_points) -> geometry::m
     return coords;
 }
 
+std::vector<box2d<double>> split_bounding_box(box2d<double> const& bbox)
+{
+    std::vector<box2d<double>> output;
+    if (bbox.minx() < 0 && bbox.maxx() > 0)
+    {
+        if (bbox.miny() < 0 && bbox.maxy() > 0)
+        {
+            output.emplace_back(bbox.minx(), bbox.miny(), 0, 0);
+            output.emplace_back(0, bbox.miny(), bbox.maxx(), 0);
+            output.emplace_back(0, 0, bbox.maxx(), bbox.maxy());
+            output.emplace_back(bbox.minx(), 0, 0, bbox.maxy());
+        }
+        else
+        {
+            output.emplace_back(bbox.minx(), bbox.miny(), 0, bbox.maxy());
+            output.emplace_back(0, bbox.miny(), bbox.maxx(), bbox.maxy());
+        }
+    }
+    else if (bbox.miny() < 0 && bbox.maxy() > 0)
+    {
+        output.emplace_back(bbox.minx(), bbox.miny(), bbox.maxx(), 0);
+        output.emplace_back(bbox.minx(), 0, bbox.maxx(), bbox.maxy());
+    }
+    else
+    {
+        output.push_back(bbox);
+    }
+    return output;
+}
+
 } // namespace
 
 proj_transform::proj_transform(projection const& source, projection const& dest)
@@ -398,7 +427,7 @@ bool proj_transform::backward(box2d<double>& box) const
 }
 
 // More robust, but expensive, bbox transform
-// in the face of proj4 out of bounds conditions.
+// in the face of libproj out of bounds conditions.
 // Can result in 20 -> 10 r/s performance hit.
 // Alternative is to provide proper clipping box
 // in the target srs by setting map 'maximum-extent'
@@ -413,18 +442,24 @@ bool proj_transform::backward(box2d<double>& env, std::size_t points) const
         return backward(env);
     }
 
-    auto coords = envelope_points(env, points); // this is always clockwise
+    box2d<double> result;
-
+    for (auto b : split_bounding_box(env))
+    {
+        auto coords = envelope_points(b, points); // this is always clockwise
         for (auto& p : coords)
         {
             double z = 0;
             if (!backward(p.x, p.y, z))
                 return false;
         }
-
+        if (!result.valid())
-    box2d<double> result;
             boost::geometry::envelope(coords, result);
-
+        else
+        {
+            box2d<double> bb;
+            boost::geometry::envelope(coords, bb);
+            result.expand_to_include(bb);
+        }
         if (is_source_longlat_ && !util::is_clockwise(coords))
         {
             // we've gone to a geographic CS, and our clockwise envelope has
@@ -435,6 +470,7 @@ bool proj_transform::backward(box2d<double>& env, std::size_t points) const
             result.expand_to_include(-180.0, miny);
             result.expand_to_include(180.0, miny);
         }
+    }
 
     env.re_center(result.center().x, result.center().y);
     env.height(result.height());
@@ -453,17 +489,31 @@ bool proj_transform::forward(box2d<double>& env, std::size_t points) const
     }
 
     auto coords = envelope_points(env, points); // this is always clockwise
+    double z = 0;
+    for (auto& p : coords)
+    {
+        if (!forward(p.x, p.y, z))
+            return false;
+    }
 
+    box2d<double> result;
+    for (auto b : split_bounding_box(env))
+    {
+        auto coords = envelope_points(b, points); // this is always clockwise
         for (auto& p : coords)
         {
             double z = 0;
             if (!forward(p.x, p.y, z))
                 return false;
         }
-
+        if (!result.valid())
-    box2d<double> result;
             boost::geometry::envelope(coords, result);
-
+        else
+        {
+            box2d<double> bb;
+            boost::geometry::envelope(coords, bb);
+            result.expand_to_include(bb);
+        }
         if (is_dest_longlat_ && !util::is_clockwise(coords))
         {
             // we've gone to a geographic CS, and our clockwise envelope has
@@ -474,11 +524,10 @@ bool proj_transform::forward(box2d<double>& env, std::size_t points) const
             result.expand_to_include(-180.0, miny);
             result.expand_to_include(180.0, miny);
         }
-
+    }
     env.re_center(result.center().x, result.center().y);
     env.height(result.height());
     env.width(result.width());
-
     return true;
 }
 

src/projection.cpp

@@ -94,7 +94,7 @@ bool projection::operator!=(const projection& other) const
     return !(*this == other);
 }
 
-void projection::init_proj() const
+void projection::init_proj()
 {
 #ifdef MAPNIK_USE_PROJ
     if (!proj_)
@@ -117,6 +117,11 @@ void projection::init_proj() const
         }
         PJ_TYPE type = proj_get_type(proj_);
         is_geographic_ = (type == PJ_TYPE_GEOGRAPHIC_2D_CRS || type == PJ_TYPE_GEOGRAPHIC_3D_CRS) ? true : false;
+        double west_lon, south_lat, east_lon, north_lat;
+        if (proj_get_area_of_use(proj_ctx_, proj_, &west_lon, &south_lat, &east_lon, &north_lat, nullptr))
+        {
+            area_of_use_ = box2d<double>{west_lon, south_lat, east_lon, north_lat};
+        }
     }
 #endif
 }
@@ -131,6 +136,11 @@ bool projection::is_geographic() const
     return is_geographic_;
 }
 
+std::optional<box2d<double>> projection::area_of_use() const
+{
+    return area_of_use_;
+}
+
 std::optional<well_known_srs_e> projection::well_known() const
 {
     return is_well_known_srs(params_);

test/unit/projection/proj_transform.cpp

@@ -16,6 +16,8 @@ TEST_CASE("projection transform")
         CHECK(proj_4326.is_geographic());
         CHECK(!proj_3857.is_geographic());
 
+        CHECK(*proj_4326.area_of_use() == mapnik::box2d<double>(-180, -90, 180, 90));
+        CHECK(*proj_3857.area_of_use() == mapnik::box2d<double>(-180, -85.06, 180, 85.06));
         mapnik::proj_transform prj_trans(proj_4326, proj_3857);
 
         double minx = -45.0;
@@ -48,6 +50,7 @@ TEST_CASE("projection transform")
     {
         mapnik::projection proj_4269("epsg:4269");
         mapnik::projection proj_3857("epsg:3857");
+
         mapnik::proj_transform prj_trans(proj_4269, proj_3857);
         mapnik::proj_transform prj_trans2(proj_3857, proj_4269);
 
@@ -200,7 +203,42 @@ TEST_CASE("projection transform")
             }
         }
     }
-#endif
+    SECTION("Test proj north/south poles bbox")
+    {
+        // WGS 84 / Arctic Polar Stereographic
+
+        mapnik::projection prj_geog("epsg:4326");
+        mapnik::projection prj_proj("epsg:3995");
+
+        mapnik::proj_transform prj_trans_fwd(prj_proj, prj_geog);
+        mapnik::proj_transform prj_trans_rev(prj_geog, prj_proj);
+
+        // bounds
+        const mapnik::box2d<double> bounds{-180.0, 60.0, 180.0, 90.0};
+        {
+            // projected bounds
+            mapnik::box2d<double> projected_bounds{-3299207.53, -3333134.03, 3299207.53, 3333134.03};
+            CHECKED_IF(prj_trans_fwd.forward(projected_bounds, PROJ_ENVELOPE_POINTS))
+            {
+                CHECK(projected_bounds.minx() == Approx(bounds.minx()));
+                CHECK(projected_bounds.miny() == Approx(48.65640)); // this is expected
+                CHECK(projected_bounds.maxx() == Approx(bounds.maxx()));
+                CHECK(projected_bounds.maxy() == Approx(bounds.maxy()));
+            }
+        }
+
+        {
+            // check the same logic works for .backward()
+            mapnik::box2d<double> projected_bounds{-3299207.53, -3333134.03, 3299207.53, 3333134.03};
+            CHECKED_IF(prj_trans_rev.backward(projected_bounds, PROJ_ENVELOPE_POINTS))
+            {
+                CHECK(projected_bounds.minx() == Approx(bounds.minx()));
+                CHECK(projected_bounds.miny() == Approx(48.65640)); // this is expected
+                CHECK(projected_bounds.maxx() == Approx(bounds.maxx()));
+                CHECK(projected_bounds.maxy() == Approx(bounds.maxy()));
+            }
+        }
+    }
 
     SECTION("proj_transform of coordinate arrays with stride > 1")
     {
@@ -237,7 +275,6 @@ TEST_CASE("projection transform")
             }
         }
 
-#ifdef MAPNIK_USE_PROJ
         SECTION("lonlat <-> New Zealand Transverse Mercator 2000")
         {
             mapnik::projection const proj_2193("epsg:2193");