GeographicLib  1.48
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-2016) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
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://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.
128  * @param[in] stdlat standard parallel (degrees), the circle of tangency.
129  * @param[in] k0 azimuthal scale on the standard parallel.
130  * @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k0 is
131  * not positive.
132  * @exception GeographicErr if \e stdlat is not in [&minus;90&deg;,
133  * 90&deg;].
134  **********************************************************************/
135  AlbersEqualArea(real a, real f, real stdlat, real k0);
136 
137  /**
138  * Constructor with two standard parallels.
139  *
140  * @param[in] a equatorial radius of ellipsoid (meters).
141  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
142  * Negative \e f gives a prolate ellipsoid.
143  * @param[in] stdlat1 first standard parallel (degrees).
144  * @param[in] stdlat2 second standard parallel (degrees).
145  * @param[in] k1 azimuthal scale on the standard parallels.
146  * @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k1 is
147  * not positive.
148  * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in
149  * [&minus;90&deg;, 90&deg;], or if \e stdlat1 and \e stdlat2 are
150  * opposite poles.
151  **********************************************************************/
152  AlbersEqualArea(real a, real f, real stdlat1, real stdlat2, real k1);
153 
154  /**
155  * Constructor with two standard parallels specified by sines and cosines.
156  *
157  * @param[in] a equatorial radius of ellipsoid (meters).
158  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
159  * Negative \e f gives a prolate ellipsoid.
160  * @param[in] sinlat1 sine of first standard parallel.
161  * @param[in] coslat1 cosine of first standard parallel.
162  * @param[in] sinlat2 sine of second standard parallel.
163  * @param[in] coslat2 cosine of second standard parallel.
164  * @param[in] k1 azimuthal scale on the standard parallels.
165  * @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k1 is
166  * not positive.
167  * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in
168  * [&minus;90&deg;, 90&deg;], or if \e stdlat1 and \e stdlat2 are
169  * opposite poles.
170  *
171  * This allows parallels close to the poles to be specified accurately.
172  * This routine computes the latitude of origin and the azimuthal scale at
173  * this latitude. If \e dlat = abs(\e lat2 &minus; \e lat1) &le; 160&deg;,
174  * then the error in the latitude of origin is less than 4.5 &times;
175  * 10<sup>&minus;14</sup>d;.
176  **********************************************************************/
177  AlbersEqualArea(real a, real f,
178  real sinlat1, real coslat1,
179  real sinlat2, real coslat2,
180  real k1);
181 
182  /**
183  * Set the azimuthal scale for the projection.
184  *
185  * @param[in] lat (degrees).
186  * @param[in] k azimuthal scale at latitude \e lat (default 1).
187  * @exception GeographicErr \e k is not positive.
188  * @exception GeographicErr if \e lat is not in (&minus;90&deg;,
189  * 90&deg;).
190  *
191  * This allows a "latitude of conformality" to be specified.
192  **********************************************************************/
193  void SetScale(real lat, real k = real(1));
194 
195  /**
196  * Forward projection, from geographic to Lambert conformal conic.
197  *
198  * @param[in] lon0 central meridian longitude (degrees).
199  * @param[in] lat latitude of point (degrees).
200  * @param[in] lon longitude of point (degrees).
201  * @param[out] x easting of point (meters).
202  * @param[out] y northing of point (meters).
203  * @param[out] gamma meridian convergence at point (degrees).
204  * @param[out] k azimuthal scale of projection at point; the radial
205  * scale is the 1/\e k.
206  *
207  * The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No
208  * false easting or northing is added and \e lat should be in the range
209  * [&minus;90&deg;, 90&deg;]. The values of \e x and \e y returned for
210  * points which project to infinity (i.e., one or both of the poles) will
211  * be large but finite.
212  **********************************************************************/
213  void Forward(real lon0, real lat, real lon,
214  real& x, real& y, real& gamma, real& k) const;
215 
216  /**
217  * Reverse projection, from Lambert conformal conic to geographic.
218  *
219  * @param[in] lon0 central meridian longitude (degrees).
220  * @param[in] x easting of point (meters).
221  * @param[in] y northing of point (meters).
222  * @param[out] lat latitude of point (degrees).
223  * @param[out] lon longitude of point (degrees).
224  * @param[out] gamma meridian convergence at point (degrees).
225  * @param[out] k azimuthal scale of projection at point; the radial
226  * scale is the 1/\e k.
227  *
228  * The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No
229  * false easting or northing is added. The value of \e lon returned is in
230  * the range [&minus;180&deg;, 180&deg;]. The value of \e lat returned is
231  * in the range [&minus;90&deg;, 90&deg;]. If the input point is outside
232  * the legal projected space the nearest pole is returned.
233  **********************************************************************/
234  void Reverse(real lon0, real x, real y,
235  real& lat, real& lon, real& gamma, real& k) const;
236 
237  /**
238  * AlbersEqualArea::Forward without returning the convergence and
239  * scale.
240  **********************************************************************/
241  void Forward(real lon0, real lat, real lon,
242  real& x, real& y) const {
243  real gamma, k;
244  Forward(lon0, lat, lon, x, y, gamma, k);
245  }
246 
247  /**
248  * AlbersEqualArea::Reverse without returning the convergence and
249  * scale.
250  **********************************************************************/
251  void Reverse(real lon0, real x, real y,
252  real& lat, real& lon) const {
253  real gamma, k;
254  Reverse(lon0, x, y, lat, lon, gamma, k);
255  }
256 
257  /** \name Inspector functions
258  **********************************************************************/
259  ///@{
260  /**
261  * @return \e a the equatorial radius of the ellipsoid (meters). This is
262  * the value used in the constructor.
263  **********************************************************************/
264  Math::real MajorRadius() const { return _a; }
265 
266  /**
267  * @return \e f the flattening of the ellipsoid. This is the value used in
268  * the constructor.
269  **********************************************************************/
270  Math::real Flattening() const { return _f; }
271 
272  /**
273  * @return latitude of the origin for the projection (degrees).
274  *
275  * This is the latitude of minimum azimuthal scale and equals the \e stdlat
276  * in the 1-parallel constructor and lies between \e stdlat1 and \e stdlat2
277  * in the 2-parallel constructors.
278  **********************************************************************/
279  Math::real OriginLatitude() const { return _lat0; }
280 
281  /**
282  * @return central scale for the projection. This is the azimuthal scale
283  * on the latitude of origin.
284  **********************************************************************/
285  Math::real CentralScale() const { return _k0; }
286  ///@}
287 
288  /**
289  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
290  * stdlat = 0, and \e k0 = 1. This degenerates to the cylindrical equal
291  * area projection.
292  **********************************************************************/
293  static const AlbersEqualArea& CylindricalEqualArea();
294 
295  /**
296  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
297  * stdlat = 90&deg;, and \e k0 = 1. This degenerates to the
298  * Lambert azimuthal equal area projection.
299  **********************************************************************/
300  static const AlbersEqualArea& AzimuthalEqualAreaNorth();
301 
302  /**
303  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
304  * stdlat = &minus;90&deg;, and \e k0 = 1. This degenerates to the
305  * Lambert azimuthal equal area projection.
306  **********************************************************************/
307  static const AlbersEqualArea& AzimuthalEqualAreaSouth();
308  };
309 
310 } // namespace GeographicLib
311 
312 #endif // GEOGRAPHICLIB_ALBERSEQUALAREA_HPP
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:91
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
static T atanh(T x)
Definition: Math.hpp:328
Albers equal area conic projection.
static T hypot(T x, T y)
Definition: Math.hpp:243
static T sq(T x)
Definition: Math.hpp:232
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