1 // Boost.Geometry - gis-projections (based on PROJ4)
3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
5 // This file was modified by Oracle on 2017, 2018, 2019.
6 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
7 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
9 // Use, modification and distribution is subject to the Boost Software License,
10 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
11 // http://www.boost.org/LICENSE_1_0.txt)
13 // This file is converted from PROJ4, http://trac.osgeo.org/proj
14 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
15 // PROJ4 is maintained by Frank Warmerdam
16 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
18 // Last updated version of proj: 5.0.0
20 // Original copyright notice:
22 // Purpose: Implementation of the aitoff (Aitoff) and wintri (Winkel Tripel)
24 // Author: Gerald Evenden (1995)
25 // Drazen Tutic, Lovro Gradiser (2015) - add inverse
26 // Thomas Knudsen (2016) - revise/add regression tests
27 // Copyright (c) 1995, Gerald Evenden
29 // Permission is hereby granted, free of charge, to any person obtaining a
30 // copy of this software and associated documentation files (the "Software"),
31 // to deal in the Software without restriction, including without limitation
32 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
33 // and/or sell copies of the Software, and to permit persons to whom the
34 // Software is furnished to do so, subject to the following conditions:
36 // The above copyright notice and this permission notice shall be included
37 // in all copies or substantial portions of the Software.
39 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
40 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
41 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
42 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
43 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
44 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
45 // DEALINGS IN THE SOFTWARE.
47 #ifndef BOOST_GEOMETRY_PROJECTIONS_AITOFF_HPP
48 #define BOOST_GEOMETRY_PROJECTIONS_AITOFF_HPP
50 #include <boost/core/ignore_unused.hpp>
52 #include <boost/geometry/srs/projections/impl/base_static.hpp>
53 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
54 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
55 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
56 #include <boost/geometry/srs/projections/impl/projects.hpp>
58 #include <boost/geometry/util/math.hpp>
60 namespace boost { namespace geometry
65 #ifndef DOXYGEN_NO_DETAIL
66 namespace detail { namespace aitoff
70 mode_winkel_tripel = 1
80 template <typename T, typename Parameters>
81 struct base_aitoff_spheroid
83 par_aitoff<T> m_proj_parm;
85 // FORWARD(s_forward) spheroid
86 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
87 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
91 if((d = acos(cos(lp_lat) * cos(c = 0.5 * lp_lon)))) {/* basic Aitoff */
92 xy_x = 2. * d * cos(lp_lat) * sin(c) * (xy_y = 1. / sin(d));
93 xy_y *= d * sin(lp_lat);
96 if (this->m_proj_parm.mode == mode_winkel_tripel) { /* Winkel Tripel */
97 xy_x = (xy_x + lp_lon * this->m_proj_parm.cosphi1) * 0.5;
98 xy_y = (xy_y + lp_lat) * 0.5;
101 /***********************************************************************************
103 * Inverse functions added by Drazen Tutic and Lovro Gradiser based on paper:
105 * I.Özbug Biklirici and Cengizhan Ipbüker. A General Algorithm for the Inverse
106 * Transformation of Map Projections Using Jacobian Matrices. In Proceedings of the
107 * Third International Symposium Mathematical & Computational Applications,
108 * pages 175{182, Turkey, September 2002.
110 * Expected accuracy is defined by epsilon = 1e-12. Should be appropriate for
111 * most applications of Aitoff and Winkel Tripel projections.
113 * Longitudes of 180W and 180E can be mixed in solution obtained.
115 * Inverse for Aitoff projection in poles is undefined, longitude value of 0 is assumed.
117 * Contact : dtutic@geof.hr
120 ************************************************************************************/
122 // INVERSE(s_inverse) sphere
123 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
124 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
126 static const T pi = detail::pi<T>();
127 static const T two_pi = detail::two_pi<T>();
128 static const T epsilon = 1e-12;
130 int iter, max_iter = 10, round = 0, max_round = 20;
131 T D, C, f1, f2, f1p, f1l, f2p, f2l, dp, dl, sl, sp, cp, cl, x, y;
133 if ((fabs(xy_x) < epsilon) && (fabs(xy_y) < epsilon )) {
134 lp_lat = 0.; lp_lon = 0.;
138 /* intial values for Newton-Raphson method */
139 lp_lat = xy_y; lp_lon = xy_x;
143 sl = sin(lp_lon * 0.5); cl = cos(lp_lon * 0.5);
144 sp = sin(lp_lat); cp = cos(lp_lat);
147 D = acos(D) / math::pow(C, T(1.5));
148 f1 = 2. * D * C * cp * sl;
150 f1p = 2.* (sl * cl * sp * cp / C - D * sp * sl);
151 f1l = cp * cp * sl * sl / C + D * cp * cl * sp * sp;
152 f2p = sp * sp * cl / C + D * sl * sl * cp;
153 f2l = 0.5 * (sp * cp * sl / C - D * sp * cp * cp * sl * cl);
154 if (this->m_proj_parm.mode == mode_winkel_tripel) { /* Winkel Tripel */
155 f1 = 0.5 * (f1 + lp_lon * this->m_proj_parm.cosphi1);
156 f2 = 0.5 * (f2 + lp_lat);
158 f1l = 0.5 * (f1l + this->m_proj_parm.cosphi1);
159 f2p = 0.5 * (f2p + 1.);
162 f1 -= xy_x; f2 -= xy_y;
163 dl = (f2 * f1p - f1 * f2p) / (dp = f1p * f2l - f2p * f1l);
164 dp = (f1 * f2l - f2 * f1l) / dp;
165 dl = fmod(dl, pi); /* set to interval [-M_PI, M_PI] */
166 lp_lat -= dp; lp_lon -= dl;
167 } while ((fabs(dp) > epsilon || fabs(dl) > epsilon) && (iter++ < max_iter));
168 if (lp_lat > two_pi) lp_lat -= 2.*(lp_lat-two_pi); /* correct if symmetrical solution for Aitoff */
169 if (lp_lat < -two_pi) lp_lat -= 2.*(lp_lat+two_pi); /* correct if symmetrical solution for Aitoff */
170 if ((fabs(fabs(lp_lat) - two_pi) < epsilon) && (!this->m_proj_parm.mode)) lp_lon = 0.; /* if pole in Aitoff, return longitude of 0 */
172 /* calculate x,y coordinates with solution obtained */
173 if((D = acos(cos(lp_lat) * cos(C = 0.5 * lp_lon))) != 0.0) {/* Aitoff */
174 x = 2. * D * cos(lp_lat) * sin(C) * (y = 1. / sin(D));
175 y *= D * sin(lp_lat);
178 if (this->m_proj_parm.mode == mode_winkel_tripel) { /* Winkel Tripel */
179 x = (x + lp_lon * this->m_proj_parm.cosphi1) * 0.5;
180 y = (y + lp_lat) * 0.5;
182 /* if too far from given values of x,y, repeat with better approximation of phi,lam */
183 } while (((fabs(xy_x-x) > epsilon) || (fabs(xy_y-y) > epsilon)) && (round++ < max_round));
185 if (iter == max_iter && round == max_round)
187 BOOST_THROW_EXCEPTION( projection_exception(error_non_convergent) );
188 //fprintf(stderr, "Warning: Accuracy of 1e-12 not reached. Last increments: dlat=%e and dlon=%e\n", dp, dl);
192 static inline std::string get_name()
194 return "aitoff_spheroid";
199 template <typename Parameters>
200 inline void setup(Parameters& par)
207 template <typename Parameters, typename T>
208 inline void setup_aitoff(Parameters& par, par_aitoff<T>& proj_parm)
210 proj_parm.mode = mode_aitoff;
215 template <typename Params, typename Parameters, typename T>
216 inline void setup_wintri(Params& params, Parameters& par, par_aitoff<T>& proj_parm)
218 static const T two_div_pi = detail::two_div_pi<T>();
222 proj_parm.mode = mode_winkel_tripel;
223 if (pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, phi1)) {
224 if ((proj_parm.cosphi1 = cos(phi1)) == 0.)
225 BOOST_THROW_EXCEPTION( projection_exception(error_lat_larger_than_90) );
226 } else /* 50d28' or phi1=acos(2/pi) */
227 proj_parm.cosphi1 = two_div_pi;
231 }} // namespace detail::aitoff
235 \brief Aitoff projection
237 \tparam Geographic latlong point type
238 \tparam Cartesian xy point type
239 \tparam Parameters parameter type
240 \par Projection characteristics
244 \image html ex_aitoff.gif
246 template <typename T, typename Parameters>
247 struct aitoff_spheroid : public detail::aitoff::base_aitoff_spheroid<T, Parameters>
249 template <typename Params>
250 inline aitoff_spheroid(Params const& , Parameters & par)
252 detail::aitoff::setup_aitoff(par, this->m_proj_parm);
257 \brief Winkel Tripel projection
259 \tparam Geographic latlong point type
260 \tparam Cartesian xy point type
261 \tparam Parameters parameter type
262 \par Projection characteristics
265 \par Projection parameters
266 - lat_1: Latitude of first standard parallel (degrees)
268 \image html ex_wintri.gif
270 template <typename T, typename Parameters>
271 struct wintri_spheroid : public detail::aitoff::base_aitoff_spheroid<T, Parameters>
273 template <typename Params>
274 inline wintri_spheroid(Params const& params, Parameters & par)
276 detail::aitoff::setup_wintri(params, par, this->m_proj_parm);
280 #ifndef DOXYGEN_NO_DETAIL
285 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_aitoff, aitoff_spheroid)
286 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_wintri, wintri_spheroid)
289 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(aitoff_entry, aitoff_spheroid)
290 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(wintri_entry, wintri_spheroid)
292 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(aitoff_init)
294 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(aitoff, aitoff_entry)
295 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(wintri, wintri_entry)
298 } // namespace detail
301 } // namespace projections
303 }} // namespace boost::geometry
305 #endif // BOOST_GEOMETRY_PROJECTIONS_AITOFF_HPP