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_STERE_HPP
41 #define BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
43 #include <boost/config.hpp>
44 #include <boost/geometry/util/math.hpp>
45 #include <boost/math/special_functions/hypot.hpp>
47 #include <boost/geometry/srs/projections/impl/base_static.hpp>
48 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
49 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
50 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
51 #include <boost/geometry/srs/projections/impl/pj_tsfn.hpp>
52 #include <boost/geometry/srs/projections/impl/projects.hpp>
54 namespace boost { namespace geometry
59 #ifndef DOXYGEN_NO_DETAIL
60 namespace detail { namespace stere
62 static const double epsilon10 = 1.e-10;
63 static const double tolerance = 1.e-8;
64 static const int n_iter = 8;
65 static const double conv_tolerance = 1.e-10;
85 inline T ssfn_(T const& phit, T sinphi, T const& eccen)
87 static const T half_pi = detail::half_pi<T>();
90 return (tan (.5 * (half_pi + phit)) *
91 math::pow((T(1) - sinphi) / (T(1) + sinphi), T(0.5) * eccen));
94 template <typename T, typename Parameters>
95 struct base_stere_ellipsoid
97 par_stere<T> m_proj_parm;
99 // FORWARD(e_forward) ellipsoid
100 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
101 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
103 static const T half_pi = detail::half_pi<T>();
105 T coslam, sinlam, sinX=0.0, cosX=0.0, X, A = 0.0, sinphi;
107 coslam = cos(lp_lon);
108 sinlam = sin(lp_lon);
109 sinphi = sin(lp_lat);
110 if (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit) {
111 sinX = sin(X = 2. * atan(ssfn_(lp_lat, sinphi, par.e)) - half_pi);
114 switch (this->m_proj_parm.mode) {
116 A = this->m_proj_parm.akm1 / (this->m_proj_parm.cosX1 * (1. + this->m_proj_parm.sinX1 * sinX +
117 this->m_proj_parm.cosX1 * cosX * coslam));
118 xy_y = A * (this->m_proj_parm.cosX1 * sinX - this->m_proj_parm.sinX1 * cosX * coslam);
119 goto xmul; /* but why not just xy.x = A * cosX; break; ? */
122 // TODO: calculate denominator once
123 /* avoid zero division */
124 if (1. + cosX * coslam == 0.0) {
127 A = this->m_proj_parm.akm1 / (1. + cosX * coslam);
140 xy_x = this->m_proj_parm.akm1 * pj_tsfn(lp_lat, sinphi, par.e);
141 xy_y = - xy_x * coslam;
145 xy_x = xy_x * sinlam;
148 // INVERSE(e_inverse) ellipsoid
149 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
150 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
152 static const T half_pi = detail::half_pi<T>();
154 T cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;
157 rho = boost::math::hypot(xy_x, xy_y);
158 switch (this->m_proj_parm.mode) {
161 cosphi = cos( tp = 2. * atan2(rho * this->m_proj_parm.cosX1 , this->m_proj_parm.akm1) );
164 phi_l = asin(cosphi * this->m_proj_parm.sinX1);
166 phi_l = asin(cosphi * this->m_proj_parm.sinX1 + (xy_y * sinphi * this->m_proj_parm.cosX1 / rho));
168 tp = tan(.5 * (half_pi + phi_l));
170 xy_y = rho * this->m_proj_parm.cosX1 * cosphi - xy_y * this->m_proj_parm.sinX1* sinphi;
178 phi_l = half_pi - 2. * atan(tp = - rho / this->m_proj_parm.akm1);
183 for (i = n_iter; i--; phi_l = lp_lat) {
184 sinphi = par.e * sin(phi_l);
185 lp_lat = T(2) * atan(tp * math::pow((T(1)+sinphi)/(T(1)-sinphi), halfe)) - halfpi;
186 if (fabs(phi_l - lp_lat) < conv_tolerance) {
187 if (this->m_proj_parm.mode == s_pole)
189 lp_lon = (xy_x == 0. && xy_y == 0.) ? 0. : atan2(xy_x, xy_y);
193 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
196 static inline std::string get_name()
198 return "stere_ellipsoid";
203 template <typename T, typename Parameters>
204 struct base_stere_spheroid
206 par_stere<T> m_proj_parm;
208 // FORWARD(s_forward) spheroid
209 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
210 inline void fwd(Parameters const& , T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
212 static const T fourth_pi = detail::fourth_pi<T>();
213 static const T half_pi = detail::half_pi<T>();
215 T sinphi, cosphi, coslam, sinlam;
217 sinphi = sin(lp_lat);
218 cosphi = cos(lp_lat);
219 coslam = cos(lp_lon);
220 sinlam = sin(lp_lon);
221 switch (this->m_proj_parm.mode) {
223 xy_y = 1. + cosphi * coslam;
226 xy_y = 1. + this->m_proj_parm.sinX1 * sinphi + this->m_proj_parm.cosX1 * cosphi * coslam;
228 if (xy_y <= epsilon10) {
229 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
231 xy_x = (xy_y = this->m_proj_parm.akm1 / xy_y) * cosphi * sinlam;
232 xy_y *= (this->m_proj_parm.mode == equit) ? sinphi :
233 this->m_proj_parm.cosX1 * sinphi - this->m_proj_parm.sinX1 * cosphi * coslam;
240 if (fabs(lp_lat - half_pi) < tolerance) {
241 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
243 xy_x = sinlam * ( xy_y = this->m_proj_parm.akm1 * tan(fourth_pi + .5 * lp_lat) );
249 // INVERSE(s_inverse) spheroid
250 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
251 inline void inv(Parameters const& par, T const& xy_x, T xy_y, T& lp_lon, T& lp_lat) const
255 sinc = sin(c = 2. * atan((rh = boost::math::hypot(xy_x, xy_y)) / this->m_proj_parm.akm1));
259 switch (this->m_proj_parm.mode) {
261 if (fabs(rh) <= epsilon10)
264 lp_lat = asin(xy_y * sinc / rh);
265 if (cosc != 0. || xy_x != 0.)
266 lp_lon = atan2(xy_x * sinc, cosc * rh);
269 if (fabs(rh) <= epsilon10)
272 lp_lat = asin(cosc * this->m_proj_parm.sinX1 + xy_y * sinc * this->m_proj_parm.cosX1 / rh);
273 if ((c = cosc - this->m_proj_parm.sinX1 * sin(lp_lat)) != 0. || xy_x != 0.)
274 lp_lon = atan2(xy_x * sinc * this->m_proj_parm.cosX1, c * rh);
280 if (fabs(rh) <= epsilon10)
283 lp_lat = asin(this->m_proj_parm.mode == s_pole ? - cosc : cosc);
284 lp_lon = (xy_x == 0. && xy_y == 0.) ? 0. : atan2(xy_x, xy_y);
289 static inline std::string get_name()
291 return "stere_spheroid";
296 template <typename Parameters, typename T>
297 inline void setup(Parameters const& par, par_stere<T>& proj_parm) /* general initialization */
299 static const T fourth_pi = detail::fourth_pi<T>();
300 static const T half_pi = detail::half_pi<T>();
304 if (fabs((t = fabs(par.phi0)) - half_pi) < epsilon10)
305 proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
307 proj_parm.mode = t > epsilon10 ? obliq : equit;
308 proj_parm.phits = fabs(proj_parm.phits);
313 switch (proj_parm.mode) {
316 if (fabs(proj_parm.phits - half_pi) < epsilon10)
317 proj_parm.akm1 = 2. * par.k0 /
318 sqrt(math::pow(T(1)+par.e,T(1)+par.e)*math::pow(T(1)-par.e,T(1)-par.e));
320 proj_parm.akm1 = cos(proj_parm.phits) /
321 pj_tsfn(proj_parm.phits, t = sin(proj_parm.phits), par.e);
323 proj_parm.akm1 /= sqrt(1. - t * t);
329 X = 2. * atan(ssfn_(par.phi0, t, par.e)) - half_pi;
331 proj_parm.akm1 = 2. * par.k0 * cos(par.phi0) / sqrt(1. - t * t);
332 proj_parm.sinX1 = sin(X);
333 proj_parm.cosX1 = cos(X);
337 switch (proj_parm.mode) {
339 proj_parm.sinX1 = sin(par.phi0);
340 proj_parm.cosX1 = cos(par.phi0);
343 proj_parm.akm1 = 2. * par.k0;
347 proj_parm.akm1 = fabs(proj_parm.phits - half_pi) >= epsilon10 ?
348 cos(proj_parm.phits) / tan(fourth_pi - .5 * proj_parm.phits) :
357 template <typename Params, typename Parameters, typename T>
358 inline void setup_stere(Params const& params, Parameters const& par, par_stere<T>& proj_parm)
360 static const T half_pi = detail::half_pi<T>();
362 if (! pj_param_r<srs::spar::lat_ts>(params, "lat_ts", srs::dpar::lat_ts, proj_parm.phits))
363 proj_parm.phits = half_pi;
365 setup(par, proj_parm);
368 // Universal Polar Stereographic
369 template <typename Params, typename Parameters, typename T>
370 inline void setup_ups(Params const& params, Parameters& par, par_stere<T>& proj_parm)
372 static const T half_pi = detail::half_pi<T>();
374 /* International Ellipsoid */
375 par.phi0 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? -half_pi: half_pi;
377 BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
382 proj_parm.phits = half_pi;
385 setup(par, proj_parm);
388 }} // namespace detail::stere
392 \brief Stereographic projection
394 \tparam Geographic latlong point type
395 \tparam Cartesian xy point type
396 \tparam Parameters parameter type
397 \par Projection characteristics
401 \par Projection parameters
402 - lat_ts: Latitude of true scale (degrees)
404 \image html ex_stere.gif
406 template <typename T, typename Parameters>
407 struct stere_ellipsoid : public detail::stere::base_stere_ellipsoid<T, Parameters>
409 template <typename Params>
410 inline stere_ellipsoid(Params const& params, Parameters const& par)
412 detail::stere::setup_stere(params, par, this->m_proj_parm);
417 \brief Stereographic projection
419 \tparam Geographic latlong point type
420 \tparam Cartesian xy point type
421 \tparam Parameters parameter type
422 \par Projection characteristics
426 \par Projection parameters
427 - lat_ts: Latitude of true scale (degrees)
429 \image html ex_stere.gif
431 template <typename T, typename Parameters>
432 struct stere_spheroid : public detail::stere::base_stere_spheroid<T, Parameters>
434 template <typename Params>
435 inline stere_spheroid(Params const& params, Parameters const& par)
437 detail::stere::setup_stere(params, par, this->m_proj_parm);
442 \brief Universal Polar Stereographic projection
444 \tparam Geographic latlong point type
445 \tparam Cartesian xy point type
446 \tparam Parameters parameter type
447 \par Projection characteristics
451 \par Projection parameters
452 - south: Denotes southern hemisphere UTM zone (boolean)
454 \image html ex_ups.gif
456 template <typename T, typename Parameters>
457 struct ups_ellipsoid : public detail::stere::base_stere_ellipsoid<T, Parameters>
459 template <typename Params>
460 inline ups_ellipsoid(Params const& params, Parameters & par)
462 detail::stere::setup_ups(params, par, this->m_proj_parm);
467 \brief Universal Polar Stereographic projection
469 \tparam Geographic latlong point type
470 \tparam Cartesian xy point type
471 \tparam Parameters parameter type
472 \par Projection characteristics
476 \par Projection parameters
477 - south: Denotes southern hemisphere UTM zone (boolean)
479 \image html ex_ups.gif
481 template <typename T, typename Parameters>
482 struct ups_spheroid : public detail::stere::base_stere_spheroid<T, Parameters>
484 template <typename Params>
485 inline ups_spheroid(Params const& params, Parameters & par)
487 detail::stere::setup_ups(params, par, this->m_proj_parm);
491 #ifndef DOXYGEN_NO_DETAIL
496 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_stere, stere_spheroid, stere_ellipsoid)
497 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_ups, ups_spheroid, ups_ellipsoid)
500 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(stere_entry, stere_spheroid, stere_ellipsoid)
501 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(ups_entry, ups_spheroid, ups_ellipsoid)
503 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(stere_init)
505 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(stere, stere_entry)
506 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(ups, ups_entry)
509 } // namespace detail
512 } // namespace projections
514 }} // namespace boost::geometry
516 #endif // BOOST_GEOMETRY_PROJECTIONS_STERE_HPP