Matplotlib Basemap tutorial 10: Shapefiles Unleached, continued

Following the comments on this year-old post I’ve had a look at pyshp which seems a little more maintained (the last line in the changelog is 1 year old…), and it is indeed a quite nice piece of code. I thus rewrote “Tutorial 7″ to use this module :

The goal:

tutorial10

The data:

http://www.gadm.org/ saved inside a new “borders/” folder !

The idea:
Opening a GADM shapefile, get region names, and plot filled regions with random color !

The power:

We first instanciate the Basemap (see below), then import and get the shapes (from the shp file) and records (from the dbf files):

import shapefile

r = shapefile.Reader(r"borders\bel_adm3")
shapes = r.shapes()
records = r.records()

then, we get the points we need to create the matplotlib polygons:

for record, shape in zip(records,shapes):
    lons,lats = zip(*shape.points)
    data = np.array(m(lons, lats)).T

within this loop, we will also check that the shape is not in multiple parts, and if yes, segment the points in different ensembles:

    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:])

and finally, we plot the polygons and give them a random facecolor :

    lines = LineCollection(segs,antialiaseds=(1,))
    lines.set_facecolors(cm.jet(np.random.rand(1)))
    lines.set_edgecolors('k')
    lines.set_linewidth(0.1)
    ax.add_collection(lines)

Voilààààà :-)

The full code is after the break:

#
# BaseMap example by geophysique.be
# tutorial 10

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

### PARAMETERS FOR MATPLOTLIB :
import matplotlib as mpl
mpl.rcParams['font.size'] = 10.
mpl.rcParams['font.family'] = 'Comic Sans MS'
mpl.rcParams['axes.labelsize'] = 8.
mpl.rcParams['xtick.labelsize'] = 6.
mpl.rcParams['ytick.labelsize'] = 6.

fig = plt.figure(figsize=(11.7,8.3))
#Custom adjust of the subplots
plt.subplots_adjust(left=0.05,right=0.95,top=0.90,bottom=0.05,wspace=0.15,hspace=0.05)
ax = plt.subplot(111)
#Let's create a basemap of Europe
x1 = -5.0
x2 = 15.
y1 = 45.
y2 = 54.

m = Basemap(resolution='i',projection='merc', llcrnrlat=y1,urcrnrlat=y2,llcrnrlon=x1,urcrnrlon=x2,lat_ts=(x1+x2)/2)
m.drawcountries(linewidth=0.5)
m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(y1,y2,2.),labels=[1,0,0,0],color='black',dashes=[1,0],labelstyle='+/-',linewidth=0.2) # draw parallels
m.drawmeridians(np.arange(x1,x2,2.),labels=[0,0,0,1],color='black',dashes=[1,0],labelstyle='+/-',linewidth=0.2) # draw meridians

from matplotlib.collections import LineCollection
from matplotlib import cm
import shapefile

r = shapefile.Reader(r"borders\bel_adm3")
shapes = r.shapes()
records = r.records()

for record, shape in zip(records,shapes):
    lons,lats = zip(*shape.points)
    data = np.array(m(lons, lats)).T

    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:])

    lines = LineCollection(segs,antialiaseds=(1,))
    lines.set_facecolors(cm.jet(np.random.rand(1)))
    lines.set_edgecolors('k')
    lines.set_linewidth(0.1)
    ax.add_collection(lines)

plt.savefig('tutorial10.png',dpi=300)
plt.show()

3 thoughts on “Matplotlib Basemap tutorial 10: Shapefiles Unleached, continued

  1. Pingback: Matplotlib Basemap tutorial 07: Shapefiles unleached | Géophysique.be

  2. Thanks for your post. I’m new to base map and was looking to this type of work. I’m have a little issue about which files you are using after downloading. Here are the files I downloaded. You will see that I have one set that I unzipped and also the zip file is in the dir or folder too.

    Your code shows the folder ‘borders’ and the filename, ‘ bel_adm3′ . No files that I downloaded is named like that. So I assume I should be importing the ‘gadm2.shp’ file is that correct? r=shapefile.Reader(r”gadm2.shp”) is that correct?

    import shapefile
    r = shapefile.Reader(r”borders\bel_adm3″)
    shapes = r.shapes()
    records = r.records()

    MY Downloaded ls:
    02/21/2014 03:16 AM 349,551,441 gadm_v2_shp.zip
    01/20/2012 10:53 PM 515,043,378 gadm2.dbf
    01/20/2012 10:52 PM 145 gadm2.prj
    01/20/2012 10:53 PM 2,037,236 gadm2.sbn
    01/20/2012 10:53 PM 39,668 gadm2.sbx
    01/20/2012 10:53 PM 993,682,208 gadm2.shp
    01/20/2012 10:53 PM 1,746,004 gadm2.shx
    12/20/2011 03:07 PM 212,369 read_me.pdf
    Thanks!

Leave a Reply

Your email address will not be published. Required fields are marked *


*

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>