boost/geometry/srs/projections/proj/tpeqd.hpp

   1 // Boost.Geometry - gis-projections (based on PROJ4)
   2
   3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
   4
   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.
   8
   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)
  12
  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
  17
  18 // Last updated version of proj: 5.0.0
  19
  20 // Original copyright notice:
  21
  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:
  28
  29 // The above copyright notice and this permission notice shall be included
  30 // in all copies or substantial portions of the Software.
  31
  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.
  39
  40 #ifndef BOOST_GEOMETRY_PROJECTIONS_TPEQD_HPP
  41 #define BOOST_GEOMETRY_PROJECTIONS_TPEQD_HPP
  42
  43 #include <boost/geometry/util/math.hpp>
  44 #include <boost/math/special_functions/hypot.hpp>
  45
  46 #include <boost/geometry/srs/projections/impl/aasincos.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/projects.hpp>
  52
  53 namespace boost { namespace geometry
  54 {
  55
  56 namespace projections
  57 {
  58     #ifndef DOXYGEN_NO_DETAIL
  59     namespace detail { namespace tpeqd
  60     {
  61             template <typename T>
  62             struct par_tpeqd
  63             {
  64                 T cp1, sp1, cp2, sp2, ccs, cs, sc, r2z0, z02, dlam2;
  65                 T hz0, thz0, rhshz0, ca, sa, lp, lamc;
  66             };
  67
  68             template <typename T, typename Parameters>
  69             struct base_tpeqd_spheroid
  70             {
  71                 par_tpeqd<T> m_proj_parm;
  72
  73                 // FORWARD(s_forward)  sphere
  74                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  75                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  76                 {
  77                     T t, z1, z2, dl1, dl2, sp, cp;
  78
  79                     sp = sin(lp_lat);
  80                     cp = cos(lp_lat);
  81                     z1 = aacos(this->m_proj_parm.sp1 * sp + this->m_proj_parm.cp1 * cp * cos(dl1 = lp_lon + this->m_proj_parm.dlam2));
  82                     z2 = aacos(this->m_proj_parm.sp2 * sp + this->m_proj_parm.cp2 * cp * cos(dl2 = lp_lon - this->m_proj_parm.dlam2));
  83                     z1 *= z1;
  84                     z2 *= z2;
  85
  86                     xy_x = this->m_proj_parm.r2z0 * (t = z1 - z2);
  87                     t = this->m_proj_parm.z02 - t;
  88                     xy_y = this->m_proj_parm.r2z0 * asqrt(4. * this->m_proj_parm.z02 * z2 - t * t);
  89                     if ((this->m_proj_parm.ccs * sp - cp * (this->m_proj_parm.cs * sin(dl1) - this->m_proj_parm.sc * sin(dl2))) < 0.)
  90                         xy_y = -xy_y;
  91                 }
  92
  93                 // INVERSE(s_inverse)  sphere
  94                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  95                 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
  96                 {
  97                     T cz1, cz2, s, d, cp, sp;
  98
  99                     cz1 = cos(boost::math::hypot(xy_y, xy_x + this->m_proj_parm.hz0));
 100                     cz2 = cos(boost::math::hypot(xy_y, xy_x - this->m_proj_parm.hz0));
 101                     s = cz1 + cz2;
 102                     d = cz1 - cz2;
 103                     lp_lon = - atan2(d, (s * this->m_proj_parm.thz0));
 104                     lp_lat = aacos(boost::math::hypot(this->m_proj_parm.thz0 * s, d) * this->m_proj_parm.rhshz0);
 105                     if ( xy_y < 0. )
 106                         lp_lat = - lp_lat;
 107                     /* lam--phi now in system relative to P1--P2 base equator */
 108                     sp = sin(lp_lat);
 109                     cp = cos(lp_lat);
 110                     lp_lat = aasin(this->m_proj_parm.sa * sp + this->m_proj_parm.ca * cp * (s = cos(lp_lon -= this->m_proj_parm.lp)));
 111                     lp_lon = atan2(cp * sin(lp_lon), this->m_proj_parm.sa * cp * s - this->m_proj_parm.ca * sp) + this->m_proj_parm.lamc;
 112                 }
 113
 114                 static inline std::string get_name()
 115                 {
 116                     return "tpeqd_spheroid";
 117                 }
 118
 119             };
 120
 121             // Two Point Equidistant
 122             template <typename Params, typename Parameters, typename T>
 123             inline void setup_tpeqd(Params const& params, Parameters& par, par_tpeqd<T>& proj_parm)
 124             {
 125                 T lam_1, lam_2, phi_1, phi_2, A12, pp;
 126
 127                 /* get control point locations */
 128                 phi_1 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
 129                 lam_1 = pj_get_param_r<T, srs::spar::lon_1>(params, "lon_1", srs::dpar::lon_1);
 130                 phi_2 = pj_get_param_r<T, srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2);
 131                 lam_2 = pj_get_param_r<T, srs::spar::lon_2>(params, "lon_2", srs::dpar::lon_2);
 132
 133                 if (phi_1 == phi_2 && lam_1 == lam_2)
 134                     BOOST_THROW_EXCEPTION( projection_exception(error_control_point_no_dist) );
 135
 136                 par.lam0 = adjlon(0.5 * (lam_1 + lam_2));
 137                 proj_parm.dlam2 = adjlon(lam_2 - lam_1);
 138
 139                 proj_parm.cp1 = cos(phi_1);
 140                 proj_parm.cp2 = cos(phi_2);
 141                 proj_parm.sp1 = sin(phi_1);
 142                 proj_parm.sp2 = sin(phi_2);
 143                 proj_parm.cs = proj_parm.cp1 * proj_parm.sp2;
 144                 proj_parm.sc = proj_parm.sp1 * proj_parm.cp2;
 145                 proj_parm.ccs = proj_parm.cp1 * proj_parm.cp2 * sin(proj_parm.dlam2);
 146                 proj_parm.z02 = aacos(proj_parm.sp1 * proj_parm.sp2 + proj_parm.cp1 * proj_parm.cp2 * cos(proj_parm.dlam2));
 147                 proj_parm.hz0 = .5 * proj_parm.z02;
 148                 A12 = atan2(proj_parm.cp2 * sin(proj_parm.dlam2),
 149                     proj_parm.cp1 * proj_parm.sp2 - proj_parm.sp1 * proj_parm.cp2 * cos(proj_parm.dlam2));
 150                 proj_parm.ca = cos(pp = aasin(proj_parm.cp1 * sin(A12)));
 151                 proj_parm.sa = sin(pp);
 152                 proj_parm.lp = adjlon(atan2(proj_parm.cp1 * cos(A12), proj_parm.sp1) - proj_parm.hz0);
 153                 proj_parm.dlam2 *= .5;
 154                 proj_parm.lamc = geometry::math::half_pi<T>() - atan2(sin(A12) * proj_parm.sp1, cos(A12)) - proj_parm.dlam2;
 155                 proj_parm.thz0 = tan(proj_parm.hz0);
 156                 proj_parm.rhshz0 = .5 / sin(proj_parm.hz0);
 157                 proj_parm.r2z0 = 0.5 / proj_parm.z02;
 158                 proj_parm.z02 *= proj_parm.z02;
 159
 160                 par.es = 0.;
 161             }
 162
 163     }} // namespace detail::tpeqd
 164     #endif // doxygen
 165
 166     /*!
 167         \brief Two Point Equidistant projection
 168         \ingroup projections
 169         \tparam Geographic latlong point type
 170         \tparam Cartesian xy point type
 171         \tparam Parameters parameter type
 172         \par Projection characteristics
 173          - Miscellaneous
 174          - Spheroid
 175         \par Projection parameters
 176          - lat_1: Latitude of first standard parallel (degrees)
 177          - lon_1 (degrees)
 178          - lat_2: Latitude of second standard parallel (degrees)
 179          - lon_2 (degrees)
 180         \par Example
 181         \image html ex_tpeqd.gif
 182     */
 183     template <typename T, typename Parameters>
 184     struct tpeqd_spheroid : public detail::tpeqd::base_tpeqd_spheroid<T, Parameters>
 185     {
 186         template <typename Params>
 187         inline tpeqd_spheroid(Params const& params, Parameters & par)
 188         {
 189             detail::tpeqd::setup_tpeqd(params, par, this->m_proj_parm);
 190         }
 191     };
 192
 193     #ifndef DOXYGEN_NO_DETAIL
 194     namespace detail
 195     {
 196
 197         // Static projection
 198         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_tpeqd, tpeqd_spheroid)
 199
 200         // Factory entry(s)
 201         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(tpeqd_entry, tpeqd_spheroid)
 202
 203         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(tpeqd_init)
 204         {
 205             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(tpeqd, tpeqd_entry)
 206         }
 207
 208     } // namespace detail
 209     #endif // doxygen
 210
 211 } // namespace projections
 212
 213 }} // namespace boost::geometry
 214
 215 #endif // BOOST_GEOMETRY_PROJECTIONS_TPEQD_HPP
 216