[geometry] projections, added Natural Earth projection

[SVN r77231]
This commit is contained in:
Barend Gehrels
2012-03-05 14:41:02 +00:00
parent ef4d45ac9a
commit ee169cb7d7
4 changed files with 191 additions and 0 deletions

View File

@@ -70,6 +70,7 @@
#include <boost/geometry/extensions/gis/projections/proj/mill.hpp>
#include <boost/geometry/extensions/gis/projections/proj/mod_ster.hpp>
#include <boost/geometry/extensions/gis/projections/proj/moll.hpp>
#include <boost/geometry/extensions/gis/projections/proj/natearth.hpp>
#include <boost/geometry/extensions/gis/projections/proj/nell.hpp>
#include <boost/geometry/extensions/gis/projections/proj/nell_h.hpp>
#include <boost/geometry/extensions/gis/projections/proj/nocol.hpp>
@@ -190,6 +191,7 @@ public:
detail::mill_init(*this);
detail::mod_ster_init(*this);
detail::moll_init(*this);
detail::natearth_init(*this);
detail::nell_init(*this);
detail::nell_h_init(*this);
detail::nocol_init(*this);

View File

@@ -42,6 +42,7 @@
#include <boost/geometry/extensions/gis/projections/impl/base_dynamic.hpp>
#include <boost/geometry/extensions/gis/projections/impl/projects.hpp>
#include <boost/geometry/extensions/gis/projections/impl/factory_entry.hpp>
#include <boost/geometry/extensions/gis/projections/impl/aasincos.hpp>
namespace boost { namespace geometry { namespace projections
{

View File

@@ -0,0 +1,186 @@
#ifndef BOOST_GEOMETRY_PROJECTIONS_NATEARTH_HPP
#define BOOST_GEOMETRY_PROJECTIONS_NATEARTH_HPP
// Boost.Geometry - extensions-gis-projections (based on PROJ4)
// This file is automatically generated. DO NOT EDIT.
// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands.
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
// This file is converted from PROJ4, http://trac.osgeo.org/proj
// PROJ4 is originally written by Gerald Evenden (then of the USGS)
// PROJ4 is maintained by Frank Warmerdam
// PROJ4 is converted to Boost.Geometry by Barend Gehrels
// Original copyright notice:
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
#include <boost/math/special_functions/hypot.hpp>
#include <boost/geometry/extensions/gis/projections/impl/base_static.hpp>
#include <boost/geometry/extensions/gis/projections/impl/base_dynamic.hpp>
#include <boost/geometry/extensions/gis/projections/impl/projects.hpp>
#include <boost/geometry/extensions/gis/projections/impl/factory_entry.hpp>
namespace boost { namespace geometry { namespace projections
{
#ifndef DOXYGEN_NO_DETAIL
namespace detail { namespace natearth{
static const double A0 = 0.8707;
static const double A1 = -0.131979;
static const double A2 = -0.013791;
static const double A3 = 0.003971;
static const double A4 = -0.001529;
static const double B0 = 1.007226;
static const double B1 = 0.015085;
static const double B2 = -0.044475;
static const double B3 = 0.028874;
static const double B4 = -0.005916;
static const double C0 = B0;
static const double C1 = (3 * B1);
static const double C2 = (7 * B2);
static const double C3 = (9 * B3);
static const double C4 = (11 * B4);
static const double EPS = 1e-11;
static const double MAX_Y = (0.8707 * 0.52 * PI);
// template class, using CRTP to implement forward/inverse
template <typename Geographic, typename Cartesian, typename Parameters>
struct base_natearth_spheroid : public base_t_fi<base_natearth_spheroid<Geographic, Cartesian, Parameters>,
Geographic, Cartesian, Parameters>
{
typedef double geographic_type;
typedef double cartesian_type;
inline base_natearth_spheroid(const Parameters& par)
: base_t_fi<base_natearth_spheroid<Geographic, Cartesian, Parameters>,
Geographic, Cartesian, Parameters>(*this, par) {}
inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
{
double phi2, phi4;
phi2 = lp_lat * lp_lat;
phi4 = phi2 * phi2;
xy_x = lp_lon * (A0 + phi2 * (A1 + phi2 * (A2 + phi4 * phi2 * (A3 + phi2 * A4))));
xy_y = lp_lat * (B0 + phi2 * (B1 + phi4 * (B2 + B3 * phi2 + B4 * phi4)));
}
inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
{
double yc, tol, y2, y4, f, fder;
/* make sure y is inside valid range */
if (xy_y > MAX_Y) {
xy_y = MAX_Y;
} else if (xy_y < -MAX_Y) {
xy_y = -MAX_Y;
}
/* latitude */
yc = xy_y;
for (;;) { /* Newton-Raphson */
y2 = yc * yc;
y4 = y2 * y2;
f = (yc * (B0 + y2 * (B1 + y4 * (B2 + B3 * y2 + B4 * y4)))) - xy_y;
fder = C0 + y2 * (C1 + y4 * (C2 + C3 * y2 + C4 * y4));
yc -= tol = f / fder;
if (fabs(tol) < EPS) {
break;
}
}
lp_lat = yc;
/* longitude */
y2 = yc * yc;
lp_lon = xy_x / (A0 + y2 * (A1 + y2 * (A2 + y2 * y2 * y2 * (A3 + y2 * A4))));
}
};
// Natural Earth
template <typename Parameters>
void setup_natearth(Parameters& par)
{
par.es = 0;
// par.inv = s_inverse;
// par.fwd = s_forward;
}
}} // namespace detail::natearth
#endif // doxygen
/*!
\brief Natural Earth projection
\ingroup projections
\tparam Geographic latlong point type
\tparam Cartesian xy point type
\tparam Parameters parameter type
\par Projection characteristics
- Pseudocylindrical
- Spheroid
\par Example
\image html ex_natearth.gif
*/
template <typename Geographic, typename Cartesian, typename Parameters = parameters>
struct natearth_spheroid : public detail::natearth::base_natearth_spheroid<Geographic, Cartesian, Parameters>
{
inline natearth_spheroid(const Parameters& par) : detail::natearth::base_natearth_spheroid<Geographic, Cartesian, Parameters>(par)
{
detail::natearth::setup_natearth(this->m_par);
}
};
#ifndef DOXYGEN_NO_DETAIL
namespace detail
{
// Factory entry(s)
template <typename Geographic, typename Cartesian, typename Parameters>
class natearth_entry : public detail::factory_entry<Geographic, Cartesian, Parameters>
{
public :
virtual projection<Geographic, Cartesian>* create_new(const Parameters& par) const
{
return new base_v_fi<natearth_spheroid<Geographic, Cartesian, Parameters>, Geographic, Cartesian, Parameters>(par);
}
};
template <typename Geographic, typename Cartesian, typename Parameters>
inline void natearth_init(detail::base_factory<Geographic, Cartesian, Parameters>& factory)
{
factory.add_to_factory("natearth", new natearth_entry<Geographic, Cartesian, Parameters>);
}
} // namespace detail
#endif // doxygen
}}} // namespace boost::geometry::projections
#endif // BOOST_GEOMETRY_PROJECTIONS_NATEARTH_HPP

View File

@@ -188,6 +188,7 @@ void test_all()
test_forward<P>("murd1", 4.897000, 52.371000, 333340.993642, 5839071.944597, "+proj=murd1 +ellps=WGS84 +units=m +lat_1=20n +lat_2=60n");
test_forward<P>("murd2", 4.897000, 52.371000, 317758.821713, 6759296.097305, "+proj=murd2 +ellps=WGS84 +units=m +lat_1=20n +lat_2=60n");
test_forward<P>("murd3", 4.897000, 52.371000, 331696.409000, 5839224.186916, "+proj=murd3 +ellps=WGS84 +units=m +lat_1=20n +lat_2=60n");
test_forward<P>("natearth", 4.897000, 52.371000, 409886.629231, 5862282.218987, "+proj=natearth +ellps=WGS84 +units=m");
test_forward<P>("nell", 4.897000, 52.371000, 454576.246081, 5355027.851999, "+proj=nell +ellps=WGS84 +units=m");
test_forward<P>("nell_h", 4.897000, 52.371000, 438979.742911, 5386970.539995, "+proj=nell_h +ellps=WGS84 +units=m");
test_forward<P>("nicol", 4.897000, 52.371000, 360493.071000, 5836451.532406, "+proj=nicol +ellps=WGS84 +units=m");
@@ -304,6 +305,7 @@ void test_all()
test_inverse<P>("murd1", 333340.993642, 5839071.944597, 4.897000, 52.371000, "+proj=murd1 +ellps=WGS84 +units=m +lat_1=20n +lat_2=60n");
test_inverse<P>("murd2", 317758.821713, 6759296.097305, 4.897000, 52.371000, "+proj=murd2 +ellps=WGS84 +units=m +lat_1=20n +lat_2=60n");
test_inverse<P>("murd3", 331696.409000, 5839224.186916, 4.897000, 52.371000, "+proj=murd3 +ellps=WGS84 +units=m +lat_1=20n +lat_2=60n");
test_inverse<P>("natearth", 409886.629231, 5862282.218987, 4.897000, 52.371000, "+proj=natearth +ellps=WGS84 +units=m");
test_inverse<P>("nell", 454576.246081, 5355027.851999, 4.897000, 52.371000, "+proj=nell +ellps=WGS84 +units=m");
test_inverse<P>("nell_h", 438979.742911, 5386970.539995, 4.897000, 52.371000, "+proj=nell_h +ellps=WGS84 +units=m");
test_inverse<P>("nsper", 0.521191, 7.919806, 4.897000, 52.371000, "+proj=nsper +ellps=WGS84 +units=m +a=10 +h=40000000");