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

Merge pull request #3935 from lightmare/fix-bbox-reprojection

Browse source 
Fix bbox reprojection
This commit is contained in:
talaj 2018-07-10 21:45:47 +02:00 committed by GitHub
commit ae3a9404b7
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
4 changed files with 148 additions and 64 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
.travis.yml
View file

@ -106,10 +106,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 || RESULT=$?
# we allow visual failures with g++ for now: https://github.com/mapnik/mapnik/issues/3567
- if [[ ${RESULT} != 0 ]] && [[ ${CXX} =~ 'clang++' ]]; then false; fi;
- make test
- enabled ${COVERAGE} coverage
- enabled ${BENCH} make bench
- ./scripts/check_glibcxx.sh

1
CHANGELOG.md
View file

@ -16,6 +16,7 @@ For a complete change history, see the git log.
#### Core
- Rewrote parser grammars using Boost.Spirit X3 (replacing Boost.Spirit Qi)
- Fixed bbox reprojection ([#3935](https://github.com/mapnik/mapnik/issues/3935))
- Fixed segfault when ShieldSymbolizer has invalid placements ([5464ae9](https://github.com/mapnik/mapnik/commit/5464ae9cdfd32983883463bcc2396dc0e51d885f), [#3604](https://github.com/mapnik/mapnik/issues/3604#issuecomment-281542148))
- Slightly improved `sql_utils::table_from_sql` ([2587bb3](https://github.com/mapnik/mapnik/commit/2587bb3a1d8db397acfa8dcc2d332da3a8a9399f))
- Added wrappers for proper quoting in SQL query construction: `sql_utils::identifier`, `sql_utils::literal` ([7b21713](https://github.com/mapnik/mapnik/commit/7b217133e2749b82c2638551045c4edbece15086))

130
src/proj_transform.cpp
View file

@ -21,13 +21,17 @@
*****************************************************************************/
// mapnik
#include <mapnik/global.hpp>
#include <mapnik/geometry/boost_adapters.hpp>
#include <mapnik/geometry/box2d.hpp>
#include <mapnik/geometry/multi_point.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;
envelope_points(coords, env, points); // this is always clockwise
auto coords = envelope_points(env, points); // this is always clockwise
double z;
for (std::vector<coord<double,2> >::iterator it = coords.begin(); it!=coords.end(); ++it) {
z = 0;
if (!backward(it->x, it->y, z)) {
for (auto & p : coords)
{
double z = 0;
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;
envelope_points(coords, env, points); // this is always clockwise
auto coords = envelope_points(env, points); // this is always clockwise
double z;
for (std::vector<coord<double,2> >::iterator it = coords.begin(); it!=coords.end(); ++it) {
z = 0;
if (!forward(it->x, it->y, z)) {
for (auto & p : coords)
{
double z = 0;
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))
{

76
test/unit/projection/proj_transform.cpp
View file

@ -119,4 +119,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()));
}
}
}