projections/proj/aea.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  // Author: Gerald Evenden (1995)
0019  //         Thomas Knudsen (2016) - revise/add regression tests
0020
0021 // Last updated version of proj: 5.0.0
0022
0023 // Original copyright notice:
0024
0025 // Purpose:  Implementation of the aea (Albers Equal Area) projection.
0026 // Author:   Gerald Evenden
0027 // Copyright (c) 1995, Gerald Evenden
0028
0029 // Permission is hereby granted, free of charge, to any person obtaining a
0030 // copy of this software and associated documentation files (the "Software"),
0031 // to deal in the Software without restriction, including without limitation
0032 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0033 // and/or sell copies of the Software, and to permit persons to whom the
0034 // Software is furnished to do so, subject to the following conditions:
0035
0036 // The above copyright notice and this permission notice shall be included
0037 // in all copies or substantial portions of the Software.
0038
0039 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0040 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0041 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0042 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0043 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0044 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0045 // DEALINGS IN THE SOFTWARE.
0046
0047 #ifndef BOOST_GEOMETRY_PROJECTIONS_AEA_HPP
0048 #define BOOST_GEOMETRY_PROJECTIONS_AEA_HPP
0049
0050 #include <boost/core/ignore_unused.hpp>
0051 #include <boost/geometry/util/math.hpp>
0052 #include <boost/math/special_functions/hypot.hpp>
0053
0054 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0055 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0056 #include <boost/geometry/srs/projections/impl/projects.hpp>
0057 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0058 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
0059 #include <boost/geometry/srs/projections/impl/pj_msfn.hpp>
0060 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0061 #include <boost/geometry/srs/projections/impl/pj_qsfn.hpp>
0062
0063
0064 namespace boost { namespace geometry
0065 {
0066
0067 namespace projections
0068 {
0069     #ifndef DOXYGEN_NO_DETAIL
0070     namespace detail { namespace aea
0071     {
0072
0073             static const double epsilon10 = 1.e-10;
0074             static const double tolerance7 = 1.e-7;
0075             static const double epsilon = 1.0e-7;
0076             static const double tolerance = 1.0e-10;
0077             static const int n_iter = 15;
0078
0079             template <typename T>
0080             struct par_aea
0081             {
0082                 T    ec;
0083                 T    n;
0084                 T    c;
0085                 T    dd;
0086                 T    n2;
0087                 T    rho0;
0088                 T    phi1;
0089                 T    phi2;
0090                 detail::en<T> en;
0091                 bool ellips;
0092             };
0093
0094             /* determine latitude angle phi-1 */
0095             template <typename T>
0096             inline T phi1_(T const& qs, T const& Te, T const& Tone_es)
0097             {
0098                 int i;
0099                 T Phi, sinpi, cospi, con, com, dphi;
0100
0101                 Phi = asin (.5 * qs);
0102                 if (Te < epsilon)
0103                     return( Phi );
0104                 i = n_iter;
0105                 do {
0106                     sinpi = sin (Phi);
0107                     cospi = cos (Phi);
0108                     con = Te * sinpi;
0109                     com = 1. - con * con;
0110                     dphi = .5 * com * com / cospi * (qs / Tone_es -
0111                        sinpi / com + .5 / Te * log ((1. - con) /
0112                        (1. + con)));
0113                     Phi += dphi;
0114                 } while (fabs(dphi) > tolerance && --i);
0115                 return( i ? Phi : HUGE_VAL );
0116             }
0117
0118             template <typename T, typename Parameters>
0119             struct base_aea_ellipsoid
0120             {
0121                 par_aea<T> m_proj_parm;
0122
0123                 // FORWARD(e_forward)  ellipsoid & spheroid
0124                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0125                 inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0126                 {
0127                     T rho = this->m_proj_parm.c - (this->m_proj_parm.ellips
0128                                                                     ? this->m_proj_parm.n * pj_qsfn(sin(lp_lat), par.e, par.one_es)
0129                                                                     : this->m_proj_parm.n2 * sin(lp_lat));
0130                     if (rho < 0.)
0131                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0132                     rho = this->m_proj_parm.dd * sqrt(rho);
0133                     xy_x = rho * sin( lp_lon *= this->m_proj_parm.n );
0134                     xy_y = this->m_proj_parm.rho0 - rho * cos(lp_lon);
0135                 }
0136
0137                 // INVERSE(e_inverse)  ellipsoid & spheroid
0138                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0139                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0140                 {
0141                     static const T half_pi = detail::half_pi<T>();
0142
0143                     T rho = 0.0;
0144                     if( (rho = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.rho0 - xy_y)) != 0.0 ) {
0145                         if (this->m_proj_parm.n < 0.) {
0146                             rho = -rho;
0147                             xy_x = -xy_x;
0148                             xy_y = -xy_y;
0149                         }
0150                         lp_lat =  rho / this->m_proj_parm.dd;
0151                         if (this->m_proj_parm.ellips) {
0152                             lp_lat = (this->m_proj_parm.c - lp_lat * lp_lat) / this->m_proj_parm.n;
0153                             if (fabs(this->m_proj_parm.ec - fabs(lp_lat)) > tolerance7) {
0154                                 if ((lp_lat = phi1_(lp_lat, par.e, par.one_es)) == HUGE_VAL)
0155                                     BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0156                             } else
0157                                 lp_lat = lp_lat < 0. ? -half_pi : half_pi;
0158                         } else if (fabs(lp_lat = (this->m_proj_parm.c - lp_lat * lp_lat) / this->m_proj_parm.n2) <= 1.)
0159                             lp_lat = asin(lp_lat);
0160                         else
0161                             lp_lat = lp_lat < 0. ? -half_pi : half_pi;
0162                         lp_lon = atan2(xy_x, xy_y) / this->m_proj_parm.n;
0163                     } else {
0164                         lp_lon = 0.;
0165                         lp_lat = this->m_proj_parm.n > 0. ? half_pi : - half_pi;
0166                     }
0167                 }
0168
0169                 static inline std::string get_name()
0170                 {
0171                     return "aea_ellipsoid";
0172                 }
0173
0174             };
0175
0176             template <typename Parameters, typename T>
0177             inline void setup(Parameters const& par, par_aea<T>& proj_parm)
0178             {
0179                 T cosphi, sinphi;
0180                 int secant;
0181
0182                 if (fabs(proj_parm.phi1 + proj_parm.phi2) < epsilon10)
0183                     BOOST_THROW_EXCEPTION( projection_exception(error_conic_lat_equal) );
0184                 proj_parm.n = sinphi = sin(proj_parm.phi1);
0185                 cosphi = cos(proj_parm.phi1);
0186                 secant = fabs(proj_parm.phi1 - proj_parm.phi2) >= epsilon10;
0187                 if( (proj_parm.ellips = (par.es > 0.))) {
0188                     T ml1, m1;
0189
0190                     proj_parm.en = pj_enfn<T>(par.es);
0191                     m1 = pj_msfn(sinphi, cosphi, par.es);
0192                     ml1 = pj_qsfn(sinphi, par.e, par.one_es);
0193                     if (secant) { /* secant cone */
0194                         T ml2, m2;
0195
0196                         sinphi = sin(proj_parm.phi2);
0197                         cosphi = cos(proj_parm.phi2);
0198                         m2 = pj_msfn(sinphi, cosphi, par.es);
0199                         ml2 = pj_qsfn(sinphi, par.e, par.one_es);
0200                         if (ml2 == ml1)
0201                             BOOST_THROW_EXCEPTION( projection_exception(0) );
0202
0203                         proj_parm.n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
0204                     }
0205                     proj_parm.ec = 1. - .5 * par.one_es * log((1. - par.e) /
0206                         (1. + par.e)) / par.e;
0207                     proj_parm.c = m1 * m1 + proj_parm.n * ml1;
0208                     proj_parm.dd = 1. / proj_parm.n;
0209                     proj_parm.rho0 = proj_parm.dd * sqrt(proj_parm.c - proj_parm.n * pj_qsfn(sin(par.phi0),
0210                         par.e, par.one_es));
0211                 } else {
0212                     if (secant) proj_parm.n = .5 * (proj_parm.n + sin(proj_parm.phi2));
0213                     proj_parm.n2 = proj_parm.n + proj_parm.n;
0214                     proj_parm.c = cosphi * cosphi + proj_parm.n2 * sinphi;
0215                     proj_parm.dd = 1. / proj_parm.n;
0216                     proj_parm.rho0 = proj_parm.dd * sqrt(proj_parm.c - proj_parm.n2 * sin(par.phi0));
0217                 }
0218             }
0219
0220
0221             // Albers Equal Area
0222             template <typename Params, typename Parameters, typename T>
0223             inline void setup_aea(Params const& params, Parameters const& par, par_aea<T>& proj_parm)
0224             {
0225                 proj_parm.phi1 = 0.0;
0226                 proj_parm.phi2 = 0.0;
0227                 bool is_phi1_set = pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, proj_parm.phi1);
0228                 bool is_phi2_set = pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, proj_parm.phi2);
0229
0230                 // Boost.Geometry specific, set default parameters manually
0231                 if (! is_phi1_set || ! is_phi2_set) {
0232                     bool const use_defaults = ! pj_get_param_b<srs::spar::no_defs>(params, "no_defs", srs::dpar::no_defs);
0233                     if (use_defaults) {
0234                         if (!is_phi1_set)
0235                             proj_parm.phi1 = 29.5;
0236                         if (!is_phi2_set)
0237                             proj_parm.phi2 = 45.5;
0238                     }
0239                 }
0240
0241                 setup(par, proj_parm);
0242             }
0243
0244             // Lambert Equal Area Conic
0245             template <typename Params, typename Parameters, typename T>
0246             inline void setup_leac(Params const& params, Parameters const& par, par_aea<T>& proj_parm)
0247             {
0248                 static const T half_pi = detail::half_pi<T>();
0249
0250                 proj_parm.phi2 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
0251                 proj_parm.phi1 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? -half_pi : half_pi;
0252                 setup(par, proj_parm);
0253             }
0254
0255     }} // namespace detail::aea
0256     #endif // doxygen
0257
0258     /*!
0259         \brief Albers Equal Area projection
0260         \ingroup projections
0261         \tparam Geographic latlong point type
0262         \tparam Cartesian xy point type
0263         \tparam Parameters parameter type
0264         \par Projection characteristics
0265          - Conic
0266          - Spheroid
0267          - Ellipsoid
0268         \par Projection parameters
0269          - lat_1: Latitude of first standard parallel (degrees)
0270          - lat_2: Latitude of second standard parallel (degrees)
0271         \par Example
0272         \image html ex_aea.gif
0273     */
0274     template <typename T, typename Parameters>
0275     struct aea_ellipsoid : public detail::aea::base_aea_ellipsoid<T, Parameters>
0276     {
0277         template <typename Params>
0278         inline aea_ellipsoid(Params const& params, Parameters const& par)
0279         {
0280             detail::aea::setup_aea(params, par, this->m_proj_parm);
0281         }
0282     };
0283
0284     /*!
0285         \brief Lambert Equal Area Conic projection
0286         \ingroup projections
0287         \tparam Geographic latlong point type
0288         \tparam Cartesian xy point type
0289         \tparam Parameters parameter type
0290         \par Projection characteristics
0291          - Conic
0292          - Spheroid
0293          - Ellipsoid
0294         \par Projection parameters
0295          - lat_1: Latitude of first standard parallel (degrees)
0296          - south: Denotes southern hemisphere UTM zone (boolean)
0297         \par Example
0298         \image html ex_leac.gif
0299     */
0300     template <typename T, typename Parameters>
0301     struct leac_ellipsoid : public detail::aea::base_aea_ellipsoid<T, Parameters>
0302     {
0303         template <typename Params>
0304         inline leac_ellipsoid(Params const& params, Parameters const& par)
0305         {
0306             detail::aea::setup_leac(params, par, this->m_proj_parm);
0307         }
0308     };
0309
0310     #ifndef DOXYGEN_NO_DETAIL
0311     namespace detail
0312     {
0313
0314         // Static projection
0315         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_aea, aea_ellipsoid)
0316         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_leac, leac_ellipsoid)
0317
0318         // Factory entry(s)
0319         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(aea_entry, aea_ellipsoid)
0320         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(leac_entry, leac_ellipsoid)
0321
0322         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(aea_init)
0323         {
0324             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(aea, aea_entry)
0325             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(leac, leac_entry)
0326         }
0327
0328     } // namespace detail
0329     #endif // doxygen
0330
0331 } // namespace projections
0332
0333 }} // namespace boost::geometry
0334
0335 #endif // BOOST_GEOMETRY_PROJECTIONS_AEA_HPP
0336