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[0]
self.ZCTA_name = _tuple[1]
self.zipcode = self.ZCTA_name[0:5]
self.pop = _tuple[2]
self.housing_units = _tuple[3]
self.land_area_metric = _tuple[4]
self.water_area_metric = _tuple[5]
self.land_area = _tuple[6]
self.water_area = _tuple[7]
self.latitude = float(_tuple[8])
self.longitude = float(_tuple[9])
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
lon1 = radians(z1.longitude)
lon2 = radians(z2.longitude)
lat1 = radians(z1.latitude)
lat2 = radians(z2.latitude)
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)
_______________________________________________
Tutor maillist - Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor
