[Tutor] Geolocating objects

[Terry Carroll]

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)