[Tutor] Geolocating objects

Terry Carroll carroll at tjc.com
Tue Feb 13 02:29:26 CET 2007

```anil maran was asking about determining distances between two US zip
codes.

I pointed him to a couple resources, but just for the fun of it, tried my
hand at it.

As far as I know, this actually works:

from math import radians, sin, asin, cos, sqrt, atan2

class ZCTA:
"""
class for each of the Zip Code Tabulation Areas (U.S. Census)
"""

def __init__(self, line):
"""
format from http://www.census.gov/geo/www/gazetteer/places2k.html:

* Columns 1-2: United States Postal Service State Abbreviation
* Columns 3-66: Name (e.g. 35004 5-Digit ZCTA - there are no
post office names)
* Columns 67-75: Total Population (2000)
* Columns 76-84: Total Housing Units (2000)
* Columns 85-98: Land Area (square meters) - Created for
statistical purposes only.
* Columns 99-112: Water Area (square meters) - Created for
statistical purposes only.
* Columns 113-124: Land Area (square miles) - Created for
statistical purposes only.
* Columns 125-136: Water Area (square miles) - Created for
statistical purposes only.
* Columns 137-146: Latitude (decimal degrees) First character
is blank or "-" denoting North or South latitude respectively
* Columns 147-157: Longitude (decimal degrees) First character
is blank or "-" denoting East or West longitude respectively
"""

import struct

format = "2s64s9s9s14s14s12s12s10s11sc"
assert len(line) == 158
assert struct.calcsize(format) == 158
# 157 chars as defined above, + one character for newline

_tuple = struct.unpack(format, line)
self.state = _tuple
self.ZCTA_name = _tuple
self.zipcode = self.ZCTA_name[0:5]
self.pop = _tuple
self.housing_units = _tuple
self.land_area_metric = _tuple
self.water_area_metric = _tuple
self.land_area = _tuple
self.water_area = _tuple
self.latitude = float(_tuple)
self.longitude = float(_tuple)

zfile = open("zcta5.txt","r")
ZCTAs = {}
for line in zfile:
z = ZCTA(line)
ZCTAs[z.zipcode] = z

R = 3956.08833 # circumference of the earth (miles)

while True:
zipcode1 = raw_input("first zip code: ")
if zipcode1 == "":
break
z1 = ZCTAs[zipcode1]
print "starting: %s, lat: %s, long: %s" % \
(z1.zipcode, z1.latitude, z1.longitude)
zipcode2 = raw_input("second zip code: ")
if zipcode2 == "":
break
z2 = ZCTAs[zipcode2]
print "starting: %s, lat: %s, long: %s" % \
(z2.zipcode, z2.latitude, z2.longitude)

# Calculate the distance:
# http://www.movable-type.co.uk/scripts/GIS-FAQ-5.1.html
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin((dlat/2.0)**2) + cos(lat1) * cos(lat2) * sin((dlon/2.0)**2)
# c = 2 * asin(min(1,sqrt(a)))    # alternate formula
c = 2 * atan2(sqrt(a), sqrt(1-a)) # distance in radians
d_miles = R * c # distance in miles
d_km = d_miles * 1.609344

print "distance from %s to %s is %d miles (%d km)\n" % \
(zipcode1, zipcode2, d_miles, d_km)

```