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