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 // Permission is hereby granted, free of charge, to any person obtaining a
23 // copy of this software and associated documentation files (the "Software"),
24 // to deal in the Software without restriction, including without limitation
25 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
26 // and/or sell copies of the Software, and to permit persons to whom the
27 // Software is furnished to do so, subject to the following conditions:
29 // The above copyright notice and this permission notice shall be included
30 // in all copies or substantial portions of the Software.
32 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
33 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
34 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
35 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
36 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
37 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
38 // DEALINGS IN THE SOFTWARE.
40 #ifndef BOOST_GEOMETRY_PROJECTIONS_MOD_STER_HPP
41 #define BOOST_GEOMETRY_PROJECTIONS_MOD_STER_HPP
43 #include <boost/geometry/util/math.hpp>
44 #include <boost/math/special_functions/hypot.hpp>
46 #include <boost/geometry/srs/projections/impl/base_static.hpp>
47 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
48 #include <boost/geometry/srs/projections/impl/projects.hpp>
49 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
50 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
51 #include <boost/geometry/srs/projections/impl/pj_zpoly1.hpp>
53 namespace boost { namespace geometry
58 #ifndef DOXYGEN_NO_DETAIL
59 namespace detail { namespace mod_ster
62 static const double epsilon = 1e-12;
68 pj_complex<T>* zcoeff;
72 /* based upon Snyder and Linck, USGS-NMD */
74 template <typename T, typename Parameters>
75 struct base_mod_ster_ellipsoid
77 par_mod_ster<T> m_proj_parm;
79 // FORWARD(e_forward) ellipsoid
80 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
81 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
83 static const T half_pi = detail::half_pi<T>();
85 T sinlon, coslon, esphi, chi, schi, cchi, s;
90 esphi = par.e * sin(lp_lat);
91 chi = 2. * atan(tan((half_pi + lp_lat) * .5) *
92 math::pow((T(1) - esphi) / (T(1) + esphi), par.e * T(0.5))) - half_pi;
95 s = 2. / (1. + this->m_proj_parm.schio * schi + this->m_proj_parm.cchio * cchi * coslon);
96 p.r = s * cchi * sinlon;
97 p.i = s * (this->m_proj_parm.cchio * schi - this->m_proj_parm.schio * cchi * coslon);
98 p = pj_zpoly1(p, this->m_proj_parm.zcoeff, this->m_proj_parm.n);
103 // INVERSE(e_inverse) ellipsoid
104 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
105 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
107 static const T half_pi = detail::half_pi<T>();
110 pj_complex<T> p, fxy, fpxy, dp;
111 T den, rh = 0, z, sinz = 0, cosz = 0, chi, phi = 0, dphi, esphi;
115 for (nn = 20; nn ;--nn) {
116 fxy = pj_zpolyd1(p, this->m_proj_parm.zcoeff, this->m_proj_parm.n, &fpxy);
119 den = fpxy.r * fpxy.r + fpxy.i * fpxy.i;
120 dp.r = -(fxy.r * fpxy.r + fxy.i * fpxy.i) / den;
121 dp.i = -(fxy.i * fpxy.r - fxy.r * fpxy.i) / den;
124 if ((fabs(dp.r) + fabs(dp.i)) <= epsilon)
128 rh = boost::math::hypot(p.r, p.i);
129 z = 2. * atan(.5 * rh);
133 if (fabs(rh) <= epsilon) {
134 /* if we end up here input coordinates were (0,0).
135 * pj_inv() adds P->lam0 to lp.lam, this way we are
136 * sure to get the correct offset */
141 chi = aasin(cosz * this->m_proj_parm.schio + p.i * sinz * this->m_proj_parm.cchio / rh);
143 for (nn = 20; nn ;--nn) {
144 esphi = par.e * sin(phi);
145 dphi = 2. * atan(tan((half_pi + chi) * .5) *
146 math::pow((T(1) + esphi) / (T(1) - esphi), par.e * T(0.5))) - half_pi - phi;
148 if (fabs(dphi) <= epsilon)
154 lp_lon = atan2(p.r * sinz, rh * this->m_proj_parm.cchio * cosz - p.i *
155 this->m_proj_parm.schio * sinz);
157 lp_lon = lp_lat = HUGE_VAL;
160 static inline std::string get_name()
162 return "mod_ster_ellipsoid";
167 template <typename Parameters, typename T>
168 inline void setup(Parameters& par, par_mod_ster<T>& proj_parm) /* general initialization */
170 static T const half_pi = detail::half_pi<T>();
175 esphi = par.e * sin(par.phi0);
176 chio = 2. * atan(tan((half_pi + par.phi0) * .5) *
177 math::pow((T(1) - esphi) / (T(1) + esphi), par.e * T(0.5))) - half_pi;
180 proj_parm.schio = sin(chio);
181 proj_parm.cchio = cos(chio);
185 /* Miller Oblated Stereographic */
186 template <typename Parameters, typename T>
187 inline void setup_mil_os(Parameters& par, par_mod_ster<T>& proj_parm)
189 static const T d2r = geometry::math::d2r<T>();
191 static pj_complex<T> AB[] = {
198 par.lam0 = d2r * 20.;
199 par.phi0 = d2r * 18.;
200 proj_parm.zcoeff = AB;
203 setup(par, proj_parm);
206 /* Lee Oblated Stereographic */
207 template <typename Parameters, typename T>
208 inline void setup_lee_os(Parameters& par, par_mod_ster<T>& proj_parm)
210 static const T d2r = geometry::math::d2r<T>();
212 static pj_complex<T> AB[] = {
215 {-0.0088162, -0.00617325}
219 par.lam0 = d2r * -165.;
220 par.phi0 = d2r * -10.;
221 proj_parm.zcoeff = AB;
224 setup(par, proj_parm);
227 // Mod. Stererographics of 48 U.S.
228 template <typename Parameters, typename T>
229 inline void setup_gs48(Parameters& par, par_mod_ster<T>& proj_parm)
231 static const T d2r = geometry::math::d2r<T>();
233 static pj_complex<T> AB[] = { /* 48 United States */
242 par.lam0 = d2r * -96.;
243 par.phi0 = d2r * -39.;
244 proj_parm.zcoeff = AB;
248 setup(par, proj_parm);
251 // Mod. Stererographics of Alaska
252 template <typename Parameters, typename T>
253 inline void setup_alsk(Parameters& par, par_mod_ster<T>& proj_parm)
255 static const T d2r = geometry::math::d2r<T>();
257 static pj_complex<T> ABe[] = { /* Alaska ellipsoid */
259 { .0052083, -.0027404},
260 { .0072721, .0048181},
261 {-.0151089, -.1932526},
262 { .0642675, -.1381226},
263 { .3582802, -.2884586}
266 static pj_complex<T> ABs[] = { /* Alaska sphere */
268 { .0052513, -.0041175},
269 { .0074606, .0048125},
270 {-.0153783, -.1968253},
271 { .0636871, -.1408027},
272 { .3660976, -.2937382}
276 par.lam0 = d2r * -152.;
277 par.phi0 = d2r * 64.;
278 if (par.es != 0.0) { /* fixed ellipsoid/sphere */
279 proj_parm.zcoeff = ABe;
281 par.e = sqrt(par.es = 0.00676866);
283 proj_parm.zcoeff = ABs;
287 setup(par, proj_parm);
290 // Mod. Stererographics of 50 U.S.
291 template <typename Parameters, typename T>
292 inline void setup_gs50(Parameters& par, par_mod_ster<T>& proj_parm)
294 static const T d2r = geometry::math::d2r<T>();
296 static pj_complex<T> ABe[] = { /* GS50 ellipsoid */
298 { .0210669, .0053804},
299 {-.1031415, -.0571664},
300 {-.0323337, -.0322847},
301 { .0502303, .1211983},
302 { .0251805, .0895678},
303 {-.0012315, -.1416121},
304 { .0072202, -.1317091},
305 {-.0194029, .0759677},
306 {-.0210072, .0834037}
308 static pj_complex<T> ABs[] = { /* GS50 sphere */
310 { .0211642, .0037608},
311 {-.1036018, -.0575102},
312 {-.0329095, -.0320119},
313 { .0499471, .1223335},
314 { .0260460, .0899805},
315 { .0007388, -.1435792},
316 { .0075848, -.1334108},
317 {-.0216473, .0776645},
318 {-.0225161, .0853673}
322 par.lam0 = d2r * -120.;
323 par.phi0 = d2r * 45.;
324 if (par.es != 0.0) { /* fixed ellipsoid/sphere */
325 proj_parm.zcoeff = ABe;
327 par.e = sqrt(par.es = 0.00676866);
329 proj_parm.zcoeff = ABs;
333 setup(par, proj_parm);
336 }} // namespace detail::mod_ster
340 \brief Miller Oblated Stereographic projection
342 \tparam Geographic latlong point type
343 \tparam Cartesian xy point type
344 \tparam Parameters parameter type
345 \par Projection characteristics
348 \image html ex_mil_os.gif
350 template <typename T, typename Parameters>
351 struct mil_os_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
353 template <typename Params>
354 inline mil_os_ellipsoid(Params const& , Parameters & par)
356 detail::mod_ster::setup_mil_os(par, this->m_proj_parm);
361 \brief Lee Oblated Stereographic projection
363 \tparam Geographic latlong point type
364 \tparam Cartesian xy point type
365 \tparam Parameters parameter type
366 \par Projection characteristics
369 \image html ex_lee_os.gif
371 template <typename T, typename Parameters>
372 struct lee_os_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
374 template <typename Params>
375 inline lee_os_ellipsoid(Params const& , Parameters & par)
377 detail::mod_ster::setup_lee_os(par, this->m_proj_parm);
382 \brief Mod. Stererographics of 48 U.S. projection
384 \tparam Geographic latlong point type
385 \tparam Cartesian xy point type
386 \tparam Parameters parameter type
387 \par Projection characteristics
390 \image html ex_gs48.gif
392 template <typename T, typename Parameters>
393 struct gs48_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
395 template <typename Params>
396 inline gs48_ellipsoid(Params const& , Parameters & par)
398 detail::mod_ster::setup_gs48(par, this->m_proj_parm);
403 \brief Mod. Stererographics of Alaska projection
405 \tparam Geographic latlong point type
406 \tparam Cartesian xy point type
407 \tparam Parameters parameter type
408 \par Projection characteristics
411 \image html ex_alsk.gif
413 template <typename T, typename Parameters>
414 struct alsk_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
416 template <typename Params>
417 inline alsk_ellipsoid(Params const& , Parameters & par)
419 detail::mod_ster::setup_alsk(par, this->m_proj_parm);
424 \brief Mod. Stererographics of 50 U.S. projection
426 \tparam Geographic latlong point type
427 \tparam Cartesian xy point type
428 \tparam Parameters parameter type
429 \par Projection characteristics
432 \image html ex_gs50.gif
434 template <typename T, typename Parameters>
435 struct gs50_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
437 template <typename Params>
438 inline gs50_ellipsoid(Params const& , Parameters & par)
440 detail::mod_ster::setup_gs50(par, this->m_proj_parm);
444 #ifndef DOXYGEN_NO_DETAIL
449 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_mil_os, mil_os_ellipsoid)
450 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lee_os, lee_os_ellipsoid)
451 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_gs48, gs48_ellipsoid)
452 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_alsk, alsk_ellipsoid)
453 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_gs50, gs50_ellipsoid)
456 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(mil_os_entry, mil_os_ellipsoid)
457 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(lee_os_entry, lee_os_ellipsoid)
458 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(gs48_entry, gs48_ellipsoid)
459 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(alsk_entry, alsk_ellipsoid)
460 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(gs50_entry, gs50_ellipsoid)
462 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(mod_ster_init)
464 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(mil_os, mil_os_entry)
465 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(lee_os, lee_os_entry)
466 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(gs48, gs48_entry)
467 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(alsk, alsk_entry)
468 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(gs50, gs50_entry)
471 } // namespace detail
474 } // namespace projections
476 }} // namespace boost::geometry
478 #endif // BOOST_GEOMETRY_PROJECTIONS_MOD_STER_HPP