So I'm coding a program (c++), and thing leads to thing, and I need to find out how many meters to a degree of longitude at a specific latitude to a very high precision. Optimally < 1 meter at all latitudes(The program should be able to use this to plot many points given in lat and lon that are sub meter apart.) pretty much anywhere where there would normally be people (this means everywhere but within a few degrees of the poles).
When researching I came across some formulas to do this. The one that was most commonly available to me produced distances that were ~50 meters off for every 15 degrees of lat compared to what Wikipedia (and other online calcs) said they should be. The second formula seems to be more accurate, but I have no idea WHY its more accurate then the first one.
code/formulas:
double MetersAtLat(double lat)
{
long double pi = 3.14159265358979323846;
long double a = 6378137.0; //equatorial radius, meters exactly
long double b = 6356752.3142;// polar radius, meters exactly
long double eccentricitySquared = ((a*a) - (b*b))/(a*a);//0.00669437999014?;
long double rlat = (lat) * (pi/180);
//method 1; off by ~ 50m for every 15 degrees away from equator. (for some reason)
long double answ = (pi * a * cos(rlat))/180 * sqrt (1 - eccentricitySquared * (sin(rlat)*sin(rlat)));
//method 2; off by ~0.5m? at the equator and gets more accurate the more northern lats given?.
long double tanPsi = (b * tan(rlat))/a;
long double psi = atan(tanPsi);
long double answ2 = (pi * a * cos(psi))/180;
return answ2;
}
Is the problem with method 1 that it really is that inaccurate or is it really an issue with my computer/compiler/language? and is the second method as accurate as it looks? (are the numbers it returns good?)
edit: I've been using a table on Wikipedia to verify the numbers. (https://en.wikipedia.org/wiki/Latitude)
-
1Add some more parentheses on method 1 to make sure the */ are being done correctly. Possibly change 180 to 180.0. I couldn't get method 2 to work.mkennedy– mkennedy2016年01月07日 17:46:04 +00:00Commented Jan 7, 2016 at 17:46
1 Answer 1
As per mkennedy's suggestion ( thank you), I made some changes in the layout of my code, and now both methods appear to be producing the same answers. I'm making a bit of a leap of faith for the "high-accuracy/good numbers" part, only comparing them to a few web sources without actually "knowing for sure."
It appears that some floating point operations were getting in the way with how method 1 was originally being calculated, resulting in a loss of precision. Breaking it into separate steps appears to have avoided the issue.
Code Follows:
double viewManager::MetersAtLat1(double lat)
{
long double pi = 3.14159265358979323846;
long double a = 6378137.0; //equatorial radius, meters exactly
long double b = 6356752.3142;//polar radius, meters exactly
long double eccentricitySquared = ((a*a) - (b*b))/(a*a);//0.00669437999014?;
long double rlat = (lat) * (pi/180); //change to radians
long double topPart = (pi * a * cos(rlat));
long double sinLatSqr = (sin(rlat)*sin(rlat));
long double bottomPart = 180.0 * sqrt (1 - eccentricitySquared * sinLatSqr);
long double answ = topPart /bottomPart;
return answ;
}
double viewManager::MetersAtLat2(double lat)
{
long double pi = 3.14159265358979323846;
long double a = 6378137.0; //equatorial radius, meters exactly
long double b = 6356752.3142;//polar radius, meters exactly
long double eccentricitySquared = ((a*a) - (b*b))/(a*a);//0.00669437999014?;
long double rlat = (lat) * (pi/180); //change to radians
long double tanPsi = (b * tan(rlat))/a;
long double psi = atan(tanPsi);
long double answ = (pi * a * cos(psi))/180;
return answ;
}
output: both methods
Meters @ lat 0; 111319.49079327358 Meters @ lat 15; 107550.48871678663 Meters @ lat 30; 96486.280251067088 Meters @ lat 45; 78846.83509425736 Meters @ lat 60; 55800.001572733083 Meters @ lat 75; 28902.005809998191 Meters @ lat 90; 6.81e-12