Calculating distances between zipcodes

Coincidentally I was playing around with calculating distances between zipcodes three Saturdays ago (having previously snarfed the zipcode, latitude & longitude info from the Census Bureaus website).

If you want an exhaustive discussion of the various ways to calculate the distance between two points on the surface of the earth, see the discussion at: http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1 or if you just want something you can feed data into and get answers out of, look at this perl script.

I opted to use the "Law of Cosines for Spherical Trigonometry," since it was the simplest to hack into SQL (on Solid). The caveat is that it doesn't work well for small distances (less than one minute of arc). However, for finding zipcodes within a certain (large) radius, etc. I didn't that this would really be much of an issue.

My zt_zipcodes table has, among other goodies courtesty of the Census Bureau, zipcode, latitude, and longitude columns.

The following query works to calculate the distance between two arbitrary zipcodes:

select
acos(
        ( sin(z1.latitude * 0.017453293)*sin(z2.latitude * 0.017453293) )
    +
        ( cos(z1.latitude * 0.017453293) *
          cos(z2.latitude * 0.017453293) *
          cos((z2.longitude*0.017453293)-(z1.longitude*0.017453293))
        )
    )
* 3956 as distance
from zt_zipcodes z1, zt_zipcodes z2
where z1.zipcode = '84041' and z2.zipcode = '02138'; -- 2,086 miles

or to find all the zipcodes (and their respective distance—in miles truncated to one decimal place) within a fixed range from a fixed zipcode:

select z2.zipcode,
truncate(
acos(
        ( sin(z1.latitude * 0.017453293)*sin(z2.latitude * 0.017453293) )
    +
        ( cos(z1.latitude * 0.017453293) *
          cos(z2.latitude * 0.017453293) *
          cos((z2.longitude*0.017453293)-(z1.longitude*0.017453293))
        )
    ) * 3956,1) as distance
from zt_zipcodes z1, zt_zipcodes z2
where z1.zipcode = '27604' and z2.zipcode <> '27604' 
and acos(
        ( sin(z1.latitude * 0.017453293)*sin(z2.latitude * 0.017453293) )
    +
        ( cos(z1.latitude * 0.017453293) *
          cos(z2.latitude * 0.017453293) *
          cos((z2.longitude*0.017453293)-(z1.longitude*0.017453293))
        )
    ) * 3956 between 50 and 100; -- finds 310 zipcodes

Regarding the "magic numbers" above, 0.017453293 is used to convert decimal degrees to radians (pi/180). 3956 represents the radius of the earth in miles. Obviously if you were to use either of the above more than just once or twice you'd probably want to create a view.

This entry was originally posted on the photo.net (later ArsDigita.com) web/db Q&A forum, now archived in read-only mode for posterity. I've moved my post here to my blog so that I can find it myself without googling... — MC, 3 Jan 2005

— Michael A. Cleverly
Saturday, March 04, 2000
30 comments | Printer friendly version