Computational Geometry in Python

David@Pixabay

If you are looking for a good Python tool to do computational geometry, look no further than Shapely! Shapely is a powerful Python package designed to help you process and analyze geometric objects. With Shapely, you can create, manipulate, and analyze 2D and 3D geometric shapes with ease. From simple points and lines to complex polygons and multipolygons, Shapely provides an intuitive and flexible interface for working with spatial data.

In this tutorial, you’ll learn the basics of Shapely and explore many of its features. You’ll discover how to create and manipulate geometric objects, perform geometric operations, and analyze spatial relationships between objects. By the end of this tutorial, you’ll be able to use Shapely to solve problems in computational geometry. Let’s get started!

Installation

Installation of Shapely is super-easy using pip . Just run

pip install shapely

and you are ready to go.

Basic Geometry Objects in Shapely

Shapely provides a rich set of geometric objects for working with spatial data in Python. Let's take a look at some of the basic geometric objects that Shapely provides.

Points

Points are the simplest geometric objects in Shapely. They represent a single point in space and are defined by their x, y (and optionally z) coordinates. Here’s how you can create a point object in Shapely and visualize it using the matplotlib library:

from shapely.geometry import Point
import matplotlib.pyplot as plt

p = Point(2, 3)

plt.scatter(p.x, p.y)
plt.show()

Alternatively, if you don’t need the axes, you can just the variable in the Jupyter cell, and it will display nicely:

p

LineStrings

LineString s are sequences of points that define a line. They can be straight or curved and can have any number of vertices. Here's how you can create a LineString object in Shapely:

from shapely.geometry import LineString
import matplotlib.pyplot as plt

line = LineString([(0, 0), (1, 1)])
plt.plot(line.xy[0], line.xy[1])
plt.show()

If you need the boundary points of the line, you can get them from the boundary attribute.

Polygons

Polygons are closed LineString objects that define a shape with an interior and an exterior. They can be simple or complex and can have any number of vertices. Here's an example:

from shapely.geometry import Polygon
import matplotlib.pyplot as plt
polygon = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
plt.plot(polygon.exterior.xy[0], polygon.exterior.xy[1])
plt.gca().set_aspect("equal")
plt.show()

The exterior attribute is used to access the outer boundary of the polygon. In Shapely, a polygon is defined by its boundary, which is composed of one or more so-called "LinearRing" objects. The outer boundary of a polygon is a LinearRing, while any holes within the polygon would also be represented as LinearRing objects. The exterior attribute of a Polygon object returns a LinearRing object representing the outer boundary of the polygon. You can access the coordinates of the points that make up the boundary by calling the xy attribute on the LinearRing. The xy attribute returns a tuple containing the x-coordinates and y-coordinates of the points that make up the boundary:

polygon.exterior.xy
(array('d', [0.0, 0.0, 1.0, 1.0, 0.0]), array('d', [0.0, 1.0, 1.0, 0.0, 0.0]))

If the polygon had any holes, we could access them using the interiors attribute of the Polygon object. The interiors attribute returns a list of LinearRing objects representing the holes in the polygon. Note that while a shape can have only one exterior boundary, it can have several disconnected interior boundaries.

In addition to these basic geometric objects, Shapely also provides more advanced objects for working with multi-part geometries, such as MultiPoints, MultiLineStrings, and MultiPolygons which we will cover in a later section of this article.

Using these basic geometric objects, you can create and manipulate spatial data in Python with ease. In the next section, we’ll explore how to perform geometric operations on these objects.

Geometric Operations

Once you’ve created geometric objects in Shapely, you can perform a wide range of geometric operations on them. In this section, we’ll explore some of the basic geometric operations that Shapely provides.

Intersection

The intersection of two geometric objects is the portion of space that they have in common. Here’s how you can compute the intersection of two LineStrings:

from shapely.geometry import LineString
line1 = LineString([(0, 0), (1, 1)])
line2 = LineString([(0, 1), (1, 0)])
intersection = line1.intersection(line2)
print(intersection)
POINT (0.5 0.5)

If we had computed the intersection of two parallel lines, we would get

from shapely.geometry import LineString
line1 = LineString([(0, 0), (1, 0)])
line2 = LineString([(0, 1), (1, 1)])
intersection = line1.intersection(line2)
print(intersection)
LINESTRING Z EMPTY

Of course, you can also do this with 2D shapes:

from shapely.geometry import Polygon
import matplotlib.pyplot as plt
triangle1 = Polygon([(0, 0), (2, 0), (1, 1)])
triangle2 = Polygon([(1, 0), (3, 0), (2, 1)])
intersection = triangle1.intersection(triangle2)
fig, ax = plt.subplots()
ax.plot(triangle1.exterior.xy[0], triangle1.exterior.xy[1], label="triangle 1")
ax.plot(triangle2.exterior.xy[0], triangle2.exterior.xy[1], label="triangle 2")
ax.plot(intersection.exterior.xy[0], intersection.exterior.xy[1], label="intersection")
ax.legend()
plt.show()

Union

The union of two geometric objects is the smallest object that contains both objects. Here’s how you can compute the union of two Polygons:

from shapely.geometry import Polygon
import matplotlib.pyplot as plt
polygon1 = Polygon([(0, 0), (0, 1), (1, 1), (1, 0)])
polygon2 = Polygon([(1, 0), (1, 1), (2, 1), (2, 0)])
union = polygon1.union(polygon2)
plt.plot(polygon1.exterior.xy[0], polygon1.exterior.xy[1], color="blue")
plt.plot(polygon2.exterior.xy[0], polygon2.exterior.xy[1], color="orange")
plt.plot(union.exterior.xy[0], union.exterior.xy[1], color="green")
plt.show()

Difference

The difference of two geometric objects is the portion of space that is contained in one object but not the other. Here’s how you can compute the difference of two Polygons:

from shapely.geometry import Polygon
import matplotlib.pyplot as plt
polygon1 = Polygon([(0, 0), (0, 2), (2, 2), (2, 0)])
polygon2 = Polygon([(1, 1), (1, 3), (3, 3), (3, 1)])
difference = polygon1.difference(polygon2)
plt.plot(polygon1.exterior.xy[0], polygon1.exterior.xy[1], ls="--")
plt.plot(polygon2.exterior.xy[0], polygon2.exterior.xy[1], ls="--")
plt.plot(difference.exterior.xy[0], difference.exterior.xy[1])
plt.show()

This will display a plot with a shape that represents the area contained in polygon1 but not in polygon2 . Conversely, we can say

difference = polygon2.difference(polygon1)
plt.plot(polygon1.exterior.xy[0], polygon1.exterior.xy[1], ls="--")
plt.plot(polygon2.exterior.xy[0], polygon2.exterior.xy[1], ls="--")
plt.plot(difference.exterior.xy[0], difference.exterior.xy[1])
plt.show()

Spatial Relationships

In addition to performing basic geometric operations, Shapely also provides a wide range of predicates for testing spatial relationships between geometric objects. These predicates return True or False depending on whether the spatial relationship holds between the two objects.

Here are a few examples:

from shapely.geometry import Point, Polygon

point = Point(1, 1)
polygon = Polygon([(0, 0), (0, 2), (2, 2), (2, 0)])

print(point.within(polygon)) # True

print(point.touches(polygon)) # True

print(point.disjoint(polygon)) # False
True
False
False

These predicates can be very useful for performing more complex spatial analyses on your geometric objects.

By combining these basic geometric operations and spatial predicates, you can perform complex spatial analyses on your geometric objects. You can use these operations to answer questions like "what is the area of the intersection between two polygons?" or "which points lie within a given polygon?"

In addition to working with individual geometric objects, Shapely provides support for working with multi-part geometries. A multi-part geometry is a collection of one or more geometric objects, where each object can be of a different type (e.g. Point, LineString, or Polygon).

Conclusion

In this tutorial, we covered the basics of Shapely and explored many of its features. We started by learning about the basic geometric objects in Shapely, including points, LineStrings, and polygons. We then went on to cover geometric operations such as intersection, union, and difference, and we also discussed spatial relationships between objects.

Throughout the tutorial, we provided code examples and plots to illustrate each concept, making it easier for readers to follow along and understand how Shapely works. We also showed how Shapely can be used for more complex spatial analyses, including working with multi-part geometries.

Overall, Shapely is a powerful tool that can greatly simplify the process of working with spatial data in Python. We hope this tutorial has given you a good introduction to the capabilities of Shapely and has inspired you to explore its many features further.