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