13

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.Generated By Code Actual File

  • Not familiar with the shapefile reader, however I can tell that you are just appending all the points in the file to one big list without separating each shape into its component parts. You need a big list of shapes to which you append each shapes points – user681 Jan 25 '15 at 23:08
  • Right. Have to find a way to separate the shapes. But that is what I am unable to do at the moment. – statBeginner Jan 25 '15 at 23:24
  • @DanPatterson Can you specify how to plot multiple shapes in the same figure after I manage to separate the shapes? If i use plt.plot(listx,listy) for every shape, it keeps generating a new figure every time, instead of using the same figure. – statBeginner Jan 26 '15 at 0:12
10

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()
  • oh.. wonder how i missed that. i do get the shapes printed in different colors though. Will have to fix that :) – statBeginner Jan 26 '15 at 5:46
  • how to get or isolate the different shapes? – FaCoffee Jun 21 '16 at 10:29
15

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.

  • 6
    I know this might be redundant information, but for those not yet familiar with the subject it would have been useful to say that import shapefile refers to the pyshp package: pypi.python.org/pypi/pyshp – FaCoffee Feb 15 '17 at 10:42
  • This is not ok when you have a bunch of islands, as these points will be connected by lines to points on the mainland, creating something similar to what the OP posted. – FaCoffee Sep 28 '17 at 11:21
  • 1
    @FaCoffee, you are right. My answer gis.stackexchange.com/a/309780/126618 should address this. – Gus Jan 24 at 17:07
7

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

enter image description here

  • That's really helpful, but can you do this in a for loop if you have multiple polygons to plot? – FaCoffee May 15 '17 at 13:10
  • Yes without problem – gene May 15 '17 at 15:39
  • I noticed that the descartes solution does not work if you try and plot two different shapefiles on two adjacent subplots using fig, ax = plt.subplots(1,2,figsize=(15, 8)) and then ax[0].add_patch(PolygonPatch(poly_geo, fc='#d3d3d3', ec='#000000', alpha=0, zorder=5)) and ax[1].add_patch(PolygonPatch(poly_geo, fc='#d3d3d3', ec='#000000', alpha=0, zorder=5)). The result is an empty image. Any idea? – FaCoffee Jan 7 '18 at 18:02
2

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

2

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.

0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.