Handle bounds reprojections which cross the anti-meridian. mapnik#2648

When doing an envelope-points reprojection to a geographic CS, check
the points stay in clockwise order. Otherwise expand the resulting bounds
to include the world.

Includes visual test. Cairo tests aren't finding differences, but the renderer
is doing the right thing.

src/proj_transform.cpp

@@ -219,6 +219,8 @@ 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();
@@ -236,15 +238,32 @@ void envelope_points(std::vector< coord<double,2> > & coords, box2d<double>& env
     double xstep = width / steps;
     double ystep = height / steps;
 
-    for (int i=0; i<=steps; i++) {
-        coords.push_back(coord<double,2>(env.minx() + i * xstep, env.miny()));
-        coords.push_back(coord<double,2>(env.minx() + i * xstep, env.maxy()));
-
+    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);
     }
-    for (int i=1; i<steps; i++) {
-        coords.push_back(coord<double,2>(env.minx(), env.miny() + i * ystep));
-        coords.push_back(coord<double,2>(env.maxx(), env.miny() + i * ystep));
+}
+
+/* determine if an ordered sequence of coordinates is in clockwise order */
+bool is_clockwise(const std::vector< coord<double,2> > & coords)
+{
+    int n = coords.size();
+    coord<double,2> c1, c2;
+    double a = 0.0;
+
+    for (int i=0; i<n; i++) {
+        c1 = coords[i];
+        c2 = coords[(i + 1) % n];
+        a += (c1.x * c2.y - c2.x * c1.y);
     }
+    return a <= 0.0;
 }
 
 box2d<double> calculate_bbox(std::vector<coord<double,2> > & points) {
@@ -276,7 +295,7 @@ bool proj_transform::backward(box2d<double>& env, int points) const
     }
 
     std::vector<coord<double,2> > coords;
-    envelope_points(coords, env, points);
+    envelope_points(coords, env, points);  // this is always clockwise
 
     double z;
     for (std::vector<coord<double,2> >::iterator it = coords.begin(); it!=coords.end(); ++it) {
@@ -287,6 +306,16 @@ bool proj_transform::backward(box2d<double>& env, int points) const
     }
 
     box2d<double> result = calculate_bbox(coords);
+    if (is_source_longlat_ && !is_clockwise(coords)) {
+        /* we've gone to a geographic CS, and our clockwise envelope has
+         * changed into an anticlockwise one. This means we've crossed the antimeridian, and
+         * need to expand the X direction to +/-180 to include all the data. Once we can deal
+         * with multiple bboxes in queries we can improve.
+         */
+         double miny = result.miny();
+         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());
@@ -306,7 +335,7 @@ bool proj_transform::forward(box2d<double>& env, int points) const
     }
 
     std::vector<coord<double,2> > coords;
-    envelope_points(coords, env, points);
+    envelope_points(coords, env, points);  // this is always clockwise
 
     double z;
     for (std::vector<coord<double,2> >::iterator it = coords.begin(); it!=coords.end(); ++it) {
@@ -318,6 +347,17 @@ bool proj_transform::forward(box2d<double>& env, int points) const
 
     box2d<double> result = calculate_bbox(coords);
 
+    if (is_dest_longlat_ && !is_clockwise(coords)) {
+        /* we've gone to a geographic CS, and our clockwise envelope has
+         * changed into an anticlockwise one. This means we've crossed the antimeridian, and
+         * need to expand the X direction to +/-180 to include all the data. Once we can deal
+         * with multiple bboxes in queries we can improve.
+         */
+         double miny = result.miny();
+         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());

tests/python_tests/projection_test.py

@@ -5,7 +5,7 @@
 import mapnik
 import random
 import math
-from utilities import run_all
+from utilities import run_all, assert_box2d_almost_equal
 
 # Tests that exercise map projections.
 
@@ -114,5 +114,39 @@ def test_proj_transform_between_init_and_literal():
             eq_(math.fabs(coord.y - lon_lat_coord4.y) < 1,True)
 
 
+# Github Issue #2648
+def test_proj_antimeridian_bbox():
+    # this is logic from feature_style_processor::prepare_layer()
+    PROJ_ENVELOPE_POINTS = 20  # include/mapnik/config.hpp
+
+    prjGeog = mapnik.Projection('+init=epsg:4326')
+    prjProj = mapnik.Projection('+init=epsg:2193')
+    prj_trans_fwd = mapnik.ProjTransform(prjProj, prjGeog)
+    prj_trans_rev = mapnik.ProjTransform(prjGeog, prjProj)
+
+    # bad = mapnik.Box2d(-177.31453250437079, -62.33374815225163, 178.02778363316355, -24.584597490955804)
+    better = mapnik.Box2d(-180.0, -62.33374815225163, 180.0, -24.584597490955804)
+
+    buffered_query_ext = mapnik.Box2d(274000, 3087000, 3327000, 7173000)
+    fwd_ext = prj_trans_fwd.forward(buffered_query_ext, PROJ_ENVELOPE_POINTS)
+    assert_box2d_almost_equal(fwd_ext, better)
+
+    # check the same logic works for .backward()
+    ext = mapnik.Box2d(274000, 3087000, 3327000, 7173000)
+    rev_ext = prj_trans_rev.backward(ext, PROJ_ENVELOPE_POINTS)
+    assert_box2d_almost_equal(rev_ext, better)
+
+    # checks for not being snapped (ie. not antimeridian)
+    normal = mapnik.Box2d(148.766759749,-60.1222810238,159.95484893,-24.9771195151)
+    buffered_query_ext = mapnik.Box2d(274000, 3087000, 276000, 7173000)
+    fwd_ext = prj_trans_fwd.forward(buffered_query_ext, PROJ_ENVELOPE_POINTS)
+    assert_box2d_almost_equal(fwd_ext, normal)
+
+    # check the same logic works for .backward()
+    ext = mapnik.Box2d(274000, 3087000, 276000, 7173000)
+    rev_ext = prj_trans_rev.backward(ext, PROJ_ENVELOPE_POINTS)
+    assert_box2d_almost_equal(rev_ext, normal)
+
+
 if __name__ == "__main__":
     exit(run_all(eval(x) for x in dir() if x.startswith("test_")))

tests/python_tests/utilities.py

@@ -2,6 +2,7 @@
 # -*- coding: utf-8 -*-
 
 from nose.plugins.errorclass import ErrorClass, ErrorClassPlugin
+from nose.tools import assert_almost_equal
 
 import os, sys, inspect, traceback
 import mapnik
@@ -86,3 +87,10 @@ def side_by_side_image(left_im, right_im):
     im.blend(0, 0, left_im, 1.0)
     im.blend(left_im.width() + 1, 0, right_im, 1.0)
     return im
+
+def assert_box2d_almost_equal(a, b, msg=None):
+    msg = msg or ("%r != %r" % (a, b))
+    assert_almost_equal(a.minx, b.minx, msg=msg)
+    assert_almost_equal(a.maxx, b.maxx, msg=msg)
+    assert_almost_equal(a.miny, b.miny, msg=msg)
+    assert_almost_equal(a.maxy, b.maxy, msg=msg)