00001 /** 00002 * \file AlbersEqualArea.hpp 00003 * \brief Header for GeographicLib::AlbersEqualArea 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_ALBERSEQUALAREA_HPP) 00011 #define GEOGRAPHICLIB_ALBERSEQUALAREA_HPP 1 00012 00013 #include <GeographicLib/Constants.hpp> 00014 00015 namespace GeographicLib { 00016 00017 /** 00018 * \brief Albers equal area 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. 101--102. 00025 * 00026 * This is a implementation of the equations in Snyder except that divided 00027 * differences will be [have been] used to transform the expressions into 00028 * ones which may be evaluated accurately. [In this implementation, the 00029 * projection correctly becomes the cylindrical equal area or the azimuthal 00030 * equal area projection when the standard latitude is the equator or a 00031 * pole.] 00032 * 00033 * The ellipsoid parameters, the standard parallels, and the scale on the 00034 * standard parallels are set in the constructor. Internally, the case with 00035 * two standard parallels is converted into a single standard parallel, the 00036 * latitude of minimum azimuthal scale, with an azimuthal scale specified on 00037 * this parallel. This latitude is also used as the latitude of origin which 00038 * is returned by AlbersEqualArea::OriginLatitude. The azimuthal scale on 00039 * the latitude of origin is given by AlbersEqualArea::CentralScale. The 00040 * case with two standard parallels at opposite poles is singular and is 00041 * disallowed. The central meridian (which is a trivial shift of the 00042 * longitude) is specified as the \e lon0 argument of the 00043 * AlbersEqualArea::Forward and AlbersEqualArea::Reverse functions. 00044 * AlbersEqualArea::Forward and AlbersEqualArea::Reverse also return the 00045 * meridian convergence, γ, and azimuthal scale, \e k. A small square 00046 * aligned with the cardinal directions is projected to a rectangle with 00047 * dimensions \e k (in the E-W direction) and 1/\e k (in the N-S direction). 00048 * The E-W sides of the rectangle are oriented γ degrees 00049 * counter-clockwise from the \e x axis. There is no provision in this class 00050 * for specifying a false easting or false northing or a different latitude 00051 * of origin. 00052 * 00053 * Example of use: 00054 * \include example-AlbersEqualArea.cpp 00055 * 00056 * <a href="ConicProj.1.html">ConicProj</a> is a command-line utility 00057 * providing access to the functionality of LambertConformalConic and 00058 * AlbersEqualArea. 00059 **********************************************************************/ 00060 class GEOGRAPHICLIB_EXPORT AlbersEqualArea { 00061 private: 00062 typedef Math::real real; 00063 real eps_, epsx_, epsx2_, tol_, tol0_; 00064 real _a, _f, _fm, _e2, _e, _e2m, _qZ, _qx; 00065 real _sign, _lat0, _k0; 00066 real _n0, _m02, _nrho0, _k2, _txi0, _scxi0, _sxi0; 00067 static const int numit_ = 5; // Newton iterations in Reverse 00068 static const int numit0_ = 20; // Newton iterations in Init 00069 static inline real hyp(real x) { return Math::hypot(real(1), x); } 00070 // atanh( e * x)/ e if f > 0 00071 // atan (sqrt(-e2) * x)/sqrt(-e2) if f < 0 00072 // x if f = 0 00073 inline real atanhee(real x) const { 00074 using std::atan2; using std::abs; 00075 return _f > 0 ? Math::atanh(_e * x)/_e : 00076 // We only invoke atanhee in txif for positive latitude. Then x is 00077 // only negative for very prolate ellipsoids (_b/_a >= sqrt(2)) and we 00078 // still need to return a positive result in this case; hence the need 00079 // for the call to atan2. 00080 (_f < 0 ? (atan2(_e * abs(x), x < 0 ? -1 : 1)/_e) : x); 00081 } 00082 // return atanh(sqrt(x))/sqrt(x) - 1, accurate for small x 00083 static real atanhxm1(real x); 00084 00085 // Divided differences 00086 // Definition: Df(x,y) = (f(x)-f(y))/(x-y) 00087 // See: 00088 // W. M. Kahan and R. J. Fateman, 00089 // Symbolic computation of divided differences, 00090 // SIGSAM Bull. 33(3), 7-28 (1999) 00091 // http://dx.doi.org/10.1145/334714.334716 00092 // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf 00093 // 00094 // General rules 00095 // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y) 00096 // h(x) = f(x)*g(x): 00097 // Dh(x,y) = Df(x,y)*g(x) + Dg(x,y)*f(y) 00098 // = Df(x,y)*g(y) + Dg(x,y)*f(x) 00099 // = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2 00100 // 00101 // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^2)) 00102 static inline real Dsn(real x, real y, real sx, real sy) { 00103 // sx = x/hyp(x) 00104 real t = x * y; 00105 return t > 0 ? (x + y) * Math::sq( (sx * sy)/t ) / (sx + sy) : 00106 (x - y != 0 ? (sx - sy) / (x - y) : 1); 00107 } 00108 // Datanhee(x,y) = atanhee((x-y)/(1-e^2*x*y))/(x-y) 00109 inline real Datanhee(real x, real y) const { 00110 real t = x - y, d = 1 - _e2 * x * y; 00111 return t ? atanhee(t / d) / t : 1 / d; 00112 } 00113 // DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x) 00114 real DDatanhee(real x, real y) const; 00115 void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1); 00116 real txif(real tphi) const; 00117 real tphif(real txi) const; 00118 00119 friend class Ellipsoid; // For access to txif, tphif, etc. 00120 public: 00121 00122 /** 00123 * Constructor with a single standard parallel. 00124 * 00125 * @param[in] a equatorial radius of ellipsoid (meters). 00126 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00127 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set 00128 * flattening to 1/\e f. 00129 * @param[in] stdlat standard parallel (degrees), the circle of tangency. 00130 * @param[in] k0 azimuthal scale on the standard parallel. 00131 * @exception GeographicErr if \e a, (1 − \e f) \e a, or \e k0 is 00132 * not positive. 00133 * @exception GeographicErr if \e stdlat is not in [−90°, 00134 * 90°]. 00135 **********************************************************************/ 00136 AlbersEqualArea(real a, real f, real stdlat, real k0); 00137 00138 /** 00139 * Constructor with two standard parallels. 00140 * 00141 * @param[in] a equatorial radius of ellipsoid (meters). 00142 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00143 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set 00144 * flattening to 1/\e f. 00145 * @param[in] stdlat1 first standard parallel (degrees). 00146 * @param[in] stdlat2 second standard parallel (degrees). 00147 * @param[in] k1 azimuthal scale on the standard parallels. 00148 * @exception GeographicErr if \e a, (1 − \e f) \e a, or \e k1 is 00149 * not positive. 00150 * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in 00151 * [−90°, 90°], or if \e stdlat1 and \e stdlat2 are 00152 * opposite poles. 00153 **********************************************************************/ 00154 AlbersEqualArea(real a, real f, real stdlat1, real stdlat2, real k1); 00155 00156 /** 00157 * Constructor with two standard parallels specified by sines and cosines. 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] sinlat1 sine of first standard parallel. 00164 * @param[in] coslat1 cosine of first standard parallel. 00165 * @param[in] sinlat2 sine of second standard parallel. 00166 * @param[in] coslat2 cosine of second standard parallel. 00167 * @param[in] k1 azimuthal scale on the standard parallels. 00168 * @exception GeographicErr if \e a, (1 − \e f) \e a, or \e k1 is 00169 * not positive. 00170 * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in 00171 * [−90°, 90°], or if \e stdlat1 and \e stdlat2 are 00172 * opposite poles. 00173 * 00174 * This allows parallels close to the poles to be specified accurately. 00175 * This routine computes the latitude of origin and the azimuthal scale at 00176 * this latitude. If \e dlat = abs(\e lat2 − \e lat1) ≤ 160°, 00177 * then the error in the latitude of origin is less than 4.5 × 00178 * 10<sup>−14</sup>d;. 00179 **********************************************************************/ 00180 AlbersEqualArea(real a, real f, 00181 real sinlat1, real coslat1, 00182 real sinlat2, real coslat2, 00183 real k1); 00184 00185 /** 00186 * Set the azimuthal scale for the projection. 00187 * 00188 * @param[in] lat (degrees). 00189 * @param[in] k azimuthal scale at latitude \e lat (default 1). 00190 * @exception GeographicErr \e k is not positive. 00191 * @exception GeographicErr if \e lat is not in (−90°, 00192 * 90°). 00193 * 00194 * This allows a "latitude of conformality" to be specified. 00195 **********************************************************************/ 00196 void SetScale(real lat, real k = real(1)); 00197 00198 /** 00199 * Forward projection, from geographic to Lambert conformal conic. 00200 * 00201 * @param[in] lon0 central meridian longitude (degrees). 00202 * @param[in] lat latitude of point (degrees). 00203 * @param[in] lon longitude of point (degrees). 00204 * @param[out] x easting of point (meters). 00205 * @param[out] y northing of point (meters). 00206 * @param[out] gamma meridian convergence at point (degrees). 00207 * @param[out] k azimuthal scale of projection at point; the radial 00208 * scale is the 1/\e k. 00209 * 00210 * The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No 00211 * false easting or northing is added and \e lat should be in the range 00212 * [−90°, 90°]; \e lon and \e lon0 should be in the 00213 * range [−540°, 540°). The values of \e x and \e y 00214 * returned for points which project to infinity (i.e., one or both of the 00215 * poles) will be large but finite. 00216 **********************************************************************/ 00217 void Forward(real lon0, real lat, real lon, 00218 real& x, real& y, real& gamma, real& k) const; 00219 00220 /** 00221 * Reverse projection, from Lambert conformal conic to geographic. 00222 * 00223 * @param[in] lon0 central meridian longitude (degrees). 00224 * @param[in] x easting of point (meters). 00225 * @param[in] y northing of point (meters). 00226 * @param[out] lat latitude of point (degrees). 00227 * @param[out] lon longitude of point (degrees). 00228 * @param[out] gamma meridian convergence at point (degrees). 00229 * @param[out] k azimuthal scale of projection at point; the radial 00230 * scale is the 1/\e k. 00231 * 00232 * The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No 00233 * false easting or northing is added. \e lon0 should be in the range 00234 * [−540°, 540°). The value of \e lon returned is in 00235 * the range [−180°, 180°). The value of \e lat 00236 * returned is in the range [−90°, 90°]. If the 00237 * input point is outside the legal projected space the nearest pole is 00238 * returned. 00239 **********************************************************************/ 00240 void Reverse(real lon0, real x, real y, 00241 real& lat, real& lon, real& gamma, real& k) const; 00242 00243 /** 00244 * AlbersEqualArea::Forward without returning the convergence and 00245 * scale. 00246 **********************************************************************/ 00247 void Forward(real lon0, real lat, real lon, 00248 real& x, real& y) const { 00249 real gamma, k; 00250 Forward(lon0, lat, lon, x, y, gamma, k); 00251 } 00252 00253 /** 00254 * AlbersEqualArea::Reverse without returning the convergence and 00255 * scale. 00256 **********************************************************************/ 00257 void Reverse(real lon0, real x, real y, 00258 real& lat, real& lon) const { 00259 real gamma, k; 00260 Reverse(lon0, x, y, lat, lon, gamma, k); 00261 } 00262 00263 /** \name Inspector functions 00264 **********************************************************************/ 00265 ///@{ 00266 /** 00267 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00268 * the value used in the constructor. 00269 **********************************************************************/ 00270 Math::real MajorRadius() const { return _a; } 00271 00272 /** 00273 * @return \e f the flattening of the ellipsoid. This is the value used in 00274 * the constructor. 00275 **********************************************************************/ 00276 Math::real Flattening() const { return _f; } 00277 00278 /// \cond SKIP 00279 /** 00280 * <b>DEPRECATED</b> 00281 * @return \e r the inverse flattening of the ellipsoid. 00282 **********************************************************************/ 00283 Math::real InverseFlattening() const { return 1/_f; } 00284 /// \endcond 00285 00286 /** 00287 * @return latitude of the origin for the projection (degrees). 00288 * 00289 * This is the latitude of minimum azimuthal scale and equals the \e stdlat 00290 * in the 1-parallel constructor and lies between \e stdlat1 and \e stdlat2 00291 * in the 2-parallel constructors. 00292 **********************************************************************/ 00293 Math::real OriginLatitude() const { return _lat0; } 00294 00295 /** 00296 * @return central scale for the projection. This is the azimuthal scale 00297 * on the latitude of origin. 00298 **********************************************************************/ 00299 Math::real CentralScale() const { return _k0; } 00300 ///@} 00301 00302 /** 00303 * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e 00304 * stdlat = 0, and \e k0 = 1. This degenerates to the cylindrical equal 00305 * area projection. 00306 **********************************************************************/ 00307 static const AlbersEqualArea& CylindricalEqualArea(); 00308 00309 /** 00310 * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e 00311 * stdlat = 90°, and \e k0 = 1. This degenerates to the 00312 * Lambert azimuthal equal area projection. 00313 **********************************************************************/ 00314 static const AlbersEqualArea& AzimuthalEqualAreaNorth(); 00315 00316 /** 00317 * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e 00318 * stdlat = −90°, and \e k0 = 1. This degenerates to the 00319 * Lambert azimuthal equal area projection. 00320 **********************************************************************/ 00321 static const AlbersEqualArea& AzimuthalEqualAreaSouth(); 00322 }; 00323 00324 } // namespace GeographicLib 00325 00326 #endif // GEOGRAPHICLIB_ALBERSEQUALAREA_HPP