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.

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!

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.