zachary.com

sixtyPercent: Cochlear Implants, Aviation, Technlology, and Philosophy 2005/01

Python Geo Programming

I found this free database of US ZIPcodes with city name, state, and latitude, longitude information. With a bit more searching, I also found these databases of US 2000 Census infomration by ZIPcode. All this information looked like a good opportunity for some Python geographic programming fun

Unfortunately, the census data is in DBase format. Fortunately, there are a couple of ways to convert this data to something useful to me. I found a Python DBase file reader, and used it plus some simple Python DB API code to populate a SQLite database with the information.

I've also been playing with SQLObject recently, which is a very nice Python Object/Relational Mapping system (ORM). It generally just does its job well, and stays out of the way the rest of the time. This fits well with my programming philosophy of simple tools that do one job. Python is excellent at both creating these simple tools, and integrating them together. So I created this simple class:

from sqlobject import *
import math

#
# The following formulas are adapted from the Aviation Formulary
# http://williams.best.vwh.net/avform.htm
#

nauticalMilePerLat = 60.00721
nauticalMilePerLongitude = 60.10793
rad = math.pi / 180.0
milesPerNauticalMile = 1.15078

def calcDistance(lat1, lon1, lat2, lon2):                      
    """
    Caclulate distance between two lat lons in NM
    """
    yDistance = (lat2 - lat1) * nauticalMilePerLat
    xDistance = (math.cos(lat1 * rad) + math.cos(lat2 * rad)) * (lon2 - lon1) * (nauticalMilePerLongitude / 2)

    distance = math.sqrt( yDistance**2 + xDistance**2 )

    return distance * milesPerNauticalMile

def milesBox( lat, lon, radius ):
    """
    Returns two lat/lon pairs as (lat1, lon2, lat2, lon2) which define a box that exactly surrounds
    a circle of radius of the given amount in miles.
    """
    latRange = radius / ( milesPerNauticalMile * 60.0 )
    lonRange = radius / ( math.cos(lat * rad) * milesPerNauticalMile * 60.0)
    
    return ( lat - latRange, lon - lonRange, lat + latRange, lon + lonRange )

#
# Zipcode persistent information class
#

class Zipcode( SQLObject ):
    """
    A collection of information about a US ZIP code. Database key is the 5 digit zipcode itself,
    as an integer. For some zipcodes, no city name is given.
    """
    _connection = 'sqlite:///tmp/zipdb.db'
    _idName = 'zip'
    
    # 'zip' is the ID column
    city = StringCol( default = 'UNKNOWN' )
    state = StringCol( length = 2 )
    latitude = FloatCol()
    longitude = FloatCol()   
    population = IntCol( default = 0 )
    area = FloatCol( default = 0.0 )            # in square miles

    def distance( self, zip ):
        """
        Return the great circle distance between self and the given zipcode.
        """
        return calcDistance( self.latitude, self.longitude,
                                zip.latitude, zip.longitude )

    def nearby( self, radius = 10.0 ):
        """
        Return a list (generator) of all zipcodes within the radius (in miles) of self --
        more or less -- see milesBox doc for true definition.
        """
        (lat1,lon1,lat2,lon2) = milesBox( self.latitude, self.longitude, radius )
        return Zipcode.select( AND (
                                    AND( Zipcode.q.latitude >= lat1,
                                         Zipcode.q.latitude <= lat2 ),
                                    AND( Zipcode.q.longitude >= lon1,
                                         Zipcode.q.longitude <= lon2 ) ) )

***
highlight file error
***

With the class built, we can start to have some fun:

% python
>>> from zipcode import Zipcode
>>> # Southern-most named town in Wyoming
>>> [z for z in Zipcode.select( Zipcode.q.state == 'WY', orderBy='latitude' )][0].city
'Jelm'

>>> # distance from Dallas to San Francisco -- note precision is only to 4 decimal places
>>> Zipcode.get( 94114 ).distance( Zipcode.get( 75248 ) )
1482.3659983019043

>>> # cities within 20 miles of Palo Alto, CA
>>> import sets
>>> sets.Set([x.city for x in Zipcode.get( 94306 ).nearby( 20.0 )])
Set(['Millbrae', 'Santa Clara', 'Los Gatos', 'Union City', 'San Leandro', 'San Gregorio', 'Newark', 'Pacifica', 'San Francisco', 'Loma Mar', 'San Mateo', 'Hayward', 'San Carlos', 'Half Moon Bay', 'Pleasanton', 'Burlingame', 'Mountain View', 'Stanford', 'Fremont', 'San Lorenzo', 'Los Altos', 'La Honda', 'Menlo Park', 'Saratoga', 'San Jose', 'Portola Valley', 'Milpitas', 'Campbell', 'Atherton', 'Sunnyvale', 'Palo Alto', 'Alviso', 'Pescadero', 'Redwood City', 'San Bruno', 'Cupertino', 'Daly City', 'Boulder Creek', 'Brisbane', 'Belmont', 'Castro Valley', 'South San Francisco', 'Sunol'])

***
highlight file error
***

by David Creemer : 2005/01/12 : Categories python : 0 trackbacks : 5 comments (permalink)