Great Circle Distance

  1. 8 months ago
    Edited 8 months ago

    Can anyone suggest a more accurate set of functions to calculate great circle distance?

    Here's what I am using, but the result doesn't seem to be very accurate.

    input data: lat:|N |33 |45 |43.00 long:|W |84 |27 |29.00 |
    ===================

    getDistance( long As String, lat As String ) As Double

    dim d1, m1, s1 As Integer
    dim d2, m2, s2 As Integer
    Dim decLat, decLong, dist As Double

    dim a, b, c As String

    a = mid( long, 3, 4 )
    b = mid( long, 8, 2 )
    c = mid( long, 12, 7 )

    d1 = val( a )
    m1 = val( b )
    s1 = val( c )

    a = mid( lat, 3, 4 )
    b = mid( lat, 8, 2 )
    c = mid( lat, 12, 7 )

    d2 = val( a )
    m2 = val( b )
    s2 = val( c )

    decLong = decimalDeg( d1, m1, s1 )
    if Left( long, 1 ) = "W" Then decLong = decLong * -1
    decLat = decimalDeg( d2, m2, s2 )
    if Left( lat, 1 ) = "S" Then decLat = decLat * -1
    dist = greatCircleDistance( decLong, decLat )
    Return dist

    ========================

    getGreatCircleDistance( long As Double, lat As Double ) As Double

    Const earthRadius = 3958.75 // miles
    Const pi = 3.1415926535897932384626433832795
    Const deg2rad = 0.01745329252
    Dim x1, y1, x2, y2, a, b, c, distance As Double

    Const myLat = 34.04527778
    Const myLong = -118.54694444

    x1 = myLong * deg2rad
    y1 = mylat * deg2rad
    x2 = long * deg2rad
    y2 = lat * deg2rad

    a = sin(( y2-y1 ) / 2.0 )^2
    b = sin(( x2-x1 ) / 2.0 )^2
    c = sqrt( a + cos( y2 ) * cos( y1 ) * b )

    distance = 2 * asin( c ) * earthRadius

    Return distance

    ====================

    Thank you all!

    Here's the final code:

    ======
    Function radians(d as double) As Double
    const pi = 3.1415926535897932384626433832795
    return pi/180*d
    End Function

    ======
    Function decimalDeg(dir As string, d As integer, m As integer, s As Double) As Double
    Dim tSeconds As Integer
    Dim fractionalPart As Double
    Dim result As Double

    tSeconds = ( m * 60 ) + s
    fractionalPart = tSeconds / 3600 result = d + fractionalPart
    if dir = "S" Then result = result * -1
    if dir = "W" Then result = result * -1
    Return result
    End Function

    ======
    Function GCdistance(Lat as double, Long as double) As double
    Const myLat = 34.04527778
    Const myLong = -118.54694444
    Dim DistPart1, Miles, km As Double

    DistPart1=Sin(Radians(myLat))*Sin(Radians(Lat))+_ Cos(Radians(myLat))*Cos(Radians(Lat))*_ Cos(Abs(Radians(Long-myLong)))

    Miles= Atan(Sqrt((1-DistPart1^2)/DistPart1^2)) * 3959.871
    km= Atan(Sqrt((1-DistPart1^2)/DistPart1^2))* 6372.795
    Return miles 'or km
    End Function

  2. Richard H

    13 Apr 2019 University of East Anglia, Nor...

    If you assume Earth is spherical, then you first need to find the solid angle (as measured from the Earth's centre) between the two points on its surface. So:

    1. Convert the lat-long for the two points to 3D cartesian coordinates, so you have V1.x, V1.y, V1.z and V2.x, V2.y, V2.z, where V1 and V2 indicate vectors or points and .x, .y and .z are their direction cosines,
    2. The Dot Product of two vectors gives the solid angle between them, i.e. V1 dot V2 = acos(angle),
    3. Then get the great circle distance from the angle you've just calculated and Earth's radius.

    If you want to take into account the slight flattening of the Earth (about 1 in 300), it gets more involved.

    Regards - Richard.

  3. Robert W

    13 Apr 2019 Western Canada

    This is the formula that I use for great circle distance.

    Function GCdistance(LatitudeLocation1 as double, LongitudeLocation1 as double, LatitudeLocation2 as double, LongitudeLocation2 as double) as double
      'Latitude and Longitude are decimal values
      DistPart1=Sin(Radians(LatitudeLocation1))*Sin(Radians(LatitudeLocation2))+_
        Cos(Radians(LatitudeLocation1))*Cos(Radians(LatitudeLocation2))*_
        Cos(Abs(Radians(LongitudeLocation2-LongitudeLocation1)))
    
      Miles= Atan(Sqrt((1-DistPart1^2)/DistPart1^2)) * 3959.871
      km= Atan(Sqrt((1-DistPart1^2)/DistPart1^2))* 6372.795
      return miles 'or km
    End Function
    
    Function Radians(d as double) as double
      const pi = 3.1415926535897932384626433832795
      return pi/180*d
    End Function 
  4. Edited 8 months ago

    Thanks.

    Robert:

    Would you be kind enough to also show a function for converting from h.m.s to decimal?

    I am using
    ===============
    decDegree( d As integer, m As integer, s As Integer ) as double

    Dim tSeconds As Integer
    Dim fractionalPart As Double
    Dim result As Double

    tSeconds = ( m * 60 ) + s

    fractionalPart = tSeconds / 3600
    result = d + fractionalPart

    Return result
    ================

  5. Robert W

    13 Apr 2019 Western Canada
    Edited 8 months ago

    Since I don't know the exact string format of your input, I can't help with that without more info. This assumes that you've already extracted the degrees minutes and seconds from the input string, and converted them to numeric values:

    Function DegMinSec2Decimal(d As Double, m As Double, s As Double, dir As String) as Double
      'Convert degrees (or hours)/minutes/seconds to decimal degrees (or hours)
      'West or South values are treated as negative. 
      'Parameters d,m,s are degrees, minutes and seconds respectively.
      'Parameter dir should contain one of N,E,S,W. Any other value will give a zero result
      static dSign() As Integer = Array(0,1,1,-1,-1)
      return (d+m/60+s/3600)*dSign(instr("NESW",dir))
    End Function
  6. Ed S

    13 Apr 2019 Answer

    Thank you all!

    Here's the final code:

    ======
    Function radians(d as double) As Double
    const pi = 3.1415926535897932384626433832795
    return pi/180*d
    End Function

    ======
    Function decimalDeg(dir As string, d As integer, m As integer, s As Double) As Double
    Dim tSeconds As Integer
    Dim fractionalPart As Double
    Dim result As Double

    tSeconds = ( m * 60 ) + s
    fractionalPart = tSeconds / 3600 result = d + fractionalPart
    if dir = "S" Then result = result * -1
    if dir = "W" Then result = result * -1
    Return result
    End Function

    ======
    Function GCdistance(Lat as double, Long as double) As double
    Const myLat = 34.04527778
    Const myLong = -118.54694444
    Dim DistPart1, Miles, km As Double

    DistPart1=Sin(Radians(myLat))*Sin(Radians(Lat))+_ Cos(Radians(myLat))*Cos(Radians(Lat))*_ Cos(Abs(Radians(Long-myLong)))

    Miles= Atan(Sqrt((1-DistPart1^2)/DistPart1^2)) * 3959.871
    km= Atan(Sqrt((1-DistPart1^2)/DistPart1^2))* 6372.795
    Return miles 'or km
    End Function

or Sign Up to reply!