1 // Boost.Geometry - gis-projections (based on PROJ4)
3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
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.
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)
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
18 // Last updated version of proj: 5.0.0
20 // Original copyright notice:
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:
29 // The above copyright notice and this permission notice shall be included
30 // in all copies or substantial portions of the Software.
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.
40 #ifndef BOOST_GEOMETRY_PROJECTIONS_LCC_HPP
41 #define BOOST_GEOMETRY_PROJECTIONS_LCC_HPP
43 #include <boost/geometry/srs/projections/impl/base_static.hpp>
44 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
45 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
46 #include <boost/geometry/srs/projections/impl/pj_msfn.hpp>
47 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
48 #include <boost/geometry/srs/projections/impl/pj_phi2.hpp>
49 #include <boost/geometry/srs/projections/impl/pj_tsfn.hpp>
50 #include <boost/geometry/srs/projections/impl/projects.hpp>
52 #include <boost/geometry/util/math.hpp>
54 #include <boost/math/special_functions/hypot.hpp>
56 namespace boost { namespace geometry
61 #ifndef DOXYGEN_NO_DETAIL
62 namespace detail { namespace lcc
64 static const double epsilon10 = 1.e-10;
77 template <typename T, typename Parameters>
78 struct base_lcc_ellipsoid
80 par_lcc<T> m_proj_parm;
82 // FORWARD(e_forward) ellipsoid & spheroid
83 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
84 inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
86 static const T fourth_pi = detail::fourth_pi<T>();
87 static const T half_pi = detail::half_pi<T>();
91 if (fabs(fabs(lp_lat) - half_pi) < epsilon10) {
92 if ((lp_lat * this->m_proj_parm.n) <= 0.) {
93 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
97 rho = this->m_proj_parm.c * (this->m_proj_parm.ellips
98 ? math::pow(pj_tsfn(lp_lat, sin(lp_lat), par.e), this->m_proj_parm.n)
99 : math::pow(tan(fourth_pi + T(0.5) * lp_lat), -this->m_proj_parm.n));
101 lp_lon *= this->m_proj_parm.n;
102 xy_x = par.k0 * (rho * sin( lp_lon) );
103 xy_y = par.k0 * (this->m_proj_parm.rho0 - rho * cos(lp_lon) );
106 // INVERSE(e_inverse) ellipsoid & spheroid
107 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
108 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
110 static const T half_pi = detail::half_pi<T>();
117 xy_y = this->m_proj_parm.rho0 - xy_y;
118 rho = boost::math::hypot(xy_x, xy_y);
120 if (this->m_proj_parm.n < 0.) {
125 if (this->m_proj_parm.ellips) {
126 lp_lat = pj_phi2(math::pow(rho / this->m_proj_parm.c, T(1)/this->m_proj_parm.n), par.e);
127 if (lp_lat == HUGE_VAL) {
128 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
131 lp_lat = 2. * atan(math::pow(this->m_proj_parm.c / rho, T(1)/this->m_proj_parm.n)) - half_pi;
132 lp_lon = atan2(xy_x, xy_y) / this->m_proj_parm.n;
135 lp_lat = this->m_proj_parm.n > 0. ? half_pi : -half_pi;
139 static inline std::string get_name()
141 return "lcc_ellipsoid";
146 // Lambert Conformal Conic
147 template <typename Params, typename Parameters, typename T>
148 inline void setup_lcc(Params const& params, Parameters& par, par_lcc<T>& proj_parm)
150 static const T fourth_pi = detail::fourth_pi<T>();
151 static const T half_pi = detail::half_pi<T>();
156 proj_parm.phi1 = 0.0;
157 proj_parm.phi2 = 0.0;
158 bool is_phi1_set = pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, proj_parm.phi1);
159 bool is_phi2_set = pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, proj_parm.phi2);
161 // Boost.Geometry specific, set default parameters manually
162 if (! is_phi1_set || ! is_phi2_set) {
163 bool const use_defaults = ! pj_get_param_b<srs::spar::no_defs>(params, "no_defs", srs::dpar::no_defs);
177 proj_parm.phi2 = proj_parm.phi1;
178 if (! pj_param_exists<srs::spar::lat_0>(params, "lat_0", srs::dpar::lat_0))
179 par.phi0 = proj_parm.phi1;
181 if (fabs(proj_parm.phi1 + proj_parm.phi2) < epsilon10)
182 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.)) ) {
190 par.e = sqrt(par.es); // TODO: Isn't it already set?
191 m1 = pj_msfn(sinphi, cosphi, par.es);
192 ml1 = pj_tsfn(proj_parm.phi1, sinphi, par.e);
193 if (secant) { /* secant cone */
194 sinphi = sin(proj_parm.phi2);
195 proj_parm.n = log(m1 / pj_msfn(sinphi, cos(proj_parm.phi2), par.es));
196 proj_parm.n /= log(ml1 / pj_tsfn(proj_parm.phi2, sinphi, par.e));
198 proj_parm.c = (proj_parm.rho0 = m1 * math::pow(ml1, -proj_parm.n) / proj_parm.n);
199 proj_parm.rho0 *= (fabs(fabs(par.phi0) - half_pi) < epsilon10) ? T(0) :
200 math::pow(pj_tsfn(par.phi0, sin(par.phi0), par.e), proj_parm.n);
203 proj_parm.n = log(cosphi / cos(proj_parm.phi2)) /
204 log(tan(fourth_pi + .5 * proj_parm.phi2) /
205 tan(fourth_pi + .5 * proj_parm.phi1));
206 proj_parm.c = cosphi * math::pow(tan(fourth_pi + T(0.5) * proj_parm.phi1), proj_parm.n) / proj_parm.n;
207 proj_parm.rho0 = (fabs(fabs(par.phi0) - half_pi) < epsilon10) ? 0. :
208 proj_parm.c * math::pow(tan(fourth_pi + T(0.5) * par.phi0), -proj_parm.n);
212 }} // namespace detail::lcc
216 \brief Lambert Conformal Conic projection
218 \tparam Geographic latlong point type
219 \tparam Cartesian xy point type
220 \tparam Parameters parameter type
221 \par Projection characteristics
225 \par Projection parameters
226 - lat_1: Latitude of first standard parallel (degrees)
227 - lat_2: Latitude of second standard parallel (degrees)
228 - lat_0: Latitude of origin
230 \image html ex_lcc.gif
232 template <typename T, typename Parameters>
233 struct lcc_ellipsoid : public detail::lcc::base_lcc_ellipsoid<T, Parameters>
235 template <typename Params>
236 inline lcc_ellipsoid(Params const& params, Parameters & par)
238 detail::lcc::setup_lcc(params, par, this->m_proj_parm);
242 #ifndef DOXYGEN_NO_DETAIL
247 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lcc, lcc_ellipsoid)
250 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(lcc_entry, lcc_ellipsoid)
252 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(lcc_init)
254 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(lcc, lcc_entry);
257 } // namespace detail
260 } // namespace projections
262 }} // namespace boost::geometry
264 #endif // BOOST_GEOMETRY_PROJECTIONS_LCC_HPP