How to extract the intersection between two polygons in Python using Shapely ?

Introduction

In this article, we will discuss how to extract the intersection between two polygons using Shapely - a Python package for geometric manipulation and analysis.

Before diving into the extraction process, let's first understand what polygons are. In simple terms, a polygon is a closed shape with three or more straight sides. It is a two-dimensional figure that is made up of line segments connected end to end.

Examples of polygons include squares, rectangles, triangles, and pentagons.

Create two polygons using Shapely

Shapely is a Python package that provides geometric objects, predicates, and operations for the manipulation and analysis of planar geometric objects.

To find the intersection between two polygons using Shapely, we must start by creating two polygon objects. This is achieved by specifying a list of coordinates that define the vertices of each polygon. Let's begin by creating the first polygon using Shapely:

from shapely.geometry import Polygon

import numpy as np

a = 2.0 # square size

pos = np.array((1.0, 2.0))

c1 = np.array((0.0, 0.0)) + pos
c2 = np.array((0.0, a)) + pos
c3 = np.array((a, a)) + pos
c4 = np.array((a, 0.0)) + pos

coords = (c1,c2,c3,c4)

rect1 = Polygon( coords )

and then a second one:

a = 2.0 # square size

pos = np.array((2.0, 3.0))

c1 = np.array((0.0, 0.0)) + pos
c2 = np.array((0.0, a)) + pos
c3 = np.array((a, a)) + pos
c4 = np.array((a, 0.0)) + pos

coords = (c1,c2,c3,c4)

rect2 = Polygon( coords )

Once we have our polygon objects, we can then use the "intersection" method to find the common area between them:

rect3 = rect1.intersection(rect2)

The "rect3" variable will now contain a new polygon object representing the common area between "rect1" and "rect2".

To obtain the coordinates of the new polygon rect3, one possible solution is to utilize the exterior.coords.xy function:

rect3.exterior.coords.xy

gives

(array('d', [2.0, 3.0, 3.0, 2.0, 2.0]), array('d', [4.0, 4.0, 3.0, 3.0, 4.0]))

Visualization

Plotting the two polygons

In a previous article, we learned about plotting polygons in Python using Shapely and Matplotlib. Now, let's utilize these libraries to plot the two polygons we have defined earlier.

from matplotlib.pyplot import figure

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

fig = figure(num=None, figsize=(12, 10), dpi=100, edgecolor='k')

ax = plt.axes()

#----------------------------------------#
# Plot rectangle 1

X,Y = rect1.exterior.coords.xy

poly = mpatches.Polygon(np.stack( (X, Y), axis=1 ), 
                        closed=True, 
                        ec='k', 
                        lw=1, 
                        color='coral',
                        alpha=0.5)

ax.add_patch(poly)

#----------------------------------------#
# Plot rectangle 2

X,Y = rect2.exterior.coords.xy

poly = mpatches.Polygon(np.stack( (X, Y), axis=1 ), 
                        closed=True, 
                        ec='k', 
                        lw=1, 
                        color='lightblue',
                        alpha=0.5)

ax.add_patch(poly)

 #----------------------------------------#

plt.xlim(0,6)
plt.ylim(0,6)

plt.title('How to extract the intersection between two polygons in Python using Shapely ?')

plt.savefig("polygons_intersection_01.png", bbox_inches='tight', facecolor='white')

plt.show()

Plotting the intersection

from matplotlib.pyplot import figure

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

fig = figure(num=None, figsize=(12, 10), dpi=100, edgecolor='k')

ax = plt.axes()

#----------------------------------------#
# Plot intersection

X,Y = rect3.exterior.coords.xy

poly = mpatches.Polygon(np.stack( (X, Y), axis=1 ), 
                        closed=True, 
                        ec='k', 
                        lw=1, 
                        color='black',
                        alpha=0.5)

ax.add_patch(poly)

#----------------------------------------#

plt.xlim(0,6)
plt.ylim(0,6)

plt.title('How to extract the intersection between two polygons in Python using Shapely ?')

plt.savefig("polygons_intersection_02.png", bbox_inches='tight', facecolor='white')

plt.show()

Plotting all together

from matplotlib.pyplot import figure

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

fig = figure(num=None, figsize=(12, 10), dpi=100, edgecolor='k')

ax = plt.axes()

#----------------------------------------#
# Plot rectangle 1

X,Y = rect1.exterior.coords.xy

poly = mpatches.Polygon(np.stack( (X, Y), axis=1 ), 
                        closed=True, 
                        ec='k', 
                        lw=1, 
                        color='coral',
                        alpha=0.5)

ax.add_patch(poly)

#----------------------------------------#
# Plot rectangle 2

X,Y = rect2.exterior.coords.xy

poly = mpatches.Polygon(np.stack( (X, Y), axis=1 ), 
                        closed=True, 
                        ec='k', 
                        lw=1, 
                        color='lightblue',
                        alpha=0.5)

ax.add_patch(poly)

#----------------------------------------#
# Plot intersection

X,Y = rect3.exterior.coords.xy

poly = mpatches.Polygon(np.stack( (X, Y), axis=1 ), 
                        closed=True, 
                        ec='k', 
                        lw=1, 
                        color='black',
                        alpha=0.5)

ax.add_patch(poly)

#----------------------------------------#

plt.xlim(0,6)
plt.ylim(0,6)

plt.title('How to extract the intersection between two polygons in Python using Shapely ?')

plt.savefig("polygons_intersection_03.png", bbox_inches='tight', facecolor='white')

plt.show()

References