The Shapely Manual

Author: Sean Gillies
Address: 
sgillies@frii.com
Revision: 1.0.5
Date: 20 May 2008
Copyright: This work is licensed under a Creative Commons Attribution 3.0 United States License.
abstract:This document describes the Shapely Python package for programming with geospatial geometries.

Contents

1   Background

Shapely is a Python package for programming with geospatial geometries. It is based on the GEOS library, a port of the Java Topology Suite. Please refer to the JTS overview for definitions and illustrations of the various geometry types, operations, and predicates. See also the Shapely wiki and Python Package Index record.

2   Geometries

The basic, standard GIS geometry model consists of single points, line strings, and polygons, homogeneous multi-point, multi-line string, and multi-polygon collections, and heterogeneous geometry collections. See the JTS illustration.

Shapely 1.0 does not provide circles, arcs, splines or the like.

2.1   Factories

Geometries can be created in the typical Python fashion, using the geometry classes themselves as factories.

Pseudo-code blocks in this section will use the following notation. Let a be a Cartesian x, y, and optional z coordinate sequence. The coordinates values must be numeric types. Let (a1, ..., aM) and (b1, ..., bN) be ordered sequences of M and N such coordinate sequences, defining lines or rings.

2.1.1   Points

The point factory Point takes a coordinate sequence parameter

>>> from shapely.geometry import Point
>>> point = Point(a)

The alternate form is to pass individual coordinate parameters

>>> point = Point(x0, y0 [, z0])

2.1.2   LineStrings

To create a line string, pass in an ordered sequence of coordinate sequences:

>>> from shapely.geometry import LineString
>>> line = LineString((a1, ..., aM))

2.1.3   Polygons

A polygon with only an exterior boundary and no holes is created by passing the sequence representation of a closed ring

>>> from shapely.geometry import Polygon
>>> polygon = Polygon((a1, ..., aM))

If a1 is not exactly equal to aM, the factory will close the ring. The following (unit square) polygons are therefore topologically equal

>>> polygon1 = Polygon(((0, 0), (0, 1), (1, 1), (1, 0), (0, 0)))
>>> polygon2 = Polygon(((0, 0), (0, 1), (1, 1), (1, 0)))

To create a polygon with interior boundaries pass a sequence of rings to the second parameter (holes)

>>> polygon = Polygon((a1, ..., aM), [(b1, ..., bN), ...])

Rings must be non-crossing, ordered coordinate sequences. The order may be clockwise or counter-clockwise. The resulting topology is independent of the order.

2.2   Multipart Geometry Factories

2.2.1   MultiPoints

An N-point geometry is created by passing an unordered sequence of coordinate sequences [c1, ..., cN]

>>> from shapely.geometry import MultiPoint
>>> points = MultiPoint([c1, ..., cN])

2.2.2   MultiLineStrings

A multi-line geometry is created by passing a sequence of representations of lines

>>> from shapely.geometry import MultiLineString
>>> lines = MultiLineString([(a1, ..., aM), (b1, ..., bN), ...])

2.2.3   MultiPolygons

A multi-polygon geometry is created by passing a sequence of exterior ring and hole list tuples

>>> from shapely.geometry import MultiPolygon
>>> lines = MultiPolygon([((a1, ..., aM), [(b1, ..., bN), ...]), ...])

More explicit notation for the exterior and interior boundaries (or shells and holes) makes usage more clear

>>> shell = (a1, ..., aM)
>>> holes = [(b1, ..., bN), ...]
>>> lines = MultiPolygon([(shell, holes), ...])

2.3   Null Geometries

Null geometries can be created by calling the factories with no arguments, but almost nothing can be done with a null geometry.

>>> line_null = LineString()
>>> line_null.length
Traceback (most recent call last):
...
ValueError: Null geometry supports no operations

The coordinates of a null geometry can be set (see Section 3), after which the geometry is no longer null.

>>> l_null.coords = [(0, 0), (1, 1)]
>>> print l_null.length
1.414...

2.4   Constructive Spatial Analysis Methods

There are methods of geometry classes that also serve as factories for new geometries. It is important to note that these are topological and not point-wise operations, and therefore may produce results that are not what one might expect from operations on Python sets.

See also the JTS illustration.

2.4.1   Example Geometries

>>> polygon = Polygon(((-1.0, -1.0), (-1.0, 1.0), (1.0, 1.0), (1.0, -1.0)))
>>> point_r = Point(-1.5, 1.2)
>>> point_g = Point(-1.0, 1.0)
>>> point_b = Point(-0.5, 0.5)
>>> line_r = LineString(((-0.5, 0.5), (0.5, 0.5)))
>>> line_g = LineString(((1.0, -1.0), (1.8, 0.5)))
>>> line_b = LineString(((-1.8, -1.2), (1.8, 0.5)))

2.4.2   Buffer

.buffer(width, quadsegs=16) : geometry
Returns a buffer region having the given width and with a specified number of segments used to approximate curves.

The default result of buffering a point is an N-gon approximation of a circle:

>>> buffered = point_r.buffer(1.0)
>>> buffered
<shapely.geometry.polygon.Polygon object at ...>
>>> buffered.length
6.2806623139097271
>>> buffered.area
3.1365484905463727
>>> len(buffered.exterior.coords)
66

2.4.3   Boundary

.boundary : geometry
Returns a lower dimension geometry. 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.
>>> polygon.boundary
<shapely.geometry.linestring.LineString object at ...>
>>> line_b.boundary
<shapely.geometry.multipoint.MultiPoint object at ...>
>>> point_r.boundary.is_empty
True

2.4.4   Centroid

.centroid : geometry
Returns the centroid, or geometric center of the polygon.
>>> centroid_point = polygon.centroid
>>> centroid_point.wkt
'POINT (-0.0000000000000000 -0.0000000000000000)'

2.4.5   Convex Hull

.convex_hull : geometry
Imagine an elastic band stretched around the geometry: that's a convex hull, more or less.

For example, collect the three points into a multi-point geometry, and get the triangular polygon that is their convex hull:

>>> multi_point = point_r.union(point_g)
>>> multi_point = multi_point.union(point_b)
>>> multi_point.convex_hull
<shapely.geometry.polygon.Polygon object at ...>

2.4.6   Difference

.difference(other) : geometry
Returns a geometry representing the points making up this geometry that do not make up other. Note that A.difference(B) is not necessarily equal to B.difference(A).
>>> hull = multi_point.convex_hull
>>> polygon.difference(hull)
<shapely.geometry.polygon.Polygon object at ...>

2.4.7   Envelope

.envelope : geometry
Returns the geometry's rectangular polygon envelope.
>>> polygon.envelope
<shapely.geometry.polygon.Polygon object at ...>

2.4.8   Intersection

.intersection(other) : geometry
Returns the intersection of one geometry and the other geometry.
>>> polygon.intersection(hull)
<shapely.geometry.polygon.Polygon object at ...>

2.4.9   Symmetric Difference

.symmetric_difference(other) : geometry
Returns a geometry combining the points in this geometry not in other, and the points in other not in this geometry.
>>> polygon.symmetric_difference(hull)
<shapely.geometry.multipolygon.MultiPolygon object at ...>

2.4.10   Union

.union(other) : geometry
Returns the union of one geometry and the other geometry.

Point unions were demonstrated above under convex hull. The union of polygons will be a polygon or a multi-polygon depending on whether they intersect or not:

>>> hull.union(polygon)
<shapely.geometry.polygon.Polygon object at ...>

2.5   Other Operations

2.5.1   Polygonization

shapely.ops.polygonize(lines) : iterator
Returns an iterator over polygons constructed from the lines iterator. The elements of lines may be Shapely geometries, objects that provide the geo interface, or Numpy arrays or Python sequences shaped like LineStrings.
>>> from shapely.ops import polygonize
>>> lines = [
...     ((0, 0), (1, 1)),
...     ((0, 0), (0, 1)),
...     ((0, 1), (1, 1)),
...     ((1, 1), (1, 0)),
...     ((1, 0), (0, 0))
...     ]
>>> result = polygonize(lines)
>>> list(result.geoms)
[<shapely.geometry.polygon.Polygon object at ...>, <shapely.geometry.polygon.Polygon object at ...>]

2.6   Unary Spatial Predicates

These are implemented as Python attributes.

2.6.1   Is Empty

.is_empty : bool
True if the set of points in this geometry is empty, else False. For more details, see http://geos.refractions.net/ro/doxygen_docs/html/classgeos_1_1geom_1_1Geometry.html#a17.

2.6.2   Is Valid

.is_valid : bool
True if the geometry is valid (definition depends on sub-class), else False. For more details, see http://geos.refractions.net/ro/doxygen_docs/html/classgeos_1_1geom_1_1Geometry.html#a16.

2.6.3   Is Ring

.is_ring : bool
True if the geometry is a closed ring, else False.

2.6.4   Has Z

.has_z : bool
True if the geometry's coordinate sequence(s) have z values (are 3-dimensional)

2.6.5   Examples

>>> polygon.is_empty
False
>>> polygon.is_valid
True
>>> polygon.is_ring
False
>>> polygon.boundary.is_ring
True
>>> polygon.has_z
False

(Note: that last return value exposes a bug in GEOS 2.2.3.)

2.7   Binary Spatial Predicates

All of these methods take a single positional argument, an other geometry. It is important to note that these are topological and not point-wise operations, and therefore may produce results that are not what one might expect from operations on Python.

2.7.1   Contains

.contains(other) : bool
True if the geometry is spatially within, without touching. Applies to all types of geometries.
>>> polygon.contains(point_b)
True

2.7.2   Crosses

.crosses(other) : bool
Only linear geometries (lines, rings, polygon boundaries) may ever cross. No geometry may ever cross a point.
>>> line_b.crosses(polygon)
True

2.7.3   Disjoint

.disjoint(other) : bool
True if geometries do not spatially relate in any way, else False. See the complementary intersects.
>>> polygon.disjoint(point_r)
True

2.7.4   Equals

.equals(other) : bool
Two geometries are topologically equal if their interiors intersect and no part of the interior or boundary of one geometry intersects the exterior of the other. Not to be confused with Python's __equals__.

2.7.5   Intersects

.intersects(other) : bool
This predicate is the complement of disjoint: geometries that do not intersect are disjoint. Intersects is the most inclusive predicate.
>>> polygon.intersects(point_b)
True

2.7.6   Touches

.touches(other) : bool
True if geometries only touch. The least inclusive predicate.
>>> polygon.touches(line_g)
True
>>> polygon.touches(line_b)
False

2.7.7   Within

.within(other): bool
The inverse of contains.

2.8   General Methods

2.8.1   Distance

.distance(other) : geometry
The minimum distance from one geometry to the other.
>>> Point(0,0).distance(Point(1,1))
1.4142135623730951

2.9   Scalar Properties

2.9.1   Area

.area : float
Area of the geometry, unitless. Non-zero only for surfaces (polygons, multi-polygons).

2.9.2   Bounds

.bounds : tuple
The geometry's (minx, miny, maxx, maxy) bounding box.

2.9.3   Length

.length : float
Length of the geometry, unitless. Non-zero only for linear geometries (line strings, rings, polygon boundaries)

2.9.4   Examples

>>> polygon.area
4.0
>>> polygon.bounds
(-1.0, -1.0, 1.0, 1.0)
>>> polygon.length
8.0
>>> line_r.length
1.0
>>> line_b.length
3.9812058474788765

3   Geometry Parts and Coordinates

3.1   Coordinate Sequences

The coordinates of points, line strings, and polygon rings can be accessed through the coords attribute of a geometry. Coords is an iterator over coordinate tuples.

>>> point_r.coords
<shapely.geometry.base.CoordinateSequence object at ...>
>>> len(point_r.coords)
1
>>> point_r.coords[0]
(-1.5, 1.2)
>>> list(point_r.coords)
[(-1.5, 1.2)]

The coordinate sequence can be modifed by assigning a sequence (a1, ..., aM) to the coords attribute.

>>> point_new = Point(0, 0)
>>> point_new.coords = (1, 1)
>>> list(point_new.coords)
[(1.0, 1.0)]

For line strings:

>>> line_new = LineString([(0,0), (1,1)])
>>> line_new.coords = [(1,1), (2,2)]
>>> list(line_new.coords)
[(1.0, 1.0), (2.0, 2.0)]

3.2   Polygon Rings

The exterior boundary of a polygon can be accessed through the exterior attribute of the geometry object.

>>> polygon.exterior
<shapely.geometry.polygon.LinearRing object at ...>
>>> list(polygon.exterior.coords)
[(-1.0, -1.0), (-1.0, 1.0), (1.0, 1.0), (1.0, -1.0), (-1.0, -1.0)]

The interior boundaries (or holes) of a polygon can be accessed through the interiors attribute, which is a list of rings.

3.3   Sub-geometries

The parts of a multi-part geometry can be accessed through the geoms attribute of the geometry object, which is an iterator over the sub-geometries:

>>> multi_point.geoms
<shapely.geometry.base.GeometrySequence object at ...>
>>> len(multi_point.geoms)
3
>>> from pprint import pprint
>>> pprint(list(multi_point.geoms))
[<shapely.geometry.point.Point object at ...>,
 <shapely.geometry.point.Point object at ...>,
 <shapely.geometry.point.Point object at ...>]

The coordinate sequences of these sub-geometries can then be accessed as described above.

3.4   Point Coordinates

For the sake of convenience the coordinate values of points can be accessed read-only via x, y, and z attributes:

>>> point = Point(1.0, 1.0)
>>> point.x
1.0
>>> point.y
1.0

4   Interoperation

Shapely provides 4 avenues for interoperation with other Python and GIS software.

4.1   Well-known Formats

4.1.1   Well-known Text (WKT)

The WKT representation of any geometry object can be had via the wkt attribute:

>>> point_r.wkt
'POINT (-1.5000000000000000 1.2000000000000000)'

Hex-encode that string and you have a value that can be conveniently inserted directly into PostGIS

>>> point_r.wkt.encode('hex')
'504f494e5420282d312e3530303030303030303030303030303020312e3230303030303030303030303030303029'

New geometries can be created from WKT representations using the shapely.wkt.loads factory (inspired by the pickle module)

>>> from shapely.wkt import loads
>>> loads('POINT (0 0)')
<shapely.geometry.point.Point object at ...>

4.1.2   Well-known Binary (WKB)

The WKB representation of any geometry object can be had via the wkb attribute. New geometries can be created from WKB data using the shapely.wkb.loads factory. Use this format to interoperate with ogr.py:

>>> import ogr
>>> from shapely.wkb import loads
>>> source = ogr.Open("/tmp/world_borders.shp")
>>> borders = source.GetLayerByName("world_borders")
>>> feature = borders.GetNextFeature()
>>> loads(feature.GetGeometryRef().ExportToWkb())
<shapely.geometry.polygon.Polygon object at ...>

4.2   Python Sequences

Python sequence data can be analyzed as Shapely geometries using the shapely.geometry.as* adapters while leaving the data in its original storage. A pair of float can be treated as a point with asPoint:

>>> from shapely.geometry import asPoint
>>> coords = [3.0, 4.0]
>>> pa = asPoint(coords)
>>> pa.wkt
'POINT (3.0000000000000000 4.0000000000000000)'

Move the coordinates and watch the geometry adapter change

>>> coords[0] = 1.0
>>> pa.wkt
'POINT (1.0000000000000000 4.0000000000000000)'

The asLineString adapter works much the same. The asPolygon adapter is used like

>>> from shapely.geometry import asPolygon
>>> coords = [[0.0, 0.0], [0.0, 1.0], [1.0, 1.0], [1.0, 0.0]]
>>> hole_coords = [((0.1,0.1), (0.1,0.2), (0.2,0.2), (0.2,0.1))]
>>> pa = asPolygon(coords, hole_coords)
>>> len(pa.exterior.coords)
5
>>> len(pa.interiors)
1
>>> len(pa.interiors[0].coords)
5

4.3   Numpy Array Interface

Shapely geometries provide the Numpy array interface which means that points, line strings, and polygon rings can be used as Numpy arrays:

>>> from numpy import array
>>> a = array(polygon.exterior)
>>> a
array([[-1., -1.],
       [-1.,  1.],
       [ 1.,  1.],
       [ 1., -1.],
       [-1., -1.]])

The numpy.asarray function does not copy coordinate values at the price of slower numpy access to coordinates.

The shapely.geometry.as* functions can also be used to wrap numpy arrays, which can then be analyzed using Shapely while maintaining their original storage. A 1 x 2 array can be adapted to a point

>>> a = array([1.0, 2.0])
>>> pa = asPoint(a)
>>> pa.wkt
'POINT (1.0000000000000000 2.0000000000000000)'

and a N x 2 array can be adapted to a line string

>>> from shapely.geometry import asLineString
>>> a = array([[1.0, 2.0], [3.0, 4.0]])
>>> la = asLineString(a)
>>> la.wkt
'LINESTRING (1.0000000000000000 2.0000000000000000, 3.0000000000000000 4.0000000000000000)'

There is no Numpy array representation of a polygon.

4.4   Python Geo Interface

Any object that provides the GeoJSON-like Python geo interface can be adapted and used as a Shapely geometry using the shapely.geometry.asShape function. For example, a dictionary:

>>> from shapely.geometry import asShape
>>> d = {"type": "Point", "coordinates": (0.0, 0.0)}
>>> shape = asShape(d)
>>> shape.geom_type
'Point'
>>> list(shape.coords)
[(0.0, 0.0)]

Or a simple placemark-type object:

>>> class GeoThing(object):
...     def __init__(self, d):
...         self.__geo_interface__ = d
>>> thing = GeoThing({"type": "Point", "coordinates": (0.0, 0.0)})
>>> shape = asShape(thing)
>>> shape.geom_type
'Point'
>>> list(shape.coords)
[(0.0, 0.0)]

If you want to copy coordinate data to a new geometry, use the shapely.geometry.shape function instead.

5   Advanced Features

5.1   Iterative Operations

Shapely provides functions for efficient operations on large sets of geometries.

5.1.1   Contains

To find the subset of points that are contained within a polygon, use shapely.iterops.contains:

>>> from shapely.geometry import Polygon
>>> from shapely.geometry import Point
>>> coords = ((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0))
>>> polygon = Polygon(coords)
>>> points = [Point(0.5, 0.5), Point(2.0, 2.0)]
>>> from shapely import iterops
>>> list(iterops.contains(polygon, points, True))
[<shapely.geometry.point.Point object at ...>]

The second parameter to iterops.contains can be any kind of iterator, even a generator of objects. If it yields tuples, then the second element of the tuple will be ultimately yielded from iterops.contains.

>>> list(iterops.contains(polygon, iter((p, p.wkt) for p in points)))
['POINT (0.5000000000000000 0.5000000000000000)']

6   Credits

Shapely is written by Sean Gillies with contributions from Howard Butler, Kai Lautaportti (Hexagon IT), Frédéric Junod (Camptocamp SA), Eric Lemoine (Camptocamp SA) and ctypes tips from Justin Bronn (GeoDjango).