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