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#
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 carruing 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).
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:
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]
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
(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 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)
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 two-dimensional figure bounded by a linear ring
|
| A polygon 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.
|
| 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 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) 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 self|value.
|
| __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)]),
| ... 1e-6
| ... )
| 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 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.
|
| 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)')
| >>> g.buffer(1.0).area # 16-gon approx of a unit radius circle
| 3.1365484905459...
| >>> g.buffer(1.0, 128).area # 128-gon 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=94359748446880)
|
| 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.
|
| 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 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.
|
| 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)
| 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
| 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 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__ = 94359748446880
|
| __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.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, multipoint]
[<shapely.geometry.polygon.Polygon at 0x7f9358d8eb00>,
<shapely.geometry.multipoint.MultiPoint at 0x7f937822f4c0>]
# 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
shapely 2.0
While we are having this course, the team developing shapely is preparing the library’s next updates. It will be a major version that breaks with some of the programming patterns that were possible with earlier versions.
When you work with shapely in the future, be sure to check out what will have changed.