diff options
Diffstat (limited to 'boost/geometry/srs/projections/impl')
29 files changed, 5100 insertions, 0 deletions
diff --git a/boost/geometry/srs/projections/impl/aasincos.hpp b/boost/geometry/srs/projections/impl/aasincos.hpp new file mode 100644 index 0000000000..4678029bc0 --- /dev/null +++ b/boost/geometry/srs/projections/impl/aasincos.hpp @@ -0,0 +1,114 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_AASINCOS_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_AASINCOS_HPP + + +#include <cmath> + +#include <boost/geometry/util/math.hpp> + + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + +namespace aasincos +{ + template <typename T> + inline T ONE_TOL() { return 1.00000000000001; } + //template <typename T> + //inline T TOL() { return 0.000000001; } + template <typename T> + inline T ATOL() { return 1e-50; } +} + +template <typename T> +inline T aasin(T const& v) +{ + T av = 0; + + if ((av = geometry::math::abs(v)) >= 1.0) + { + if (av > aasincos::ONE_TOL<T>()) + { + BOOST_THROW_EXCEPTION( projection_exception(-19) ); + } + return (v < 0.0 ? -geometry::math::half_pi<T>() : geometry::math::half_pi<T>()); + } + + return asin(v); +} + +template <typename T> +inline T aacos(T const& v) +{ + T av = 0; + + if ((av = geometry::math::abs(v)) >= 1.0) + { + if (av > aasincos::ONE_TOL<T>()) + { + BOOST_THROW_EXCEPTION( projection_exception(-19) ); + } + return (v < 0.0 ? geometry::math::pi<T>() : 0.0); + } + + return acos(v); +} + +template <typename T> +inline T asqrt(T const& v) +{ + return ((v <= 0) ? 0 : sqrt(v)); +} + +template <typename T> +inline T aatan2(T const& n, T const& d) +{ + return ((geometry::math::abs(n) < aasincos::ATOL<T>() + && geometry::math::abs(d) < aasincos::ATOL<T>()) ? 0.0 : atan2(n, d)); +} + + +} // namespace detail + + +}}} // namespace boost::geometry::projections + + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_AASINCOS_HPP diff --git a/boost/geometry/srs/projections/impl/adjlon.hpp b/boost/geometry/srs/projections/impl/adjlon.hpp new file mode 100644 index 0000000000..7f0c9ac4ca --- /dev/null +++ b/boost/geometry/srs/projections/impl/adjlon.hpp @@ -0,0 +1,70 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_ADJLON_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_ADJLON_HPP + +#include <boost/math/constants/constants.hpp> +#include <boost/geometry/util/math.hpp> + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + +/* reduce argument to range +/- PI */ +template <typename T> +inline T adjlon (T lon) +{ + if (geometry::math::abs(lon) <= boost::math::constants::pi<T>()) + { + return lon; + } + + /* adjust to 0..2pi rad */ + lon += boost::math::constants::pi<T>(); + /* remove integral # of 'revolutions'*/ + lon -= boost::math::constants::two_pi<T>() * + std::floor(lon / boost::math::constants::two_pi<T>()); + /* adjust back to -pi..pi rad */ + lon -= boost::math::constants::pi<T>(); + + return lon; +} + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_ADJLON_HPP diff --git a/boost/geometry/srs/projections/impl/base_dynamic.hpp b/boost/geometry/srs/projections/impl/base_dynamic.hpp new file mode 100644 index 0000000000..3979c34300 --- /dev/null +++ b/boost/geometry/srs/projections/impl/base_dynamic.hpp @@ -0,0 +1,151 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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) + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_BASE_DYNAMIC_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_BASE_DYNAMIC_HPP + +#include <string> + +#include <boost/geometry/srs/projections/impl/projects.hpp> + +namespace boost { namespace geometry { namespace projections +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +/*! + \brief projection virtual base class + \details class containing virtual methods + \ingroup projection + \tparam CT calculation type + \tparam P parameters type +*/ +template <typename CT, typename P> +class base_v +{ +public : + /// Forward projection, from Latitude-Longitude to Cartesian + template <typename LL, typename XY> + inline bool forward(LL const& lp, XY& xy) const + { + try + { + pj_fwd(*this, this->params(), lp, xy); + return true; + } + catch (...) + { + return false; + } + } + + /// Inverse projection, from Cartesian to Latitude-Longitude + template <typename LL, typename XY> + inline bool inverse(XY const& xy, LL& lp) const + { + try + { + pj_inv(*this, this->params(), xy, lp); + return true; + } + catch (projection_not_invertible_exception &) + { + BOOST_RETHROW + } + catch (...) + { + return false; + } + } + + /// Forward projection using lon / lat and x / y separately + virtual void fwd(CT& lp_lon, CT& lp_lat, CT& xy_x, CT& xy_y) const = 0; + + /// Inverse projection using x / y and lon / lat + virtual void inv(CT& xy_x, CT& xy_y, CT& lp_lon, CT& lp_lat) const = 0; + + /// Returns name of projection + virtual std::string name() const = 0; + + /// Returns parameters of projection + virtual P const& params() const = 0; + + /// Returns mutable parameters of projection + virtual P& mutable_params() = 0; + + virtual ~base_v() {} +}; + +// Base-virtual-forward +template <typename Prj, typename CT, typename P> +class base_v_f : public base_v<CT, P> +{ +public: + base_v_f(P const& params) + : m_proj(params) + {} + + template <typename ProjP> + base_v_f(P const& params, ProjP const& proj_params) + : m_proj(params, proj_params) + {} + + virtual void fwd(CT& lp_lon, CT& lp_lat, CT& xy_x, CT& xy_y) const + { + m_proj.fwd(lp_lon, lp_lat, xy_x, xy_y); + } + + virtual void inv(CT& , CT& , CT& , CT& ) const + { + BOOST_THROW_EXCEPTION(projection_not_invertible_exception(params().name)); + } + + virtual std::string name() const { return m_proj.name(); } + + virtual P const& params() const { return m_proj.params(); } + + virtual P& mutable_params() { return m_proj.mutable_params(); } + +protected: + Prj m_proj; +}; + +// Base-virtual-forward/inverse +template <typename Prj, typename CT, typename P> +class base_v_fi : public base_v_f<Prj, CT, P> +{ + typedef base_v_f<Prj, CT, P> base_t; + +public: + base_v_fi(P const& params) + : base_t(params) + {} + + template <typename ProjP> + base_v_fi(P const& params, ProjP const& proj_params) + : base_t(params, proj_params) + {} + + virtual void inv(CT& xy_x, CT& xy_y, CT& lp_lon, CT& lp_lat) const + { + this->m_proj.inv(xy_x, xy_y, lp_lon, lp_lat); + } +}; + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_BASE_DYNAMIC_HPP diff --git a/boost/geometry/srs/projections/impl/base_static.hpp b/boost/geometry/srs/projections/impl/base_static.hpp new file mode 100644 index 0000000000..aa205bcee5 --- /dev/null +++ b/boost/geometry/srs/projections/impl/base_static.hpp @@ -0,0 +1,139 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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) + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_BASE_STATIC_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_BASE_STATIC_HPP + +#if defined(_MSC_VER) +// For CRTP, *this is acceptable in constructor -> turn warning off +#pragma warning( disable : 4355 ) +#endif // defined(_MSC_VER) + + +#include <string> + +#include <boost/geometry/core/tags.hpp> + +#include <boost/geometry/srs/projections/impl/pj_fwd.hpp> +#include <boost/geometry/srs/projections/impl/pj_inv.hpp> + +#include <boost/mpl/assert.hpp> + + +namespace boost { namespace geometry { namespace projections +{ + + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +template <typename Prj, typename CSTag, typename BGP, typename CT, typename P> +struct static_projection_type +{ + BOOST_MPL_ASSERT_MSG((false), + NOT_IMPLEMENTED_FOR_THIS_PROJECTION_OR_CSTAG, + (Prj, CSTag)); +}; + +#define BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(PROJ, P_SPHERE, P_SPHEROID) \ +template <typename BGP, typename CT, typename P> \ +struct static_projection_type<PROJ, srs_sphere_tag, BGP, CT, P> \ +{ \ + typedef P_SPHERE<CT, P> type; \ +}; \ +template <typename BGP, typename CT, typename P> \ +struct static_projection_type<PROJ, srs_spheroid_tag, BGP, CT, P> \ +{ \ + typedef P_SPHEROID<CT, P> type; \ +}; \ + +// Base-template-forward +template <typename Prj, typename CT, typename P> +struct base_t_f +{ +public: + + inline base_t_f(Prj const& prj, P const& params) + : m_par(params), m_prj(prj) + {} + + inline P const& params() const { return m_par; } + + inline P& mutable_params() { return m_par; } + + template <typename LL, typename XY> + inline bool forward(LL const& lp, XY& xy) const + { + try + { + pj_fwd(m_prj, m_par, lp, xy); + return true; + } + catch(...) + { + return false; + } + } + + template <typename XY, typename LL> + inline bool inverse(XY const& , LL& ) const + { + BOOST_MPL_ASSERT_MSG((false), + PROJECTION_IS_NOT_INVERTABLE, + (Prj)); + return false; + } + + inline std::string name() const + { + return this->m_par.name; + } + +protected: + + P m_par; + const Prj& m_prj; +}; + +// Base-template-forward/inverse +template <typename Prj, typename CT, typename P> +struct base_t_fi : public base_t_f<Prj, CT, P> +{ +public : + inline base_t_fi(Prj const& prj, P const& params) + : base_t_f<Prj, CT, P>(prj, params) + {} + + template <typename XY, typename LL> + inline bool inverse(XY const& xy, LL& lp) const + { + try + { + pj_inv(this->m_prj, this->m_par, xy, lp); + return true; + } + catch(...) + { + return false; + } + } +}; + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + + +}}} // namespace boost::geometry::projections + + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_BASE_STATIC_HPP diff --git a/boost/geometry/srs/projections/impl/dms_parser.hpp b/boost/geometry/srs/projections/impl/dms_parser.hpp new file mode 100644 index 0000000000..bbecc9b1a2 --- /dev/null +++ b/boost/geometry/srs/projections/impl/dms_parser.hpp @@ -0,0 +1,278 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017, 2018. +// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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) + +#ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_DMS_PARSER_HPP +#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_DMS_PARSER_HPP + +// This file is totally revised from PROJ4 dmstor.c + +// PROJ4 is originally written by Gerald Evenden (then of the USGS) +// PROJ4 is maintained by Frank Warmerdam +// PROJ4 is converted to Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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 <string> + +#include <boost/static_assert.hpp> + +#if !defined(BOOST_GEOMETRY_NO_LEXICAL_CAST) +#include <boost/lexical_cast.hpp> +#endif // !defined(BOOST_GEOMETRY_NO_LEXICAL_CAST) + +#include <boost/algorithm/string.hpp> + +#include <boost/config.hpp> + +#include <boost/geometry/core/cs.hpp> + +#include <boost/geometry/util/math.hpp> + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + +template <typename T> +struct dms_result +{ + enum axis_selector {axis_lat = 1, axis_lon = 0}; + + private : + T m_angle; + axis_selector m_axis; + + public : + + explicit dms_result(T const& v, axis_selector ax) + : m_angle(v) + , m_axis(ax) + {} + + inline axis_selector axis() const { return m_axis; } + + inline T angle() const { return m_angle; } + + template <typename CH, typename TR> + inline friend std::basic_ostream<CH, TR>& operator<<(std::basic_ostream<CH, TR>& os, + const dms_result& d) + { + os << d.m_angle; + return os; + } + +}; + + +template <typename T + , bool as_radian = true + , char N = 'N', char E = 'E', char S = 'S', char W = 'W' // translatable + , char MIN = '\'', char SEC = '"' // other char's possible + , char D = 'D', char R = 'R' // degree sign might be small o + > +struct dms_parser +{ + + + // Question from Barend: can we compile-time select that it is case-sensitive/case-insensitive? + // We have to change the switch then -> specializations + + // For now: make it (compile-time) case sensitive + static const int diff = 'a' - 'A'; +#ifndef __GNUC__ + BOOST_STATIC_ASSERT((diff > 0)); // make sure we've the right assumption. GCC does not accept this here. +#endif + static const char n_alter = N <= 'Z' ? N + diff : N - diff; + static const char e_alter = E <= 'Z' ? E + diff : E - diff; + static const char s_alter = S <= 'Z' ? S + diff : S - diff; + static const char w_alter = W <= 'Z' ? W + diff : W - diff; + + static const char r_alter = R <= 'Z' ? R + diff : R - diff; + + // degree is normally D (proj4) but might be superscript o + // Note d_alter is not correct then, so map it to NULL now, guarded by the while + static const char d_alter = + ((D >= 'A' && D <= 'Z') || (D >= 'a' && D <= 'z')) ? (D <= 'Z' ? D + diff : D - diff) : '\0'; + + + struct dms_value + { + T dms[3]; + bool has_dms[3]; + + dms_value() +#ifndef BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX + : dms{0, 0, 0} + , has_dms{false, false, false} +#endif + { +#ifdef BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX + std::fill(dms, dms + 3, T(0)); + std::fill(has_dms, has_dms + 3, false); +#endif + } + }; + + + template <size_t I> + static inline void assign_dms(dms_value& dms, std::string& value, bool& has_value) + { +#if !defined(BOOST_GEOMETRY_NO_LEXICAL_CAST) + dms.dms[I] = boost::lexical_cast<T>(value.c_str()); +#else // !defined(BOOST_GEOMETRY_NO_LEXICAL_CAST) + dms.dms[I] = std::atof(value.c_str()); +#endif // !defined(BOOST_GEOMETRY_NO_LEXICAL_CAST) + dms.has_dms[I] = true; + has_value = false; + value.clear(); + } + + static inline void process(dms_value& dms, std::string& value, bool& has_value) + { + if (has_value) + { + // Assign last one, sequentially + if (! dms.has_dms[0]) assign_dms<0>(dms, value, has_value); + else if (! dms.has_dms[1]) assign_dms<1>(dms, value, has_value); + else if (! dms.has_dms[2]) assign_dms<2>(dms, value, has_value); + } + } + + dms_result<T> apply(std::string const& is) const + { + return apply(is.c_str()); + } + + dms_result<T> apply(const char* is) const + { + dms_value dms; + bool has_value = false; + std::string value; + + T factor = 1.0; // + denotes N/E values, -1 denotes S/W values + typename dms_result<T>::axis_selector axis = dms_result<T>::axis_lon; // true denotes N/S values + bool in_radian = false; // true denotes values as "0.1R" + + while(*is) + { + switch(*is) + { + case '-' : + if (! has_value && ! dms.has_dms[0]) + { + factor = -factor; + } + break; + case N : + case n_alter : + axis = dms_result<T>::axis_lat; + break; + case S : + case s_alter : + axis = dms_result<T>::axis_lat; + factor = -factor; + break; + case E : + case e_alter : + axis = dms_result<T>::axis_lon; + break; + case W : + case w_alter : + axis = dms_result<T>::axis_lon; + factor = -factor; + break; + case D : + case d_alter : + if (! dms.has_dms[0] && has_value) + { + assign_dms<0>(dms, value, has_value); + } + break; + case R : + case r_alter : + if (! dms.has_dms[0] && has_value) + { + // specified value is in radian! + in_radian = true; + assign_dms<0>(dms, value, has_value); + } + break; + case MIN: + if (! dms.has_dms[1] && has_value) + { + assign_dms<1>(dms, value, has_value); + } + break; + case SEC : + if (! dms.has_dms[2] && has_value) + { + assign_dms<2>(dms, value, has_value); + } + break; + case ' ' : + case '\t' : + case '\n' : + process(dms, value, has_value); + break; + default : + value += *is; + has_value = true; + break; + } + is++; + } + + // Assign last one, if any + process(dms, value, has_value); + + T const d2r = math::d2r<T>(); + T const r2d = math::r2d<T>(); + + return dms_result<T>(factor * + (in_radian && as_radian + ? dms.dms[0] + : in_radian && ! as_radian + ? dms.dms[0] * r2d + : ! in_radian && as_radian + ? dms.dms[0] * d2r + dms.dms[1] * d2r / 60.0 + dms.dms[2] * d2r / 3600.0 + : dms.dms[0] + dms.dms[1] / 60.0 + dms.dms[2] / 3600.0) + , axis); + } +}; + + +} // namespace detail + + +}}} // namespace boost::geometry::projections + + +#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_DMS_PARSER_HPP diff --git a/boost/geometry/srs/projections/impl/factory_entry.hpp b/boost/geometry/srs/projections/impl/factory_entry.hpp new file mode 100644 index 0000000000..7c587af5b6 --- /dev/null +++ b/boost/geometry/srs/projections/impl/factory_entry.hpp @@ -0,0 +1,51 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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) + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_FACTORY_ENTRY_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_FACTORY_ENTRY_HPP + +#include <string> + +#include <boost/geometry/srs/projections/impl/base_dynamic.hpp> + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + +// forward declaration needed by some projections +template <typename CT, typename Parameters> +class factory; + +template <typename CT, typename P> +class factory_entry +{ +public: + + virtual ~factory_entry() {} + virtual base_v<CT, P>* create_new(P const& par) const = 0; +}; + +template <typename CT, typename P> +class base_factory +{ +public: + + virtual ~base_factory() {} + virtual void add_to_factory(std::string const& name, factory_entry<CT, P>* sub) = 0; +}; + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_FACTORY_ENTRY_HPP diff --git a/boost/geometry/srs/projections/impl/function_overloads.hpp b/boost/geometry/srs/projections/impl/function_overloads.hpp new file mode 100644 index 0000000000..43d33c6bea --- /dev/null +++ b/boost/geometry/srs/projections/impl/function_overloads.hpp @@ -0,0 +1,46 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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) + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_FUNCTION_OVERLOADS_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_FUNCTION_OVERLOADS_HPP + +#include <cmath> + +namespace boost { namespace geometry { namespace projections +{ + +// Functions to resolve ambiguity when compiling with coordinates of different types +/* +template <typename T> +inline T atan2(T const& a, T const& b) +{ + using std::atan2; + return atan2(a, b); +} +template <typename T> +inline T pow(T const& a, T const& b) +{ + using std::pow; + return pow(a, b); +} +*/ + +template <typename T> +inline int int_floor(T const& f) +{ + using std::floor; + return int(floor(f)); +} + +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_FUNCTION_OVERLOADS_HPP diff --git a/boost/geometry/srs/projections/impl/geocent.hpp b/boost/geometry/srs/projections/impl/geocent.hpp new file mode 100644 index 0000000000..8ae2f393ce --- /dev/null +++ b/boost/geometry/srs/projections/impl/geocent.hpp @@ -0,0 +1,487 @@ +// Boost.Geometry +// This file is manually converted from PROJ4 + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 +// This file was converted to Geometry Library by Adam Wulkiewicz + +// Original copyright notice: + +/***************************************************************************/ +/* RSC IDENTIFIER: GEOCENTRIC + * + * ABSTRACT + * + * This component provides conversions between Geodetic coordinates (latitude, + * longitude in radians and height in meters) and Geocentric coordinates + * (X, Y, Z) in meters. + * + * ERROR HANDLING + * + * This component checks parameters for valid values. If an invalid value + * is found, the error code is combined with the current error code using + * the bitwise or. This combining allows multiple error codes to be + * returned. The possible error codes are: + * + * GEOCENT_NO_ERROR : No errors occurred in function + * GEOCENT_LAT_ERROR : Latitude out of valid range + * (-90 to 90 degrees) + * GEOCENT_LON_ERROR : Longitude out of valid range + * (-180 to 360 degrees) + * GEOCENT_A_ERROR : Semi-major axis lessthan or equal to zero + * GEOCENT_B_ERROR : Semi-minor axis lessthan or equal to zero + * GEOCENT_A_LESS_B_ERROR : Semi-major axis less than semi-minor axis + * + * + * REUSE NOTES + * + * GEOCENTRIC is intended for reuse by any application that performs + * coordinate conversions between geodetic coordinates and geocentric + * coordinates. + * + * + * REFERENCES + * + * An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion, + * Ralph Toms, February 1996 UCRL-JC-123138. + * + * Further information on GEOCENTRIC can be found in the Reuse Manual. + * + * GEOCENTRIC originated from : U.S. Army Topographic Engineering Center + * Geospatial Information Division + * 7701 Telegraph Road + * Alexandria, VA 22310-3864 + * + * LICENSES + * + * None apply to this component. + * + * RESTRICTIONS + * + * GEOCENTRIC has no restrictions. + * + * ENVIRONMENT + * + * GEOCENTRIC was tested and certified in the following environments: + * + * 1. Solaris 2.5 with GCC version 2.8.1 + * 2. Windows 95 with MS Visual C++ version 6 + * + * MODIFICATIONS + * + * Date Description + * ---- ----------- + * 25-02-97 Original Code + * + */ + + +#ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_GEOCENT_HPP +#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_GEOCENT_HPP + + +#include <boost/geometry/util/math.hpp> + + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + +/***************************************************************************/ +/* + * DEFINES + */ +static const long GEOCENT_NO_ERROR = 0x0000; +static const long GEOCENT_LAT_ERROR = 0x0001; +static const long GEOCENT_LON_ERROR = 0x0002; +static const long GEOCENT_A_ERROR = 0x0004; +static const long GEOCENT_B_ERROR = 0x0008; +static const long GEOCENT_A_LESS_B_ERROR = 0x0010; + +template <typename T> +struct GeocentricInfo +{ + T Geocent_a; /* Semi-major axis of ellipsoid in meters */ + T Geocent_b; /* Semi-minor axis of ellipsoid */ + T Geocent_a2; /* Square of semi-major axis */ + T Geocent_b2; /* Square of semi-minor axis */ + T Geocent_e2; /* Eccentricity squared */ + T Geocent_ep2; /* 2nd eccentricity squared */ +}; + +template <typename T> +inline T COS_67P5() +{ + /*return 0.38268343236508977*/; + return cos(T(67.5) * math::d2r<T>()); /* cosine of 67.5 degrees */ +} +template <typename T> +inline T AD_C() +{ + return 1.0026000; /* Toms region 1 constant */ +} + + +/***************************************************************************/ +/* + * FUNCTIONS + */ + +template <typename T> +inline long pj_Set_Geocentric_Parameters (GeocentricInfo<T> & gi, T const& a, T const& b) + +{ /* BEGIN Set_Geocentric_Parameters */ +/* + * The function Set_Geocentric_Parameters receives the ellipsoid parameters + * as inputs and sets the corresponding state variables. + * + * a : Semi-major axis, in meters. (input) + * b : Semi-minor axis, in meters. (input) + */ + long Error_Code = GEOCENT_NO_ERROR; + + if (a <= 0.0) + Error_Code |= GEOCENT_A_ERROR; + if (b <= 0.0) + Error_Code |= GEOCENT_B_ERROR; + if (a < b) + Error_Code |= GEOCENT_A_LESS_B_ERROR; + if (!Error_Code) + { + gi.Geocent_a = a; + gi.Geocent_b = b; + gi.Geocent_a2 = a * a; + gi.Geocent_b2 = b * b; + gi.Geocent_e2 = (gi.Geocent_a2 - gi.Geocent_b2) / gi.Geocent_a2; + gi.Geocent_ep2 = (gi.Geocent_a2 - gi.Geocent_b2) / gi.Geocent_b2; + } + return (Error_Code); +} /* END OF Set_Geocentric_Parameters */ + + +template <typename T> +inline void pj_Get_Geocentric_Parameters (GeocentricInfo<T> const& gi, + T & a, + T & b) +{ /* BEGIN Get_Geocentric_Parameters */ +/* + * The function Get_Geocentric_Parameters returns the ellipsoid parameters + * to be used in geocentric coordinate conversions. + * + * a : Semi-major axis, in meters. (output) + * b : Semi-minor axis, in meters. (output) + */ + + a = gi.Geocent_a; + b = gi.Geocent_b; +} /* END OF Get_Geocentric_Parameters */ + + +template <typename T> +inline long pj_Convert_Geodetic_To_Geocentric (GeocentricInfo<T> const& gi, + T Longitude, T Latitude, T Height, + T & X, T & Y, T & Z) +{ /* BEGIN Convert_Geodetic_To_Geocentric */ +/* + * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates + * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z), + * according to the current ellipsoid parameters. + * + * Latitude : Geodetic latitude in radians (input) + * Longitude : Geodetic longitude in radians (input) + * Height : Geodetic height, in meters (input) + * X : Calculated Geocentric X coordinate, in meters (output) + * Y : Calculated Geocentric Y coordinate, in meters (output) + * Z : Calculated Geocentric Z coordinate, in meters (output) + * + */ + long Error_Code = GEOCENT_NO_ERROR; + T Rn; /* Earth radius at location */ + T Sin_Lat; /* sin(Latitude) */ + T Sin2_Lat; /* Square of sin(Latitude) */ + T Cos_Lat; /* cos(Latitude) */ + + static const T PI = math::pi<T>(); + static const T PI_OVER_2 = math::half_pi<T>(); + + /* + ** Don't blow up if Latitude is just a little out of the value + ** range as it may just be a rounding issue. Also removed longitude + ** test, it should be wrapped by cos() and sin(). NFW for PROJ.4, Sep/2001. + */ + if( Latitude < -PI_OVER_2 && Latitude > -1.001 * PI_OVER_2 ) + Latitude = -PI_OVER_2; + else if( Latitude > PI_OVER_2 && Latitude < 1.001 * PI_OVER_2 ) + Latitude = PI_OVER_2; + else if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2)) + { /* Latitude out of range */ + Error_Code |= GEOCENT_LAT_ERROR; + } + + if (!Error_Code) + { /* no errors */ + if (Longitude > PI) + Longitude -= (2*PI); + Sin_Lat = sin(Latitude); + Cos_Lat = cos(Latitude); + Sin2_Lat = Sin_Lat * Sin_Lat; + Rn = gi.Geocent_a / (sqrt(1.0e0 - gi.Geocent_e2 * Sin2_Lat)); + X = (Rn + Height) * Cos_Lat * cos(Longitude); + Y = (Rn + Height) * Cos_Lat * sin(Longitude); + Z = ((Rn * (1 - gi.Geocent_e2)) + Height) * Sin_Lat; + } + return (Error_Code); +} /* END OF Convert_Geodetic_To_Geocentric */ + +/* + * The function Convert_Geocentric_To_Geodetic converts geocentric + * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude, + * and height), according to the current ellipsoid parameters. + * + * X : Geocentric X coordinate, in meters. (input) + * Y : Geocentric Y coordinate, in meters. (input) + * Z : Geocentric Z coordinate, in meters. (input) + * Latitude : Calculated latitude value in radians. (output) + * Longitude : Calculated longitude value in radians. (output) + * Height : Calculated height value, in meters. (output) + */ + +#define BOOST_GEOMETRY_PROJECTIONS_USE_ITERATIVE_METHOD + +template <typename T> +inline void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo<T> const& gi, + T X, T Y, T Z, + T & Longitude, T & Latitude, T & Height) +{ /* BEGIN Convert_Geocentric_To_Geodetic */ + + static const T PI_OVER_2 = math::half_pi<T>(); + +#if !defined(BOOST_GEOMETRY_PROJECTIONS_USE_ITERATIVE_METHOD) + + static const T COS_67P5 = detail::COS_67P5<T>(); + static const T AD_C = detail::AD_C<T>(); + +/* + * The method used here is derived from 'An Improved Algorithm for + * Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996 + */ + +/* Note: Variable names follow the notation used in Toms, Feb 1996 */ + + T W; /* distance from Z axis */ + T W2; /* square of distance from Z axis */ + T T0; /* initial estimate of vertical component */ + T T1; /* corrected estimate of vertical component */ + T S0; /* initial estimate of horizontal component */ + T S1; /* corrected estimate of horizontal component */ + T Sin_B0; /* sin(B0), B0 is estimate of Bowring aux variable */ + T Sin3_B0; /* cube of sin(B0) */ + T Cos_B0; /* cos(B0) */ + T Sin_p1; /* sin(phi1), phi1 is estimated latitude */ + T Cos_p1; /* cos(phi1) */ + T Rn; /* Earth radius at location */ + T Sum; /* numerator of cos(phi1) */ + bool At_Pole; /* indicates location is in polar region */ + + At_Pole = false; + if (X != 0.0) + { + Longitude = atan2(Y,X); + } + else + { + if (Y > 0) + { + Longitude = PI_OVER_2; + } + else if (Y < 0) + { + Longitude = -PI_OVER_2; + } + else + { + At_Pole = true; + Longitude = 0.0; + if (Z > 0.0) + { /* north pole */ + Latitude = PI_OVER_2; + } + else if (Z < 0.0) + { /* south pole */ + Latitude = -PI_OVER_2; + } + else + { /* center of earth */ + Latitude = PI_OVER_2; + Height = -Geocent_b; + return; + } + } + } + W2 = X*X + Y*Y; + W = sqrt(W2); + T0 = Z * AD_C; + S0 = sqrt(T0 * T0 + W2); + Sin_B0 = T0 / S0; + Cos_B0 = W / S0; + Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0; + T1 = Z + gi.Geocent_b * gi.Geocent_ep2 * Sin3_B0; + Sum = W - gi.Geocent_a * gi.Geocent_e2 * Cos_B0 * Cos_B0 * Cos_B0; + S1 = sqrt(T1*T1 + Sum * Sum); + Sin_p1 = T1 / S1; + Cos_p1 = Sum / S1; + Rn = gi.Geocent_a / sqrt(1.0 - gi.Geocent_e2 * Sin_p1 * Sin_p1); + if (Cos_p1 >= COS_67P5) + { + Height = W / Cos_p1 - Rn; + } + else if (Cos_p1 <= -COS_67P5) + { + Height = W / -Cos_p1 - Rn; + } + else + { + Height = Z / Sin_p1 + Rn * (gi.Geocent_e2 - 1.0); + } + if (At_Pole == false) + { + Latitude = atan(Sin_p1 / Cos_p1); + } +#else /* defined(BOOST_GEOMETRY_PROJECTIONS_USE_ITERATIVE_METHOD) */ +/* +* Reference... +* ============ +* Wenzel, H.-G.(1985): Hochauflösende Kugelfunktionsmodelle für +* das Gravitationspotential der Erde. Wiss. Arb. Univ. Hannover +* Nr. 137, p. 130-131. + +* Programmed by GGA- Leibniz-Institute of Applied Geophysics +* Stilleweg 2 +* D-30655 Hannover +* Federal Republic of Germany +* Internet: www.gga-hannover.de +* +* Hannover, March 1999, April 2004. +* see also: comments in statements +* remarks: +* Mathematically exact and because of symmetry of rotation-ellipsoid, +* each point (X,Y,Z) has at least two solutions (Latitude1,Longitude1,Height1) and +* (Latitude2,Longitude2,Height2). Is point=(0.,0.,Z) (P=0.), so you get even +* four solutions, every two symmetrical to the semi-minor axis. +* Here Height1 and Height2 have at least a difference in order of +* radius of curvature (e.g. (0,0,b)=> (90.,0.,0.) or (-90.,0.,-2b); +* (a+100.)*(sqrt(2.)/2.,sqrt(2.)/2.,0.) => (0.,45.,100.) or +* (0.,225.,-(2a+100.))). +* The algorithm always computes (Latitude,Longitude) with smallest |Height|. +* For normal computations, that means |Height|<10000.m, algorithm normally +* converges after to 2-3 steps!!! +* But if |Height| has the amount of length of ellipsoid's axis +* (e.g. -6300000.m), algorithm needs about 15 steps. +*/ + +/* local definitions and variables */ +/* end-criterium of loop, accuracy of sin(Latitude) */ +static const T genau = 1.E-12; +static const T genau2 = (genau*genau); +static const int maxiter = 30; + + T P; /* distance between semi-minor axis and location */ + T RR; /* distance between center and location */ + T CT; /* sin of geocentric latitude */ + T ST; /* cos of geocentric latitude */ + T RX; + T RK; + T RN; /* Earth radius at location */ + T CPHI0; /* cos of start or old geodetic latitude in iterations */ + T SPHI0; /* sin of start or old geodetic latitude in iterations */ + T CPHI; /* cos of searched geodetic latitude */ + T SPHI; /* sin of searched geodetic latitude */ + T SDPHI; /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */ + int iter; /* # of continuous iteration, max. 30 is always enough (s.a.) */ + + P = sqrt(X*X+Y*Y); + RR = sqrt(X*X+Y*Y+Z*Z); + +/* special cases for latitude and longitude */ + if (P/gi.Geocent_a < genau) { + +/* special case, if P=0. (X=0., Y=0.) */ + Longitude = 0.; + +/* if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis + * of ellipsoid (=center of mass), Latitude becomes PI/2 */ + if (RR/gi.Geocent_a < genau) { + Latitude = PI_OVER_2; + Height = -gi.Geocent_b; + return ; + + } + } + else { +/* ellipsoidal (geodetic) longitude + * interval: -PI < Longitude <= +PI */ + Longitude=atan2(Y,X); + } + +/* -------------------------------------------------------------- + * Following iterative algorithm was developed by + * "Institut für Erdmessung", University of Hannover, July 1988. + * Internet: www.ife.uni-hannover.de + * Iterative computation of CPHI,SPHI and Height. + * Iteration of CPHI and SPHI to 10**-12 radian resp. + * 2*10**-7 arcsec. + * -------------------------------------------------------------- + */ + CT = Z/RR; + ST = P/RR; + RX = 1.0/sqrt(1.0-gi.Geocent_e2*(2.0-gi.Geocent_e2)*ST*ST); + CPHI0 = ST*(1.0-gi.Geocent_e2)*RX; + SPHI0 = CT*RX; + iter = 0; + +/* loop to find sin(Latitude) resp. Latitude + * until |sin(Latitude(iter)-Latitude(iter-1))| < genau */ + do + { + iter++; + RN = gi.Geocent_a/sqrt(1.0-gi.Geocent_e2*SPHI0*SPHI0); + +/* ellipsoidal (geodetic) height */ + Height = P*CPHI0+Z*SPHI0-RN*(1.0-gi.Geocent_e2*SPHI0*SPHI0); + + RK = gi.Geocent_e2*RN/(RN+Height); + RX = 1.0/sqrt(1.0-RK*(2.0-RK)*ST*ST); + CPHI = ST*(1.0-RK)*RX; + SPHI = CT*RX; + SDPHI = SPHI*CPHI0-CPHI*SPHI0; + CPHI0 = CPHI; + SPHI0 = SPHI; + } + while (SDPHI*SDPHI > genau2 && iter < maxiter); + +/* ellipsoidal (geodetic) latitude */ + Latitude=atan(SPHI/fabs(CPHI)); + + return; +#endif /* defined(BOOST_GEOMETRY_PROJECTIONS_USE_ITERATIVE_METHOD) */ +} /* END OF Convert_Geocentric_To_Geodetic */ + + +} // namespace detail + + +}}} // namespace boost::geometry::projections + + +#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_GEOCENT_HPP diff --git a/boost/geometry/srs/projections/impl/pj_auth.hpp b/boost/geometry/srs/projections/impl/pj_auth.hpp new file mode 100644 index 0000000000..899d0b25b1 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_auth.hpp @@ -0,0 +1,95 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_AUTH_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_AUTH_HPP + +#include <cassert> +#include <cmath> + +#include <boost/geometry/core/assert.hpp> + +namespace boost { namespace geometry { namespace projections { + +namespace detail { + +static const int APA_SIZE = 3; + +/* determine latitude from authalic latitude */ +template <typename T> +inline bool pj_authset(T const& es, T* APA) +{ + BOOST_GEOMETRY_ASSERT(0 != APA); + + static const T P00 = .33333333333333333333; + static const T P01 = .17222222222222222222; + static const T P02 = .10257936507936507936; + static const T P10 = .06388888888888888888; + static const T P11 = .06640211640211640211; + static const T P20 = .01641501294219154443; + + T t = 0; + + // if (APA = (double *)pj_malloc(APA_SIZE * sizeof(double))) + { + APA[0] = es * P00; + t = es * es; + APA[0] += t * P01; + APA[1] = t * P10; + t *= es; + APA[0] += t * P02; + APA[1] += t * P11; + APA[2] = t * P20; + } + return true; +} + +template <typename T> +inline T pj_authlat(T const& beta, const T* APA) +{ + BOOST_GEOMETRY_ASSERT(0 != APA); + + T const t = beta + beta; + + return(beta + APA[0] * sin(t) + APA[1] * sin(t + t) + APA[2] * sin(t + t + t)); +} + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_AUTH_HPP diff --git a/boost/geometry/srs/projections/impl/pj_datum_set.hpp b/boost/geometry/srs/projections/impl/pj_datum_set.hpp new file mode 100644 index 0000000000..5301dc7cfe --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_datum_set.hpp @@ -0,0 +1,213 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017, 2018. +// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_DATUM_SET_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_DATUM_SET_HPP + + +#include <string> +#include <vector> + +#include <boost/algorithm/string.hpp> + +#include <boost/geometry/srs/projections/exception.hpp> +#include <boost/geometry/srs/projections/impl/projects.hpp> +#include <boost/geometry/srs/projections/impl/pj_datums.hpp> +#include <boost/geometry/srs/projections/impl/pj_param.hpp> +#include <boost/geometry/srs/projections/par4.hpp> +#include <boost/geometry/srs/projections/proj4.hpp> + + +namespace boost { namespace geometry { namespace projections { + +namespace detail { + + +/* SEC_TO_RAD = Pi/180/3600 */ +template <typename T> +inline T SEC_TO_RAD() { return 4.84813681109535993589914102357e-6; } + +template <typename BGParams, typename T> +inline void pj_datum_add_defn(BGParams const& , std::vector<pvalue<T> >& pvalues) +{ + /* -------------------------------------------------------------------- */ + /* Is there a datum definition in the parameter list? If so, */ + /* add the defining values to the parameter list. Note that */ + /* this will append the ellipse definition as well as the */ + /* towgs84= and related parameters. It should also be pointed */ + /* out that the addition is permanent rather than temporary */ + /* like most other keyword expansion so that the ellipse */ + /* definition will last into the pj_ell_set() function called */ + /* after this one. */ + /* -------------------------------------------------------------------- */ + std::string name = pj_param(pvalues, "sdatum").s; + if(! name.empty()) + { + /* find the datum definition */ + const int n = sizeof(pj_datums) / sizeof(pj_datums[0]); + int index = -1; + for (int i = 0; i < n && index == -1; i++) + { + if(pj_datums[i].id == name) + { + index = i; + } + } + + if (index == -1) + { + BOOST_THROW_EXCEPTION( projection_exception(-9) ); + } + + if(! pj_datums[index].ellipse_id.empty()) + { + std::string entry("ellps="); + entry +=pj_datums[index].ellipse_id; + pvalues.push_back(pj_mkparam<T>(entry)); + } + + if(! pj_datums[index].defn.empty()) + { + pvalues.push_back(pj_mkparam<T>(pj_datums[index].defn)); + } + } +} + +template <BOOST_GEOMETRY_PROJECTIONS_DETAIL_TYPENAME_PX, typename T> +inline void pj_datum_add_defn(srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAIL_PX> const& /*bg_params*/, + std::vector<pvalue<T> >& pvalues) +{ + typedef srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAIL_PX> bg_parameters_type; + typedef typename srs::par4::detail::tuples_find_if + < + bg_parameters_type, + //srs::par4::detail::is_datum + srs::par4::detail::is_param_t<srs::par4::datum>::pred + >::type datum_type; + typedef typename srs::par4::detail::datum_traits + < + datum_type + > datum_traits; + + // is unknown if datum parameter found but traits are not specialized + static const bool not_set_or_known = boost::is_same<datum_type, void>::value + || ! boost::is_same<typename datum_traits::ellps_type, void>::value; + BOOST_MPL_ASSERT_MSG((not_set_or_known), UNKNOWN_DATUM, (bg_parameters_type)); + + std::string defn = datum_traits::definition(); + + if (! defn.empty()) + { + pvalues.push_back(pj_mkparam<T>(defn)); + } +} + +/************************************************************************/ +/* pj_datum_set() */ +/************************************************************************/ + +template <typename BGParams, typename T> +inline void pj_datum_set(BGParams const& bg_params, std::vector<pvalue<T> >& pvalues, parameters<T>& projdef) +{ + static const T SEC_TO_RAD = detail::SEC_TO_RAD<T>(); + + projdef.datum_type = PJD_UNKNOWN; + + pj_datum_add_defn(bg_params, pvalues); + +/* -------------------------------------------------------------------- */ +/* Check for nadgrids parameter. */ +/* -------------------------------------------------------------------- */ + std::string nadgrids = pj_param(pvalues, "snadgrids").s; + std::string towgs84 = pj_param(pvalues, "stowgs84").s; + if(! nadgrids.empty()) + { + /* We don't actually save the value separately. It will continue + to exist int he param list for use in pj_apply_gridshift.c */ + + projdef.datum_type = PJD_GRIDSHIFT; + } + +/* -------------------------------------------------------------------- */ +/* Check for towgs84 parameter. */ +/* -------------------------------------------------------------------- */ + else if(! towgs84.empty()) + { + int parm_count = 0; + + int n = sizeof(projdef.datum_params) / sizeof(projdef.datum_params[0]); + + /* parse out the pvalues */ + std::vector<std::string> parm; + boost::split(parm, towgs84, boost::is_any_of(" ,")); + for (std::vector<std::string>::const_iterator it = parm.begin(); + it != parm.end() && parm_count < n; + ++it) + { + projdef.datum_params[parm_count++] = atof(it->c_str()); + } + + if( projdef.datum_params[3] != 0.0 + || projdef.datum_params[4] != 0.0 + || projdef.datum_params[5] != 0.0 + || projdef.datum_params[6] != 0.0 ) + { + projdef.datum_type = PJD_7PARAM; + + /* transform from arc seconds to radians */ + projdef.datum_params[3] *= SEC_TO_RAD; + projdef.datum_params[4] *= SEC_TO_RAD; + projdef.datum_params[5] *= SEC_TO_RAD; + /* transform from parts per million to scaling factor */ + projdef.datum_params[6] = + (projdef.datum_params[6]/1000000.0) + 1; + } + else + { + projdef.datum_type = PJD_3PARAM; + } + + /* Note that pj_init() will later switch datum_type to + PJD_WGS84 if shifts are all zero, and ellipsoid is WGS84 or GRS80 */ + } +} + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_DATUM_SET_HPP diff --git a/boost/geometry/srs/projections/impl/pj_datums.hpp b/boost/geometry/srs/projections/impl/pj_datums.hpp new file mode 100644 index 0000000000..55da24a2ca --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_datums.hpp @@ -0,0 +1,112 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_DATUMS_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_DATUMS_HPP + +#include <boost/geometry/srs/projections/impl/projects.hpp> + +namespace boost { namespace geometry { namespace projections { + +namespace detail { + +/* + * The ellipse code must match one from pj_ellps.c. The datum id should + * be kept to 12 characters or less if possible. Use the official OGC + * datum name for the comments if available. + */ + +static const PJ_DATUMS pj_datums[] = +{ + /* id definition ellipse comments */ + /* -- ---------- ------- -------- */ + {"WGS84", "towgs84=0,0,0", + "WGS84", ""}, + + {"GGRS87", "towgs84=-199.87,74.79,246.62", + "GRS80", "Greek_Geodetic_Reference_System_1987"}, + + {"NAD83", "towgs84=0,0,0", + "GRS80", "North_American_Datum_1983"}, + + {"NAD27", "nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat", + "clrk66", "North_American_Datum_1927"}, + + {"potsdam", "towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7", + "bessel", "Potsdam Rauenberg 1950 DHDN"}, + + {"carthage", "towgs84=-263.0,6.0,431.0", + "clrk80ign", "Carthage 1934 Tunisia"}, + + {"hermannskogel", "towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232", + "bessel", "Hermannskogel"}, + + {"ire65", "towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15", + "mod_airy", "Ireland 1965"}, + + {"nzgd49", "towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993", + "intl", "New Zealand Geodetic Datum 1949"}, + + {"OSGB36", "towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894", + "airy", "Airy 1830"} +}; + + +static const PJ_PRIME_MERIDIANS pj_prime_meridians[] = +{ + /* id definition */ + /* -- ---------- */ + { "greenwich", "0dE" }, + { "lisbon", "9d07'54.862\"W" }, + { "paris", "2d20'14.025\"E" }, + { "bogota", "74d04'51.3\"W" }, + { "madrid", "3d41'16.58\"W" }, + { "rome", "12d27'8.4\"E" }, + { "bern", "7d26'22.5\"E" }, + { "jakarta", "106d48'27.79\"E" }, + { "ferro", "17d40'W" }, + { "brussels", "4d22'4.71\"E" }, + { "stockholm", "18d3'29.8\"E" }, + { "athens", "23d42'58.815\"E" }, + { "oslo", "10d43'22.5\"E" } +}; + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_DATUMS_HPP diff --git a/boost/geometry/srs/projections/impl/pj_ell_set.hpp b/boost/geometry/srs/projections/impl/pj_ell_set.hpp new file mode 100644 index 0000000000..6f5f14b780 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_ell_set.hpp @@ -0,0 +1,205 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017, 2018. +// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELL_SET_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELL_SET_HPP + +#include <string> +#include <vector> + +#include <boost/geometry/formulas/eccentricity_sqr.hpp> +#include <boost/geometry/util/math.hpp> + +#include <boost/geometry/srs/projections/impl/pj_ellps.hpp> +#include <boost/geometry/srs/projections/impl/pj_param.hpp> +#include <boost/geometry/srs/projections/proj4.hpp> + + +namespace boost { namespace geometry { namespace projections { + +namespace detail { + +/* set ellipsoid parameters a and es */ +template <typename T> +inline T SIXTH() { return .1666666666666666667; } /* 1/6 */ +template <typename T> +inline T RA4() { return .04722222222222222222; } /* 17/360 */ +template <typename T> +inline T RA6() { return .02215608465608465608; } /* 67/3024 */ +template <typename T> +inline T RV4() { return .06944444444444444444; } /* 5/72 */ +template <typename T> +inline T RV6() { return .04243827160493827160; } /* 55/1296 */ + +/* initialize geographic shape parameters */ +template <typename BGParams, typename T> +inline void pj_ell_set(BGParams const& /*bg_params*/, std::vector<pvalue<T> >& parameters, T &a, T &es) +{ + T b = 0.0; + T e = 0.0; + std::string name; + + /* check for varying forms of ellipsoid input */ + a = es = 0.; + + /* R takes precedence */ + if (pj_param(parameters, "tR").i) + a = pj_param(parameters, "dR").f; + else { /* probable elliptical figure */ + + /* check if ellps present and temporarily append its values to pl */ + name = pj_param(parameters, "sellps").s; + if (! name.empty()) + { + const int n = sizeof(pj_ellps) / sizeof(pj_ellps[0]); + int index = -1; + for (int i = 0; i < n && index == -1; i++) + { + if(pj_ellps[i].id == name) + { + index = i; + } + } + + if (index == -1) { + BOOST_THROW_EXCEPTION( projection_exception(-9) ); + } + + parameters.push_back(pj_mkparam<T>(pj_ellps[index].major)); + parameters.push_back(pj_mkparam<T>(pj_ellps[index].ell)); + } + a = pj_param(parameters, "da").f; + if (pj_param(parameters, "tes").i) /* eccentricity squared */ + es = pj_param(parameters, "des").f; + else if (pj_param(parameters, "te").i) { /* eccentricity */ + e = pj_param(parameters, "de").f; + es = e * e; + } else if (pj_param(parameters, "trf").i) { /* recip flattening */ + es = pj_param(parameters, "drf").f; + if (!es) { + BOOST_THROW_EXCEPTION( projection_exception(-10) ); + } + es = 1./ es; + es = es * (2. - es); + } else if (pj_param(parameters, "tf").i) { /* flattening */ + es = pj_param(parameters, "df").f; + es = es * (2. - es); + } else if (pj_param(parameters, "tb").i) { /* minor axis */ + b = pj_param(parameters, "db").f; + es = 1. - (b * b) / (a * a); + } /* else es == 0. and sphere of radius a */ + if (!b) + b = a * sqrt(1. - es); + /* following options turn ellipsoid into equivalent sphere */ + if (pj_param(parameters, "bR_A").i) { /* sphere--area of ellipsoid */ + a *= 1. - es * (SIXTH<T>() + es * (RA4<T>() + es * RA6<T>())); + es = 0.; + } else if (pj_param(parameters, "bR_V").i) { /* sphere--vol. of ellipsoid */ + a *= 1. - es * (SIXTH<T>() + es * (RV4<T>() + es * RV6<T>())); + es = 0.; + } else if (pj_param(parameters, "bR_a").i) { /* sphere--arithmetic mean */ + a = .5 * (a + b); + es = 0.; + } else if (pj_param(parameters, "bR_g").i) { /* sphere--geometric mean */ + a = sqrt(a * b); + es = 0.; + } else if (pj_param(parameters, "bR_h").i) { /* sphere--harmonic mean */ + a = 2. * a * b / (a + b); + es = 0.; + } else { + int i = pj_param(parameters, "tR_lat_a").i; + if (i || /* sphere--arith. */ + pj_param(parameters, "tR_lat_g").i) { /* or geom. mean at latitude */ + T tmp; + + tmp = sin(pj_param(parameters, i ? "rR_lat_a" : "rR_lat_g").f); + if (geometry::math::abs(tmp) > geometry::math::half_pi<T>()) { + BOOST_THROW_EXCEPTION( projection_exception(-11) ); + } + tmp = 1. - es * tmp * tmp; + a *= i ? .5 * (1. - es + tmp) / ( tmp * sqrt(tmp)) : + sqrt(1. - es) / tmp; + es = 0.; + } + } + } + + /* some remaining checks */ + if (es < 0.) { + BOOST_THROW_EXCEPTION( projection_exception(-12) ); + } + if (a <= 0.) { + BOOST_THROW_EXCEPTION( projection_exception(-13) ); + } +} + +template <BOOST_GEOMETRY_PROJECTIONS_DETAIL_TYPENAME_PX, typename T> +inline void pj_ell_set(srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAIL_PX> const& bg_params, + std::vector<pvalue<T> >& /*parameters*/, T &a, T &es) +{ + typedef srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAIL_PX> static_parameters_type; + typedef typename srs::par4::detail::pick_ellps + < + static_parameters_type + > pick_ellps; + + typename pick_ellps::model_type model = pick_ellps::model(bg_params); + + a = geometry::get_radius<0>(model); + T b = geometry::get_radius<2>(model); + es = 0.; + if (a != b) + { + es = formula::eccentricity_sqr<T>(model); + + // Ignore all other parameters passed in string, at least for now + } + + /* some remaining checks */ + if (es < 0.) { + BOOST_THROW_EXCEPTION( projection_exception(-12) ); + } + if (a <= 0.) { + BOOST_THROW_EXCEPTION( projection_exception(-13) ); + } +} + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELL_SET_HPP diff --git a/boost/geometry/srs/projections/impl/pj_ellps.hpp b/boost/geometry/srs/projections/impl/pj_ellps.hpp new file mode 100644 index 0000000000..586802778c --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_ellps.hpp @@ -0,0 +1,98 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELLPS_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELLPS_HPP + +#include <boost/geometry/srs/projections/impl/projects.hpp> + +namespace boost { namespace geometry { namespace projections { + +namespace detail { + +static const PJ_ELLPS pj_ellps[] = +{ + {"MERIT", "a=6378137.0", "rf=298.257", "MERIT 1983"}, + {"SGS85", "a=6378136.0", "rf=298.257", "Soviet Geodetic System 85"}, + {"GRS80", "a=6378137.0", "rf=298.257222101", "GRS 1980(IUGG, 1980)"}, + {"IAU76", "a=6378140.0", "rf=298.257", "IAU 1976"}, + {"airy", "a=6377563.396", "b=6356256.910", "Airy 1830"}, + {"APL4.9", "a=6378137.0.", "rf=298.25", "Appl. Physics. 1965"}, + {"NWL9D", "a=6378145.0.", "rf=298.25", "Naval Weapons Lab., 1965"}, + {"mod_airy", "a=6377340.189", "b=6356034.446", "Modified Airy"}, + {"andrae", "a=6377104.43", "rf=300.0", "Andrae 1876 (Den., Iclnd.)"}, + {"aust_SA", "a=6378160.0", "rf=298.25", "Australian Natl & S. Amer. 1969"}, + {"GRS67", "a=6378160.0", "rf=298.2471674270", "GRS 67(IUGG 1967)"}, + {"bessel", "a=6377397.155", "rf=299.1528128", "Bessel 1841"}, + {"bess_nam", "a=6377483.865", "rf=299.1528128", "Bessel 1841 (Namibia)"}, + {"clrk66", "a=6378206.4", "b=6356583.8", "Clarke 1866"}, + {"clrk80", "a=6378249.145", "rf=293.4663", "Clarke 1880 mod."}, + {"clrk80ign", "a=6378249.2", "rf=293.4660212936269", "Clarke 1880 (IGN)."}, + {"CPM", "a=6375738.7", "rf=334.29", "Comm. des Poids et Mesures 1799"}, + {"delmbr", "a=6376428.", "rf=311.5", "Delambre 1810 (Belgium)"}, + {"engelis", "a=6378136.05", "rf=298.2566", "Engelis 1985"}, + {"evrst30", "a=6377276.345", "rf=300.8017", "Everest 1830"}, + {"evrst48", "a=6377304.063", "rf=300.8017", "Everest 1948"}, + {"evrst56", "a=6377301.243", "rf=300.8017", "Everest 1956"}, + {"evrst69", "a=6377295.664", "rf=300.8017", "Everest 1969"}, + {"evrstSS", "a=6377298.556", "rf=300.8017", "Everest (Sabah & Sarawak)"}, + {"fschr60", "a=6378166.", "rf=298.3", "Fischer (Mercury Datum) 1960"}, + {"fschr60m", "a=6378155.", "rf=298.3", "Modified Fischer 1960"}, + {"fschr68", "a=6378150.", "rf=298.3", "Fischer 1968"}, + {"helmert", "a=6378200.", "rf=298.3", "Helmert 1906"}, + {"hough", "a=6378270.0", "rf=297.", "Hough"}, + {"intl", "a=6378388.0", "rf=297.", "International 1909 (Hayford)"}, + {"krass", "a=6378245.0", "rf=298.3", "Krassovsky, 1942"}, + {"kaula", "a=6378163.", "rf=298.24", "Kaula 1961"}, + {"lerch", "a=6378139.", "rf=298.257", "Lerch 1979"}, + {"mprts", "a=6397300.", "rf=191.", "Maupertius 1738"}, + {"new_intl", "a=6378157.5", "b=6356772.2", "New International 1967"}, + {"plessis", "a=6376523.", "b=6355863.", "Plessis 1817 (France)"}, + {"SEasia", "a=6378155.0", "b=6356773.3205", "Southeast Asia"}, + {"walbeck", "a=6376896.0", "b=6355834.8467", "Walbeck"}, + {"WGS60", "a=6378165.0", "rf=298.3", "WGS 60"}, + {"WGS66", "a=6378145.0", "rf=298.25", "WGS 66"}, + {"WGS72", "a=6378135.0", "rf=298.26", "WGS 72"}, + {"WGS84", "a=6378137.0", "rf=298.257223563", "WGS 84"}, + {"sphere", "a=6370997.0", "b=6370997.0", "Normal Sphere (r=6370997)"} +}; + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELLPS_HPP diff --git a/boost/geometry/srs/projections/impl/pj_fwd.hpp b/boost/geometry/srs/projections/impl/pj_fwd.hpp new file mode 100644 index 0000000000..73101b7f40 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_fwd.hpp @@ -0,0 +1,101 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_FWD_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_FWD_HPP + +#include <boost/geometry/core/radian_access.hpp> +#include <boost/geometry/util/math.hpp> + +#include <boost/geometry/srs/projections/impl/adjlon.hpp> +#include <boost/geometry/srs/projections/impl/projects.hpp> + +#include <boost/math/constants/constants.hpp> + +/* general forward projection */ + +namespace boost { namespace geometry { namespace projections { + +namespace detail { + +/* forward projection entry */ +template <typename Prj, typename LL, typename XY, typename P> +inline void pj_fwd(Prj const& prj, P const& par, LL const& ll, XY& xy) +{ + typedef typename P::type calc_t; + static const calc_t EPS = 1.0e-12; + + using namespace detail; + + calc_t lp_lon = geometry::get_as_radian<0>(ll); + calc_t lp_lat = geometry::get_as_radian<1>(ll); + calc_t const t = geometry::math::abs(lp_lat) - geometry::math::half_pi<calc_t>(); + + /* check for forward and latitude or longitude overange */ + if (t > EPS || geometry::math::abs(lp_lon) > 10.) + { + BOOST_THROW_EXCEPTION( projection_exception(-14) ); + } + + if (geometry::math::abs(t) <= EPS) + { + lp_lat = lp_lat < 0. ? -geometry::math::half_pi<calc_t>() : geometry::math::half_pi<calc_t>(); + } + else if (par.geoc) + { + lp_lat = atan(par.rone_es * tan(lp_lat)); + } + + lp_lon -= par.lam0; /* compute del lp.lam */ + if (! par.over) + { + lp_lon = adjlon(lp_lon); /* post_forward del longitude */ + } + + calc_t x = 0; + calc_t y = 0; + + prj.fwd(lp_lon, lp_lat, x, y); + geometry::set<0>(xy, par.fr_meter * (par.a * x + par.x0)); + geometry::set<1>(xy, par.fr_meter * (par.a * y + par.y0)); +} + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_FWD_HPP diff --git a/boost/geometry/srs/projections/impl/pj_gauss.hpp b/boost/geometry/srs/projections/impl/pj_gauss.hpp new file mode 100644 index 0000000000..2c90870434 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_gauss.hpp @@ -0,0 +1,142 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_GAUSS_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_GAUSS_HPP + + +#include <boost/geometry/util/math.hpp> + + +namespace boost { namespace geometry { namespace projections { + +namespace detail { namespace gauss { + + +static const int MAX_ITER = 20; + +template <typename T> +struct GAUSS +{ + T C; + T K; + T e; + T ratexp; +}; + +template <typename T> +inline T srat(T const& esinp, T const& exp) +{ + return (pow((1.0 - esinp) / (1.0 + esinp), exp)); +} + +template <typename T> +inline GAUSS<T> gauss_ini(T const& e, T const& phi0, T& chi, T& rc) +{ + static const T FORTPI = detail::FORTPI<T>(); + + using std::asin; + using std::cos; + using std::sin; + using std::sqrt; + using std::tan; + + T sphi = 0; + T cphi = 0; + T es = 0; + + GAUSS<T> en; + es = e * e; + en.e = e; + sphi = sin(phi0); + cphi = cos(phi0); + cphi *= cphi; + + rc = sqrt(1.0 - es) / (1.0 - es * sphi * sphi); + en.C = sqrt(1.0 + es * cphi * cphi / (1.0 - es)); + chi = asin(sphi / en.C); + en.ratexp = 0.5 * en.C * e; + en.K = tan(0.5 * chi + FORTPI) + / (pow(tan(0.5 * phi0 + FORTPI), en.C) * srat(en.e * sphi, en.ratexp)); + + return en; +} + +template <typename T> +inline void gauss(GAUSS<T> const& en, T& lam, T& phi) +{ + static const T FORTPI = detail::FORTPI<T>(); + + phi = 2.0 * atan(en.K * pow(tan(0.5 * phi + FORTPI), en.C) + * srat(en.e * sin(phi), en.ratexp) ) - geometry::math::half_pi<T>(); + + lam *= en.C; +} + +template <typename T> +inline void inv_gauss(GAUSS<T> const& en, T& lam, T& phi) +{ + static const T FORTPI = detail::FORTPI<T>(); + static const T DEL_TOL = 1e-14; + + lam /= en.C; + const T num = pow(tan(0.5 * phi + FORTPI) / en.K, 1.0 / en.C); + + int i = 0; + for (i = MAX_ITER; i; --i) + { + const T elp_phi = 2.0 * atan(num * srat(en.e * sin(phi), - 0.5 * en.e)) - geometry::math::half_pi<T>(); + + if (geometry::math::abs(elp_phi - phi) < DEL_TOL) + { + break; + } + phi = elp_phi; + } + + /* convergence failed */ + if (!i) + { + BOOST_THROW_EXCEPTION( projection_exception(-17) ); + } +} + +}} // namespace detail::gauss +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_GAUSS_HPP diff --git a/boost/geometry/srs/projections/impl/pj_init.hpp b/boost/geometry/srs/projections/impl/pj_init.hpp new file mode 100644 index 0000000000..8b93343bcd --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_init.hpp @@ -0,0 +1,395 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017, 2018. +// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_INIT_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_INIT_HPP + +#include <cstdlib> +#include <string> +#include <vector> + +#include <boost/algorithm/string.hpp> +#include <boost/lexical_cast.hpp> +#include <boost/range.hpp> +#include <boost/type_traits/is_same.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/condition.hpp> + +#include <boost/geometry/srs/projections/impl/dms_parser.hpp> +#include <boost/geometry/srs/projections/impl/pj_datum_set.hpp> +#include <boost/geometry/srs/projections/impl/pj_datums.hpp> +#include <boost/geometry/srs/projections/impl/pj_ell_set.hpp> +#include <boost/geometry/srs/projections/impl/pj_param.hpp> +#include <boost/geometry/srs/projections/impl/pj_units.hpp> +#include <boost/geometry/srs/projections/impl/projects.hpp> +#include <boost/geometry/srs/projections/proj4.hpp> + + +namespace boost { namespace geometry { namespace projections +{ + + +namespace detail +{ + +template <typename BGParams, typename T> +inline void pj_push_defaults(BGParams const& /*bg_params*/, parameters<T>& pin) +{ + pin.params.push_back(pj_mkparam<T>("ellps=WGS84")); + + if (pin.name == "aea") + { + pin.params.push_back(pj_mkparam<T>("lat_1=29.5")); + pin.params.push_back(pj_mkparam<T>("lat_2=45.5 ")); + } + else if (pin.name == "lcc") + { + pin.params.push_back(pj_mkparam<T>("lat_1=33")); + pin.params.push_back(pj_mkparam<T>("lat_2=45")); + } + else if (pin.name == "lagrng") + { + pin.params.push_back(pj_mkparam<T>("W=2")); + } +} + +template <BOOST_GEOMETRY_PROJECTIONS_DETAIL_TYPENAME_PX, typename T> +inline void pj_push_defaults(srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAIL_PX> const& /*bg_params*/, + parameters<T>& pin) +{ + typedef srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAIL_PX> static_parameters_type; + typedef typename srs::par4::detail::pick_proj_tag + < + static_parameters_type + >::type proj_tag; + + // statically defaulting to WGS84 + //pin.params.push_back(pj_mkparam("ellps=WGS84")); + + if (BOOST_GEOMETRY_CONDITION((boost::is_same<proj_tag, srs::par4::aea>::value))) + { + pin.params.push_back(pj_mkparam<T>("lat_1=29.5")); + pin.params.push_back(pj_mkparam<T>("lat_2=45.5 ")); + } + else if (BOOST_GEOMETRY_CONDITION((boost::is_same<proj_tag, srs::par4::lcc>::value))) + { + pin.params.push_back(pj_mkparam<T>("lat_1=33")); + pin.params.push_back(pj_mkparam<T>("lat_2=45")); + } + else if (BOOST_GEOMETRY_CONDITION((boost::is_same<proj_tag, srs::par4::lagrng>::value))) + { + pin.params.push_back(pj_mkparam<T>("W=2")); + } +} + +template <typename T> +inline void pj_init_units(std::vector<pvalue<T> > const& params, + std::string const& sunits, + std::string const& sto_meter, + T & to_meter, + T & fr_meter, + T const& default_to_meter, + T const& default_fr_meter) +{ + std::string s; + std::string units = pj_param(params, sunits).s; + if (! units.empty()) + { + const int n = sizeof(pj_units) / sizeof(pj_units[0]); + int index = -1; + for (int i = 0; i < n && index == -1; i++) + { + if(pj_units[i].id == units) + { + index = i; + } + } + + if (index == -1) { + BOOST_THROW_EXCEPTION( projection_exception(-7) ); + } + s = pj_units[index].to_meter; + } + + if (s.empty()) + { + s = pj_param(params, sto_meter).s; + } + + if (! s.empty()) + { + std::size_t const pos = s.find('/'); + if (pos == std::string::npos) + { + to_meter = lexical_cast<T>(s); + } + else + { + T const numerator = lexical_cast<T>(s.substr(0, pos)); + T const denominator = lexical_cast<T>(s.substr(pos + 1)); + if (numerator == 0.0 || denominator == 0.0) + { + BOOST_THROW_EXCEPTION( projection_exception(-99) ); + } + to_meter = numerator / denominator; + } + if (to_meter == 0.0) + { + BOOST_THROW_EXCEPTION( projection_exception(-99) ); + } + fr_meter = 1. / to_meter; + } + else + { + to_meter = default_to_meter; + fr_meter = default_fr_meter; + } +} + +/************************************************************************/ +/* pj_init() */ +/* */ +/* Main entry point for initialing a PJ projections */ +/* definition. Note that the projection specific function is */ +/* called to do the initial allocation so it can be created */ +/* large enough to hold projection specific parameters. */ +/************************************************************************/ +template <typename T, typename BGParams, typename R> +inline parameters<T> pj_init(BGParams const& bg_params, R const& arguments, bool use_defaults = true) +{ + parameters<T> pin; + for (std::vector<std::string>::const_iterator it = boost::begin(arguments); + it != boost::end(arguments); it++) + { + pin.params.push_back(pj_mkparam<T>(*it)); + } + + /* check if +init present */ + if (pj_param(pin.params, "tinit").i) + { + // maybe TODO: handle "init" parameter + //if (!(curr = get_init(&arguments, curr, pj_param(pin.params, "sinit").s))) + } + + // find projection -> implemented in projection factory + pin.name = pj_param(pin.params, "sproj").s; + // exception thrown in projection<> + // TODO: consider throwing here both projection_unknown_id_exception and + // projection_not_named_exception in order to throw before other exceptions + //if (pin.name.empty()) + //{ BOOST_THROW_EXCEPTION( projection_not_named_exception() ); } + + // set defaults, unless inhibited + // GL-Addition, if use_defaults is false then defaults are ignored + if (use_defaults && ! pj_param(pin.params, "bno_defs").i) + { + // proj4 gets defaults from "proj_def.dat", file of 94/02/23 with a few defaults. + // Here manually + pj_push_defaults(bg_params, pin); + //curr = get_defaults(&arguments, curr, name); + } + + /* allocate projection structure */ + // done by BGParams constructor: + // pin.is_latlong = 0; + // pin.is_geocent = 0; + // pin.long_wrap_center = 0.0; + // pin.long_wrap_center = 0.0; + pin.is_long_wrap_set = false; + + /* set datum parameters */ + pj_datum_set(bg_params, pin.params, pin); + + /* set ellipsoid/sphere parameters */ + pj_ell_set(bg_params, pin.params, pin.a, pin.es); + + pin.a_orig = pin.a; + pin.es_orig = pin.es; + + pin.e = sqrt(pin.es); + pin.ra = 1. / pin.a; + pin.one_es = 1. - pin.es; + if (pin.one_es == 0.) { + BOOST_THROW_EXCEPTION( projection_exception(-6) ); + } + pin.rone_es = 1./pin.one_es; + + /* Now that we have ellipse information check for WGS84 datum */ + if( pin.datum_type == PJD_3PARAM + && pin.datum_params[0] == 0.0 + && pin.datum_params[1] == 0.0 + && pin.datum_params[2] == 0.0 + && pin.a == 6378137.0 + && geometry::math::abs(pin.es - 0.006694379990) < 0.000000000050 )/*WGS84/GRS80*/ + { + pin.datum_type = PJD_WGS84; + } + + /* set pin.geoc coordinate system */ + pin.geoc = (pin.es && pj_param(pin.params, "bgeoc").i); + + /* over-ranging flag */ + pin.over = pj_param(pin.params, "bover").i; + + /* longitude center for wrapping */ + pin.is_long_wrap_set = pj_param(pin.params, "tlon_wrap").i != 0; + if (pin.is_long_wrap_set) + pin.long_wrap_center = pj_param(pin.params, "rlon_wrap").f; + + /* central meridian */ + pin.lam0 = pj_param(pin.params, "rlon_0").f; + + /* central latitude */ + pin.phi0 = pj_param(pin.params, "rlat_0").f; + + /* false easting and northing */ + pin.x0 = pj_param(pin.params, "dx_0").f; + pin.y0 = pj_param(pin.params, "dy_0").f; + + /* general scaling factor */ + if (pj_param(pin.params, "tk_0").i) + pin.k0 = pj_param(pin.params, "dk_0").f; + else if (pj_param(pin.params, "tk").i) + pin.k0 = pj_param(pin.params, "dk").f; + else + pin.k0 = 1.; + if (pin.k0 <= 0.) { + BOOST_THROW_EXCEPTION( projection_exception(-31) ); + } + + /* set units */ + pj_init_units(pin.params, "sunits", "sto_meter", + pin.to_meter, pin.fr_meter, 1., 1.); + pj_init_units(pin.params, "svunits", "svto_meter", + pin.vto_meter, pin.vfr_meter, pin.to_meter, pin.fr_meter); + + /* prime meridian */ + std::string pm = pj_param(pin.params, "spm").s; + if (! pm.empty()) + { + std::string value; + + int n = sizeof(pj_prime_meridians) / sizeof(pj_prime_meridians[0]); + for (int i = 0; i < n ; i++) + { + if(pj_prime_meridians[i].id == pm) + { + value = pj_prime_meridians[i].defn; + break; + } + } + + dms_parser<T, true> parser; + + // TODO: Handle case when lexical_cast is not used consistently. + // This should probably be done in dms_parser. + BOOST_TRY + { + if (value.empty()) { + pin.from_greenwich = parser.apply(pm).angle(); + } else { + pin.from_greenwich = parser.apply(value).angle(); + } + } + BOOST_CATCH(boost::bad_lexical_cast const&) + { + BOOST_THROW_EXCEPTION( projection_exception(-46) ); + } + BOOST_CATCH_END + } + else + { + pin.from_greenwich = 0.0; + } + + return pin; +} + +/************************************************************************/ +/* pj_init_plus() */ +/* */ +/* Same as pj_init() except it takes one argument string with */ +/* individual arguments preceeded by '+', such as "+proj=utm */ +/* +zone=11 +ellps=WGS84". */ +/************************************************************************/ +template <typename T, typename BGParams> +inline parameters<T> pj_init_plus(BGParams const& bg_params, std::string const& definition, bool use_defaults = true) +{ + const char* sep = " +"; + + /* split into arguments based on '+' and trim white space */ + + // boost::split splits on one character, here it should be on " +", so implementation below + // todo: put in different routine or sort out + std::vector<std::string> arguments; + std::string def = boost::trim_copy(definition); + boost::trim_left_if(def, boost::is_any_of(sep)); + + std::string::size_type loc = def.find(sep); + while (loc != std::string::npos) + { + std::string par = def.substr(0, loc); + boost::trim(par); + if (! par.empty()) + { + arguments.push_back(par); + } + + def.erase(0, loc); + boost::trim_left_if(def, boost::is_any_of(sep)); + loc = def.find(sep); + } + + if (! def.empty()) + { + arguments.push_back(def); + } + + /*boost::split(arguments, definition, boost::is_any_of("+")); + for (std::vector<std::string>::iterator it = arguments.begin(); it != arguments.end(); it++) + { + boost::trim(*it); + }*/ + return pj_init<T>(bg_params, arguments, use_defaults); +} + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_INIT_HPP diff --git a/boost/geometry/srs/projections/impl/pj_inv.hpp b/boost/geometry/srs/projections/impl/pj_inv.hpp new file mode 100644 index 0000000000..91cf4c50f6 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_inv.hpp @@ -0,0 +1,82 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_PJ_INV_HPP +#define BOOST_GEOMETRY_PROJECTIONS_PJ_INV_HPP + + + +#include <boost/geometry/srs/projections/impl/adjlon.hpp> +#include <boost/geometry/core/radian_access.hpp> +#include <boost/geometry/util/math.hpp> + +/* general inverse projection */ + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + + /* inverse projection entry */ +template <typename PRJ, typename LL, typename XY, typename PAR> +inline void pj_inv(PRJ const& prj, PAR const& par, XY const& xy, LL& ll) +{ + typedef typename PAR::type calc_t; + static const calc_t EPS = 1.0e-12; + + /* can't do as much preliminary checking as with forward */ + /* descale and de-offset */ + calc_t xy_x = (geometry::get<0>(xy) * par.to_meter - par.x0) * par.ra; + calc_t xy_y = (geometry::get<1>(xy) * par.to_meter - par.y0) * par.ra; + calc_t lon = 0, lat = 0; + prj.inv(xy_x, xy_y, lon, lat); /* inverse project */ + lon += par.lam0; /* reduce from del lp.lam */ + if (!par.over) + lon = adjlon(lon); /* adjust longitude to CM */ + if (par.geoc && geometry::math::abs(geometry::math::abs(lat)-geometry::math::half_pi<calc_t>()) > EPS) + lat = atan(par.one_es * tan(lat)); + + geometry::set_from_radian<0>(ll, lon); + geometry::set_from_radian<1>(ll, lat); +} + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif diff --git a/boost/geometry/srs/projections/impl/pj_mlfn.hpp b/boost/geometry/srs/projections/impl/pj_mlfn.hpp new file mode 100644 index 0000000000..1317dd2e2f --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_mlfn.hpp @@ -0,0 +1,121 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_PJ_MLFN_HPP +#define BOOST_GEOMETRY_PROJECTIONS_PJ_MLFN_HPP + + + +#include <boost/geometry/util/math.hpp> + + +namespace boost { namespace geometry { namespace projections { + +namespace detail { + +static const int EN_SIZE = 5; + +template <typename T> +inline bool pj_enfn(T const& es, T* en) +{ + static const T C00 = 1.; + static const T C02 = .25; + static const T C04 = .046875; + static const T C06 = .01953125; + static const T C08 = .01068115234375; + static const T C22 = .75; + static const T C44 = .46875; + static const T C46 = .01302083333333333333; + static const T C48 = .00712076822916666666; + static const T C66 = .36458333333333333333; + static const T C68 = .00569661458333333333; + static const T C88 = .3076171875; + + T t; //, *en; + + //if (en = (double *)pj_malloc(EN_SIZE * sizeof(double))) + { + en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08))); + en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08))); + en[2] = (t = es * es) * (C44 - es * (C46 + es * C48)); + en[3] = (t *= es) * (C66 - es * C68); + en[4] = t * es * C88; + } + // return en; + return true; +} + +template <typename T> +inline T pj_mlfn(T const& phi, T sphi, T cphi, const T *en) +{ + cphi *= sphi; + sphi *= sphi; + return(en[0] * phi - cphi * (en[1] + sphi*(en[2] + + sphi*(en[3] + sphi*en[4])))); +} + +template <typename T> +inline T pj_inv_mlfn(T const& arg, T const& es, const T *en) +{ + /* meridinal distance for ellipsoid and inverse + ** 8th degree - accurate to < 1e-5 meters when used in conjuction + ** with typical major axis values. + ** Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds. + */ + static const T EPS = 1e-11; + static const int MAX_ITER = 10; + + T s, t, phi, k = 1./(1.-es); + int i; + + phi = arg; + for (i = MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */ + s = sin(phi); + t = 1. - es * s * s; + phi -= t = (pj_mlfn(phi, s, cos(phi), en) - arg) * (t * sqrt(t)) * k; + if (geometry::math::abs(t) < EPS) + return phi; + } + BOOST_THROW_EXCEPTION( projection_exception(-17) ); + return phi; +} + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif diff --git a/boost/geometry/srs/projections/impl/pj_msfn.hpp b/boost/geometry/srs/projections/impl/pj_msfn.hpp new file mode 100644 index 0000000000..f99cc2c843 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_msfn.hpp @@ -0,0 +1,59 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_PJ_MSFN_HPP +#define BOOST_GEOMETRY_PROJECTIONS_PJ_MSFN_HPP + + +namespace boost { namespace geometry { namespace projections { + +namespace detail { + + +/* determine constant small m */ +template <typename T> +inline T pj_msfn(T const& sinphi, T const& cosphi, T const& es) +{ + return (cosphi / sqrt (1. - es * sinphi * sinphi)); +} + + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif diff --git a/boost/geometry/srs/projections/impl/pj_param.hpp b/boost/geometry/srs/projections/impl/pj_param.hpp new file mode 100644 index 0000000000..4f33ad837f --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_param.hpp @@ -0,0 +1,161 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_PJ_PARAM_HPP +#define BOOST_GEOMETRY_PROJECTIONS_PJ_PARAM_HPP + + +#include <string> +#include <vector> + +#include <boost/geometry/srs/projections/impl/dms_parser.hpp> +#include <boost/geometry/srs/projections/impl/projects.hpp> + + +namespace boost { namespace geometry { namespace projections { + +namespace detail { + + + +/* create pvalue list entry */ +template <typename T> +inline pvalue<T> pj_mkparam(std::string const& str) +{ + std::string name = str; + std::string value; + boost::trim_left_if(name, boost::is_any_of("+")); + std::string::size_type loc = name.find("="); + if (loc != std::string::npos) + { + value = name.substr(loc + 1); + name.erase(loc); + } + + + pvalue<T> newitem; + newitem.param = name; + newitem.s = value; + newitem.used = 0; + newitem.i = atoi(value.c_str()); + newitem.f = atof(value.c_str()); + return newitem; +} + +/************************************************************************/ +/* pj_param() */ +/* */ +/* Test for presence or get pvalue value. The first */ +/* character in `opt' is a pvalue type which can take the */ +/* values: */ +/* */ +/* `t' - test for presence, return TRUE/FALSE in pvalue.i */ +/* `i' - integer value returned in pvalue.i */ +/* `d' - simple valued real input returned in pvalue.f */ +/* `r' - degrees (DMS translation applied), returned as */ +/* radians in pvalue.f */ +/* `s' - string returned in pvalue.s */ +/* `b' - test for t/T/f/F, return in pvalue.i */ +/* */ +/************************************************************************/ + +template <typename T> +inline pvalue<T> pj_param(std::vector<pvalue<T> > const& pl, std::string opt) +{ + char type = opt[0]; + opt.erase(opt.begin()); + + pvalue<T> value; + + /* simple linear lookup */ + typedef typename std::vector<pvalue<T> >::const_iterator iterator; + for (iterator it = pl.begin(); it != pl.end(); it++) + { + if (it->param == opt) + { + //it->used = 1; + switch (type) + { + case 't': + value.i = 1; + break; + case 'i': /* integer input */ + value.i = atoi(it->s.c_str()); + break; + case 'd': /* simple real input */ + value.f = atof(it->s.c_str()); + break; + case 'r': /* degrees input */ + { + dms_parser<T, true> parser; + value.f = parser.apply(it->s.c_str()).angle(); + } + break; + case 's': /* char string */ + value.s = it->s; + break; + case 'b': /* boolean */ + switch (it->s[0]) + { + case 'F': case 'f': + value.i = 0; + break; + case '\0': case 'T': case 't': + value.i = 1; + break; + default: + value.i = 0; + break; + } + break; + } + return value; + } + + } + + value.i = 0; + value.f = 0.0; + value.s = ""; + return value; +} + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif diff --git a/boost/geometry/srs/projections/impl/pj_phi2.hpp b/boost/geometry/srs/projections/impl/pj_phi2.hpp new file mode 100644 index 0000000000..71f0cf1249 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_phi2.hpp @@ -0,0 +1,73 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_PHI2_HPP +#define BOOST_GEOMETRY_PROJECTIONS_PHI2_HPP + +#include <boost/geometry/util/math.hpp> + +namespace boost { namespace geometry { namespace projections { +namespace detail { + +template <typename T> +inline T pj_phi2(T const& ts, T const& e) +{ + static const T TOL = 1.0e-10; + static const int N_ITER = 15; + + T eccnth, Phi, con, dphi; + int i; + + eccnth = .5 * e; + Phi = geometry::math::half_pi<T>() - 2. * atan (ts); + i = N_ITER; + do { + con = e * sin (Phi); + dphi = geometry::math::half_pi<T>() - 2. * atan (ts * pow((1. - con) / + (1. + con), eccnth)) - Phi; + Phi += dphi; + } while ( geometry::math::abs(dphi) > TOL && --i); + if (i <= 0) + BOOST_THROW_EXCEPTION( projection_exception(-18) ); + return Phi; +} + +} // namespace detail +}}} // namespace boost::geometry::projections + +#endif diff --git a/boost/geometry/srs/projections/impl/pj_qsfn.hpp b/boost/geometry/srs/projections/impl/pj_qsfn.hpp new file mode 100644 index 0000000000..606fab5e4e --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_qsfn.hpp @@ -0,0 +1,95 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_PJ_QSFN_HPP +#define BOOST_GEOMETRY_PROJECTIONS_PJ_QSFN_HPP + + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail { + +/* determine small q */ +template <typename T> +inline T pj_qsfn(T const& sinphi, T const& e, T const& one_es) +{ + static const T EPSILON = 1.0e-7; + + if (e >= EPSILON) + { + T con = e * sinphi; + return (one_es * (sinphi / (1. - con * con) - + (.5 / e) * log ((1. - con) / (1. + con)))); + } else + return (sinphi + sinphi); +} + + +static const int MAX_C = 9; + +template <typename T> +struct AUTHALIC +{ + T C[MAX_C], CP[MAX_C], CQ[MAX_C]; +}; + +/** + * @brief determine authalic latitude + * @param[in] phi geographic latitude + * @param[in] a initialized structure pointer + * @return authalic latitude + */ +template <typename T> +inline T proj_qsfn(T const& phi, AUTHALIC<T> const& a) +{ + T s, s2, sum; + int i = MAX_C; + + s = sin(phi); + s2 = s * s; + sum = a.CQ[MAX_C - 1]; + while (--i) sum = a.CQ[i] + s2 * sum; + return(s * sum); +} + +} // namespace detail + +}}} // namespace boost::geometry::projections + +#endif diff --git a/boost/geometry/srs/projections/impl/pj_strerrno.hpp b/boost/geometry/srs/projections/impl/pj_strerrno.hpp new file mode 100644 index 0000000000..22e0c48af4 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_strerrno.hpp @@ -0,0 +1,144 @@ +// Boost.Geometry +// This file is manually converted from PROJ4 + +// This file was modified by Oracle on 2017, 2018. +// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 +// This file was converted to Geometry Library by Adam Wulkiewicz + +// Original copyright notice: + +// None + +/* list of projection system pj_errno values */ + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_STRERRNO_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_STRERRNO_HPP + +#include <cerrno> +#include <cstring> +#include <sstream> +#include <string> + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + + static const char * +pj_err_list[] = { + "no arguments in initialization list", /* -1 */ + "no options found in 'init' file", /* -2 */ + "no colon in init= string", /* -3 */ + "projection not named", /* -4 */ + "unknown projection id", /* -5 */ + "effective eccentricity = 1.", /* -6 */ + "unknown unit conversion id", /* -7 */ + "invalid boolean param argument", /* -8 */ + "unknown elliptical parameter name", /* -9 */ + "reciprocal flattening (1/f) = 0", /* -10 */ + "|radius reference latitude| > 90", /* -11 */ + "squared eccentricity < 0", /* -12 */ + "major axis or radius = 0 or not given", /* -13 */ + "latitude or longitude exceeded limits", /* -14 */ + "invalid x or y", /* -15 */ + "improperly formed DMS value", /* -16 */ + "non-convergent inverse meridional dist", /* -17 */ + "non-convergent inverse phi2", /* -18 */ + "acos/asin: |arg| >1.+1e-14", /* -19 */ + "tolerance condition error", /* -20 */ + "conic lat_1 = -lat_2", /* -21 */ + "lat_1 >= 90", /* -22 */ + "lat_1 = 0", /* -23 */ + "lat_ts >= 90", /* -24 */ + "no distance between control points", /* -25 */ + "projection not selected to be rotated", /* -26 */ + "W <= 0 or M <= 0", /* -27 */ + "lsat not in 1-5 range", /* -28 */ + "path not in range", /* -29 */ + "h <= 0", /* -30 */ + "k <= 0", /* -31 */ + "lat_0 = 0 or 90 or alpha = 90", /* -32 */ + "lat_1=lat_2 or lat_1=0 or lat_2=90", /* -33 */ + "elliptical usage required", /* -34 */ + "invalid UTM zone number", /* -35 */ + "arg(s) out of range for Tcheby eval", /* -36 */ + "failed to find projection to be rotated", /* -37 */ + "failed to load datum shift file", /* -38 */ + "both n & m must be spec'd and > 0", /* -39 */ + "n <= 0, n > 1 or not specified", /* -40 */ + "lat_1 or lat_2 not specified", /* -41 */ + "|lat_1| == |lat_2|", /* -42 */ + "lat_0 is pi/2 from mean lat", /* -43 */ + "unparseable coordinate system definition", /* -44 */ + "geocentric transformation missing z or ellps", /* -45 */ + "unknown prime meridian conversion id", /* -46 */ + "illegal axis orientation combination", /* -47 */ + "point not within available datum shift grids", /* -48 */ + "invalid sweep axis, choose x or y", /* -49 */ + "malformed pipeline", /* -50 */ +}; + +inline std::string pj_generic_strerrno(std::string const& msg, int err) +{ + std::stringstream ss; + ss << msg << " (" << err << ")"; + return ss.str(); +} + +inline std::string pj_strerrno(int err) { + if (0==err) + { + return ""; + } + else if (err > 0) + { + // std::strerror function may be not thread-safe + //return std::strerror(err); + + switch(err) + { +#ifdef EINVAL + case EINVAL: + return "Invalid argument"; +#endif +#ifdef EDOM + case EDOM: + return "Math argument out of domain of func"; +#endif +#ifdef ERANGE + case ERANGE: + return "Math result not representable"; +#endif + default: + return pj_generic_strerrno("system error", err); + } + } + else /*if (err < 0)*/ + { + size_t adjusted_err = - err - 1; + if (adjusted_err < (sizeof(pj_err_list) / sizeof(char *))) + { + return(pj_err_list[adjusted_err]); + } + else + { + return pj_generic_strerrno("invalid projection system error", err); + } + } +} + +} // namespace detail + +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_STRERRNO_HPP diff --git a/boost/geometry/srs/projections/impl/pj_transform.hpp b/boost/geometry/srs/projections/impl/pj_transform.hpp new file mode 100644 index 0000000000..8c2095642f --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_transform.hpp @@ -0,0 +1,1009 @@ +// Boost.Geometry +// This file is manually converted from PROJ4 + +// This file was modified by Oracle on 2017, 2018. +// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 +// This file was converted to Geometry Library by Adam Wulkiewicz + +// Original copyright notice: + +// Copyright (c) 2000, Frank Warmerdam + +// 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. + +#ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP +#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP + + +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/radian_access.hpp> + +#include <boost/geometry/srs/projections/impl/geocent.hpp> +#include <boost/geometry/srs/projections/impl/projects.hpp> +#include <boost/geometry/srs/projections/invalid_point.hpp> + +#include <boost/geometry/util/range.hpp> + +#include <cstring> +#include <cmath> + + +namespace boost { namespace geometry { namespace projections +{ + +namespace detail +{ + +// ----------------------------------------------------------- +// Boost.Geometry helpers begin +// ----------------------------------------------------------- + +template +< + typename Point, + bool HasCoord2 = (dimension<Point>::value > 2) +> +struct z_access +{ + typedef typename coordinate_type<Point>::type type; + static inline type get(Point const& point) + { + return geometry::get<2>(point); + } + static inline void set(Point & point, type const& h) + { + return geometry::set<2>(point, h); + } +}; + +template <typename Point> +struct z_access<Point, false> +{ + typedef typename coordinate_type<Point>::type type; + static inline type get(Point const& ) + { + return type(0); + } + static inline void set(Point & , type const& ) + {} +}; + +template <typename XYorXYZ> +inline typename z_access<XYorXYZ>::type +get_z(XYorXYZ const& xy_or_xyz) +{ + return z_access<XYorXYZ>::get(xy_or_xyz); +} + +template <typename XYorXYZ> +inline void set_z(XYorXYZ & xy_or_xyz, + typename z_access<XYorXYZ>::type const& z) +{ + return z_access<XYorXYZ>::set(xy_or_xyz, z); +} + +template +< + typename Range, + bool AddZ = (dimension<typename boost::range_value<Range>::type>::value < 3) +> +struct range_wrapper +{ + typedef Range range_type; + typedef typename boost::range_value<Range>::type point_type; + typedef typename coordinate_type<point_type>::type coord_t; + + range_wrapper(Range & range) + : m_range(range) + {} + + range_type & get_range() { return m_range; } + + coord_t get_z(std::size_t i) { return detail::get_z(range::at(m_range, i)); } + void set_z(std::size_t i, coord_t const& v) { return detail::set_z(range::at(m_range, i), v); } + +private: + Range & m_range; +}; + +template <typename Range> +struct range_wrapper<Range, true> +{ + typedef Range range_type; + typedef typename boost::range_value<Range>::type point_type; + typedef typename coordinate_type<point_type>::type coord_t; + + range_wrapper(Range & range) + : m_range(range) + , m_zs(boost::size(range), coord_t(0)) + {} + + range_type & get_range() { return m_range; } + + coord_t get_z(std::size_t i) { return m_zs[i]; } + void set_z(std::size_t i, coord_t const& v) { m_zs[i] = v; } + +private: + Range & m_range; + std::vector<coord_t> m_zs; +}; + +// ----------------------------------------------------------- +// Boost.Geometry helpers end +// ----------------------------------------------------------- + +/*#ifndef SRS_WGS84_SEMIMAJOR +#define SRS_WGS84_SEMIMAJOR 6378137.0 +#endif + +#ifndef SRS_WGS84_ESQUARED +#define SRS_WGS84_ESQUARED 0.0066943799901413165 +#endif*/ + +template <typename Par> +inline typename Par::type Dx_BF(Par const& defn) { return defn.datum_params[0]; } +template <typename Par> +inline typename Par::type Dy_BF(Par const& defn) { return defn.datum_params[1]; } +template <typename Par> +inline typename Par::type Dz_BF(Par const& defn) { return defn.datum_params[2]; } +template <typename Par> +inline typename Par::type Rx_BF(Par const& defn) { return defn.datum_params[3]; } +template <typename Par> +inline typename Par::type Ry_BF(Par const& defn) { return defn.datum_params[4]; } +template <typename Par> +inline typename Par::type Rz_BF(Par const& defn) { return defn.datum_params[5]; } +template <typename Par> +inline typename Par::type M_BF(Par const& defn) { return defn.datum_params[6]; } + +/* +** This table is intended to indicate for any given error code in +** the range 0 to -44, whether that error will occur for all locations (ie. +** it is a problem with the coordinate system as a whole) in which case the +** value would be 0, or if the problem is with the point being transformed +** in which case the value is 1. +** +** At some point we might want to move this array in with the error message +** list or something, but while experimenting with it this should be fine. +*/ + +static const int transient_error[50] = { + /* 0 1 2 3 4 5 6 7 8 9 */ + /* 0 to 9 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, + /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, + /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }; + + +template <typename T, typename Range> +inline int pj_geocentric_to_geodetic( T const& a, T const& es, + Range & range ); +template <typename T, typename Range> +inline int pj_geodetic_to_geocentric( T const& a, T const& es, + Range & range ); + +/************************************************************************/ +/* pj_transform() */ +/* */ +/* Currently this function doesn't recognise if two projections */ +/* are identical (to short circuit reprojection) because it is */ +/* difficult to compare PJ structures (since there are some */ +/* projection specific components). */ +/************************************************************************/ + +template <typename SrcPrj, typename DstPrj2, typename Par, typename Range> +inline bool pj_transform(SrcPrj const& srcprj, Par const& srcdefn, + DstPrj2 const& dstprj, Par const& dstdefn, + Range & range) + +{ + typedef typename boost::range_value<Range>::type point_type; + typedef typename coordinate_type<point_type>::type coord_t; + static const std::size_t dimension = geometry::dimension<point_type>::value; + std::size_t point_count = boost::size(range); + bool result = true; + +/* -------------------------------------------------------------------- */ +/* Transform unusual input coordinate axis orientation to */ +/* standard form if needed. */ +/* -------------------------------------------------------------------- */ + // Ignored + +/* -------------------------------------------------------------------- */ +/* Transform Z to meters if it isn't already. */ +/* -------------------------------------------------------------------- */ + if( srcdefn.vto_meter != 1.0 && dimension > 2 ) + { + for( std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = geometry::range::at(range, i); + set_z(point, get_z(point) * srcdefn.vto_meter); + } + } + +/* -------------------------------------------------------------------- */ +/* Transform geocentric source coordinates to lat/long. */ +/* -------------------------------------------------------------------- */ + if( srcdefn.is_geocent ) + { + // Point should be cartesian 3D (ECEF) + if (dimension < 3) + BOOST_THROW_EXCEPTION( projection_exception(PJD_ERR_GEOCENTRIC) ); + //return PJD_ERR_GEOCENTRIC; + + if( srcdefn.to_meter != 1.0 ) + { + for(std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(range, i); + if( ! is_invalid_point(point) ) + { + set<0>(point, get<0>(point) * srcdefn.to_meter); + set<1>(point, get<1>(point) * srcdefn.to_meter); + } + } + } + + range_wrapper<Range, false> rng(range); + int err = pj_geocentric_to_geodetic( srcdefn.a_orig, srcdefn.es_orig, + rng ); + if( err != 0 ) + BOOST_THROW_EXCEPTION( projection_exception(err) ); + //return err; + + // NOTE: here 3D cartesian ECEF is converted into 3D geodetic LLH + } + +/* -------------------------------------------------------------------- */ +/* Transform source points to lat/long, if they aren't */ +/* already. */ +/* -------------------------------------------------------------------- */ + else if( !srcdefn.is_latlong ) + { + // Point should be cartesian 2D or 3D (map projection) + + /* Check first if projection is invertible. */ + /*if( (srcdefn->inv3d == NULL) && (srcdefn->inv == NULL)) + { + pj_ctx_set_errno( pj_get_ctx(srcdefn), -17 ); + pj_log( pj_get_ctx(srcdefn), PJ_LOG_ERROR, + "pj_transform(): source projection not invertible" ); + return -17; + }*/ + + /* If invertible - First try inv3d if defined */ + //if (srcdefn->inv3d != NULL) + //{ + // /* Three dimensions must be defined */ + // if ( z == NULL) + // { + // pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC); + // return PJD_ERR_GEOCENTRIC; + // } + + // for (i=0; i < point_count; i++) + // { + // XYZ projected_loc; + // XYZ geodetic_loc; + + // projected_loc.u = x[point_offset*i]; + // projected_loc.v = y[point_offset*i]; + // projected_loc.w = z[point_offset*i]; + + // if (projected_loc.u == HUGE_VAL) + // continue; + + // geodetic_loc = pj_inv3d(projected_loc, srcdefn); + // if( srcdefn->ctx->last_errno != 0 ) + // { + // if( (srcdefn->ctx->last_errno != 33 /*EDOM*/ + // && srcdefn->ctx->last_errno != 34 /*ERANGE*/ ) + // && (srcdefn->ctx->last_errno > 0 + // || srcdefn->ctx->last_errno < -44 || point_count == 1 + // || transient_error[-srcdefn->ctx->last_errno] == 0 ) ) + // return srcdefn->ctx->last_errno; + // else + // { + // geodetic_loc.u = HUGE_VAL; + // geodetic_loc.v = HUGE_VAL; + // geodetic_loc.w = HUGE_VAL; + // } + // } + + // x[point_offset*i] = geodetic_loc.u; + // y[point_offset*i] = geodetic_loc.v; + // z[point_offset*i] = geodetic_loc.w; + + // } + + //} + //else + { + /* Fallback to the original PROJ.4 API 2d inversion - inv */ + for( std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(range, i); + + if( is_invalid_point(point) ) + continue; + + try + { + pj_inv(srcprj, srcdefn, point, point); + } + catch(projection_exception const& e) + { + if( (e.code() != 33 /*EDOM*/ + && e.code() != 34 /*ERANGE*/ ) + && (e.code() > 0 + || e.code() < -44 /*|| point_count == 1*/ + || transient_error[-e.code()] == 0) ) { + BOOST_RETHROW + } else { + set_invalid_point(point); + result = false; + if (point_count == 1) + return result; + } + } + } + } + } + +/* -------------------------------------------------------------------- */ +/* But if they are already lat long, adjust for the prime */ +/* meridian if there is one in effect. */ +/* -------------------------------------------------------------------- */ + if( srcdefn.from_greenwich != 0.0 ) + { + for( std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(range, i); + + if( ! is_invalid_point(point) ) + set<0>(point, get<0>(point) + srcdefn.from_greenwich); + } + } + +/* -------------------------------------------------------------------- */ +/* Do we need to translate from geoid to ellipsoidal vertical */ +/* datum? */ +/* -------------------------------------------------------------------- */ + /*if( srcdefn->has_geoid_vgrids && z != NULL ) + { + if( pj_apply_vgridshift( srcdefn, "sgeoidgrids", + &(srcdefn->vgridlist_geoid), + &(srcdefn->vgridlist_geoid_count), + 0, point_count, point_offset, x, y, z ) != 0 ) + return pj_ctx_get_errno(srcdefn->ctx); + }*/ + +/* -------------------------------------------------------------------- */ +/* Convert datums if needed, and possible. */ +/* -------------------------------------------------------------------- */ + if ( ! pj_datum_transform( srcdefn, dstdefn, range ) ) + { + result = false; + } + +/* -------------------------------------------------------------------- */ +/* Do we need to translate from ellipsoidal to geoid vertical */ +/* datum? */ +/* -------------------------------------------------------------------- */ + /*if( dstdefn->has_geoid_vgrids && z != NULL ) + { + if( pj_apply_vgridshift( dstdefn, "sgeoidgrids", + &(dstdefn->vgridlist_geoid), + &(dstdefn->vgridlist_geoid_count), + 1, point_count, point_offset, x, y, z ) != 0 ) + return dstdefn->ctx->last_errno; + }*/ + +/* -------------------------------------------------------------------- */ +/* But if they are staying lat long, adjust for the prime */ +/* meridian if there is one in effect. */ +/* -------------------------------------------------------------------- */ + if( dstdefn.from_greenwich != 0.0 ) + { + for( std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(range, i); + + if( ! is_invalid_point(point) ) + set<0>(point, get<0>(point) - dstdefn.from_greenwich); + } + } + +/* -------------------------------------------------------------------- */ +/* Transform destination latlong to geocentric if required. */ +/* -------------------------------------------------------------------- */ + if( dstdefn.is_geocent ) + { + // Point should be cartesian 3D (ECEF) + if (dimension < 3) + BOOST_THROW_EXCEPTION( projection_exception(PJD_ERR_GEOCENTRIC) ); + //return PJD_ERR_GEOCENTRIC; + + // NOTE: In the original code the return value of the following + // function is not checked + range_wrapper<Range, false> rng(range); + int err = pj_geodetic_to_geocentric( dstdefn.a_orig, dstdefn.es_orig, + rng ); + if( err == -14 ) + result = false; + else + BOOST_THROW_EXCEPTION( projection_exception(err) ); + + if( dstdefn.fr_meter != 1.0 ) + { + for( std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(range, i); + if( ! is_invalid_point(point) ) + { + set<0>(point, get<0>(point) * dstdefn.fr_meter); + set<1>(point, get<1>(point) * dstdefn.fr_meter); + } + } + } + } + +/* -------------------------------------------------------------------- */ +/* Transform destination points to projection coordinates, if */ +/* desired. */ +/* -------------------------------------------------------------------- */ + else if( !dstdefn.is_latlong ) + { + + //if( dstdefn->fwd3d != NULL) + //{ + // for( i = 0; i < point_count; i++ ) + // { + // XYZ projected_loc; + // LPZ geodetic_loc; + + // geodetic_loc.u = x[point_offset*i]; + // geodetic_loc.v = y[point_offset*i]; + // geodetic_loc.w = z[point_offset*i]; + + // if (geodetic_loc.u == HUGE_VAL) + // continue; + + // projected_loc = pj_fwd3d( geodetic_loc, dstdefn); + // if( dstdefn->ctx->last_errno != 0 ) + // { + // if( (dstdefn->ctx->last_errno != 33 /*EDOM*/ + // && dstdefn->ctx->last_errno != 34 /*ERANGE*/ ) + // && (dstdefn->ctx->last_errno > 0 + // || dstdefn->ctx->last_errno < -44 || point_count == 1 + // || transient_error[-dstdefn->ctx->last_errno] == 0 ) ) + // return dstdefn->ctx->last_errno; + // else + // { + // projected_loc.u = HUGE_VAL; + // projected_loc.v = HUGE_VAL; + // projected_loc.w = HUGE_VAL; + // } + // } + + // x[point_offset*i] = projected_loc.u; + // y[point_offset*i] = projected_loc.v; + // z[point_offset*i] = projected_loc.w; + // } + + //} + //else + { + for(std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(range, i); + + if( is_invalid_point(point) ) + continue; + + try { + pj_fwd(dstprj, dstdefn, point, point); + } catch (projection_exception const& e) { + + if( (e.code() != 33 /*EDOM*/ + && e.code() != 34 /*ERANGE*/ ) + && (e.code() > 0 + || e.code() < -44 /*|| point_count == 1*/ + || transient_error[-e.code()] == 0) ) { + BOOST_RETHROW + } else { + set_invalid_point(point); + result = false; + if (point_count == 1) + return result; + } + } + } + } + } + +/* -------------------------------------------------------------------- */ +/* If a wrapping center other than 0 is provided, rewrap around */ +/* the suggested center (for latlong coordinate systems only). */ +/* -------------------------------------------------------------------- */ + else if( dstdefn.is_latlong && dstdefn.is_long_wrap_set ) + { + for( std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(range, i); + coord_t x = get_as_radian<0>(point); + + if( is_invalid_point(point) ) + continue; + + // TODO - units-dependant constants could be used instead + while( x < dstdefn.long_wrap_center - math::pi<coord_t>() ) + x += math::two_pi<coord_t>(); + while( x > dstdefn.long_wrap_center + math::pi<coord_t>() ) + x -= math::two_pi<coord_t>(); + + set_from_radian<0>(point, x); + } + } + +/* -------------------------------------------------------------------- */ +/* Transform Z from meters if needed. */ +/* -------------------------------------------------------------------- */ + if( dstdefn.vto_meter != 1.0 && dimension > 2 ) + { + for( std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = geometry::range::at(range, i); + set_z(point, get_z(point) * dstdefn.vfr_meter); + } + } + +/* -------------------------------------------------------------------- */ +/* Transform normalized axes into unusual output coordinate axis */ +/* orientation if needed. */ +/* -------------------------------------------------------------------- */ + // Ignored + + return result; +} + +/************************************************************************/ +/* pj_geodetic_to_geocentric() */ +/************************************************************************/ + +template <typename T, typename Range, bool AddZ> +inline int pj_geodetic_to_geocentric( T const& a, T const& es, + range_wrapper<Range, AddZ> & range_wrapper ) + +{ + //typedef typename boost::range_iterator<Range>::type iterator; + typedef typename boost::range_value<Range>::type point_type; + //typedef typename coordinate_type<point_type>::type coord_t; + + Range & rng = range_wrapper.get_range(); + std::size_t point_count = boost::size(rng); + + int ret_errno = 0; + + T const b = (es == 0.0) ? a : a * sqrt(1-es); + + GeocentricInfo<T> gi; + if( pj_Set_Geocentric_Parameters( gi, a, b ) != 0 ) + { + return PJD_ERR_GEOCENTRIC; + } + + for( std::size_t i = 0 ; i < point_count ; ++i ) + { + point_type & point = range::at(rng, i); + + if( is_invalid_point(point) ) + continue; + + T X = 0, Y = 0, Z = 0; + if( pj_Convert_Geodetic_To_Geocentric( gi, + get_as_radian<0>(point), + get_as_radian<1>(point), + range_wrapper.get_z(i), // Height + X, Y, Z ) != 0 ) + { + ret_errno = -14; + set_invalid_point(point); + /* but keep processing points! */ + } + else + { + set<0>(point, X); + set<1>(point, Y); + range_wrapper.set_z(i, Z); + } + } + + return ret_errno; +} + +/************************************************************************/ +/* pj_geodetic_to_geocentric() */ +/************************************************************************/ + +template <typename T, typename Range, bool AddZ> +inline int pj_geocentric_to_geodetic( T const& a, T const& es, + range_wrapper<Range, AddZ> & range_wrapper ) + +{ + //typedef typename boost::range_iterator<Range>::type iterator; + typedef typename boost::range_value<Range>::type point_type; + //typedef typename coordinate_type<point_type>::type coord_t; + + Range & rng = range_wrapper.get_range(); + std::size_t point_count = boost::size(rng); + + T const b = (es == 0.0) ? a : a * sqrt(1-es); + + GeocentricInfo<T> gi; + if( pj_Set_Geocentric_Parameters( gi, a, b ) != 0 ) + { + return PJD_ERR_GEOCENTRIC; + } + + for( std::size_t i = 0 ; i < point_count ; ++i ) + { + point_type & point = range::at(rng, i); + + if( is_invalid_point(point) ) + continue; + + T Longitude = 0, Latitude = 0, Height = 0; + pj_Convert_Geocentric_To_Geodetic( gi, + get<0>(point), + get<1>(point), + range_wrapper.get_z(i), // z + Longitude, Latitude, Height ); + + set_from_radian<0>(point, Longitude); + set_from_radian<1>(point, Latitude); + range_wrapper.set_z(i, Height); // Height + } + + return 0; +} + +/************************************************************************/ +/* pj_compare_datums() */ +/* */ +/* Returns TRUE if the two datums are identical, otherwise */ +/* FALSE. */ +/************************************************************************/ + +template <typename Par> +inline bool pj_compare_datums( Par & srcdefn, Par & dstdefn ) +{ + if( srcdefn.datum_type != dstdefn.datum_type ) + { + return false; + } + else if( srcdefn.a_orig != dstdefn.a_orig + || math::abs(srcdefn.es_orig - dstdefn.es_orig) > 0.000000000050 ) + { + /* the tolerance for es is to ensure that GRS80 and WGS84 are + considered identical */ + return false; + } + else if( srcdefn.datum_type == PJD_3PARAM ) + { + return (srcdefn.datum_params[0] == dstdefn.datum_params[0] + && srcdefn.datum_params[1] == dstdefn.datum_params[1] + && srcdefn.datum_params[2] == dstdefn.datum_params[2]); + } + else if( srcdefn.datum_type == PJD_7PARAM ) + { + return (srcdefn.datum_params[0] == dstdefn.datum_params[0] + && srcdefn.datum_params[1] == dstdefn.datum_params[1] + && srcdefn.datum_params[2] == dstdefn.datum_params[2] + && srcdefn.datum_params[3] == dstdefn.datum_params[3] + && srcdefn.datum_params[4] == dstdefn.datum_params[4] + && srcdefn.datum_params[5] == dstdefn.datum_params[5] + && srcdefn.datum_params[6] == dstdefn.datum_params[6]); + } + else if( srcdefn.datum_type == PJD_GRIDSHIFT ) + { + return pj_param(srcdefn.params,"snadgrids").s + == pj_param(dstdefn.params,"snadgrids").s; + } + else + return true; +} + +/************************************************************************/ +/* pj_geocentic_to_wgs84() */ +/************************************************************************/ + +template <typename Par, typename Range, bool AddZ> +inline int pj_geocentric_to_wgs84( Par const& defn, + range_wrapper<Range, AddZ> & range_wrapper ) + +{ + typedef typename boost::range_value<Range>::type point_type; + typedef typename coordinate_type<point_type>::type coord_t; + + Range & rng = range_wrapper.get_range(); + std::size_t point_count = boost::size(rng); + + if( defn.datum_type == PJD_3PARAM ) + { + for(std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(rng, i); + + if( is_invalid_point(point) ) + continue; + + set<0>(point, get<0>(point) + Dx_BF(defn)); + set<1>(point, get<1>(point) + Dy_BF(defn)); + range_wrapper.set_z(i, range_wrapper.get_z(i) + Dz_BF(defn)); + } + } + else if( defn.datum_type == PJD_7PARAM ) + { + for(std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(rng, i); + + if( is_invalid_point(point) ) + continue; + + coord_t x = get<0>(point); + coord_t y = get<1>(point); + coord_t z = range_wrapper.get_z(i); + + coord_t x_out, y_out, z_out; + + x_out = M_BF(defn)*( x - Rz_BF(defn)*y + Ry_BF(defn)*z) + Dx_BF(defn); + y_out = M_BF(defn)*( Rz_BF(defn)*x + y - Rx_BF(defn)*z) + Dy_BF(defn); + z_out = M_BF(defn)*(-Ry_BF(defn)*x + Rx_BF(defn)*y + z) + Dz_BF(defn); + + set<0>(point, x_out); + set<1>(point, y_out); + range_wrapper.set_z(i, z_out); + } + } + + return 0; +} + +/************************************************************************/ +/* pj_geocentic_from_wgs84() */ +/************************************************************************/ + +template <typename Par, typename Range, bool AddZ> +inline int pj_geocentric_from_wgs84( Par const& defn, + range_wrapper<Range, AddZ> & range_wrapper ) + +{ + typedef typename boost::range_value<Range>::type point_type; + typedef typename coordinate_type<point_type>::type coord_t; + + Range & rng = range_wrapper.get_range(); + std::size_t point_count = boost::size(rng); + + if( defn.datum_type == PJD_3PARAM ) + { + for(std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(rng, i); + + if( is_invalid_point(point) ) + continue; + + set<0>(point, get<0>(point) - Dx_BF(defn)); + set<1>(point, get<1>(point) - Dy_BF(defn)); + range_wrapper.set_z(i, range_wrapper.get_z(i) - Dz_BF(defn)); + } + } + else if( defn.datum_type == PJD_7PARAM ) + { + for(std::size_t i = 0; i < point_count; i++ ) + { + point_type & point = range::at(rng, i); + + if( is_invalid_point(point) ) + continue; + + coord_t x = get<0>(point); + coord_t y = get<1>(point); + coord_t z = range_wrapper.get_z(i); + + coord_t x_tmp = (x - Dx_BF(defn)) / M_BF(defn); + coord_t y_tmp = (y - Dy_BF(defn)) / M_BF(defn); + coord_t z_tmp = (z - Dz_BF(defn)) / M_BF(defn); + + x = x_tmp + Rz_BF(defn)*y_tmp - Ry_BF(defn)*z_tmp; + y = -Rz_BF(defn)*x_tmp + y_tmp + Rx_BF(defn)*z_tmp; + z = Ry_BF(defn)*x_tmp - Rx_BF(defn)*y_tmp + z_tmp; + + set<0>(point, x); + set<1>(point, y); + range_wrapper.set_z(i, z); + } + } + + return 0; +} + + +inline bool pj_datum_check_error(int err) +{ + return err != 0 && (err > 0 || transient_error[-err] == 0); +} + +/************************************************************************/ +/* pj_datum_transform() */ +/* */ +/* The input should be long/lat/z coordinates in radians in the */ +/* source datum, and the output should be long/lat/z */ +/* coordinates in radians in the destination datum. */ +/************************************************************************/ + +template <typename Par, typename Range> +inline bool pj_datum_transform( Par const& srcdefn, Par const& dstdefn, + Range & range ) + +{ + typedef typename Par::type calc_t; + bool result = true; + + calc_t src_a, src_es, dst_a, dst_es; + +/* -------------------------------------------------------------------- */ +/* We cannot do any meaningful datum transformation if either */ +/* the source or destination are of an unknown datum type */ +/* (ie. only a +ellps declaration, no +datum). This is new */ +/* behavior for PROJ 4.6.0. */ +/* -------------------------------------------------------------------- */ + if( srcdefn.datum_type == PJD_UNKNOWN + || dstdefn.datum_type == PJD_UNKNOWN ) + return result; + +/* -------------------------------------------------------------------- */ +/* Short cut if the datums are identical. */ +/* -------------------------------------------------------------------- */ + if( pj_compare_datums( srcdefn, dstdefn ) ) + return result; + + src_a = srcdefn.a_orig; + src_es = srcdefn.es_orig; + + dst_a = dstdefn.a_orig; + dst_es = dstdefn.es_orig; + +/* -------------------------------------------------------------------- */ +/* Create a temporary Z array if one is not provided. */ +/* -------------------------------------------------------------------- */ + + range_wrapper<Range> z_range(range); + +/* -------------------------------------------------------------------- */ +/* If this datum requires grid shifts, then apply it to geodetic */ +/* coordinates. */ +/* -------------------------------------------------------------------- */ + /*if( srcdefn.datum_type == PJD_GRIDSHIFT ) + { + try { + pj_apply_gridshift_2( srcdefn, 0, point_count, point_offset, x, y, z ); + } catch (projection_exception const& e) { + if (pj_datum_check_error(e.code())) { + BOOST_RETHROW + } + } + + src_a = SRS_WGS84_SEMIMAJOR; + src_es = SRS_WGS84_ESQUARED; + } + + if( dstdefn.datum_type == PJD_GRIDSHIFT ) + { + dst_a = SRS_WGS84_SEMIMAJOR; + dst_es = SRS_WGS84_ESQUARED; + }*/ + +/* ==================================================================== */ +/* Do we need to go through geocentric coordinates? */ +/* ==================================================================== */ + if( src_es != dst_es || src_a != dst_a + || srcdefn.datum_type == PJD_3PARAM + || srcdefn.datum_type == PJD_7PARAM + || dstdefn.datum_type == PJD_3PARAM + || dstdefn.datum_type == PJD_7PARAM) + { +/* -------------------------------------------------------------------- */ +/* Convert to geocentric coordinates. */ +/* -------------------------------------------------------------------- */ + int err = pj_geodetic_to_geocentric( src_a, src_es, z_range ); + if (pj_datum_check_error(err)) + BOOST_THROW_EXCEPTION( projection_exception(err) ); + else if (err != 0) + result = false; + +/* -------------------------------------------------------------------- */ +/* Convert between datums. */ +/* -------------------------------------------------------------------- */ + if( srcdefn.datum_type == PJD_3PARAM + || srcdefn.datum_type == PJD_7PARAM ) + { + try { + pj_geocentric_to_wgs84( srcdefn, z_range ); + } catch (projection_exception const& e) { + if (pj_datum_check_error(e.code())) { + BOOST_RETHROW + } + } + } + + if( dstdefn.datum_type == PJD_3PARAM + || dstdefn.datum_type == PJD_7PARAM ) + { + try { + pj_geocentric_from_wgs84( dstdefn, z_range ); + } catch (projection_exception const& e) { + if (pj_datum_check_error(e.code())) { + BOOST_RETHROW + } + } + } + +/* -------------------------------------------------------------------- */ +/* Convert back to geodetic coordinates. */ +/* -------------------------------------------------------------------- */ + err = pj_geocentric_to_geodetic( dst_a, dst_es, z_range ); + if (pj_datum_check_error(err)) + BOOST_THROW_EXCEPTION( projection_exception(err) ); + else if (err != 0) + result = false; + } + +/* -------------------------------------------------------------------- */ +/* Apply grid shift to destination if required. */ +/* -------------------------------------------------------------------- */ + /*if( dstdefn.datum_type == PJD_GRIDSHIFT ) + { + try { + pj_apply_gridshift_2( dstdefn, 1, point_count, point_offset, x, y, z ); + } catch (projection_exception const& e) { + if (pj_datum_check_error(e.code())) + BOOST_RETHROW + } + }*/ + + return result; +} + +} // namespace detail + +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP diff --git a/boost/geometry/srs/projections/impl/pj_tsfn.hpp b/boost/geometry/srs/projections/impl/pj_tsfn.hpp new file mode 100644 index 0000000000..d39a0fb2fc --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_tsfn.hpp @@ -0,0 +1,58 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_PJ_TSFN_HPP +#define BOOST_GEOMETRY_PROJECTIONS_PJ_TSFN_HPP + +#include <boost/geometry/util/math.hpp> + +namespace boost { namespace geometry { namespace projections { +namespace detail { + + /* determine small t */ + template <typename T> + inline T pj_tsfn(T const& phi, T sinphi, T const& e) + { + sinphi *= e; + return (tan (.5 * (geometry::math::half_pi<T>() - phi)) / + pow((1. - sinphi) / (1. + sinphi), .5 * e)); + } + +} // namespace detail +}}} // namespace boost::geometry::projections +#endif diff --git a/boost/geometry/srs/projections/impl/pj_units.hpp b/boost/geometry/srs/projections/impl/pj_units.hpp new file mode 100644 index 0000000000..269b8ff92e --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_units.hpp @@ -0,0 +1,79 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_UNITS_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_UNITS_HPP + +#include <boost/geometry/srs/projections/impl/projects.hpp> + +namespace boost { namespace geometry { namespace projections { +namespace detail { + +/* Field 2 that contains the multiplier to convert named units to meters +** may be expressed by either a simple floating point constant or a +** numerator/denomenator values (e.g. 1/1000) */ + +static const PJ_UNITS pj_units[] = +{ + { "km", "1000.", "Kilometer" }, + { "m", "1.", "Meter" }, + { "dm", "1/10", "Decimeter" }, + { "cm", "1/100", "Centimeter" }, + { "mm", "1/1000", "Millimeter" }, + { "kmi", "1852.0", "International Nautical Mile" }, + { "in", "0.0254", "International Inch" }, + { "ft", "0.3048", "International Foot" }, + { "yd", "0.9144", "International Yard" }, + { "mi", "1609.344", "International Statute Mile" }, + { "fath", "1.8288", "International Fathom" }, + { "ch", "20.1168", "International Chain" }, + { "link", "0.201168", "International Link" }, + { "us-in", "1./39.37", "U.S. Surveyor's Inch" }, + { "us-ft", "0.304800609601219", "U.S. Surveyor's Foot" }, + { "us-yd", "0.914401828803658", "U.S. Surveyor's Yard" }, + { "us-ch", "20.11684023368047", "U.S. Surveyor's Chain" }, + { "us-mi", "1609.347218694437", "U.S. Surveyor's Statute Mile" }, + { "ind-yd", "0.91439523", "Indian Yard" }, + { "ind-ft", "0.30479841", "Indian Foot" }, + { "ind-ch", "20.11669506", "Indian Chain" } +}; + +} // detail +}}} // namespace boost::geometry::projections + +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_UNITS_HPP diff --git a/boost/geometry/srs/projections/impl/pj_zpoly1.hpp b/boost/geometry/srs/projections/impl/pj_zpoly1.hpp new file mode 100644 index 0000000000..613b6d2b35 --- /dev/null +++ b/boost/geometry/srs/projections/impl/pj_zpoly1.hpp @@ -0,0 +1,106 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_ZPOLY1_HPP +#define BOOST_GEOMETRY_PROJECTIONS_ZPOLY1_HPP + + +#include <boost/geometry/srs/projections/impl/projects.hpp> + + +namespace boost { namespace geometry { namespace projections { namespace detail { + + /* evaluate complex polynomial */ + + /* note: coefficients are always from C_1 to C_n + ** i.e. C_0 == (0., 0) + ** n should always be >= 1 though no checks are made + */ + template <typename T> + inline COMPLEX<T> + pj_zpoly1(COMPLEX<T> z, const COMPLEX<T> *C, int n) + { + COMPLEX<T> a; + T t; + + a = *(C += n); + while (n-- > 0) + { + a.r = (--C)->r + z.r * (t = a.r) - z.i * a.i; + a.i = C->i + z.r * a.i + z.i * t; + } + a.r = z.r * (t = a.r) - z.i * a.i; + a.i = z.r * a.i + z.i * t; + return a; + } + + /* evaluate complex polynomial and derivative */ + template <typename T> + inline COMPLEX<T> + pj_zpolyd1(COMPLEX<T> z, const COMPLEX<T> *C, int n, COMPLEX<T> *der) + { + T t; + bool first = true; + + COMPLEX<T> a = *(C += n); + COMPLEX<T> b = a; + while (n-- > 0) + { + if (first) + { + first = false; + } + else + { + b.r = a.r + z.r * (t = b.r) - z.i * b.i; + b.i = a.i + z.r * b.i + z.i * t; + } + a.r = (--C)->r + z.r * (t = a.r) - z.i * a.i; + a.i = C->i + z.r * a.i + z.i * t; + } + b.r = a.r + z.r * (t = b.r) - z.i * b.i; + b.i = a.i + z.r * b.i + z.i * t; + a.r = z.r * (t = a.r) - z.i * a.i; + a.i = z.r * a.i + z.i * t; + *der = b; + return a; + } + +}}}} // namespace boost::geometry::projections::detail + +#endif diff --git a/boost/geometry/srs/projections/impl/proj_mdist.hpp b/boost/geometry/srs/projections/impl/proj_mdist.hpp new file mode 100644 index 0000000000..1c325df1c6 --- /dev/null +++ b/boost/geometry/srs/projections/impl/proj_mdist.hpp @@ -0,0 +1,144 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_PROJ_MDIST_HPP +#define BOOST_GEOMETRY_PROJECTIONS_PROJ_MDIST_HPP + + +#include <boost/geometry/util/math.hpp> + + +namespace boost { namespace geometry { namespace projections +{ +namespace detail +{ + static const int MDIST_MAX_ITER = 20; + + template <typename T> + struct MDIST + { + int nb; + T es; + T E; + T b[MDIST_MAX_ITER]; + }; + + template <typename CT> + inline bool proj_mdist_ini(CT const& es, MDIST<CT>& b) + { + CT numf, numfi, twon1, denf, denfi, ens, T, twon; + CT den, El, Es; + CT E[MDIST_MAX_ITER]; + int i, j; + + /* generate E(e^2) and its terms E[] */ + ens = es; + numf = twon1 = denfi = 1.; + denf = 1.; + twon = 4.; + Es = El = E[0] = 1.; + for (i = 1; i < MDIST_MAX_ITER ; ++i) + { + numf *= (twon1 * twon1); + den = twon * denf * denf * twon1; + T = numf/den; + Es -= (E[i] = T * ens); + ens *= es; + twon *= 4.; + denf *= ++denfi; + twon1 += 2.; + if (Es == El) /* jump out if no change */ + break; + El = Es; + } + b.nb = i - 1; + b.es = es; + b.E = Es; + /* generate b_n coefficients--note: collapse with prefix ratios */ + b.b[0] = Es = 1. - Es; + numf = denf = 1.; + numfi = 2.; + denfi = 3.; + for (j = 1; j < i; ++j) + { + Es -= E[j]; + numf *= numfi; + denf *= denfi; + b.b[j] = Es * numf / denf; + numfi += 2.; + denfi += 2.; + } + return true; + } + + template <typename T> + inline T proj_mdist(T const& phi, T const& sphi, T const& cphi, MDIST<T> const& b) + { + T sc, sum, sphi2, D; + int i; + + sc = sphi * cphi; + sphi2 = sphi * sphi; + D = phi * b.E - b.es * sc / sqrt(1. - b.es * sphi2); + sum = b.b[i = b.nb]; + while (i) sum = b.b[--i] + sphi2 * sum; + return(D + sc * sum); + } + + template <typename T> + inline T proj_inv_mdist(T const& dist, MDIST<T> const& b) + { + static const T TOL = 1e-14; + T s, t, phi, k; + int i; + + k = 1./(1.- b.es); + i = MDIST_MAX_ITER; + phi = dist; + while ( i-- ) { + s = sin(phi); + t = 1. - b.es * s * s; + phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) * + (t * sqrt(t)) * k; + if (geometry::math::abs(t) < TOL) /* that is no change */ + return phi; + } + /* convergence failed */ + BOOST_THROW_EXCEPTION( projection_exception(-17) ); + } +} // namespace detail + +}}} // namespace boost::geometry::projections + +#endif diff --git a/boost/geometry/srs/projections/impl/projects.hpp b/boost/geometry/srs/projections/impl/projects.hpp new file mode 100644 index 0000000000..232ae67ae9 --- /dev/null +++ b/boost/geometry/srs/projections/impl/projects.hpp @@ -0,0 +1,272 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// This file is manually converted from PROJ4 (projects.h) + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// 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 Geometry Library by Barend Gehrels (Geodan, Amsterdam) + +// 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. + +#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PROJECTS_HPP +#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PROJECTS_HPP + + +#include <cstring> +#include <string> +#include <vector> + +#include <boost/geometry/srs/projections/exception.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/mpl/if.hpp> +#include <boost/type_traits/is_pod.hpp> + + +namespace boost { namespace geometry { namespace projections +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +/* some useful constants */ +template <typename T> +inline T ONEPI() { return boost::math::constants::pi<T>(); } +template <typename T> +inline T HALFPI() { return boost::math::constants::half_pi<T>(); } +template <typename T> +inline T FORTPI() { return boost::math::constants::pi<T>() / T(4); } +template <typename T> +inline T TWOPI() { return boost::math::constants::two_pi<T>(); } +template <typename T> +inline T TWO_D_PI() { return boost::math::constants::two_div_pi<T>(); } +template <typename T> +inline T HALFPI_SQR() { return 2.4674011002723396547086227499689; } +template <typename T> +inline T PI_SQR() { return boost::math::constants::pi_sqr<T>(); } +template <typename T> +inline T THIRD() { return 0.3333333333333333333333333333333; } +template <typename T> +inline T TWOTHIRD() { return 0.6666666666666666666666666666666; } +template <typename T> +inline T PI_HALFPI() { return 4.7123889803846898576939650749193; } +template <typename T> +inline T TWOPI_HALFPI() { return 7.8539816339744830961566084581988; } +template <typename T> +inline T PI_DIV_3() { return 1.0471975511965977461542144610932; } + +/* datum_type values */ +static const int PJD_UNKNOWN = 0; +static const int PJD_3PARAM = 1; +static const int PJD_7PARAM = 2; +static const int PJD_GRIDSHIFT = 3; +static const int PJD_WGS84 = 4; /* WGS84 (or anything considered equivelent) */ + +/* library errors */ +static const int PJD_ERR_GEOCENTRIC = -45; +static const int PJD_ERR_AXIS = -47; +static const int PJD_ERR_GRID_AREA = -48; +static const int PJD_ERR_CATALOG = -49; + +template <typename T> +struct pvalue +{ + std::string param; + int used; + + int i; + T f; + std::string s; +}; + +template <typename T> +struct pj_const_pod +{ + int over; /* over-range flag */ + int geoc; /* geocentric latitude flag */ + int is_latlong; /* proj=latlong ... not really a projection at all */ + int is_geocent; /* proj=geocent ... not really a projection at all */ + T + a, /* major axis or radius if es==0 */ + a_orig, /* major axis before any +proj related adjustment */ + es, /* e ^ 2 */ + es_orig, /* es before any +proj related adjustment */ + e, /* eccentricity */ + ra, /* 1/A */ + one_es, /* 1 - e^2 */ + rone_es, /* 1/one_es */ + lam0, phi0, /* central longitude, latitude */ + x0, y0, /* easting and northing */ + k0, /* general scaling factor */ + to_meter, fr_meter, /* cartesian scaling */ + vto_meter, vfr_meter; /* Vertical scaling. Internal unit [m] */ + + int datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */ + T datum_params[7]; + T from_greenwich; /* prime meridian offset (in radians) */ + T long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/ + bool is_long_wrap_set; + + // Initialize all variables to zero + pj_const_pod() + { + std::memset(this, 0, sizeof(pj_const_pod)); + } +}; + +template <typename T> +struct pj_const_non_pod +{ + int over; /* over-range flag */ + int geoc; /* geocentric latitude flag */ + int is_latlong; /* proj=latlong ... not really a projection at all */ + int is_geocent; /* proj=geocent ... not really a projection at all */ + T + a, /* major axis or radius if es==0 */ + a_orig, /* major axis before any +proj related adjustment */ + es, /* e ^ 2 */ + es_orig, /* es before any +proj related adjustment */ + e, /* eccentricity */ + ra, /* 1/A */ + one_es, /* 1 - e^2 */ + rone_es, /* 1/one_es */ + lam0, phi0, /* central longitude, latitude */ + x0, y0, /* easting and northing */ + k0, /* general scaling factor */ + to_meter, fr_meter, /* cartesian scaling */ + vto_meter, vfr_meter; /* Vertical scaling. Internal unit [m] */ + + int datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */ + T datum_params[7]; + T from_greenwich; /* prime meridian offset (in radians) */ + T long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/ + bool is_long_wrap_set; + + // Initialize all variables to zero + pj_const_non_pod() + : over(0), geoc(0), is_latlong(0), is_geocent(0) + , a(0), a_orig(0), es(0), es_orig(0), e(0), ra(0) + , one_es(0), rone_es(0), lam0(0), phi0(0), x0(0), y0(0), k0(0) + , to_meter(0), fr_meter(0), vto_meter(0), vfr_meter(0) + , datum_type(PJD_UNKNOWN) + , from_greenwich(0), long_wrap_center(0), is_long_wrap_set(false) + { + datum_params[0] = 0; + datum_params[1] = 0; + datum_params[2] = 0; + datum_params[3] = 0; + datum_params[4] = 0; + datum_params[5] = 0; + datum_params[6] = 0; + } +}; + +template <typename T> +struct pj_const + : boost::mpl::if_c + < + boost::is_pod<T>::value, + pj_const_pod<T>, + pj_const_non_pod<T> + >::type +{}; + +// PROJ4 complex. Might be replaced with std::complex +template <typename T> +struct COMPLEX { T r, i; }; + +struct PJ_ELLPS +{ + std::string id; /* ellipse keyword name */ + std::string major; /* a= value */ + std::string ell; /* elliptical parameter */ + std::string name; /* comments */ +}; + +struct PJ_DATUMS +{ + std::string id; /* datum keyword */ + std::string defn; /* ie. "to_wgs84=..." */ + std::string ellipse_id; /* ie from ellipse table */ + std::string comments; /* EPSG code, etc */ +}; + +struct PJ_PRIME_MERIDIANS +{ + std::string id; /* prime meridian keyword */ + std::string defn; /* offset from greenwich in DMS format. */ +}; + +struct PJ_UNITS +{ + std::string id; /* units keyword */ + std::string to_meter; /* multiply by value to get meters */ + std::string name; /* comments */ +}; + +template <typename T> +struct DERIVS +{ + T x_l, x_p; /* derivatives of x for lambda-phi */ + T y_l, y_p; /* derivatives of y for lambda-phi */ +}; + +template <typename T> +struct FACTORS +{ + DERIVS<T> der; + T h, k; /* meridinal, parallel scales */ + T omega, thetap; /* angular distortion, theta prime */ + T conv; /* convergence */ + T s; /* areal scale factor */ + T a, b; /* max-min scale error */ + int code; /* info as to analytics, see following */ +}; + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + +/*! + \brief parameters, projection parameters + \details This structure initializes all projections + \ingroup projection +*/ +template <typename T> +struct parameters : public detail::pj_const<T> +{ + typedef T type; + + std::string name; + std::vector<detail::pvalue<T> > params; +}; + +}}} // namespace boost::geometry::projections +#endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PROJECTS_HPP |