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()

4 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!

  3. Hi Doglus,
    Dont know if you are still dealing with that problem. The files is a zip it conatins all the files needed to build the shapes on the map you must just specify the folder with all thoses files in them. Make sure you dont doubleup on names as windows sometimes does.
    Regards Kaspar

    ps. I really enjoyed this tutorial it has helped me a lot in my work, I am investigating wheter basemap or cartography is better if the author has any comments that would be helpful.

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>