Geoid.hpp

Go to the documentation of this file.

#if !defined(GEOGRAPHICLIB_GEOID_HPP)
#define GEOGRAPHICLIB_GEOID_HPP 1

#include <vector>
#include <fstream>
#include <GeographicLib/Constants.hpp>

#if defined(_MSC_VER)
// Squelch warnings about dll vs vector and constant conditional expressions
#  pragma warning (push)
#  pragma warning (disable: 4251 4127)
#endif

#if !defined(GEOGRAPHICLIB_GEOID_PGM_PIXEL_WIDTH)

#  define GEOGRAPHICLIB_GEOID_PGM_PIXEL_WIDTH 2
#endif

namespace GeographicLib {

  class GEOGRAPHICLIB_EXPORT Geoid {
  private:
    typedef Math::real real;
#if GEOGRAPHICLIB_GEOID_PGM_PIXEL_WIDTH != 4
    typedef unsigned short pixel_t;
    static const unsigned pixel_size_ = 2;
    static const unsigned pixel_max_ = 0xffffu;
#else
    typedef unsigned pixel_t;
    static const unsigned pixel_size_ = 4;
    static const unsigned pixel_max_ = 0xffffffffu;
#endif
    static const unsigned stencilsize_ = 12;
    static const unsigned nterms_ = ((3 + 1) * (3 + 2))/2; // for a cubic fit
    static const int c0_;
    static const int c0n_;
    static const int c0s_;
    static const int c3_[stencilsize_ * nterms_];
    static const int c3n_[stencilsize_ * nterms_];
    static const int c3s_[stencilsize_ * nterms_];

    std::string _name, _dir, _filename;
    const bool _cubic;
    const real _a, _e2, _degree, _eps;
    mutable std::ifstream _file;
    real _rlonres, _rlatres;
    std::string _description, _datetime;
    real _offset, _scale, _maxerror, _rmserror;
    int _width, _height;
    unsigned long long _datastart, _swidth;
    bool _threadsafe;
    // Area cache
    mutable std::vector< std::vector<pixel_t> > _data;
    mutable bool _cache;
    // NE corner and extent of cache
    mutable int _xoffset, _yoffset, _xsize, _ysize;
    // Cell cache
    mutable int _ix, _iy;
    mutable real _v00, _v01, _v10, _v11;
    mutable real _t[nterms_];
    void filepos(int ix, int iy) const {
      _file.seekg(
#if !(defined(__GNUC__) && __GNUC__ < 4)
                  // g++ 3.x doesn't know about the cast to streamoff.
                  std::ios::streamoff
#endif
                  (_datastart +
                   pixel_size_ * (unsigned(iy)*_swidth + unsigned(ix))));
    }
    real rawval(int ix, int iy) const {
      if (ix < 0)
        ix += _width;
      else if (ix >= _width)
        ix -= _width;
      if (_cache && iy >= _yoffset && iy < _yoffset + _ysize &&
          ((ix >= _xoffset && ix < _xoffset + _xsize) ||
           (ix + _width >= _xoffset && ix + _width < _xoffset + _xsize))) {
        return real(_data[iy - _yoffset]
                    [ix >= _xoffset ? ix - _xoffset : ix + _width - _xoffset]);
      } else {
        if (iy < 0 || iy >= _height) {
          iy = iy < 0 ? -iy : 2 * (_height - 1) - iy;
          ix += (ix < _width/2 ? 1 : -1) * _width/2;
        }
        try {
          filepos(ix, iy);
          // initial values to suppress warnings in case get fails
          char a = 0, b = 0;
          _file.get(a);
          _file.get(b);
          unsigned r = ((unsigned char)(a) << 8) | (unsigned char)(b);
          if (pixel_size_ == 4) {
            _file.get(a);
            _file.get(b);
            r = (r << 16) | ((unsigned char)(a) << 8) | (unsigned char)(b);
          }
          return real(r);
        }
        catch (const std::exception& e) {
          // throw GeographicErr("Error reading " + _filename + ": "
          //                      + e.what());
          // triggers complaints about the "binary '+'" under Visual Studio.
          // So use '+=' instead.
          std::string err("Error reading ");
          err += _filename;
          err += ": ";
          err += e.what();
          throw GeographicErr(err);
        }
      }
    }
    real height(real lat, real lon) const;
    Geoid(const Geoid&);            // copy constructor not allowed
    Geoid& operator=(const Geoid&); // copy assignment not allowed
  public:

    enum convertflag {
      ELLIPSOIDTOGEOID = -1,
      NONE = 0,
      GEOIDTOELLIPSOID = 1,
    };


    explicit Geoid(const std::string& name, const std::string& path = "",
                   bool cubic = true, bool threadsafe = false);

    void CacheArea(real south, real west, real north, real east) const;

    void CacheAll() const { CacheArea(real(-90), real(0),
                                      real(90), real(360)); }

    void CacheClear() const;



    Math::real operator()(real lat, real lon) const {
      return height(lat, lon);
    }

    Math::real ConvertHeight(real lat, real lon, real h,
                             convertflag d) const {
      return h + real(d) * height(lat, lon);
    }



    const std::string& Description() const { return _description; }

    const std::string& DateTime() const { return _datetime; }

    const std::string& GeoidFile() const { return _filename; }

    const std::string& GeoidName() const { return _name; }

    const std::string& GeoidDirectory() const { return _dir; }

    const std::string Interpolation() const
    { return std::string(_cubic ? "cubic" : "bilinear"); }

    Math::real MaxError() const { return _maxerror; }

    Math::real RMSError() const { return _rmserror; }

    Math::real Offset() const { return _offset; }

    Math::real Scale() const { return _scale; }

    bool ThreadSafe() const { return _threadsafe; }

    bool Cache() const { return _cache; }

    Math::real CacheWest() const {
      return _cache ? ((_xoffset + (_xsize == _width ? 0 : _cubic)
                        + _width/2) % _width - _width/2) / _rlonres :
        0;
    }

    Math::real CacheEast() const {
      return  _cache ?
        CacheWest() +
        (_xsize - (_xsize == _width ? 0 : 1 + 2 * _cubic)) / _rlonres :
        0;
    }

    Math::real CacheNorth() const {
      return _cache ? 90 - (_yoffset + _cubic) / _rlatres : 0;
    }

    Math::real CacheSouth() const {
      return _cache ? 90 - ( _yoffset + _ysize - 1 - _cubic) / _rlatres : 0;
    }

    Math::real MajorRadius() const
    { return Constants::WGS84_a(); }

    Math::real Flattening() const { return Constants::WGS84_f(); }

    static std::string DefaultGeoidPath();

    static std::string DefaultGeoidName();

  };

} // namespace GeographicLib

#if defined(_MSC_VER)
#  pragma warning (pop)
#endif

#endif  // GEOGRAPHICLIB_GEOID_HPP