From 78835c9d87f2cc99dc15ba6318ad52ec8305c52e Mon Sep 17 00:00:00 2001 From: Blake Thompson Date: Wed, 1 Apr 2015 16:18:05 -0500 Subject: [PATCH] Initial commit for new method to provide geometry reprojection --- include/mapnik/geometry_reprojection.hpp | 43 ++++ include/mapnik/proj_transform.hpp | 9 + include/mapnik/well_known_srs.hpp | 31 +++ src/build.py | 1 + src/geometry_reprojection.cpp | 257 +++++++++++++++++++++++ src/proj_transform.cpp | 70 ++++++ 6 files changed, 411 insertions(+) create mode 100644 include/mapnik/geometry_reprojection.hpp create mode 100644 src/geometry_reprojection.cpp diff --git a/include/mapnik/geometry_reprojection.hpp b/include/mapnik/geometry_reprojection.hpp new file mode 100644 index 000000000..6102b6537 --- /dev/null +++ b/include/mapnik/geometry_reprojection.hpp @@ -0,0 +1,43 @@ +/***************************************************************************** + * + * This file is part of Mapnik (c++ mapping toolkit) + * + * Copyright (C) 2014 Artem Pavlenko + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + * + *****************************************************************************/ + +#ifndef MAPNIK_GEOMETRY_REPROJECTION_HPP +#define MAPNIK_GEOMETRY_REPROJECTION_HPP + +// mapnik +#include +#include + +// std +#include + +namespace mapnik { + +MAPNIK_DECL geometry::geometry reproject(geometry::geometry const& geom, proj_transform const& proj_trans, unsigned int & n_err, bool backwards = false); + +MAPNIK_DECL geometry::geometry reproject(geometry::geometry const& geom, projection const& source, projection const& dest, unsigned int & n_err); + +MAPNIK_DECL geometry::geometry reproject(geometry::geometry const& geom, std::string const& source, std::string const& dest, unsigned int & n_err); + +} // end mapnik ns + +#endif // MAPNIK_GEOMETRY_REPROJECTION_HPP diff --git a/include/mapnik/proj_transform.hpp b/include/mapnik/proj_transform.hpp index aa32a1aa0..da26d242c 100644 --- a/include/mapnik/proj_transform.hpp +++ b/include/mapnik/proj_transform.hpp @@ -29,6 +29,10 @@ namespace mapnik { +namespace geometry { +struct point; +struct line_string; +} class projection; template class box2d; @@ -39,10 +43,15 @@ public: projection const& dest); bool equal() const; + bool is_known() const; bool forward (double& x, double& y , double& z) const; bool backward (double& x, double& y , double& z) const; bool forward (double *x, double *y , double *z, int point_count) const; bool backward (double *x, double *y , double *z, int point_count) const; + bool forward (geometry::point & p) const; + bool backward (geometry::point & p) const; + unsigned int forward (geometry::line_string & ls) const; + unsigned int backward (geometry::line_string & ls) const; bool forward (box2d & box) const; bool backward (box2d & box) const; bool forward (box2d & box, int points) const; diff --git a/include/mapnik/well_known_srs.hpp b/include/mapnik/well_known_srs.hpp index 2c0d40e1e..13251e88e 100644 --- a/include/mapnik/well_known_srs.hpp +++ b/include/mapnik/well_known_srs.hpp @@ -26,6 +26,7 @@ // mapnik #include // for M_PI on windows #include +#include // boost #include @@ -89,6 +90,36 @@ static inline bool merc2lonlat(double * x, double * y , int point_count) return true; } +static inline bool lonlat2merc(geometry::line_string & ls) +{ + for(auto p: ls) + { + if (p.x > 180) p.x = 180; + else if (p.x < -180) p.x = -180; + if (p.y > MAX_LATITUDE) p.y = MAX_LATITUDE; + else if (p.y < -MAX_LATITUDE) p.y = -MAX_LATITUDE; + p.x = p.x * MAXEXTENTby180; + p.y = std::log(std::tan((90 + p.y) * M_PIby360)) * R2D; + p.y = p.y * MAXEXTENTby180; + } + return true; +} + +static inline bool merc2lonlat(geometry::line_string & ls) +{ + for (auto p : ls) + { + if (p.x > MAXEXTENT) p.x = MAXEXTENT; + else if (p.x < -MAXEXTENT) p.x = -MAXEXTENT; + if (p.y > MAXEXTENT) p.y = MAXEXTENT; + else if (p.y < -MAXEXTENT) p.y = -MAXEXTENT; + p.x = (p.x / MAXEXTENT) * 180; + p.y = (p.y / MAXEXTENT) * 180; + p.y = R2D * (2 * std::atan(std::exp(p.y * D2R)) - M_PI_by2); + } + return true; +} + } #endif // MAPNIK_WELL_KNOWN_SRS_HPP diff --git a/src/build.py b/src/build.py index 0d9ebcaf9..e6179fb66 100644 --- a/src/build.py +++ b/src/build.py @@ -159,6 +159,7 @@ source = Split( datasource_cache.cpp datasource_cache_static.cpp debug.cpp + geometry_reprojection.cpp expression_node.cpp expression_string.cpp expression.cpp diff --git a/src/geometry_reprojection.cpp b/src/geometry_reprojection.cpp new file mode 100644 index 000000000..bdddc0046 --- /dev/null +++ b/src/geometry_reprojection.cpp @@ -0,0 +1,257 @@ +/***************************************************************************** + * + * This file is part of Mapnik (c++ mapping toolkit) + * + * Copyright (C) 2014 Artem Pavlenko + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + * + *****************************************************************************/ + +// mapnik +#include + +namespace mapnik { + +namespace detail { + +struct geom_reproj_visitor { + + geom_reproj_visitor(proj_transform const & proj_trans, unsigned int & n_err, bool reverse) + : proj_trans_(proj_trans), n_err_(n_err), reverse_(reverse) {} + + geometry::geometry operator() (geometry::geometry_empty const&) + { + return geometry::geometry(geometry::geometry_empty()); + } + + bool trans_point(geometry::point & new_p) + { + if (!reverse_) + { + if (!proj_trans_.forward(new_p)) + { + n_err_++; + return false; + } + } + else + { + if (!proj_trans_.backward(new_p)) + { + n_err_++; + return false; + } + } + return true; + } + + geometry::geometry operator() (geometry::point const & p) + { + geometry::point new_p(p); + if (!trans_point(new_p)) + { + return geometry::geometry(geometry::geometry_empty()); + } + return geometry::geometry(std::move(new_p)); + } + + bool trans_ls(geometry::line_string & new_ls) + { + unsigned int err = 0; + if (!reverse_) + { + err = proj_trans_.forward(new_ls); + } + else + { + err = proj_trans_.backward(new_ls); + } + if (err > 0) + { + n_err_ += err; + return false; + } + return true; + } + + geometry::geometry operator() (geometry::line_string const & ls) + { + geometry::line_string new_ls(ls); + if (!trans_ls(new_ls)) + { + return geometry::geometry(geometry::geometry_empty()); + } + return geometry::geometry(std::move(new_ls)); + } + + bool trans_poly(geometry::polygon & new_poly, geometry::polygon const & poly) + { + geometry::linear_ring new_ext(poly.exterior_ring); + int err = 0; + if (!reverse_) + { + err = proj_trans_.forward(new_ext); + } + else + { + err = proj_trans_.backward(new_ext); + } + // If the exterior ring doesn't transform don't bother with the holes. + if (err > 0) + { + n_err_ += err; + return false; + } + new_poly.set_exterior_ring(std::move(new_ext)); + new_poly.interior_rings.reserve(poly.interior_rings.size()); + + for (auto lr : poly.interior_rings) + { + geometry::linear_ring new_lr(lr); + if (!reverse_) + { + err = proj_trans_.forward(new_lr); + } + else + { + err = proj_trans_.backward(new_lr); + } + if (err > 0) + { + n_err_ += err; + // If there is an error in interior ring drop + // it from polygon. + continue; + } + new_poly.add_hole(std::move(new_lr)); + } + return true; + } + + geometry::geometry operator() (geometry::polygon const & poly) + { + geometry::polygon new_poly; + if (!trans_poly(new_poly, poly)) + { + return std::move(geometry::geometry_empty()); + } + return std::move(new_poly); + } + + geometry::geometry operator() (geometry::multi_point const & mp) + { + if (proj_trans_.is_known()) + { + geometry::multi_point new_mp(mp); + if (!trans_ls(new_mp)) + { + // should be impossible to reach this currently + /* LCOV_EXCL_START */ + return std::move(geometry::geometry_empty()); + /* LCOV_EXCL_END */ + } + return std::move(new_mp); + } + geometry::multi_point new_mp; + new_mp.reserve(mp.size()); + for (auto p : mp) + { + geometry::point new_p(p); + if (trans_point(new_p)) + { + new_mp.emplace_back(new_p); + } + } + if (new_mp.empty()) + { + return std::move(geometry::geometry_empty()); + } + return std::move(new_mp); + } + + geometry::geometry operator() (geometry::multi_line_string const & mls) + { + geometry::multi_line_string new_mls; + new_mls.reserve(mls.size()); + for (auto ls : mls) + { + geometry::line_string new_ls(ls); + if (trans_ls(new_ls)) + { + new_mls.emplace_back(new_ls); + } + } + if (new_mls.empty()) + { + return std::move(geometry::geometry_empty()); + } + return std::move(new_mls); + } + + geometry::geometry operator() (geometry::multi_polygon const & mpoly) + { + geometry::multi_polygon new_mpoly; + new_mpoly.reserve(mpoly.size()); + for (auto poly : mpoly) + { + geometry::polygon new_poly; + if (trans_poly(new_poly, poly)) + { + new_mpoly.emplace_back(new_poly); + } + } + if (new_mpoly.empty()) + { + return std::move(geometry::geometry_empty()); + } + return std::move(new_mpoly); + } + + geometry::geometry operator() (geometry::geometry_collection const & c) + { + geometry::geometry_collection new_c; + new_c.reserve(c.size()); + for (auto g : c) + { + geometry::geometry new_g(std::move(reproject(g, proj_trans_, n_err_, reverse_))); + if (!new_g.is()) + { + new_c.emplace_back(new_g); + } + } + if (new_c.empty()) + { + return std::move(geometry::geometry_empty()); + } + return std::move(new_c); + } + + private: + proj_transform const& proj_trans_; + unsigned int & n_err_; + bool reverse_; + +}; + +} // end detail ns + +geometry::geometry reproject(geometry::geometry const& geom, proj_transform const& proj_trans, unsigned int & n_err, bool backwards) +{ + detail::geom_reproj_visitor visit(proj_trans, n_err, backwards); + return util::apply_visitor(visit, geom); +} + +} // end mapnik ns diff --git a/src/proj_transform.cpp b/src/proj_transform.cpp index 37ca8d703..0735b4daf 100644 --- a/src/proj_transform.cpp +++ b/src/proj_transform.cpp @@ -87,11 +87,49 @@ bool proj_transform::equal() const return is_source_equal_dest_; } +bool proj_transform::is_known() const +{ + return merc_to_wgs84_ || wgs84_to_merc_; +} + bool proj_transform::forward (double & x, double & y , double & z) const { return forward(&x, &y, &z, 1); } +bool proj_transform::forward (geometry::point & p) const +{ + double z = 0; + return forward(&(p.x), &(p.y), &z, 1); +} + +unsigned int proj_transform::forward (geometry::line_string & ls) const +{ + if (is_source_equal_dest_) + return 0; + + if (wgs84_to_merc_) + { + lonlat2merc(ls); + return 0; + } + else if (merc_to_wgs84_) + { + merc2lonlat(ls); + return 0; + } + + unsigned int err_n = 0; + for (auto p : ls) + { + if (!forward(p)) + { + err_n++; + } + } + return err_n; +} + bool proj_transform::forward (double * x, double * y , double * z, int point_count) const { @@ -182,6 +220,38 @@ bool proj_transform::backward (double & x, double & y , double & z) const return backward(&x, &y, &z, 1); } +bool proj_transform::backward (geometry::point & p) const +{ + double z = 0; + return backward(&(p.x), &(p.y), &z, 1); +} + +unsigned int proj_transform::backward (geometry::line_string & ls) const +{ + if (is_source_equal_dest_) + return 0; + + if (wgs84_to_merc_) + { + merc2lonlat(ls); + return 0; + } + else if (merc_to_wgs84_) + { + lonlat2merc(ls); + return 0; + } + + unsigned int err_n = 0; + for (auto p : ls) + { + if (!backward(p)) + { + err_n++; + } + } + return err_n; +} bool proj_transform::forward (box2d & box) const {