Volunteered Geographic Information

Thanks to Sean Gillies for commenting on my last post, he put me onto a couple of Python packages that he’s been involved in creating that allow you to do some really excellent geospatial things. The shapely package is a great implementation of a lot of spatial analyses that you can do on projected (i.e. flattened) datasets, including topological operations and a full set of object types. The descartes package allows better integration of matplotlib with spatial data, particularly in terms of not having to use the “fill” plotting function repeatedly, but creating a more efficient set of “patches” which can then be added to the figure plot. The overal impression I got from descartes is that it wasn’t spectacularly different from the method detailed in my previous post, but it gives you more control and stability over the map plotting process; whereas using raw matplotlib you are inclined to hope that the map outputs correctly (it all seems a bit up to chance), using descartes you have a more robust and easily manipulable output.

In order to test this I rewrote my previous thematic map script to: firstly convert the shapefile geometries into shapely polygons, and secondly to pass those shapely polygons to descartes and draw a map plot using descartes-matplotlib. The only slightly odd piece of functionality that I found was that you can’t pass the shapely polygon object a list of shapely points in order to create the polygon, rather you have to pass a list of x,y tuples – much less satisfying!

Nonetheless, the changes were easy to implement, and with the previous script as given basically include:

from shapely.geometry import Polygon

points = []
for i in range(0,number of points in shapefile):
 tempx = float(x coord of point in shapefile polygon)
 tempy = float(y coord of point in shapefile polygon)

 points.append((tempx,tempy))
polygon = Polygon(points)

The above method creates a simple polygon without holes, shapely can accomodate this is need be though. Having created the shapely polygons, all that remains is to create a patch.

from descartes import PolygonPatch

patch = PolygonPatch(polygon, plus colour and line considerations)

Then you simply add the patch to the matplotlib figure you have already created so:

from matplotlib import pyplot

fig = pyplot.figure(1, figsize = [10,10], dpi = 300)   #create 10x10 figure
ax = fig.addsubplot(111)    #Add the map frame (single plot)

# here you create all the polygons and patches

ax.addpatch(patch)   # simply add the patch to the subplot
# set plot vars
ax.set_xlim(get xmin and xmax values from data)
ax.set_ylim(get ymin and ymax values from data)
ax.set_aspect(1)

pyplot.show()

Using these basics I was able to create a basic OAC map using Welsh OAs as an example: