root/Shapely/trunk/manual/manual.txt

==================
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`__.

.. __: http://creativecommons.org/licenses/by/3.0/us/

:abstract: This document describes the Shapely Python package for programming
  with geospatial geometries.

.. sectnum::

.. contents::


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_.

.. _GEOS: http://geos.refractions.net
.. _Java Topology Suite: http://www.jump-project.org/project.php?PID=JTS&SID=OVER
.. _record: http://pypi.python.org/pypi/Shapely
.. _wiki: http://trac.gispython.org/projects/PCL/wiki/Shapely


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_.

.. _illustration: http://www.jump-project.org/project.php?PID=JTS&SID=OVER#spatialdatatypes

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


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 (**a**\ 1, ..., **a**\ M) and (**b**\ 1, ...,
**b**\ N) be ordered sequences of *M* and *N* such coordinate sequences,
defining lines or rings.

Points
++++++

The point factory *Point* takes a coordinate sequence parameter

.. code-block:: python

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

The alternate form is to pass individual coordinate parameters

.. code-block:: python

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

LineStrings
+++++++++++

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

.. code-block:: python

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

Polygons
++++++++

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

.. code-block:: python

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

If **a**\ 1 is not exactly equal to **a**\ M, the factory will close the ring.
The following (unit square) polygons are therefore topologically equal

.. code-block:: python

  >>> 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*)

.. code-block:: python

  >>> 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.


Multipart Geometry Factories
----------------------------

MultiPoints
+++++++++++

An *N*\ -point geometry is created by passing an unordered sequence of
coordinate sequences [**c**\ 1, ..., **c**\ N]

.. code-block:: python

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

MultiLineStrings
++++++++++++++++

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

.. code-block:: python

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

MultiPolygons
+++++++++++++

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

.. code-block:: python

  >>> 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

.. code-block:: python

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


Null Geometries
---------------

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

.. code-block:: python

  >>> 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.

.. code-block:: python

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


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 |illustration2|_.

.. |illustration2| replace:: illustration
.. _illustration2: http://www.jump-project.org/project.php?PID=JTS&SID=OVER#spatialanalysismethods

Example Geometries
++++++++++++++++++

.. code-block:: python

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

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:

.. code-block:: python

  >>> 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

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.

.. code-block:: python

  >>> polygon.boundary
  <shapely.geometry.linestring.LineString object at ...>
  >>> line_b.boundary
  <shapely.geometry.multipoint.MultiPoint object at ...>
  >>> point_r.boundary.is_empty
  True

Centroid
++++++++

.centroid : geometry
  Returns the centroid, or geometric center of the polygon.

.. code-block:: python

  >>> centroid_point = polygon.centroid
  >>> centroid_point.wkt
  'POINT (-0.0000000000000000 -0.0000000000000000)'

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:

.. code-block:: python

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

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).

.. code-block:: python

  >>> hull = multi_point.convex_hull
  >>> polygon.difference(hull)
  <shapely.geometry.polygon.Polygon object at ...>

Envelope
++++++++

.envelope : geometry
  Returns the geometry's rectangular polygon envelope.

.. code-block:: python

  >>> polygon.envelope
  <shapely.geometry.polygon.Polygon object at ...>

Intersection
++++++++++++

.intersection(other) : geometry
  Returns the intersection of one geometry and the *other* geometry.

.. code-block:: python

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

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.

.. code-block:: python

  >>> polygon.symmetric_difference(hull)
  <shapely.geometry.multipolygon.MultiPolygon object at ...>

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:

.. code-block:: python

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


Other Operations
----------------

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.

.. code-block:: python

  >>> 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 ...>]


Unary Spatial Predicates
------------------------

These are implemented as Python attributes.

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.  

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.

Is Ring
+++++++

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

Has Z
+++++

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

Examples
++++++++

.. code-block:: python

  >>> 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.)


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.

Contains
++++++++

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

.. code-block:: python

  >>> polygon.contains(point_b)
  True

Crosses
+++++++

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

.. code-block:: python

  >>> line_b.crosses(polygon)
  True

Disjoint
++++++++

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

.. code-block:: python

  >>> polygon.disjoint(point_r)
  True

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__*.

Intersects
++++++++++

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

.. code-block:: python

  >>> polygon.intersects(point_b)
  True

Touches
+++++++

.touches(other) : bool
  True if geometries *only* touch. The least inclusive predicate.

.. code-block:: python

  >>> polygon.touches(line_g)
  True
  >>> polygon.touches(line_b)
  False

Within
++++++

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


General Methods
---------------

Distance
++++++++

.distance(other) : geometry
  The minimum distance from one geometry to the other.

.. code-block:: python

  >>> Point(0,0).distance(Point(1,1))
  1.4142135623730951


Scalar Properties
-----------------

Area
++++

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

Bounds
++++++

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

Length
++++++

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

Examples
++++++++

.. code-block:: python

  >>> 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


Geometry Parts and Coordinates
==============================

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.

.. code-block:: python

  >>> 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 (**a**\ 1, ...,
**a**\ M) to the coords attribute.

.. code-block:: python

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

For line strings:

.. code-block:: python

  >>> 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)]


Polygon Rings
-------------

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

.. code-block:: python

  >>> 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.


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:

.. code-block:: python

  >>> 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.


Point Coordinates
-----------------

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

.. code-block:: python

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


Interoperation
==============

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

Well-known Formats
------------------

Well-known Text (WKT)
+++++++++++++++++++++

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

.. code-block:: python

  >>> 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

.. code-block:: python

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

.. code-block:: python

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

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:

.. code-block:: python

  >>> 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 ...>


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**:

.. code-block:: python

  >>> 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

.. code-block:: python

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

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

.. code-block:: python

  >>> 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


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:

.. code-block:: python

  >>> 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

.. code-block:: python

  >>> 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

.. code-block:: python

  >>> 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.


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:

.. code-block:: python

  >>> 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:

.. code-block:: python

  >>> 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.

.. _Python geo interface: http://trac.gispython.org/projects/PCL/wiki/PythonGeoInterface


Advanced Features
=================

Iterative Operations
--------------------

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

Contains
++++++++

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

.. code-block:: python

  >>> 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*.

.. code-block:: python

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


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).