I recently started working in crime analysis and started to need to use latitudes and longitudes and calculate the distances between them. For this calculation I used the Haversine formula which calculates the geodesic distance between two points using the Spherical Law of Cosines.

Warning: this formula assumes the Earth is a perfect sphere,which is not the case and results in inaccuracies compared to more rigorous methods. For my purposes it works fine.

First, since the function I’m going to create uses radians, and latitudes and longitudes are basically a measurement of degrees, we need to create some helper functions to translate back and forth. Here they are:

#function to convert degrees to radians deg2rad <- function(deg) return(deg*pi/180)

Now we’re going to create the actual function. One of the things that messed with my head a little bit with this whole process is the fact that every R function that deals with GIS uses longitude first, then latitude. As a former navigator, I did the opposite for years. But the reason R does it this way, is because when you’re plotting latitude and longitude, longitude is going to end up in the X plane, which is always first in plotting functions.

Also take note, that this function is going to return the distance in statute miles. This is easy to change, just change the Earth mean radius part to your measurement of choice, and the function will return the distance in that metric. Here’s the function:

gcd.slc <- function(long1, lat1, long2, lat2) { # Convert degrees to radians long1<-deg2rad (long1) lat1<-deg2rad(lat1) long2<-deg2rad(long2) lat2<-deg2rad(lat2) R <- 3959 # Earth mean radius [miles] d <- acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(long2-long1)) * R return(d) # Distance in miles }

And that’s it. To call it:

gcd.slc(-112.0738, 33.448266, -111.9322, 33.477) [1] 8.400443

Notice the WGS84 coordinates, that’s the system that Google uses in it’s maps, so it’s convenient to flip back and forth. For this answer you should get 8.4 statute miles.