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 /***************************************************************************** 41 42 Lambert Conformal Conic Alternative 43 ----------------------------------- 44 45 This is Gerald Evenden's 2003 implementation of an alternative 46 "almost" LCC, which has been in use historically, but which 47 should NOT be used for new projects - i.e: use this implementation 48 if you need interoperability with old data represented in this 49 projection, but not in any other case. 50 51 The code was originally discussed on the PROJ.4 mailing list in 52 a thread archived over at 53 54 http://lists.maptools.org/pipermail/proj/2003-March/000644.html 55 56 It was discussed again in the thread starting at 57 58 http://lists.maptools.org/pipermail/proj/2017-October/007828.html 59 and continuing at 60 http://lists.maptools.org/pipermail/proj/2017-November/007831.html 61 62 which prompted Clifford J. Mugnier to add these clarifying notes: 63 64 The French Army Truncated Cubic Lambert (partially conformal) Conic 65 projection is the Legal system for the projection in France between 66 the late 1800s and 1948 when the French Legislature changed the law 67 to recognize the fully conformal version. 68 69 It was (might still be in one or two North African prior French 70 Colonies) used in North Africa in Algeria, Tunisia, & Morocco, as 71 well as in Syria during the Levant. 72 73 Last time I have seen it used was about 30+ years ago in 74 Algeria when it was used to define Lease Block boundaries for 75 Petroleum Exploration & Production. 76 77 (signed) 78 79 Clifford J. Mugnier, c.p., c.m.s. 80 Chief of Geodesy 81 LSU Center for GeoInformatics 82 Dept. of Civil Engineering 83 LOUISIANA STATE UNIVERSITY 84 85 *****************************************************************************/ 86 87 #ifndef BOOST_GEOMETRY_PROJECTIONS_LCCA_HPP 88 #define BOOST_GEOMETRY_PROJECTIONS_LCCA_HPP 89 90 #include <boost/geometry/srs/projections/impl/base_static.hpp> 91 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> 92 #include <boost/geometry/srs/projections/impl/projects.hpp> 93 #include <boost/geometry/srs/projections/impl/factory_entry.hpp> 94 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp> 95 96 namespace boost { namespace geometry 97 { 98 99 namespace projections 100 { 101 #ifndef DOXYGEN_NO_DETAIL 102 namespace detail { namespace lcca 103 { 104 105 static const int max_iter = 10; 106 static const double del_tol = 1e-12; 107 108 template <typename T> 109 struct par_lcca 110 { 111 detail::en<T> en; 112 T r0, l, M0; 113 T C; 114 }; 115 116 template <typename T> /* func to compute dr */ fS(T const & S,T const & C)117 inline T fS(T const& S, T const& C) 118 { 119 return(S * ( 1. + S * S * C)); 120 } 121 122 template <typename T> /* deriv of fs */ fSp(T const & S,T const & C)123 inline T fSp(T const& S, T const& C) 124 { 125 return(1. + 3.* S * S * C); 126 } 127 128 template <typename T, typename Parameters> 129 struct base_lcca_ellipsoid 130 { 131 par_lcca<T> m_proj_parm; 132 133 // FORWARD(e_forward) ellipsoid 134 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::lcca::base_lcca_ellipsoid135 inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 136 { 137 T S, r, dr; 138 139 S = pj_mlfn(lp_lat, sin(lp_lat), cos(lp_lat), this->m_proj_parm.en) - this->m_proj_parm.M0; 140 dr = fS(S, this->m_proj_parm.C); 141 r = this->m_proj_parm.r0 - dr; 142 xy_x = par.k0 * (r * sin( lp_lon *= this->m_proj_parm.l ) ); 143 xy_y = par.k0 * (this->m_proj_parm.r0 - r * cos(lp_lon) ); 144 } 145 146 // INVERSE(e_inverse) ellipsoid & spheroid 147 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::lcca::base_lcca_ellipsoid148 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const 149 { 150 T theta, dr, S, dif; 151 int i; 152 153 xy_x /= par.k0; 154 xy_y /= par.k0; 155 theta = atan2(xy_x , this->m_proj_parm.r0 - xy_y); 156 dr = xy_y - xy_x * tan(0.5 * theta); 157 lp_lon = theta / this->m_proj_parm.l; 158 S = dr; 159 for (i = max_iter; i ; --i) { 160 S -= (dif = (fS(S, this->m_proj_parm.C) - dr) / fSp(S, this->m_proj_parm.C)); 161 if (fabs(dif) < del_tol) break; 162 } 163 if (!i) { 164 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 165 } 166 lp_lat = pj_inv_mlfn(S + this->m_proj_parm.M0, par.es, this->m_proj_parm.en); 167 } 168 get_nameboost::geometry::projections::detail::lcca::base_lcca_ellipsoid169 static inline std::string get_name() 170 { 171 return "lcca_ellipsoid"; 172 } 173 174 }; 175 176 // Lambert Conformal Conic Alternative 177 template <typename Parameters, typename T> setup_lcca(Parameters const & par,par_lcca<T> & proj_parm)178 inline void setup_lcca(Parameters const& par, par_lcca<T>& proj_parm) 179 { 180 T s2p0, N0, R0, tan0; 181 182 proj_parm.en = pj_enfn<T>(par.es); 183 184 if (par.phi0 == 0.) { 185 BOOST_THROW_EXCEPTION( projection_exception(error_lat_0_is_zero) ); 186 } 187 proj_parm.l = sin(par.phi0); 188 proj_parm.M0 = pj_mlfn(par.phi0, proj_parm.l, cos(par.phi0), proj_parm.en); 189 s2p0 = proj_parm.l * proj_parm.l; 190 R0 = 1. / (1. - par.es * s2p0); 191 N0 = sqrt(R0); 192 R0 *= par.one_es * N0; 193 tan0 = tan(par.phi0); 194 proj_parm.r0 = N0 / tan0; 195 proj_parm.C = 1. / (6. * R0 * N0); 196 } 197 198 }} // namespace detail::lcca 199 #endif // doxygen 200 201 /*! 202 \brief Lambert Conformal Conic Alternative projection 203 \ingroup projections 204 \tparam Geographic latlong point type 205 \tparam Cartesian xy point type 206 \tparam Parameters parameter type 207 \par Projection characteristics 208 - Conic 209 - Spheroid 210 - Ellipsoid 211 \par Projection parameters 212 - lat_0: Latitude of origin 213 \par Example 214 \image html ex_lcca.gif 215 */ 216 template <typename T, typename Parameters> 217 struct lcca_ellipsoid : public detail::lcca::base_lcca_ellipsoid<T, Parameters> 218 { 219 template <typename Params> lcca_ellipsoidboost::geometry::projections::lcca_ellipsoid220 inline lcca_ellipsoid(Params const& , Parameters const& par) 221 { 222 detail::lcca::setup_lcca(par, this->m_proj_parm); 223 } 224 }; 225 226 #ifndef DOXYGEN_NO_DETAIL 227 namespace detail 228 { 229 230 // Static projection BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lcca,lcca_ellipsoid)231 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lcca, lcca_ellipsoid) 232 233 // Factory entry(s) 234 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(lcca_entry, lcca_ellipsoid) 235 236 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(lcca_init) 237 { 238 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(lcca, lcca_entry) 239 } 240 241 } // namespace detail 242 #endif // doxygen 243 244 } // namespace projections 245 246 }} // namespace boost::geometry 247 248 #endif // BOOST_GEOMETRY_PROJECTIONS_LCCA_HPP 249 250