Merge pull request #3940 from lightmare/fix-bbox-reprojection-v3.0.x

Backport Fix bbox reprojection to v3.0.x

.travis.yml

@@ -103,10 +103,7 @@ script:
  # (and might work) for the next build
  - DURATION=2400
  - scripts/travis-command-wrapper.py -s "date" -i 120 --deadline=$(( $(date +%s) + ${DURATION} )) make
- - RESULT=0
+ - make test
- - make test || RESULT=$?
- # we allow visual failures with g++ for now: https://github.com/mapnik/mapnik/issues/3567
- - if [[ ${RESULT} != 0 ]] && [[ ${CXX} =~ 'clang++' ]]; then false; fi;
  - enabled ${COVERAGE} coverage
  - enabled ${BENCH} make bench
  - ./scripts/check_glibcxx.sh

src/proj_transform.cpp

@@ -21,13 +21,17 @@
  *****************************************************************************/
 
 // mapnik
-#include <mapnik/global.hpp>
 #include <mapnik/box2d.hpp>
+#include <mapnik/geometry.hpp>
+#include <mapnik/geometry_adapters.hpp>
 #include <mapnik/projection.hpp>
 #include <mapnik/proj_transform.hpp>
 #include <mapnik/coord.hpp>
 #include <mapnik/util/is_clockwise.hpp>
 
+// boost
+#include <boost/geometry/algorithms/envelope.hpp>
+
 #ifdef MAPNIK_USE_PROJ4
 // proj4
 #include <proj_api.h>
@@ -39,6 +43,56 @@
 
 namespace mapnik {
 
+namespace { // (local)
+
+// Returns points in clockwise order. This allows us to do anti-meridian checks.
+template <typename T>
+auto envelope_points(box2d<T> const& env, std::size_t num_points)
+    -> geometry::multi_point<T>
+{
+    auto width = env.width();
+    auto height = env.height();
+
+    geometry::multi_point<T> coords;
+    coords.reserve(num_points);
+
+    // top side: left >>> right
+    // gets extra point if (num_points % 4 >= 1)
+    for (std::size_t i = 0, n = (num_points + 3) / 4; i < n; ++i)
+    {
+        auto x = env.minx() + (i * width) / n;
+        coords.emplace_back(x, env.maxy());
+    }
+
+    // right side: top >>> bottom
+    // gets extra point if (num_points % 4 >= 3)
+    for (std::size_t i = 0, n = (num_points + 1) / 4; i < n; ++i)
+    {
+        auto y = env.maxy() - (i * height) / n;
+        coords.emplace_back(env.maxx(), y);
+    }
+
+    // bottom side: right >>> left
+    // gets extra point if (num_points % 4 >= 2)
+    for (std::size_t i = 0, n = (num_points + 2) / 4; i < n; ++i)
+    {
+        auto x = env.maxx() - (i * width) / n;
+        coords.emplace_back(x, env.miny());
+    }
+
+    // left side: bottom >>> top
+    // never gets extra point
+    for (std::size_t i = 0, n = (num_points + 0) / 4; i < n; ++i)
+    {
+        auto y = env.miny() + (i * height) / n;
+        coords.emplace_back(env.minx(), y);
+    }
+
+    return coords;
+}
+
+} // namespace mapnik::(local)
+
 proj_transform::proj_transform(projection const& source,
                                projection const& dest)
     : source_(source),
@@ -334,49 +388,6 @@ bool proj_transform::backward (box2d<double> & box) const
     return true;
 }
 
-// Returns points in clockwise order. This allows us to do anti-meridian checks.
-void envelope_points(std::vector< coord<double,2> > & coords, box2d<double>& env, int points)
-{
-    double width = env.width();
-    double height = env.height();
-
-    int steps;
-
-    if (points <= 4) {
-        steps = 0;
-    } else {
-        steps = static_cast<int>(std::ceil((points - 4) / 4.0));
-    }
-
-    steps += 1;
-    double xstep = width / steps;
-    double ystep = height / steps;
-
-    coords.resize(points);
-    for (int i=0; i<steps; i++) {
-        // top: left>right
-        coords[i] = coord<double, 2>(env.minx() + i * xstep, env.maxy());
-        // right: top>bottom
-        coords[i + steps] = coord<double, 2>(env.maxx(), env.maxy() - i * ystep);
-        // bottom: right>left
-        coords[i + steps * 2] = coord<double, 2>(env.maxx() - i * xstep, env.miny());
-        // left: bottom>top
-        coords[i + steps * 3] = coord<double, 2>(env.minx(), env.miny() + i * ystep);
-    }
-}
-
-box2d<double> calculate_bbox(std::vector<coord<double,2> > & points) {
-    std::vector<coord<double,2> >::iterator it = points.begin();
-    std::vector<coord<double,2> >::iterator it_end = points.end();
-
-    box2d<double> env(*it, *(++it));
-    for (; it!=it_end; ++it) {
-        env.expand_to_include(*it);
-    }
-    return env;
-}
-
-
 // More robust, but expensive, bbox transform
 // in the face of proj4 out of bounds conditions.
 // Can result in 20 -> 10 r/s performance hit.
@@ -393,18 +404,18 @@ bool proj_transform::backward(box2d<double>& env, int points) const
         return backward(env);
     }
 
-    std::vector<coord<double,2> > coords;
+    auto coords = envelope_points(env, points);  // this is always clockwise
-    envelope_points(coords, env, points);  // this is always clockwise
 
-    double z;
+    for (auto & p : coords)
-    for (std::vector<coord<double,2> >::iterator it = coords.begin(); it!=coords.end(); ++it) {
+    {
-        z = 0;
+        double z = 0;
-        if (!backward(it->x, it->y, z)) {
+        if (!backward(p.x, p.y, z))
             return false;
-        }
     }
 
-    box2d<double> result = calculate_bbox(coords);
+    box2d<double> result;
+    boost::geometry::envelope(coords, result);
+
     if (is_source_longlat_ && !util::is_clockwise(coords))
     {
         // we've gone to a geographic CS, and our clockwise envelope has
@@ -432,18 +443,17 @@ bool proj_transform::forward(box2d<double>& env, int points) const
         return forward(env);
     }
 
-    std::vector<coord<double,2> > coords;
+    auto coords = envelope_points(env, points);  // this is always clockwise
-    envelope_points(coords, env, points);  // this is always clockwise
 
-    double z;
+    for (auto & p : coords)
-    for (std::vector<coord<double,2> >::iterator it = coords.begin(); it!=coords.end(); ++it) {
+    {
-        z = 0;
+        double z = 0;
-        if (!forward(it->x, it->y, z)) {
+        if (!forward(p.x, p.y, z))
             return false;
-        }
     }
 
-    box2d<double> result = calculate_bbox(coords);
+    box2d<double> result;
+    boost::geometry::envelope(coords, result);
 
     if (is_dest_longlat_ && !util::is_clockwise(coords))
     {

test/unit/projection/proj_transform.cpp

@@ -120,4 +120,80 @@ SECTION("test pj_transform failure behavior")
 
 #endif
 
+// Github Issue https://github.com/mapnik/mapnik/issues/2648
+SECTION("Test proj antimeridian bbox")
+{
+    mapnik::projection prj_geog("+init=epsg:4326");
+    mapnik::projection prj_proj("+init=epsg:2193");
+
+    mapnik::proj_transform prj_trans_fwd(prj_proj, prj_geog);
+    mapnik::proj_transform prj_trans_rev(prj_geog, prj_proj);
+
+    // reference values taken from proj4 command line tool:
+    // (non-corner points assume PROJ_ENVELOPE_POINTS == 20)
+    //
+    //  cs2cs -Ef %.10f +init=epsg:2193 +to +init=epsg:4326 <<END
+    //        2105800 3087000 # left-most
+    //        1495200 3087000 # bottom-most
+    //        2105800 7173000 # right-most
+    //        3327000 7173000 # top-most
+    //  END
+    //
+    // wrong = mapnik.Box2d(-177.3145325044, -62.3337481525,
+    //                       178.0277836332, -24.5845974912)
+    const mapnik::box2d<double> better(-180.0, -62.3337481525,
+                                        180.0, -24.5845974912);
+
+    {
+        mapnik::box2d<double> ext(274000, 3087000, 3327000, 7173000);
+        prj_trans_fwd.forward(ext, PROJ_ENVELOPE_POINTS);
+        CHECK(ext.minx() == Approx(better.minx()));
+        CHECK(ext.miny() == Approx(better.miny()));
+        CHECK(ext.maxx() == Approx(better.maxx()));
+        CHECK(ext.maxy() == Approx(better.maxy()));
+    }
+
+    {
+        // check the same logic works for .backward()
+        mapnik::box2d<double> ext(274000, 3087000, 3327000, 7173000);
+        prj_trans_rev.backward(ext, PROJ_ENVELOPE_POINTS);
+        CHECK(ext.minx() == Approx(better.minx()));
+        CHECK(ext.miny() == Approx(better.miny()));
+        CHECK(ext.maxx() == Approx(better.maxx()));
+        CHECK(ext.maxy() == Approx(better.maxy()));
+    }
+
+    // reference values taken from proj4 command line tool:
+    //
+    //  cs2cs -Ef %.10f +init=epsg:2193 +to +init=epsg:4326 <<END
+    //        274000 3087000 # left-most
+    //        276000 3087000 # bottom-most
+    //        276000 7173000 # right-most
+    //        274000 7173000 # top-most
+    //  END
+    //
+    const mapnik::box2d<double> normal(148.7667597489, -60.1222810241,
+                                       159.9548489296, -24.9771195155);
+
+    {
+        // checks for not being snapped (ie. not antimeridian)
+        mapnik::box2d<double> ext(274000, 3087000, 276000, 7173000);
+        prj_trans_fwd.forward(ext, PROJ_ENVELOPE_POINTS);
+        CHECK(ext.minx() == Approx(normal.minx()));
+        CHECK(ext.miny() == Approx(normal.miny()));
+        CHECK(ext.maxx() == Approx(normal.maxx()));
+        CHECK(ext.maxy() == Approx(normal.maxy()));
+    }
+
+    {
+        // check the same logic works for .backward()
+        mapnik::box2d<double> ext(274000, 3087000, 276000, 7173000);
+        prj_trans_rev.backward(ext, PROJ_ENVELOPE_POINTS);
+        CHECK(ext.minx() == Approx(normal.minx()));
+        CHECK(ext.miny() == Approx(normal.miny()));
+        CHECK(ext.maxx() == Approx(normal.maxx()));
+        CHECK(ext.maxy() == Approx(normal.maxy()));
+    }
+}
+
 }