Rash thoughts about .NET, C#, F# and Dynamics NAV.

"Every solution will only lead to new problems."

Friday, 5. December 2008

The distance of two cities on the surface of the earth in F#

Filed under: English posts,F#,Informatik — Steffen Forkmann at 12:19 Uhr

If we want to calculate the distance between two objects, we need some sort of distance function. I recently showed how we can compute the “edit distance” between two strings (see Damerau-Levenshtein-Distance in F# – part I, part II and part III). This time I will consider the distance between two points in the plane and on the surface of the earth.

First of all we need some sort of location type, which stores information about the coordinates (I used this type before). I will use the “Units of Measure”-feature in F# to distinct between <degree> vs. <radian> and <m> vs. <km>. If you are not familiar with this concept see Andrew Kennedy’s excellent blog series about Units of Measures in F#.

open Microsoft.FSharp.Math.SI
type degree

type radian =
  static member per_degree = Math.PI/180.0<degree/radian>
type km =  
  static member per_meter = 1000.<m/km>
type Location =
   Longitude:float<degree> }
  member x.YPos = x.Longitude
  member x.XPos = x.Latitude        
  member x.Longitude_radian = x.Longitude * radian.per_degree  
  member x.Latitude_radian = x.Latitude * radian.per_degree   
  member x.Coordinates = (x.Longitude,x.Latitude)  
  static member CreateLocation(latitude, longitude) = 
    {new Location with
      Latitude = latitude and
      Longitude = longitude}

Now we can easily compute the Euclidean distance for two points in the plane.

let sq x = x * x
let CalcEuclideanDistance (s:Location) (t:Location) =

let d1 = (t.Longitude - s.Longitude) / 1.<degree> let d2 = (t.Latitude - s.Latitude) / 1.<degree> ((sq d1 + sq d2) |> Math.Sqrt) * 1.<m>

If we want to approximate the distance of two cities on the surface of earth, we can use the “Great circle distance“.

“The great-circle distance is the shortest distance between any two points on the surface of a sphere measured along a path on the surface of the sphere (as opposed to going through the sphere’s interior)” – Wikipedia

let kmFactor = 6372.795<km>
let sin_rad x = Math.Sin(x / 1.<radian>)
let cos_rad x = Math.Cos(x / 1.<radian>)

/// Great Circle Distance 
let CalcGreatCircleDistance (s:Location) (t:Location)=
  if s = t then
    let factor = kmFactor * km.per_meter
    let b =
         sin_rad s.Latitude_radian * sin_rad t.Latitude_radian +
         cos_rad s.Latitude_radian * cos_rad t.Latitude_radian *
           cos_rad (t.Longitude_radian - s.Longitude_radian) 
    Math.Acos(b) * factor
Tags: , , ,