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