zip code formula/algorithm question [Solved]

[ice_531]

Hey all.

Just curious about this as I've never come across the need.

Is there any way to calculate distance between different zip codes, mathematically or vb functions?

Example.

A listbox containing 6 different zipcodes within a state.

Is there any decent way to sort the listbox of zip codes based on distance to a zipcode in a textbox the user inputs.

Nearest zip code at top of listbox and furthest at bottom.

any code ideas? other than using a 3rd party website functions.

[LaVolpe]

Not sure it is possible without a database with geo-location information regarding each zip code. Maybe some exist. I'd imagine a simplified method would be to know (i.e., database), the lat/long of the center of each zipcode and calculate the distance between the two. There are formulas for calculating distances based on lat/long.

[MartinLiss]

You can get a database of zip codes with latitude and longitude for $30 here. (Never used it)

[MartinLiss]

And this may be free including non-VB code.

[LaVolpe]

You can get a database of zip codes with latitude and longitude for $30 here. (Never used it)

If I recall, I think you can actually get one from the US Post Office website (usps.com), but I'm no longer positive about that.

Edited: The latest one Martinliss posted sounds promising; though it doesn't indicate how current the info is.

[ice_531]

Thanks for the replies.

I had done some research and found the haversine forumla (used to calculate distance between 2 points on a sphere.) which I read is the most accurate to use for this. You can apparently get a free database from US Census Bureau of zip codes to use.

Didn't know if there was an easier way or not. Looks like I'll be going that route, havent seen anything but C# or C code implementing it though.

[MartinLiss]

If you post in one of the C forums someone can probably covert the code for you.

[ice_531]

Ok. After finding different code samples for this problem, none of which are in VB6 unfortunately....ive put together some code in a test app. I'm just not getting accurate results. I mostly went off javascript code.. which presented an issue being there is an atan2 function whereas VB does not have this. also js calls a float function.

Here is my code from test app, and also some sample results.

Public Function Atan2(ByVal y As Double, ByVal x As Double) _
    As Double
  Dim theta As Double
 
  If (Abs(x) < 0.0000001) Then
    If (Abs(y) < 0.0000001) Then
      theta = 0#
    ElseIf (y > 0#) Then
      theta = 1.5707963267949
    Else
      theta = -1.5707963267949
    End If
  Else
    theta = Atn(y / x)
  
    If (x < 0) Then
      If (y >= 0#) Then
        theta = 3.14159265358979 + theta
      Else
        theta = theta - 3.14159265358979
      End If
    End If
  End If
    
  Atan2 = theta
End Function
 
Private Sub Calculate()
Dim R As Double  '' Earth's mean radius in km
R = 6371
Const Pi = 3.14159265358979 '...32384626433832795028841971693993751...
Dim dLat As Double
Dim dLon As Double
Dim ylatFirst As Double
Dim ylatSecond As Double
Dim ylonFirst As Double
Dim ylonSecond As Double
Dim ylatFirst_floated As Double
Dim ylatSecond_floated As Double
Dim lat1 As Double
Dim lat2 As Double
Dim lon1 As Double
Dim lon2 As Double
Dim i As Double
Dim z As Double
Dim a As Double
Dim c As Double
Dim d As Double
 
ylatFirst = Text1.Text
ylonFirst = Text2.Text
 
ylatSecond = Text3.Text
ylonSecond = Text4.Text
 
 
dLat = (ylatSecond - ylatFirst) * Pi / 180 '' convert this to radians
dLon = (ylonSecond - ylonFirst) * Pi / 180 '' convert this to radians
ylatFirst_floated = (ylatFirst) * Pi / 180  '' float ( )
ylatSecond_floated = (ylatSecond) * Pi / 180 '' float ( )
a = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) + Math.Cos(ylatFirst_floated) * Math.Cos(lat2) * Math.Sin(dLon / 2) * Math.Sin(dLon / 2)
c = 2 * Atan2(Math.Sqr(a), Math.Sqr(1 - a))
d = R * c  '' d equals to the distance between the two points in kilometers.
 
Text5.Text = d
End Sub

Here are actual results from 2 different examples.

36.066545  - lat
-80.30733  - long



36.116854  - lat
-79.88291  - long


Actual Distance is 23.9 miles.

My app comes out to 26.6 miles though.

2.7 mile difference

--------------------------------------------

36.066545  - lat
-80.30733  - long



28.945269  - lat
-95.935707 - long


Actual Distance is 1,033.5 miles.

My app comes out to 1089 miles.

55.5 mile difference

I'm sure my calculations are off somewhere obviously. However, almost every code snippet i saw that used the haversine formula, was different in some way. even the ones written in Jscript. Here is the link to the Javascript code I based the above from: http://www.developer.nokia.com/Commu..._in_JavaScript

Any ideas as to what I'm doing wrong?

[LaVolpe]

I'm finding this returns pretty close results

Haversine from http://www.xtremevbtalk.com/showthread.php?t=290297
atan2 from http://forums.devx.com/archive/index.php/t-69063.html

Private Function atan2(y As Double, x As Double) As Double

If x > 0 Then
    atan2 = Atn(y / x)
ElseIf x < 0 Then
    atan2 = Sgn(y) * (Pi - Atn(Abs(y / x)))
ElseIf y = 0 Then
    atan2 = 0
Else
    atan2 = Sgn(y) * Pi / 2
End If

End Function

Private Function Haversine(ByVal Lat1 As Double, _
                            ByVal Long1 As Double, _
                            ByVal Lat2 As Double, _
                            ByVal Long2 As Double) As Double
    
    Const R As Integer = 6371 * 0.621371192  'earth radius in miles
        
    Dim DeltaLat As Double, DeltaLong As Double
    Dim a As Double, c As Double
    Dim Pi As Double
    
    On Error GoTo ErrorExit
    
    Pi = 4 * Atn(1)
    
    'convert Lat1, Lat2, Long1, Long2 from degrees into radians
    Lat1 = Lat1 * Pi / 180
    Lat2 = Lat2 * Pi / 180
    Long1 = Long1 * Pi / 180
    Long2 = Long2 * Pi / 180
    
    'calculate change in Latitude and Longitude
    DeltaLat = Lat2 - Lat1
    DeltaLong = Long2 - Long1
    
    a = ((Sin(DeltaLat / 2)) ^ 2) + (Cos(Lat1) * Cos(Lat2) * ((Sin(DeltaLong / 2)) ^ 2))
    c = 2 * atan2(Sqr(a), Sqr(1 - a))
    Haversine = R * c
ErrorExit:
End Function

Edited: Change Const R to Double vs Integer or instead convert Haversine from kilometers to miles after it returns

Took me a bit to see it, but your error in your haversine formula was simply this. In the calculation of variable "a":
change Math.Cos(Lat2) to Math.Cos(ylatSecond_floated)

[Milk]

I wish I'd seen this thread earlier...

This is my version of Vincenty's formulae.
It's more complicated than Haversine but is supposed to address all distances better, Haversine errors with antipodal points... kind of irrelevant here.

Private Function GetDistance(ByVal Latt1 As Single, _
                             ByVal Long1 As Single, _
                             ByVal Latt2 As Single, _
                             ByVal Long2 As Single, _
                             Optional ByVal InputRadians As Boolean) As Single
'Const EARTHRADIUS As Single = 3438 'NMi (Average)
Const EARTHRADIUS As Single = 3959 'mi (Average)
'Const EARTHRADIUS As Single = 6371 'KM (Average)
'Const EARTHRADIUS As Single = 6371000 'M (Average)
Const PI As Single = 3.14159265
Dim n As Single, d As Single, cxd As Single, sxd As Single
Dim sy1 As Single, sy2 As Single, cy1 As Single, cy2 As Single

  If InputRadians = False Then
    d = PI / 180!
    Latt1 = Latt1 * d
    Long1 = Long1 * d
    Latt2 = Latt2 * d
    Long2 = Long2 * d
  End If

  d = Long2 - Long1
  cxd = Cos(d)
  sxd = Sin(d)
  cy1 = Cos(Latt1)
  cy2 = Cos(Latt2)
  sy1 = Sin(Latt1)
  sy2 = Sin(Latt2)
  
  d = (cy2 * sxd)
  n = (cy1 * sy2) - (sy1 * cy2 * cxd)
  n = Sqr(d * d + n * n)
  d = (sy1 * sy2) + (cy1 * cy2 * cxd)
  
  If d <> 0 Then
    GetDistance = Atn(n / d) * EARTHRADIUS
  ElseIf n <> 0 Then
    GetDistance = PI * 0.5! * EARTHRADIUS
  Else
    GetDistance = 0
  End If
  
End Function

I can see I used atan rather than atan2 although I think that was intentional because I've an earlier version using atan2 which I've commented out (I can't remember my thinking it was ages ago)

[ice_531]

Oddly enough, Milk's code gets me the most accurate number it seems.


With Milks code and the first set of test results i posted...i get 23.95163 miles distance. using the following coordinates.

36.066545  - lat
-80.30733  - long

36.116854  - lat
-79.88291  - long

Which is pretty close to the 23.92 i'm getting from different websites. For my purpose I only need accuracy to the first decimal place anyways.


I do like the atan2 function. i think the only difference between the atan and atan2 is (Y,X) instead of (x/y). but anywho.

Thanks for the help. Once i get code finished to pull all the zip codes and their coordinates into lists and sort that out ill post full code in case anyone would like to see it.

[jmsrickland]

If the Java script code worked for you then why don't you just use the Java script code in your VB project

[jmsrickland]

Thanks for the help. Once i get code finished to pull all the zip codes and their coordinates into lists and sort that out ill post full code in case anyone would like to see it.
__________________
:::`DISCLAIMER`:::
Do NOT take anything i have posted to be truthful in any way, shape or form.
Thank You!

So, then, you might not

[ice_531]

Making new topic at this point. Milk's code worked out perfectly for this threads problem.

[ice_531]

Thanks for the help. Once i get code finished to pull all the zip codes and their coordinates into lists and sort that out ill post full code in case anyone would like to see it.
__________________

I don't think I understand either of your posts. Care to clarify?

[jmsrickland]

I don't think I understand either of your posts. Care to clarify?

Well, my first post was in reference to your statement in post #8. I figured if you had the Java script code and it worked for you then why not use it in your VB code so you wouldn't have to write a bunch of VB code that may or may not work as well as the Java.

My second post was just a joke. Your disclamer Do NOT take anything i have posted to be truthful in any way, shape or form + your comment ...ill post full code...
prompted me to say So, then, you might not

[ice_531]

I don't think I understand either of your posts. Care to clarify?

Well, my first post was in reference to your statement in post #8. I figured if you had the Java script code and it worked for you then why not use it in your VB code so you wouldn't have to write a bunch of VB code that may or may not work as well as the Java.

My second post was just a joke. Your disclamer Do NOT take anything i have posted to be truthful in any way, shape or form + your comment ...ill post full code...
prompted me to say So, then, you might not

Ah. gotcha. lol


Well the javascript code I hadnt tested, just assumed it worked and tested after I re-wrote it in VB. I didn't really want to use other code whether it be javascript,c# etc. trying to get back into VB again, so need to stay within it for the most part

It's been a lonnnng time since I even read that signature. Think I put that back in 2003 or something. Been here a long time, but inactive for many years after switching majors away from programming. The irony is with my current job, I do a lot of side projects by request of different department supervisors and most of them involve me having to code in some aspect. Whether it be simple batch scripting or database help...and more recently actual applications in Vb.

It's fun and scary at the same time, getting back into things. I have a lot of source code from apps i did years ago, and sadly 50% of it is above my head at this point. So will be drifting about these forums more actively now.

[DrUnicode]

Try Function distVincenty which is supposed to be quite accurate:

Download site appears to be offline:
http://www.windreader.com/geodesy/Do.../GeodesyVB.zip

Code for Vincenty distance is also here:
http://lost-species.livejournal.com/38453.htm

Option Explicit

Private Const PI = 3.14159265358979
Private Const EPSILON As Double = 1e-12

Public Sub test()
    MsgBox distVincenty(52.874, 4.389, 45.001, 15.716)
End Sub

Private Function toRad(ByVal degrees As Double) As Double
    toRad = degrees * (PI / 180)
End Function

Private Function Atan2(ByVal X As Double, ByVal Y As Double) As Double
 ' code nicked from:
 ' http://en.wikibooks.org/wiki/Programming:Visual_Basic_Classic/Simple_Arithmetic#Trigonometrical_Functions
 ' If you re-use this watch out the x and y have been reversed.
    If Y > 0 Then
        If X >= Y Then
            Atan2 = Atn(Y / X)
        ElseIf X <= -Y Then
            Atan2 = Atn(Y / X) + PI
        Else
        Atan2 = PI / 2 - Atn(X / Y)
    End If
        Else
            If X >= -Y Then
            Atan2 = Atn(Y / X)
        ElseIf X <= Y Then
            Atan2 = Atn(Y / X) - PI
        Else
            Atan2 = -Atn(X / Y) - PI / 2
        End If
    End If
End Function

Public Function distVincenty(ByVal lat1 As Double, ByVal lon1 As Double, ByVal lat2 As Double, ByVal lon2 As Double) As Double
'=================================================================================
' Calculate geodesic distance (in m) between two points specified by latitude/longitude (in numeric degrees)
' using Vincenty inverse formula for ellipsoids
'=================================================================================
' Code has been ported by lost_species from www.aliencoffee.co.uk to VBA from javascript published at:
' http://www.movable-type.co.uk/scripts/latlong-vincenty.html
' * from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the
' *       Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975
' *       http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
'=================================================================================
' Copyright lost_species 2008 LGPL http://www.fsf.org/licensing/licenses/lgpl.html
'=================================================================================

    Dim low_a As Double
    Dim low_b As Double
    Dim f As Double
    Dim L As Double
    Dim U1 As Double
    Dim U2 As Double
    Dim sinU1 As Double
    Dim sinU2 As Double
    Dim cosU1 As Double
    Dim cosU2 As Double
    Dim lambda As Double
    Dim lambdaP As Double
    Dim iterLimit As Integer
    Dim sinLambda As Double
    Dim cosLambda As Double
    Dim sinSigma As Double
    Dim cosSigma As Double
    Dim sigma As Double
    
    Dim sinAlpha As Double
    Dim cosSqAlpha As Double
    Dim cos2SigmaM As Double
    Dim C As Double
    Dim uSq As Double
    Dim upper_A As Double
    Dim upper_B As Double
    Dim deltaSigma As Double
    Dim s As Double
    
    low_a = 6378137
    low_b = 6356752.3142
    f = 1 / 298.257223563 'WGS-84 ellipsiod
    
    L = toRad(lon2 - lon1)
    U1 = Atn((1 - f) * Tan(toRad(lat1)))
    U2 = Atn((1 - f) * Tan(toRad(lat2)))
    sinU1 = Sin(U1)
    cosU1 = Cos(U1)
    sinU2 = Sin(U2)
    cosU2 = Cos(U2)
    
    lambda = L
    lambdaP = 2 * PI
    iterLimit = 20
    
     While (Abs(lambda - lambdaP) > EPSILON) And (iterLimit > 0)
        iterLimit = iterLimit - 1
        
        sinLambda = Sin(lambda)
        cosLambda = Cos(lambda)
        sinSigma = Sqr(((cosU2 * sinLambda) ^ 2) + ((cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ 2))
        If sinSigma = 0 Then
            distVincenty = 0  'co-incident points
            Exit Function
        End If
        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
        sigma = Atan2(cosSigma, sinSigma)
        sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
        cosSqAlpha = 1 - sinAlpha * sinAlpha

        If cosSqAlpha = 0 Then 'check for a divide by zero
            cos2SigmaM = 0 '2 points on the equator
        Else
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
        End If
        
        C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
        lambdaP = lambda
        lambda = L + (1 - C) * f * sinAlpha * _
            (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * (cos2SigmaM ^ 2))))
     Wend
    
    If iterLimit < 1 Then
        MsgBox "iteration limit has been reached, something didn't work."
        Exit Function
    End If

    uSq = cosSqAlpha * (low_a ^ 2 - low_b ^ 2) / (low_b ^ 2)
    upper_A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    upper_B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
    deltaSigma = upper_B * sinSigma * (cos2SigmaM + upper_B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM ^ 2) _
        - upper_B / 6 * cos2SigmaM * (-3 + 4 * sinSigma ^ 2) * (-3 + 4 * cos2SigmaM ^ 2)))
    
    s = low_b * upper_A * (sigma - deltaSigma)
    distVincenty = s

End Function

[Milk]

Ahh, I posted a special case Vincenty which treats both axes as equal which simplifies loads (third formulae down). DrUnicode has posted the classic version which accounts for the earth being fatter around the equator.

[ice_531]

Haven't found any downsides to your formula on all that I've tested thus far Milk. It's much simpler than the original above. For my understanding at least

also before, do not require super accuracy, as the first decimal place is fine.

Does make me wonder which is faster though, out of all the different variations i've seen of the haversine and Vincenty formula's.

[Jonney]

Read http://williams.best.vwh.net/avform.htm#LL
attach http://www.windreader.com/geodesy/Do.../GeodesyVB.zip