The matplotlib basemap toolkit is a library for plotting 2D data on maps in Python. The Basemap constructor takes as input points specifying the latitude and longitude of the location of interest, along with an optional argument, resolution, that specifies the desired quality of the rendered map. This Basemap instance exposes the readshapefile method, which contains the necessary details to render 2D maps with various types of boundaries. The Basemap constructor can accept many additional arguments, but they will not be covered here. Instead, I’ll demonstrate how to render a map of Chicago, Illinois that represents variation in population density by zip code using matplotlib’s default colormaps. The end result will be rendered as follows:
In this post, we reference the shapefiles available for download
here. The archive will
contain a number of files which are then consumed by the readshapefile method of
the basemap instance. The Chicago population data by zipcode can be downloaded here.
To begin, we first demonstrate how to create an un-shaded 2D map of Chicago with zipcode boundaries rendered. Chicago’s coordinates are (41.8781N, 87.6298W), and we can use this knowledge to construct a 4-point window that surrounds the geographic location of interest:
"""
Render a 2D map of Chicago with zip code boundaries.
"""
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import os
import os.path
us_shape_file_dir = "cb_2017_us_zcta510_500k"
os.chdir(us_shape_file_dir)
# Chicago coordinates.
lowerlon = -88.2
upperlon = -87.2
lowerlat = 41.62
upperlat = 42.05
m = Basemap(
llcrnrlon=lowerlon,
llcrnrlat=lowerlat,
urcrnrlon=upperlon,
urcrnrlat=upperlat,
resolution='c',
projection='lcc',
lat_0=lowerlat,
lat_1=upperlat,
lon_0=lowerlon,
lon_1=upperlon
)
shp_info = m.readshapefile(os.path.basename(us_shape_file_dir), 'state')
plt.gca().axis("off")
plt.show()
Running the above produces the following:
The Basemap constructor resolution argument can be set to one of the following:
- “c” (crude)
- “l” (low)
- “i” (intermediate)
- “h” (high)
- “f” (full)
- None
To represent population density by zip code graphically is slightly more involved. Instead of describing step-by-step how to go about it, I’ll present the code, then will comment on a few sections after the rendering is generated:
from mpl_toolkits.basemap import Basemap
from matplotlib.collections import LineCollection
import matplotlib as mpl
from matplotlib.colors import rgb2hex
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import os
# Read in population data.
pop_path = "ChicagoPopulation.csv"
DF = read_csv("Sheet1")
colormap = plt.cm.Purples
us_shape_file_dir = "cb_2017_us_zcta510_500k"
os.chdir(us_shape_file_dir)
# Chicago coordinates.
lowerlon = -88.2
upperlon = -87.2
lowerlat = 41.62
upperlat = 42.05
m = Basemap(
llcrnrlon=lowerlon,
llcrnrlat=lowerlat,
urcrnrlon=upperlon,
urcrnrlat=upperlat,
projection="lcc",
resolution="c",
lat_0=lowerlat,
lat_1=upperlat,
lon_0=lowerlon,
lon_1=upperlon
)
shp_info = m.readshapefile(
os.path.basename(us_shape_file_dir), 'states', drawbounds=True
)
# Convert integer ZIP5 field to character dtype.
DF['ZIP5'] = DF['ZIP5'].astype(str)
# Read population density info into popdens dict. Take square root of
# actual density for better color mapping.
popdens = {
str(ii): np.sqrt(jj) for (ii, jj) in zip(DF.ZIP5.values, DF.POPULATION.values)
}
# Choose a color for each state based on population density. Range
# vmin-vmax has arbitrarily been set to 0-6.
ziplist = []
colors = {}
vmin = 0.
vmax = 6.
# Filter m.states_info to only Chicago zipcodes.
zip_info = m.states_info
popdiv = (max(popdens.values()) / (vmax-vmin))
popdensscl = {ii: (jj / popdiv) for (ii ,jj) in popdens.items()}
for d in zip_info:
iterzip = d["ZCTA5CE10"]
if iterzip in popdensscl.keys():
iterpop = popdensscl.get(iterzip, 0)
colors[iterzip] = colormap(iterpop / vmax)[:3]
ziplist.append(iterzip)
for nshape,seg in enumerate(m.states):
i, j = zip(*seg)
if ziplist[nshape] in popdensscl.keys():
color = rgb2hex(colors[ziplist[nshape]])
edgecolor = "#000000"
plt.fill(i,j,color,edgecolor=edgecolor)
# (Optional) include colorbar.
sm = plt.cm.ScalarMappable(
cmap=colormap, norm=mpl.colors.Normalize(vmin=vmin, vmax=vmax)
)
mm = plt.cm.ScalarMappable(cmap=colormap)
mm.set_array([vmin, vmax])
plt.colorbar(mm, ticks=np.arange(vmin, vmax + 1, 1), orientation="vertical")
plt.title("Chicago Population by ZIP5")
plt.gca().axis("off")
plt.show()
Running the code above will produce the following: