Plot shapefile with matplotlib

Question

I am trying to read a shapefile and plot it using matplotlib. Here is the code:

import matplotlib.pyplot as plt
import shapefile   

shpFilePath = "D:\test.shp"  
listx=[]
listy=[]
test = shapefile.Reader(shpFilePath)
for sr in test.shapeRecords():
    for xNew,yNew in sr.shape.points:
        listx.append(xNew)
        listy.append(yNew)
plt.plot(listx,listy)
plt.show()

However, i get lines connecting my polygons. How can I draw the polygons such that they are the way in the shapefile. Here are screenshots of the plot and the shapefile when it is opened with ArcGIS.

Answer 1

For future references, here is the solution I have came to after following the advices above.

import shapefile as shp  # Requires the pyshp package
import matplotlib.pyplot as plt

sf = shp.Reader("test.shp")

plt.figure()
for shape in sf.shapeRecords():
    x = [i[0] for i in shape.shape.points[:]]
    y = [i[1] for i in shape.shape.points[:]]
    plt.plot(x,y)
plt.show()

The resulting figure will be very colorful, but then, you just need to adjust the plot keywords.

Answer 2

I will leave it to you how to collect the shapes but this is the principle

import numpy as np
from matplotlib import pyplot as p  #contains both numpy and pyplot
x1 = [-1,-1,10,10,-1]; y1 = [-1,10,10,-1,-1]
x2 = [21,21,29,29,21]; y2 = [21,29,29,21,21]
shapes = [[x1,y1],[x2,y2]]
for shape in shapes:
  x,y = shape
  p.plot(x,y)
p.show()

Answer 3

You need to use matplotlib paths and patches and there is a Python module dedicated to plot polygons from shapefiles using these functions Descartes.

As Pyshp (shapefile) has the geo_interface (New geo_interface for PyShp) convention, you can use it.

polys  = shapefile.Reader("polygon")
# first polygon
poly = polys.iterShapes().next().__geo_interface__
print poly
{'type': 'Polygon', 'coordinates': (((151116.87238259654, 135890.8706318218), (153492.19971554304, 134793.3055883224), (153934.50204650551, 133892.31935858406), (152623.97662143156, 131811.86024627919), (150903.91200102202, 130894.49244872745), (149347.66305874675, 132991.33312884573), (149151.08424498566, 134383.76639298678), (151116.87238259654, 135890.8706318218)),)}

The result is the GeoJSON representation of the geometry and you can use the solution of How to plot geo-data using matplotlib/python

import matplotlib.pyplot as plt 
from descartes import PolygonPatch
BLUE = '#6699cc'
fig = plt.figure() 
ax = fig.gca() 
ax.add_patch(PolygonPatch(poly, fc=BLUE, ec=BLUE, alpha=0.5, zorder=2 ))
ax.axis('scaled')
plt.show()

Answer 4

In addition to ldocao answer, and responding to FaCoffee question. When you have isolated islands and they are part of the same feature, you can try next:

import shapefile as shp
import matplotlib.pyplot as plt

sf = shp.Reader("test.shp")

plt.figure()
for shape in sf.shapeRecords():
    for i in range(len(shape.shape.parts)):
        i_start = shape.shape.parts[i]
        if i==len(shape.shape.parts)-1:
            i_end = len(shape.shape.points)
        else:
            i_end = shape.shape.parts[i+1]
        x = [i[0] for i in shape.shape.points[i_start:i_end]]
        y = [i[1] for i in shape.shape.points[i_start:i_end]]
        plt.plot(x,y)
plt.show()

This makes it work for me. The propertie "parts" of a shape returns the starting indexes of differents geometries inside a feature.

Answer 5

It can be done using either geopandas or pyshp as discussed in this answer. Geopandas use matplotlib at its backend for plotting.

Answer 6

Still, in one shapefile shape, there may be multiple parts. This will plot each part within one shape, separately.

import matplotlib.pyplot as plt
import shapefile
import numpy as np

this_shapefile = shapefile.Reader(map_file_base) # whichever file
shape = this_shapefile.shape(i) # whichever shape
points = np.array(shape.points)

intervals = list(shape.parts) + [len(shape.points)]

ax = plt.gca()
ax.set_aspect(1)

for (i, j) in zip(intervals[:-1], intervals[1:]):
    ax.plot(*zip(*points[i:j]))

Answer 7

The most simple answer is using geopandas itself.

import geopandas as gpd
df=gpd.read_file("filename.shp")
df.plot()