Python GEOS Module

Update (2008-01-24): See http://pypi.python.org/pypi/Shapely/1.0.

Update (2007-06-12): http://sgillies.net/blog/478/shapely-geometries-for-python/.

Update (2007-04-10): http://sgillies.net/blog/429/geometries-for-python-update/.

Update (2006-06-01): I'm no longer using SWIG to make an interface for GEOS. Instead, I'm working on a better GEOS-based geometry module for Python.

GEOS is the Geometry Engine Open Source, a library for programming with 2D geometries and spatial predicates. Just before the 2.0 release, I committed a SWIG interface file and Python distutils setup script for generating a Python 'geos' module on Linux. Hobu has put together a geos (from GEOS 2.1.1) module for Python 2.4 and is distributing it from http://hobu.biz/software/PyGEOS. This is the module that he and I will be using in our OSG05 workshop.

It's easy to get started with the Python geos module. There's no specific documentation for it, but you can readily adapt the GEOS docs to use with Python; the API is almost exactly the same. Below are examples of using the geos GeometryFactory, WKTReader, and WKTWriter classes, and also of the Union and intersection methods of the geos Geometry class.

For kicks, I'll render the results using PCL's cartography module. First, an example of two triangle polygons rendered separately using the default polygon symbolizer and a simple text symbolizer:

from cartography import data, mapping, rendering, spatial
import geos

# Setup the renderer and styling
# ----------------------------------------------------------------------
FONTSET = '/home/sean/projects/ms_45/mapserver/tests/fonts.txt'
mapper = rendering.MapRenderer('MAPSERVER', fontset=FONTSET)

# Text and Polygon symbolizers
text_symb = mapping.TextSymbolizer('FNAME', {'font-family': 'Vera',
                                             'font-size': 10.0},
                                   placement='AUTO')
poly_symb = mapping.PolygonSymbolizer()

# Style and two rules
pe_style = mapping.ExtrinsicStyle('pe_style')
pe_style.rules.append(mapping.Rule('catch-all', text_symb))
pe_style.rules.append(mapping.Rule('catch-all', poly_symb))

# Create two intersecting triangles using PCL
# ----------------------------------------------------------------------
# We're treating everything as WGS 1984
srs = spatial.SpatialReference(epsg=4326)

t1wkt = 'POLYGON ((-1.0 50.5, -0.5 51.2, 0.3 50.9, -1 50.5))'
triangle1 = spatial.Polygon(wkt=t1wkt)
atts1 = ('A Polygon',)

t2wkt = 'POLYGON((-0.7 50.3, 0.1 51.0, 0.6 50.1, -0.7 50.3))'
triangle2 = spatial.Polygon(wkt=t2wkt)
atts2 = ('B Polygon',)

# Render the two shapes separately
# -----------------------------------------------------------------------
store = data.MemoryFeatureStore([('FNAME', 'string')],
            zip([triangle1, triangle2], [atts1, atts2]), srs=srs)
layer = mapping.Layer(store, 'POLYGON')
view = mapping.View(srs, spatial.Point(-1.2, 51.2), spatial.Point(1, 50))
size = [300, 300]
image = mapper.render([(layer, pe_style)], view,
                      'image/png', size, transparent=1)
f = open('separate.png', 'wb')
f.write(image.getdata())
f.close()

the resulting image, separate.png:

Use the read() method of WKTWriter to convert the PCL polygons to geos Geometries and then call the intersection method of one geometry with the other as the single argument. A new geometry is returned, and can be rendered by stuffing it into a PCL data store and wrapping with a layer:

# Next, render intersection of the polygons
# -----------------------------------------------------------------------
geom_factory = geos.GeometryFactory()
wkt_reader = geos.WKTReader(geom_factory)
wkt_writer = geos.WKTWriter()

geom1 = wkt_reader.read(triangle1.toWKT())
geom2 = wkt_reader.read(triangle2.toWKT())
g = geom1.intersection(geom2)

# make new PCL geometry
inter = spatial.Polygon(wkt=wkt_writer.write(g))

inter_store = data.MemoryFeatureStore([('FNAME', 'string')],
                                      [(inter, ('Intersection',))], srs=srs)
inter_layer = mapping.Layer(inter_store, 'POLYGON')
image = mapper.render([(inter_layer, pe_style)], view,
                      'image/png', size, transparent=1)
f = open('intersection.png', 'wb')
f.write(image.getdata())
f.close()

the resulting image, intersection.png:

The union of the two triangle polygons is just as easily obtained:

# Now, render the union of the polygons
# -----------------------------------------------------------------------
geom_factory = geos.GeometryFactory()
wkt_reader = geos.WKTReader(geom_factory)
wkt_writer = geos.WKTWriter()

geom1 = wkt_reader.read(triangle1.toWKT())
geom2 = wkt_reader.read(triangle2.toWKT())
g = geom1.Union(geom2)  # <-- watch out for capitalized union

# make new PCL geometry
union = spatial.Polygon(wkt=wkt_writer.write(g))

union_store = data.MemoryFeatureStore([('FNAME', 'string')],
                                      [(union, ('Union',))], srs=srs)
union_layer = mapping.Layer(union_store, 'POLYGON')
image = mapper.render([(union_layer, pe_style)], view,
                      'image/png', size, transparent=1)
f = open('union.png', 'wb')
f.write(image.getdata())
f.close()

resulting image, union.png: