Working around a mapping bug in cartopy

Posted in hacks and kludges on Thursday, September 01 2016

I recently discovered a bug in how I was generating choropleth maps with cartopy. Basically the add_geometries() function for adding shapes to matplotlib maps was not mapping the correct face colours to the correct geometry. As an example consider the two maps below (of where university students live in Edmonton).

University students in Edmonton, wrong

This first map was done using the following:

ax.add_geometries(df['geometry'], ..., facecolour=df['colour'])

How convenient! A one line way of making a choropleth map! Except the neighbourhoods are all coloured incorrectly. Note it is not that the colour field for each row is wrong, it's that the colour on the map is not the colour in the dataframe for that neighbourhood. The correct map should look like this:

University students in Edmonton, correct

So what is going on? As far as I can tell (without actually going into any sourcecode or otherwised doing anything helpful) is that cartopy is flattening the list of geometries then applying colours. Unfortunately some neighbourhoods are single Polygons and others are MultiPolygons, so flattening this all down into a list of polygons is quite destructive. It no longer lines up with the list of colours (and I assume there is some default behaviour as to why it isn't crying foul when it runs off the end of the list of colours).

If I iterate through the dataframe row by row and add geometries that way, everything works out fine. This is how I produced the correct map.

for index, row in df.iterrows():
    # the add_geometries function expects an iterator, so better give it one
    except TypeError:
        geom = (row['geometry'],)
        geom = row['geometry']

    ax.add_geometries(geom, ccrs.PlateCarree(), facecolor=row['colour'],
                      edgecolor='#555555', alpha=0.8)