# Geometries (shapely)#

Last updated: 2023-02-25 13:43:42


## Introduction#

In the previous chapters we covered the basics of working with Python, starting with the working environment setup (see Setting up the environment) to advanced methods of working with tables (Table reshaping and joins). Now, we move on to the main topic of this book, working with spatial data.

Spatial data can be divided into two categories:

• Vector layers—points, lines, and polygons, associated with attributes

• Rasters—numeric arrays representing a regular grid over a rectangular area

In the next three chapters, we go through the basics of working with the first type, vector layers, in Python. In this chapter, we cover the shapely package, which is used to represent and work with individual vector geometries. The next two chapters (see Vector layers (geopandas) and Geometric operations) deal with another package, called geopandas, which is used to represent an work with vector layers. The individual geometries within a vector layer are stored as shapely geometries. Therefore it is essential to be familiar with the shapely package, first, before moving on geopandas.

Specifically, in this chapter we are going to cover the following topics:

## What is shapely?#

shapely is a Python package for working with vector geometries, that is, the geometric component of vector layers (the other component being non-spatial attributes). shapely includes functions for creating geometries, as well as functions for applying geometric operations on geometries, such as calculating the centroid of a polygon.

Technically, shapely is a Python interface to the Geometry Engine Open Source (GEOS) software. GEOS is used in numerous open-source programs other than shapely, such as QGIS, and interfaced in many programming languages, such as the sf package in R.

By design, shapely only deals with individual geometries, their creation, their derived properties, and spatial operation applied to them. shapely does not include any functions to read or write geometries to or from files, or more complex data structures to represent collections of geometries with or without non-spatial attributes (i.e., vector layers). These are accomlished with geopandas, which builds on top of pandas and shapely (Fig. 5), to provide more high-level methods to work with complete vector layers.

Note

The shapely documentation (Fig. 36) is quite detailed, and has many illustrations which are useful to understand how the various functions work.

## What are Simple Features?#

The GEOS software, and accordingly also the shapely package, implements the Simple Features specification of geometry types. The geometry type specification and functions may be familiar to you, since both the GEOS software and the Simple Features standard are widely used in open-source GIS software and in other programming languages, such as R’s sf package and the PostGIS spatial database.

The Simple Features defines at least 17 geometry types, but only seven types are commonly used in spatial analysis. These are shown in Fig. 37. Note that there are separate geometry types for points ("Point"), lines ("LineString"), and polygons ("Polygon"), as well are their multi-part versions ("MultiPoint", "MultiLineString", "MultiPolygon"), which sums up to the first six geometry types. The multi-part versions are intended for geometries that contain more than one single-part geometry at once. The seventh type ("GeometryCollection") is intended for any combination of the other six types.

Well Known Text (WKT) is a plain text format for representing Simple Feature geometries, and part of the Simple Feature standard. WKT is understood by many programs, including shapely, the sf package in R, and PostGIS. WKT strings of the geometries shown in Fig. 37 are given in Table 17 1.

Table 17 Simple feature geometry types and WKT examples#

Type

WKT example

"Point"

POINT (30 10)

"LineString"

LINESTRING (30 10, 10 30, 40 40)

"Polygon"

POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10),(20 30, 35 35, 30 20, 20 30))

"MultiPoint"

MULTIPOINT (10 40, 40 30, 20 20, 30 10)

"MultiLineString"

MULTILINESTRING ((10 10, 20 20, 10 40),(40 40, 30 30, 40 20, 30 10))

"MultiPolygon"

MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)),((20 35, 10 30, 10 10, 30 5, 45 20, 20 35),(30 20, 20 15, 20 25, 30 20)))

"GeometryCollection"

GEOMETRYCOLLECTION (POINT (40 10),LINESTRING (10 10, 20 20, 10 40),POLYGON ((40 40, 20 45, 45 30, 40 40)))

Do not worry if you are not familiar with the exact syntax of the WKT format. You will have an opportunity to see various WKT strings in the examples in this chapter. For now, just note that these are strings specifying the geometry type and the coordinates comprising the geometry. Parentheses (() are used the indicate the nested structure of coordinate sequences withing the geometry, such as the separate polygons comprising a "MultiPolygon" geometry, or the exterior and interior (i.e., holes) in a "Polygon". For example:

• POINT (30 10)—A point at $$(30,10)$$

• LINESTRING (0 0, 0 1)—A straight line going from $$(0,0)$$ to $$(0,1)$$

## Creating geometries#

### shapely from WKT#

Before starting to work with shapely, we need to load the package. The shapely package is organized into sub-packages, of which we need two. Specifically, we load the sub-package shapely.geometry which contains all geometry-related functions we will be working with. We also load the shapely.wkt sub-package, which is necessary to convert Well Known Text (WKT) strings to shapely geometries:

import shapely.geometry
import shapely.wkt


One way to create a shapely geometry object is to import a WKT string. To convert a WKT string into a shapely geometry, we pass the WKT string to the shapely.wkt.loads function (short for “load string”). In our first example, we create a polygon named pol1. Note that this is a rectange, composed of five coordinates. The fifth coordinate is identical to the first, used to “close” the shape, as required by the WKT syntax:

pol1 = shapely.wkt.loads('POLYGON ((0 0, 0 -1, 7.5 -1, 7.5 0, 0 0))')


Typing a shapely geometry in a Jupyter notebook conveniently draws a small figure:

pol1


We can also print the geometry, which prints its WKT string representation:

print(pol1)

POLYGON ((0 0, 0 -1, 7.5 -1, 7.5 0, 0 0))


The type of a shapely geometry consists of several components, starting from the general class, shapely.geometry, and ending with the specific geometry type, such as "Polygon":

type(pol1)

shapely.geometry.polygon.Polygon


Let us create another polygon geometry, named pol2:

pol2 = shapely.wkt.loads('POLYGON ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1))')
print(pol2)

POLYGON ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1))


This time we have more coordinates, and, accordingly, the shape of the polygon is more complex:

pol2


Finally, let us create another type of geometry, a "MultiPolygon", named pol3. Note how the "MultiPolygon" has several sequences of shapes, representing the outline and hole(s) of one or more polygons:

pol3 = shapely.wkt.loads('''
MULTIPOLYGON
(((40 40, 20 45, 45 30, 40 40)),
((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20)))
''')


Note

The above expression uses multi-line string input in Python, using the ''' (or """) character at the start and end of the string. This tells the Python interpreter that the string spans several lines. Otherwise we would get an error since a “newline” character ends the expression.

In pol3, we have one polygon without any holes, and another polygon that has one hole:

pol3


Exercise 07-a

• Go to the WKT article on Wikipedia.

• Locate examples of other types of geometries, such as "Point", "LineString", "MultiPoint", and "MultiLineString". Create the corresponding shapely geometries.

• Print the shapes to inspect them graphically. Compare to the figures on Wikipedia.

### shapely from type-specific functions#

Another method of creating shapely geometries is to pass lists of coordinates to specific functions, according to the required geometry type (Table 18).

Table 18 shapely functions to create geometries#

Function

Geometry type

shapely.geometry.Point

"Point"

shapely.geometry.MultiPoint

"MultiPoint"

shapely.geometry.LineString

"LineString"

shapely.geometry.MultiLineString

"MultiLineString"

shapely.geometry.Polygon

"Polygon"

shapely.geometry.MultiPolygon

"MultiPolygon"

shapely.geometry.GeometryCollection

"GeometryCollection"

Each of the function takes input of a different form, according to the geometry type. For example, to create a "Point" geometry, we need to pass a tuple of the form (x,y) to shapely.geometry.Point:

pnt1 = shapely.geometry.Point((2, 0.5))
pnt1


A "MultiPoint", however, can be created from a list of tuples, where each list elements represents an individual point, using shapely.geometry.MultiPoint:

coords = [(2, 0.5), (1, 1), (-1, 0), (1, 0)]
pnt2 = shapely.geometry.MultiPoint(coords)
pnt2


A "LineString" is also created from a list of tuples, using shapely.geometry.LineString:

coords = [(2, 0.5), (1, 1), (-1, 0), (1, 0)]
line1 = shapely.geometry.LineString(coords)
line1


A "MultiLineString" can be created from a list of "LineString" geometries, using shapely.geometry.MultiLineString:

l1 = shapely.geometry.LineString([(2, 0.5), (1, 1), (-1, 0), (1, -1)])
l2 = shapely.geometry.LineString([(-2, 1), (2, -0.5), (3, 1)])
line2 = shapely.geometry.MultiLineString([l1, l2])
line2


Exercise 07-b

• The shapely.geometry.GeometryCollection function accepts a list of geometries, and returns a "GeometryCollection" composed of them.

• Use the above objects to construct and plot a "GeometryCollection" that contains a point, a line, and a polygon.

A "Polygon" can be created using the shapely.geometry.Polygon function and using a coordinate sequence where the first and last points are identical:

coords = [(0, 0), (0, -1), (7.5, -1), [7.5, 0], (0, 0)]
shapely.geometry.Polygon(coords)


Optionally, a second list may used to pass a list of sequences comprising polygon holes. For example, the following modified expressions create a polygon with one hole:

coords_exterior = [(0, 0), (0, -1), (7.5, -1), [7.5, 0], (0, 0)]
coords_interiors = [[(0.4, -0.3), (5, -0.3), (5, -0.7), (0.4, -0.7), (0.4, -0.7)]]
shapely.geometry.Polygon(coords_exterior, coords_interiors)


A "MultiPolygon" can be created from a list of individual "Polygon" geometries, using shapely.geometry.MultiPolygon:

shapely.geometry.MultiPolygon([pol1, pol2])


Note that what we got in the last expression is an individual "MultiPolygon" geometry, not two separate polygon geometries:

print(shapely.geometry.MultiPolygon([pol1, pol2]))

MULTIPOLYGON (((0 0, 0 -1, 7.5 -1, 7.5 0, 0 0)), ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1)))


Note

The above "MultiPolygon" geometry is considered invalid, i.e., not in agreement with the Simple Features specification (see What are Simple Features?). Validity can be examined using the .is_valid property. It is also suggested by the green/red color of the geometry visualization. Specifically, a valid "MultiPolygon" geometry cannot overlap in an infinite number of points.

### shapely from GeoJSON-like dict#

Our third method to create shapely geometries is the shapely.geometry.shape function. This function accepts a GeoJSON-like dict with the geometry definition. In general, the dict needs to have two properties:

• The "type" property contains the geometry type

• The "coordinates" property contains the sequences of geometries/parts/coordinates, as lists or tuples

For example, the following expression uses shapely.geometry.shape to create a "Point" geometry:

d = {'type': 'Point', 'coordinates': (0, 1)}
shapely.geometry.shape(d)


The following expression uses shapely.geometry.shape to create a "MultiPolygon" geometry. Accordingly, the input dict is more complex:

d = {
'type': 'MultiPolygon',
'coordinates': [
[
[[40, 40], [20, 45], [45, 30], [40, 40]]
],
[
[[20, 35], [10, 30], [10, 10], [30, 5], [45, 20],
[20, 35]],
[[30, 20], [20, 15], [20, 25], [30, 20]]
]
]
}
shapely.geometry.shape(d)


## Geometry type#

As demonstrated above, the information contained in a shapely geometry is basically composed of two items:

• The geometry type, such as "Point"

• The coordinate sequence(s) of all shapes comprising the geometry, such as (0,1)

Using shapely methods, we can access the two items. In this section we show how to access geometry type, which is simple enough. In the next section, we discuss the methods to access the coordinates, which is more complex (see Geometry coordinates).

The .geom_type property of a shapely geometry contains the geometry type:

pol1.geom_type

'Polygon'


This is a string:

type(pol1.geom_type)

str


## Geometry coordinates#

### Overview#

The other component of a shape (in addition to type, see Geometry type) is the shape coordinates. In the following sub-sections, we demonstrate the ways we can access the coordinates of the five most commonly used geometry types (Table 19).

Table 19 Accessing shapely geometry coordinates#

Geometry type

Access method

Point coordinates

.coords

LineString coordinates

.coords

MultiLineString coordinates

.geoms[0].coords

Polygon coordinates

.exterior.coords, .interiors[0].coords

MultiPolygon coordinates

.geoms[0].exterior.coords, .geoms[0].interiors[0].coords

As we will see, the access method is specific per geometry types, due to their different levels of complexity.

In practice, it is rarely necessary to access the coordinates directly, because there are predefined functions for all common operations involving the coordinates (such as calculating the centroid, see Centroids (shapely)). However, being aware of the coordinate access methods gives you the freedom, if necessary, to develop any type of geometric calculation or algorithm, even if it is not defined in shapely, or any other software. Additionally, through the access methods we will become more familiar with the structure and differences between the geometry types.

### Point coordinates#

"Point" geometry coordiantes are accessible thought the .coords property. The returned object is of a special type called CoordinateSequence:

pnt1.coords

<shapely.coords.CoordinateSequence at 0x7fcd58186cb0>


It can be converted to a more convenient form, of a list, using the list function:

list(pnt1.coords)

[(2.0, 0.5)]


As we can see, coordinates are organized as a list of tuples. In this case, since with are talking about a Point geometry, the list is of length 1, containing a tuple of the point coordinates.

### LineString coordinates#

The coordinates of a "LineString" are accessible in exactly the same way as those of "Point" (see Point coordinates).

list(line1.coords)

[(2.0, 0.5), (1.0, 1.0), (-1.0, 0.0), (1.0, 0.0)]


However, the list is not of length 1, but at least of length 2. For line1, the list is of lenght 4, because the line is composed of four points.

### MultiLineString coordinates#

For "MultiLineString" geometries, there is no direct .coords method to access the coordinates of all parts at once. Instead, the parts have to be accessed one by one, using the .geoms property, combined with an index (starting from 0). Any particular “part” of a "MultiLineString" geometry behaves the same way as a "LineString" geometry. Once we are “inside” a particular "MultiLineString" part, we can use any of the properties or methods applicable to a "LineString" geometry, such as coords.

For example, here is how we can access the two parts of the "MultiLineString" named line2:

line2.geoms[0]

line2.geoms[1]


Note that each part of a "MultiLineString" is a "LineString":

type(line2.geoms[0])

shapely.geometry.linestring.LineString


In case we do not know in advance how many parts there are, we can use the len function to find out:

len(line2.geoms)

2


Going further, we can access the coordinates of each part using .coords (see LineString coordinates):

list(line2.geoms[0].coords)

[(2.0, 0.5), (1.0, 1.0), (-1.0, 0.0), (1.0, -1.0)]

list(line2.geoms[1].coords)

[(-2.0, 1.0), (2.0, -0.5), (3.0, 1.0)]


### Polygon coordinates#

"Polygon" coordinate sequences are associated with the concepts of exterior and interiors coordinate sequences:

• The exterior is the outer bound of the polygon. There is only one exterior sequence per polygon. The .exterior property gives direct access to the exterior geometry.

• The interiors refers to polygon holes. There can be zero or more holes in each polygon. Accordingly, the .interiors property is a sequence of length zero or more.

For example, pol1 is a "Polygon" with no holes. We can access the exterior coordinates using .exterior, and then .coords:

pol1.exterior

list(pol1.exterior.coords)

[(0.0, 0.0), (0.0, -1.0), (7.5, -1.0), (7.5, 0.0), (0.0, 0.0)]


However, the .interiors sequence is empty:

len(pol1.interiors)

0


You will get a chance to see a non-empty .interiors sequence in the next example (see MultiPolygon coordinates).

### MultiPolygon coordinates#

"MultiPolygon" coordinates are accessed by first extracting one of the polygon part, using .geoms, similarly to the way we extract "MultiLineString" parts (see MultiLineString coordinates). Each of the "MultiPolygon" parts is a "Polygon" geometry. Therefore, we can use the exterior and interiors methods to get the coordinates of each part, just like with a standalone "Polygon" (see Polygon coordinates).

For example, pol3 is a "MultiPolygon" geometry with two parts, whereas one of the parts has one hole:

pol3


The .geoms method lets us access the parts one by one:

pol3.geoms[0]

pol3.geoms[1]


Note that each part of a "MultiPolygon" geometry is a "Polygon" geometry:

type(pol3.geoms[0])

shapely.geometry.polygon.Polygon


Starting with the first part, here is the exterior:

pol3.geoms[0].exterior

list(pol3.geoms[0].exterior.coords)

[(40.0, 40.0), (20.0, 45.0), (45.0, 30.0), (40.0, 40.0)]


As expected, there are no interiors:

len(pol3.geoms[0].interiors)

0


Now, here is the second part of the "MultiPolygon":

pol3.geoms[1]


Here is the exterior:

pol3.geoms[1].exterior


and its coordinates:

list(pol3.geoms[1].exterior.coords)

[(20.0, 35.0),
(10.0, 30.0),
(10.0, 10.0),
(30.0, 5.0),
(45.0, 20.0),
(20.0, 35.0)]


This part does contain interiors, specifically one “interior” (i.e., hole):

len(pol3.geoms[1].interiors)

1


Here is how we can access the interior:

pol3.geoms[1].interiors[0]


and its coordinates:

list(pol3.geoms[1].interiors[0].coords)

[(30.0, 20.0), (20.0, 15.0), (20.0, 25.0), (30.0, 20.0)]


## Derived properties#

### Overview#

The geometry type (Geometry type), and the geometry coordinates along with their division to parts and interior/exterior sequences (Geometry coordinates), are basic “internal” geometry properties. Additionally, there are derived useful properties that can be calculated for each geometry.

In the next few section we cover several useful derived geometry properties, namely:

### Bounds (shapely)#

The .bounds property of a shapely geometry contains its bounds, i.e., the coordinates of the bottom-left $$(x_{min}, y_{min})$$ and top-right $$(x_{max}, y_{max})$$ corner. The four coordinates are returned in the form of a tuple of length four $$(x_{min}, y_{min}, x_{max}, y_{max})$$. For example:

pol1.bounds

(0.0, -1.0, 7.5, 0.0)

pol2.bounds

(0.0, -0.5, 7.0, 1.0)

pol3.bounds

(10.0, 5.0, 45.0, 45.0)


Note that even a "Point" geometry has bounds, where $$x_{min}$$ is equal to $$x_{max}$$, and $$y_{min}$$ is equal to $$y_{max}$$:

pnt1.bounds

(2.0, 0.5, 2.0, 0.5)


The shapely.geometry.box function can be used to create a "Polygon" geometry reflecting the bounding box. For example, the following expression creates the bounding box geometry of pol3:

shapely.geometry.box(*pol3.bounds)


Note that the shapely.geometry.box function requires four separate positional arguments (i.e., box(xmin,ymin,xmax,ymax)), rather than a tuple. Therefore, in the last expression, we used the positional argument technique (see Variable-length positional arguments) to “unpack” the pol3.bounds tuple into four arguments passed to box. The equivalent expression without using the positional argument technique would be much longer:

shapely.geometry.box(pol3.bounds[0], pol3.bounds[1], pol3.bounds[2], pol3.bounds[3])


Note

The .envelope property (Table 20) can be used to directly access the bounding box geometry.

### Length (shapely)#

The .length property gives the length of the geometry. For example, here we calculate the length of line1, a "LineString" geometry:

line1

line1.length

5.354101966249685


and here we calculate the length of line2, which is a "MultiLineString" geometry:

line2

line2.length

11.664947454140234


Note that the length (or area, see Area (shapely)) of a multi-part geometry is the sum of lengths (or areas) of all its parts. For example, using the access methods shown above (see Geometry coordinates), we can calculate the length of the two parts that comprise the "MultiLineString" named line2:

line2.geoms[0].length

5.5901699437494745

line2.geoms[1].length

6.07477751039076


and demonstrate that the sum of lengths of the two parts is equal to the total geometry length:

line2.geoms[0].length + line2.geoms[1].length == line2.length

True


Keep in mind that:

• The .length of "Point" or "MultiPoint" geometries is equal to 0

• The .length of "Polygon" or "MultiPolygon" geometries refers to the summed perimeters of the exterior and all interiors

Exercise 07-c

• pol1 is a "Polygon" geometry representing a rectangle in parallel to the x- and y-axes.

• Calculate the rectangle perimeter “manually” using the formula $$(|x_{max} - x_{min}|)*2 + (|y_{max} - y_{min}|)*2$$.

• Demonstrate that the result is equal to pol1.length.

### Area (shapely)#

The .area property gives the area of a given geometry. For example, here is how we can get the areas of pol1, pol2, and pol3:

pol1.area

7.5

pol2.area

6.25

pol3.area

712.5


Keep in mind that the area of "Point", "MultiPoint", "LineString", and "MultiLineString" geometries is equal to 0!

Exercise 07-d

• Write expressions to calculate the area of each of the two parts that comprise the "MultiPolygon" named pol3.

• Write an expression that demonstrates that the sum of the two areas is equal to pol3.area (similarly to what we did for length of line2, see Length (shapely)).

## New geometries 1#

### Overview#

The shapely package contains many commonly used methods to create new geometries based on existing ones (Table 20). These are defined as properties, or methods, of an existing geometry. In the following sections, we go over three most common ones: .centroid, .buffer(), and .convex_hull.

Table 20 shapely functions to create a new geometry from an input individual geometry#

Method

Meaning

.centroid

Centroids (shapely)

.buffer()

Buffer (shapely)

.convex_hull

Convex hull

.envelope

Bounding box

.minimum_rotated_rectangle

Minimum rotated rectangle

.simplify

Simplified geometry

### Centroids (shapely)#

The centroid is defined as the center of mass of a geometry. The centroid of a given geometry is obtained through the .centroid property. For example:

pol2.centroid


To visualize the centroid along with the original geometry, we may take advantage of the "GeometryCollection" geometry type, as follows:

shapely.geometry.GeometryCollection([pol2, pol2.centroid])


### Buffer (shapely)#

A buffer is a geometry that encompasses all points within a given distance from the input geometry. In shapely, a buffer is calculated using the .buffer method, given the specified distance. For example, here is a buffer 5 units of distance around pol1:

pol1.buffer(5)


Here is another example, pol2 buffered to a distance of 0.75:

pol2.buffer(0.75)


Keep in mind that shapely does not have a notion of Coordinate Reference Systems (CRS). The buffer distance, as well as the distances in any other distance-related calculation, are specified in the same (undefined) units as the geometry coordinates. CRS are defined “on top” of shapely geometries in higher-level packages such as geopandas (see Coordinate Reference System (CRS) and CRS and reprojection).

Also note that the .buffer distance can be negative, in which case the buffer is “internal” rather than “external”.

Exercise 07-e

• What do you think will happen when we use a negative buffer which is large enough so that the geometry “disappears”?

• Use the plot and print methods to see how shapely represents such an empty polygon.

### Convex hull#

A Convex Hull is the minimal convex polygon that encompasses the given set of points (or coordinates of a non-point geometry). The convex hull can be obtained through the convex_hull property of any shapely geometry. For example, here is the convex hull of pol3:

pol3.convex_hull


## Boolean functions#

shapely includes the standard geometric functions to evaluate relations between geometries. These are defined as methods of one geometry, whereas the argument is another geometry (Table 21).

Table 21 shapely boolean methods#

Method

Meaning

.equals

Equals

.almost_equals

Approximately equal

.contains

Contains

.covers

Covers

.covered_by

Covered by

.crosses

Crosses

.disjoint

Disjoint

.intersects

Intersects

.overlaps

Overlaps

.touches

Touches

.within

Within

For example, the .intersects method evaluates whether one geometry intersects with another. Here, we examine whether pol1 intersects with pol3:

pol1.intersects(pol3)

False


The answer is False, meaning that the two geometries do not intersect, as evident in the following figure:

shapely.geometry.GeometryCollection([pol1, pol3])


However, pol1 does intersect with pol2:

pol1.intersects(pol2)

True


as demonstrated below:

shapely.geometry.GeometryCollection([pol1, pol2])


The .intersects method is probably the most useful one, since we are usually interested in simply finding out whether the given two geometries intersect or not. There are other functions referring to more specific relations, such as .equals, .contains, .crosses, and .touches (Table 21).

Note

See the Binary Predicates section in the shapely documentation for a complete list and more details on boolean functions for evaluating geometry relations.

## New geometries 2#

### Overview#

Earlier, we dealt with methods to create a new geometry given an individual input geometry, such as calculating a centroid, a buffer, or a convex hull (see New geometries 1). shapely also defines the standard methods to create an new geometry based on a pair of input geometries, which is what we cover in this section (shapely functions to create a new geometry from an input pair of geometries).

Table 22 shapely functions to create a new geometry from an input pair of geometries#

Method

Meaning

.intersection

Intersection

.union

Union

.difference

Difference

.symmetric_difference

Symmetric difference

For example, given two partially overlapping circles, x (left) and y (right):

x = shapely.geometry.Point((0, 0)).buffer(1)
y = shapely.geometry.Point((1, 0)).buffer(1)
shapely.geometry.GeometryCollection([x, y])


Here is the intersection:

x.intersection(y)


Here is the difference of x from y:

x.difference(y)


And here is the union:

x.union(y)


In the next three sub-sections we elaborate on each of those methods.

Note

Another function of the same “group”, which we do not cover here, is .symmetric_difference. See the documentation: https://shapely.readthedocs.io/en/stable/manual.html#object.symmetric_difference.

### Intersection#

The intersection of two geometries is defined, in plain terms, as their shared area. Given two geometries x and y, the intersection is given by the .intersection method, as in x.intersection(y). Note that, since the intersection operator is symmetric, y.intersection(x) gives the same result.

As an example, let us calculate the area of intersection of pol1 and pol2:

pol1.intersection(pol2)


These two polygons were intentionally 2 constructed in such a way that their intersection results in three parts: a "Point", a "LineString", and a "Polygon". Since the result must still be contained in a single geometry object, this must be a "GeometryCollection" (Table 17):

print(pol1.intersection(pol2))

GEOMETRYCOLLECTION (POLYGON ((7 0, 7 -0.5, 6 -0.5, 5.5 0, 7 0)), LINESTRING (4 0, 3 0), POINT (1 0))


Note that the intersection of two geometries may be an empty geometry, in case the geometries do not have even one point in common:

print(pol1.intersection(pol3))

POLYGON EMPTY


### Union#

The .union method, as in x.union(y), returns the area of covered by at least one of the geometries x and y (or both). For example, recall pol1 and pol2:

shapely.geometry.GeometryCollection([pol1, pol2])


Here is their union:

pol1.union(pol2)


Like .intersection (see Intersection), the .union operation is symmetric. For example, pol2.union(pol1) gives the same result as above:

pol2.union(pol1)


### Difference#

The third commonly used pairwise operator to create a new geometry is the .difference method. The operation x.difference(y) returns the area of x minus the area of y. For example:

x.difference(y)  # The area of 'x' minus the area of 'y'


Unlike .intersection and .union, the .difference method is asymmetric. Namely, y.difference(x) is not the same as x.difference(y):

y.difference(x)  # The area of 'y' minus the area of 'x'


Exercise 07-f

• Calculate the difference between the bounding box of pol3, and pol3 (Fig. 38).

• Create a polygon where the pol1 rectangular polygon comprises a polygon hole. To do that, apply a buffer on pol1, then “erase” the the original pol1 from the buffered one using .difference (Fig. 39).

## Distance (shapely)#

Another very common operation in spatial analysis is the calculation of distance. As part of our analysis we often need to figure out how far are the given geometries from each other. For example, we may need to calculate what is the distance of a particular city from a source of air pollution such as a power plants, and so on.

In shapely, the distance between two geometries x and y is calculated using the .distance method, as in x.distance(y). For example, the following expression returns the distance between pol1 and pol3:

pol1.distance(pol3)

10.307764064044152


Exercise 07-g

• Do you think the .distance operator is symmetric?

• Execute the reverse pol3.distance(pol1) to check your answer.

Keep in mind that distance is defnied as “shortest” distance, i.e., the minimal distance between any pair of points on the two geometries.

As can be expected, the distance between two geometries that intersect is zero. For example:

pol1.distance(pol2)

0.0


## Points to line (shapely)#

Geometry transformations, sometimes called “geometry casting”, are operations where we convert a geometry from one type to another. Not all combinations make sense, and not all are useful. For example, there is no sensible way to transform a "Point" to a "Polygon", because a valid polygon needs to be composed of at least three coordinates (a triangle). A "Polygon", however, can always be transformed to a "MultiPoint", which represents the the polygon coordinates as separate points.

One of the most useful geometry casting operations, which we will come back to later on in the book (see Points to line (geopandas)) is the transformation of a "MultiPoint" geometry, or a set of "Point" geometries, to a "LineString" geometry which “connects” the points according to their given order. For example, a series of chronologically ordered GPS readings, may be connected to a "LineString" representing the recorded route. This type of casting facilitates further analysis, such as visualizing the route, calculating its total length, figuring out which areas the route goes through, and so on.

We have already seen that the shapely.geometry.LineString function accepts a list of coordinates (shapely from type-specific functions). Alternatively, the shapely.geometry.LineString function can accept list of "Point" geometries, and “transform” it into a "LineString". To demonstrate, let us create five "Point" geometries named p1, p2, p3, p4, and p5, using wkt.loads (shapely from WKT):

p1 = shapely.wkt.loads('POINT (0 0)')


Next, we combine the geometries into a list named p:

p = [p1, p2, p3, p4, p5]
p

[<shapely.geometry.point.Point at 0x7fcd58186920>,
<shapely.geometry.point.Point at 0x7fcd58185d20>,
<shapely.geometry.point.Point at 0x7fcd58186770>,
<shapely.geometry.point.Point at 0x7fcd58186e30>,
<shapely.geometry.point.Point at 0x7fcd58184c40>]


Now we can use the shapely.geometry.LineString function to transform the list of point geometries to a "LineString" geometry:

shapely.geometry.LineString(p)


We will come back to this technque for a more realistic example (see Points to line (geopandas)), where we are going to calculate a "LineString" geometryies of public transport routes based on the GTFS dataset (see What is GTFS?).

## More exercises#

Exercise 07-h

• Create two geometries—a square and a circle—in such a way that they partially overlap, so that $$1\over4$$ of the circle area overlaps with the square (Fig. 40). Hint: use the buffer function applied on a point geometry to create the circle.

• Calculate and plot the geometry representing the union of the square and the circle (Fig. 41).

• Calculate and plot the geometry representing the intersection of the square and the circle (Fig. 42).

• Calculate the area of the intersection and of the the union geometries.

Exercise 07-i

• Execute the code section given below to get a DataFrame with the location of bus stops along bus line 24 in Beer-Sheva, the route which goes from Ramot to Beer-Sheva Central Bus Station. The bus stop coordinates are in the "stop_lon" and "stop_lat" columns.

• Calculate and plot a MultiPoint geometry with all stop locations (Fig. 43). Hint: extract the contents of the "stop_lon" and "stop_lat" columns as a list of lists/tuples, then use appropriate function to create a MulpiPoint geometry (see shapely from type-specific functions).

• Calculate and plot a LineString geometry with of the bus route (Fig. 44).

• Calculate and plot a LineString geometry with the straight line between the first and last bus stops (Fig. 45).

• Calculate and plot the convex hull polygon (see Convex hull) of the bus stops (Fig. 46).

import pandas as pd
stops = stops[["stop_id", "stop_name", "stop_lon", "stop_lat"]]
stop_times = stop_times[["trip_id", "stop_id", "stop_sequence"]]
sel = stop_times["trip_id"] == "55745843_180421"
stop_times = stop_times[sel]
stops = pd.merge(stop_times, stops, on="stop_id", how="left")
stops = stops.sort_values(by="stop_sequence")


1

WKT examples are from Wikipedia (https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry).

2

The polygon coordinates are from an example in the sf R package documentation: https://cran.r-project.org/web/packages/sf/vignettes/sf1.html.