14 #if !defined(GEOGRAPHICLIB_DATA)
16 # define GEOGRAPHICLIB_DATA "C:/ProgramData/GeographicLib"
18 # define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib"
22 #if !defined(GEOGRAPHICLIB_GEOID_DEFAULT_NAME)
23 # define GEOGRAPHICLIB_GEOID_DEFAULT_NAME "egm96-5"
28 # pragma warning (disable: 4996)
99 const int Geoid::c0_ = 240;
100 const int Geoid::c3_[stencilsize_ * nterms_] = {
101 9, -18, -88, 0, 96, 90, 0, 0, -60, -20,
102 -9, 18, 8, 0, -96, 30, 0, 0, 60, -20,
103 9, -88, -18, 90, 96, 0, -20, -60, 0, 0,
104 186, -42, -42, -150, -96, -150, 60, 60, 60, 60,
105 54, 162, -78, 30, -24, -90, -60, 60, -60, 60,
106 -9, -32, 18, 30, 24, 0, 20, -60, 0, 0,
107 -9, 8, 18, 30, -96, 0, -20, 60, 0, 0,
108 54, -78, 162, -90, -24, 30, 60, -60, 60, -60,
109 -54, 78, 78, 90, 144, 90, -60, -60, -60, -60,
110 9, -8, -18, -30, -24, 0, 20, 60, 0, 0,
111 -9, 18, -32, 0, 24, 30, 0, 0, -60, 20,
112 9, -18, -8, 0, -24, -30, 0, 0, 60, 20,
149 const int Geoid::c0n_ = 372;
150 const int Geoid::c3n_[stencilsize_ * nterms_] = {
151 0, 0, -131, 0, 138, 144, 0, 0, -102, -31,
152 0, 0, 7, 0, -138, 42, 0, 0, 102, -31,
153 62, 0, -31, 0, 0, -62, 0, 0, 0, 31,
154 124, 0, -62, 0, 0, -124, 0, 0, 0, 62,
155 124, 0, -62, 0, 0, -124, 0, 0, 0, 62,
156 62, 0, -31, 0, 0, -62, 0, 0, 0, 31,
157 0, 0, 45, 0, -183, -9, 0, 93, 18, 0,
158 0, 0, 216, 0, 33, 87, 0, -93, 12, -93,
159 0, 0, 156, 0, 153, 99, 0, -93, -12, -93,
160 0, 0, -45, 0, -3, 9, 0, 93, -18, 0,
161 0, 0, -55, 0, 48, 42, 0, 0, -84, 31,
162 0, 0, -7, 0, -48, -42, 0, 0, 84, 31,
183 const int Geoid::c0s_ = 372;
184 const int Geoid::c3s_[stencilsize_ * nterms_] = {
185 18, -36, -122, 0, 120, 135, 0, 0, -84, -31,
186 -18, 36, -2, 0, -120, 51, 0, 0, 84, -31,
187 36, -165, -27, 93, 147, -9, 0, -93, 18, 0,
188 210, 45, -111, -93, -57, -192, 0, 93, 12, 93,
189 162, 141, -75, -93, -129, -180, 0, 93, -12, 93,
190 -36, -21, 27, 93, 39, 9, 0, -93, -18, 0,
191 0, 0, 62, 0, 0, 31, 0, 0, 0, -31,
192 0, 0, 124, 0, 0, 62, 0, 0, 0, -62,
193 0, 0, 124, 0, 0, 62, 0, 0, 0, -62,
194 0, 0, 62, 0, 0, 31, 0, 0, 0, -31,
195 -18, 36, -64, 0, 66, 51, 0, 0, -102, 31,
196 18, -36, 2, 0, -66, -51, 0, 0, 102, 31,
199 Geoid::Geoid(
const std::string& name,
const std::string& path,
bool cubic,
206 , _degree(
Math::degree() )
207 , _eps( sqrt(numeric_limits<real>::epsilon()) )
210 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(pixel_t) == pixel_size_,
211 "pixel_t has the wrong size");
214 _filename = _dir +
"/" + _name + (pixel_size_ != 4 ?
".pgm" :
".pgm4");
215 _file.open(_filename.c_str(), ios::binary);
219 if (!(getline(_file, s) && s ==
"P5"))
221 _offset = numeric_limits<real>::max();
223 _maxerror = _rmserror = -1;
224 _description =
"NONE";
225 _datetime =
"UNKNOWN";
226 while (getline(_file, s)) {
231 string commentid, key;
232 if (!(is >> commentid >> key) || commentid !=
"#")
234 if (key ==
"Description" || key ==
"DateTime") {
235 string::size_type p =
236 s.find_first_not_of(
" \t",
unsigned(is.tellg()));
237 if (p != string::npos)
238 (key ==
"Description" ? _description : _datetime) = s.substr(p);
239 }
else if (key ==
"Offset") {
240 if (!(is >> _offset))
242 }
else if (key ==
"Scale") {
245 }
else if (key == (_cubic ?
"MaxCubicError" :
"MaxBilinearError")) {
248 }
else if (key == (_cubic ?
"RMSCubicError" :
"RMSBilinearError")) {
254 if (!(is >> _width >> _height))
255 throw GeographicErr(
"Error reading raster size " + _filename);
261 if (!(_file >> maxval))
263 if (maxval != pixel_max_)
264 throw GeographicErr(
"Incorrect value of maxval " + _filename);
266 _datastart = (
unsigned long long)(_file.tellg()) + 1ULL;
267 _swidth = (
unsigned long long)(_width);
269 if (_offset == numeric_limits<real>::max())
275 if (_height < 2 || _width < 2)
284 _file.seekg(0, ios::end);
286 _datastart + pixel_size_ * _swidth * (
unsigned long long)(_height) !=
287 (
unsigned long long)(_file.tellg()))
290 throw GeographicErr(
"File has the wrong length " + _filename);
291 _rlonres = _width / real(360);
292 _rlatres = (_height - 1) / real(180);
297 _file.exceptions(ifstream::eofbit | ifstream::failbit | ifstream::badbit);
314 fy = -lat * _rlatres;
317 iy = min((_height - 1)/2 - 1,
int(floor(fy)));
320 iy += (_height - 1)/2;
321 ix += ix < 0 ? _width : (ix >= _width ? -_width : 0);
322 real v00 = 0, v01 = 0, v10 = 0, v11 = 0;
325 if (_threadsafe || !(ix == _ix && iy == _iy)) {
327 v00 = rawval(ix , iy );
328 v01 = rawval(ix + 1, iy );
329 v10 = rawval(ix , iy + 1);
330 v11 = rawval(ix + 1, iy + 1);
332 real v[stencilsize_];
334 v[k++] = rawval(ix , iy - 1);
335 v[k++] = rawval(ix + 1, iy - 1);
336 v[k++] = rawval(ix - 1, iy );
337 v[k++] = rawval(ix , iy );
338 v[k++] = rawval(ix + 1, iy );
339 v[k++] = rawval(ix + 2, iy );
340 v[k++] = rawval(ix - 1, iy + 1);
341 v[k++] = rawval(ix , iy + 1);
342 v[k++] = rawval(ix + 1, iy + 1);
343 v[k++] = rawval(ix + 2, iy + 1);
344 v[k++] = rawval(ix , iy + 2);
345 v[k++] = rawval(ix + 1, iy + 2);
347 const int* c3x = iy == 0 ? c3n_ : (iy == _height - 2 ? c3s_ : c3_);
348 int c0x = iy == 0 ? c0n_ : (iy == _height - 2 ? c0s_ : c0_);
349 for (
unsigned i = 0; i < nterms_; ++i) {
351 for (
unsigned j = 0; j < stencilsize_; ++j)
352 t[i] += v[j] * c3x[nterms_ * j + i];
363 copy(_t, _t + nterms_, t);
367 a = (1 - fx) * v00 + fx * v01,
368 b = (1 - fx) * v10 + fx * v11,
369 c = (1 - fy) * a + fy * b,
370 h = _offset + _scale * c;
376 n = 1/sqrt(1 - _e2 * sinphi * sinphi);
377 gradn = ((1 - fx) * (v00 - v10) + fx * (v01 - v11)) *
378 _rlatres / (_degree * _a * (1 - _e2) * n * n * n);
379 grade = (cosphi > _eps ?
380 ((1 - fy) * (v01 - v00) + fy * (v11 - v10)) / cosphi :
381 (sinphi > 0 ? v11 - v10 : v01 - v00) *
382 _rlatres / _degree ) *
383 _rlonres / (_degree * _a * n);
397 real h = t[0] + fx * (t[1] + fx * (t[3] + fx * t[6])) +
398 fy * (t[2] + fx * (t[4] + fx * t[7]) +
399 fy * (t[5] + fx * t[8] + fy * t[9]));
400 h = _offset + _scale * h;
403 lat = min(lat, 90 - 1/(100 * _rlatres));
404 lat = max(lat, -90 + 1/(100 * _rlatres));
405 fy = (90 - lat) * _rlatres;
411 n = 1/sqrt(1 - _e2 * sinphi * sinphi);
412 gradn = t[2] + fx * (t[4] + fx * t[7]) +
413 fy * (2 * t[5] + fx * 2 * t[8] + 3 * fy * t[9]);
414 grade = t[1] + fx * (2 * t[3] + fx * 3 * t[6]) +
415 fy * (t[4] + fx * 2 * t[7] + fy * t[8]);
416 gradn *= - _rlatres / (_degree * _a * (1 - _e2) * n * n * n) * _scale;
417 grade *= _rlonres / (_degree * _a * n * cosphi) * _scale;
422 copy(t, t + nterms_, _t);
434 vector< vector<pixel_t> >().swap(_data);
436 catch (
const exception&) {
443 throw GeographicErr(
"Attempt to change cache of threadsafe Geoid");
453 iw = int(floor(west * _rlonres)),
454 ie = int(floor(east * _rlonres)),
455 in = int(floor(-north * _rlatres)) + (_height - 1)/2,
456 is =
int(floor(-south * _rlatres)) + (_height - 1)/2;
457 in = max(0, min(_height - 2, in));
458 is = max(0, min(_height - 2, is));
467 if (ie - iw >= _width - 1) {
472 ie += iw < 0 ? _width : (iw >= _width ? -_width : 0);
473 iw += iw < 0 ? _width : (iw >= _width ? -_width : 0);
475 int oysize = int(_data.size());
476 _xsize = ie - iw + 1;
477 _ysize = is - in + 1;
482 _data.resize(_ysize, vector<pixel_t>(_xsize));
483 for (
int iy = min(oysize, _ysize); iy--;)
484 _data[iy].resize(_xsize);
486 catch (
const bad_alloc&) {
488 throw GeographicErr(
"Insufficient memory for caching " + _filename);
492 for (
int iy = in; iy <= is; ++iy) {
493 int iy1 = iy, iw1 = iw;
494 if (iy < 0 || iy >= _height) {
496 iy1 = iy1 < 0 ? -iy1 : 2 * (_height - 1) - iy1;
501 int xs1 = min(_width - iw1, _xsize);
503 Utility::readarray<pixel_t, pixel_t, true>
504 (_file, &(_data[iy - in][0]), xs1);
508 Utility::readarray<pixel_t, pixel_t, true>
509 (_file, &(_data[iy - in][xs1]), _xsize - xs1);
514 catch (
const exception& e) {
516 throw GeographicErr(
string(
"Error filling cache ") + e.what());
522 char* geoidpath = getenv(
"GEOGRAPHICLIB_GEOID_PATH");
524 path = string(geoidpath);
527 char* datapath = getenv(
"GEOGRAPHICLIB_DATA");
529 path = string(datapath);
535 char* geoidname = getenv(
"GEOGRAPHICLIB_GEOID_NAME");
537 name = string(geoidname);
static T AngNormalize(T x)
GeographicLib::Math::real real
Header for GeographicLib::Utility class.
Mathematical functions needed by GeographicLib.
static std::string DefaultGeoidPath()
#define GEOGRAPHICLIB_GEOID_DEFAULT_NAME
void CacheArea(real south, real west, real north, real east) const
#define GEOGRAPHICLIB_DATA
Namespace for GeographicLib.
Constants needed by GeographicLib
Exception handling for GeographicLib.
static std::string DefaultGeoidName()
Header for GeographicLib::Geoid class.