[GIS] Calculating the earth radius at latitude

algorithmcoordinate systemjavascript

This is an extension to this question: How do you compute the earth's radius at a given geodetic latitude? and to the Wikipedia page: https://en.wikipedia.org/wiki/Earth_radius#Geocentric_radius

I tried implementing the function (in Javascript) to calculate the horizontal radius for a given latitude. But nobody specified the units and I am at a loss why my function simply outputs enormous values:

var WGS_ELLIPSOID = { a: 6378137.0, b: 6356752.314 }; // meter

function get_radius_at_lat(lat) {

    var f1 = Math.pow((Math.pow(WGS_ELLIPSOID.a, 2) * Math.cos(lat)), 2);
    var f2 = Math.pow((Math.pow(WGS_ELLIPSOID.b, 2) * Math.sin(lat)), 2);
    var f3 = Math.pow((WGS_ELLIPSOID.a * Math.cos(lat)), 2);
    var f4 = Math.pow((WGS_ELLIPSOID.b * Math.sin(lat)), 2);

    var radius =  Math.sqrt((f1 + f2) / (f3 / f4));

    return radius;
}

For something like 0 (equator) it returns 0 (incorrect). For something like 45 degree north, it returns an astronomial value of 65354531512255.51, which I am not sure what unit it has. Also, it's completely wrong – the radius of earth has a maximum of 6378137 m.

Meanwhile, if I do: 6371000.0 * Math.cos(lat), I get 3346826.3907577563 (supposedly meter?) for 45 degree north. Which is much more plausible, if the earth at the equator has a radius of 6000 km, 3300 for 45° north is OK, but this does not account for the earths flatness.

What am I doing wrong? Is there an error in the formula or some Javascript weirdness or are my units off? The input is in degrees, not radians (I assumed this because of the cos() functions).

The formula I used:

radius(lat) = sqrt( (a^2 cos(lat))^2 + (b^2 sin(lat))^2 ) / 
                    ( (a cos(lat))^2   + (b sin(lat))^2 ) )

Best Answer

You made an error in the line defining the radius in your function. Instead of f3 / f4 it should be f3 + f4. The corrected line should be as follows -

var radius =  Math.sqrt((f1 + f2) / (f3 + f4));

Also the Math.sin() and Math.cos() functions take the inputs a radians. So if your input for lat is in degrees then you would need to convert that. Add a line before defining f1 in the function -

lat  = lat * (Math.PI / 180);
Related Question