This notebook is the result of random experiments working with Matplotlib/Basemap and Shapefiles. It reuses chuncks of code from several blog posts, stack overflow answers, and other sources.
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm
from matplotlib.collections import LineCollection
from matplotlib.patches import Polygon
from matplotlib.patches import Polygon
from mpl_toolkits.basemap import Basemap
%matplotlib inline
import shapefile
def draw_portugal(islands = True):
"""
This functions draws and returns a map of Portugal, either just of the mainland or including the Azores and Madeira islands.
"""
fig = plt.figure(figsize=(15.7,12.3))
projection='merc'
llcrnrlat=-80
urcrnrlat=90
llcrnrlon=-180
urcrnrlon=180
resolution='i'
if islands:
x1 = -32.
x2 = 5.
y1 = 30.
y2 = 45.
else:
x1 = -12.
x2 = -5.
y1 = 36.
y2 = 44.
m = Basemap(projection=projection, llcrnrlat=y1, urcrnrlat=y2, llcrnrlon=x1,
urcrnrlon=x2, resolution=resolution)
m.drawcoastlines()
m.drawmapboundary()
m.drawcountries()
m.fillcontinents(color = '#C0C0C0')
return m
draw_portugal()
plt.show()
Alright, now that we set up our imports and support function, let's start using the shapefiles. For this part I used an official Shapefile provided by the Portuguese authorities. You can download the files here.
Let's load and inspect the information:
r = shapefile.Reader(r"maps\prt_adm1")
shapes = r.shapes()
records = r.records()
print r.fields
print "\n----------------------------------------------------------------------------------\n"
for record in records:
print record
[('DeletionFlag', 'C', 1, 0), ['ID_0', 'N', 9, 0], ['ISO', 'C', 3, 0], ['NAME_0', 'C', 75, 0], ['ID_1', 'N', 9, 0], ['NAME_1', 'C', 75, 0], ['NL_NAME_1', 'C', 50, 0], ['VARNAME_1', 'C', 150, 0], ['TYPE_1', 'C', 50, 0], ['ENGTYPE_1', 'C', 50, 0]] ---------------------------------------------------------------------------------- [180, 'PRT', 'Portugal', 1, 'Aveiro', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 2, 'Azores', ' ', 'A\xe7ores', 'Regi\xf5es aut\xf4noma', 'Autonomous region'] [180, 'PRT', 'Portugal', 3, 'Beja', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 4, 'Braga', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 5, 'Bragan\xe7a', ' ', 'Braganza', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 6, 'Castelo Branco', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 7, 'Coimbra', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 8, '\xc9vora', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 9, 'Faro', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 10, 'Guarda', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 11, 'Leiria', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 12, 'Lisboa', ' ', 'Lisbon|Lisbona|Lisbonne|Lissabon', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 13, 'Madeira', ' ', ' ', 'Regi\xf5es aut\xf4noma', 'Autonomous region'] [180, 'PRT', 'Portugal', 14, 'Portalegre', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 15, 'Porto', ' ', 'Oporto', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 16, 'Santar\xe9m', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 17, 'Set\xfabal', ' ', ' ', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 18, 'Viana do Castelo', ' ', 'Vianna do Castello', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 19, 'Vila Real', ' ', 'Villa Real', 'Distrito', 'District'] [180, 'PRT', 'Portugal', 20, 'Viseu', ' ', 'Vizeu', 'Distrito', 'District']
Each record, which corresponds to a shape, represents a Portuguese "Distrito". Let's append the population of each district to the record as a number between 0 and 1 (1 being the population of the most inhabited district).
populacao = {"Aveiro": 735790,
"Azores": 246746,
"Beja": 150287,
"Braga": 848185,
"Bragança": 139344,
"Castelo Branco": 225916,
"Coimbra": 429714,
"Évora": 168034,
"Faro": 434023,
"Guarda": 168898,
"Leiria": 470895,
"Lisboa": 2244799,
"Madeira": 267785,
"Portalegre": 118952,
"Porto": 2027191,
"Santarém": 465701,
"Setúbal": 866794,
"Viana do Castelo": 250390,
"Vila Real": 213775,
"Viseu": 391215
}
max_pop = np.max(populacao.values())+0.0
for record in records:
record.append(populacao[record[4].decode('Latin1').encode('UTF-8')]/max_pop)
#print a record to make sure that the information was added
print records[0]
[180, 'PRT', 'Portugal', 1, 'Aveiro', ' ', ' ', 'Distrito', 'District', 0.32777544893774452]
m = draw_portugal(islands=False)
ax = plt.subplot(111)
for record, shape in zip(records,shapes):
#read shape
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
#each shape may have different segments
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
#draws the segments, and sets its properties. A colormap is used to get the gradient effect.
lines = LineCollection(segs,antialiaseds=(1,))
lines.set_facecolors(cm.YlGn(record[-1]))
lines.set_edgecolors('k')
lines.set_linewidth(0.1)
ax.add_collection(lines)
plt.show()
Basemap comes with its own Shapefile reader. In this case, we'll use a shapefile with world countries (download it from Geocommons - direct link).
This method is easier to use when you know what you want to do, but it's harder to explore since it doesn't separate the information from the map itself.
fig = plt.figure(figsize=(15.7,12.3))
ax = plt.subplot(111)
map = Basemap(projection='cyl',llcrnrlat=35,urcrnrlat=72,\
llcrnrlon=-17,urcrnrlon=42,resolution='c')
map.drawmapboundary(fill_color='aqua')
map.readshapefile('maps\world_by_country', 'countries_sf', drawbounds=True)
#the readshapefile method allow you to call the shapefile's shapes and info.
#Both are lists, the first containing a list of tuples (coordinates), and the second containig a dictionary with associated metadata
print len(map.countries_sf) #number of shapes
print map.countries_sf_info[0].keys() #metadata
#plots the shapes as Polygons with a random facecolor (if European), or gray (if not)
for shape, country in zip(map.countries_sf, map.countries_sf_info):
if country['REGION'] == 'Europe':
poly = Polygon(shape, facecolor=cm.YlGn(np.random.rand()))
plt.gca().add_patch(poly)
else:
poly = Polygon(shape, facecolor="gray")
plt.gca().add_patch(poly)
#adds the colormap legend
cmleg = np.zeros((1,20))
for i in range(20):
cmleg[0,i] = (i*5)/100.0
plt.imshow(cmleg, cmap=plt.get_cmap("YlGn"), aspect=1)
plt.colorbar()
plt.show()
982 ['NAME_1', 'NAME', 'AREA', 'REGION', 'RINGNUM', 'ISO_2_CODE', 'ISO_3_CODE', 'GMI_CNTRY', 'SHAPENUM', 'NAME_12', 'POP2005']