AlbersEqualArea.hpp

Go to the documentation of this file.
00001 /**
00002  * \file AlbersEqualArea.hpp
00003  * \brief Header for GeographicLib::AlbersEqualArea class
00004  *
00005  * Copyright (c) Charles Karney (2010) <charles@karney.com> and licensed under
00006  * the LGPL.  For more information, see http://geographiclib.sourceforge.net/
00007  **********************************************************************/
00008 
00009 #if !defined(GEOGRAPHICLIB_ALBERSEQUALAREA_HPP)
00010 #define GEOGRAPHICLIB_ALBERSEQUALAREA_HPP "$Id: AlbersEqualArea.hpp 6919 2010-12-21 13:23:47Z karney $"
00011 
00012 #include "GeographicLib/Constants.hpp"
00013 #include <algorithm>
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&ndash;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, \e gamma, 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 \e gamma 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   class AlbersEqualArea {
00054   private:
00055     typedef Math::real real;
00056     const real _a, _r, _f, _fm, _e2, _e, _e2m, _qZ, _qx;
00057     real _sign, _lat0, _k0;
00058     real _n0, _m02, _nrho0, _k2, _txi0, _scxi0, _sxi0;
00059     static const real eps, epsx, epsx2, tol, tol0, ahypover;
00060     static const int numit = 5;   // Newton iterations in Reverse
00061     static const int numit0 = 20; // Newton iterations in Init
00062     static inline real sq(real x) throw() { return x * x; }
00063     static inline real hyp(real x) throw() { return Math::hypot(real(1), x); }
00064     // atanh(      e   * x)/      e   if f > 0
00065     // atan (sqrt(-e2) * x)/sqrt(-e2) if f < 0
00066     // x                              if f = 0
00067     inline real atanhee(real x) const throw() {
00068       return _f > 0 ? Math::atanh(_e * x)/_e :
00069         (_f < 0 ? std::atan(_e * x)/_e : x);
00070     }
00071     // return atanh(sqrt(x))/sqrt(x) - 1, accurate for small x
00072     static real atanhxm1(real x) throw();
00073 
00074     // Divided differences
00075     // Definition: Df(x,y) = (f(x)-f(y))/(x-y)
00076     // See: W. M. Kahan and R. J. Fateman,
00077     // Symbolic computation of divided differences,
00078     // SIGSAM Bull. 33(3), 7-28 (1999)
00079     // http://doi.acm.org/10.1145/334714.334716
00080     // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf
00081     //
00082     // General rules
00083     // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y)
00084     // h(x) = f(x)*g(x):
00085     //        Dh(x,y) = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2
00086     //
00087     // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^2))
00088     static inline real Dsn(real x, real y, real sx, real sy) throw() {
00089       // sx = x/hyp(x)
00090       real t = x * y;
00091       return t > 0 ? (x + y) * sq( (sx * sy)/t ) / (sx + sy) :
00092         (x - y != 0 ? (sx - sy) / (x - y) : 1);
00093     }
00094     // Datanhee(x,y) = atanhee((x-y)/(1-e^2*x*y))/(x-y)
00095     inline real Datanhee(real x, real y) const throw() {
00096       real t = x - y, d = 1 - _e2 * x * y;
00097       return t != 0 ? atanhee(t / d) / t : 1 / d;
00098     }
00099     // DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x)
00100     real DDatanhee(real x, real y) const throw();
00101     void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1) throw();
00102     real txif(real tphi) const throw();
00103     real tphif(real txi) const throw();
00104   public:
00105 
00106     /**
00107      * Constructor with a single standard parallel.
00108      *
00109      * @param[in] a equatorial radius of ellipsoid (meters)
00110      * @param[in] r reciprocal flattening of ellipsoid.  Setting \e r = 0
00111      *   implies \e r = inf or flattening = 0 (i.e., a sphere).  Negative \e r
00112      *   indicates a prolate ellipsoid.
00113      * @param[in] stdlat standard parallel (degrees), the circle of tangency.
00114      * @param[in] k0 azimuthal scale on the standard parallel.
00115      *
00116      * An exception is thrown if \e a or \e k0 is not positive or if \e stdlat
00117      * is not in the range [-90, 90].
00118      **********************************************************************/
00119     AlbersEqualArea(real a, real r, real stdlat, real k0);
00120 
00121     /**
00122      * Constructor with two standard parallels.
00123      *
00124      * @param[in] a equatorial radius of ellipsoid (meters)
00125      * @param[in] r reciprocal flattening of ellipsoid.  Setting \e r = 0
00126      *   implies \e r = inf or flattening = 0 (i.e., a sphere).  Negative \e r
00127      *   indicates a prolate ellipsoid.
00128      * @param[in] stdlat1 first standard parallel (degrees).
00129      * @param[in] stdlat2 second standard parallel (degrees).
00130      * @param[in] k1 azimuthal scale on the standard parallels.
00131      *
00132      * An exception is thrown if \e a or \e k0 is not positive or if \e stdlat1
00133      * or \e stdlat2 is not in the range [-90, 90].  In addition, an exception
00134      * is thrown if \e stdlat1 and \e stdlat2 are opposite poles.
00135      **********************************************************************/
00136     AlbersEqualArea(real a, real r, real stdlat1, real stdlat2, real k1);
00137 
00138     /**
00139      * Constructor with two standard parallels specified by sines and cosines.
00140      *
00141      * @param[in] a equatorial radius of ellipsoid (meters)
00142      * @param[in] r reciprocal flattening of ellipsoid.  Setting \e r = 0
00143      *   implies \e r = inf or flattening = 0 (i.e., a sphere).  Negative \e r
00144      *   indicates a prolate ellipsoid.
00145      * @param[in] sinlat1 sine of first standard parallel.
00146      * @param[in] coslat1 cosine of first standard parallel.
00147      * @param[in] sinlat2 sine of second standard parallel.
00148      * @param[in] coslat2 cosine of second standard parallel.
00149      * @param[in] k1 azimuthal scale on the standard parallels.
00150      *
00151      * This allows parallels close to the poles to be specified accurately.
00152      * This routine computes the latitude of origin and the azimuthal scale at
00153      * this latitude.  If \e dlat = abs(\e lat2 - \e lat1) <= 160<sup>o</sup>,
00154      * then the error in the latitude of origin is less than
00155      * 4.5e-14<sup>o</sup>.
00156      **********************************************************************/
00157     AlbersEqualArea(real a, real r,
00158                     real sinlat1, real coslat1,
00159                     real sinlat2, real coslat2,
00160                     real k1);
00161 
00162     /**
00163      * Set the azimuthal scale for the projection.
00164      *
00165      * @param[in] lat (degrees).
00166      * @param[in] k azimuthal scale at latitude \e lat (default 1).
00167      *
00168      * This allows a "latitude of conformality" to be specified.  An exception
00169      * is thrown if \e k is not positive or if \e lat is not in the range (-90,
00170      * 90).
00171      **********************************************************************/
00172     void SetScale(real lat, real k = real(1));
00173 
00174     /**
00175      * Forward projection, from geographic to Lambert conformal conic.
00176      *
00177      * @param[in] lon0 central meridian longitude (degrees).
00178      * @param[in] lat latitude of point (degrees).
00179      * @param[in] lon longitude of point (degrees).
00180      * @param[out] x easting of point (meters).
00181      * @param[out] y northing of point (meters).
00182      * @param[out] gamma meridian convergence at point (degrees).
00183      * @param[out] k azimuthal scale of projection at point; the radial
00184      *   scale is the 1/\e k.
00185      *
00186      * The latitude origin is given by AlbersEqualArea::LatitudeOrigin().  No
00187      * false easting or northing is added and \e lat should be in the range
00188      * [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360].  The
00189      * values of \e x and \e y returned for points which project to infinity
00190      * (i.e., one or both of the poles) will be large but finite.
00191      **********************************************************************/
00192     void Forward(real lon0, real lat, real lon,
00193                  real& x, real& y, real& gamma, real& k) const throw();
00194 
00195     /**
00196      * Reverse projection, from Lambert conformal conic to geographic.
00197      *
00198      * @param[in] lon0 central meridian longitude (degrees).
00199      * @param[in] x easting of point (meters).
00200      * @param[in] y northing of point (meters).
00201      * @param[out] lat latitude of point (degrees).
00202      * @param[out] lon longitude of point (degrees).
00203      * @param[out] gamma meridian convergence at point (degrees).
00204      * @param[out] k azimuthal scale of projection at point; the radial
00205      *   scale is the 1/\e k.
00206      *
00207      * The latitude origin is given by AlbersEqualArea::LatitudeOrigin().  No
00208      * false easting or northing is added.  \e lon0 should be in the range
00209      * [-180, 360].  The value of \e lon returned is in the range [-180, 180).
00210      * The value of \e lat returned is in the range [-90,90].  If the input
00211      * point is outside the legal projected space the nearest pole is returned.
00212      **********************************************************************/
00213     void Reverse(real lon0, real x, real y,
00214                  real& lat, real& lon, real& gamma, real& k) const throw();
00215 
00216     /**
00217      * AlbersEqualArea::Forward without returning the convergence and
00218      * scale.
00219      **********************************************************************/
00220     void Forward(real lon0, real lat, real lon,
00221                  real& x, real& y) const throw() {
00222       real gamma, k;
00223       Forward(lon0, lat, lon, x, y, gamma, k);
00224     }
00225 
00226     /**
00227      * AlbersEqualArea::Reverse without returning the convergence and
00228      * scale.
00229      **********************************************************************/
00230     void Reverse(real lon0, real x, real y,
00231                  real& lat, real& lon) const throw() {
00232       real gamma, k;
00233       Reverse(lon0, x, y, lat, lon, gamma, k);
00234     }
00235 
00236     /** \name Inspector functions
00237      **********************************************************************/
00238     ///@{
00239     /**
00240      * @return \e a the equatorial radius of the ellipsoid (meters).  This is
00241      *   the value used in the constructor.
00242      **********************************************************************/
00243     Math::real MajorRadius() const throw() { return _a; }
00244 
00245     /**
00246      * @return \e r the inverse flattening of the ellipsoid.  This is the
00247      *   value used in the constructor.  A value of 0 is returned for a sphere
00248      *   (infinite inverse flattening).
00249      **********************************************************************/
00250     Math::real InverseFlattening() const throw() { return _r; }
00251 
00252     /**
00253      * @return latitude of the origin for the projection (degrees).
00254      *
00255      * This is the latitude of minimum azimuthal scale and equals the \e stdlat
00256      * in the 1-parallel constructor and lies between \e stdlat1 and \e stdlat2
00257      * in the 2-parallel constructors.
00258      **********************************************************************/
00259     Math::real OriginLatitude() const throw() { return _lat0; }
00260 
00261     /**
00262      * @return central scale for the projection.  This is the azimuthal scale
00263      *   on the latitude of origin.
00264      **********************************************************************/
00265     Math::real CentralScale() const throw() { return _k0; }
00266     ///@}
00267 
00268     /**
00269      * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
00270      * stdlat = 0, and \e k0 = 1.  This degenerates to the cylindrical equal
00271      * area projection.
00272      **********************************************************************/
00273     static const AlbersEqualArea CylindricalEqualArea;
00274 
00275     /**
00276      * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
00277      * stdlat = 90<sup>o</sup>, and \e k0 = 1.  This degenerates to the
00278      * Lambert azimuthal equal area projection.
00279      **********************************************************************/
00280     static const AlbersEqualArea AzimuthalEqualAreaNorth;
00281 
00282     /**
00283      * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
00284      * stdlat = -90<sup>o</sup>, and \e k0 = 1.  This degenerates to the
00285      * Lambert azimuthal equal area projection.
00286      **********************************************************************/
00287     static const AlbersEqualArea AzimuthalEqualAreaSouth;
00288   };
00289 
00290 } // namespace GeographicLib
00291 
00292 #endif