GeographicLib  1.42
Geoid.hpp
Go to the documentation of this file.
1 /**
2  * \file Geoid.hpp
3  * \brief Header for GeographicLib::Geoid class
4  *
5  * Copyright (c) Charles Karney (2009-2015) <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_GEOID_HPP)
11 #define GEOGRAPHICLIB_GEOID_HPP 1
12 
13 #include <vector>
14 #include <fstream>
16 
17 #if defined(_MSC_VER)
18 // Squelch warnings about dll vs vector and constant conditional expressions
19 # pragma warning (push)
20 # pragma warning (disable: 4251 4127)
21 #endif
22 
23 #if !defined(GEOGRAPHICLIB_GEOID_PGM_PIXEL_WIDTH)
24 /**
25  * The size of the pixel data in the pgm data files for the geoids. 2 is the
26  * standard size corresponding to a maxval 2<sup>16</sup>&minus;1. Setting it
27  * to 4 uses a maxval of 2<sup>32</sup>&minus;1 and changes the extension for
28  * the data files from .pgm to .pgm4. Note that the format of these pgm4 files
29  * is a non-standard extension of the pgm format.
30  **********************************************************************/
31 # define GEOGRAPHICLIB_GEOID_PGM_PIXEL_WIDTH 2
32 #endif
33 
34 namespace GeographicLib {
35 
36  /**
37  * \brief Looking up the height of the geoid
38  *
39  * This class evaluated the height of one of the standard geoids, EGM84,
40  * EGM96, or EGM2008 by bilinear or cubic interpolation into a rectangular
41  * grid of data. These geoid models are documented in
42  * - EGM84:
43  * http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html
44  * - EGM96:
45  * http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html
46  * - EGM2008:
47  * http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008
48  *
49  * The geoids are defined in terms of spherical harmonics. However in order
50  * to provide a quick and flexible method of evaluating the geoid heights,
51  * this class evaluates the height by interpolation into a grid of
52  * precomputed values.
53  *
54  * The geoid height, \e N, can be used to convert a height above the
55  * ellipsoid, \e h, to the corresponding height above the geoid (roughly the
56  * height above mean sea level), \e H, using the relations
57  *
58  * &nbsp;&nbsp;&nbsp;\e h = \e N + \e H;
59  * &nbsp;&nbsp;\e H = &minus;\e N + \e h.
60  *
61  * See \ref geoid for details of how to install the data sets, the data
62  * format, estimates of the interpolation errors, and how to use caching.
63  *
64  * In addition to returning the geoid height, the gradient of the geoid can
65  * be calculated. The gradient is defined as the rate of change of the geoid
66  * as a function of position on the ellipsoid. This uses the parameters for
67  * the WGS84 ellipsoid. The gradient defined in terms of the interpolated
68  * heights. As a result of the way that the geoid data is stored, the
69  * calculation of gradients can result in large quantization errors. This is
70  * particularly acute for fine grids, at high latitudes, and for the easterly
71  * gradient. For this reason, the use of this facility is <b>DEPRECATED</b>.
72  * Instead, use the GravityModel class to evaluate the gravity vector.
73  *
74  * This class is typically \e not thread safe in that a single instantiation
75  * cannot be safely used by multiple threads because of the way the object
76  * reads the data set and because it maintains a single-cell cache. If
77  * multiple threads need to calculate geoid heights they should all construct
78  * thread-local instantiations. Alternatively, set the optional \e
79  * threadsafe parameter to true in the constructor. This causes the
80  * constructor to read all the data into memory and to turn off the
81  * single-cell caching which results in a Geoid object which \e is thread
82  * safe.
83  *
84  * Example of use:
85  * \include example-Geoid.cpp
86  *
87  * <a href="GeoidEval.1.html">GeoidEval</a> is a command-line utility
88  * providing access to the functionality of Geoid.
89  **********************************************************************/
90 
92  private:
93  typedef Math::real real;
94 #if GEOGRAPHICLIB_GEOID_PGM_PIXEL_WIDTH != 4
95  typedef unsigned short pixel_t;
96  static const unsigned pixel_size_ = 2;
97  static const unsigned pixel_max_ = 0xffffu;
98 #else
99  typedef unsigned pixel_t;
100  static const unsigned pixel_size_ = 4;
101  static const unsigned pixel_max_ = 0xffffffffu;
102 #endif
103  static const unsigned stencilsize_ = 12;
104  static const unsigned nterms_ = ((3 + 1) * (3 + 2))/2; // for a cubic fit
105  static const int c0_;
106  static const int c0n_;
107  static const int c0s_;
108  static const int c3_[stencilsize_ * nterms_];
109  static const int c3n_[stencilsize_ * nterms_];
110  static const int c3s_[stencilsize_ * nterms_];
111 
112  std::string _name, _dir, _filename;
113  const bool _cubic;
114  const real _a, _e2, _degree, _eps;
115  mutable std::ifstream _file;
116  real _rlonres, _rlatres;
117  std::string _description, _datetime;
118  real _offset, _scale, _maxerror, _rmserror;
119  int _width, _height;
120  unsigned long long _datastart, _swidth;
121  bool _threadsafe;
122  // Area cache
123  mutable std::vector< std::vector<pixel_t> > _data;
124  mutable bool _cache;
125  // NE corner and extent of cache
126  mutable int _xoffset, _yoffset, _xsize, _ysize;
127  // Cell cache
128  mutable int _ix, _iy;
129  mutable real _v00, _v01, _v10, _v11;
130  mutable real _t[nterms_];
131  void filepos(int ix, int iy) const {
132  _file.seekg(
133 #if !(defined(__GNUC__) && __GNUC__ < 4)
134  // g++ 3.x doesn't know about the cast to streamoff.
135  std::ios::streamoff
136 #endif
137  (_datastart +
138  pixel_size_ * (unsigned(iy)*_swidth + unsigned(ix))));
139  }
140  real rawval(int ix, int iy) const {
141  if (ix < 0)
142  ix += _width;
143  else if (ix >= _width)
144  ix -= _width;
145  if (_cache && iy >= _yoffset && iy < _yoffset + _ysize &&
146  ((ix >= _xoffset && ix < _xoffset + _xsize) ||
147  (ix + _width >= _xoffset && ix + _width < _xoffset + _xsize))) {
148  return real(_data[iy - _yoffset]
149  [ix >= _xoffset ? ix - _xoffset : ix + _width - _xoffset]);
150  } else {
151  if (iy < 0 || iy >= _height) {
152  iy = iy < 0 ? -iy : 2 * (_height - 1) - iy;
153  ix += (ix < _width/2 ? 1 : -1) * _width/2;
154  }
155  try {
156  filepos(ix, iy);
157  // initial values to suppress warnings in case get fails
158  char a = 0, b = 0;
159  _file.get(a);
160  _file.get(b);
161  unsigned r = ((unsigned char)(a) << 8) | (unsigned char)(b);
162  if (pixel_size_ == 4) {
163  _file.get(a);
164  _file.get(b);
165  r = (r << 16) | ((unsigned char)(a) << 8) | (unsigned char)(b);
166  }
167  return real(r);
168  }
169  catch (const std::exception& e) {
170  // throw GeographicErr("Error reading " + _filename + ": "
171  // + e.what());
172  // triggers complaints about the "binary '+'" under Visual Studio.
173  // So use '+=' instead.
174  std::string err("Error reading ");
175  err += _filename;
176  err += ": ";
177  err += e.what();
178  throw GeographicErr(err);
179  }
180  }
181  }
182  real height(real lat, real lon, bool gradp,
183  real& grade, real& gradn) const;
184  Geoid(const Geoid&); // copy constructor not allowed
185  Geoid& operator=(const Geoid&); // copy assignment not allowed
186  public:
187 
188  /**
189  * Flags indicating conversions between heights above the geoid and heights
190  * above the ellipsoid.
191  **********************************************************************/
192  enum convertflag {
193  /**
194  * The multiplier for converting from heights above the geoid to heights
195  * above the ellipsoid.
196  **********************************************************************/
197  ELLIPSOIDTOGEOID = -1,
198  /**
199  * No conversion.
200  **********************************************************************/
201  NONE = 0,
202  /**
203  * The multiplier for converting from heights above the ellipsoid to
204  * heights above the geoid.
205  **********************************************************************/
206  GEOIDTOELLIPSOID = 1,
207  };
208 
209  /** \name Setting up the geoid
210  **********************************************************************/
211  ///@{
212  /**
213  * Construct a geoid.
214  *
215  * @param[in] name the name of the geoid.
216  * @param[in] path (optional) directory for data file.
217  * @param[in] cubic (optional) interpolation method; false means bilinear,
218  * true (the default) means cubic.
219  * @param[in] threadsafe (optional), if true, construct a thread safe
220  * object. The default is false
221  * @exception GeographicErr if the data file cannot be found, is
222  * unreadable, or is corrupt.
223  * @exception GeographicErr if \e threadsafe is true but the memory
224  * necessary for caching the data can't be allocated.
225  *
226  * The data file is formed by appending ".pgm" to the name. If \e path is
227  * specified (and is non-empty), then the file is loaded from directory, \e
228  * path. Otherwise the path is given by DefaultGeoidPath(). If the \e
229  * threadsafe parameter is true, the data set is read into memory, the data
230  * file is closed, and single-cell caching is turned off; this results in a
231  * Geoid object which \e is thread safe.
232  **********************************************************************/
233  explicit Geoid(const std::string& name, const std::string& path = "",
234  bool cubic = true, bool threadsafe = false);
235 
236  /**
237  * Set up a cache.
238  *
239  * @param[in] south latitude (degrees) of the south edge of the cached area.
240  * @param[in] west longitude (degrees) of the west edge of the cached area.
241  * @param[in] north latitude (degrees) of the north edge of the cached area.
242  * @param[in] east longitude (degrees) of the east edge of the cached area.
243  * @exception GeographicErr if the memory necessary for caching the data
244  * can't be allocated (in this case, you will have no cache and can try
245  * again with a smaller area).
246  * @exception GeographicErr if there's a problem reading the data.
247  * @exception GeographicErr if this is called on a threadsafe Geoid.
248  *
249  * Cache the data for the specified "rectangular" area bounded by the
250  * parallels \e south and \e north and the meridians \e west and \e east.
251  * \e east is always interpreted as being east of \e west, if necessary by
252  * adding 360&deg; to its value. \e south and \e north should be in
253  * the range [&minus;90&deg;, 90&deg;]; \e west and \e east should
254  * be in the range [&minus;540&deg;, 540&deg;).
255  **********************************************************************/
256  void CacheArea(real south, real west, real north, real east) const;
257 
258  /**
259  * Cache all the data.
260  *
261  * @exception GeographicErr if the memory necessary for caching the data
262  * can't be allocated (in this case, you will have no cache and can try
263  * again with a smaller area).
264  * @exception GeographicErr if there's a problem reading the data.
265  * @exception GeographicErr if this is called on a threadsafe Geoid.
266  *
267  * On most computers, this is fast for data sets with grid resolution of 5'
268  * or coarser. For a 1' grid, the required RAM is 450MB; a 2.5' grid needs
269  * 72MB; and a 5' grid needs 18MB.
270  **********************************************************************/
271  void CacheAll() const { CacheArea(real(-90), real(0),
272  real(90), real(360)); }
273 
274  /**
275  * Clear the cache. This never throws an error. (This does nothing with a
276  * thread safe Geoid.)
277  **********************************************************************/
278  void CacheClear() const;
279 
280  ///@}
281 
282  /** \name Compute geoid heights
283  **********************************************************************/
284  ///@{
285  /**
286  * Compute the geoid height at a point
287  *
288  * @param[in] lat latitude of the point (degrees).
289  * @param[in] lon longitude of the point (degrees).
290  * @exception GeographicErr if there's a problem reading the data; this
291  * never happens if (\e lat, \e lon) is within a successfully cached area.
292  * @return geoid height (meters).
293  *
294  * The latitude should be in [&minus;90&deg;, 90&deg;] and
295  * longitude should be in [&minus;540&deg;, 540&deg;).
296  **********************************************************************/
297  Math::real operator()(real lat, real lon) const {
298  real gradn, grade;
299  return height(lat, lon, false, gradn, grade);
300  }
301 
302  /**
303  * Compute the geoid height and gradient at a point
304  *
305  * @param[in] lat latitude of the point (degrees).
306  * @param[in] lon longitude of the point (degrees).
307  * @param[out] gradn northerly gradient (dimensionless).
308  * @param[out] grade easterly gradient (dimensionless).
309  * @exception GeographicErr if there's a problem reading the data; this
310  * never happens if (\e lat, \e lon) is within a successfully cached area.
311  * @return geoid height (meters).
312  *
313  * The latitude should be in [&minus;90&deg;, 90&deg;] and longitude should
314  * be in [&minus;540&deg;, 540&deg;). As a result of the way that the
315  * geoid data is stored, the calculation of gradients can result in large
316  * quantization errors. This is particularly acute for fine grids, at high
317  * latitudes, and for the easterly gradient. For this reason, the
318  * computation of the gradient is <b>DEPRECATED</b>. If you need to
319  * compute the direction of the acceleration due to gravity accurately, you
320  * should use GravityModel::Gravity.
321  **********************************************************************/
322  Math::real operator()(real lat, real lon, real& gradn, real& grade) const {
323  return height(lat, lon, true, gradn, grade);
324  }
325 
326  /**
327  * Convert a height above the geoid to a height above the ellipsoid and
328  * vice versa.
329  *
330  * @param[in] lat latitude of the point (degrees).
331  * @param[in] lon longitude of the point (degrees).
332  * @param[in] h height of the point (degrees).
333  * @param[in] d a Geoid::convertflag specifying the direction of the
334  * conversion; Geoid::GEOIDTOELLIPSOID means convert a height above the
335  * geoid to a height above the ellipsoid; Geoid::ELLIPSOIDTOGEOID means
336  * convert a height above the ellipsoid to a height above the geoid.
337  * @exception GeographicErr if there's a problem reading the data; this
338  * never happens if (\e lat, \e lon) is within a successfully cached area.
339  * @return converted height (meters).
340  **********************************************************************/
341  Math::real ConvertHeight(real lat, real lon, real h,
342  convertflag d) const {
343  real gradn, grade;
344  return h + real(d) * height(lat, lon, true, gradn, grade);
345  }
346 
347  ///@}
348 
349  /** \name Inspector functions
350  **********************************************************************/
351  ///@{
352  /**
353  * @return geoid description, if available, in the data file; if
354  * absent, return "NONE".
355  **********************************************************************/
356  const std::string& Description() const { return _description; }
357 
358  /**
359  * @return date of the data file; if absent, return "UNKNOWN".
360  **********************************************************************/
361  const std::string& DateTime() const { return _datetime; }
362 
363  /**
364  * @return full file name used to load the geoid data.
365  **********************************************************************/
366  const std::string& GeoidFile() const { return _filename; }
367 
368  /**
369  * @return "name" used to load the geoid data (from the first argument of
370  * the constructor).
371  **********************************************************************/
372  const std::string& GeoidName() const { return _name; }
373 
374  /**
375  * @return directory used to load the geoid data.
376  **********************************************************************/
377  const std::string& GeoidDirectory() const { return _dir; }
378 
379  /**
380  * @return interpolation method ("cubic" or "bilinear").
381  **********************************************************************/
382  const std::string Interpolation() const
383  { return std::string(_cubic ? "cubic" : "bilinear"); }
384 
385  /**
386  * @return estimate of the maximum interpolation and quantization error
387  * (meters).
388  *
389  * This relies on the value being stored in the data file. If the value is
390  * absent, return &minus;1.
391  **********************************************************************/
392  Math::real MaxError() const { return _maxerror; }
393 
394  /**
395  * @return estimate of the RMS interpolation and quantization error
396  * (meters).
397  *
398  * This relies on the value being stored in the data file. If the value is
399  * absent, return &minus;1.
400  **********************************************************************/
401  Math::real RMSError() const { return _rmserror; }
402 
403  /**
404  * @return offset (meters).
405  *
406  * This in used in converting from the pixel values in the data file to
407  * geoid heights.
408  **********************************************************************/
409  Math::real Offset() const { return _offset; }
410 
411  /**
412  * @return scale (meters).
413  *
414  * This in used in converting from the pixel values in the data file to
415  * geoid heights.
416  **********************************************************************/
417  Math::real Scale() const { return _scale; }
418 
419  /**
420  * @return true if the object is constructed to be thread safe.
421  **********************************************************************/
422  bool ThreadSafe() const { return _threadsafe; }
423 
424  /**
425  * @return true if a data cache is active.
426  **********************************************************************/
427  bool Cache() const { return _cache; }
428 
429  /**
430  * @return west edge of the cached area; the cache includes this edge.
431  **********************************************************************/
433  return _cache ? ((_xoffset + (_xsize == _width ? 0 : _cubic)
434  + _width/2) % _width - _width/2) / _rlonres :
435  0;
436  }
437 
438  /**
439  * @return east edge of the cached area; the cache excludes this edge.
440  **********************************************************************/
442  return _cache ?
443  CacheWest() +
444  (_xsize - (_xsize == _width ? 0 : 1 + 2 * _cubic)) / _rlonres :
445  0;
446  }
447 
448  /**
449  * @return north edge of the cached area; the cache includes this edge.
450  **********************************************************************/
452  return _cache ? 90 - (_yoffset + _cubic) / _rlatres : 0;
453  }
454 
455  /**
456  * @return south edge of the cached area; the cache excludes this edge
457  * unless it's the south pole.
458  **********************************************************************/
460  return _cache ? 90 - ( _yoffset + _ysize - 1 - _cubic) / _rlatres : 0;
461  }
462 
463  /**
464  * @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
465  *
466  * (The WGS84 value is returned because the supported geoid models are all
467  * based on this ellipsoid.)
468  **********************************************************************/
470  { return Constants::WGS84_a(); }
471 
472  /**
473  * @return \e f the flattening of the WGS84 ellipsoid.
474  *
475  * (The WGS84 value is returned because the supported geoid models are all
476  * based on this ellipsoid.)
477  **********************************************************************/
479  ///@}
480 
481  /// \cond SKIP
482  /**
483  * <b>DEPRECATED</b>
484  * @return \e r the inverse flattening of the WGS84 ellipsoid.
485  **********************************************************************/
486  Math::real InverseFlattening() const
487  { return 1/Constants::WGS84_f(); }
488  /// \endcond
489 
490  /**
491  * @return the default path for geoid data files.
492  *
493  * This is the value of the environment variable GEOGRAPHICLIB_GEOID_PATH,
494  * if set; otherwise, it is $GEOGRAPHICLIB_DATA/geoids if the environment
495  * variable GEOGRAPHICLIB_DATA is set; otherwise, it is a compile-time
496  * default (/usr/local/share/GeographicLib/geoids on non-Windows systems
497  * and C:/ProgramData/GeographicLib/geoids on Windows systems).
498  **********************************************************************/
499  static std::string DefaultGeoidPath();
500 
501  /**
502  * @return the default name for the geoid.
503  *
504  * This is the value of the environment variable GEOGRAPHICLIB_GEOID_NAME,
505  * if set; otherwise, it is "egm96-5". The Geoid class does not use this
506  * function; it is just provided as a convenience for a calling program
507  * when constructing a Geoid object.
508  **********************************************************************/
509  static std::string DefaultGeoidName();
510 
511  };
512 
513 } // namespace GeographicLib
514 
515 #if defined(_MSC_VER)
516 # pragma warning (pop)
517 #endif
518 
519 #endif // GEOGRAPHICLIB_GEOID_HPP
const std::string & GeoidDirectory() const
Definition: Geoid.hpp:377
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:90
Math::real MaxError() const
Definition: Geoid.hpp:392
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
Math::real CacheWest() const
Definition: Geoid.hpp:432
Math::real CacheNorth() const
Definition: Geoid.hpp:451
const std::string & DateTime() const
Definition: Geoid.hpp:361
Math::real MajorRadius() const
Definition: Geoid.hpp:469
const std::string & GeoidFile() const
Definition: Geoid.hpp:366
bool Cache() const
Definition: Geoid.hpp:427
Math::real ConvertHeight(real lat, real lon, real h, convertflag d) const
Definition: Geoid.hpp:341
const std::string Interpolation() const
Definition: Geoid.hpp:382
Math::real Offset() const
Definition: Geoid.hpp:409
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
Math::real RMSError() const
Definition: Geoid.hpp:401
bool ThreadSafe() const
Definition: Geoid.hpp:422
Math::real Flattening() const
Definition: Geoid.hpp:478
Exception handling for GeographicLib.
Definition: Constants.hpp:382
Header for GeographicLib::Constants class.
Math::real CacheSouth() const
Definition: Geoid.hpp:459
Math::real CacheEast() const
Definition: Geoid.hpp:441
const std::string & Description() const
Definition: Geoid.hpp:356
const std::string & GeoidName() const
Definition: Geoid.hpp:372
Math::real Scale() const
Definition: Geoid.hpp:417
void CacheAll() const
Definition: Geoid.hpp:271
Math::real operator()(real lat, real lon) const
Definition: Geoid.hpp:297
Math::real operator()(real lat, real lon, real &gradn, real &grade) const
Definition: Geoid.hpp:322
Looking up the height of the geoid.
Definition: Geoid.hpp:91