]>
Commit | Line | Data |
---|---|---|
92f5a8d4 | 1 | // Boost.Geometry - gis-projections (based on PROJ4) |
11fdf7f2 TL |
2 | |
3 | // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. | |
4 | ||
92f5a8d4 TL |
5 | // This file was modified by Oracle on 2017, 2018, 2019. |
6 | // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates. | |
11fdf7f2 TL |
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 | ||
92f5a8d4 | 18 | // Last updated version of proj: 5.0.0 |
11fdf7f2 TL |
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 | ||
92f5a8d4 TL |
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 | |
11fdf7f2 TL |
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 | ||
11fdf7f2 TL |
99 | namespace projections |
100 | { | |
101 | #ifndef DOXYGEN_NO_DETAIL | |
102 | namespace detail { namespace lcca | |
103 | { | |
104 | ||
92f5a8d4 TL |
105 | static const int max_iter = 10; |
106 | static const double del_tol = 1e-12; | |
11fdf7f2 TL |
107 | |
108 | template <typename T> | |
109 | struct par_lcca | |
110 | { | |
92f5a8d4 | 111 | detail::en<T> en; |
11fdf7f2 TL |
112 | T r0, l, M0; |
113 | T C; | |
114 | }; | |
115 | ||
116 | template <typename T> /* func to compute dr */ | |
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 */ | |
123 | inline T fSp(T const& S, T const& C) | |
124 | { | |
125 | return(1. + 3.* S * S * C); | |
126 | } | |
127 | ||
92f5a8d4 TL |
128 | template <typename T, typename Parameters> |
129 | struct base_lcca_ellipsoid | |
11fdf7f2 | 130 | { |
92f5a8d4 | 131 | par_lcca<T> m_proj_parm; |
11fdf7f2 TL |
132 | |
133 | // FORWARD(e_forward) ellipsoid | |
134 | // Project coordinates from geographic (lon, lat) to cartesian (x, y) | |
92f5a8d4 | 135 | inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const |
11fdf7f2 | 136 | { |
92f5a8d4 | 137 | T S, r, dr; |
11fdf7f2 TL |
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; | |
92f5a8d4 TL |
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) ); | |
11fdf7f2 TL |
144 | } |
145 | ||
146 | // INVERSE(e_inverse) ellipsoid & spheroid | |
147 | // Project coordinates from cartesian (x, y) to geographic (lon, lat) | |
92f5a8d4 | 148 | inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const |
11fdf7f2 | 149 | { |
92f5a8d4 | 150 | T theta, dr, S, dif; |
11fdf7f2 TL |
151 | int i; |
152 | ||
92f5a8d4 TL |
153 | xy_x /= par.k0; |
154 | xy_y /= par.k0; | |
11fdf7f2 TL |
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; | |
92f5a8d4 | 159 | for (i = max_iter; i ; --i) { |
11fdf7f2 | 160 | S -= (dif = (fS(S, this->m_proj_parm.C) - dr) / fSp(S, this->m_proj_parm.C)); |
92f5a8d4 TL |
161 | if (fabs(dif) < del_tol) break; |
162 | } | |
163 | if (!i) { | |
164 | BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); | |
11fdf7f2 | 165 | } |
92f5a8d4 | 166 | lp_lat = pj_inv_mlfn(S + this->m_proj_parm.M0, par.es, this->m_proj_parm.en); |
11fdf7f2 TL |
167 | } |
168 | ||
169 | 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> | |
92f5a8d4 | 178 | inline void setup_lcca(Parameters const& par, par_lcca<T>& proj_parm) |
11fdf7f2 | 179 | { |
92f5a8d4 TL |
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 | } | |
11fdf7f2 TL |
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); | |
11fdf7f2 TL |
194 | proj_parm.r0 = N0 / tan0; |
195 | proj_parm.C = 1. / (6. * R0 * N0); | |
11fdf7f2 TL |
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 | */ | |
92f5a8d4 TL |
216 | template <typename T, typename Parameters> |
217 | struct lcca_ellipsoid : public detail::lcca::base_lcca_ellipsoid<T, Parameters> | |
11fdf7f2 | 218 | { |
92f5a8d4 TL |
219 | template <typename Params> |
220 | inline lcca_ellipsoid(Params const& , Parameters const& par) | |
11fdf7f2 | 221 | { |
92f5a8d4 | 222 | detail::lcca::setup_lcca(par, this->m_proj_parm); |
11fdf7f2 TL |
223 | } |
224 | }; | |
225 | ||
226 | #ifndef DOXYGEN_NO_DETAIL | |
227 | namespace detail | |
228 | { | |
229 | ||
230 | // Static projection | |
92f5a8d4 | 231 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lcca, lcca_ellipsoid) |
11fdf7f2 TL |
232 | |
233 | // Factory entry(s) | |
92f5a8d4 TL |
234 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(lcca_entry, lcca_ellipsoid) |
235 | ||
236 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(lcca_init) | |
11fdf7f2 | 237 | { |
92f5a8d4 | 238 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(lcca, lcca_entry) |
11fdf7f2 TL |
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 |