Shapely and geometry objects

In this lesson, you will learn how to create and manipulate geometries in Python using the Shapely Python Package.

Sources: Parts of this chapter are based on shapely’s documentation and on chapter 3 of Westra E. (2013).

Spatial data model

Fundamental geometric objects (‘simple features’) that can be used in Python with shapely. (Figures by M. W. Toews; cf.Wikipedia’s article on GeoJSON) ::

All geospatial vector data can be described by a combining a set of fundamental geometric objects: points, lines, and polygons are the basic ingredients of working with spatial data.

In Python, the library package shapely has become the standard tool for working with geo-spatial objects, and for carrying out a range of geometric operations. A basic understanding of how shapely works is paramount for using higher-level tools, such as, for instance, geopandas (see lesson 2) that handles entire data sets of geographic information (a.k.a. ‘layers’).

Shapely, like the vast majority of geospatial software, follows the data model set forward in the Open Geospatial Consortium’s Simple Feature Access standard. In this chapter, we use the terminology used by shapely, but the general concepts are applicable much more widely.

Geometric objects are described and constructed by coordinate tuples

Even more fundamental than that, coordinate tuples are what is used to define the geometric fundamentals point, line, and polygon. Coordinate tuples typically are either two-dimensional (x, y), or three-dimensional (x, y, z).

Info: Tuples

A tuple is a Python data structure that consists of a number of values separated by commas. Coordinate pairs are often represented as a tuple. For example, (60.192059, 24.945831) is a tuple with two values, (1,) a tuple with one value (note the comma that distinguishes a one-value tuple from a mathematical expression in parentheses), and tuple([1, 2, 3]) converts (‘casts’) a list into a tuple.

Tuples belong to the sequence data types in Python. Other sequence data types are lists and ranges. Tuples have many similarities with lists and ranges, but they are often used for different purposes.

The main difference between tuples and lists is that tuples are immutable, which means that the contents of a tuple cannot be altered (while lists are mutable; you can, for example, add and remove values from lists).

Points

Point geometries represent a singular point (in two- or three-dimensional Euclidean space). Points are defined by a single coordinate tuple.

LineStrings

LineString geometries (and the related LinearRings) represent lines. They are defined by a sequence of points. By connecting the points in order, a line is formed, if the first and last point are the same, a linear ring. Consequently, to describe a LineString, at least two points are necessary, for a LinearRing at least three (first and last being identical).

Polygons

Finally, Polygon geometries represent an area. A polygon is defined by exactly one LinearRing as its circumference, and any number of additional LinearRings representing holes that are cut out. As such, at minimum, three coordinate tuples are needed to define a Polygon (but it becomes more difficult quickly, as, naturally, the coordinates that define the holes have to lie within the exterior LinearRing, and also must not intersect each other).

Collections of geometric objects

All of the fundamental geometric types can be combined to represent more complex geographic features, such as, for instance, administrative areas that consist of multiple discontinuous regions:

Multipolygons

Tyrol, a federal state of Austria, is a MultiPolygon. (Data:Statistics Austria)

• MultiPoint geometries represent collections of points.

• MultiLineString geometries represent collections of lines.

• MultiPolygon geometries represent collections of polygons.

• GeometryCollection geometries are collections of points, lines, and polygons, as well as multi-points, multi-lines, and multi-polygons.

Using shapely to create, access, manipulate and analyse geometric objects

Shapely can perform many useful operations on geometries, and provides a range of attributes of geometries. For instance, you can: - create lines and polygons from a set of points - calculate the area, length, perimeter, etc., of geometries - perform geometric operations on a set of geometries, for instance, merging (union), subtracting (difference), or calculating the distance between geometries. - query the relationship between geometries, i.e., find out whether two geometries intersect, touch, cross, etc.

Creating Point geometries and accessing their properties

Creating a point geometry object is easy: simply pass coordinates (x, y, and possibly z) to its constructor.

# Import `shapely.geometry.Point` class
from shapely.geometry import Point

# Create `Point` objects:
point1 = Point(2.2, 4.2)
point2 = Point(7.2, -25.1)
point3 = Point(9.26, -2.456)
point4_3D = Point(9.26, -2.456, 0.57)

Let’s see what these variables now contain:

point1

As we can see, Jupyter notebook is able to display the shape directly on the screen.

Alternatively, use a print statement to show the text representation of a shapely geometry object:

print(point1)
print(point4_3D)

POINT (2.2 4.2)
POINT Z (9.26 -2.456 0.57)

This text representation is in ‘Well-Known Text’ (WKT) format, a standard set forward in the Open Geospatial Consortium’s (OGC) Simple Feature Access (see above). This includes the additional letter ‘Z’ that marks the three-dimensional version of a geometry (e.g., our point point4_3D).

Let’s also check the data type of a point:

type(point1)

shapely.geometry.point.Point

We can see that the type of the point is a shapely.geometry.point.Point (which is equivalent to shapely.geometry.Point, the class we used to instantiate the point object).

Under the hood, shapely uses GEOS to handle geometry objects. GEOS is a C++ library (much faster than Python code), and is one of the fundamental pillars of the open source GIS world, powering geospatial processing for many projects, including QGIS.

Point properties and methods

Points and other shapely geometry objects have useful built-in properties and methods. Using the available attributes, we can for example extract the coordinate values of a Point and calculate the Euclidian distance between points.

The geom_type property contains information about the geometry type of a shapely geometry:

point1.geom_type

'Point'

There are multiple ways to access the coordinates of geometry object. For instance, coords is a spapely.coords.CoordinateSequence. It is an Iterator, an efficient Python data structure to iterate over lists of items, but for now we can simply convert it into a list of the (one pair of) coordinates:

# Get coordinate tuple(s)
list(point1.coords)

[(2.2, 4.2)]

However, since points, by definition, only contain one coordinate tuple, shapely.geometry.Points have properties to directly access its coordinate values: the properties x, y, and (possibly) z, which are basic float type decimal numbers.

# Read x and y coordinates separately
x = point1.x
y = point1.y

print(x, y)

2.2 4.2

# Read x, y, and z
x = point4_3D.x
y = point4_3D.y
z = point4_3D.z

print(x, y, z)

9.26 -2.456 0.57

It is also possible to calculate the distance between two objects using the distance method.

In our example the distance is calculated in a cartesian coordinate system. When working with real GIS data the distance is based on the used coordinate reference system. always check what is the unit of measurement (for example, meters) in the coordinate reference system you are using.

Let’s calculate the distance between point1 and point2:

# Check input data
print(point1)
print(point2)

POINT (2.2 4.2)
POINT (7.2 -25.1)

# Calculate the distance between point1 and point2
dist = point1.distance(point2)

# Print out a nicely formatted info message
print(f"Distance between the points is {dist:.2f} units")

Distance between the points is 29.72 units

Caution:

Shapely geometries are, by design, agnostic (unaware) of the reference system used to represent them. Distances and surface area calculated using the built-in shapely methods will always: a) assume a flat, Cartesian, Euclidean space, and b) return the calculated value in the unit of the coordinates (e.g., meters, or degrees).

This is perfectly fine for small-scale geo-spatial operations, if you keep yourself aware of the expected output unit. Most packages built on top of shapely, for instance GeoPandas, which we will get to know in lesson 2, bring their own functions and take the coordinate reference system into consideration.

Lines

Creating LineString objects is similar to creating Points. Instead of a single coordinate tuple, we pass a list of coordinate tuples, or a list of points, that make up the line:

# import the LineString class
from shapely.geometry import LineString

# Create a LineString from our Point objects
line = LineString([point1, point2, point3])

# Create a LineString from a list of coordinates:
# (with the same coordinate values as the points, so results should be identical)
line2 = LineString([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])

# Check if the lines are, indeed, identical:
line == line2

True

Let’s see how our line looks like:

line

print(line)

LINESTRING (2.2 4.2, 7.2 -25.1, 9.26 -2.456)

Again, the text representation is in WKT format. WKT is convenient as it is a human-readable text format that also most GIS tools can readily use.

It’s not surprising, but we can see that a LineString is constituted of multiple coordinate tuples. In fact, the value(s) of a WKT LINESTRING are made up of the values of multiple WKT POINTS, joined together with a comma.

Check also the data type:

# Check data type of the line object
type(line)

shapely.geometry.linestring.LineString

# Check geometry type of the line object
line.geom_type

'LineString'

LineString properties and methods

Linear geometries in their shapely representations (LineString, LinearRing, MultiLineString) have a variety of properties and methods that expose useful functionality. For instance, it is possible to access a geometry’s coordinates, calculate its lengths, find its centre point, create points along the line at a specified interval, or compute the closest distance between a line an another geometry.

Tip:

Consult the Shapely user manual for a complete list of geometry attributes and operations on one or more geometries.

Fundamentally, accessing the coordinates of a line is very similar as accessing the ones of a point:

# Get coordinate tuples
list(line.coords)

[(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)]

Because a line has to have at least two coordinate tuples, the list now contains more than the one value we saw earlier with points.

If you would need to access all x-coordinates or all y-coordinates of the line, you can use its xy attribute (an iterator, but, again, for now, we can use them as lists):

# Obtain x and y coordinates
xcoords = list(line.xy[0])
ycoords = list(line.xy[1])

print(xcoords)
print(ycoords)

[2.2, 7.2, 9.26]
[4.2, -25.1, -2.456]

Warning: 3D-LineStrings

Note that the xy property of shapely geometries does not return z values for three-dimensional geometries.

Other properties of lines that are useful for GIS analyses include the length and the centre point (centroid) of lines:

# Get the length of the line
line_length = line.length
print(f"Length of our line: {line_length:.1f} units")

Length of our line: 52.5 units

# Get the centre point of the line
print(line.centroid)

POINT (6.229961354035622 -11.892411157572392)

The centroid (or centre point) of a line (or any other shapely geometry) is a shapely.geometry.Point object.

Polygon

Creating a polygon geometry follows the same logic as creating a point or line geometry. However, as discussed above the rules for what constitutes a polygon are more complex: It is constructed of exactly one linear ring forming its exterior (perimeter), and any number of additional linear rings forming holes that are cut out of the exterior shell.

Consequently, the shapely.geometry.Polygon constructor function accepts two parameter: the first one, shell, is a list of coordinate tuples, a list of points, or a LinearRing, and will form the outer hull of the new polygon. The second, optional, parameter holes can be a list of holes to cut out of shell (the items in the list can be the same data types as shell).

For now, let’s create a simple polygon without any holes. The first example uses (at least three) coordinate tuples (three points are required to form a surface):

from shapely.geometry import Polygon

# Create a Polygon from the coordinates
polygon1 = Polygon([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])

We can also construct the polygon directly from a list of points:

polygon2 = Polygon([point1, point2, point3])

… or from a `LinearRing <https://shapely.readthedocs.io/en/stable/manual.html#linearrings>`__ (which has an almost identical behaviour as a LineString, except that it is closed, i.e., the first and last point are identical):

from shapely.geometry import LinearRing

shell = LinearRing([point1, point2, point3, point1])
polygon3 = Polygon(shell)

(When constructing a shapely.geometry.LinearRing, you can omit listing the first point again at the end; shapely will then implicitely add the first point another time to the end of the list of points)

We used different methods to construct the three polygons, but we used the same values. Let’s see whether they ended up describing identical geometries:

polygon1 == polygon2 == polygon3

np.True_

Let’s also see how the polygon looks like drawn, and what its text representation is:

polygon1

print(polygon1)

POLYGON ((2.2 4.2, 7.2 -25.1, 9.26 -2.456, 2.2 4.2))

Just like with points and lines, the text representation of a shapely.geometry.Polygon is in the Well-Known Text format. Note how a WKT POLYGON is made up of the values of one or more WKT LINEARRING’s values (closed line strings), in parentheses, and joined together by commas. The first linear ring represents the exterior, all following ones holes. (Our example polygon consists of one linear ring, only, so no need for the comma).

Check also the data type:

# Data type
type(polygon1)

shapely.geometry.polygon.Polygon

# Geometry type
polygon1.geom_type

'Polygon'

Tip:

You can always use the built-in help() function to find out how a function or class works, which parameters it expects, and what properties and methods you can use.

# Check the help for Polygon objects:
help(Polygon)

Help on class Polygon in module shapely.geometry.polygon:

class Polygon(shapely.geometry.base.BaseGeometry)
 |  Polygon(shell=None, holes=None)
 |
 |  A geometry type representing an area that is enclosed by a linear ring.
 |
 |  A polygon is a two-dimensional feature and has a non-zero area. It may
 |  have one or more negative-space "holes" which are also bounded by linear
 |  rings. If any rings cross each other, the feature is invalid and
 |  operations on it may fail.
 |
 |  Parameters
 |  ----------
 |  shell : sequence
 |      A sequence of (x, y [,z]) numeric coordinate pairs or triples, or
 |      an array-like with shape (N, 2) or (N, 3).
 |      Also can be a sequence of Point objects.
 |  holes : sequence
 |      A sequence of objects which satisfy the same requirements as the
 |      shell parameters above
 |
 |  Attributes
 |  ----------
 |  exterior : LinearRing
 |      The ring which bounds the positive space of the polygon.
 |  interiors : sequence
 |      A sequence of rings which bound all existing holes.
 |
 |  Examples
 |  --------
 |  Create a square polygon with no holes
 |
 |  >>> coords = ((0., 0.), (0., 1.), (1., 1.), (1., 0.), (0., 0.))
 |  >>> polygon = Polygon(coords)
 |  >>> polygon.area
 |  1.0
 |
 |  Method resolution order:
 |      Polygon
 |      shapely.geometry.base.BaseGeometry
 |      shapely.lib.Geometry
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  __eq__(self, other)
 |      Return self==value.
 |
 |  __hash__(self)
 |      Return hash(self).
 |
 |  svg(self, scale_factor=1.0, fill_color=None, opacity=None)
 |      Returns SVG path element for the Polygon geometry.
 |
 |      Parameters
 |      ==========
 |      scale_factor : float
 |          Multiplication factor for the SVG stroke-width.  Default is 1.
 |      fill_color : str, optional
 |          Hex string for fill color. Default is to use "#66cc99" if
 |          geometry is valid, and "#ff3333" if invalid.
 |      opacity : float
 |          Float number between 0 and 1 for color opacity. Default value is 0.6
 |
 |  ----------------------------------------------------------------------
 |  Class methods defined here:
 |
 |  from_bounds(xmin, ymin, xmax, ymax)
 |      Construct a `Polygon()` from spatial bounds.
 |
 |  ----------------------------------------------------------------------
 |  Static methods defined here:
 |
 |  __new__(self, shell=None, holes=None)
 |      Create and return a new object.  See help(type) for accurate signature.
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties defined here:
 |
 |  __geo_interface__
 |      Dictionary representation of the geometry
 |
 |  coords
 |      Access to geometry's coordinates (CoordinateSequence)
 |
 |  exterior
 |
 |  interiors
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from shapely.geometry.base.BaseGeometry:
 |
 |  __and__(self, other)
 |
 |  __bool__(self)
 |
 |  __format__(self, format_spec)
 |      Format a geometry using a format specification.
 |
 |  __ne__(self, other)
 |      Return self!=value.
 |
 |  __nonzero__(self)
 |
 |  __or__(self, other)
 |      Return self|value.
 |
 |  __reduce__(self)
 |      Helper for pickle.
 |
 |  __repr__(self)
 |      Return repr(self).
 |
 |  __str__(self)
 |      Return str(self).
 |
 |  __sub__(self, other)
 |
 |  __xor__(self, other)
 |
 |  almost_equals(self, other, decimal=6)
 |      True if geometries are equal at all coordinates to a
 |      specified decimal place.
 |
 |      .. deprecated:: 1.8.0
 |          The 'almost_equals()' method is deprecated
 |          and will be removed in Shapely 2.1 because the name is
 |          confusing. The 'equals_exact()' method should be used
 |          instead.
 |
 |      Refers to approximate coordinate equality, which requires
 |      coordinates to be approximately equal and in the same order for
 |      all components of a geometry.
 |
 |      Because of this it is possible for "equals()" to be True for two
 |      geometries and "almost_equals()" to be False.
 |
 |      Examples
 |      --------
 |      >>> LineString(
 |      ...     [(0, 0), (2, 2)]
 |      ... ).equals_exact(
 |      ...     LineString([(0, 0), (1, 1), (2, 2)]),
 |      ...     1e-6
 |      ... )
 |      False
 |
 |      Returns
 |      -------
 |      bool
 |
 |  buffer(self, distance, quad_segs=16, cap_style='round', join_style='round', mitre_limit=5.0, single_sided=False, **kwargs)
 |      Get a geometry that represents all points within a distance
 |      of this geometry.
 |
 |      A positive distance produces a dilation, a negative distance an
 |      erosion. A very small or zero distance may sometimes be used to
 |      "tidy" a polygon.
 |
 |      Parameters
 |      ----------
 |      distance : float
 |          The distance to buffer around the object.
 |      resolution : int, optional
 |          The resolution of the buffer around each vertex of the
 |          object.
 |      quad_segs : int, optional
 |          Sets the number of line segments used to approximate an
 |          angle fillet.
 |      cap_style : shapely.BufferCapStyle or {'round', 'square', 'flat'}, default 'round'
 |          Specifies the shape of buffered line endings. BufferCapStyle.round ('round')
 |          results in circular line endings (see ``quad_segs``). Both BufferCapStyle.square
 |          ('square') and BufferCapStyle.flat ('flat') result in rectangular line endings,
 |          only BufferCapStyle.flat ('flat') will end at the original vertex,
 |          while BufferCapStyle.square ('square') involves adding the buffer width.
 |      join_style : shapely.BufferJoinStyle or {'round', 'mitre', 'bevel'}, default 'round'
 |          Specifies the shape of buffered line midpoints. BufferJoinStyle.ROUND ('round')
 |          results in rounded shapes. BufferJoinStyle.bevel ('bevel') results in a beveled
 |          edge that touches the original vertex. BufferJoinStyle.mitre ('mitre') results
 |          in a single vertex that is beveled depending on the ``mitre_limit`` parameter.
 |      mitre_limit : float, optional
 |          The mitre limit ratio is used for very sharp corners. The
 |          mitre ratio is the ratio of the distance from the corner to
 |          the end of the mitred offset corner. When two line segments
 |          meet at a sharp angle, a miter join will extend the original
 |          geometry. To prevent unreasonable geometry, the mitre limit
 |          allows controlling the maximum length of the join corner.
 |          Corners with a ratio which exceed the limit will be beveled.
 |      single_side : bool, optional
 |          The side used is determined by the sign of the buffer
 |          distance:
 |
 |              a positive distance indicates the left-hand side
 |              a negative distance indicates the right-hand side
 |
 |          The single-sided buffer of point geometries is the same as
 |          the regular buffer.  The End Cap Style for single-sided
 |          buffers is always ignored, and forced to the equivalent of
 |          CAP_FLAT.
 |      quadsegs : int, optional
 |          Deprecated alias for `quad_segs`.
 |
 |      Returns
 |      -------
 |      Geometry
 |
 |      Notes
 |      -----
 |      The return value is a strictly two-dimensional geometry. All
 |      Z coordinates of the original geometry will be ignored.
 |
 |      Examples
 |      --------
 |      >>> from shapely.wkt import loads
 |      >>> g = loads('POINT (0.0 0.0)')
 |
 |      16-gon approx of a unit radius circle:
 |
 |      >>> g.buffer(1.0).area  # doctest: +ELLIPSIS
 |      3.1365484905459...
 |
 |      128-gon approximation:
 |
 |      >>> g.buffer(1.0, 128).area  # doctest: +ELLIPSIS
 |      3.141513801144...
 |
 |      triangle approximation:
 |
 |      >>> g.buffer(1.0, 3).area
 |      3.0
 |      >>> list(g.buffer(1.0, cap_style=BufferCapStyle.square).exterior.coords)
 |      [(1.0, 1.0), (1.0, -1.0), (-1.0, -1.0), (-1.0, 1.0), (1.0, 1.0)]
 |      >>> g.buffer(1.0, cap_style=BufferCapStyle.square).area
 |      4.0
 |
 |  contains(self, other)
 |      Returns True if the geometry contains the other, else False
 |
 |  contains_properly(self, other)
 |      Returns True if the geometry completely contains the other, with no
 |      common boundary points, else False
 |
 |      Refer to `shapely.contains_properly` for full documentation.
 |
 |  covered_by(self, other)
 |      Returns True if the geometry is covered by the other, else False
 |
 |  covers(self, other)
 |      Returns True if the geometry covers the other, else False
 |
 |  crosses(self, other)
 |      Returns True if the geometries cross, else False
 |
 |  difference(self, other, grid_size=None)
 |      Returns the difference of the geometries.
 |
 |      Refer to `shapely.difference` for full documentation.
 |
 |  disjoint(self, other)
 |      Returns True if geometries are disjoint, else False
 |
 |  distance(self, other)
 |      Unitless distance to other geometry (float)
 |
 |  dwithin(self, other, distance)
 |      Returns True if geometry is within a given distance from the other, else False.
 |
 |      Refer to `shapely.dwithin` for full documentation.
 |
 |  equals(self, other)
 |      Returns True if geometries are equal, else False.
 |
 |      This method considers point-set equality (or topological
 |      equality), and is equivalent to (self.within(other) &
 |      self.contains(other)).
 |
 |      Examples
 |      --------
 |      >>> LineString(
 |      ...     [(0, 0), (2, 2)]
 |      ... ).equals(
 |      ...     LineString([(0, 0), (1, 1), (2, 2)])
 |      ... )
 |      True
 |
 |      Returns
 |      -------
 |      bool
 |
 |  equals_exact(self, other, tolerance)
 |      True if geometries are equal to within a specified
 |      tolerance.
 |
 |      Parameters
 |      ----------
 |      other : BaseGeometry
 |          The other geometry object in this comparison.
 |      tolerance : float
 |          Absolute tolerance in the same units as coordinates.
 |
 |      This method considers coordinate equality, which requires
 |      coordinates to be equal and in the same order for all components
 |      of a geometry.
 |
 |      Because of this it is possible for "equals()" to be True for two
 |      geometries and "equals_exact()" to be False.
 |
 |      Examples
 |      --------
 |      >>> LineString(
 |      ...     [(0, 0), (2, 2)]
 |      ... ).equals_exact(
 |      ...     LineString([(0, 0), (1, 1), (2, 2)]),
 |      ...     1e-6
 |      ... )
 |      False
 |
 |      Returns
 |      -------
 |      bool
 |
 |  geometryType(self)
 |
 |  hausdorff_distance(self, other)
 |      Unitless hausdorff distance to other geometry (float)
 |
 |  interpolate(self, distance, normalized=False)
 |      Return a point at the specified distance along a linear geometry
 |
 |      Negative length values are taken as measured in the reverse
 |      direction from the end of the geometry. Out-of-range index
 |      values are handled by clamping them to the valid range of values.
 |      If the normalized arg is True, the distance will be interpreted as a
 |      fraction of the geometry's length.
 |
 |      Alias of `line_interpolate_point`.
 |
 |  intersection(self, other, grid_size=None)
 |      Returns the intersection of the geometries.
 |
 |      Refer to `shapely.intersection` for full documentation.
 |
 |  intersects(self, other)
 |      Returns True if geometries intersect, else False
 |
 |  line_interpolate_point(self, distance, normalized=False)
 |      Return a point at the specified distance along a linear geometry
 |
 |      Negative length values are taken as measured in the reverse
 |      direction from the end of the geometry. Out-of-range index
 |      values are handled by clamping them to the valid range of values.
 |      If the normalized arg is True, the distance will be interpreted as a
 |      fraction of the geometry's length.
 |
 |      Alias of `interpolate`.
 |
 |  line_locate_point(self, other, normalized=False)
 |      Returns the distance along this geometry to a point nearest the
 |      specified point
 |
 |      If the normalized arg is True, return the distance normalized to the
 |      length of the linear geometry.
 |
 |      Alias of `project`.
 |
 |  normalize(self)
 |      Converts geometry to normal form (or canonical form).
 |
 |      This method orders the coordinates, rings of a polygon and parts of
 |      multi geometries consistently. Typically useful for testing purposes
 |      (for example in combination with `equals_exact`).
 |
 |      Examples
 |      --------
 |      >>> from shapely import MultiLineString
 |      >>> line = MultiLineString([[(0, 0), (1, 1)], [(3, 3), (2, 2)]])
 |      >>> line.normalize()
 |      <MULTILINESTRING ((2 2, 3 3), (0 0, 1 1))>
 |
 |  overlaps(self, other)
 |      Returns True if geometries overlap, else False
 |
 |  point_on_surface(self)
 |      Returns a point guaranteed to be within the object, cheaply.
 |
 |      Alias of `representative_point`.
 |
 |  project(self, other, normalized=False)
 |      Returns the distance along this geometry to a point nearest the
 |      specified point
 |
 |      If the normalized arg is True, return the distance normalized to the
 |      length of the linear geometry.
 |
 |      Alias of `line_locate_point`.
 |
 |  relate(self, other)
 |      Returns the DE-9IM intersection matrix for the two geometries
 |      (string)
 |
 |  relate_pattern(self, other, pattern)
 |      Returns True if the DE-9IM string code for the relationship between
 |      the geometries satisfies the pattern, else False
 |
 |  representative_point(self)
 |      Returns a point guaranteed to be within the object, cheaply.
 |
 |      Alias of `point_on_surface`.
 |
 |  reverse(self)
 |      Returns a copy of this geometry with the order of coordinates reversed.
 |
 |      If the geometry is a polygon with interior rings, the interior rings are also
 |      reversed.
 |
 |      Points are unchanged.
 |
 |      See also
 |      --------
 |      is_ccw : Checks if a geometry is clockwise.
 |
 |      Examples
 |      --------
 |      >>> from shapely import LineString, Polygon
 |      >>> LineString([(0, 0), (1, 2)]).reverse()
 |      <LINESTRING (1 2, 0 0)>
 |      >>> Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]).reverse()
 |      <POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))>
 |
 |  segmentize(self, max_segment_length)
 |      Adds vertices to line segments based on maximum segment length.
 |
 |      Additional vertices will be added to every line segment in an input geometry
 |      so that segments are no longer than the provided maximum segment length. New
 |      vertices will evenly subdivide each segment.
 |
 |      Only linear components of input geometries are densified; other geometries
 |      are returned unmodified.
 |
 |      Parameters
 |      ----------
 |      max_segment_length : float or array_like
 |          Additional vertices will be added so that all line segments are no
 |          longer this value.  Must be greater than 0.
 |
 |      Examples
 |      --------
 |      >>> from shapely import LineString, Polygon
 |      >>> LineString([(0, 0), (0, 10)]).segmentize(max_segment_length=5)
 |      <LINESTRING (0 0, 0 5, 0 10)>
 |      >>> Polygon([(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)]).segmentize(max_segment_length=5)
 |      <POLYGON ((0 0, 5 0, 10 0, 10 5, 10 10, 5 10, 0 10, 0 5, 0 0))>
 |
 |  simplify(self, tolerance, preserve_topology=True)
 |      Returns a simplified geometry produced by the Douglas-Peucker
 |      algorithm
 |
 |      Coordinates of the simplified geometry will be no more than the
 |      tolerance distance from the original. Unless the topology preserving
 |      option is used, the algorithm may produce self-intersecting or
 |      otherwise invalid geometries.
 |
 |  symmetric_difference(self, other, grid_size=None)
 |      Returns the symmetric difference of the geometries.
 |
 |      Refer to `shapely.symmetric_difference` for full documentation.
 |
 |  touches(self, other)
 |      Returns True if geometries touch, else False
 |
 |  union(self, other, grid_size=None)
 |      Returns the union of the geometries.
 |
 |      Refer to `shapely.union` for full documentation.
 |
 |  within(self, other)
 |      Returns True if geometry is within the other, else False
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from shapely.geometry.base.BaseGeometry:
 |
 |  area
 |      Unitless area of the geometry (float)
 |
 |  boundary
 |      Returns a lower dimension geometry that bounds the object
 |
 |      The boundary of a polygon is a line, the boundary of a line is a
 |      collection of points. The boundary of a point is an empty (null)
 |      collection.
 |
 |  bounds
 |      Returns minimum bounding region (minx, miny, maxx, maxy)
 |
 |  centroid
 |      Returns the geometric center of the object
 |
 |  convex_hull
 |      Imagine an elastic band stretched around the geometry: that's a
 |      convex hull, more or less
 |
 |      The convex hull of a three member multipoint, for example, is a
 |      triangular polygon.
 |
 |  envelope
 |      A figure that envelopes the geometry
 |
 |  geom_type
 |      Name of the geometry's type, such as 'Point'
 |
 |  has_z
 |      True if the geometry's coordinate sequence(s) have z values (are
 |      3-dimensional)
 |
 |  is_closed
 |      True if the geometry is closed, else False
 |
 |      Applicable only to 1-D geometries.
 |
 |  is_empty
 |      True if the set of points in this geometry is empty, else False
 |
 |  is_ring
 |      True if the geometry is a closed ring, else False
 |
 |  is_simple
 |      True if the geometry is simple, meaning that any self-intersections
 |      are only at boundary points, else False
 |
 |  is_valid
 |      True if the geometry is valid (definition depends on sub-class),
 |      else False
 |
 |  length
 |      Unitless length of the geometry (float)
 |
 |  minimum_clearance
 |      Unitless distance by which a node could be moved to produce an invalid geometry (float)
 |
 |  minimum_rotated_rectangle
 |      Returns the oriented envelope (minimum rotated rectangle) that
 |      encloses the geometry.
 |
 |      Unlike `envelope` this rectangle is not constrained to be parallel to the
 |      coordinate axes. If the convex hull of the object is a degenerate (line
 |      or point) this degenerate is returned.
 |
 |      Alias of `oriented_envelope`.
 |
 |  oriented_envelope
 |      Returns the oriented envelope (minimum rotated rectangle) that
 |      encloses the geometry.
 |
 |      Unlike envelope this rectangle is not constrained to be parallel to the
 |      coordinate axes. If the convex hull of the object is a degenerate (line
 |      or point) this degenerate is returned.
 |
 |      Alias of `minimum_rotated_rectangle`.
 |
 |  type
 |
 |  wkb
 |      WKB representation of the geometry
 |
 |  wkb_hex
 |      WKB hex representation of the geometry
 |
 |  wkt
 |      WKT representation of the geometry
 |
 |  xy
 |      Separate arrays of X and Y coordinate values
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from shapely.lib.Geometry:
 |
 |  __ge__(self, value, /)
 |      Return self>=value.
 |
 |  __gt__(self, value, /)
 |      Return self>value.
 |
 |  __le__(self, value, /)
 |      Return self<=value.
 |
 |  __lt__(self, value, /)
 |      Return self<value.
 |
 |  __setstate__(...)
 |      For unpickling pre-shapely 2.0 pickles

Let’s still see how to create a polygon with a hole:

# define the exterior
outer = LinearRing([(-180, 90), (-180, -90), (180, -90), (180, 90)])

# define a hole:
hole = LinearRing([(-170, 80), (-100, -80), (100, -80), (170, 80)])

Let’s see how the exterior shell and the hole look like on their own:

outer

hole

A polygon using only the exterior shell:

polygon_without_hole = Polygon(outer)
polygon_without_hole

And, finally, a polygon defined by the exterior shell, and one hole (note that holes need to be specified as a list):

polygon_with_hole = Polygon(outer, [hole])
polygon_with_hole

print(polygon_without_hole)
print(polygon_with_hole)

POLYGON ((-180 90, -180 -90, 180 -90, 180 90, -180 90))
POLYGON ((-180 90, -180 -90, 180 -90, 180 90, -180 90), (-170 80, -100 -80, 100 -80, 170 80, -170 80))

Polygon properties and methods

Very similar to lines and points, also shapely.geometry.Polygons expose a number of properties and methods that can be useful for spatial analysis tasks. Consult the shapely user manual for a complete list, and see a few examples here:

print(f"Polygon centroid: {polygon_with_hole.centroid}")
print(f"Polygon area: {polygon_with_hole.area}")
print(f"Polygon bounding box: {polygon_with_hole.bounds}")
print(f"Polygon exterior ring: {polygon_with_hole.exterior}")
print(f"Polygon circumference: {polygon_with_hole.exterior.length}")

Polygon centroid: POINT (0 -13.827160493827162)
Polygon area: 21600.0
Polygon bounding box: (-180.0, -90.0, 180.0, 90.0)
Polygon exterior ring: LINEARRING (-180 90, -180 -90, 180 -90, 180 90, -180 90)
Polygon circumference: 1080.0

As we can see above, it is again fairly straightforward to access different attributes of Polygon objects. Note that distance metrics will make more sense when we start working with data in projected coordinate systems.

Check your understanding

Plot these shapes using shapely!

• Pentagon, example coordinates: (30, 2.01), (31.91, 0.62), (31.18, -1.63), (28.82, -1.63), (28.09, 0.62)

• Triangle

• Square

• Circle

# Pentagon
Polygon([(30, 2.01), (31.91, 0.62), (31.18, -1.63), (28.82, -1.63), (28.09, 0.62)])

# Triangle
Polygon([(0,0), (2,4), (4,0)])

# Square
Polygon([(0,0), (0,4), (4,4), (4,0)])

# Circle (using a buffer around a point)
point = Point((0,0))
point.buffer(1)

Geometry collections (optional)

Sometimes, it can be useful to store multiple geometries (for example, several points or several polygons) in a single feature. See, for instance, the example of Tyrol above. Its two parts share the same attributes, and together are equivalent to any of the other provinces of Austria. Semantically, it would not make sense to store the two polygons in separate rows. If expressed as a MultiPolygon, Tyrol would represent one row of information in the attribute table with multiple polygons attached.

Caution:

By convention, data sets should always consist of either single- or multi-geometries. Some file formats enforce this, and many GIS tools refuse operation on data sets with mixed single- and multi-geometries.

If one feature in a data set is a MultiGeometry, all other features should be converted, too. All single geometries can be expressed as a collection of one item.

In shapely, collections of points are implemented as MultiPoint geometries, collections of lines as MultiLineString geometries, and collections of polygons as MultiPolygon geometries.

from shapely.geometry import MultiPoint, MultiLineString, MultiPolygon

# Create a MultiPoint object of our points 1,2 and 3
multipoint = MultiPoint([point1, point2, point3])

# We can also create a MultiLineString with two lines
line1 = LineString([point1, point2])
line2 = LineString([point2, point3])
multiline = MultiLineString([line1, line2])

print(multipoint)
print(multiline)

MULTIPOINT (2.2 4.2, 7.2 -25.1, 9.26 -2.456)
MULTILINESTRING ((2.2 4.2, 7.2 -25.1), (7.2 -25.1, 9.26 -2.456))

multipoint

multiline

MultiPolygons are constructed in a similar manner. Let’s create a bounding box for ‘the world’ by combining two separate polygons that represent the western and eastern hemispheres.

# Let’s create the exterior of the western part of the world
western_hemisphere = Polygon([(-180, 90), (-180, -90), (0, -90), (0, 90)])
print(western_hemisphere)
western_hemisphere

POLYGON ((-180 90, -180 -90, 0 -90, 0 90, -180 90))

Shapely has a short-hand function for creating rectangular polygons, ‘boxes’. It can be used, for instance, to create bounding boxes using minimum and maximum x and y coordinates.

Let’s use shapely.geometry.box() for creating the polygon representing the the eastern hemisphere:

from shapely.geometry import box
min_x = 0
max_x = 180
min_y = -90
max_y = 90

eastern_hemisphere = box(min_x, min_y, max_x, max_y)

print(eastern_hemisphere)
eastern_hemisphere

POLYGON ((180 -90, 180 90, 0 90, 0 -90, 180 -90))

Finally, we can combine the two polygons into a MultiPolygon:

# Let’s create our MultiPolygon.
# Pass multiple Polygon objects as a list
multipolygon = MultiPolygon([western_hemisphere, eastern_hemisphere])

print(multipolygon)
multipolygon

MULTIPOLYGON (((-180 90, -180 -90, 0 -90, 0 90, -180 90)), ((180 -90, 180 90, 0 90, 0 -90, 180 -90)))

Multi-geometries are in many ways similar to simple geometries, like the ones we created earlier. The main difference is that they can combine multiple geometric primitives into one feature.

Convex hull and envelope

A ‘convex hull’ refers to the smallest possible convex polygon that can contain a geometry or a set of geometries. Alongside bounding boxes, convex hulls are useful to describe the extent of data sets.

# Check input geometry
multipoint

# Convex Hull
multipoint.convex_hull

# Envelope (smalles rectangular polygon around a geometry/set of geometries):
multipoint.envelope

Validity of geometries

As discussed on the very top of this page, already the geometric primitives have certain requirements. For instance, a LineString must consist of at least two points, and a Polygon’s exterior shell and holes must not intersect.

Each shapely geometry has a built-in check that can be of great help, for instance, finding topological errors:

print(f"Is polygon valid?: {polygon_with_hole.is_valid}")

Is polygon valid?: True