I just did this the other day. I doubt I set this up the best way but... it's something. This gives you the distance in miles, if you read the macro it should be easy to figure out how to make it in kilometers. Also note that if you have a ton of points, the TempDistance table is goign to be huge. I haven't really figured out hte best way to use this, I usually just do something like select a to_city and then search for distance <=50 or something. First you have to do something like this:
TempDistance:
NOCONCATENATE LOAD DISTINCT
city,
city.latitude,
city.longitude
RESIDENT CityGPS;
JOIN (TempDistance) LOAD DISTINCT
city as 'to_city',
city.latitude as 'to_latitude',
city.longitude as 'to_longitude'
RESIDENT CityGPS;
That join gives you all combinations that you need the distance from. In my case, I wanted the distance from each city, to each city. Then do something like this:
Distance:
LOAD DISTINCT
city,
to_city,
round(haversine(city.latitude,city.longitude,to_latitude,to_longitude),'.1') as 'distance'
RESIDENT TempDistance;
Add the two macros below:
Function haversine(lat1,long1,lat2,long2)
pi=4*Atn(1)
r = 6371 'radius of the earth or something
lat1 = lat1 * pi / 180
lat2 = lat2 * pi / 180
long1 = long1 * pi / 180
long2 = long2 * pi/180
dLat = (lat2 - lat1) '* (pi / 180)
dLon = (long2 - long1)' * (pi / 180)
a = (sin(dLat / 2) * sin(dLat / 2)) + cos(lat1) * cos(lat2)* (sin(dLon / 2) * sin(dLon / 2))
c = 2 * atan2(sqr(a),sqr(1-a))
d = r * c
d = d * .621371192 '.621371192 converts km to miles
haversine = d
end Function
Function atan2(ys,xs)
If xs <> 0 Then
theta = Atn(ys / xs)
If xs < 0 Then
theta = theta + pi
End If
Else
If ys < 0 Then
theta = 3 * pi / 2 '90
Else
theta = pi / 2 '270
End If
End If
atan2 = theta
End Function
