Great circle distance

Compute the great circle distance of 2 points

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
42: 
43: 
44: 
45: 
46: 
47: 
48: 
49: 
50: 
51: 
52: 
53: 
54:

open System

module Earth =
    type Point(lon,lat) = 
        let toRad deg = deg * (Math.PI / 180.0)
        
        member this.Lon = lon
        member this.Lat = lat
        
        member this.LonRad = toRad lon
        member this.LatRad = toRad lat

    
    let greatCircleDistance (p1:Point) (p2:Point) =
        // code adapted from 
        // http://www.codeproject.com/Articles/12269/Distance-between-locations-using-latitude-and-long
        (*
            The Haversine formula according to Dr. Math.
            http://mathforum.org/library/drmath/view/51879.html
                
            dlon = lon2 - lon1
            dlat = lat2 - lat1
            a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
            c = 2 * atan2(sqrt(a), sqrt(1-a)) 
            d = R * c
                
            Where
                * dlon is the change in longitude
                * dlat is the change in latitude
                * c is the great circle distance in Radians.
                * R is the radius of a spherical Earth.
                * The locations of the two points in 
                    spherical coordinates (longitude and 
                    latitude) are lon1,lat1 and lon2, lat2.
        *)
        
        let dlon = p2.LonRad - p1.LonRad;
        let dlat = p2.LatRad - p1.LatRad;

        // Intermediate result a.
        let a = (sin (dlat / 2.0)) ** 2.0 + ((cos p1.LatRad) * (cos p2.LatRad) * (sin (dlon / 2.0)) ** 2.0);

        // Intermediate result c (great circle distance in Radians).
        let c = 2.0 * (asin (sqrt a));

        // Distance.
        let earthRadiusKms = 6371.0;
        let distance = earthRadiusKms * c;

        distance

    let test = 
        let d = greatCircleDistance (new Point(5.0, -32.0)) (new Point(-3.0, 4.0))
        printfn "%f" d // 4091 km

namespace System

module Earth

from Script

Multiple items
type Point =
  new : lon:float * lat:float -> Point
  member Lat : float
  member LatRad : float
  member Lon : float
  member LonRad : float

Full name: Script.Earth.Point

--------------------
new : lon:float * lat:float -> Point

val lon : float

val lat : float

val toRad : (float -> float)

val deg : float

type Math =
  static val PI : float
  static val E : float
  static member Abs : value:sbyte -> sbyte + 6 overloads
  static member Acos : d:float -> float
  static member Asin : d:float -> float
  static member Atan : d:float -> float
  static member Atan2 : y:float * x:float -> float
  static member BigMul : a:int * b:int -> int64
  static member Ceiling : d:decimal -> decimal + 1 overload
  static member Cos : d:float -> float
  ...

Full name: System.Math

field Math.PI = 3.14159265359

val this : Point

member Point.Lon : float

Full name: Script.Earth.Point.Lon

member Point.Lat : float

Full name: Script.Earth.Point.Lat

member Point.LonRad : float

Full name: Script.Earth.Point.LonRad

member Point.LatRad : float

Full name: Script.Earth.Point.LatRad

val greatCircleDistance : p1:Point -> p2:Point -> float

Full name: Script.Earth.greatCircleDistance

val p1 : Point

val p2 : Point

val dlon : float

property Point.LonRad: float

val dlat : float

property Point.LatRad: float

val a : float

val sin : value:'T -> 'T (requires member Sin)

Full name: Microsoft.FSharp.Core.Operators.sin

val cos : value:'T -> 'T (requires member Cos)

Full name: Microsoft.FSharp.Core.Operators.cos

val c : float

val asin : value:'T -> 'T (requires member Asin)

Full name: Microsoft.FSharp.Core.Operators.asin

val sqrt : value:'T -> 'U (requires member Sqrt)

Full name: Microsoft.FSharp.Core.Operators.sqrt

val earthRadiusKms : float

val distance : float

val test : unit

Full name: Script.Earth.test

val d : float

val printfn : format:Printf.TextWriterFormat<'T> -> 'T

Full name: Microsoft.FSharp.Core.ExtraTopLevelOperators.printfn

Raw view Test code New version

More information

Link:	http://fssnip.net/jp
Posted:	10 years ago
Author:	Samuel Bosch
Tags:	math , spatial