00001 /** 00002 * \file LambertConformalConic.hpp 00003 * \brief Header for GeographicLib::LambertConformalConic class 00004 * 00005 * Copyright (c) Charles Karney (2010-2014) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP) 00011 #define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP 1 00012 00013 #include <GeographicLib/Constants.hpp> 00014 00015 namespace GeographicLib { 00016 00017 /** 00018 * \brief Lambert conformal conic projection 00019 * 00020 * Implementation taken from the report, 00021 * - J. P. Snyder, 00022 * <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projections: A 00023 * Working Manual</a>, USGS Professional Paper 1395 (1987), 00024 * pp. 107--109. 00025 * 00026 * This is a implementation of the equations in Snyder except that divided 00027 * differences have been used to transform the expressions into ones which 00028 * may be evaluated accurately and that Newton's method is used to invert the 00029 * projection. In this implementation, the projection correctly becomes the 00030 * Mercator projection or the polar stereographic projection when the 00031 * standard latitude is the equator or a pole. The accuracy of the 00032 * projections is about 10 nm (10 nanometers). 00033 * 00034 * The ellipsoid parameters, the standard parallels, and the scale on the 00035 * standard parallels are set in the constructor. Internally, the case with 00036 * two standard parallels is converted into a single standard parallel, the 00037 * latitude of tangency (also the latitude of minimum scale), with a scale 00038 * specified on this parallel. This latitude is also used as the latitude of 00039 * origin which is returned by LambertConformalConic::OriginLatitude. The 00040 * scale on the latitude of origin is given by 00041 * LambertConformalConic::CentralScale. The case with two distinct standard 00042 * parallels where one is a pole is singular and is disallowed. The central 00043 * meridian (which is a trivial shift of the longitude) is specified as the 00044 * \e lon0 argument of the LambertConformalConic::Forward and 00045 * LambertConformalConic::Reverse functions. There is no provision in this 00046 * class for specifying a false easting or false northing or a different 00047 * latitude of origin. However these are can be simply included by the 00048 * calling function. For example the Pennsylvania South state coordinate 00049 * system (<a href="http://www.spatialreference.org/ref/epsg/3364/"> 00050 * EPSG:3364</a>) is obtained by: 00051 * \include example-LambertConformalConic.cpp 00052 * 00053 * <a href="ConicProj.1.html">ConicProj</a> is a command-line utility 00054 * providing access to the functionality of LambertConformalConic and 00055 * AlbersEqualArea. 00056 **********************************************************************/ 00057 class GEOGRAPHICLIB_EXPORT LambertConformalConic { 00058 private: 00059 typedef Math::real real; 00060 real eps_, epsx_, tol_, ahypover_; 00061 real _a, _f, _fm, _e2, _e, _e2m; 00062 real _sign, _n, _nc, _t0nm1, _scale, _lat0, _k0; 00063 real _scbet0, _tchi0, _scchi0, _psi0, _nrho0, _drhomax; 00064 static const int numit_ = 5; 00065 static inline real hyp(real x) { return Math::hypot(real(1), x); } 00066 // e * atanh(e * x) = log( ((1 + e*x)/(1 - e*x))^(e/2) ) if f >= 0 00067 // - sqrt(-e2) * atan( sqrt(-e2) * x) if f < 0 00068 inline real eatanhe(real x) const { 00069 using std::atan; 00070 return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * atan(_e * x); 00071 } 00072 // Divided differences 00073 // Definition: Df(x,y) = (f(x)-f(y))/(x-y) 00074 // See: 00075 // W. M. Kahan and R. J. Fateman, 00076 // Symbolic computation of divided differences, 00077 // SIGSAM Bull. 33(3), 7-28 (1999) 00078 // http://dx.doi.org/10.1145/334714.334716 00079 // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf 00080 // 00081 // General rules 00082 // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y) 00083 // h(x) = f(x)*g(x): 00084 // Dh(x,y) = Df(x,y)*g(x) + Dg(x,y)*f(y) 00085 // = Df(x,y)*g(y) + Dg(x,y)*f(x) 00086 // = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2 00087 // 00088 // hyp(x) = sqrt(1+x^2): Dhyp(x,y) = (x+y)/(hyp(x)+hyp(y)) 00089 static inline real Dhyp(real x, real y, real hx, real hy) 00090 // hx = hyp(x) 00091 { return (x + y) / (hx + hy); } 00092 // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^2)) 00093 static inline real Dsn(real x, real y, real sx, real sy) { 00094 // sx = x/hyp(x) 00095 real t = x * y; 00096 return t > 0 ? (x + y) * Math::sq( (sx * sy)/t ) / (sx + sy) : 00097 (x - y != 0 ? (sx - sy) / (x - y) : 1); 00098 } 00099 // Dlog1p(x,y) = log1p((x-y)/(1+y)/(x-y) 00100 static inline real Dlog1p(real x, real y) { 00101 real t = x - y; if (t < 0) { t = -t; y = x; } 00102 return t ? Math::log1p(t / (1 + y)) / t : 1 / (1 + x); 00103 } 00104 // Dexp(x,y) = exp((x+y)/2) * 2*sinh((x-y)/2)/(x-y) 00105 static inline real Dexp(real x, real y) { 00106 using std::sinh; using std::exp; 00107 real t = (x - y)/2; 00108 return (t ? sinh(t)/t : 1) * exp((x + y)/2); 00109 } 00110 // Dsinh(x,y) = 2*sinh((x-y)/2)/(x-y) * cosh((x+y)/2) 00111 // cosh((x+y)/2) = (c+sinh(x)*sinh(y)/c)/2 00112 // c=sqrt((1+cosh(x))*(1+cosh(y))) 00113 // cosh((x+y)/2) = sqrt( (sinh(x)*sinh(y) + cosh(x)*cosh(y) + 1)/2 ) 00114 static inline real Dsinh(real x, real y, real sx, real sy, real cx, real cy) 00115 // sx = sinh(x), cx = cosh(x) 00116 { 00117 // real t = (x - y)/2, c = sqrt((1 + cx) * (1 + cy)); 00118 // return (t ? sinh(t)/t : real(1)) * (c + sx * sy / c) /2; 00119 using std::sinh; using std::sqrt; 00120 real t = (x - y)/2; 00121 return (t ? sinh(t)/t : 1) * sqrt((sx * sy + cx * cy + 1) /2); 00122 } 00123 // Dasinh(x,y) = asinh((x-y)*(x+y)/(x*sqrt(1+y^2)+y*sqrt(1+x^2)))/(x-y) 00124 // = asinh((x*sqrt(1+y^2)-y*sqrt(1+x^2)))/(x-y) 00125 static inline real Dasinh(real x, real y, real hx, real hy) { 00126 // hx = hyp(x) 00127 real t = x - y; 00128 return t ? 00129 Math::asinh(x*y > 0 ? t * (x+y) / (x*hy + y*hx) : x*hy - y*hx) / t : 00130 1/hx; 00131 } 00132 // Deatanhe(x,y) = eatanhe((x-y)/(1-e^2*x*y))/(x-y) 00133 inline real Deatanhe(real x, real y) const { 00134 real t = x - y, d = 1 - _e2 * x * y; 00135 return t ? eatanhe(t / d) / t : _e2 / d; 00136 } 00137 void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1); 00138 public: 00139 00140 /** 00141 * Constructor with a single standard parallel. 00142 * 00143 * @param[in] a equatorial radius of ellipsoid (meters). 00144 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00145 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set 00146 * flattening to 1/\e f. 00147 * @param[in] stdlat standard parallel (degrees), the circle of tangency. 00148 * @param[in] k0 scale on the standard parallel. 00149 * @exception GeographicErr if \e a, (1 − \e f) \e a, or \e k0 is 00150 * not positive. 00151 * @exception GeographicErr if \e stdlat is not in [−90°, 00152 * 90°]. 00153 **********************************************************************/ 00154 LambertConformalConic(real a, real f, real stdlat, real k0); 00155 00156 /** 00157 * Constructor with two standard parallels. 00158 * 00159 * @param[in] a equatorial radius of ellipsoid (meters). 00160 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00161 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set 00162 * flattening to 1/\e f. 00163 * @param[in] stdlat1 first standard parallel (degrees). 00164 * @param[in] stdlat2 second standard parallel (degrees). 00165 * @param[in] k1 scale on the standard parallels. 00166 * @exception GeographicErr if \e a, (1 − \e f) \e a, or \e k1 is 00167 * not positive. 00168 * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in 00169 * [−90°, 90°], or if either \e stdlat1 or \e 00170 * stdlat2 is a pole and \e stdlat1 is not equal \e stdlat2. 00171 **********************************************************************/ 00172 LambertConformalConic(real a, real f, real stdlat1, real stdlat2, real k1); 00173 00174 /** 00175 * Constructor with two standard parallels specified by sines and cosines. 00176 * 00177 * @param[in] a equatorial radius of ellipsoid (meters). 00178 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00179 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set 00180 * flattening to 1/\e f. 00181 * @param[in] sinlat1 sine of first standard parallel. 00182 * @param[in] coslat1 cosine of first standard parallel. 00183 * @param[in] sinlat2 sine of second standard parallel. 00184 * @param[in] coslat2 cosine of second standard parallel. 00185 * @param[in] k1 scale on the standard parallels. 00186 * @exception GeographicErr if \e a, (1 − \e f) \e a, or \e k1 is 00187 * not positive. 00188 * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in 00189 * [−90°, 90°], or if either \e stdlat1 or \e 00190 * stdlat2 is a pole and \e stdlat1 is not equal \e stdlat2. 00191 * 00192 * This allows parallels close to the poles to be specified accurately. 00193 * This routine computes the latitude of origin and the scale at this 00194 * latitude. In the case where \e lat1 and \e lat2 are different, the 00195 * errors in this routines are as follows: if \e dlat = abs(\e lat2 − 00196 * \e lat1) ≤ 160° and max(abs(\e lat1), abs(\e lat2)) ≤ 90 00197 * − min(0.0002, 2.2 × 10<sup>−6</sup>(180 − \e 00198 * dlat), 6 × 10<sup>−8</sup> <i>dlat</i><sup>2</sup>) (in 00199 * degrees), then the error in the latitude of origin is less than 4.5 00200 * × 10<sup>−14</sup>d and the relative error in the scale is 00201 * less than 7 × 10<sup>−15</sup>. 00202 **********************************************************************/ 00203 LambertConformalConic(real a, real f, 00204 real sinlat1, real coslat1, 00205 real sinlat2, real coslat2, 00206 real k1); 00207 00208 /** 00209 * Set the scale for the projection. 00210 * 00211 * @param[in] lat (degrees). 00212 * @param[in] k scale at latitude \e lat (default 1). 00213 * @exception GeographicErr \e k is not positive. 00214 * @exception GeographicErr if \e lat is not in [−90°, 00215 * 90°]. 00216 **********************************************************************/ 00217 void SetScale(real lat, real k = real(1)); 00218 00219 /** 00220 * Forward projection, from geographic to Lambert conformal conic. 00221 * 00222 * @param[in] lon0 central meridian longitude (degrees). 00223 * @param[in] lat latitude of point (degrees). 00224 * @param[in] lon longitude of point (degrees). 00225 * @param[out] x easting of point (meters). 00226 * @param[out] y northing of point (meters). 00227 * @param[out] gamma meridian convergence at point (degrees). 00228 * @param[out] k scale of projection at point. 00229 * 00230 * The latitude origin is given by LambertConformalConic::LatitudeOrigin(). 00231 * No false easting or northing is added and \e lat should be in the range 00232 * [−90°, 90°]; \e lon and \e lon0 should be in the 00233 * range [−540°, 540°). The error in the projection 00234 * is less than about 10 nm (10 nanometers), true distance, and the errors 00235 * in the meridian convergence and scale are consistent with this. The 00236 * values of \e x and \e y returned for points which project to infinity 00237 * (i.e., one or both of the poles) will be large but finite. 00238 **********************************************************************/ 00239 void Forward(real lon0, real lat, real lon, 00240 real& x, real& y, real& gamma, real& k) const; 00241 00242 /** 00243 * Reverse projection, from Lambert conformal conic to geographic. 00244 * 00245 * @param[in] lon0 central meridian longitude (degrees). 00246 * @param[in] x easting of point (meters). 00247 * @param[in] y northing of point (meters). 00248 * @param[out] lat latitude of point (degrees). 00249 * @param[out] lon longitude of point (degrees). 00250 * @param[out] gamma meridian convergence at point (degrees). 00251 * @param[out] k scale of projection at point. 00252 * 00253 * The latitude origin is given by LambertConformalConic::LatitudeOrigin(). 00254 * No false easting or northing is added. \e lon0 should be in the range 00255 * [−540°, 540°). The value of \e lon returned is in 00256 * the range [−180°, 180°). The error in the 00257 * projection is less than about 10 nm (10 nanometers), true distance, and 00258 * the errors in the meridian convergence and scale are consistent with 00259 * this. 00260 **********************************************************************/ 00261 void Reverse(real lon0, real x, real y, 00262 real& lat, real& lon, real& gamma, real& k) const; 00263 00264 /** 00265 * LambertConformalConic::Forward without returning the convergence and 00266 * scale. 00267 **********************************************************************/ 00268 void Forward(real lon0, real lat, real lon, 00269 real& x, real& y) const { 00270 real gamma, k; 00271 Forward(lon0, lat, lon, x, y, gamma, k); 00272 } 00273 00274 /** 00275 * LambertConformalConic::Reverse without returning the convergence and 00276 * scale. 00277 **********************************************************************/ 00278 void Reverse(real lon0, real x, real y, 00279 real& lat, real& lon) const { 00280 real gamma, k; 00281 Reverse(lon0, x, y, lat, lon, gamma, k); 00282 } 00283 00284 /** \name Inspector functions 00285 **********************************************************************/ 00286 ///@{ 00287 /** 00288 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00289 * the value used in the constructor. 00290 **********************************************************************/ 00291 Math::real MajorRadius() const { return _a; } 00292 00293 /** 00294 * @return \e f the flattening of the ellipsoid. This is the 00295 * value used in the constructor. 00296 **********************************************************************/ 00297 Math::real Flattening() const { return _f; } 00298 00299 /// \cond SKIP 00300 /** 00301 * <b>DEPRECATED</b> 00302 * @return \e r the inverse flattening of the ellipsoid. 00303 **********************************************************************/ 00304 Math::real InverseFlattening() const { return 1/_f; } 00305 /// \endcond 00306 00307 /** 00308 * @return latitude of the origin for the projection (degrees). 00309 * 00310 * This is the latitude of minimum scale and equals the \e stdlat in the 00311 * 1-parallel constructor and lies between \e stdlat1 and \e stdlat2 in the 00312 * 2-parallel constructors. 00313 **********************************************************************/ 00314 Math::real OriginLatitude() const { return _lat0; } 00315 00316 /** 00317 * @return central scale for the projection. This is the scale on the 00318 * latitude of origin. 00319 **********************************************************************/ 00320 Math::real CentralScale() const { return _k0; } 00321 ///@} 00322 00323 /** 00324 * A global instantiation of LambertConformalConic with the WGS84 00325 * ellipsoid, \e stdlat = 0, and \e k0 = 1. This degenerates to the 00326 * Mercator projection. 00327 **********************************************************************/ 00328 static const LambertConformalConic& Mercator(); 00329 }; 00330 00331 } // namespace GeographicLib 00332 00333 #endif // GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP