Plotting shapely polygon in matplotlib

Question

When working with geometric objects I like working with the shapely module (Shapely - PyPI). In here it is easy to create polygons and define union, intersection and differences. Like so:

from shapely.geometry import Polygon
x1 = [3, 3, 6, 6, 3]
y1 = [4, 8, 8, 4, 4]
x2 = [1, 8, 8, 1, 1]
y2 = [1, 1, 6, 6, 1]
rectangle_1 = Polygon([*zip(x1, y1)])
rectangle_2 = Polygon([*zip(x2, y2)])
intersection = rectangle_1.intersection(rectangle_2)
union = rectangle_1.union(rectangle_2)
difference_1_2 = rectangle_1.difference(rectangle_2)
difference_2_1 = rectangle_2.difference(rectangle_1)

However, when plotting this polygons in Matplotlib I could not find a direct method where I can plot the exterior and interior paths that can exist in shapely Polygons. In particular being able to plot 'holes' in a bigger polygon created by differences of smaller polygons fully embedded in the bigger one.

enter image description here

I came up with the following solution whereby the exterior and interior paths of the polygon are drawn using the the mpl Path module and convert this to a patch.

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from shapely.geometry import Polygon
class PatchPolygon:
 def __init__(self, polygon, **kwargs):
 polygon_path = self.pathify(polygon)
 self._patch = PathPatch(polygon_path, **kwargs)
 @property
 def patch(self):
 return self._patch
 @staticmethod
 def pathify(polygon):
 ''' Convert coordinates to path vertices. Objects produced by Shapely's
 analytic methods have the proper coordinate order, no need to sort.
 The codes will be all "LINETO" commands, except for "MOVETO"s at the
 beginning of each subpath
 '''
 vertices = list(polygon.exterior.coords)
 codes = [Path.MOVETO if i == 0 else Path.LINETO
 for i in range(len(polygon.exterior.coords))]
 for interior in polygon.interiors:
 vertices += list(interior.coords)
 codes += [Path.MOVETO if i == 0 else Path.LINETO
 for i in range(len(interior.coords))]
 return Path(vertices, codes)
x1 = [3, 6, 6, 3]
y1 = [2, 2, 4, 4]
x2 = [1, 1, 2, 2]
y2 = [1, 2, 2, 1]
x3 = [0, 9, 9, 0]
y3 = [0, 0, 6, 6]
fig, ax1 = plt.subplots(figsize=(6, 4))
ax1.set_xlim(-1, 10)
ax1.set_ylim(-1, 10)
rectangle_1 = Polygon([*zip(x1, y1)])
rectangle_2 = Polygon([*zip(x2, y2)])
rectangle_3 = Polygon([*zip(x3, y3)])
difference = rectangle_3.difference(rectangle_2).difference(rectangle_1)
ax1.add_patch(PatchPolygon(difference, facecolor='blue', edgecolor='red').patch)
plt.show()

Credit goes to Sean Gillies (painting-punctured-polygons-with-matplotlib) who inspired me to the solution. However we are now almost 10 years further, so wondered if there is a better way!

• Mast: Please do not update the code in your question to incorporate feedback from answers, doing so goes against the Question + Answer style of Code Review. This is not a forum where you should keep the most updated version in your question. Please see what you may and may not do after receiving answers .

Answer 1

Unfortunately I don't know of a better way (but have also not looked into one, tbh). Here are some comments on your pathify method instead.

You are currently checking if i == 0 on every loop iteration, but it is true only once at the beginning. Just special case that one out and do

codes = [Path.MOVETO]
codes += [Path.LINETO for _ in range(1, len(polygon.exterior.coords))]

Even better, lists can be multiplied. Since Path.LINETO seems to be a constant, you can just do

codes = [Path.MOVETO] + [Path.LINETO] * (len(polygon.exterior.coords) - 1)

For the vertices I would use list.extend instead of list addition. This way you don't need to cast it to a list yourself but can let extend handle it. It might even be implemented in such a way that it just consumes an iterable. CPython actually does this, but it guesses the size by which the original list needs to be extended by asking the iterable object. If the object does not return a guess, it uses 8 instead (which might not be the most efficient for e.g. a long generator).

for interior in polygon.interiors:
 vertices.extend(interior.coords)
 codes.extend([Path.MOVETO] + [Path.LINETO] * (len(interior.coords) - 1))

At this point it might make sense to put generating the commands into a function, since we already had to use it twice:

@staticmethod
def generate_codes(n):
 """ The first command needs to be a "MOVETO" command,
 all following commands are "LINETO" commands.
 """
 return [Path.MOVETO] + [Path.LINETO] * (n - 1)
@staticmethod
def pathify(polygon):
 ''' Convert coordinates to path vertices. Objects produced by Shapely's
 analytic methods have the proper coordinate order, no need to sort.
 The codes will be all "LINETO" commands, except for "MOVETO"s at the
 beginning of each subpath
 '''
 vertices = list(polygon.exterior.coords)
 codes = self.generate_codes(len(polygon.exterior.coords))
 for interior in polygon.interiors:
 vertices.extend(interior.coords)
 codes.extend(self.generate_codes(len(interior.coords)))
 return Path(vertices, codes)

You should also probably put your calling code under a if __name__ == "__main__": guard to allow importing from this module without running the example code.

• Bruno Vermeulen: Thanks for the suggestions. I will have a look at extend that I had forgotten about! you are right on the if name == 'main' guard which I will edit in the text as it is a bit of a distraction from the real issue.
• Graipher: @BrunoVermeulen: Please don't edit your code after having received an answer, it would invalidate it! Have a look at what you are allowed and not allowed to do after you have received an answer. I rolled back your edit. In any case, since I mentioned it, any other reviewers might just ignore it, or put it in anyway (in which case you can just ignore it).