Computational Geometry in Python
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.