Geoid.cpp

Go to the documentation of this file.
00001 /**
00002  * \file Geoid.cpp
00003  * \brief Implementation 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 #include "GeographicLib/Geoid.hpp"
00011 #include <sstream>
00012 #include <cstdlib>
00013 
00014 #define GEOGRAPHICLIB_GEOID_CPP "$Id: Geoid.cpp 6785 2010-01-05 22:15:42Z karney $"
00015 
00016 RCSID_DECL(GEOGRAPHICLIB_GEOID_CPP)
00017 RCSID_DECL(GEOGRAPHICLIB_GEOID_HPP)
00018 
00019 #if !defined(GEOID_DEFAULT_PATH)
00020 #if defined(_MSC_VER)
00021 #define GEOID_DEFAULT_PATH "C:/cygwin/usr/local/share/GeographicLib/geoids"
00022 #else
00023 #define GEOID_DEFAULT_PATH "/usr/local/share/GeographicLib/geoids"
00024 #endif
00025 #endif
00026 
00027 #if defined(_MSC_VER)
00028 // Squelch warnings about unsafe use of getenv
00029 #pragma warning (disable: 4996)
00030 #endif
00031 
00032 namespace GeographicLib {
00033 
00034   using namespace std;
00035 
00036   // This is the transfer matrix for a 3rd order fit with a 12-point stencil
00037   // with weights
00038   //
00039   //  \ x -1  0  1  2
00040   //  y
00041   //  -1   .  1  1  .
00042   //   0   1  2  2  1
00043   //   1   1  2  2  1
00044   //   2   .  1  1  .
00045   //
00046   // A algorithm for n-dimensional polynomial fits is described in
00047   //   F. H. Lesh,
00048   //   Multi-dimensional least-squares polynomial curve fitting,
00049   //   CACM 2, 29-30 (1959).
00050   //
00051   // Here's the Maxima code to generate this matrix:
00052   //
00053   // /* The stencil and the weights */
00054   // xarr:[
00055   //     0, 1,
00056   // -1, 0, 1, 2,
00057   // -1, 0, 1, 2,
00058   //     0, 1]$
00059   // yarr:[
00060   //   -1,-1,
00061   // 0, 0, 0, 0,
00062   // 1, 1, 1, 1,
00063   //    2, 2]$
00064   // warr:[
00065   //    1, 1,
00066   // 1, 2, 2, 1,
00067   // 1, 2, 2, 1,
00068   //    1, 1]$
00069   //
00070   // /* [x exponent, y exponent] for cubic fit */
00071   // pows:[
00072   // [0,0],
00073   // [1,0],[0,1],
00074   // [2,0],[1,1],[0,2],
00075   // [3,0],[2,1],[1,2],[0,3]]$
00076   //
00077   // basisvec(x,y,pows):=map(lambda([ex],(if ex[1]=0 then 1 else x^ex[1])*
00078   //     (if ex[2]=0 then 1 else y^ex[2])),pows)$
00079   // addterm(x,y,f,w,pows):=block([a,b,bb:basisvec(x,y,pows)],
00080   //   a:w*(transpose(bb).bb),
00081   //   b:(w*f) * bb,
00082   //   [a,b])$
00083   //
00084   // c3row(k):=block([a,b,c,pows:pows,n],
00085   //   n:length(pows),
00086   //   a:zeromatrix(n,n),
00087   //   b:copylist(part(a,1)),
00088   //   c:[a,b],
00089   //   for i:1 thru length(xarr) do
00090   //   c:c+addterm(xarr[i],yarr[i],if i=k then 1 else 0,warr[i],pows),
00091   //   a:c[1],b:c[2],
00092   //   part(transpose( a^^-1 . transpose(b)),1))$
00093   // c3:[]$
00094   // for k:1 thru length(warr) do c3:endcons(c3row(k),c3)$
00095   // c3:apply(matrix,c3)$
00096   // c0:part(ratsimp(
00097   // genmatrix(yc,1,length(warr)).abs(c3).genmatrix(yd,length(pows),1)),2)$
00098   // c3:c0*c3$
00099 
00100   const Math::real Geoid::c0 = 240; // Common denominator
00101   const Math::real Geoid::c3[stencilsize * nterms] = {
00102       9, -18, -88,    0,  96,   90,   0,   0, -60, -20,
00103      -9,  18,   8,    0, -96,   30,   0,   0,  60, -20,
00104       9, -88, -18,   90,  96,    0, -20, -60,   0,   0,
00105     186, -42, -42, -150, -96, -150,  60,  60,  60,  60,
00106      54, 162, -78,   30, -24,  -90, -60,  60, -60,  60,
00107      -9, -32,  18,   30,  24,    0,  20, -60,   0,   0,
00108      -9,   8,  18,   30, -96,    0, -20,  60,   0,   0,
00109      54, -78, 162,  -90, -24,   30,  60, -60,  60, -60,
00110     -54,  78,  78,   90, 144,   90, -60, -60, -60, -60,
00111       9,  -8, -18,  -30, -24,    0,  20,  60,   0,   0,
00112      -9,  18, -32,    0,  24,   30,   0,   0, -60,  20,
00113       9, -18,  -8,    0, -24,  -30,   0,   0,  60,  20,
00114   };
00115 
00116   // Like c3, but with the coeffs of x, x^2, and x^3 constrained to be zero.
00117   // Use this at the N pole so that the height in independent of the longitude
00118   // there.
00119   //
00120   // Here's the Maxima code to generate this matrix (continued from above).
00121   //
00122   // /* figure which terms to exclude so that fit is indep of x at y=0 */
00123   // mask:part(zeromatrix(1,length(pows)),1)+1$
00124   // for i:1 thru length(pows) do
00125   // if pows[i][1]>0 and pows[i][2]=0 then mask[i]:0$
00126   //
00127   // /* Same as c3row but with masked pows. */
00128   // c3nrow(k):=block([a,b,c,powsa:[],n,d,e],
00129   //   for i:1 thru length(mask) do if mask[i]>0 then
00130   //   powsa:endcons(pows[i],powsa),
00131   //   n:length(powsa),
00132   //   a:zeromatrix(n,n),
00133   //   b:copylist(part(a,1)),
00134   //   c:[a,b],
00135   //   for i:1 thru length(xarr) do
00136   //   c:c+addterm(xarr[i],yarr[i],if i=k then 1 else 0,warr[i],powsa),
00137   //   a:c[1],b:c[2],
00138   //   d:part(transpose( a^^-1 . transpose(b)),1),
00139   //   e:[],
00140   //   for i:1 thru length(mask) do
00141   //   if mask[i]>0 then (e:endcons(first(d),e),d:rest(d)) else e:endcons(0,e),
00142   //   e)$
00143   // c3n:[]$
00144   // for k:1 thru length(warr) do c3n:endcons(c3nrow(k),c3n)$
00145   // c3n:apply(matrix,c3n)$
00146   // c0n:part(ratsimp(
00147   //     genmatrix(yc,1,length(warr)).abs(c3n).genmatrix(yd,length(pows),1)),2)$
00148   // c3n:c0n*c3n$
00149 
00150   const Math::real Geoid::c0n = 372; // Common denominator
00151   const Math::real Geoid::c3n[stencilsize * nterms] = {
00152       0, 0, -131, 0,  138,  144, 0,   0, -102, -31,
00153       0, 0,    7, 0, -138,   42, 0,   0,  102, -31,
00154      62, 0,  -31, 0,    0,  -62, 0,   0,    0,  31,
00155     124, 0,  -62, 0,    0, -124, 0,   0,    0,  62,
00156     124, 0,  -62, 0,    0, -124, 0,   0,    0,  62,
00157      62, 0,  -31, 0,    0,  -62, 0,   0,    0,  31,
00158       0, 0,   45, 0, -183,   -9, 0,  93,   18,   0,
00159       0, 0,  216, 0,   33,   87, 0, -93,   12, -93,
00160       0, 0,  156, 0,  153,   99, 0, -93,  -12, -93,
00161       0, 0,  -45, 0,   -3,    9, 0,  93,  -18,   0,
00162       0, 0,  -55, 0,   48,   42, 0,   0,  -84,  31,
00163       0, 0,   -7, 0,  -48,  -42, 0,   0,   84,  31,
00164   };
00165 
00166   // Like c3n, but y -> 1-y so that h is independent of x at y = 1.  Use this
00167   // at the S pole so that the height in independent of the longitude there.
00168   //
00169   // Here's the Maxima code to generate this matrix (continued from above).
00170   //
00171   // /* Transform c3n to c3s by transforming y -> 1-y */
00172   // vv:[
00173   //      v[11],v[12],
00174   // v[7],v[8],v[9],v[10],
00175   // v[3],v[4],v[5],v[6],
00176   //      v[1],v[2]]$
00177   // poly:expand(vv.(c3n/c0n).transpose(basisvec(x,1-y,pows)))$
00178   // c3sf[i,j]:=coeff(coeff(coeff(poly,v[i]),x,pows[j][1]),y,pows[j][2])$
00179   // c3s:genmatrix(c3sf,length(vv),length(pows))$
00180   // c0s:part(ratsimp(
00181   //     genmatrix(yc,1,length(warr)).abs(c3s).genmatrix(yd,length(pows),1)),2)$
00182   // c3s:c0s*c3s$
00183 
00184   const Math::real Geoid::c0s = 372; // Common denominator
00185   const Math::real Geoid::c3s[stencilsize * nterms] = {
00186      18,  -36, -122,   0,  120,  135, 0,   0,  -84, -31,
00187     -18,   36,   -2,   0, -120,   51, 0,   0,   84, -31,
00188      36, -165,  -27,  93,  147,   -9, 0, -93,   18,   0,
00189     210,   45, -111, -93,  -57, -192, 0,  93,   12,  93,
00190     162,  141,  -75, -93, -129, -180, 0,  93,  -12,  93,
00191     -36,  -21,   27,  93,   39,    9, 0, -93,  -18,   0,
00192       0,    0,   62,   0,    0,   31, 0,   0,    0, -31,
00193       0,    0,  124,   0,    0,   62, 0,   0,    0, -62,
00194       0,    0,  124,   0,    0,   62, 0,   0,    0, -62,
00195       0,    0,   62,   0,    0,   31, 0,   0,    0, -31,
00196     -18,   36,  -64,   0,   66,   51, 0,   0, -102,  31,
00197      18,  -36,    2,   0,  -66,  -51, 0,   0,  102,  31,
00198   };
00199 
00200   Geoid::Geoid(const std::string& name, const std::string& path, bool cubic)
00201     : _name(name)
00202     , _dir(path)
00203     , _cubic(cubic)
00204     , _a( Constants::WGS84_a() )
00205     , _e2( (2 - 1/Constants::WGS84_r())/Constants::WGS84_r() )
00206     , _degree( Constants::degree() )
00207     , _eps( sqrt(numeric_limits<real>::epsilon()) ) {
00208     if (_dir.empty())
00209       _dir = GeoidPath();
00210     if (_dir.empty())
00211       _dir = DefaultPath();
00212     _filename = _dir + "/" + _name + ".pgm";
00213     _file.open(_filename.c_str(), ios::binary);
00214     if (!(_file.good()))
00215       throw GeographicErr("File not readable " + _filename);
00216     string s;
00217     if (!(getline(_file, s) && s == "P5"))
00218       throw GeographicErr("File not in PGM format " + _filename);
00219     _offset = numeric_limits<real>::max();
00220     _scale = 0;
00221     _maxerror = _rmserror = -1;
00222     _description = "NONE";
00223     _datetime = "UNKNOWN";
00224     while (getline(_file, s)) {
00225       if (s.empty())
00226         continue;
00227       if (s[0] == '#') {
00228         istringstream is(s);
00229         string commentid, key;
00230         if (!(is >> commentid >> key) || commentid != "#")
00231           continue;
00232         if (key == "Description" || key =="DateTime") {
00233           string::size_type p = s.find_first_not_of(" \t", is.tellg());
00234           if (p != string::npos)
00235             (key == "Description" ? _description : _datetime) = s.substr(p);
00236         } else if (key == "Offset") {
00237           if (!(is >> _offset))
00238             throw GeographicErr("Error reading offset " + _filename);
00239         } else if (key == "Scale") {
00240           if (!(is >> _scale))
00241             throw GeographicErr("Error reading scale " + _filename);
00242         } else if (key == (_cubic ? "MaxCubicError" : "MaxBilinearError")) {
00243           // It's not an error if the error can't be read
00244           is >> _maxerror;
00245         } else if (key == (_cubic ? "RMSCubicError" : "RMSBilinearError")) {
00246           // It's not an error if the error can't be read
00247           is >> _rmserror;
00248         }
00249       } else {
00250         istringstream is(s);
00251         if (!(is >> _width >> _height))
00252           throw GeographicErr("Error reading raster size " + _filename);
00253         break;
00254       }
00255     }
00256     {
00257       unsigned maxval;
00258       if (!(_file >> maxval))
00259         throw GeographicErr("Error reading maxval " + _filename);
00260       if (maxval != 0xffffu)
00261         throw GeographicErr("Maxval not equal to 2^16-1 " + _filename);
00262       // Add 1 for whitespace after maxval
00263       _datastart = (unsigned long long)(_file.tellg()) + 1ULL;
00264       _swidth = (unsigned long long)(_width);
00265     }
00266     if (_offset == numeric_limits<real>::max())
00267       throw GeographicErr("Offset not set " + _filename);
00268     if (_scale == 0)
00269       throw GeographicErr("Scale not set " + _filename);
00270     if (_scale < 0)
00271       throw GeographicErr("Scale must be positive " + _filename);
00272     if (_height < 2 || _width < 2)
00273       // Coarsest grid spacing is 180deg.
00274       throw GeographicErr("Raster size too small " + _filename);
00275     if (_width & 1)
00276       // This is so that longitude grids can be extended thru the poles.
00277       throw GeographicErr("Raster width is odd " + _filename);
00278     if (!(_height & 1))
00279       // This is so that latitude grid includes the equator.
00280       throw GeographicErr("Raster height is even " + _filename);
00281     _file.seekg(0, ios::end);
00282     if (!_file.good() ||
00283         _datastart + 2ULL * _swidth * (unsigned long long)(_height) !=
00284         (unsigned long long)(_file.tellg()))
00285       // Possibly this test should be "<" because the file contains, e.g., a
00286       // second image.  However, for now we are more strict.
00287       throw GeographicErr("File has the wrong length " + _filename);
00288     _rlonres = _width / real(360);
00289     _rlatres = (_height - 1) / real(180);
00290     _cache = false;
00291     _ix = _width;
00292     _iy = _height;
00293     // Ensure that file errors throw exceptions
00294     _file.exceptions(ifstream::eofbit | ifstream::failbit | ifstream::badbit);
00295   }
00296 
00297   Math::real Geoid::height(real lat, real lon, bool gradp,
00298                            real& gradn, real& grade) const {
00299     real
00300       fx =  lon * _rlonres,
00301       fy = -lat * _rlatres;
00302     int
00303       ix = int(floor(fx)),
00304       iy = min((_height - 1)/2 - 1, int(floor(fy)));
00305     fx -= ix;
00306     fy -= iy;
00307     iy += (_height - 1)/2;
00308     ix += ix < 0 ? _width : ix >= _width ? -_width : 0;
00309     if (!(ix == _ix && iy == _iy)) {
00310       _ix = ix;
00311       _iy = iy;
00312       if (!_cubic) {
00313         _v00 = rawval(ix    , iy    );
00314         _v01 = rawval(ix + 1, iy    );
00315         _v10 = rawval(ix    , iy + 1);
00316         _v11 = rawval(ix + 1, iy + 1);
00317       } else {
00318         real v[stencilsize];
00319         int k = 0;
00320         v[k++] = rawval(ix    , iy - 1);
00321         v[k++] = rawval(ix + 1, iy - 1);
00322         v[k++] = rawval(ix - 1, iy    );
00323         v[k++] = rawval(ix    , iy    );
00324         v[k++] = rawval(ix + 1, iy    );
00325         v[k++] = rawval(ix + 2, iy    );
00326         v[k++] = rawval(ix - 1, iy + 1);
00327         v[k++] = rawval(ix    , iy + 1);
00328         v[k++] = rawval(ix + 1, iy + 1);
00329         v[k++] = rawval(ix + 2, iy + 1);
00330         v[k++] = rawval(ix    , iy + 2);
00331         v[k++] = rawval(ix + 1, iy + 2);
00332 
00333         const real* c3x = iy == 0 ? c3n : iy == _height - 2 ? c3s : c3;
00334         real c0x = iy == 0 ? c0n : iy == _height - 2 ? c0s : c0;
00335         for (unsigned i = 0; i < nterms; ++i) {
00336           _t[i] = 0;
00337           for (unsigned j = 0; j < stencilsize; ++j)
00338             _t[i] += v[j] * c3x[nterms * j + i];
00339           _t[i] /= c0x;
00340         }
00341       }
00342     }
00343     if (!_cubic) {
00344       real
00345         a = (1 - fx) * _v00 + fx * _v01,
00346         b = (1 - fx) * _v10 + fx * _v11,
00347         c = (1 - fy) * a + fy * b,
00348         h = _offset + _scale * c;
00349       if (gradp) {
00350         real
00351           phi = lat * _degree,
00352           cosphi = cos(phi),
00353           sinphi = sin(phi),
00354           n = 1/sqrt(1 - _e2 * sinphi * sinphi);
00355         gradn = ((1 - fx) * (_v00 - _v10) + fx * (_v01 - _v11)) *
00356           _rlatres / (_degree * _a * (1 - _e2) * n * n * n);
00357         grade = (cosphi > _eps ?
00358                  ((1 - fy) * (_v01 - _v00) + fy * (_v11 - _v10)) /   cosphi :
00359                  (sinphi > 0 ? _v11 - _v10 : _v01 - _v00) *
00360                  _rlatres / _degree ) *
00361           _rlonres / (_degree * _a * n);
00362         gradn *= _scale;
00363         grade *= _scale;
00364       }
00365       return h;
00366     } else {
00367       real h = _t[0] + fx * (_t[1] + fx * (_t[3] + fx * _t[6])) +
00368         fy * (_t[2] + fx * (_t[4] + fx * _t[7]) +
00369              fy * (_t[5] + fx * _t[8] + fy * _t[9]));
00370       h = _offset + _scale * h;
00371       if (gradp) {
00372         // Avoid 0/0 at the poles by backing off 1/100 of a cell size
00373         lat = min(lat,  90 - 1/(100 * _rlatres));
00374         lat = max(lat, -90 + 1/(100 * _rlatres));
00375         fy = (90 - lat) * _rlatres;
00376         fy -=  int(fy);
00377         real
00378           phi = lat * _degree,
00379           cosphi = cos(phi),
00380           sinphi = sin(phi),
00381           n = 1/sqrt(1 - _e2 * sinphi * sinphi);
00382         gradn = _t[2] + fx * (_t[4] + fx * _t[7]) +
00383           fy * (2 * _t[5] + fx * 2 * _t[8] + 3 * fy * _t[9]);
00384         grade = _t[1] + fx * (2 * _t[3] + fx * 3 * _t[6]) +
00385           fy * (_t[4] + fx * 2 * _t[7] + fy * _t[8]);
00386         gradn *= - _rlatres / (_degree * _a * (1 - _e2) * n * n * n) * _scale;
00387         grade *= _rlonres / (_degree * _a * n * cosphi) * _scale;
00388       }
00389       return h;
00390     }
00391   }
00392 
00393   void Geoid::CacheClear() const throw() {
00394     _cache = false;
00395     try {
00396       _data.clear();
00397       // Use swap to release memory back to system
00398       vector< vector<unsigned short> >().swap(_data);
00399     }
00400     catch (const exception&) {
00401     }
00402   }
00403 
00404   void Geoid::CacheArea(real south, real west, real north, real east) const {
00405     if (south > north) {
00406       CacheClear();
00407       return;
00408     }
00409     if (west >= 180)
00410       west -= 360;              // west in [-180, 180)
00411     if (east >= 180)
00412       east -= 360;
00413     if (east <= west)
00414       east += 360;              // east - west in (0, 360]
00415     int
00416       iw = int(floor(west * _rlonres)),
00417       ie = int(floor(east * _rlonres)),
00418       in = int(floor(-north * _rlatres)) + (_height - 1)/2,
00419       is = int(floor(-south * _rlatres)) + (_height - 1)/2;
00420     in = max(0, min(_height - 2, in));
00421     is = max(0, min(_height - 2, is));
00422     is += 1;
00423     ie += 1;
00424     if (_cubic) {
00425       in -= 1;
00426       is += 1;
00427       iw -= 1;
00428       ie += 1;
00429     }
00430     if (ie - iw >= _width - 1) {
00431       // Include entire longitude range
00432       iw = 0;
00433       ie = _width - 1;
00434     } else {
00435       ie += iw < 0 ? _width : iw >= _width ? -_width : 0;
00436       iw += iw < 0 ? _width : iw >= _width ? -_width : 0;
00437     }
00438     int oysize = int(_data.size());
00439     _xsize = ie - iw + 1;
00440     _ysize = is - in + 1;
00441     _xoffset = iw;
00442     _yoffset = in;
00443 
00444     try {
00445       _data.resize(_ysize, vector<unsigned short>(_xsize));
00446       for (int iy = min(oysize, _ysize); iy--;)
00447         _data[iy].resize(_xsize);
00448     }
00449     catch (const bad_alloc&) {
00450       CacheClear();
00451       throw GeographicErr("Insufficient memory for caching " + _filename);
00452     }
00453 
00454     try {
00455       vector<char> buf(2 * _xsize);
00456       for (int iy = in; iy <= is; ++iy) {
00457         int iy1 = iy, iw1 = iw;
00458         if (iy < 0 || iy >= _height) {
00459           iy1 = iy1 < 0 ? -iy1 : 2 * (_height - 1) - iy1;
00460           iw1 += _width/2;
00461           if (iw1 >= _width)
00462             iw1 -= _width;
00463         }
00464         int xs1 = min(_width - iw1, _xsize);
00465         filepos(iw1, iy1);
00466         _file.read(&(buf[0]), 2 * xs1);
00467         if (xs1 < _xsize) {
00468           // Wrap around longitude = 0
00469           filepos(0, iy1);
00470           _file.read(&(buf[2 * xs1]), 2 * (_xsize - xs1));
00471         }
00472         for (int ix = 0; ix < _xsize; ++ix)
00473           _data[iy - in][ix] =
00474             (unsigned short)((unsigned char)buf[2 * ix] * 256u +
00475                              (unsigned char)buf[2 * ix + 1]);
00476       }
00477       _cache = true;
00478     }
00479     catch (const exception& e) {
00480       CacheClear();
00481       throw GeographicErr(string("Error filling cache ") + e.what());
00482     }
00483   }
00484 
00485   std::string Geoid::DefaultPath() {
00486     return string(GEOID_DEFAULT_PATH);
00487   }
00488 
00489   std::string Geoid::GeoidPath() {
00490     string path;
00491     char* geoidpath = getenv("GEOID_PATH");
00492     if (geoidpath)
00493       path = string(geoidpath);
00494     return path;
00495   }
00496 }

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1