GeographicLib  1.35
AlbersEqualArea.hpp
Go to the documentation of this file.
1 /**
2  * \file AlbersEqualArea.hpp
3  * \brief Header for GeographicLib::AlbersEqualArea class
4  *
5  * Copyright (c) Charles Karney (2010-2012) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
10 #if !defined(GEOGRAPHICLIB_ALBERSEQUALAREA_HPP)
11 #define GEOGRAPHICLIB_ALBERSEQUALAREA_HPP 1
12 
14 
15 namespace GeographicLib {
16 
17  /**
18  * \brief Albers equal area conic projection
19  *
20  * Implementation taken from the report,
21  * - J. P. Snyder,
22  * <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projections: A
23  * Working Manual</a>, USGS Professional Paper 1395 (1987),
24  * pp. 101--102.
25  *
26  * This is a implementation of the equations in Snyder except that divided
27  * differences will be [have been] used to transform the expressions into
28  * ones which may be evaluated accurately. [In this implementation, the
29  * projection correctly becomes the cylindrical equal area or the azimuthal
30  * equal area projection when the standard latitude is the equator or a
31  * pole.]
32  *
33  * The ellipsoid parameters, the standard parallels, and the scale on the
34  * standard parallels are set in the constructor. Internally, the case with
35  * two standard parallels is converted into a single standard parallel, the
36  * latitude of minimum azimuthal scale, with an azimuthal scale specified on
37  * this parallel. This latitude is also used as the latitude of origin which
38  * is returned by AlbersEqualArea::OriginLatitude. The azimuthal scale on
39  * the latitude of origin is given by AlbersEqualArea::CentralScale. The
40  * case with two standard parallels at opposite poles is singular and is
41  * disallowed. The central meridian (which is a trivial shift of the
42  * longitude) is specified as the \e lon0 argument of the
43  * AlbersEqualArea::Forward and AlbersEqualArea::Reverse functions.
44  * AlbersEqualArea::Forward and AlbersEqualArea::Reverse also return the
45  * meridian convergence, &gamma;, and azimuthal scale, \e k. A small square
46  * aligned with the cardinal directions is projected to a rectangle with
47  * dimensions \e k (in the E-W direction) and 1/\e k (in the N-S direction).
48  * The E-W sides of the rectangle are oriented &gamma; degrees
49  * counter-clockwise from the \e x axis. There is no provision in this class
50  * for specifying a false easting or false northing or a different latitude
51  * of origin.
52  *
53  * Example of use:
54  * \include example-AlbersEqualArea.cpp
55  *
56  * <a href="ConicProj.1.html">ConicProj</a> is a command-line utility
57  * providing access to the functionality of LambertConformalConic and
58  * AlbersEqualArea.
59  **********************************************************************/
61  private:
62  typedef Math::real real;
63  real _a, _f, _fm, _e2, _e, _e2m, _qZ, _qx;
64  real _sign, _lat0, _k0;
65  real _n0, _m02, _nrho0, _k2, _txi0, _scxi0, _sxi0;
66  static const real eps_;
67  static const real epsx_;
68  static const real epsx2_;
69  static const real tol_;
70  static const real tol0_;
71  static const real ahypover_;
72  static const int numit_ = 5; // Newton iterations in Reverse
73  static const int numit0_ = 20; // Newton iterations in Init
74  static inline real hyp(real x) throw() { return Math::hypot(real(1), x); }
75  // atanh( e * x)/ e if f > 0
76  // atan (sqrt(-e2) * x)/sqrt(-e2) if f < 0
77  // x if f = 0
78  inline real atanhee(real x) const throw() {
79  return _f > 0 ? Math::atanh(_e * x)/_e :
80  // We only invoke atanhee in txif for positive latitude. Then x is
81  // only negative for very prolate ellipsoids (_b/_a >= sqrt(2)) and we
82  // still need to return a positive result in this case; hence the need
83  // for the call to atan2.
84  (_f < 0 ? (std::atan2(_e * std::abs(x), x < 0 ? -1 : 1)/_e) : x);
85  }
86  // return atanh(sqrt(x))/sqrt(x) - 1, accurate for small x
87  static real atanhxm1(real x) throw();
88 
89  // Divided differences
90  // Definition: Df(x,y) = (f(x)-f(y))/(x-y)
91  // See:
92  // W. M. Kahan and R. J. Fateman,
93  // Symbolic computation of divided differences,
94  // SIGSAM Bull. 33(3), 7-28 (1999)
95  // http://dx.doi.org/10.1145/334714.334716
96  // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf
97  //
98  // General rules
99  // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y)
100  // h(x) = f(x)*g(x):
101  // Dh(x,y) = Df(x,y)*g(x) + Dg(x,y)*f(y)
102  // = Df(x,y)*g(y) + Dg(x,y)*f(x)
103  // = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2
104  //
105  // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^2))
106  static inline real Dsn(real x, real y, real sx, real sy) throw() {
107  // sx = x/hyp(x)
108  real t = x * y;
109  return t > 0 ? (x + y) * Math::sq( (sx * sy)/t ) / (sx + sy) :
110  (x - y != 0 ? (sx - sy) / (x - y) : 1);
111  }
112  // Datanhee(x,y) = atanhee((x-y)/(1-e^2*x*y))/(x-y)
113  inline real Datanhee(real x, real y) const throw() {
114  real t = x - y, d = 1 - _e2 * x * y;
115  return t != 0 ? atanhee(t / d) / t : 1 / d;
116  }
117  // DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x)
118  real DDatanhee(real x, real y) const throw();
119  void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1) throw();
120  real txif(real tphi) const throw();
121  real tphif(real txi) const throw();
122 
123  friend class Ellipsoid; // For access to txif, tphif, etc.
124  public:
125 
126  /**
127  * Constructor with a single standard parallel.
128  *
129  * @param[in] a equatorial radius of ellipsoid (meters).
130  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
131  * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flattening
132  * to 1/\e f.
133  * @param[in] stdlat standard parallel (degrees), the circle of tangency.
134  * @param[in] k0 azimuthal scale on the standard parallel.
135  * @exception GeographicErr if \e a, (1 &minus; \e f ) \e a, or \e k0 is
136  * not positive.
137  * @exception GeographicErr if \e stdlat is not in [&minus;90&deg;,
138  * 90&deg;].
139  **********************************************************************/
140  AlbersEqualArea(real a, real f, real stdlat, real k0);
141 
142  /**
143  * Constructor with two standard parallels.
144  *
145  * @param[in] a equatorial radius of ellipsoid (meters).
146  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
147  * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flattening
148  * to 1/\e f.
149  * @param[in] stdlat1 first standard parallel (degrees).
150  * @param[in] stdlat2 second standard parallel (degrees).
151  * @param[in] k1 azimuthal scale on the standard parallels.
152  * @exception GeographicErr if \e a, (1 &minus; \e f ) \e a, or \e k1 is
153  * not positive.
154  * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in
155  * [&minus;90&deg;, 90&deg;], or if \e stdlat1 and \e stdlat2 are
156  * opposite poles.
157  **********************************************************************/
158  AlbersEqualArea(real a, real f, real stdlat1, real stdlat2, real k1);
159 
160  /**
161  * Constructor with two standard parallels specified by sines and cosines.
162  *
163  * @param[in] a equatorial radius of ellipsoid (meters).
164  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
165  * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flattening
166  * to 1/\e f.
167  * @param[in] sinlat1 sine of first standard parallel.
168  * @param[in] coslat1 cosine of first standard parallel.
169  * @param[in] sinlat2 sine of second standard parallel.
170  * @param[in] coslat2 cosine of second standard parallel.
171  * @param[in] k1 azimuthal scale on the standard parallels.
172  * @exception GeographicErr if \e a, (1 &minus; \e f ) \e a, or \e k1 is
173  * not positive.
174  * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in
175  * [&minus;90&deg;, 90&deg;], or if \e stdlat1 and \e stdlat2 are
176  * opposite poles.
177  *
178  * This allows parallels close to the poles to be specified accurately.
179  * This routine computes the latitude of origin and the azimuthal scale at
180  * this latitude. If \e dlat = abs(\e lat2 &minus; \e lat1) &le; 160&deg;,
181  * then the error in the latitude of origin is less than 4.5 &times;
182  * 10<sup>&minus;14</sup>d;.
183  **********************************************************************/
184  AlbersEqualArea(real a, real f,
185  real sinlat1, real coslat1,
186  real sinlat2, real coslat2,
187  real k1);
188 
189  /**
190  * Set the azimuthal scale for the projection.
191  *
192  * @param[in] lat (degrees).
193  * @param[in] k azimuthal scale at latitude \e lat (default 1).
194  * @exception GeographicErr \e k is not positive.
195  * @exception GeographicErr if \e lat is not in (&minus;90&deg;,
196  * 90&deg;).
197  *
198  * This allows a "latitude of conformality" to be specified.
199  **********************************************************************/
200  void SetScale(real lat, real k = real(1));
201 
202  /**
203  * Forward projection, from geographic to Lambert conformal conic.
204  *
205  * @param[in] lon0 central meridian longitude (degrees).
206  * @param[in] lat latitude of point (degrees).
207  * @param[in] lon longitude of point (degrees).
208  * @param[out] x easting of point (meters).
209  * @param[out] y northing of point (meters).
210  * @param[out] gamma meridian convergence at point (degrees).
211  * @param[out] k azimuthal scale of projection at point; the radial
212  * scale is the 1/\e k.
213  *
214  * The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No
215  * false easting or northing is added and \e lat should be in the range
216  * [&minus;90&deg;, 90&deg;]; \e lon and \e lon0 should be in the
217  * range [&minus;540&deg;, 540&deg;). The values of \e x and \e y
218  * returned for points which project to infinity (i.e., one or both of the
219  * poles) will be large but finite.
220  **********************************************************************/
221  void Forward(real lon0, real lat, real lon,
222  real& x, real& y, real& gamma, real& k) const throw();
223 
224  /**
225  * Reverse projection, from Lambert conformal conic to geographic.
226  *
227  * @param[in] lon0 central meridian longitude (degrees).
228  * @param[in] x easting of point (meters).
229  * @param[in] y northing of point (meters).
230  * @param[out] lat latitude of point (degrees).
231  * @param[out] lon longitude of point (degrees).
232  * @param[out] gamma meridian convergence at point (degrees).
233  * @param[out] k azimuthal scale of projection at point; the radial
234  * scale is the 1/\e k.
235  *
236  * The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No
237  * false easting or northing is added. \e lon0 should be in the range
238  * [&minus;540&deg;, 540&deg;). The value of \e lon returned is in
239  * the range [&minus;180&deg;, 180&deg;). The value of \e lat
240  * returned is in the range [&minus;90&deg;, 90&deg;]. If the
241  * input point is outside the legal projected space the nearest pole is
242  * returned.
243  **********************************************************************/
244  void Reverse(real lon0, real x, real y,
245  real& lat, real& lon, real& gamma, real& k) const throw();
246 
247  /**
248  * AlbersEqualArea::Forward without returning the convergence and
249  * scale.
250  **********************************************************************/
251  void Forward(real lon0, real lat, real lon,
252  real& x, real& y) const throw() {
253  real gamma, k;
254  Forward(lon0, lat, lon, x, y, gamma, k);
255  }
256 
257  /**
258  * AlbersEqualArea::Reverse without returning the convergence and
259  * scale.
260  **********************************************************************/
261  void Reverse(real lon0, real x, real y,
262  real& lat, real& lon) const throw() {
263  real gamma, k;
264  Reverse(lon0, x, y, lat, lon, gamma, k);
265  }
266 
267  /** \name Inspector functions
268  **********************************************************************/
269  ///@{
270  /**
271  * @return \e a the equatorial radius of the ellipsoid (meters). This is
272  * the value used in the constructor.
273  **********************************************************************/
274  Math::real MajorRadius() const throw() { return _a; }
275 
276  /**
277  * @return \e f the flattening of the ellipsoid. This is the value used in
278  * the constructor.
279  **********************************************************************/
280  Math::real Flattening() const throw() { return _f; }
281 
282  /// \cond SKIP
283  /**
284  * <b>DEPRECATED</b>
285  * @return \e r the inverse flattening of the ellipsoid.
286  **********************************************************************/
287  Math::real InverseFlattening() const throw() { return 1/_f; }
288  /// \endcond
289 
290  /**
291  * @return latitude of the origin for the projection (degrees).
292  *
293  * This is the latitude of minimum azimuthal scale and equals the \e stdlat
294  * in the 1-parallel constructor and lies between \e stdlat1 and \e stdlat2
295  * in the 2-parallel constructors.
296  **********************************************************************/
297  Math::real OriginLatitude() const throw() { return _lat0; }
298 
299  /**
300  * @return central scale for the projection. This is the azimuthal scale
301  * on the latitude of origin.
302  **********************************************************************/
303  Math::real CentralScale() const throw() { return _k0; }
304  ///@}
305 
306  /**
307  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
308  * stdlat = 0, and \e k0 = 1. This degenerates to the cylindrical equal
309  * area projection.
310  **********************************************************************/
312 
313  /**
314  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
315  * stdlat = 90&deg;, and \e k0 = 1. This degenerates to the
316  * Lambert azimuthal equal area projection.
317  **********************************************************************/
319 
320  /**
321  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
322  * stdlat = &minus;90&deg;, and \e k0 = 1. This degenerates to the
323  * Lambert azimuthal equal area projection.
324  **********************************************************************/
326  };
327 
328 } // namespace GeographicLib
329 
330 #endif // GEOGRAPHICLIB_ALBERSEQUALAREA_HPP
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:52
static const AlbersEqualArea AzimuthalEqualAreaSouth
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
static const AlbersEqualArea CylindricalEqualArea
static T atanh(T x)
Definition: Math.hpp:315
static const AlbersEqualArea AzimuthalEqualAreaNorth
void Reverse(real lon0, real x, real y, real &lat, real &lon) const
Albers equal area conic projection.
static T hypot(T x, T y)
Definition: Math.hpp:165
static T sq(T x)
Definition: Math.hpp:153
Properties of an ellipsoid.
Definition: Ellipsoid.hpp:39
Header for GeographicLib::Constants class.
Math::real OriginLatitude() const
void Forward(real lon0, real lat, real lon, real &x, real &y) const