sij/mapnik
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Handle bounds reprojections which cross the anti-meridian.

Browse source 
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.

Backport of #2657 from master to 2.3.x. Fixes #2648.
This commit is contained in:
Robert Coup 2015-01-28 13:20:33 +13:00
parent 7c4325fd76
commit a27b886e19
9 changed files with 229 additions and 10 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

58
src/proj_transform.cpp
View file

@ -216,6 +216,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();
@ -233,15 +235,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) {
@ -268,7 +287,7 @@ bool proj_transform::backward(box2d<double>& env, int points) const
return true;
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) {
@ -279,6 +298,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());
@ -293,7 +322,7 @@ bool proj_transform::forward(box2d<double>& env, int points) const
return true;
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) {
@ -305,6 +334,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());

36
tests/python_tests/projection_test.py
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_")))

8
tests/python_tests/utilities.py
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)

44
tests/visual_tests/grids/antimeridian-bbox-500-100-1.0-grid-reference.json Normal file
View file

@ -0,0 +1,44 @@
{
"keys": [
"",
"47",
"39",
"24",
"31",
"17",
"11",
"48",
"6",
"54",
"1",
"62"
],
"data": {},
"grid": [
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" !#$ ",
" % &' ( ",
" ) * ",
" + , ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" ",
" "
]
}

BIN
tests/visual_tests/images/antimeridian-bbox-500-100-1.0-agg-reference.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 202 B

BIN
tests/visual_tests/images/antimeridian-bbox-500-100-1.0-cairo-reference.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 448 B

BIN
tests/visual_tests/images/antimeridian-bbox-500-100-2.0-agg-reference.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 363 B

BIN
tests/visual_tests/images/antimeridian-bbox-500-100-2.0-cairo-reference.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 465 B

93
tests/visual_tests/styles/antimeridian-bbox.xml Normal file
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>