// calculation of geoid bby usnig EGM08 by refer to the solution of openterrain based on NGA // using 2.5m x 2.5m interval and indicae by interpolation // using 116 value of geoid #include #include #include #include // for nan #include #include #ifndef NAN # define NAN (std::numeric_limits::quiet_NaN()) #endif // check whether of not a floating point value is not a number inline int isNAN(const double v) { return(v != v); } //convert from deg min sec to decimal degree float convert(float deg, float min, float sec); // bi-linearily sample data matrix at normalized position (x,y) inline float interpolate(float x, float y, const float *data, unsigned int width, unsigned int height); float geoid_heigt(double latitude, double longitude); // main function int main() { float dlat, mlat, slat, dlon, mlon, slon; std::cout << "\nplease insert the latitude ( ex: 8 10 20) : "; std::cin >> dlat>>mlat>>slat; std::cout << "\nplease insert the longitude ( ex: 100 10 10) : "; std::cin >> dlon>>mlon>>slon; float flat = convert(dlat, mlat, slat); float flon = convert(dlon, mlon, slon); float lat = abs(flat); float lon = abs(flon); if ((dlat >= 0) && (dlon >= 0)) { printf("\nThe geoid height at %g N %g E", lat, lon); } else if ((dlat >= 0) && (dlon < 0)) { printf("\nThe geoid height at %g N %g W", lat, lon); } else if ((dlat < 0) && (dlon < 0)) { printf("\nThe geoid height at %g S %g W", lat, lon); } else { printf("\nThe geoid height at %g S %g E", lat, lon); } printf(" is %g meter \n", geoid_heigt(lat, lon)); system("pause"); return(0); } //convert from deg min sec to decimal degree float convert(float deg, float min, float sec){ /* deg min sec -> dd d=deg m=min/60 s=sec/3600 dd=d+m+s */ float d, m, s, dd; d = deg; m = min / 60; s = sec / 3600; dd = d + m + s; return dd; } // bi-linearily sample data matrix at normalized position (x,y) inline float interpolate(float x, float y, const float *data, unsigned int width, unsigned int height) { unsigned int i, j; float val1, val2, val3, val4; const float *floatptr; if (x >= 0.0f && x <= 1.0f && y >= 0.0f && y <= 1.0f) { x *= width - 1; y *= height - 1; i = (unsigned int)(x + 0.5f); j = (unsigned int)(y + 0.5f); x -= i; y -= j; if (i == width - 1 && width > 1) { i = width - 2; x = 1.0f; } if (j == height - 1 && height > 1) { j = height - 2; y = 1.0f; } floatptr = &data[i + j * width]; val1 = floatptr[0]; if (isNAN(val1)) return(NAN); if (width < 2) return(val1); else { val2 = floatptr[1]; if (isNAN(val2)) return(NAN); if (height < 2) return((1.0f - x)*val1 + x * val2); else { val3 = floatptr[width]; val4 = floatptr[width + 1]; if (isNAN(val3) || isNAN(val4)) return(NAN); return((1.0f - y)*((1.0f - x)*val1 + x * val2) + y * ((1.0f - x)*val3 + x * val4)); } } } return(NAN); } // value of geoid static const int size_x=11; static const int size_y=6; static const double LLWGS84_swx=0; static const double LLWGS84_swy=-90; static const double LLWGS84_nwx=0; static const double LLWGS84_nwy=90; static const double LLWGS84_nex=360; static const double LLWGS84_ney=90; static const double LLWGS84_sex=360; static const double LLWGS84_sey=-90; static const float data[size_y][size_x]={ -30.15,-30.15,-30.15,-30.15,-30.15,-30.15,-30.15,-30.15,-30.15,-30.15,-30.15, 22.8511,40.7621,22.669,-11.7373,-11.8084,-7.656,-12.5149,-8.93637,0.880717,0.00731354,22.8511, 17.0403,-33.2788,-73.6445,15.6048,45.5325,27.1284,5.36011,-8.95659,-6.25443,-17.5287,17.0403, 18.9033,-25.0505,-76.6696,13.3805,39.3692,20.335,-0.984197,-13.8498,-0.883649,-20.8845,18.9033, 47.2155,10.0123,-37.2059,-3.50886,2.8773,-1.61806,-7.87575,-27.4857,-22.7701,37.2808,47.2155, 14.8985,14.8985,14.8985,14.8985,14.8985,14.8985,14.8985,14.8985,14.8985,14.8985,14.8985 }; // sample the geoid height at a geographic position float geoid_heigt(double latitude, double longitude) { // calculate normalized coordinates from geographic position double nx = (longitude - LLWGS84_swx) / (LLWGS84_sex - LLWGS84_swx); double ny = (latitude - LLWGS84_swy) / (LLWGS84_nwy - LLWGS84_swy); // adjust wraparound longitudes nx = nx - floor(nx); // clamp out-of-range latitudes if (ny < 0.0) ny = 0.0; else if (ny > 1.0) ny = 1.0; // interpolate geoid height return(interpolate(nx, ny, (const float *)data, size_x, size_y)); }

calculator geoid height