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) <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 6785 2010-01-05 22:15:42Z karney $"
00012 
00013 #include "GeographicLib/Constants.hpp"
00014 #include <vector>
00015 #include <fstream>
00016 
00017 namespace GeographicLib {
00018 
00019   /**
00020    * \brief Computing the height of the geoid
00021    *
00022    * This class evaluated the height of one of the standard geoids, EGM84,
00023    * EGM96, or EGM2008 by bilinear or cubic interpolation into a rectangular
00024    * grid of data.  These geoid models are documented in
00025    * - EGM84:
00026    *   http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html
00027    * - EGM96:
00028    *   http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html
00029    * - EGM2008:
00030    *   http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008
00031    *
00032    * The geoids are defined in terms of spherical harmonics.  However in order
00033    * to provide a quick and flexible method of evaluating the geoid heights,
00034    * this class evaluates the height by interpolation inot a grid of
00035    * precomputed values.
00036    *
00037    * See \ref geoid for details of how to install the data sets, the data
00038    * format, estimates of the interpolation errors, and how to use caching.
00039    *
00040    * In addition to returning the geoid height, the gradient of the geoid can
00041    * be calculated.  The gradient is defined as the rate of change of the geoid
00042    * as a function of position on the ellipsoid.  This uses the parameters for
00043    * the WGS84 ellipsoid.  The gradient defined in terms of the interpolated
00044    * heights.
00045    *
00046    * This class is \e not thread safe in that a single instantiation cannot be
00047    * safely used by multiple threads.  If multiple threads need to calculate
00048    * geoid heights they should all construct thread-local instantiations.
00049    **********************************************************************/
00050 
00051   class Geoid {
00052   private:
00053     typedef Math::real real;
00054     static const unsigned stencilsize = 12;
00055     static const unsigned nterms = ((3 + 1) * (3 + 2))/2; // for a cubic fit
00056     static const real c0, c0n, c0s;
00057     static const real c3[stencilsize * nterms];
00058     static const real c3n[stencilsize * nterms];
00059     static const real c3s[stencilsize * nterms];
00060 
00061     std::string _name, _dir, _filename;
00062     const bool _cubic;
00063     const real _a, _e2, _degree, _eps;
00064     mutable std::ifstream _file;
00065     real _rlonres, _rlatres;
00066     std::string _description, _datetime;
00067     real _offset, _scale, _maxerror, _rmserror;
00068     int _width, _height;
00069     unsigned long long _datastart, _swidth;
00070     // Area cache
00071     mutable std::vector< std::vector<unsigned short> > _data;
00072     mutable bool _cache;
00073     // NE corner and extent of cache
00074     mutable int _xoffset, _yoffset, _xsize, _ysize;
00075     // Cell cache
00076     mutable int _ix, _iy;
00077     mutable real _v00, _v01, _v10, _v11;
00078     mutable real _t[nterms];
00079     void filepos(int ix, int iy) const {
00080       _file.seekg(std::ios::streamoff(_datastart +
00081                                       2ULL * (unsigned(iy) * _swidth +
00082                                               unsigned(ix))));
00083     }
00084     real rawval(int ix, int iy) const {
00085       if (ix < 0)
00086         ix += _width;
00087       else if (ix >= _width)
00088         ix -= _width;
00089       if (_cache && iy >= _yoffset && iy < _yoffset + _ysize &&
00090           ((ix >= _xoffset && ix < _xoffset + _xsize) ||
00091            (ix + _width >= _xoffset && ix + _width < _xoffset + _xsize))) {
00092         return real(_data[iy - _yoffset]
00093                     [ix >= _xoffset ? ix - _xoffset : ix + _width - _xoffset]);
00094       } else {
00095         if (iy < 0 || iy >= _height) {
00096           iy = iy < 0 ? -iy : 2 * (_height - 1) - iy;
00097           ix += (ix < _width/2 ? 1 : -1) * _width/2;
00098         }
00099         try {
00100           filepos(ix, iy);
00101           char a, b;
00102           _file.get(a);
00103           _file.get(b);
00104           return real((unsigned char)(a) * 256u + (unsigned char)(b));
00105         }
00106         catch (const std::exception& e) {
00107           // throw GeographicErr("Error reading " + _filename + ": "
00108           //                      + e.what());
00109           // triggers complaints about the "binary '+'" under Visual Studio.
00110           // So use '+=' instead.
00111           std::string err("Error reading ");
00112           err += _filename;
00113           err += ": ";
00114           err += e.what();
00115           throw GeographicErr(err);
00116 
00117         }
00118       }
00119     }
00120     real height(real lat, real lon, bool gradp,
00121                 real& grade, real& gradn) const;
00122     Geoid(const Geoid&);            // copy constructor not allowed
00123     Geoid& operator=(const Geoid&); // copy assignment not allowed
00124   public:
00125 
00126     /**
00127      * Create a Geoid loading the data for geoid \e name.  The data file is
00128      * formed by appending ".pgm" to the name.  If \e path is specified (and is
00129      * non-empty), then the file is loaded from directory, \e path.  Otherwise
00130      * the path is given by the GEOID_PATH environment variable.  If that is
00131      * undefined, a compile-time default path is used
00132      * (/usr/local/share/GeographicLib/geoids on non-Windows systems and
00133      * C:/cygwin/usr/local/share/GeographicLib/geoids on Windows systems).  The
00134      * final \e cubic argument specifies whether to use bilinear (\e cubic =
00135      * false) or cubic (\e cubic = true, the default) interpolation.  This may
00136      * throw an error because the file does not exist, is unreadable, or is
00137      * corrupt.
00138      **********************************************************************/
00139     explicit Geoid(const std::string& name, const std::string& path = "",
00140                    bool cubic = true);
00141 
00142     /**
00143      * Cache the data for the rectangular area defined by the four arguments \e
00144      * south, \e west, \e north, \e east (all in degrees).  \e east is always
00145      * interpreted as being east of \e west, if necessary by adding
00146      * 360<sup>o</sup> to its value.  This may throw an error because of
00147      * insufficent memory or because of an error reading the data from the
00148      * file.  In this case, you can catch the error and either do nothing (you
00149      * will have no cache in this case) or try again with a smaller area.  \e
00150      * south and \e north should be in the range [-90, 90]; \e west and \e east
00151      * should be in the range [-180, 360].
00152      **********************************************************************/
00153     void CacheArea(real south, real west, real north, real east) const;
00154 
00155     /**
00156      * Cache all the data.  On most computers, this is fast for data sets with
00157      * grid resolution of 5' or coarser.  For a 1' grid, the required RAM is
00158      * 450MB; a 2.5' grid needs 72MB; and a 5' grid needs 18MB.  This may throw
00159      * an error because of insufficent memory or because of an error reading
00160      * the data from the file.  In this case, you can catch the error and
00161      * either do nothing (you will have no cache in this case) or try using
00162      * Geoid::CacheArea on a specific area.
00163      **********************************************************************/
00164     void CacheAll() const { CacheArea(real(-90), real(0),
00165                                       real(90), real(360)); }
00166 
00167     /**
00168      * Clear the cache.  This never throws an error.
00169      **********************************************************************/
00170     void CacheClear() const throw();
00171 
00172     /**
00173      * Return the geoid height in meters for latitude \e lat (in [-90, 90]) and
00174      * longitude \e lon (in [-180,360]), both in degrees.  This may throw an
00175      * error because of an error reading data from disk.  However, it will not
00176      * throw if (\e lat, \e lon) is within a successfully cached area.
00177      **********************************************************************/
00178     Math::real operator()(real lat, real lon) const {
00179       real gradn, grade;
00180       return height(lat, lon, false, gradn, grade);
00181     }
00182     /**
00183      * Return the geoid height in meters for latitude \e lat (in [-90, 90]) and
00184      * longitude \e lon (in [-180,360]), both in degrees.  In addition compute
00185      * the gradient of the geoid height in the northerly \e gradn and easterly
00186      * \e grade directions.  This may throw an error because of an error
00187      * reading data from disk.  However, it will not throw if (\e lat, \e lon)
00188      * is within a successfully cached area.
00189      **********************************************************************/
00190     Math::real operator()(real lat, real lon, real& gradn, real& grade) const {
00191       return height(lat, lon, true, gradn, grade);
00192     }
00193 
00194     /**
00195      * Return the geoid description if available in the data file.  If absent,
00196      * return "NONE".
00197      **********************************************************************/
00198     const std::string& Description() const throw() { return _description; }
00199 
00200     /**
00201      * Return the date of the data file.  If absent, return "UNKNOWN".
00202      **********************************************************************/
00203     const std::string& DateTime() const throw() { return _datetime; }
00204 
00205     /**
00206      * Return the full file name used to load the geoid data.
00207      **********************************************************************/
00208     const std::string& GeoidFile() const throw() { return _filename; }
00209 
00210     /**
00211      * Return the "name" used to load the geoid data (from the first argument
00212      * of the constructor).
00213      **********************************************************************/
00214     const std::string& GeoidName() const throw() { return _name; }
00215 
00216     /**
00217      * Return the directory used to load the geoid data.
00218      **********************************************************************/
00219     const std::string& GeoidDirectory() const throw() { return _dir; }
00220 
00221     /**
00222      * Return the interpolation method (cubic or bilinear).
00223      **********************************************************************/
00224     const std::string Interpolation() const
00225     { return std::string(_cubic ? "cubic" : "bilinear"); }
00226 
00227     /**
00228      * Return a estimate of the maximum interpolation and quantization error
00229      * (meters).  This relies on the value being stored in the data file.  If
00230      * the value is absent, return -1.
00231      **********************************************************************/
00232     Math::real MaxError() const throw() { return _maxerror; }
00233 
00234     /**
00235      * Return a estimate of the RMS interpolation and quantization error
00236      * (meters).  This relies on the value being stored in the data file.  If
00237      * the value is absent, return -1.
00238      **********************************************************************/
00239     Math::real RMSError() const throw() { return _rmserror; }
00240 
00241     /**
00242      * Return offset (meters) for converting pixel values to geoid heights.
00243      **********************************************************************/
00244     Math::real Offset() const throw() { return _offset; }
00245 
00246     /**
00247      * Return scale (meters) for converting pixel values to geoid
00248      * heights.
00249      **********************************************************************/
00250     Math::real Scale() const throw() { return _scale; }
00251 
00252     /**
00253      * Is a data cache active?
00254      **********************************************************************/
00255     bool Cache() const throw() { return _cache; }
00256 
00257     /**
00258      * Return the west edge of the cached area.  The cache includes this edge.
00259      **********************************************************************/
00260     Math::real CacheWest() const throw() {
00261       return _cache ? ((_xoffset + (_xsize == _width ? 0 : _cubic)
00262                         + _width/2) % _width - _width/2) / _rlonres :
00263         0;
00264     }
00265 
00266     /**
00267      * Return the east edge of the cached area.  The cache excludes this edge.
00268      **********************************************************************/
00269     Math::real CacheEast() const throw() {
00270       return  _cache ?
00271         CacheWest() +
00272         (_xsize - (_xsize == _width ? 0 : 1 + 2 * _cubic)) / _rlonres :
00273         0;
00274     }
00275 
00276     /**
00277      * Return the north edge of the cached area.  The cache includes this edge.
00278      **********************************************************************/
00279     Math::real CacheNorth() const throw() {
00280       return _cache ? 90 - (_yoffset + _cubic) / _rlatres : 0;
00281     }
00282 
00283     /**
00284      * Return the south edge of the cached area.  The cache excludes this edge
00285      * unless it's the south pole.
00286      **********************************************************************/
00287     Math::real CacheSouth() const throw() {
00288       return _cache ? 90 - ( _yoffset + _ysize - 1 - _cubic) / _rlatres : 0;
00289     }
00290 
00291     /**
00292      * The major radius of the ellipsoid (meters).  This is the value for the
00293      * WGS84 ellipsoid because the supported geoid models are all based on this
00294      * ellipsoid.
00295      **********************************************************************/
00296     Math::real MajorRadius() const throw() { return Constants::WGS84_a(); }
00297 
00298     /**
00299      * The inverse flattening of the ellipsoid.  This is the value for the
00300      * WGS84 ellipsoid because the supported geoid models are all based on this
00301      * ellipsoid.
00302      **********************************************************************/
00303     Math::real InverseFlattening() const throw()
00304     { return Constants::WGS84_r(); }
00305 
00306     /**
00307      * Return the compile-time default path for the geoid data files.
00308      **********************************************************************/
00309     static std::string DefaultPath();
00310 
00311     /**
00312      * Return the value of the environment variable GEOID_PATH.
00313      **********************************************************************/
00314     static std::string GeoidPath();
00315   };
00316 
00317 } // namespace GeographicLib
00318 #endif

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1