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