Rendering Shapely Geometries in Matplotlib

Update: I just received a reminder that numpy.asarray doesn't copy data.

Here's an example of using Shapely to connect OGR data sources to Matplotlib:

import ogr
import pylab
from numpy import asarray
from shapely.wkb import loads

source = ogr.Open("/tmp/world_borders.shp")
borders = source.GetLayerByName("world_borders")

fig = pylab.figure(1, figsize=(4,2), dpi=300)

while 1:
    feature = borders.GetNextFeature()
    if not feature:
        break

    # Make a Shapely geometry from exported WKB
    geom = loads(feature.GetGeometryRef().ExportToWkb())

    # Wrap the geometry in a Numpy array, slice out lat/long, and plot
    a = asarray(geom)
    pylab.plot(a[:,0], a[:,1])

pylab.show()

The result:

I hope you'll agree that this is considerably simpler than the code I used at the 2005 Open Source Geospatial workshop. An even more direct solution for ogr.py fans would be to provide the Numpy array interface directly from OGR geometries.