Deterministic spatial analysis is an important component of computational
approaches to problems in agriculture, ecology, epidemiology, sociology, and
many other fields. What is the surveyed perimeter/area ratio of these patches
of animal habitat? Which properties in this town intersect with the 50-year
flood contour from this new flooding model? What are the extents of findspots
for ancient ceramic wares with maker’s marks “A” and “B”, and where do the
extents overlap? What’s the path from home to office that best skirts
identified zones of location based spam? These are just a few of the possible
questions addressable using non-statistical spatial analysis, and more
specifically, computational geometry.

Shapely is a Python package for set-theoretic analysis and manipulation of
planar features using (via Python’s ctypes module) functions from the
well known and widely deployed GEOS library. GEOS, a port of the Java
Topology Suite (JTS), is the geometry engine of the PostGIS spatial
extension for the PostgreSQL RDBMS. The designs of JTS and GEOS are largely
guided by the Open Geospatial Consortium‘s Simple Features Access
Specification [1] and Shapely adheres mainly to the same set of standard
classes and operations. Shapely is thereby deeply rooted in the conventions of
the geographic information systems (GIS) world, but aspires to be equally
useful to programmers working on non-conventional problems.

The first premise of Shapely is that Python programmers should be able to
perform PostGIS type geometry operations outside of an RDBMS. Not all
geographic data originate or reside in a RDBMS or are best processed using SQL.
We can load data into a spatial RDBMS to do work, but if there’s no mandate to
manage (the “M” in “RDBMS”) the data over time in the database we’re using the
wrong tool for the job. The second premise is that the persistence,
serialization, and map projection of features are significant, but orthogonal
problems. You may not need a hundred GIS format readers and writers or the
multitude of State Plane projections, and Shapely doesn’t burden you with them.
The third premise is that Python idioms trump GIS (or Java, in this case, since
the GEOS library is derived from JTS, a Java project) idioms.

If you enjoy and profit from idiomatic Python, appreciate packages that do one
thing well, and agree that a spatially enabled RDBMS is often enough the wrong
tool for your computational geometry job, Shapely might be for you.

The fundamental types of geometric objects implemented by Shapely are points,
curves, and surfaces. Each is associated with three sets of (possibly infinite)
points in the plane. The interior, boundary, and exterior sets of a
feature are mutually exclusive and their union coincides with the entire plane
[2].

A Point has an interior set of exactly one point, a boundary set of
exactly no points, and an exterior set of all other points. A Point has
a topological dimension of 0.

A Curve has an interior set consisting of the infinitely many points
along its length (imagine a Point dragged in space), a boundary set
consisting of its two end points, and an exterior set of all other points.
A Curve has a topological dimension of 1.

A Surface has an interior set consisting of the infinitely many points
within (imagine a Curve dragged in space to cover an area), a boundary
set consisting of one or more Curves, and an exterior set of all other
points including those within holes that might exist in the surface. A
Surface has a topological dimension of 2.

That may seem a bit esoteric, but will help clarify the meanings of Shapely’s
spatial predicates, and it’s as deep into theory as this manual will go.
Consequences of point-set theory, including some that manifest themselves as
“gotchas”, for different classes will be discussed later in this manual.

The point type is implemented by a Point class; curve by the LineString and
LinearRing classes; and surface by a Polygon class. Shapely implements no
smooth (i.e. having continuous tangents) curves. All curves must be
approximated by linear splines. All rounded patches must be approximated by
regions bounded by linear splines.

Collections of points are implemented by a MultiPoint class, collections of
curves by a MultiLineString class, and collections of surfaces by a
MultiPolygon class. These collections aren’t computationally significant, but
are useful for modeling certain kinds of features. A Y-shaped line feature, for
example, is well modeled as a whole by a MultiLineString.

The standard data model has additional constraints specific to certain types
of geometric objects that will be discussed in following sections of this
manual.

The spatial data model is accompanied by a group of natural language
relationships between geometric objects – contains, intersects, overlaps,
touches, etc. – and a theoretical framework for understanding them using the
3x3 matrix of the mutual intersections of their component point sets [2]: the
DE-9IM. A comprehensive review of the relationships in terms of the DE-9IM is
found in [4] and will not be reiterated in this manual.

Following the JTS technical specs [5], this manual will make a distinction
between constructive (buffer, convex hull) and set-theoretic operations
(intersection, union, etc.). The individual operations will be fully
described in a following section of the manual.

Even though the Earth is not flat – and for that matter not exactly spherical –
there are many analytic problems that can be approached by transforming Earth
features to a Cartesian plane, applying tried and true algorithms, and then
transforming the results back to geographic coordinates. This practice is as
old as the tradition of accurate paper maps.

Shapely does not support coordinate system transformations. All operations on
two or more features presume that the features exist in the same Cartesian
plane.

Geometric objects are created in the typical Python fashion, using the classes
themselves as instance factories. A few of their intrinsic properties will be
discussed in this sections, others in the following sections on operations and
serializations.

Instances of Point, LineString, and LinearRing have as their most
important attribute a finite sequence of coordinates that determines their
interior, boundary, and exterior point sets. A line string can be determined by
as few as 2 points, but contains an infinite number of points. Coordinate
sequences are immutable. A third z coordinate value may be used when
constructing instances, but has no effect on geometric analysis. All
operations are performed in the x-y plane.

In all constructors, numeric values are converted to type float. In other
words, Point(0,0) and Point(0.0,0.0) produce geometrically equivalent
instances. Shapely does not check the topological simplicity or validity of
instances when they are constructed as the cost is unwarranted in most cases.
Validating factories are trivially implemented, using the is_valid
predicate, by users that require them.

The LineString constructor takes an ordered sequence of 2 or more
(x,y[,z]) point tuples.

The constructed LineString object represents one or more connected linear
splines between the points. Repeated points in the ordered sequence are
allowed, but may incur performance penalties and should be avoided. A
LineString may cross itself (i.e. be complex and not simple).

The LinearRing constructor takes an ordered sequence of (x,y[,z])
point tuples.

The sequence may be explicitly closed by passing identical values in the first
and last indices. Otherwise, the sequence will be implicitly closed by copying
the first tuple to the last index. As with a LineString, repeated points in
the ordered sequence are allowed, but may incur performance penalties and
should be avoided. A LinearRing may not cross itself, and may not touch
itself at a single point.

The Polygon constructor takes two positional parameters. The first is an
ordered sequence of (x,y[,z]) point tuples and is treated exactly as in
the LinearRing case. The second is an optional unordered sequence of
ring-like sequences specifying the interior boundaries or “holes” of the
feature.

Rings of a validPolygon may not cross each other, but may touch at a
single point only. Again, Shapely will not prevent the creation of invalid
features, but exceptions will be raised when they are operated on.

Figure 3. On the left, a valid Polygon with one interior ring that touches
the exterior ring at one point, and on the right a Polygon that is invalid
because its interior ring touches the exterior ring at more than one point. The
points that describe the rings are shown in grey.

Returns a properly oriented copy of the given polygon. The signed area of the
result will have the given sign. A sign of 1.0 means that the coordinates of
the product’s exterior ring will be oriented counter-clockwise.

Heterogeneous collections of geometric objects may result from some Shapely
operations. For example, two LineStrings may intersect along a line and at a
point. To represent these kind of results, Shapely provides frozenset-like,
immutable collections of geometric objects. The collections may be homogeneous
(MultiPoint etc.) or heterogeneous.

Figure 6. On the left, a simple, disconnected MultiLineString, and on the
right, a non-simple MultiLineString. The points defining the objects are
shown in gray, the boundaries of the objects in black.

An “empty” feature is one with a point set that coincides with the empty set;
not None, but like set([]). Empty features can be created by calling
the various constructors with no arguments. Almost no operations are supported
by empty features.

It can be useful to specify position along linear features such as LineStrings
and MultiLineStrings with a 1-dimensional referencing system. Shapely
supports linear referencing based on length or distance, evaluating the
distance along a geometric object to the projection of a given point, or the
point at a given distance along the object.

For example, the linear referencing methods might be used to cut lines at a
specified distance.

defcut(line,distance):# Cuts a line in two at a distance from its starting pointifdistance<=0.0ordistance>=line.length:return[LineString(line)]coords=list(line.coords)fori,pinenumerate(coords):pd=line.project(Point(p))ifpd==distance:return[LineString(coords[:i+1]),LineString(coords[i:])]ifpd>distance:cp=line.interpolate(distance)return[LineString(coords[:i]+[(cp.x,cp.y)]),LineString([(cp.x,cp.y)]+coords[i:])]

A valid LinearRing may not cross itself or touch itself at a single point. A
valid Polygon may not possess any overlapping exterior or interior rings. A
valid MultiPolygon may not collect any overlapping polygons. Operations on
invalid features may fail.

The two points above are close enough that the polygons resulting from the
buffer operations (explained in a following section) overlap.

Note

The is_valid predicate can be used to write a validating decorator that
could ensure that only valid objects are returned from a constructor
function.

fromfunctoolsimportwrapsdefvalidate(func):@wraps(func)defwrapper(*args,**kwargs):ob=func(*args,**kwargs)ifnotob.is_valid:raiseTopologicalError("Given arguments do not determine a valid geometric object")returnobreturnwrapper

Standard binary predicates are implemented as methods. These predicates
evaluate topological, set-theoretic relationships. In a few cases the results
may not be what one might expect starting from different assumptions. All take
another geometric object as argument and return True or False.

Returns True if the set-theoretic boundary, interior, and exterior
of the object coincide with those of the other.

The coordinates passed to the object constructors are of these sets, and
determine them, but are not the entirety of the sets. This is a potential
“gotcha” for new users. Equivalent lines, for example, can be constructed
differently.

Used in a sorted()key, within() makes it easy to spatially sort
objects. Let’s say we have 4 stereotypic features: a point that is contained by
a polygon which is itself contained by another polygon, and a free spirited
point contained by none

that we’d prefer to have ordered as [d,c,c,b,a] in reverse containment
order. As explained in the Python Sorting HowTo, we can define a key
function that operates on each list element and returns a value for comparison.
Our key function will be a wrapper class that implements __lt__() using
Shapely’s binary within() predicate.

As the howto says, the less than comparison is guaranteed to be used in
sorting. That’s what we’ll rely on to spatially sort, and the reason why we use
within() in reverse instead of contains(). Trying it out on features
d and c, we see that it works.

Returns a string representation of the DE-9IM matrix of relationships
between an object’s interior, boundary, exterior and those of another
geometric object.

The named relationship predicates (contains(), etc.) are typically
implemented as wrappers around relate().

Two different points have mainly F (false) values in their matrix; the
intersection of their external sets (the 9th element) is a 2 dimensional
object (the rest of the plane). The intersection of the interior of one with
the exterior of the other is a 0 dimensional object (3rd and 7th elements
of the matrix).

>>> Point(0,0).relate(Point(1,1))'FF0FFF0F2'

The matrix for a line and a point on the line has more “true” (not F)
elements.

Shapely can not represent the difference between an object and a lower
dimensional object (such as the difference between a polygon and a line or
point) as a single object, and in these cases the difference method returns a
copy of the object named self.

Returns an approximate representation of all points within a given distance
of the this geometric object.

The styles of caps are specified by integer values: 1 (round), 2 (flat),
3 (square). These values are also enumerated by the object
shapely.geometry.CAP_STYLE (see below).

The styles of joins between offset segments are specified by integer values:
1 (round), 2 (mitre), and 3 (bevel). These values are also enumerated by the
object shapely.geometry.JOIN_STYLE (see below).

Returns a representation of the smallest convex Polygon containing all the
points in the object unless the number of points in the object is less than
three. For two points, the convex hull collapses to a LineString; for 1, a
Point.

Returns a LineString or MultiLineString geometry at a distance from the
object on its right or its left side.

Distance must be a positive float value. The side parameter may be ‘left’ or
‘right’. The resolution of the offset around each vertex of the object is
parameterized as in the buffer method.

The join style is for outside corners between line segments. Accepted integer
values are 1 (round), 2 (mitre), and 3 (bevel). See also
shapely.geometry.JOIN_STYLE.

Severely mitered corners can be controlled by the mitre_limit parameter
(spelled in British English, en-gb). The ratio of the distance from the
corner to the end of the mitred offset corner is the miter ratio. Corners
with a ratio which exceed the limit will be beveled.

All points in the simplified object will be within the tolerance distance of
the original geometry. By default a slower algorithm is used that preserves
topology. If preserve topology is set to False the much quicker
Douglas-Peucker algorithm [6] is used.

A collection of affine transform functions are in the shapely.affinity
module, which return transformed geometries by either directly supplying
coefficients to an affine transformation matrix, or by using a specific, named
transform (rotate, scale, etc.). The functions can be used with all
geometry types (except GeometryCollection), and 3D types are either
preserved or supported by 3D affine transformations.

Dangles are edges which have one or both ends which are not incident on
another edge endpoint. Cut edges are connected at both ends but do not
form part of polygon. Invalid ring lines form rings which are invalid
(bowties, etc).

Prepared geometries instances have the following methods: contains,
contains_properly, covers, and intersects. All have exactly the
same arguments and usage as their counterparts in non-prepared geometric
objects.

A Well Known Text (WKT) or Well Known Binary (WKB) representation [1] of
any geometric object can be had via its wkt or wkb attribute.
These representations allow interchange with many GIS programs. PostGIS, for
example, trades in hex-encoded WKB.

The shapely.wkt and shapely.wkb modules provide dumps() and loads()
functions that work almost exactly as their pickle and simplejson module
counterparts. To serialize a geometric object to a binary or text string, use
dumps(). To deserialize a string and get a new geometric object of the
appropriate type, use loads().

The shapely.geometry.asShape() family of functions can be used to wrap
Numpy coordinate arrays so that they can then be analyzed using Shapely while
maintaining their original storage. A 1 x 2 array can be adapted to a point

Shapely uses the GEOS library for all operations. GEOS is written in C++ and
used in many applications and you can expect that all operations are highly
optimized. The creation of new geometries with many coordinates, however,
involves some overhead that might slow down your code.

New in version 1.2.10.

The shapely.speedups module contains performance enhancements written in
C. They are automaticaly installed when Python has access to a compiler and
GEOS development headers during installation.

You can check if the speedups are installed with the available
attribute. The constructor speedups are disabled by default. To enable the
speedups call enable(). You can revert to the default implementation with
disable().

David H. Douglas and Thomas K. Peucker,
“Algorithms for the Reduction of the Number of Points Required to Represent
a Digitized Line or its Caricature,” Cartographica: The International
Journal for Geographic Information and Geovisualization, vol. 10, Dec.
1973, pp. 112-122.