2

I'm trying to calculate distance in kilometers between two geographical coordinates using the haversine formula.

Code:

Dim dbl_dLat As Double
Dim dbl_dLon As Double
Dim dbl_a As Double

dbl_P = WorksheetFunction.Pi / 180
dbl_dLat = dbl_P * (dbl_Latitude2 - dbl_Latitude1)
dbl_dLon = dbl_P * (dbl_Longitude2 - dbl_Longitude1)

dbl_a = Sin(dbl_dLat / 2) * Sin(dbl_dLat / 2) + Cos(dbl_Latitude1 * dbl_P) * Cos(dbl_Latitude2 * dbl_P) * Sin(dbl_dLon / 2) * Sin(dbl_dLon / 2)

dbl_Distance_KM = 6371 * 2 * WorksheetFunction.Atan2(Sqr(dbl_a), Sqr(1 - dbl_a))

I'm testing with these coordinates:

dbl_Longitude1 = 55.629178
dbl_Longitude2 = 29.846686
dbl_Latitude1 = 37.659466
dbl_Latitude2 = 30.24441

And the code returns 20015.09, which is obviously wrong. It should be 642 km according to Yandex maps.

Where am I wrong? Are the longitude and latitude in wrong format?

cxw
  • 16,685
  • 2
  • 45
  • 81
User24112017
  • 65
  • 1
  • 4
  • Why not `Lon` rather than `dbl_dLon`? Do those pseudo-Hungarian prefixes really help with code reliability? It definitely doesn't help in the readability department. – John Coleman Dec 28 '17 at 13:14
  • Welcome to the site! Check out the [tour](https://stackoverflow.com/tour) for more about how the site runs. You can [edit your question](https://stackoverflow.com/posts/48008116/edit) to include more information. In this case, the coordinates given are [2512 km apart](http://www.gcmap.com/mapui?P=37.659466%2C55.629178-30.24441%2C29.846686&MS=wls&DU=km). Would you please double-check your input coordinates? – cxw Dec 28 '17 at 13:16
  • 1
    @ cxw sorry, dbl_Longitude2 should read as 59.846686 (then the points are 642 km away). Anyway, as you pointed out, the trouble was because of the wrong argument placement for Atan2 function. – User24112017 Dec 28 '17 at 14:02
  • @User24112017 Thanks for the update! By the way, please don't put a space after the at sign --- `@cxw` should notify me, but `@ cxw` did not. – cxw Dec 28 '17 at 14:51

2 Answers2

2

As far as I can tell, the issue is that the order of arguments to atan2() varies by language. The following works* for me:

Option Explicit

Public Sub Distance()
    Dim dbl_Longitude1 As Double, dbl_Longitude2 As Double, dbl_Latitude1 As Double, dbl_Latitude2 As Double

    dbl_Longitude1 = 55.629178
    dbl_Longitude2 = 29.846686
    dbl_Latitude1 = 37.659466
    dbl_Latitude2 = 30.24441

    Dim dbl_dLat As Double
    Dim dbl_dLon As Double
    Dim dbl_a As Double
    Dim dbl_P As Double

    dbl_P = WorksheetFunction.Pi / 180
    dbl_dLat = dbl_P * (dbl_Latitude2 - dbl_Latitude1)      'to radians
    dbl_dLon = dbl_P * (dbl_Longitude2 - dbl_Longitude1)    'to radians

    dbl_a = Sin(dbl_dLat / 2) * Sin(dbl_dLat / 2) + _
            Cos(dbl_Latitude1 * dbl_P) * Cos(dbl_Latitude2 * dbl_P) * Sin(dbl_dLon / 2) * Sin(dbl_dLon / 2)

    Dim c As Double
    Dim dbl_Distance_KM As Double
    c = 2 * WorksheetFunction.Atan2(Sqr(1 - dbl_a), Sqr(dbl_a))  ' *** swapped arguments to Atan2
    dbl_Distance_KM = 6371 * c

    Debug.Print dbl_Distance_KM
End Sub

*Output: 2507.26205401321, although gcmap.com says the answer is 2512 km. This might be a precision issue --- I think it's close enough to count as working. (Edit it might also be that gcmap uses local earth radii rather than the mean radius; I am not sure.)

Explanation

I found this description of the haversine formula for great-circle distance, which is what you are implementing. The JavaScript implementation on that page gives this computation for c:

var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));

In JavaScript, atan2() takes parameters y, x. However, in Excel VBA, WorksheetFunction.Atan2 takes parameters x, y. Your original code passed Sqr(dbl_a) as the first parameter, as it would be in JavaScript. However, Sqr(dbl_a) needs to be the second parameter in Excel VBA.

A comment on naming

Building on @JohnColeman's point, there are lots of ways to name variables. In this case, I would recommend using the prefixes for unit rather than for type: e.g., deg_Latitude1, RadPerDeg = Pi/180, and rad_dLat = RadPerDeg * (deg_Latitude2 - deg_Latitude1). I personally think that helps avoid unit-conversion mishaps.

cxw
  • 16,685
  • 2
  • 45
  • 81
1

My VBA code that returns the answer in feet; However 'd' is the answer in kilometers.

Imports System.Math
Module Haversine
Public Function GlobalAddressDistance(sLat1 As String, sLon1 As String, sLat2 As String, sLon2 As String) As String
    Const R As Integer = 6371
    Const cMetersToFeet As Single = 3.2808399
    Const cKiloMetersToMeters As Integer = 1000
    Dim a As Double = 0, c As Double = 0, d As Double = 0

    'Convert strings to numberic double values
    Dim dLat1 As Double = Val(sLat1)
    Dim dLat2 As Double = Val(sLat2)
    Dim dLatDiff As Double = DegreesToRadians(CDbl(sLat2) - CDbl(sLat1))
    Dim dLonDiff As Double = DegreesToRadians(CDbl(sLon2) - CDbl(sLon1))

    a = Pow(Sin(dLatDiff / 2), 2) + Cos(DegreesToRadians(dLat1)) * Cos(DegreesToRadians(dLat2)) * Pow(Sin(dLonDiff / 2), 2)
    c = 2 * Atan2(Sqrt(a), Sqrt(1 - a))
    d = R * c

    'Convert kilometers to feet
    Return Format((d * cKiloMetersToMeters * cMetersToFeet), "0.##").ToString
End Function

Private Function DegreesToRadians(ByVal dDegrees As Double) As Double
    Return (dDegrees * PI) / 180
End Function

End Module

MarkHL
  • 11
  • 2