projections/proj/tpeqd.hpp

0001 // Boost.Geometry - gis-projections (based on PROJ4)
0002
0003 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2017, 2018, 2019.
0006 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0014 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0015 // PROJ4 is maintained by Frank Warmerdam
0016 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0017
0018 // Last updated version of proj: 5.0.0
0019
0020 // Original copyright notice:
0021
0022 // Permission is hereby granted, free of charge, to any person obtaining a
0023 // copy of this software and associated documentation files (the "Software"),
0024 // to deal in the Software without restriction, including without limitation
0025 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0026 // and/or sell copies of the Software, and to permit persons to whom the
0027 // Software is furnished to do so, subject to the following conditions:
0028
0029 // The above copyright notice and this permission notice shall be included
0030 // in all copies or substantial portions of the Software.
0031
0032 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0033 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0034 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0035 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0036 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0037 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0038 // DEALINGS IN THE SOFTWARE.
0039
0040 #ifndef BOOST_GEOMETRY_PROJECTIONS_TPEQD_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_TPEQD_HPP
0042
0043 #include <boost/geometry/util/math.hpp>
0044 #include <boost/math/special_functions/hypot.hpp>
0045
0046 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
0047 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0048 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0049 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0050 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0051 #include <boost/geometry/srs/projections/impl/projects.hpp>
0052
0053 namespace boost { namespace geometry
0054 {
0055
0056 namespace projections
0057 {
0058     #ifndef DOXYGEN_NO_DETAIL
0059     namespace detail { namespace tpeqd
0060     {
0061             template <typename T>
0062             struct par_tpeqd
0063             {
0064                 T cp1, sp1, cp2, sp2, ccs, cs, sc, r2z0, z02, dlam2;
0065                 T hz0, thz0, rhshz0, ca, sa, lp, lamc;
0066             };
0067
0068             template <typename T, typename Parameters>
0069             struct base_tpeqd_spheroid
0070             {
0071                 par_tpeqd<T> m_proj_parm;
0072
0073                 // FORWARD(s_forward)  sphere
0074                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0075                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0076                 {
0077                     T t, z1, z2, dl1, dl2, sp, cp;
0078
0079                     sp = sin(lp_lat);
0080                     cp = cos(lp_lat);
0081                     z1 = aacos(this->m_proj_parm.sp1 * sp + this->m_proj_parm.cp1 * cp * cos(dl1 = lp_lon + this->m_proj_parm.dlam2));
0082                     z2 = aacos(this->m_proj_parm.sp2 * sp + this->m_proj_parm.cp2 * cp * cos(dl2 = lp_lon - this->m_proj_parm.dlam2));
0083                     z1 *= z1;
0084                     z2 *= z2;
0085
0086                     xy_x = this->m_proj_parm.r2z0 * (t = z1 - z2);
0087                     t = this->m_proj_parm.z02 - t;
0088                     xy_y = this->m_proj_parm.r2z0 * asqrt(4. * this->m_proj_parm.z02 * z2 - t * t);
0089                     if ((this->m_proj_parm.ccs * sp - cp * (this->m_proj_parm.cs * sin(dl1) - this->m_proj_parm.sc * sin(dl2))) < 0.)
0090                         xy_y = -xy_y;
0091                 }
0092
0093                 // INVERSE(s_inverse)  sphere
0094                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0095                 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0096                 {
0097                     T cz1, cz2, s, d, cp, sp;
0098
0099                     cz1 = cos(boost::math::hypot(xy_y, xy_x + this->m_proj_parm.hz0));
0100                     cz2 = cos(boost::math::hypot(xy_y, xy_x - this->m_proj_parm.hz0));
0101                     s = cz1 + cz2;
0102                     d = cz1 - cz2;
0103                     lp_lon = - atan2(d, (s * this->m_proj_parm.thz0));
0104                     lp_lat = aacos(boost::math::hypot(this->m_proj_parm.thz0 * s, d) * this->m_proj_parm.rhshz0);
0105                     if ( xy_y < 0. )
0106                         lp_lat = - lp_lat;
0107                     /* lam--phi now in system relative to P1--P2 base equator */
0108                     sp = sin(lp_lat);
0109                     cp = cos(lp_lat);
0110                     lp_lat = aasin(this->m_proj_parm.sa * sp + this->m_proj_parm.ca * cp * (s = cos(lp_lon -= this->m_proj_parm.lp)));
0111                     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;
0112                 }
0113
0114                 static inline std::string get_name()
0115                 {
0116                     return "tpeqd_spheroid";
0117                 }
0118
0119             };
0120
0121             // Two Point Equidistant
0122             template <typename Params, typename Parameters, typename T>
0123             inline void setup_tpeqd(Params const& params, Parameters& par, par_tpeqd<T>& proj_parm)
0124             {
0125                 T lam_1, lam_2, phi_1, phi_2, A12, pp;
0126
0127                 /* get control point locations */
0128                 phi_1 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
0129                 lam_1 = pj_get_param_r<T, srs::spar::lon_1>(params, "lon_1", srs::dpar::lon_1);
0130                 phi_2 = pj_get_param_r<T, srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2);
0131                 lam_2 = pj_get_param_r<T, srs::spar::lon_2>(params, "lon_2", srs::dpar::lon_2);
0132
0133                 if (phi_1 == phi_2 && lam_1 == lam_2)
0134                     BOOST_THROW_EXCEPTION( projection_exception(error_control_point_no_dist) );
0135
0136                 par.lam0 = adjlon(0.5 * (lam_1 + lam_2));
0137                 proj_parm.dlam2 = adjlon(lam_2 - lam_1);
0138
0139                 proj_parm.cp1 = cos(phi_1);
0140                 proj_parm.cp2 = cos(phi_2);
0141                 proj_parm.sp1 = sin(phi_1);
0142                 proj_parm.sp2 = sin(phi_2);
0143                 proj_parm.cs = proj_parm.cp1 * proj_parm.sp2;
0144                 proj_parm.sc = proj_parm.sp1 * proj_parm.cp2;
0145                 proj_parm.ccs = proj_parm.cp1 * proj_parm.cp2 * sin(proj_parm.dlam2);
0146                 proj_parm.z02 = aacos(proj_parm.sp1 * proj_parm.sp2 + proj_parm.cp1 * proj_parm.cp2 * cos(proj_parm.dlam2));
0147                 proj_parm.hz0 = .5 * proj_parm.z02;
0148                 A12 = atan2(proj_parm.cp2 * sin(proj_parm.dlam2),
0149                     proj_parm.cp1 * proj_parm.sp2 - proj_parm.sp1 * proj_parm.cp2 * cos(proj_parm.dlam2));
0150                 proj_parm.ca = cos(pp = aasin(proj_parm.cp1 * sin(A12)));
0151                 proj_parm.sa = sin(pp);
0152                 proj_parm.lp = adjlon(atan2(proj_parm.cp1 * cos(A12), proj_parm.sp1) - proj_parm.hz0);
0153                 proj_parm.dlam2 *= .5;
0154                 proj_parm.lamc = geometry::math::half_pi<T>() - atan2(sin(A12) * proj_parm.sp1, cos(A12)) - proj_parm.dlam2;
0155                 proj_parm.thz0 = tan(proj_parm.hz0);
0156                 proj_parm.rhshz0 = .5 / sin(proj_parm.hz0);
0157                 proj_parm.r2z0 = 0.5 / proj_parm.z02;
0158                 proj_parm.z02 *= proj_parm.z02;
0159
0160                 par.es = 0.;
0161             }
0162
0163     }} // namespace detail::tpeqd
0164     #endif // doxygen
0165
0166     /*!
0167         \brief Two Point Equidistant projection
0168         \ingroup projections
0169         \tparam Geographic latlong point type
0170         \tparam Cartesian xy point type
0171         \tparam Parameters parameter type
0172         \par Projection characteristics
0173          - Miscellaneous
0174          - Spheroid
0175         \par Projection parameters
0176          - lat_1: Latitude of first standard parallel (degrees)
0177          - lon_1 (degrees)
0178          - lat_2: Latitude of second standard parallel (degrees)
0179          - lon_2 (degrees)
0180         \par Example
0181         \image html ex_tpeqd.gif
0182     */
0183     template <typename T, typename Parameters>
0184     struct tpeqd_spheroid : public detail::tpeqd::base_tpeqd_spheroid<T, Parameters>
0185     {
0186         template <typename Params>
0187         inline tpeqd_spheroid(Params const& params, Parameters & par)
0188         {
0189             detail::tpeqd::setup_tpeqd(params, par, this->m_proj_parm);
0190         }
0191     };
0192
0193     #ifndef DOXYGEN_NO_DETAIL
0194     namespace detail
0195     {
0196
0197         // Static projection
0198         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_tpeqd, tpeqd_spheroid)
0199
0200         // Factory entry(s)
0201         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(tpeqd_entry, tpeqd_spheroid)
0202
0203         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(tpeqd_init)
0204         {
0205             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(tpeqd, tpeqd_entry)
0206         }
0207
0208     } // namespace detail
0209     #endif // doxygen
0210
0211 } // namespace projections
0212
0213 }} // namespace boost::geometry
0214
0215 #endif // BOOST_GEOMETRY_PROJECTIONS_TPEQD_HPP
0216