GeographicLib  1.43
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-2014) <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 eps_, epsx_, epsx2_, tol_, tol0_;
64  real _a, _f, _fm, _e2, _e, _e2m, _qZ, _qx;
65  real _sign, _lat0, _k0;
66  real _n0, _m02, _nrho0, _k2, _txi0, _scxi0, _sxi0;
67  static const int numit_ = 5; // Newton iterations in Reverse
68  static const int numit0_ = 20; // Newton iterations in Init
69  static inline real hyp(real x) { return Math::hypot(real(1), x); }
70  // atanh( e * x)/ e if f > 0
71  // atan (sqrt(-e2) * x)/sqrt(-e2) if f < 0
72  // x if f = 0
73  inline real atanhee(real x) const {
74  using std::atan2; using std::abs;
75  return _f > 0 ? Math::atanh(_e * x)/_e :
76  // We only invoke atanhee in txif for positive latitude. Then x is
77  // only negative for very prolate ellipsoids (_b/_a >= sqrt(2)) and we
78  // still need to return a positive result in this case; hence the need
79  // for the call to atan2.
80  (_f < 0 ? (atan2(_e * abs(x), x < 0 ? -1 : 1)/_e) : x);
81  }
82  // return atanh(sqrt(x))/sqrt(x) - 1, accurate for small x
83  static real atanhxm1(real x);
84 
85  // Divided differences
86  // Definition: Df(x,y) = (f(x)-f(y))/(x-y)
87  // See:
88  // W. M. Kahan and R. J. Fateman,
89  // Symbolic computation of divided differences,
90  // SIGSAM Bull. 33(3), 7-28 (1999)
91  // https://dx.doi.org/10.1145/334714.334716
92  // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf
93  //
94  // General rules
95  // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y)
96  // h(x) = f(x)*g(x):
97  // Dh(x,y) = Df(x,y)*g(x) + Dg(x,y)*f(y)
98  // = Df(x,y)*g(y) + Dg(x,y)*f(x)
99  // = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2
100  //
101  // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^2))
102  static inline real Dsn(real x, real y, real sx, real sy) {
103  // sx = x/hyp(x)
104  real t = x * y;
105  return t > 0 ? (x + y) * Math::sq( (sx * sy)/t ) / (sx + sy) :
106  (x - y != 0 ? (sx - sy) / (x - y) : 1);
107  }
108  // Datanhee(x,y) = atanhee((x-y)/(1-e^2*x*y))/(x-y)
109  inline real Datanhee(real x, real y) const {
110  real t = x - y, d = 1 - _e2 * x * y;
111  return t ? atanhee(t / d) / t : 1 / d;
112  }
113  // DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x)
114  real DDatanhee(real x, real y) const;
115  void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1);
116  real txif(real tphi) const;
117  real tphif(real txi) const;
118 
119  friend class Ellipsoid; // For access to txif, tphif, etc.
120  public:
121 
122  /**
123  * Constructor with a single standard parallel.
124  *
125  * @param[in] a equatorial radius of ellipsoid (meters).
126  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
127  * Negative \e f gives a prolate ellipsoid. If \e f &gt; 1, set
128  * flattening to 1/\e f.
129  * @param[in] stdlat standard parallel (degrees), the circle of tangency.
130  * @param[in] k0 azimuthal scale on the standard parallel.
131  * @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k0 is
132  * not positive.
133  * @exception GeographicErr if \e stdlat is not in [&minus;90&deg;,
134  * 90&deg;].
135  **********************************************************************/
136  AlbersEqualArea(real a, real f, real stdlat, real k0);
137 
138  /**
139  * Constructor with two standard parallels.
140  *
141  * @param[in] a equatorial radius of ellipsoid (meters).
142  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
143  * Negative \e f gives a prolate ellipsoid. If \e f &gt; 1, set
144  * flattening to 1/\e f.
145  * @param[in] stdlat1 first standard parallel (degrees).
146  * @param[in] stdlat2 second standard parallel (degrees).
147  * @param[in] k1 azimuthal scale on the standard parallels.
148  * @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k1 is
149  * not positive.
150  * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in
151  * [&minus;90&deg;, 90&deg;], or if \e stdlat1 and \e stdlat2 are
152  * opposite poles.
153  **********************************************************************/
154  AlbersEqualArea(real a, real f, real stdlat1, real stdlat2, real k1);
155 
156  /**
157  * Constructor with two standard parallels specified by sines and cosines.
158  *
159  * @param[in] a equatorial radius of ellipsoid (meters).
160  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
161  * Negative \e f gives a prolate ellipsoid. If \e f &gt; 1, set
162  * flattening to 1/\e f.
163  * @param[in] sinlat1 sine of first standard parallel.
164  * @param[in] coslat1 cosine of first standard parallel.
165  * @param[in] sinlat2 sine of second standard parallel.
166  * @param[in] coslat2 cosine of second standard parallel.
167  * @param[in] k1 azimuthal scale on the standard parallels.
168  * @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k1 is
169  * not positive.
170  * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in
171  * [&minus;90&deg;, 90&deg;], or if \e stdlat1 and \e stdlat2 are
172  * opposite poles.
173  *
174  * This allows parallels close to the poles to be specified accurately.
175  * This routine computes the latitude of origin and the azimuthal scale at
176  * this latitude. If \e dlat = abs(\e lat2 &minus; \e lat1) &le; 160&deg;,
177  * then the error in the latitude of origin is less than 4.5 &times;
178  * 10<sup>&minus;14</sup>d;.
179  **********************************************************************/
180  AlbersEqualArea(real a, real f,
181  real sinlat1, real coslat1,
182  real sinlat2, real coslat2,
183  real k1);
184 
185  /**
186  * Set the azimuthal scale for the projection.
187  *
188  * @param[in] lat (degrees).
189  * @param[in] k azimuthal scale at latitude \e lat (default 1).
190  * @exception GeographicErr \e k is not positive.
191  * @exception GeographicErr if \e lat is not in (&minus;90&deg;,
192  * 90&deg;).
193  *
194  * This allows a "latitude of conformality" to be specified.
195  **********************************************************************/
196  void SetScale(real lat, real k = real(1));
197 
198  /**
199  * Forward projection, from geographic to Lambert conformal conic.
200  *
201  * @param[in] lon0 central meridian longitude (degrees).
202  * @param[in] lat latitude of point (degrees).
203  * @param[in] lon longitude of point (degrees).
204  * @param[out] x easting of point (meters).
205  * @param[out] y northing of point (meters).
206  * @param[out] gamma meridian convergence at point (degrees).
207  * @param[out] k azimuthal scale of projection at point; the radial
208  * scale is the 1/\e k.
209  *
210  * The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No
211  * false easting or northing is added and \e lat should be in the range
212  * [&minus;90&deg;, 90&deg;]; \e lon and \e lon0 should be in the
213  * range [&minus;540&deg;, 540&deg;). The values of \e x and \e y
214  * returned for points which project to infinity (i.e., one or both of the
215  * poles) will be large but finite.
216  **********************************************************************/
217  void Forward(real lon0, real lat, real lon,
218  real& x, real& y, real& gamma, real& k) const;
219 
220  /**
221  * Reverse projection, from Lambert conformal conic to geographic.
222  *
223  * @param[in] lon0 central meridian longitude (degrees).
224  * @param[in] x easting of point (meters).
225  * @param[in] y northing of point (meters).
226  * @param[out] lat latitude of point (degrees).
227  * @param[out] lon longitude of point (degrees).
228  * @param[out] gamma meridian convergence at point (degrees).
229  * @param[out] k azimuthal scale of projection at point; the radial
230  * scale is the 1/\e k.
231  *
232  * The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No
233  * false easting or northing is added. \e lon0 should be in the range
234  * [&minus;540&deg;, 540&deg;). The value of \e lon returned is in
235  * the range [&minus;180&deg;, 180&deg;). The value of \e lat
236  * returned is in the range [&minus;90&deg;, 90&deg;]. If the
237  * input point is outside the legal projected space the nearest pole is
238  * returned.
239  **********************************************************************/
240  void Reverse(real lon0, real x, real y,
241  real& lat, real& lon, real& gamma, real& k) const;
242 
243  /**
244  * AlbersEqualArea::Forward without returning the convergence and
245  * scale.
246  **********************************************************************/
247  void Forward(real lon0, real lat, real lon,
248  real& x, real& y) const {
249  real gamma, k;
250  Forward(lon0, lat, lon, x, y, gamma, k);
251  }
252 
253  /**
254  * AlbersEqualArea::Reverse without returning the convergence and
255  * scale.
256  **********************************************************************/
257  void Reverse(real lon0, real x, real y,
258  real& lat, real& lon) const {
259  real gamma, k;
260  Reverse(lon0, x, y, lat, lon, gamma, k);
261  }
262 
263  /** \name Inspector functions
264  **********************************************************************/
265  ///@{
266  /**
267  * @return \e a the equatorial radius of the ellipsoid (meters). This is
268  * the value used in the constructor.
269  **********************************************************************/
270  Math::real MajorRadius() const { return _a; }
271 
272  /**
273  * @return \e f the flattening of the ellipsoid. This is the value used in
274  * the constructor.
275  **********************************************************************/
276  Math::real Flattening() const { return _f; }
277 
278  /// \cond SKIP
279  /**
280  * <b>DEPRECATED</b>
281  * @return \e r the inverse flattening of the ellipsoid.
282  **********************************************************************/
283  Math::real InverseFlattening() const { return 1/_f; }
284  /// \endcond
285 
286  /**
287  * @return latitude of the origin for the projection (degrees).
288  *
289  * This is the latitude of minimum azimuthal scale and equals the \e stdlat
290  * in the 1-parallel constructor and lies between \e stdlat1 and \e stdlat2
291  * in the 2-parallel constructors.
292  **********************************************************************/
293  Math::real OriginLatitude() const { return _lat0; }
294 
295  /**
296  * @return central scale for the projection. This is the azimuthal scale
297  * on the latitude of origin.
298  **********************************************************************/
299  Math::real CentralScale() const { return _k0; }
300  ///@}
301 
302  /**
303  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
304  * stdlat = 0, and \e k0 = 1. This degenerates to the cylindrical equal
305  * area projection.
306  **********************************************************************/
307  static const AlbersEqualArea& CylindricalEqualArea();
308 
309  /**
310  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
311  * stdlat = 90&deg;, and \e k0 = 1. This degenerates to the
312  * Lambert azimuthal equal area projection.
313  **********************************************************************/
314  static const AlbersEqualArea& AzimuthalEqualAreaNorth();
315 
316  /**
317  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
318  * stdlat = &minus;90&deg;, and \e k0 = 1. This degenerates to the
319  * Lambert azimuthal equal area projection.
320  **********************************************************************/
321  static const AlbersEqualArea& AzimuthalEqualAreaSouth();
322  };
323 
324 } // namespace GeographicLib
325 
326 #endif // GEOGRAPHICLIB_ALBERSEQUALAREA_HPP
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:90
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
static T atanh(T x)
Definition: Math.hpp:340
Albers equal area conic projection.
static T hypot(T x, T y)
Definition: Math.hpp:255
static T sq(T x)
Definition: Math.hpp:244
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
void Forward(real lon0, real lat, real lon, real &x, real &y) const
Properties of an ellipsoid.
Definition: Ellipsoid.hpp:39
Header for GeographicLib::Constants class.
void Reverse(real lon0, real x, real y, real &lat, real &lon) const