GeographicLib  1.21
Geoid.hpp
Go to the documentation of this file.
00001 /**
00002  * \file Geoid.hpp
00003  * \brief Header for GeographicLib::Geoid class
00004  *
00005  * Copyright (c) Charles Karney (2009-2012) <charles@karney.com> and licensed
00006  * under the MIT/X11 License.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #if !defined(GEOGRAPHICLIB_GEOID_HPP)
00011 #define GEOGRAPHICLIB_GEOID_HPP \
00012   "$Id: 4e4eb5941d16ad00416798703d246a6f7ef5fe46 $"
00013 
00014 #include <string>
00015 #include <vector>
00016 #include <fstream>
00017 #include <GeographicLib/Constants.hpp>
00018 
00019 #if defined(_MSC_VER)
00020 // Squelch warnings about dll vs vector
00021 #pragma warning (push)
00022 #pragma warning (disable: 4251)
00023 #endif
00024 
00025 #if !defined(PGM_PIXEL_WIDTH)
00026 /**
00027  * The size of the pixel data in the pgm data files for the geoids.  2
00028  * is the standard size corresponding to a maxval 2^16-1.  Setting it
00029  * to 4 uses a maxval of 2^32-1 and changes the extension for the data
00030  * files from .pgm to .pgm4.  Note that the format of these pgm4 files
00031  * is a non-standard extension of the pgm format.
00032  **********************************************************************/
00033 #define PGM_PIXEL_WIDTH 2
00034 #endif
00035 
00036 namespace GeographicLib {
00037 
00038   /**
00039    * \brief Looking up the height of the geoid
00040    *
00041    * This class evaluated the height of one of the standard geoids, EGM84,
00042    * EGM96, or EGM2008 by bilinear or cubic interpolation into a rectangular
00043    * grid of data.  These geoid models are documented in
00044    * - EGM84:
00045    *   http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html
00046    * - EGM96:
00047    *   http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html
00048    * - EGM2008:
00049    *   http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008
00050    *
00051    * The geoids are defined in terms of spherical harmonics.  However in order
00052    * to provide a quick and flexible method of evaluating the geoid heights,
00053    * this class evaluates the height by interpolation into a grid of
00054    * precomputed values.
00055    *
00056    * See \ref geoid for details of how to install the data sets, the data
00057    * format, estimates of the interpolation errors, and how to use caching.
00058    *
00059    * In addition to returning the geoid height, the gradient of the geoid can
00060    * be calculated.  The gradient is defined as the rate of change of the geoid
00061    * as a function of position on the ellipsoid.  This uses the parameters for
00062    * the WGS84 ellipsoid.  The gradient defined in terms of the interpolated
00063    * heights.  As a result of the way that the geoid data is stored, the
00064    * calculation of gradients can result in large quantization errors.  This is
00065    * particularly acute for fine grids, at high latitudes, and for the easterly
00066    * gradient.
00067    *
00068    * This class is typically \e not thread safe in that a single instantiation
00069    * cannot be safely used by multiple threads because of the way the object
00070    * reads the data set and because it maintains a single-cell cache.  If
00071    * multiple threads need to calculate geoid heights they should all construct
00072    * thread-local instantiations.  Alternatively, set the optional \e
00073    * threadsafe parameter to true in the constructor.  This causes the
00074    * constructor to read all the data into memory and to turn off the
00075    * single-cell caching which results in a Geoid object which \e is thread
00076    * safe.
00077    *
00078    * Example of use:
00079    * \include example-Geoid.cpp
00080    *
00081    * <a href="GeoidEval.1.html">GeoidEval</a> is a command-line utility
00082    * providing access to the functionality of Geoid.
00083    **********************************************************************/
00084 
00085   class GEOGRAPHIC_EXPORT Geoid {
00086   private:
00087     typedef Math::real real;
00088 #if PGM_PIXEL_WIDTH != 4
00089     typedef unsigned short pixel_t;
00090     static const unsigned pixel_size_ = 2;
00091     static const unsigned pixel_max_ = 0xffffu;
00092 #else
00093     typedef unsigned pixel_t;
00094     static const unsigned pixel_size_ = 4;
00095     static const unsigned pixel_max_ = 0xffffffffu;
00096 #endif
00097     static const unsigned stencilsize_ = 12;
00098     static const unsigned nterms_ = ((3 + 1) * (3 + 2))/2; // for a cubic fit
00099     static const real c0_;
00100     static const real c0n_;
00101     static const real c0s_;
00102     static const real c3_[stencilsize_ * nterms_];
00103     static const real c3n_[stencilsize_ * nterms_];
00104     static const real c3s_[stencilsize_ * nterms_];
00105 
00106     std::string _name, _dir, _filename;
00107     const bool _cubic;
00108     const real _a, _e2, _degree, _eps;
00109     mutable std::ifstream _file;
00110     real _rlonres, _rlatres;
00111     std::string _description, _datetime;
00112     real _offset, _scale, _maxerror, _rmserror;
00113     int _width, _height;
00114     unsigned long long _datastart, _swidth;
00115     bool _threadsafe;
00116     // Area cache
00117     mutable std::vector< std::vector<pixel_t> > _data;
00118     mutable bool _cache;
00119     // NE corner and extent of cache
00120     mutable int _xoffset, _yoffset, _xsize, _ysize;
00121     // Cell cache
00122     mutable int _ix, _iy;
00123     mutable real _v00, _v01, _v10, _v11;
00124     mutable real _t[nterms_];
00125     void filepos(int ix, int iy) const {
00126       _file.seekg(
00127 #if !(defined(__GNUC__) && __GNUC__ < 4)
00128                   // g++ 3.x doesn't know about the cast to streamoff.
00129                   std::ios::streamoff
00130 #endif
00131                   (_datastart +
00132                    pixel_size_ * (unsigned(iy)*_swidth + unsigned(ix))));
00133     }
00134     real rawval(int ix, int iy) const {
00135       if (ix < 0)
00136         ix += _width;
00137       else if (ix >= _width)
00138         ix -= _width;
00139       if (_cache && iy >= _yoffset && iy < _yoffset + _ysize &&
00140           ((ix >= _xoffset && ix < _xoffset + _xsize) ||
00141            (ix + _width >= _xoffset && ix + _width < _xoffset + _xsize))) {
00142         return real(_data[iy - _yoffset]
00143                     [ix >= _xoffset ? ix - _xoffset : ix + _width - _xoffset]);
00144       } else {
00145         if (iy < 0 || iy >= _height) {
00146           iy = iy < 0 ? -iy : 2 * (_height - 1) - iy;
00147           ix += (ix < _width/2 ? 1 : -1) * _width/2;
00148         }
00149         try {
00150           filepos(ix, iy);
00151           char a, b;
00152           _file.get(a);
00153           _file.get(b);
00154           unsigned r = ((unsigned char)(a) << 8) | (unsigned char)(b);
00155           if (pixel_size_ == 4) {
00156             _file.get(a);
00157             _file.get(b);
00158             r = (r << 16) | ((unsigned char)(a) << 8) | (unsigned char)(b);
00159           }
00160           return real(r);
00161         }
00162         catch (const std::exception& e) {
00163           // throw GeographicErr("Error reading " + _filename + ": "
00164           //                      + e.what());
00165           // triggers complaints about the "binary '+'" under Visual Studio.
00166           // So use '+=' instead.
00167           std::string err("Error reading ");
00168           err += _filename;
00169           err += ": ";
00170           err += e.what();
00171           throw GeographicErr(err);
00172         }
00173       }
00174     }
00175     real height(real lat, real lon, bool gradp,
00176                 real& grade, real& gradn) const;
00177     Geoid(const Geoid&);            // copy constructor not allowed
00178     Geoid& operator=(const Geoid&); // copy assignment not allowed
00179   public:
00180 
00181     /**
00182      * Flags indicating conversions between heights above the geoid and heights
00183      * above the ellipsoid.
00184      **********************************************************************/
00185     enum convertflag {
00186       /**
00187        * The multiplier for converting from heights above the geoid to heights
00188        * above the ellipsoid.
00189        **********************************************************************/
00190       ELLIPSOIDTOGEOID = -1,
00191       /**
00192        * No conversion.
00193        **********************************************************************/
00194       NONE = 0,
00195       /**
00196        * The multiplier for converting from heights above the ellipsoid to
00197        * heights above the geoid.
00198        **********************************************************************/
00199       GEOIDTOELLIPSOID = 1,
00200     };
00201 
00202     /** \name Setting up the geoid
00203      **********************************************************************/
00204     ///@{
00205     /**
00206      * Construct a geoid.
00207      *
00208      * @param[in] name the name of the geoid.
00209      * @param[in] path (optional) directory for data file.
00210      * @param[in] cubic (optional) interpolation method; false means bilinear,
00211      *   true (the default) means cubic.
00212      * @param[in] threadsafe (optional), if true, construct a thread safe
00213      *   object.  The default is false
00214      *
00215      * The data file is formed by appending ".pgm" to the name.  If \e path is
00216      * specified (and is non-empty), then the file is loaded from directory, \e
00217      * path.  Otherwise the path is given by DefaultGeoidPath().  This may
00218      * throw an exception because the file does not exist, is unreadable, or is
00219      * corrupt.  If the \e threadsafe parameter is true, the data set is read
00220      * into memory (which this may also cause an exception to be thrown), the
00221      * data file is closed, and single-cell caching is turned off; this results
00222      * in a Geoid object which \e is thread safe.
00223      **********************************************************************/
00224     explicit Geoid(const std::string& name, const std::string& path = "",
00225                    bool cubic = true, bool threadsafe = false);
00226 
00227     /**
00228      * Set up a cache.
00229      *
00230      * @param[in] south latitude (degrees) of the south edge of the cached area.
00231      * @param[in] west longitude (degrees) of the west edge of the cached area.
00232      * @param[in] north latitude (degrees) of the north edge of the cached area.
00233      * @param[in] east longitude (degrees) of the east edge of the cached area.
00234      *
00235      * Cache the data for the specified "rectangular" area bounded by the
00236      * parallels \e south and \e north and the meridians \e west and \e east.
00237      * \e east is always interpreted as being east of \e west, if necessary by
00238      * adding 360<sup>o</sup> to its value.  This may throw an error because of
00239      * insufficient memory or because of an error reading the data from the
00240      * file.  In this case, you can catch the error and either do nothing (you
00241      * will have no cache in this case) or try again with a smaller area.  \e
00242      * south and \e north should be in the range [-90, 90]; \e west and \e east
00243      * should be in the range [-180, 360].  An exception is thrown if this
00244      * routine is called on a thread safe Geoid.
00245      **********************************************************************/
00246     void CacheArea(real south, real west, real north, real east) const;
00247 
00248     /**
00249      * Cache all the data.  On most computers, this is fast for data sets with
00250      * grid resolution of 5' or coarser.  For a 1' grid, the required RAM is
00251      * 450MB; a 2.5' grid needs 72MB; and a 5' grid needs 18MB.  This may throw
00252      * an error because of insufficient memory or because of an error reading
00253      * the data from the file.  In this case, you can catch the error and
00254      * either do nothing (you will have no cache in this case) or try using
00255      * Geoid::CacheArea on a specific area.  An exception is thrown if this
00256      * routine is called on a thread safe Geoid.
00257      **********************************************************************/
00258     void CacheAll() const { CacheArea(real(-90), real(0),
00259                                       real(90), real(360)); }
00260 
00261     /**
00262      * Clear the cache.  This never throws an error.  (This does nothing with a
00263      * thread safe Geoid.)
00264      **********************************************************************/
00265     void CacheClear() const throw();
00266 
00267     ///@}
00268 
00269     /** \name Compute geoid heights
00270      **********************************************************************/
00271     ///@{
00272     /**
00273      * Compute the geoid height at a point
00274      *
00275      * @param[in] lat latitude of the point (degrees).
00276      * @param[in] lon longitude of the point (degrees).
00277      * @return geoid height (meters).
00278      *
00279      * The latitude should be in [-90, 90] and longitude should be in
00280      * [-180,360].  This may throw an error because of an error reading data
00281      * from disk.  However, it will not throw if (\e lat, \e lon) is within a
00282      * successfully cached area.
00283      **********************************************************************/
00284     Math::real operator()(real lat, real lon) const {
00285       real gradn, grade;
00286       return height(lat, lon, false, gradn, grade);
00287     }
00288 
00289     /**
00290      * Compute the geoid height and gradient at a point
00291      *
00292      * @param[in] lat latitude of the point (degrees).
00293      * @param[in] lon longitude of the point (degrees).
00294      * @param[out] gradn northerly gradient (dimensionless).
00295      * @param[out] grade easterly gradient (dimensionless).
00296      * @return geoid height (meters).
00297      *
00298      * The latitude should be in [-90, 90] and longitude should be in [-180,
00299      * 360].  This may throw an error because of an error reading data from
00300      * disk.  However, it will not throw if (\e lat, \e lon) is within a
00301      * successfully cached area.  As a result of the way that the geoid data is
00302      * stored, the calculation of gradients can result in large quantization
00303      * errors.  This is particularly acute for fine grids, at high latitudes,
00304      * and for the easterly gradient.  If you need to compute the direction of
00305      * the acceleration due to gravity accurately, you should use
00306      * GravityModel::Gravity.
00307      **********************************************************************/
00308     Math::real operator()(real lat, real lon, real& gradn, real& grade) const {
00309       return height(lat, lon, true, gradn, grade);
00310     }
00311 
00312     /**
00313      * Convert a height above the geoid to a height above the ellipsoid and
00314      * vice versa.
00315      *
00316      * @param[in] lat latitude of the point (degrees).
00317      * @param[in] lon longitude of the point (degrees).
00318      * @param[in] h height of the point (degrees).
00319      * @param[in] d a Geoid::convertflag specifying the direction of the
00320      *   conversion; Geoid::GEOIDTOELLIPSOID means convert a height above the
00321      *   geoid to a height above the ellipsoid; Geoid::ELLIPSOIDTOGEOID means
00322      *   convert a height above the ellipsoid to a height above the geoid.
00323      * @return converted height (meters).
00324      **********************************************************************/
00325     Math::real ConvertHeight(real lat, real lon, real h,
00326                              convertflag d) const {
00327       real gradn, grade;
00328       return h + real(d) * height(lat, lon, true, gradn, grade);
00329     }
00330 
00331     ///@}
00332 
00333     /** \name Inspector functions
00334      **********************************************************************/
00335     ///@{
00336     /**
00337      * @return geoid description, if available, in the data file; if
00338      *   absent, return "NONE".
00339      **********************************************************************/
00340     const std::string& Description() const throw() { return _description; }
00341 
00342     /**
00343      * @return date of the data file; if absent, return "UNKNOWN".
00344      **********************************************************************/
00345     const std::string& DateTime() const throw() { return _datetime; }
00346 
00347     /**
00348      * @return full file name used to load the geoid data.
00349      **********************************************************************/
00350     const std::string& GeoidFile() const throw() { return _filename; }
00351 
00352     /**
00353      * @return "name" used to load the geoid data (from the first argument of
00354      *   the constructor).
00355      **********************************************************************/
00356     const std::string& GeoidName() const throw() { return _name; }
00357 
00358     /**
00359      * @return directory used to load the geoid data.
00360      **********************************************************************/
00361     const std::string& GeoidDirectory() const throw() { return _dir; }
00362 
00363     /**
00364      * @return interpolation method ("cubic" or "bilinear").
00365      **********************************************************************/
00366     const std::string Interpolation() const
00367     { return std::string(_cubic ? "cubic" : "bilinear"); }
00368 
00369     /**
00370      * @return estimate of the maximum interpolation and quantization error
00371      *   (meters).
00372      *
00373      * This relies on the value being stored in the data file.  If the value is
00374      * absent, return -1.
00375      **********************************************************************/
00376     Math::real MaxError() const throw() { return _maxerror; }
00377 
00378     /**
00379      * @return estimate of the RMS interpolation and quantization error
00380      *   (meters).
00381      *
00382      * This relies on the value being stored in the data file.  If the value is
00383      * absent, return -1.
00384      **********************************************************************/
00385     Math::real RMSError() const throw() { return _rmserror; }
00386 
00387     /**
00388      * @return offset (meters).
00389      *
00390      * This in used in converting from the pixel values in the data file to
00391      * geoid heights.
00392      **********************************************************************/
00393     Math::real Offset() const throw() { return _offset; }
00394 
00395     /**
00396      * @return scale (meters).
00397      *
00398      * This in used in converting from the pixel values in the data file to
00399      * geoid heights.
00400      **********************************************************************/
00401     Math::real Scale() const throw() { return _scale; }
00402 
00403     /**
00404      * @return true if the object is constructed to be thread safe.
00405      **********************************************************************/
00406     bool ThreadSafe() const throw() { return _threadsafe; }
00407 
00408     /**
00409      * @return true if a data cache is active.
00410      **********************************************************************/
00411     bool Cache() const throw() { return _cache; }
00412 
00413     /**
00414      * @return west edge of the cached area; the cache includes this edge.
00415      **********************************************************************/
00416     Math::real CacheWest() const throw() {
00417       return _cache ? ((_xoffset + (_xsize == _width ? 0 : _cubic)
00418                         + _width/2) % _width - _width/2) / _rlonres :
00419         0;
00420     }
00421 
00422     /**
00423      * @return east edge of the cached area; the cache excludes this edge.
00424      **********************************************************************/
00425     Math::real CacheEast() const throw() {
00426       return  _cache ?
00427         CacheWest() +
00428         (_xsize - (_xsize == _width ? 0 : 1 + 2 * _cubic)) / _rlonres :
00429         0;
00430     }
00431 
00432     /**
00433      * @return north edge of the cached area; the cache includes this edge.
00434      **********************************************************************/
00435     Math::real CacheNorth() const throw() {
00436       return _cache ? 90 - (_yoffset + _cubic) / _rlatres : 0;
00437     }
00438 
00439     /**
00440      * @return south edge of the cached area; the cache excludes this edge
00441      *   unless it's the south pole.
00442      **********************************************************************/
00443     Math::real CacheSouth() const throw() {
00444       return _cache ? 90 - ( _yoffset + _ysize - 1 - _cubic) / _rlatres : 0;
00445     }
00446 
00447     /**
00448      * @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
00449      *
00450      * (The WGS84 value is returned because the supported geoid models are all
00451      * based on this ellipsoid.)
00452      **********************************************************************/
00453     Math::real MajorRadius() const throw()
00454     { return Constants::WGS84_a<real>(); }
00455 
00456     /**
00457      * @return \e f the flattening of the WGS84 ellipsoid.
00458      *
00459      * (The WGS84 value is returned because the supported geoid models are all
00460      * based on this ellipsoid.)
00461      **********************************************************************/
00462     Math::real Flattening() const throw() { return Constants::WGS84_f<real>(); }
00463     ///@}
00464 
00465     /// \cond SKIP
00466     /**
00467      * <b>DEPRECATED</b>
00468      * @return \e r the inverse flattening of the WGS84 ellipsoid.
00469      **********************************************************************/
00470     Math::real InverseFlattening() const throw()
00471     { return 1/Constants::WGS84_f<real>(); }
00472     /// \endcond
00473 
00474     /**
00475      * @return the default path for geoid data files.
00476      *
00477      * This is the value of the environment variable GEOID_PATH, if set;
00478      * otherwise, it is $GEOGRAPHICLIB_DATA/geoids if the environment variable
00479      * GEOGRAPHICLIB_DATA is set; otherwise, it is a compile-time default
00480      * (/usr/local/share/GeographicLib/geoids on non-Windows systems and
00481      * C:/Documents and Settings/All Users/Application
00482      * Data/GeographicLib/geoids on Windows systems).
00483      **********************************************************************/
00484     static std::string DefaultGeoidPath();
00485 
00486     /**
00487      * @return the default name for the geoid.
00488      *
00489      * This is the value of the environment variable GEOID_NAME, if set,
00490      * otherwise, it is "egm96-5".  The Geoid class does not use this function;
00491      * it is just provided as a convenience for a calling program when
00492      * constructing a Geoid object.
00493      **********************************************************************/
00494     static std::string DefaultGeoidName();
00495 
00496   };
00497 
00498 } // namespace GeographicLib
00499 
00500 #if defined(_MSC_VER)
00501 #pragma warning (pop)
00502 #endif
00503 
00504 #endif  // GEOGRAPHICLIB_GEOID_HPP