3

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)

  • 1
    Add 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 Jan 7 '16 at 17:46
2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.