GeographicLib
1.21
|
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