Plotting GIS shapes

Here's another installment in my series (iterators, revenge of the iterators, features) that considers different Python GIS APIs and environments: plotting geometries with Matplotlib. Matplotlib has a number of functions for plotting and graphing, and they typically take as their first arguments sequences of x and y coordinate values. These sequences are immediately adapted to Numpy arrays. In some sense what we're actually considering here is how well different Python APIs integrate with Numpy.

In every code snippet below we import the Matplotlib pylab module and then obtained a native geometry object geom or feature object in some way (see previous post for details).

ESRI:

x = []
y = []
pnt = geom.Next()
while pnt is not None:
    x.append(pnt.x)
    y.append(pnt.y)
    pnt = geom.Next()
pylab.plot(x, y)

Wrapping the geometry up in a point or vertex iterator would help tidy this up.

FME:

x, y = zip(*feature.getCoordinates())
pylab.plot(x, y)

No actual geometry object in FME, but this isn't bad at all in combination with zip. Remember that zip(*x) unzips the sequence of items x. Within a function call, the * operator explodes a sequence.

OGR:

x, y = zip(*((geom.GetX(i), geom.GetY(i)) for i in range(geom.GetGeometryRef(0).GetPointCount()))
pylab.plot(x, y)

Note that using zip gets this done in one pass over the geometry instead of two as in the example linked through OGR above. I'm telling you: idiommatic Python wins.

QGIS:

x, y = zip(*((p.x, p.y) for p in geom.asPolyline()))
pylab.plot(x, y)

geojson:

x, y = zip(*geom.coordinates)
pylab.plot(x, y)

Easier and easier.

Shapely (1.2a6):

x, y = geom.xy
pylab.plot(x, y)

Easiest yet.

Where possible, I've used a Python idiom to compress the code to two lines for each API example. There are a lot of calls in the OGR and QGIS samples and efficiency and clarity are at odds. Since a list is being built no matter what, there's no harm in doing it with more lines of code for those APIs. With Shapely, I've deliberately made it hard to do it in more than two lines. In fact, it's easier to do it in one:

pylab.plot(*geom.xy)