GeoAlchemy Example

This example neighborhood.p.y shows how to get GeoAlchemy <http://geoalchemy.org> to recognize an existing PostGIS table, read the table’s contents, and then call a couple of PostGIS functions on that data, all from python.

Setting up PostGIS

First we assume a working PostGIS installation (see <http://postgis.refractions.net/docs/ch02.html#PGInstall>). The code expects to connect to a database called postgres, with a user called postgres (without a password).

Using another database name, user and/or passoword?
That's cool...just change the following line in neighborhood.py:

database = "postgresql://postgres@localhost:5432/postgres"

Example data

The example is built on a table called ‘neighborhood’, which is the result of a ‘shp2pgsql’ on a file we’ll download from Zillow. First choose a neighborhood .shp file from here: http://www.zillow.com/howto/api/neighborhood-boundaries.htm. Then load that into your database.

wget http://www.zillow.com/static/shp/ZillowNeighborhoods-MA.zip
unzip ZillowNeighborhoods-MA.zip
shp2pgsql ZillowNeighborhoods-MA.shp neighborhood > neighborhood.sql
psql -U postgres -f neighborhood.sql

Run neighborhood.p.y

Fairly simple here, just type python neighborhood.py, and you’re running the example, and should see something like the following (i.e., since I have more than MA data loaded, your result will differ a bit...thinking Sixteen Acres will be the result of MA only data)

python neighborhood.py
....
DB Columns ['neighborhood.gid', 'neighborhood.state', 'neighborhood.county', 'neighborhood.city', 'neighborhood.name', 'neighborhood.regionid', 'neighborhood.the_geom']
GeoColumn Type ST_MultiPolygon
Salem centriod is at: 44.9177349523 , -123.014934057 and this big 0.0116839510301

Summary

Please take neighborhood.p.y apart, to see how it works. I’m hoping that this helps (more than hurts) any novice python / SqlAlchemy / GeoAlchemy folks out there. Creating this has helped me start to come up to speed with these technologies. In the process of learning, I never found a stand-alone example to read existing PostGIS tables, so I’m hoping this fills that void.

Good luck... http://frankpurcell.com

neighborhood.py

#
# Example GeoAlchemy query from existing PostGIS table 
# Frank Purcell
# June 1, 2010
#
import sys
from sqlalchemy     import *
from sqlalchemy.orm import *
from geoalchemy     import *
from sqlalchemy.sql import func

#
# DB connection
# NOTE: set echo=True to see all that *Alchemy is doing to talk to the database
#
database = "postgresql://postgres@localhost:5432/postgres"
engine   = create_engine(database, echo=False)
metadata = MetaData()


#
# define relational object, mapping to db table
# tell sqlalchemy about the special type of the PostGIS geom column
# other columns will be mapped automatically from the table's metadata
#
source = Table('neighborhood', 
               metadata,
               GeometryExtensionColumn('the_geom', Geometry(2)),
               autoload=True,
               autoload_with=engine
         )
# NOTE that the Did not recognize type 'geometry' of column 'the_geom' warning is happening here...not sure why
# sys.exit(1)


#
# open a bound session on the database
#
Session = sessionmaker(bind=engine)
session = Session()


#
# query for the largest polygon by doing an order_by on the st_area(geom), 
# then choosing the first row from the result set
#
# NOTE: The order_by content is simply a string that gets appended to geoalchemy's sql query being 
#       sent to the db.  The contents of the order_by (as shown here), is any valid sql applicable to the  
#       given (postgis) db.  The negative to this approach would be lack of portability (eg: sqllite).
#
q = session.query(source)
#row = q..first()
#row = q.order_by(desc("name")).first()
#row = q.order_by("st_area(the_geom)").first()         # smallest polygon
row = q.order_by(desc("st_area(the_geom)")).first()   # largest polygon


#
# print out information about our table
#
print "DB Columns", source.columns
print "GeoColumn Type", session.scalar(functions.geometry_type(row.the_geom))


#
# Call an ST_ function on our geom column, and print out information about the first row in the table  
#
center = DBSpatialElement(row.the_geom.centroid())
area   = session.scalar(DBSpatialElement(row.the_geom.area()))
lon    = session.scalar(center.y)
lat    = session.scalar(center.x)
print row.name, "centriod is at:", lon, ",", lat, "and this big",area