I'm trying to plot filled polygons of countries on the world map with matplotlib in python.
I've got a shapefile with country boundary coordinates of every country. Now, I want to convert these coordinates (for each country) into a polygon with matplotlib. Without using Basemap. Unfortunately, the parts are crossing or overlapping. Is there a workarund, maybe using the distance from point to point.. or reordering them ?
Ha! I found out, how.. I completely neglected, the sf.shapes[i].parts information! Then it comes down to:
# -- import --
import shapefile
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
# -- input --
sf = shapefile.Reader("./shapefiles/world_countries_boundary_file_world_2002")
recs = sf.records()
shapes = sf.shapes()
Nshp = len(shapes)
cns = []
for nshp in xrange(Nshp):
cns.append(recs[nshp][1])
cns = array(cns)
cm = get_cmap('Dark2')
cccol = cm(1.*arange(Nshp)/Nshp)
# -- plot --
fig = plt.figure()
ax = fig.add_subplot(111)
for nshp in xrange(Nshp):
ptchs = []
pts = array(shapes[nshp].points)
prt = shapes[nshp].parts
par = list(prt) + [pts.shape[0]]
for pij in xrange(len(prt)):
ptchs.append(Polygon(pts[par[pij]:par[pij+1]]))
ax.add_collection(PatchCollection(ptchs,facecolor=cccol[nshp,:],edgecolor='k', linewidths=.1))
ax.set_xlim(-180,+180)
ax.set_ylim(-90,90)
fig.savefig('test.png')
Then it will look like this:
Here is another piece of code I used to plot polygon shapefiles. It uses GDAL/OGR to read shapefile and plots correctly donut shape polygons:
from osgeo import ogr
import numpy as np
import matplotlib.path as mpath
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
# Extract first layer of features from shapefile using OGR
ds = ogr.Open('world_countries_boundary_file_world_2002.shp')
nlay = ds.GetLayerCount()
lyr = ds.GetLayer(0)
# Get extent and calculate buffer size
ext = lyr.GetExtent()
xoff = (ext[1]-ext[0])/50
yoff = (ext[3]-ext[2])/50
# Prepare figure
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlim(ext[0]-xoff,ext[1]+xoff)
ax.set_ylim(ext[2]-yoff,ext[3]+yoff)
paths = []
lyr.ResetReading()
# Read all features in layer and store as paths
for feat in lyr:
geom = feat.geometry()
codes = []
all_x = []
all_y = []
for i in range(geom.GetGeometryCount()):
# Read ring geometry and create path
r = geom.GetGeometryRef(i)
x = [r.GetX(j) for j in range(r.GetPointCount())]
y = [r.GetY(j) for j in range(r.GetPointCount())]
# skip boundary between individual rings
codes += [mpath.Path.MOVETO] + \
(len(x)-1)*[mpath.Path.LINETO]
all_x += x
all_y += y
path = mpath.Path(np.column_stack((all_x,all_y)), codes)
paths.append(path)
# Add paths as patches to axes
for path in paths:
patch = mpatches.PathPatch(path, \
facecolor='blue', edgecolor='black')
ax.add_patch(patch)
ax.set_aspect(1.0)
plt.show()
from fiona import collection
import matplotlib.pyplot as plt
from descartes import PolygonPatch
from matplotlib.collections import PatchCollection
from itertools import imap
from matplotlib.cm import get_cmap
cm = get_cmap('Dark2')
figure, axes = plt.subplots(1)
source_path = "./shapefiles/world_countries_boundary_file_world_2002"
with collection(source_path, 'r') as source:
patches = imap(PolygonPatch, (record['geometry'] for record in source)
axes.add_collection( PatchCollection ( patches, cmap=cm, linewidths=0.1 ) )
axes.set_xlim(-180,+180)
axes.set_ylim(-90,90)
plt.show()
Note this assumes polygons, MultiPolygons can be handles in a similar manner with
map(PolygonPatch, MultiPolygon(record['geometry']))
Regarding to @hannesk's answer, you should add the following imports: from numpy import array and import matplotlib and replace the line cm = get_cmap('Dark2') by cm = matplotlib.cm.get_cmap('Dark2')
(I'm not so famous to add a comment to the noticed post.)
