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