boost/geometry/srs/projections/proj/aea.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  // Author: Gerald Evenden (1995)
  19  //         Thomas Knudsen (2016) - revise/add regression tests
  20
  21 // Last updated version of proj: 5.0.0
  22
  23 // Original copyright notice:
  24
  25 // Purpose:  Implementation of the aea (Albers Equal Area) projection.
  26 // Author:   Gerald Evenden
  27 // Copyright (c) 1995, Gerald Evenden
  28
  29 // Permission is hereby granted, free of charge, to any person obtaining a
  30 // copy of this software and associated documentation files (the "Software"),
  31 // to deal in the Software without restriction, including without limitation
  32 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  33 // and/or sell copies of the Software, and to permit persons to whom the
  34 // Software is furnished to do so, subject to the following conditions:
  35
  36 // The above copyright notice and this permission notice shall be included
  37 // in all copies or substantial portions of the Software.
  38
  39 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  40 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  41 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  42 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  43 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  44 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  45 // DEALINGS IN THE SOFTWARE.
  46
  47 #ifndef BOOST_GEOMETRY_PROJECTIONS_AEA_HPP
  48 #define BOOST_GEOMETRY_PROJECTIONS_AEA_HPP
  49
  50 #include <boost/core/ignore_unused.hpp>
  51 #include <boost/geometry/util/math.hpp>
  52 #include <boost/math/special_functions/hypot.hpp>
  53
  54 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  55 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  56 #include <boost/geometry/srs/projections/impl/projects.hpp>
  57 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  58 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
  59 #include <boost/geometry/srs/projections/impl/pj_msfn.hpp>
  60 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  61 #include <boost/geometry/srs/projections/impl/pj_qsfn.hpp>
  62
  63
  64 namespace boost { namespace geometry
  65 {
  66
  67 namespace projections
  68 {
  69     #ifndef DOXYGEN_NO_DETAIL
  70     namespace detail { namespace aea
  71     {
  72
  73             static const double epsilon10 = 1.e-10;
  74             static const double tolerance7 = 1.e-7;
  75             static const double epsilon = 1.0e-7;
  76             static const double tolerance = 1.0e-10;
  77             static const int n_iter = 15;
  78
  79             template <typename T>
  80             struct par_aea
  81             {
  82                 T    ec;
  83                 T    n;
  84                 T    c;
  85                 T    dd;
  86                 T    n2;
  87                 T    rho0;
  88                 T    phi1;
  89                 T    phi2;
  90                 detail::en<T> en;
  91                 bool ellips;
  92             };
  93
  94             /* determine latitude angle phi-1 */
  95             template <typename T>
  96             inline T phi1_(T const& qs, T const& Te, T const& Tone_es)
  97             {
  98                 int i;
  99                 T Phi, sinpi, cospi, con, com, dphi;
 100
 101                 Phi = asin (.5 * qs);
 102                 if (Te < epsilon)
 103                     return( Phi );
 104                 i = n_iter;
 105                 do {
 106                     sinpi = sin (Phi);
 107                     cospi = cos (Phi);
 108                     con = Te * sinpi;
 109                     com = 1. - con * con;
 110                     dphi = .5 * com * com / cospi * (qs / Tone_es -
 111                        sinpi / com + .5 / Te * log ((1. - con) /
 112                        (1. + con)));
 113                     Phi += dphi;
 114                 } while (fabs(dphi) > tolerance && --i);
 115                 return( i ? Phi : HUGE_VAL );
 116             }
 117
 118             template <typename T, typename Parameters>
 119             struct base_aea_ellipsoid
 120             {
 121                 par_aea<T> m_proj_parm;
 122
 123                 // FORWARD(e_forward)  ellipsoid & spheroid
 124                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 125                 inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
 126                 {
 127                     T rho = this->m_proj_parm.c - (this->m_proj_parm.ellips
 128                                                                     ? this->m_proj_parm.n * pj_qsfn(sin(lp_lat), par.e, par.one_es)
 129                                                                     : this->m_proj_parm.n2 * sin(lp_lat));
 130                     if (rho < 0.)
 131                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 132                     rho = this->m_proj_parm.dd * sqrt(rho);
 133                     xy_x = rho * sin( lp_lon *= this->m_proj_parm.n );
 134                     xy_y = this->m_proj_parm.rho0 - rho * cos(lp_lon);
 135                 }
 136
 137                 // INVERSE(e_inverse)  ellipsoid & spheroid
 138                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 139                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 140                 {
 141                     static const T half_pi = detail::half_pi<T>();
 142
 143                     T rho = 0.0;
 144                     if( (rho = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.rho0 - xy_y)) != 0.0 ) {
 145                         if (this->m_proj_parm.n < 0.) {
 146                             rho = -rho;
 147                             xy_x = -xy_x;
 148                             xy_y = -xy_y;
 149                         }
 150                         lp_lat =  rho / this->m_proj_parm.dd;
 151                         if (this->m_proj_parm.ellips) {
 152                             lp_lat = (this->m_proj_parm.c - lp_lat * lp_lat) / this->m_proj_parm.n;
 153                             if (fabs(this->m_proj_parm.ec - fabs(lp_lat)) > tolerance7) {
 154                                 if ((lp_lat = phi1_(lp_lat, par.e, par.one_es)) == HUGE_VAL)
 155                                     BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 156                             } else
 157                                 lp_lat = lp_lat < 0. ? -half_pi : half_pi;
 158                         } else if (fabs(lp_lat = (this->m_proj_parm.c - lp_lat * lp_lat) / this->m_proj_parm.n2) <= 1.)
 159                             lp_lat = asin(lp_lat);
 160                         else
 161                             lp_lat = lp_lat < 0. ? -half_pi : half_pi;
 162                         lp_lon = atan2(xy_x, xy_y) / this->m_proj_parm.n;
 163                     } else {
 164                         lp_lon = 0.;
 165                         lp_lat = this->m_proj_parm.n > 0. ? half_pi : - half_pi;
 166                     }
 167                 }
 168
 169                 static inline std::string get_name()
 170                 {
 171                     return "aea_ellipsoid";
 172                 }
 173
 174             };
 175
 176             template <typename Parameters, typename T>
 177             inline void setup(Parameters const& par, par_aea<T>& proj_parm)
 178             {
 179                 T cosphi, sinphi;
 180                 int secant;
 181
 182                 if (fabs(proj_parm.phi1 + proj_parm.phi2) < epsilon10)
 183                     BOOST_THROW_EXCEPTION( projection_exception(error_conic_lat_equal) );
 184                 proj_parm.n = sinphi = sin(proj_parm.phi1);
 185                 cosphi = cos(proj_parm.phi1);
 186                 secant = fabs(proj_parm.phi1 - proj_parm.phi2) >= epsilon10;
 187                 if( (proj_parm.ellips = (par.es > 0.))) {
 188                     T ml1, m1;
 189
 190                     proj_parm.en = pj_enfn<T>(par.es);
 191                     m1 = pj_msfn(sinphi, cosphi, par.es);
 192                     ml1 = pj_qsfn(sinphi, par.e, par.one_es);
 193                     if (secant) { /* secant cone */
 194                         T ml2, m2;
 195
 196                         sinphi = sin(proj_parm.phi2);
 197                         cosphi = cos(proj_parm.phi2);
 198                         m2 = pj_msfn(sinphi, cosphi, par.es);
 199                         ml2 = pj_qsfn(sinphi, par.e, par.one_es);
 200                         if (ml2 == ml1)
 201                             BOOST_THROW_EXCEPTION( projection_exception(0) );
 202
 203                         proj_parm.n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
 204                     }
 205                     proj_parm.ec = 1. - .5 * par.one_es * log((1. - par.e) /
 206                         (1. + par.e)) / par.e;
 207                     proj_parm.c = m1 * m1 + proj_parm.n * ml1;
 208                     proj_parm.dd = 1. / proj_parm.n;
 209                     proj_parm.rho0 = proj_parm.dd * sqrt(proj_parm.c - proj_parm.n * pj_qsfn(sin(par.phi0),
 210                         par.e, par.one_es));
 211                 } else {
 212                     if (secant) proj_parm.n = .5 * (proj_parm.n + sin(proj_parm.phi2));
 213                     proj_parm.n2 = proj_parm.n + proj_parm.n;
 214                     proj_parm.c = cosphi * cosphi + proj_parm.n2 * sinphi;
 215                     proj_parm.dd = 1. / proj_parm.n;
 216                     proj_parm.rho0 = proj_parm.dd * sqrt(proj_parm.c - proj_parm.n2 * sin(par.phi0));
 217                 }
 218             }
 219
 220
 221             // Albers Equal Area
 222             template <typename Params, typename Parameters, typename T>
 223             inline void setup_aea(Params const& params, Parameters const& par, par_aea<T>& proj_parm)
 224             {
 225                 proj_parm.phi1 = 0.0;
 226                 proj_parm.phi2 = 0.0;
 227                 bool is_phi1_set = pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, proj_parm.phi1);
 228                 bool is_phi2_set = pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, proj_parm.phi2);
 229
 230                 // Boost.Geometry specific, set default parameters manually
 231                 if (! is_phi1_set || ! is_phi2_set) {
 232                     bool const use_defaults = ! pj_get_param_b<srs::spar::no_defs>(params, "no_defs", srs::dpar::no_defs);
 233                     if (use_defaults) {
 234                         if (!is_phi1_set)
 235                             proj_parm.phi1 = 29.5;
 236                         if (!is_phi2_set)
 237                             proj_parm.phi2 = 45.5;
 238                     }
 239                 }
 240
 241                 setup(par, proj_parm);
 242             }
 243
 244             // Lambert Equal Area Conic
 245             template <typename Params, typename Parameters, typename T>
 246             inline void setup_leac(Params const& params, Parameters const& par, par_aea<T>& proj_parm)
 247             {
 248                 static const T half_pi = detail::half_pi<T>();
 249
 250                 proj_parm.phi2 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
 251                 proj_parm.phi1 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? -half_pi : half_pi;
 252                 setup(par, proj_parm);
 253             }
 254
 255     }} // namespace detail::aea
 256     #endif // doxygen
 257
 258     /*!
 259         \brief Albers Equal Area projection
 260         \ingroup projections
 261         \tparam Geographic latlong point type
 262         \tparam Cartesian xy point type
 263         \tparam Parameters parameter type
 264         \par Projection characteristics
 265          - Conic
 266          - Spheroid
 267          - Ellipsoid
 268         \par Projection parameters
 269          - lat_1: Latitude of first standard parallel (degrees)
 270          - lat_2: Latitude of second standard parallel (degrees)
 271         \par Example
 272         \image html ex_aea.gif
 273     */
 274     template <typename T, typename Parameters>
 275     struct aea_ellipsoid : public detail::aea::base_aea_ellipsoid<T, Parameters>
 276     {
 277         template <typename Params>
 278         inline aea_ellipsoid(Params const& params, Parameters const& par)
 279         {
 280             detail::aea::setup_aea(params, par, this->m_proj_parm);
 281         }
 282     };
 283
 284     /*!
 285         \brief Lambert Equal Area Conic projection
 286         \ingroup projections
 287         \tparam Geographic latlong point type
 288         \tparam Cartesian xy point type
 289         \tparam Parameters parameter type
 290         \par Projection characteristics
 291          - Conic
 292          - Spheroid
 293          - Ellipsoid
 294         \par Projection parameters
 295          - lat_1: Latitude of first standard parallel (degrees)
 296          - south: Denotes southern hemisphere UTM zone (boolean)
 297         \par Example
 298         \image html ex_leac.gif
 299     */
 300     template <typename T, typename Parameters>
 301     struct leac_ellipsoid : public detail::aea::base_aea_ellipsoid<T, Parameters>
 302     {
 303         template <typename Params>
 304         inline leac_ellipsoid(Params const& params, Parameters const& par)
 305         {
 306             detail::aea::setup_leac(params, par, this->m_proj_parm);
 307         }
 308     };
 309
 310     #ifndef DOXYGEN_NO_DETAIL
 311     namespace detail
 312     {
 313
 314         // Static projection
 315         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_aea, aea_ellipsoid)
 316         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_leac, leac_ellipsoid)
 317
 318         // Factory entry(s)
 319         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(aea_entry, aea_ellipsoid)
 320         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(leac_entry, leac_ellipsoid)
 321
 322         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(aea_init)
 323         {
 324             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(aea, aea_entry)
 325             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(leac, leac_entry)
 326         }
 327
 328     } // namespace detail
 329     #endif // doxygen
 330
 331 } // namespace projections
 332
 333 }} // namespace boost::geometry
 334
 335 #endif // BOOST_GEOMETRY_PROJECTIONS_AEA_HPP
 336