123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370 |
- #ifndef BOOST_GEOMETRY_PROJECTIONS_GEOS_HPP
- #define BOOST_GEOMETRY_PROJECTIONS_GEOS_HPP
- #include <boost/math/special_functions/hypot.hpp>
- #include <boost/geometry/srs/projections/impl/base_static.hpp>
- #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
- #include <boost/geometry/srs/projections/impl/projects.hpp>
- #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
- #include <boost/geometry/srs/projections/impl/pj_param.hpp>
- namespace boost { namespace geometry
- {
- namespace projections
- {
- #ifndef DOXYGEN_NO_DETAIL
- namespace detail { namespace geos
- {
- template <typename T>
- struct par_geos
- {
- T h;
- T radius_p;
- T radius_p2;
- T radius_p_inv2;
- T radius_g;
- T radius_g_1;
- T C;
- bool flip_axis;
- };
- template <typename T, typename Parameters>
- struct base_geos_ellipsoid
- {
- par_geos<T> m_proj_parm;
-
-
- inline void fwd(Parameters const& , T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
- {
- T r, Vx, Vy, Vz, tmp;
-
- lp_lat = atan (this->m_proj_parm.radius_p2 * tan (lp_lat));
-
- r = (this->m_proj_parm.radius_p) / boost::math::hypot(this->m_proj_parm.radius_p * cos (lp_lat), sin (lp_lat));
- Vx = r * cos (lp_lon) * cos (lp_lat);
- Vy = r * sin (lp_lon) * cos (lp_lat);
- Vz = r * sin (lp_lat);
-
- if (((this->m_proj_parm.radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * this->m_proj_parm.radius_p_inv2) < 0.) {
- BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
- }
-
- tmp = this->m_proj_parm.radius_g - Vx;
- if(this->m_proj_parm.flip_axis) {
- xy_x = this->m_proj_parm.radius_g_1 * atan (Vy / boost::math::hypot (Vz, tmp));
- xy_y = this->m_proj_parm.radius_g_1 * atan (Vz / tmp);
- } else {
- xy_x = this->m_proj_parm.radius_g_1 * atan (Vy / tmp);
- xy_y = this->m_proj_parm.radius_g_1 * atan (Vz / boost::math::hypot (Vy, tmp));
- }
- }
-
-
- inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
- {
- T Vx, Vy, Vz, a, b, det, k;
-
- Vx = -1.0;
- if(this->m_proj_parm.flip_axis) {
- Vz = tan (xy_y / this->m_proj_parm.radius_g_1);
- Vy = tan (xy_x / this->m_proj_parm.radius_g_1) * boost::math::hypot(1.0, Vz);
- } else {
- Vy = tan (xy_x / this->m_proj_parm.radius_g_1);
- Vz = tan (xy_y / this->m_proj_parm.radius_g_1) * boost::math::hypot(1.0, Vy);
- }
-
- a = Vz / this->m_proj_parm.radius_p;
- a = Vy * Vy + a * a + Vx * Vx;
- b = 2 * this->m_proj_parm.radius_g * Vx;
- if ((det = (b * b) - 4 * a * this->m_proj_parm.C) < 0.) {
- BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
- }
-
- k = (-b - sqrt(det)) / (2. * a);
- Vx = this->m_proj_parm.radius_g + k * Vx;
- Vy *= k;
- Vz *= k;
-
- lp_lon = atan2 (Vy, Vx);
- lp_lat = atan (Vz * cos (lp_lon) / Vx);
- lp_lat = atan (this->m_proj_parm.radius_p_inv2 * tan (lp_lat));
- }
- static inline std::string get_name()
- {
- return "geos_ellipsoid";
- }
- };
- template <typename T, typename Parameters>
- struct base_geos_spheroid
- {
- par_geos<T> m_proj_parm;
-
-
- inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
- {
- T Vx, Vy, Vz, tmp;
-
- tmp = cos(lp_lat);
- Vx = cos (lp_lon) * tmp;
- Vy = sin (lp_lon) * tmp;
- Vz = sin (lp_lat);
-
-
- if (((this->m_proj_parm.radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.)
- BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
-
- tmp = this->m_proj_parm.radius_g - Vx;
- if(this->m_proj_parm.flip_axis) {
- xy_x = this->m_proj_parm.radius_g_1 * atan(Vy / boost::math::hypot(Vz, tmp));
- xy_y = this->m_proj_parm.radius_g_1 * atan(Vz / tmp);
- } else {
- xy_x = this->m_proj_parm.radius_g_1 * atan(Vy / tmp);
- xy_y = this->m_proj_parm.radius_g_1 * atan(Vz / boost::math::hypot(Vy, tmp));
- }
- }
-
-
- inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
- {
- T Vx, Vy, Vz, a, b, det, k;
-
- Vx = -1.0;
- if(this->m_proj_parm.flip_axis) {
- Vz = tan (xy_y / (this->m_proj_parm.radius_g - 1.0));
- Vy = tan (xy_x / (this->m_proj_parm.radius_g - 1.0)) * sqrt (1.0 + Vz * Vz);
- } else {
- Vy = tan (xy_x / (this->m_proj_parm.radius_g - 1.0));
- Vz = tan (xy_y / (this->m_proj_parm.radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);
- }
-
- a = Vy * Vy + Vz * Vz + Vx * Vx;
- b = 2 * this->m_proj_parm.radius_g * Vx;
- if ((det = (b * b) - 4 * a * this->m_proj_parm.C) < 0.) {
- BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
- }
-
- k = (-b - sqrt(det)) / (2 * a);
- Vx = this->m_proj_parm.radius_g + k * Vx;
- Vy *= k;
- Vz *= k;
-
- lp_lon = atan2 (Vy, Vx);
- lp_lat = atan (Vz * cos (lp_lon) / Vx);
- }
- static inline std::string get_name()
- {
- return "geos_spheroid";
- }
- };
- inline bool geos_flip_axis(srs::detail::proj4_parameters const& params)
- {
- std::string sweep_axis = pj_get_param_s(params, "sweep");
- if (sweep_axis.empty())
- return false;
- else {
- if (sweep_axis[1] != '\0' || (sweep_axis[0] != 'x' && sweep_axis[0] != 'y'))
- BOOST_THROW_EXCEPTION( projection_exception(error_invalid_sweep_axis) );
- if (sweep_axis[0] == 'x')
- return true;
- else
- return false;
- }
- }
- template <typename T>
- inline bool geos_flip_axis(srs::dpar::parameters<T> const& params)
- {
- typename srs::dpar::parameters<T>::const_iterator
- it = pj_param_find(params, srs::dpar::sweep);
- if (it == params.end()) {
- return false;
- } else {
- srs::dpar::value_sweep s = static_cast<srs::dpar::value_sweep>(it->template get_value<int>());
- return s == srs::dpar::sweep_x;
- }
- }
-
- template <typename Params, typename Parameters, typename T>
- inline void setup_geos(Params const& params, Parameters& par, par_geos<T>& proj_parm)
- {
- std::string sweep_axis;
- if ((proj_parm.h = pj_get_param_f<T, srs::spar::h>(params, "h", srs::dpar::h)) <= 0.)
- BOOST_THROW_EXCEPTION( projection_exception(error_h_less_than_zero) );
- if (par.phi0 != 0.0)
- BOOST_THROW_EXCEPTION( projection_exception(error_unknown_prime_meridian) );
- proj_parm.flip_axis = geos_flip_axis(params);
- proj_parm.radius_g_1 = proj_parm.h / par.a;
- proj_parm.radius_g = 1. + proj_parm.radius_g_1;
- proj_parm.C = proj_parm.radius_g * proj_parm.radius_g - 1.0;
- if (par.es != 0.0) {
- proj_parm.radius_p = sqrt (par.one_es);
- proj_parm.radius_p2 = par.one_es;
- proj_parm.radius_p_inv2 = par.rone_es;
- } else {
- proj_parm.radius_p = proj_parm.radius_p2 = proj_parm.radius_p_inv2 = 1.0;
- }
- }
- }}
- #endif
-
- template <typename T, typename Parameters>
- struct geos_ellipsoid : public detail::geos::base_geos_ellipsoid<T, Parameters>
- {
- template <typename Params>
- inline geos_ellipsoid(Params const& params, Parameters const& par)
- {
- detail::geos::setup_geos(params, par, this->m_proj_parm);
- }
- };
-
- template <typename T, typename Parameters>
- struct geos_spheroid : public detail::geos::base_geos_spheroid<T, Parameters>
- {
- template <typename Params>
- inline geos_spheroid(Params const& params, Parameters const& par)
- {
- detail::geos::setup_geos(params, par, this->m_proj_parm);
- }
- };
- #ifndef DOXYGEN_NO_DETAIL
- namespace detail
- {
-
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_geos, geos_spheroid, geos_ellipsoid)
-
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(geos_entry, geos_spheroid, geos_ellipsoid)
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(geos_init)
- {
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(geos, geos_entry);
- }
- }
- #endif
- }
- }}
- #endif
|