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 aeqd (Azimuthal Equidistant) projection.
23 // Author: Gerald Evenden
24 // Copyright (c) 1995, Gerald Evenden
26 // Permission is hereby granted, free of charge, to any person obtaining a
27 // copy of this software and associated documentation files (the "Software"),
28 // to deal in the Software without restriction, including without limitation
29 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
30 // and/or sell copies of the Software, and to permit persons to whom the
31 // Software is furnished to do so, subject to the following conditions:
33 // The above copyright notice and this permission notice shall be included
34 // in all copies or substantial portions of the Software.
36 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
37 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
38 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
39 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
40 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
41 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
42 // DEALINGS IN THE SOFTWARE.
44 #ifndef BOOST_GEOMETRY_PROJECTIONS_AEQD_HPP
45 #define BOOST_GEOMETRY_PROJECTIONS_AEQD_HPP
47 #include <boost/config.hpp>
49 #include <boost/geometry/formulas/vincenty_direct.hpp>
50 #include <boost/geometry/formulas/vincenty_inverse.hpp>
52 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
53 #include <boost/geometry/srs/projections/impl/base_static.hpp>
54 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
55 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
56 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
57 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
58 #include <boost/geometry/srs/projections/impl/projects.hpp>
60 #include <boost/geometry/util/math.hpp>
62 #include <boost/math/special_functions/hypot.hpp>
64 #include <boost/type_traits/is_same.hpp>
66 namespace boost { namespace geometry
71 #ifndef DOXYGEN_NO_DETAIL
72 namespace detail { namespace aeqd
75 static const double epsilon10 = 1.e-10;
76 static const double tolerance = 1.e-14;
99 template <typename T, typename Par, typename ProjParm>
100 inline void e_forward(T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y, Par const& par, ProjParm const& proj_parm)
102 T coslam, cosphi, sinphi, rho;
104 //T lam1, phi1, lam2, phi2;
106 coslam = cos(lp_lon);
107 cosphi = cos(lp_lat);
108 sinphi = sin(lp_lat);
109 switch (proj_parm.mode) {
114 xy_x = (rho = fabs(proj_parm.Mp - pj_mlfn(lp_lat, sinphi, cosphi, proj_parm.en))) *
120 if (fabs(lp_lon) < epsilon10 && fabs(lp_lat - par.phi0) < epsilon10) {
125 //phi1 = par.phi0; lam1 = par.lam0;
126 //phi2 = lp_lat; lam2 = lp_lon + par.lam0;
128 formula::result_inverse<T> const inv =
129 formula::vincenty_inverse
132 >::apply(par.lam0, par.phi0, lp_lon + par.lam0, lp_lat, srs::spheroid<T>(par.a, proj_parm.b));
133 //azi1 = inv.azimuth; s12 = inv.distance;
134 xy_x = inv.distance * sin(inv.azimuth) / par.a;
135 xy_y = inv.distance * cos(inv.azimuth) / par.a;
140 template <typename T, typename Par, typename ProjParm>
141 inline void e_inverse(T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat, Par const& par, ProjParm const& proj_parm)
145 if ((c = boost::math::hypot(xy_x, xy_y)) < epsilon10) {
150 if (proj_parm.mode == obliq || proj_parm.mode == equit) {
151 T const x2 = xy_x * par.a;
152 T const y2 = xy_y * par.a;
153 //T const lat1 = par.phi0;
154 //T const lon1 = par.lam0;
155 T const azi1 = atan2(x2, y2);
156 T const s12 = sqrt(x2 * x2 + y2 * y2);
157 formula::result_direct<T> const dir =
158 formula::vincenty_direct
161 >::apply(par.lam0, par.phi0, s12, azi1, srs::spheroid<T>(par.a, proj_parm.b));
166 lp_lat = pj_inv_mlfn(proj_parm.mode == n_pole ? proj_parm.Mp - c : proj_parm.Mp + c,
167 par.es, proj_parm.en);
168 lp_lon = atan2(xy_x, proj_parm.mode == n_pole ? -xy_y : xy_y);
172 template <typename T, typename Par, typename ProjParm>
173 inline void e_guam_fwd(T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y, Par const& par, ProjParm const& proj_parm)
177 cosphi = cos(lp_lat);
178 sinphi = sin(lp_lat);
179 t = 1. / sqrt(1. - par.es * sinphi * sinphi);
180 xy_x = lp_lon * cosphi * t;
181 xy_y = pj_mlfn(lp_lat, sinphi, cosphi, proj_parm.en) - proj_parm.M1 +
182 .5 * lp_lon * lp_lon * cosphi * sinphi * t;
185 template <typename T, typename Par, typename ProjParm>
186 inline void e_guam_inv(T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat, Par const& par, ProjParm const& proj_parm)
191 x2 = 0.5 * xy_x * xy_x;
193 for (i = 0; i < 3; ++i) {
194 t = par.e * sin(lp_lat);
195 lp_lat = pj_inv_mlfn(proj_parm.M1 + xy_y -
196 x2 * tan(lp_lat) * (t = sqrt(1. - t * t)), par.es, proj_parm.en);
198 lp_lon = xy_x * t / cos(lp_lat);
201 template <typename T, typename Par, typename ProjParm>
202 inline void s_forward(T const& lp_lon, T lp_lat, T& xy_x, T& xy_y, Par const& /*par*/, ProjParm const& proj_parm)
204 static const T half_pi = detail::half_pi<T>();
206 T coslam, cosphi, sinphi;
208 sinphi = sin(lp_lat);
209 cosphi = cos(lp_lat);
210 coslam = cos(lp_lon);
211 switch (proj_parm.mode) {
213 xy_y = cosphi * coslam;
216 xy_y = proj_parm.sinph0 * sinphi + proj_parm.cosph0 * cosphi * coslam;
218 if (fabs(fabs(xy_y) - 1.) < tolerance)
220 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
226 xy_x = xy_y * cosphi * sin(lp_lon);
227 xy_y *= (proj_parm.mode == equit) ? sinphi :
228 proj_parm.cosph0 * sinphi - proj_parm.sinph0 * cosphi * coslam;
236 if (fabs(lp_lat - half_pi) < epsilon10)
237 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
238 xy_x = (xy_y = (half_pi + lp_lat)) * sin(lp_lon);
244 template <typename T, typename Par, typename ProjParm>
245 inline void s_inverse(T xy_x, T xy_y, T& lp_lon, T& lp_lat, Par const& par, ProjParm const& proj_parm)
247 static const T pi = detail::pi<T>();
248 static const T half_pi = detail::half_pi<T>();
252 if ((c_rh = boost::math::hypot(xy_x, xy_y)) > pi) {
253 if (c_rh - epsilon10 > pi)
254 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
256 } else if (c_rh < epsilon10) {
261 if (proj_parm.mode == obliq || proj_parm.mode == equit) {
264 if (proj_parm.mode == equit) {
265 lp_lat = aasin(xy_y * sinc / c_rh);
269 lp_lat = aasin(cosc * proj_parm.sinph0 + xy_y * sinc * proj_parm.cosph0 /
271 xy_y = (cosc - proj_parm.sinph0 * sin(lp_lat)) * c_rh;
272 xy_x *= sinc * proj_parm.cosph0;
274 lp_lon = xy_y == 0. ? 0. : atan2(xy_x, xy_y);
275 } else if (proj_parm.mode == n_pole) {
276 lp_lat = half_pi - c_rh;
277 lp_lon = atan2(xy_x, -xy_y);
279 lp_lat = c_rh - half_pi;
280 lp_lon = atan2(xy_x, xy_y);
284 // Azimuthal Equidistant
285 template <typename Params, typename Parameters, typename T>
286 inline void setup_aeqd(Params const& params, Parameters& par, par_aeqd<T>& proj_parm, bool is_sphere, bool is_guam)
288 static const T half_pi = detail::half_pi<T>();
290 par.phi0 = pj_get_param_r<T, srs::spar::lat_0>(params, "lat_0", srs::dpar::lat_0);
291 if (fabs(fabs(par.phi0) - half_pi) < epsilon10) {
292 proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
293 proj_parm.sinph0 = par.phi0 < 0. ? -1. : 1.;
294 proj_parm.cosph0 = 0.;
295 } else if (fabs(par.phi0) < epsilon10) {
296 proj_parm.mode = equit;
297 proj_parm.sinph0 = 0.;
298 proj_parm.cosph0 = 1.;
300 proj_parm.mode = obliq;
301 proj_parm.sinph0 = sin(par.phi0);
302 proj_parm.cosph0 = cos(par.phi0);
307 proj_parm.en = pj_enfn<T>(par.es);
309 proj_parm.M1 = pj_mlfn(par.phi0, proj_parm.sinph0, proj_parm.cosph0, proj_parm.en);
311 switch (proj_parm.mode) {
313 proj_parm.Mp = pj_mlfn<T>(half_pi, 1., 0., proj_parm.en);
316 proj_parm.Mp = pj_mlfn<T>(-half_pi, -1., 0., proj_parm.en);
320 //proj_parm.N1 = 1. / sqrt(1. - par.es * proj_parm.sinph0 * proj_parm.sinph0);
321 //proj_parm.G = proj_parm.sinph0 * (proj_parm.He = par.e / sqrt(par.one_es));
322 //proj_parm.He *= proj_parm.cosph0;
325 // Boost.Geometry specific, in proj4 geodesic is initialized at the beginning
326 proj_parm.b = math::sqrt(math::sqr(par.a) * (1. - par.es));
331 template <typename T, typename Parameters>
334 par_aeqd<T> m_proj_parm;
336 // FORWARD(e_forward) elliptical
337 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
338 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
340 e_forward(lp_lon, lp_lat, xy_x, xy_y, par, this->m_proj_parm);
343 // INVERSE(e_inverse) elliptical
344 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
345 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
347 e_inverse(xy_x, xy_y, lp_lon, lp_lat, par, this->m_proj_parm);
350 static inline std::string get_name()
357 template <typename T, typename Parameters>
358 struct base_aeqd_e_guam
360 par_aeqd<T> m_proj_parm;
362 // FORWARD(e_guam_fwd) Guam elliptical
363 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
364 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
366 e_guam_fwd(lp_lon, lp_lat, xy_x, xy_y, par, this->m_proj_parm);
369 // INVERSE(e_guam_inv) Guam elliptical
370 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
371 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
373 e_guam_inv(xy_x, xy_y, lp_lon, lp_lat, par, this->m_proj_parm);
376 static inline std::string get_name()
378 return "aeqd_e_guam";
383 template <typename T, typename Parameters>
386 par_aeqd<T> m_proj_parm;
388 // FORWARD(s_forward) spherical
389 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
390 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
392 s_forward(lp_lon, lp_lat, xy_x, xy_y, par, this->m_proj_parm);
395 // INVERSE(s_inverse) spherical
396 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
397 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
399 s_inverse(xy_x, xy_y, lp_lon, lp_lat, par, this->m_proj_parm);
402 static inline std::string get_name()
409 }} // namespace detail::aeqd
413 \brief Azimuthal Equidistant projection
415 \tparam Geographic latlong point type
416 \tparam Cartesian xy point type
417 \tparam Parameters parameter type
418 \par Projection characteristics
422 \par Projection parameters
423 - lat_0: Latitude of origin (degrees)
426 \image html ex_aeqd.gif
428 template <typename T, typename Parameters>
429 struct aeqd_e : public detail::aeqd::base_aeqd_e<T, Parameters>
431 template <typename Params>
432 inline aeqd_e(Params const& params, Parameters & par)
434 detail::aeqd::setup_aeqd(params, par, this->m_proj_parm, false, false);
439 \brief Azimuthal Equidistant projection
441 \tparam Geographic latlong point type
442 \tparam Cartesian xy point type
443 \tparam Parameters parameter type
444 \par Projection characteristics
448 \par Projection parameters
449 - lat_0: Latitude of origin (degrees)
452 \image html ex_aeqd.gif
454 template <typename T, typename Parameters>
455 struct aeqd_e_guam : public detail::aeqd::base_aeqd_e_guam<T, Parameters>
457 template <typename Params>
458 inline aeqd_e_guam(Params const& params, Parameters & par)
460 detail::aeqd::setup_aeqd(params, par, this->m_proj_parm, false, true);
465 \brief Azimuthal Equidistant projection
467 \tparam Geographic latlong point type
468 \tparam Cartesian xy point type
469 \tparam Parameters parameter type
470 \par Projection characteristics
474 \par Projection parameters
475 - lat_0: Latitude of origin (degrees)
478 \image html ex_aeqd.gif
480 template <typename T, typename Parameters>
481 struct aeqd_s : public detail::aeqd::base_aeqd_s<T, Parameters>
483 template <typename Params>
484 inline aeqd_s(Params const& params, Parameters & par)
486 detail::aeqd::setup_aeqd(params, par, this->m_proj_parm, true, false);
490 #ifndef DOXYGEN_NO_DETAIL
495 template <typename BGP, typename CT, typename P>
496 struct static_projection_type<srs::spar::proj_aeqd, srs_sphere_tag, BGP, CT, P>
498 typedef static_wrapper_fi<aeqd_s<CT, P>, P> type;
500 template <typename BGP, typename CT, typename P>
501 struct static_projection_type<srs::spar::proj_aeqd, srs_spheroid_tag, BGP, CT, P>
503 typedef static_wrapper_fi
505 typename boost::mpl::if_c
509 typename geometry::tuples::find_if
512 //srs::par4::detail::is_guam
513 srs::spar::detail::is_param<srs::spar::guam>::pred
524 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_BEGIN(aeqd_entry)
526 bool const guam = pj_get_param_b<srs::spar::guam>(params, "guam", srs::dpar::guam);
528 if (parameters.es && ! guam)
529 return new dynamic_wrapper_fi<aeqd_e<T, Parameters>, T, Parameters>(params, parameters);
530 else if (parameters.es && guam)
531 return new dynamic_wrapper_fi<aeqd_e_guam<T, Parameters>, T, Parameters>(params, parameters);
533 return new dynamic_wrapper_fi<aeqd_s<T, Parameters>, T, Parameters>(params, parameters);
535 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_END
537 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(aeqd_init)
539 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(aeqd, aeqd_entry)
542 } // namespace detail
545 } // namespace projections
547 }} // namespace boost::geometry
549 #endif // BOOST_GEOMETRY_PROJECTIONS_AEQD_HPP