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