Initial commit for new method to provide geometry reprojection

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 <mapnik/proj_transform.hpp>
+#include <mapnik/geometry.hpp>
+
+// std
+#include <string>
+
+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

include/mapnik/proj_transform.hpp

@@ -29,6 +29,10 @@
 
 namespace mapnik {
 
+namespace geometry {
+struct point;
+struct line_string;
+}
 class projection;
 template <typename T> 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<double> & box) const;
     bool backward (box2d<double> & box) const;
     bool forward (box2d<double> & box, int points) const;

include/mapnik/well_known_srs.hpp

@@ -26,6 +26,7 @@
 // mapnik
 #include <mapnik/global.hpp> // for M_PI on windows
 #include <mapnik/enumeration.hpp>
+#include <mapnik/geometry.hpp>
 
 // boost
 #include <boost/optional.hpp>
@@ -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

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

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 <mapnik/geometry_reprojection.hpp>
+
+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<geometry::geometry_empty>())
+            {
+                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

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<double> & box) const
 {