00001
00002
00003
00004
00005
00006
00007
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
00029 #pragma warning (disable: 4996)
00030 #endif
00031
00032 namespace GeographicLib {
00033
00034 using namespace std;
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100 const Math::real Geoid::c0 = 240;
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
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150 const Math::real Geoid::c0n = 372;
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
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184 const Math::real Geoid::c0s = 372;
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
00244 is >> _maxerror;
00245 } else if (key == (_cubic ? "RMSCubicError" : "RMSBilinearError")) {
00246
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
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
00274 throw GeographicErr("Raster size too small " + _filename);
00275 if (_width & 1)
00276
00277 throw GeographicErr("Raster width is odd " + _filename);
00278 if (!(_height & 1))
00279
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
00286
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
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
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
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;
00411 if (east >= 180)
00412 east -= 360;
00413 if (east <= west)
00414 east += 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
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
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 }