GeographicLib  2.1.2
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-2022) <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 real hyp(real x) {
70  using std::hypot;
71  return hypot(real(1), x);
72  }
73  // atanh( e * x)/ e if f > 0
74  // atan (sqrt(-e2) * x)/sqrt(-e2) if f < 0
75  // x if f = 0
76  real atanhee(real x) const {
77  using std::atan; using std::atanh;
78  return _f > 0 ? atanh(_e * x)/_e : (_f < 0 ? (atan(_e * x)/_e) : x);
79  }
80  // return atanh(sqrt(x))/sqrt(x) - 1, accurate for small x
81  static real atanhxm1(real x);
82 
83  // Divided differences
84  // Definition: Df(x,y) = (f(x)-f(y))/(x-y)
85  // See:
86  // W. M. Kahan and R. J. Fateman,
87  // Symbolic computation of divided differences,
88  // SIGSAM Bull. 33(2), 7-28 (1999)
89  // https://doi.org/10.1145/334714.334716
90  // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf
91  //
92  // General rules
93  // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y)
94  // h(x) = f(x)*g(x):
95  // Dh(x,y) = Df(x,y)*g(x) + Dg(x,y)*f(y)
96  // = Df(x,y)*g(y) + Dg(x,y)*f(x)
97  // = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2
98  //
99  // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^2))
100  static real Dsn(real x, real y, real sx, real sy) {
101  // sx = x/hyp(x)
102  real t = x * y;
103  return t > 0 ? (x + y) * Math::sq( (sx * sy)/t ) / (sx + sy) :
104  (x - y != 0 ? (sx - sy) / (x - y) : 1);
105  }
106  // Datanhee(x,y) = (atanee(x)-atanee(y))/(x-y)
107  // = atanhee((x-y)/(1-e^2*x*y))/(x-y)
108  real Datanhee(real x, real y) const {
109  real t = x - y, d = 1 - _e2 * x * y;
110  return t == 0 ? 1 / d :
111  (x*y < 0 ? atanhee(x) - atanhee(y) : atanhee(t / d)) / t;
112  }
113  // DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x)
114  real DDatanhee(real x, real y) const;
115  real DDatanhee0(real x, real y) const;
116  real DDatanhee1(real x, real y) const;
117  real DDatanhee2(real x, real y) const;
118  void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1);
119  real txif(real tphi) const;
120  real tphif(real txi) const;
121 
122  friend class Ellipsoid; // For access to txif, tphif, etc.
123  public:
124 
125  /**
126  * Constructor with a single standard parallel.
127  *
128  * @param[in] a equatorial radius of ellipsoid (meters).
129  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
130  * Negative \e f gives a prolate ellipsoid.
131  * @param[in] stdlat standard parallel (degrees), the circle of tangency.
132  * @param[in] k0 azimuthal scale on the standard parallel.
133  * @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k0 is
134  * not positive.
135  * @exception GeographicErr if \e stdlat is not in [&minus;90&deg;,
136  * 90&deg;].
137  **********************************************************************/
138  AlbersEqualArea(real a, real f, real stdlat, real k0);
139 
140  /**
141  * Constructor with two standard parallels.
142  *
143  * @param[in] a equatorial radius of ellipsoid (meters).
144  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
145  * Negative \e f gives a prolate ellipsoid.
146  * @param[in] stdlat1 first standard parallel (degrees).
147  * @param[in] stdlat2 second standard parallel (degrees).
148  * @param[in] k1 azimuthal scale on the standard parallels.
149  * @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k1 is
150  * not positive.
151  * @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in
152  * [&minus;90&deg;, 90&deg;], or if \e stdlat1 and \e stdlat2 are
153  * opposite poles.
154  **********************************************************************/
155  AlbersEqualArea(real a, real f, real stdlat1, real stdlat2, real k1);
156 
157  /**
158  * Constructor with two standard parallels specified by sines and cosines.
159  *
160  * @param[in] a equatorial radius of ellipsoid (meters).
161  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
162  * Negative \e f gives a prolate ellipsoid.
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;]. The values of \e x and \e y returned for
213  * points which project to infinity (i.e., one or both of the poles) will
214  * be large but finite.
215  **********************************************************************/
216  void Forward(real lon0, real lat, real lon,
217  real& x, real& y, real& gamma, real& k) const;
218 
219  /**
220  * Reverse projection, from Lambert conformal conic to geographic.
221  *
222  * @param[in] lon0 central meridian longitude (degrees).
223  * @param[in] x easting of point (meters).
224  * @param[in] y northing of point (meters).
225  * @param[out] lat latitude of point (degrees).
226  * @param[out] lon longitude of point (degrees).
227  * @param[out] gamma meridian convergence at point (degrees).
228  * @param[out] k azimuthal scale of projection at point; the radial
229  * scale is the 1/\e k.
230  *
231  * The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No
232  * false easting or northing is added. The value of \e lon returned is in
233  * the range [&minus;180&deg;, 180&deg;]. The value of \e lat returned is
234  * in the range [&minus;90&deg;, 90&deg;]. If the input point is outside
235  * the legal projected space the nearest pole is returned.
236  **********************************************************************/
237  void Reverse(real lon0, real x, real y,
238  real& lat, real& lon, real& gamma, real& k) const;
239 
240  /**
241  * AlbersEqualArea::Forward without returning the convergence and
242  * scale.
243  **********************************************************************/
244  void Forward(real lon0, real lat, real lon,
245  real& x, real& y) const {
246  real gamma, k;
247  Forward(lon0, lat, lon, x, y, gamma, k);
248  }
249 
250  /**
251  * AlbersEqualArea::Reverse without returning the convergence and
252  * scale.
253  **********************************************************************/
254  void Reverse(real lon0, real x, real y,
255  real& lat, real& lon) const {
256  real gamma, k;
257  Reverse(lon0, x, y, lat, lon, gamma, k);
258  }
259 
260  /** \name Inspector functions
261  **********************************************************************/
262  ///@{
263  /**
264  * @return \e a the equatorial radius of the ellipsoid (meters). This is
265  * the value used in the constructor.
266  **********************************************************************/
267  Math::real EquatorialRadius() const { return _a; }
268 
269  /**
270  * @return \e f the flattening of the ellipsoid. This is the value used in
271  * the constructor.
272  **********************************************************************/
273  Math::real Flattening() const { return _f; }
274 
275  /**
276  * @return latitude of the origin for the projection (degrees).
277  *
278  * This is the latitude of minimum azimuthal scale and equals the \e stdlat
279  * in the 1-parallel constructor and lies between \e stdlat1 and \e stdlat2
280  * in the 2-parallel constructors.
281  **********************************************************************/
282  Math::real OriginLatitude() const { return _lat0; }
283 
284  /**
285  * @return central scale for the projection. This is the azimuthal scale
286  * on the latitude of origin.
287  **********************************************************************/
288  Math::real CentralScale() const { return _k0; }
289  ///@}
290 
291  /**
292  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
293  * stdlat = 0, and \e k0 = 1. This degenerates to the cylindrical equal
294  * area projection.
295  **********************************************************************/
296  static const AlbersEqualArea& CylindricalEqualArea();
297 
298  /**
299  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
300  * stdlat = 90&deg;, and \e k0 = 1. This degenerates to the
301  * Lambert azimuthal equal area projection.
302  **********************************************************************/
303  static const AlbersEqualArea& AzimuthalEqualAreaNorth();
304 
305  /**
306  * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
307  * stdlat = &minus;90&deg;, and \e k0 = 1. This degenerates to the
308  * Lambert azimuthal equal area projection.
309  **********************************************************************/
310  static const AlbersEqualArea& AzimuthalEqualAreaSouth();
311  };
312 
313 } // namespace GeographicLib
314 
315 #endif // GEOGRAPHICLIB_ALBERSEQUALAREA_HPP
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:67
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
Albers equal area conic projection.
void Reverse(real lon0, real x, real y, real &lat, real &lon) const
void Forward(real lon0, real lat, real lon, real &x, real &y) const
Properties of an ellipsoid.
Definition: Ellipsoid.hpp:39
static T sq(T x)
Definition: Math.hpp:212
Namespace for GeographicLib.
Definition: Accumulator.cpp:12