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