Shapefile and matplotlib: choosing a polygonal collection of shapefile coordinates

Question

I am trying to build flooded polygons of countries on a world map using matplotlib in python.

I have a shapefile with the coordinates of the borders of the country of each country. Now I want to convert these coordinates (for each country) into a polygon with matplotlib. Without using basemap. Unfortunately, the details overlap or overlap. Is there any employment, perhaps the use of distance from point to point .. or their reordering?

Answer 1

Ha! I found out how .. I completely neglected the info sf.shapes [i] .parts! 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:

Answer 2

Here is another piece of code that I used to create polygonal shapefiles. It uses GDAL / OGR to read the shapefile and graphics for the correct creation of 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() 

Answer 3

 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 that this suggests that polygons, MultiPolygons can be processed in a similar way with

 map(PolygonPatch, MultiPolygon(record['geometry'])) 

Answer 4

As for @hannesk's answer, you should add the following imports: from numpy import array and import matplotlib and replace the line cm = get_cmap('Dark2') with cm = matplotlib.cm.get_cmap('Dark2')

(I'm not so famous as to add a comment to a post I noticed.)