Computing the distance between two locations from coordinates in F#

The following code returns the distance between to locations based on each point’s longitude and latitude. The distance returned is relative to Earth’s radius. To get the distance in miles, multiply by 3960. To get the distance in kilometers, multiply by 6373.

Latitude is measured in degrees north of the equator; southern locations have negative latitude. Similarly, longitude is measured in degrees east of the Prime Meridian. A location 10° west of the Prime Meridian, for example, could be expressed as either 350° east or as -10° east.

The following F# implementation comes from Kevin Hazzard. See also Python implementation.

[<Measure>] type radian
[<Measure>] type degree
[<Measure>] type km
[<Measure>] type mi
type LatLong = { Lat : float<degree>; Long : float<degree> }
 
// a unit-of-measure (UoM) generic F# function based on 
// John D. Cook's Python version - https://www.johndcook.com
// a UoM generic suface distance 
let SphereSurfaceDistance[<Measure>] 't>
        (radius : float<'t>) (loc1 : LatLong) (loc2 : LatLong) =
    // Convert latitude and longitude to
    // spherical coordinates in radians.
    let degrees_to_radians (d : float<degree>) =
        d * System.Math.PI / 180.0<degree/radian>
 
    // phi = 90 - latitude (coded F# pipeline style)
    let phi1 = (90.0<degree> - loc1.Lat) |> degrees_to_radians
    let phi2 = (90.0<degree> - loc2.Lat) |> degrees_to_radians
 
    // theta = longitude (coded F# function call style)
    let theta1 = degrees_to_radians loc1.Long
    let theta2 = degrees_to_radians loc2.Long
 
    // Compute spherical distance from spherical coordinates.
    
    // For two locations in spherical coordinates      
    // (1, theta, phi) and (1, theta, phi)     
    // cosine( arc length ) =      
    //    sin phi sin phi' cos(theta-theta') + cos phi cos phi'     
    // distance = rho * arc length
    
    radius * acos(sin(float(phi1)) * sin(float(phi2)) *
        cos(float(theta1 - theta2)) + 
        cos(float(phi1)) * cos(float(phi2)))
 
// instantiate Earth-relative functions in km and miles
let SurfaceDistanceOnEarthKm = SphereSurfaceDistance 6371.0<km>
let SurfaceDistanceOnEarthMile = SphereSurfaceDistance 3959.0<mi>
 
let test1_SurfaceDistanceOnEarthMile =
    let richmond = { Lat = 37.542979<degree>; Long = -77.469092<degree> }
    let saopaulo = { Lat = -23.548943<degree>; Long = -46.638818<degree> }
 
    // distance between Richmond, Virginia USA and S o Paulo, Brazil
    // Google Earth says this should be about 4677.562 miles

    let tst_dist = SurfaceDistanceOnEarthMile richmond saopaulo
    let std_dist = 4677.562<mi>
    let delta = abs (tst_dist - std_dist) / std_dist
    // assert delta from standard < 1/4 of a percent
    delta < 0.0025
 
let test1_SurfaceDistanceOnEarthKm =
    let richmond = { Lat = 37.542979<degree>; Long = -77.469092<degree> }
    let saopaulo = { Lat = -23.548943<degree>; Long = -46.638818<degree> }
 
    // distance between Richmond, Virginia USA and S o Paulo, Brazil
    // Google Earth says this should be about 7527.806 km
    let tst_dist = SurfaceDistanceOnEarthKm richmond saopaulo
    let std_dist = 7527.806<km>
    let delta = abs (tst_dist - std_dist) / std_dist
    // assert delta from standard < 1/4 of a percent
    delta < 0.0025

The code above assumes the earth is perfectly spherical. For a discussion of how accurate this assumption is, see my blog post on the shape of the Earth.

The algorithm used to calculate distances is described in detail here.

A web page to calculate the distance between to cities based on longitude and latitude is available here.

This code is in the public domain. Do whatever you want with it, no strings attached.

More stand-alone numerical code

Click to find out more about consulting for numerical computing