Merge pull request #2657 from rcoup/antimeridian-box-projection-2648

Handle bounds reprojections which cross the anti-meridian. Fixes #2648
This commit is contained in:
Robert Coup 2015-01-28 11:44:50 +13:00
commit 93f67f9bc9
10 changed files with 777 additions and 10 deletions

View file

@ -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);
}
}
/* 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);
}
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));
}
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());

View file

@ -5,7 +5,7 @@ from nose.tools import *
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_")))

View file

@ -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)

View file

@ -0,0 +1,306 @@
{
"keys": [
"",
"42",
"41",
"40",
"24",
"44",
"26",
"38",
"23",
"46",
"36",
"22",
"28",
"21",
"47",
"34",
"29",
"33",
"30",
"19",
"15",
"14",
"12",
"18",
"16",
"11",
"31",
"49",
"48",
"10",
"50",
"51",
"9",
"52",
"8",
"54",
"6",
"56",
"5",
"58",
"4",
"60",
"3",
"62",
"1",
"63",
"66",
"65",
"64"
],
"data": {},
"grid": [
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" !# $ % ",
" & ' ",
" ( ) ",
" * + , ",
" - . ",
" / 0 1 ",
" 2 3 4 56 7 ",
" 8 9 : ",
" ; <= ",
" > ? ",
" @ ",
" A B ",
" C D ",
" ",
" E F ",
" ",
" G H ",
" I J ",
" K L ",
" ",
" M N ",
" OP ",
" Q ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" "
]
}

View file

@ -0,0 +1,286 @@
{
"keys": [
"",
"42",
"40",
"24",
"37",
"46",
"22",
"28",
"35",
"20",
"30",
"15",
"12",
"32",
"17",
"48",
"10",
"50",
"9",
"52",
"54",
"7",
"57",
"5",
"60",
"3",
"1",
"63",
"66"
],
"data": {},
"grid": [
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ! # $ ",
" ",
" % ",
" & ' ",
" ( ",
" ) * ",
" + , - ",
" . / ",
" 0 ",
" 1 2 ",
" ",
" 3 4 ",
" 5 ",
" 6 ",
" ",
" 7 ",
" 8 ",
" 9 ",
" : ",
" ",
" ; < ",
" = ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" "
]
}

Binary file not shown.

After

Width:  |  Height:  |  Size: 938 B

Binary file not shown.

After

Width:  |  Height:  |  Size: 2.2 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.9 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 2.4 KiB

View file

@ -0,0 +1,93 @@
<Map srs="+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs " maximum-extent="274000,3087000,3327000,7173000">
<Parameters>
<Parameter name="sizes">1000,1000</Parameter>
<Parameter name="bbox">200000.0,1087000.0,7007000.0,7173000.0</Parameter>
</Parameters>
<Style name="mystyle" filter-mode="first">
<Rule>
<PointSymbolizer/>
</Rule>
</Style>
<Layer name="src_layer" srs="+init=epsg:4326">
<StyleName>mystyle</StyleName>
<Datasource>
<Parameter name="extent">-180,-47.2754736667,180,-34.0036943</Parameter>
<Parameter name="type">csv</Parameter>
<Parameter name="inline">
WKT,id,
"POINT (-176.662157383264883 -43.096109795145765)",
"POINT (-176.854671182715663 -42.903595931039007)",
"POINT (-177.252533072573982 -42.646910782400411)",
"POINT (-177.611892228793693 -42.415894147047872)",
"POINT (-178.099593973519461 -42.159208996768882)",
"POINT (-178.56162722744844 -41.735678488584703)",
"POINT (-178.72847257314487 -41.555998877995457)",
"POINT (-178.985157723427307 -41.222308171291367)",
"POINT (-179.152003072780786 -40.965623011559309)",
"POINT (-179.511362290676971 -40.490755465927869)",
"POINT (-179.80655022231997 -40.054390693250525)",
"POINT (179.923930360184784 -39.784871275110177)",
"POINT (179.75708500598833 -39.643694437033751)",
"POINT (179.346388748199729 -39.772037019858125)",
"POINT (179.025532296262895 -39.861876827647727)",
"POINT (178.563499004992451 -40.067224958370574)",
"POINT (178.281145326853419 -40.208401797959127)",
"POINT (178.02446016489057 -40.054390701181269)",
"POINT (177.677935196187519 -39.900379604399049)",
"POINT (177.485421324683159 -39.669362958813636)",
"POINT (177.331410227479353 -39.425512055075657)",
"POINT (177.215901904577578 -39.181661151299444)",
"POINT (177.023388033074042 -38.86080469896423)",
"POINT (176.76670287107433 -38.604119537111721)",
"POINT (176.535686225282944 -38.565616762876971)",
"POINT (176.317503837598906 -38.783799150566402)",
"POINT (176.176326998513417 -38.989147280139598)",
"POINT (175.932476094645182 -39.425512055467188)",
"POINT (175.701459448881849 -39.759202766009246)",
"POINT (175.52177983551357 -40.080059218446856)",
"POINT (175.149586350685809 -40.37524715469236)",
"POINT (174.880066930641647 -40.298241606113201)",
"POINT (174.700387317279876 -40.131396250851466)",
"POINT (174.469370671529759 -39.861876830812001)",
"POINT (174.418033639140958 -39.707865733646173)",
"POINT (174.212685509585867 -39.374175023120273)",
"POINT (174.148514219099951 -39.143158377371385)",
"POINT (174.05867441241972 -38.937810247816806)",
"POINT (173.930331831447972 -38.719627860165083)",
"POINT (173.699315185698993 -38.527113988707669)",
"POINT (173.327121700881207 -38.552782504902005)",
"POINT (172.993430990354909 -38.655456569679323)",
"POINT (172.839419893188989 -38.732462118262291)",
"POINT (172.659740279828696 -38.899307473525418)",
"POINT (172.45439215027406 -39.027650054497059)",
"POINT (172.197706988330737 -39.284335216440311)",
"POINT (172.018027374970359 -39.669362959355176)",
"POINT (-170.809757996124091 -39.258667337666161)",
"POINT (-171.092108440755396 -39.489683796643384)",
"POINT (-171.374459297976472 -39.72070029102813)",
"POINT (-171.669644665725599 -39.95171681078768)",
"POINT (-171.964830399850541 -40.195567616117053)",
"POINT (-172.208679740446428 -40.426584197598665)",
"POINT (-172.542368540167502 -40.631932262228219)",
"POINT (-172.696378898013165 -40.785943335635686)",
"POINT (-173.055736583802087 -41.068296973781855)",
"POINT (-173.273918197572641 -41.299313602634044)",
"POINT (-173.658944763471993 -41.645838554189496)",
"POINT (-173.504934101572132 -41.491827462434976)",
"POINT (-173.979800429072611 -41.992363516983445)",
"POINT (-174.108142742586438 -42.172043129406539)",
"POINT (-174.313490466537729 -42.377391260144812)",
"POINT (-174.608677910682047 -42.736750491079746)",
"POINT (-176.058947562102475 -43.429800480947527)",
"POINT (-174.955202379135329 -42.967767146816229)",
"POINT (-175.481406377455045 -43.160281038130286)",
"POINT (-175.802262571951218 -43.352794922023079)",
"POINT (-176.328466827030041 -43.339960685516779)",
"POINT (-176.508146353344472 -43.250120886092468)",
</Parameter>
</Datasource>
</Layer>
</Map>