Shapely and geometry objects
Contents
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#
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 geospatial objects, and for carrying out a range of geometric operations. A basic understanding of how shapely works is paramount for using higherlevel 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 twodimensional (x, y)
, or threedimensional (x, y, z)
.
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 onevalue 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 threedimensional 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:
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 multipoints, multilines, and multipolygons.
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 ‘WellKnown 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 threedimensional 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 builtin 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.Point
s 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 builtin 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 smallscale geospatial 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 humanreadable 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 xcoordinates or all ycoordinates 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]
3DLineStrings
Note that the xy
property of shapely geometries does not return z
values
for threedimensional 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
(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
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 WellKnown 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 builtin 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)
Show code cell output
Help on class Polygon in module shapely.geometry.polygon:
class Polygon(shapely.geometry.base.BaseGeometry)
 Polygon(shell=None, holes=None)

 A twodimensional figure bounded by a linear ring

 A polygon has a nonzero area. It may have one or more negativespace
 "holes" which are also bounded by linear rings. If any rings cross each
 other, the feature is invalid and operations on it may fail.

 Attributes
 
 exterior : LinearRing
 The ring which bounds the positive space of the polygon.
 interiors : sequence
 A sequence of rings which bound all existing holes.

 Method resolution order:
 Polygon
 shapely.geometry.base.BaseGeometry
 builtins.object

 Methods defined here:

 __eq__(self, other)
 Return self==value.

 __init__(self, shell=None, holes=None)
 Parameters
 
 shell : sequence
 A sequence of (x, y [,z]) numeric coordinate pairs or triples.
 Also can be a sequence of Point objects.
 holes : sequence
 A sequence of objects which satisfy the same requirements as the
 shell parameters above

 Example
 
 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

 __ne__(self, other)
 Return self!=value.

 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 strokewidth. 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) from builtins.type
 Construct a `Polygon()` from spatial bounds.

 
 Readonly properties defined here:

 __array_interface__
 Provide the Numpy array protocol.

 __geo_interface__
 Dictionary representation of the geometry

 coords
 Access to geometry's coordinates (CoordinateSequence)

 exterior

 interiors

 
 Data and other attributes defined here:

 __hash__ = None

 
 Methods inherited from shapely.geometry.base.BaseGeometry:

 __and__(self, other)

 __bool__(self)

 __del__(self)

 __nonzero__(self)

 __or__(self, other)
 Return selfvalue.

 __reduce__(self)
 Helper for pickle.

 __setattr__(self, name, value)
 Implement setattr(self, name, value).

 __setstate__(self, state)

 __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.0 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)]),
 ... 1e6
 ... )
 False

 Returns
 
 bool

 buffer(self, distance, resolution=16, quadsegs=None, cap_style=1, join_style=1, mitre_limit=5.0, single_sided=False)
 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.
 quadsegs : int, optional
 Sets the number of line segments used to approximate an
 angle fillet. Note: the use of a `quadsegs` parameter is
 deprecated and will be gone from the next major release.
 cap_style : int, optional
 The styles of caps are: CAP_STYLE.round (1), CAP_STYLE.flat
 (2), and CAP_STYLE.square (3).
 join_style : int, optional
 The styles of joins between offset segments are:
 JOIN_STYLE.round (1), JOIN_STYLE.mitre (2), and
 JOIN_STYLE.bevel (3).
 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 lefthand side
 a negative distance indicates the righthand side

 The singlesided buffer of point geometries is the same as
 the regular buffer. The End Cap Style for singlesided
 buffers is always ignored, and forced to the equivalent of
 CAP_FLAT.

 Returns
 
 Geometry

 Notes
 
 The return value is a strictly twodimensional geometry. All
 Z coordinates of the original geometry will be ignored.

 Examples
 
 >>> from shapely.wkt import loads
 >>> g = loads('POINT (0.0 0.0)')
 >>> g.buffer(1.0).area # 16gon approx of a unit radius circle
 3.1365484905459...
 >>> g.buffer(1.0, 128).area # 128gon approximation
 3.141513801144...
 >>> round(g.buffer(1.0, 3).area, 10) # triangle approximation
 3.0
 >>> list(g.buffer(1.0, cap_style=CAP_STYLE.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=CAP_STYLE.square).area
 4.0

 contains(self, other)
 Returns True if the geometry contains the other, else False

 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)
 Returns the difference of the geometries

 disjoint(self, other)
 Returns True if geometries are disjoint, else False

 distance(self, other)
 Unitless distance to other geometry (float)

 empty(self, val=94823010514608)

 equals(self, other)
 Returns True if geometries are equal, else False.

 This method considers pointset 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)]),
 ... 1e6
 ... )
 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. Outofrange 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.

 intersection(self, other)
 Returns the intersection of the geometries

 intersects(self, other)
 Returns True if geometries intersect, else False

 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.wkt import loads
 >>> p = loads("MULTILINESTRING((0 0, 1 1), (3 3, 2 2))")
 >>> p.normalize().wkt
 'MULTILINESTRING ((2 2, 3 3), (0 0, 1 1))'

 overlaps(self, other)
 Returns True if geometries overlap, else False

 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.

 relate(self, other)
 Returns the DE9IM intersection matrix for the two geometries
 (string)

 relate_pattern(self, other, pattern)
 Returns True if the DE9IM 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.

 simplify(self, tolerance, preserve_topology=True)
 Returns a simplified geometry produced by the DouglasPeucker
 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 selfintersecting or
 otherwise invalid geometries.

 symmetric_difference(self, other)
 Returns the symmetric difference of the geometries
 (Shapely geometry)

 touches(self, other)
 Returns True if geometries touch, else False

 union(self, other)
 Returns the union of the geometries (Shapely geometry)

 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)

 array_interface_base

 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.

 ctypes
 Return ctypes buffer

 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
 3dimensional)

 is_closed
 True if the geometry is closed, else False

 Applicable only to 1D 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 selfintersections
 are only at boundary points, else False

 is_valid
 True if the geometry is valid (definition depends on subclass),
 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 general minimum bounding rectangle of
 the geometry. Can possibly be rotated. If the convex hull
 of the object is a degenerate (line or point) this same degenerate
 is returned.

 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

 
 Data descriptors inherited from shapely.geometry.base.BaseGeometry:

 __dict__
 dictionary for instance variables (if defined)

 __weakref__
 list of weak references to the object (if defined)

 
 Data and other attributes inherited from shapely.geometry.base.BaseGeometry:

 __geom__ = 94823010514608

 __p__ = None

 impl = <GEOSImpl object: GEOS C API version (1, 13, 0)>
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.Polygon
s 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 multigeometries. Some file formats enforce this, and many GIS tools refuse operation on data sets with mixed single and multigeometries.
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 shorthand 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)))
Multigeometries 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 builtin 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