How to Make a Map with Python and Matplotlib

Intro

One thing I like to do, it’s to put things on a map. I feel like this is a way to reattach them to life. I see many maps of COVID19 everywhere on the internet, but few show the regional distribution… At least here! So, this post will show you how to do your own in Python. Don’t worry it won’t be so hard.

Required Libraries

The classics

Libraries that deals with maps exist already in python. But, for this demo, we will directly use Matplotlib. This will be a good opportunity to see some functions that aren’t used much. We will need an array from Numpy one, we won’t really use it for something else this time.

Shapefile

I wish we could have used no libraries at all to get the geometries of the regions, but we will need pyshp. Usually, I like to work with GEOJSON, which is simply JSON file with standardized layout. But my province administration only makes shapefile available. They are strange collection of separated binary files, so dealing with the data would require a bit too much of reading for this quick demo… A good thing, pyshp make working with those data files easy. You can install it with your regular pip command.

#/usr/bin/python3
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import shapefile as sf

The code

The Data

I use a dictionary at the start of my code for this demo. But you could think of an automatized solution if your state/province is large (it will be, it’s small here!). If you need it, you might want to look at urllib and Beautifulsoup to get your data. I use this website for Quebec data.

casQc={
        "Bas-Saint-Laurent":4,
        "Saguenay - Lac-Saint-Jean":5,
        "Capitale-Nationale":79,
        "Mauricie":32,
        "Centre-du-Québec":32,
        "Estrie":130,
        "Montréal":439,
        "Outaouais":12,
        "Abitibi-Témiscamingue":5,
        "Côte-Nord":1,
        "Nord-du-Québec":1,
        "Gaspésie - îles-de-la-Madeleine":7,
        "Chaudière-Appalaches":14,
        "Laval":60,
        "Lanaudière":55,
        "Laurentides":44,
        "Montérégie":125
    }

Problems With French Characters (éôîè)

The government is French here. Which is not really a problem… But the government is slow here. Which is a problem. The file distributed by the province is encoded with ‘latin1’, which is now offset by the much more beautiful UTF-8 and UTF-16. We will have to deal with it.

def fixlatin(inStr):
    toFix = [
                ('\x93','ô'),
                ('\x82','é'),
                ('×','î'),
                ('\x8a','è')
            ]
    for c,o in toFix:
        inStr = inStr.replace(c,o)
    return inStr

Prepare the data and display

We will use Patches from the Matplotlib Library. There are few shapes, but we will only need Polygons. Matplotlib also has a tool to directly use data to colours the patches with a colormap. This tool is PatchCollection.

So, to produce our regional map, we will get the geographic borders of the regions, generate the polygons, add them to a collection, then display the whole set!

# Load the shapes file
shapefile = sf.Reader('region_admin_poly',encoding='latin1')

# Prepare the patch using the shapes containeds in region_admin_poly
# and also prepare the data
patches = []
rdata = []

# I use this gradient of color
# I think "Gradients" are better for this kind of data
colmap = cm.get_cmap('YlOrRd')
for (rec,shape) in zip(shapefile.records(),shapefile.shapes()):
    label = fixlatin(rec[7])
    patches.append(Polygon(shape.points))
    rdata.append(casQc[label])

# Prepare the display
# We will need direct control of the axes
# to get them you could also use (f,a=plt.subplots())
# but I find this more explicit
figure = plt.figure()
axe = plt.axes()
figure.patch.set_visible(False)
collection = PatchCollection(patches,edgecolors='k')

collection.set_cmap(colmap)
collection.set_array(np.array(colors))

axe.axis('off')

# This is the "drawing" method
axe.add_collection(collection)
axe.set_aspect('equal') 

plt.colorbar(collection)
plt.grid(False)
plt.xlim([-81,-56])
plt.ylim([44,63])
plt.show()

Gabriel St-Pierre-LemieuxMarch 25, 2020