Menu
▾
▴
SF.net SVN: matplotlib: [4226] trunk/toolkits/basemap-testing/lib/shapely
|
From: <js...@us...> - 2007-11-12 16:55:01
|
Revision: 4226
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4226&view=rev
Author: jswhit
Date: 2007-11-12 08:54:58 -0800 (Mon, 12 Nov 2007)
Log Message:
-----------
Add Shapely
Added Paths:
-----------
trunk/toolkits/basemap-testing/lib/shapely/geometry/
trunk/toolkits/basemap-testing/lib/shapely/geometry/__init__.py
trunk/toolkits/basemap-testing/lib/shapely/geometry/base.py
trunk/toolkits/basemap-testing/lib/shapely/geometry/collection.py
trunk/toolkits/basemap-testing/lib/shapely/geometry/geo.py
trunk/toolkits/basemap-testing/lib/shapely/geometry/linestring.py
trunk/toolkits/basemap-testing/lib/shapely/geometry/multilinestring.py
trunk/toolkits/basemap-testing/lib/shapely/geometry/multipoint.py
trunk/toolkits/basemap-testing/lib/shapely/geometry/multipolygon.py
trunk/toolkits/basemap-testing/lib/shapely/geometry/point.py
trunk/toolkits/basemap-testing/lib/shapely/geometry/polygon.py
Added: trunk/toolkits/basemap-testing/lib/shapely/geometry/__init__.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/shapely/geometry/__init__.py (rev 0)
+++ trunk/toolkits/basemap-testing/lib/shapely/geometry/__init__.py 2007-11-12 16:54:58 UTC (rev 4226)
@@ -0,0 +1,9 @@
+from geo import asShape
+from point import Point, asPoint
+from linestring import LineString, asLineString
+from polygon import Polygon, asPolygon
+from multipoint import MultiPoint, asMultiPoint
+from multilinestring import MultiLineString, asMultiLineString
+from multipolygon import MultiPolygon, asMultiPolygon
+from collection import GeometryCollection
+
Added: trunk/toolkits/basemap-testing/lib/shapely/geometry/base.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/shapely/geometry/base.py (rev 0)
+++ trunk/toolkits/basemap-testing/lib/shapely/geometry/base.py 2007-11-12 16:54:58 UTC (rev 4226)
@@ -0,0 +1,358 @@
+"""
+Base geometry class and utilities.
+"""
+
+from ctypes import string_at, byref, c_int, c_size_t, c_char_p, c_double
+import sys
+
+from shapely.geos import lgeos
+from shapely.predicates import BinaryPredicate, UnaryPredicate
+from shapely.topology import BinaryTopologicalOp, UnaryTopologicalOp
+
+
+GEOMETRY_TYPES = [
+ 'Point',
+ 'LineString',
+ 'LinearRing',
+ 'Polygon',
+ 'MultiPoint',
+ 'MultiLineString',
+ 'MultiPolygon',
+ 'GeometryCollection'
+ ]
+
+def geometry_type_name(g):
+ return GEOMETRY_TYPES[lgeos.GEOSGeomTypeId(g)]
+
+# Abstract geometry factory for use with topological methods below
+
+def geom_factory(g):
+ if not g:
+ raise ValueError, "No Shapely geometry can be created from this null value"
+ ob = BaseGeometry()
+ geom_type = geometry_type_name(g)
+ # TODO: check cost of dynamic import by profiling
+ mod = __import__(
+ 'shapely.geometry',
+ globals(),
+ locals(),
+ [geom_type],
+ )
+ ob.__class__ = getattr(mod, geom_type)
+ ob._geom = g
+ ob._ndim = 2 # callers should be all from 2D worlds
+ return ob
+
+
+class CoordinateSequence(object):
+
+ _cseq = None
+ _ndim = None
+ _length = 0
+ index = 0
+
+ def __init__(self, geom):
+ self._cseq = lgeos.GEOSGeom_getCoordSeq(geom._geom)
+ self._ndim = geom._ndim
+
+ def __iter__(self):
+ self.index = 0
+ self._length = self.__len__()
+ return self
+
+ def next(self):
+ dx = c_double()
+ dy = c_double()
+ dz = c_double()
+ i = self.index
+ if i < self._length:
+ lgeos.GEOSCoordSeq_getX(self._cseq, i, byref(dx))
+ lgeos.GEOSCoordSeq_getY(self._cseq, i, byref(dy))
+ if self._ndim == 3: # TODO: use hasz
+ lgeos.GEOSCoordSeq_getZ(self._cseq, i, byref(dz))
+ self.index += 1
+ return (dx.value, dy.value, dz.value)
+ else:
+ self.index += 1
+ return (dx.value, dy.value)
+ else:
+ raise StopIteration
+
+ def __len__(self):
+ cs_len = c_int(0)
+ lgeos.GEOSCoordSeq_getSize(self._cseq, byref(cs_len))
+ return cs_len.value
+
+ def __getitem__(self, i):
+ M = self.__len__()
+ if i + M < 0 or i >= M:
+ raise IndexError, "index out of range"
+ if i < 0:
+ ii = M + i
+ else:
+ ii = i
+ dx = c_double()
+ dy = c_double()
+ dz = c_double()
+ lgeos.GEOSCoordSeq_getX(self._cseq, ii, byref(dx))
+ lgeos.GEOSCoordSeq_getY(self._cseq, ii, byref(dy))
+ if self._ndim == 3: # TODO: use hasz
+ lgeos.GEOSCoordSeq_getZ(self._cseq, ii, byref(dz))
+ return (dx.value, dy.value, dz.value)
+ else:
+ return (dx.value, dy.value)
+
+
+class GeometrySequence(object):
+
+ _factory = None
+ _geom = None
+ _ndim = None
+ _index = 0
+ _length = 0
+
+ def __init__(self, geom, type):
+ self._factory = type
+ self._geom = geom._geom
+ self._ndim = geom._ndim
+
+ def __iter__(self):
+ self._index = 0
+ self._length = self.__len__()
+ return self
+
+ def next(self):
+ if self._index < self.__len__():
+ g = self._factory()
+ g._owned = True
+ g._geom = lgeos.GEOSGetGeometryN(self._geom, self._index)
+ self._index += 1
+ return g
+ else:
+ raise StopIteration
+
+ def __len__(self):
+ return lgeos.GEOSGetNumGeometries(self._geom)
+
+ def __getitem__(self, i):
+ M = self.__len__()
+ if i + M < 0 or i >= M:
+ raise IndexError, "index out of range"
+ if i < 0:
+ ii = M + i
+ else:
+ ii = i
+ g = self._factory()
+ g._owned = True
+ g._geom = lgeos.GEOSGetGeometryN(self._geom, ii)
+ return g
+
+ @property
+ def _longest(self):
+ max = 0
+ for g in iter(self):
+ l = len(g.coords)
+ if l > max:
+ max = l
+
+
+class BaseGeometry(object):
+
+ """Provides GEOS spatial predicates and topological operations.
+ """
+
+ _geom = None
+ _ctypes_data = None
+ _ndim = None
+ _crs = None
+ _owned = False
+
+ def __init__(self):
+ self._geom = lgeos.GEOSGeomFromWKT(c_char_p('GEOMETRYCOLLECTION EMPTY'))
+
+ def __del__(self):
+ if self._geom is not None and not self._owned:
+ lgeos.GEOSGeom_destroy(self._geom)
+
+ def __str__(self):
+ return self.to_wkt()
+
+ def __eq__(self, other):
+ return self.equals(other)
+
+ def __ne__(self, other):
+ return not self.equals(other)
+
+ # To support pickling
+
+ def __reduce__(self):
+ return (self.__class__, (), self.to_wkb())
+
+ def __setstate__(self, state):
+ self._geom = lgeos.GEOSGeomFromWKB_buf(
+ c_char_p(state),
+ c_size_t(len(state))
+ )
+
+ # Array and ctypes interfaces
+
+ @property
+ def ctypes(self):
+ """Return a ctypes representation.
+
+ To be overridden by extension classes."""
+ raise NotImplementedError
+
+ @property
+ def array_interface_base(self):
+ if sys.byteorder == 'little':
+ typestr = '<f8'
+ elif sys.byteorder == 'big':
+ typestr = '>f8'
+ else:
+ raise ValueError, "Unsupported byteorder: neither little nor big-endian"
+ return {
+ 'version': 3,
+ 'typestr': typestr,
+ 'data': self.ctypes,
+ }
+
+ @property
+ def __array_interface__(self):
+ """Provide the Numpy array protocol."""
+ raise NotImplementedError
+
+ def get_coords(self):
+ return CoordinateSequence(self)
+
+ def set_coords(self, ob):
+ raise NotImplementedError, \
+ "set_coords must be provided by derived classes"
+
+ coords = property(get_coords, set_coords)
+
+ # Python feature protocol
+
+ @property
+ def __geo_interface__(self):
+ raise NotImplementedError
+
+ @property
+ def type(self):
+ return self.geometryType()
+
+ # Type of geometry and its representations
+
+ def geometryType(self):
+ """Returns a string representing the geometry type, e.g. 'Polygon'."""
+ return geometry_type_name(self._geom)
+
+ def to_wkb(self):
+ """Returns a WKB byte string representation of the geometry."""
+ size = c_int()
+ bytes = lgeos.GEOSGeomToWKB_buf(self._geom, byref(size))
+ return string_at(bytes, size.value)
+
+ def to_wkt(self):
+ """Returns a WKT string representation of the geometry."""
+ return string_at(lgeos.GEOSGeomToWKT(self._geom))
+
+ geom_type = property(geometryType)
+ wkt = property(to_wkt)
+ wkb = property(to_wkb)
+
+ # Basic geometry properties
+
+ @property
+ def area(self):
+ a = c_double()
+ retval = lgeos.GEOSArea(self._geom, byref(a))
+ return a.value
+
+ @property
+ def length(self):
+ len = c_double()
+ retval = lgeos.GEOSLength(self._geom, byref(len))
+ return len.value
+
+ def distance(self, other):
+ d = c_double()
+ retval = lgeos.GEOSDistance(self._geom, other._geom, byref(d))
+ return d.value
+
+ # Topology operations
+ #
+ # These use descriptors to reduce the amount of boilerplate.
+
+ envelope = UnaryTopologicalOp(lgeos.GEOSEnvelope, geom_factory)
+ intersection = BinaryTopologicalOp(lgeos.GEOSIntersection, geom_factory)
+ convex_hull = UnaryTopologicalOp(lgeos.GEOSConvexHull, geom_factory)
+ difference = BinaryTopologicalOp(lgeos.GEOSDifference, geom_factory)
+ symmetric_difference = BinaryTopologicalOp(lgeos.GEOSSymDifference,
+ geom_factory)
+ boundary = UnaryTopologicalOp(lgeos.GEOSBoundary, geom_factory)
+ union = BinaryTopologicalOp(lgeos.GEOSUnion, geom_factory)
+ centroid = UnaryTopologicalOp(lgeos.GEOSGetCentroid, geom_factory)
+
+ # Buffer has a unique distance argument, so not a descriptor
+ def buffer(self, distance, quadsegs=16):
+ return geom_factory(
+ lgeos.GEOSBuffer(self._geom, c_double(distance), c_int(quadsegs))
+ )
+
+ # Relate has a unique string return value
+ def relate(self, other):
+ func = lgeos.GEOSRelate
+ func.restype = c_char_p
+ return lgeos.GEOSRelate(self._geom, other._geom)
+
+ # Binary predicates
+ #
+ # These use descriptors to reduce the amount of boilerplate.
+
+ # TODO: Relate Pattern?
+ disjoint = BinaryPredicate(lgeos.GEOSDisjoint)
+ touches = BinaryPredicate(lgeos.GEOSTouches)
+ intersects = BinaryPredicate(lgeos.GEOSIntersects)
+ crosses = BinaryPredicate(lgeos.GEOSCrosses)
+ within = BinaryPredicate(lgeos.GEOSWithin)
+ contains = BinaryPredicate(lgeos.GEOSContains)
+ overlaps = BinaryPredicate(lgeos.GEOSOverlaps)
+ equals = BinaryPredicate(lgeos.GEOSEquals)
+
+ # Unary predicates
+ #
+ # These use descriptors to reduce the amount of boilerplate.
+
+ is_empty = UnaryPredicate(lgeos.GEOSisEmpty)
+ is_valid = UnaryPredicate(lgeos.GEOSisValid)
+ is_simple = UnaryPredicate(lgeos.GEOSisSimple)
+ is_ring = UnaryPredicate(lgeos.GEOSisRing)
+ has_z = UnaryPredicate(lgeos.GEOSHasZ)
+
+ @property
+ def bounds(self):
+ env = self.envelope
+ if env.geom_type != 'Polygon':
+ raise ValueError, env.wkt
+ cs = lgeos.GEOSGeom_getCoordSeq(env.exterior._geom)
+ cs_len = c_int(0)
+ lgeos.GEOSCoordSeq_getSize(cs, byref(cs_len))
+
+ minx = 1.e+20
+ maxx = -1e+20
+ miny = 1.e+20
+ maxy = -1e+20
+ temp = c_double()
+ for i in xrange(cs_len.value):
+ lgeos.GEOSCoordSeq_getX(cs, i, byref(temp))
+ x = temp.value
+ if x < minx: minx = x
+ if x > maxx: maxx = x
+ lgeos.GEOSCoordSeq_getY(cs, i, byref(temp))
+ y = temp.value
+ if y < miny: miny = y
+ if y > maxy: maxy = y
+
+ return (minx, miny, maxx, maxy)
+
Added: trunk/toolkits/basemap-testing/lib/shapely/geometry/collection.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/shapely/geometry/collection.py (rev 0)
+++ trunk/toolkits/basemap-testing/lib/shapely/geometry/collection.py 2007-11-12 16:54:58 UTC (rev 4226)
@@ -0,0 +1,24 @@
+"""
+Geometry collections.
+"""
+
+from shapely.geometry.base import BaseGeometry
+
+class GeometryCollection(BaseGeometry):
+
+ """A geometry collection.
+ """
+
+ def __init__(self):
+ BaseGeometry.__init__(self)
+
+
+# Test runner
+def _test():
+ import doctest
+ doctest.testmod()
+
+
+if __name__ == "__main__":
+ _test()
+
Added: trunk/toolkits/basemap-testing/lib/shapely/geometry/geo.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/shapely/geometry/geo.py (rev 0)
+++ trunk/toolkits/basemap-testing/lib/shapely/geometry/geo.py 2007-11-12 16:54:58 UTC (rev 4226)
@@ -0,0 +1,35 @@
+
+from point import PointAdapter
+from linestring import LineStringAdapter
+from polygon import PolygonAdapter
+from multipoint import MultiPointAdapter
+from multilinestring import MultiLineStringAdapter
+from multipolygon import MultiPolygonAdapter
+
+
+def asShape(context):
+ """Adapts the context to a geometry interface. The coordinates remain
+ stored in the context.
+ """
+ if hasattr(context, "__geo_interface__"):
+ ob = context.__geo_interface__
+ else:
+ ob = context
+
+ geom_type = ob.get("type").lower()
+
+ if geom_type == "point":
+ return PointAdapter(ob["coordinates"])
+ elif geom_type == "linestring":
+ return LineStringAdapter(ob["coordinates"])
+ elif geom_type == "polygon":
+ return PolygonAdapter(ob["coordinates"][0], ob["coordinates"][1:])
+ elif geom_type == "multipoint":
+ return MultiPointAdapter(ob["coordinates"])
+ elif geom_type == "multilinestring":
+ return MultiLineStringAdapter(ob["coordinates"])
+ elif geom_type == "multipolygon":
+ return MultiPolygonAdapter(ob["coordinates"])
+ else:
+ raise ValueError, "Unknown geometry type: %s" % geom_type
+
Added: trunk/toolkits/basemap-testing/lib/shapely/geometry/linestring.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/shapely/geometry/linestring.py (rev 0)
+++ trunk/toolkits/basemap-testing/lib/shapely/geometry/linestring.py 2007-11-12 16:54:58 UTC (rev 4226)
@@ -0,0 +1,228 @@
+"""
+Line strings.
+"""
+
+from ctypes import byref, c_double, c_int, cast, POINTER, pointer
+from ctypes import ArgumentError
+
+from shapely.geos import lgeos
+from shapely.geometry.base import BaseGeometry
+
+
+def geos_linestring_from_py(ob, update_geom=None, update_ndim=0):
+ try:
+ # From array protocol
+ array = ob.__array_interface__
+ assert len(array['shape']) == 2
+ m = array['shape'][0]
+ n = array['shape'][1]
+ assert m >= 2
+ assert n == 2 or n == 3
+
+ # Make pointer to the coordinate array
+ try:
+ cp = cast(array['data'][0], POINTER(c_double))
+ except ArgumentError:
+ cp = array['data']
+
+ # Create a coordinate sequence
+ if update_geom is not None:
+ cs = lgeos.GEOSGeom_getCoordSeq(update_geom)
+ if n != update_ndim:
+ raise ValueError, \
+ "Wrong coordinate dimensions; this geometry has dimensions: %d" \
+ % update_ndim
+ else:
+ cs = lgeos.GEOSCoordSeq_create(m, n)
+
+ # add to coordinate sequence
+ for i in xrange(m):
+ dx = c_double(cp[n*i])
+ dy = c_double(cp[n*i+1])
+ dz = None
+ if n == 3:
+ dz = c_double(cp[n*i+2])
+
+ # Because of a bug in the GEOS C API,
+ # always set X before Y
+ lgeos.GEOSCoordSeq_setX(cs, i, dx)
+ lgeos.GEOSCoordSeq_setY(cs, i, dy)
+ if n == 3:
+ lgeos.GEOSCoordSeq_setZ(cs, i, dz)
+
+ except AttributeError:
+ # Fall back on list
+ m = len(ob)
+ n = len(ob[0])
+ assert m >= 2
+ assert n == 2 or n == 3
+
+ # Create a coordinate sequence
+ if update_geom is not None:
+ cs = lgeos.GEOSGeom_getCoordSeq(update_geom)
+ if n != update_ndim:
+ raise ValueError, \
+ "Wrong coordinate dimensions; this geometry has dimensions: %d" \
+ % update_ndim
+ else:
+ cs = lgeos.GEOSCoordSeq_create(m, n)
+
+ # add to coordinate sequence
+ for i in xrange(m):
+ coords = ob[i]
+ dx = c_double(coords[0])
+ dy = c_double(coords[1])
+ dz = None
+ if n == 3:
+ dz = c_double(coords[2])
+
+ # Because of a bug in the GEOS C API,
+ # always set X before Y
+ lgeos.GEOSCoordSeq_setX(cs, i, dx)
+ lgeos.GEOSCoordSeq_setY(cs, i, dy)
+ if n == 3:
+ lgeos.GEOSCoordSeq_setZ(cs, i, dz)
+
+ if update_geom is not None:
+ return None
+ else:
+ return (lgeos.GEOSGeom_createLineString(cs), n)
+
+def update_linestring_from_py(geom, ob):
+ geos_linestring_from_py(ob, geom._geom, geom._ndim)
+
+
+class LineString(BaseGeometry):
+
+ """A line string, also known as a polyline.
+
+ """
+
+ def __init__(self, coordinates=None):
+ """Initialize.
+
+ Parameters
+ ----------
+
+ coordinates : sequence or array
+ This may be an object that satisfies the numpy array protocol,
+ providing an M x 2 or M x 3 (with z) array, or it may be a sequence
+ of x, y (,z) coordinate sequences.
+
+ Example
+ -------
+
+ >>> line = LineString([[0.0, 0.0], [1.0, 2.0]])
+ >>> line = LineString(array([[0.0, 0.0], [1.0, 2.0]]))
+
+ Each result in a line string from (0.0, 0.0) to (1.0, 2.0).
+ """
+ BaseGeometry.__init__(self)
+
+ if coordinates is None:
+ # allow creation of null lines, to support unpickling
+ pass
+ else:
+ self._geom, self._ndim = geos_linestring_from_py(coordinates)
+
+ @property
+ def __geo_interface__(self):
+ return {
+ 'type': 'LineString',
+ 'coordinates': tuple(self.coords)
+ }
+
+ @property
+ def ctypes(self):
+ if not self._ctypes_data:
+ cs = lgeos.GEOSGeom_getCoordSeq(self._geom)
+ cs_len = c_int(0)
+ lgeos.GEOSCoordSeq_getSize(cs, byref(cs_len))
+ temp = c_double()
+ n = self._ndim
+ m = cs_len.value
+ array_type = c_double * (m * n)
+ data = array_type()
+ for i in xrange(m):
+ lgeos.GEOSCoordSeq_getX(cs, i, byref(temp))
+ data[n*i] = temp.value
+ lgeos.GEOSCoordSeq_getY(cs, i, byref(temp))
+ data[n*i+1] = temp.value
+ if n == 3: # TODO: use hasz
+ lgeos.GEOSCoordSeq_getZ(cs, i, byref(temp))
+ data[n*i+2] = temp.value
+ self._ctypes_data = data
+ return self._ctypes_data
+
+ def array_interface(self):
+ """Provide the Numpy array protocol."""
+ ai = self.array_interface_base
+ ai.update({'shape': (len(self.coords), self._ndim)})
+ return ai
+ __array_interface__ = property(array_interface)
+
+ # Coordinate access
+
+ def set_coords(self, coordinates):
+ update_linestring_from_py(self, coordinates)
+
+ coords = property(BaseGeometry.get_coords, set_coords)
+
+
+class LineStringAdapter(LineString):
+
+ """Adapts a Python coordinate pair or a numpy array to the line string
+ interface.
+ """
+
+ context = None
+
+ def __init__(self, context):
+ self.context = context
+
+ # Override base class __del__
+ def __del__(self):
+ pass
+
+ @property
+ def _ndim(self):
+ try:
+ # From array protocol
+ array = self.context.__array_interface__
+ n = array['shape'][1]
+ assert n == 2 or n == 3
+ return n
+ except AttributeError:
+ # Fall back on list
+ return len(self.context[0])
+
+ @property
+ def _geom(self):
+ """Keeps the GEOS geometry in synch with the context."""
+ return geos_linestring_from_py(self.context)[0]
+
+ @property
+ def __array_interface__(self):
+ """Provide the Numpy array protocol."""
+ try:
+ return self.context.__array_interface__
+ except AttributeError:
+ return self.array_interface()
+
+ coords = property(BaseGeometry.get_coords)
+
+
+def asLineString(context):
+ """Factory for PointAdapter instances."""
+ return LineStringAdapter(context)
+
+
+# Test runner
+def _test():
+ import doctest
+ doctest.testmod()
+
+
+if __name__ == "__main__":
+ _test()
+
Added: trunk/toolkits/basemap-testing/lib/shapely/geometry/multilinestring.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/shapely/geometry/multilinestring.py (rev 0)
+++ trunk/toolkits/basemap-testing/lib/shapely/geometry/multilinestring.py 2007-11-12 16:54:58 UTC (rev 4226)
@@ -0,0 +1,155 @@
+"""
+Multi-part collection of linestrings.
+"""
+
+from ctypes import byref, c_double, c_int, c_void_p, cast, POINTER, pointer
+
+from shapely.geos import lgeos
+from shapely.geometry.base import BaseGeometry, GeometrySequence
+from shapely.geometry.linestring import LineString, geos_linestring_from_py
+
+
+def geos_multilinestring_from_py(ob):
+ """ob must be either a sequence or array of sequences or arrays."""
+ try:
+ # From array protocol
+ array = ob.__array_interface__
+ assert len(array['shape']) == 1
+ L = array['shape'][0]
+ assert L >= 1
+
+ # Make pointer to the coordinate array
+ cp = cast(array['data'][0], POINTER(c_double))
+
+ # Array of pointers to sub-geometries
+ subs = (c_void_p * L)()
+
+ for l in xrange(L):
+ geom, ndims = geos_linestring_from_py(array['data'][l])
+ subs[i] = cast(geom, c_void_p)
+ N = lgeos.GEOSGeom_getDimensions(subs[0])
+
+ except AttributeError:
+ # Fall back on list
+ L = len(ob)
+ N = len(ob[0][0])
+ assert L >= 1
+ assert N == 2 or N == 3
+
+ # Array of pointers to point geometries
+ subs = (c_void_p * L)()
+
+ # add to coordinate sequence
+ for l in xrange(L):
+ geom, ndims = geos_linestring_from_py(ob[l])
+ subs[l] = cast(geom, c_void_p)
+
+ return (lgeos.GEOSGeom_createCollection(5, subs, L), N)
+
+
+class MultiLineString(BaseGeometry):
+
+ """a multiple linestring geometry.
+ """
+
+ def __init__(self, coordinates=None):
+ """Initialize.
+
+ Parameters
+ ----------
+
+ coordinates : sequence
+ Contains coordinate sequences or objects that provide the numpy
+ array protocol, providing an M x 2 or M x 3 (with z) array.
+
+ Example
+ -------
+
+ >>> geom = MultiLineString( [[[0.0, 0.0], [1.0, 2.0]]] )
+ >>> geom = MultiLineString( [ array([[0.0, 0.0], [1.0, 2.0]]) ] )
+
+ Each result in a collection containing one line string.
+ """
+ BaseGeometry.__init__(self)
+
+ if coordinates is None:
+ # allow creation of null lines, to support unpickling
+ pass
+ else:
+ self._geom, self._ndim = geos_multilinestring_from_py(coordinates)
+
+ @property
+ def __geo_interface__(self):
+ return {
+ 'type': 'MultiLineString',
+ 'coordinates': tuple(tuple(c for c in g.coords) for g in self.geoms)
+ }
+
+ @property
+ def ctypes(self):
+ raise NotImplementedError, \
+ "Multi-part geometries have no ctypes representations"
+
+ @property
+ def __array_interface__(self):
+ """Provide the Numpy array protocol."""
+ raise NotImplementedError, \
+ "Multi-part geometries do not themselves provide the array interface"
+
+ @property
+ def coords(self):
+ raise NotImplementedError, \
+ "Multi-part geometries do not provide a coordinate sequence"
+
+ @property
+ def geoms(self):
+ return GeometrySequence(self, LineString)
+
+
+class MultiLineStringAdapter(MultiLineString):
+
+ """Adapts sequences of sequences or numpy arrays to the multilinestring
+ interface.
+ """
+
+ context = None
+
+ def __init__(self, context):
+ self.context = context
+
+ # Override base class __del__
+ def __del__(self):
+ pass
+
+ @property
+ def _ndim(self):
+ try:
+ # From array protocol
+ array = self.context[0].__array_interface__
+ n = array['shape'][1]
+ assert n == 2 or n == 3
+ return n
+ except AttributeError:
+ # Fall back on list
+ return len(self.context[0][0])
+
+ @property
+ def _geom(self):
+ """Keeps the GEOS geometry in synch with the context."""
+ return geos_multilinestring_from_py(self.context)[0]
+
+
+def asMultiLineString(context):
+ """Factory for MultiLineStringAdapter instances."""
+ return MultiLineStringAdapter(context)
+
+
+# Test runner
+def _test():
+ import doctest
+ doctest.testmod()
+
+
+if __name__ == "__main__":
+ _test()
+
Added: trunk/toolkits/basemap-testing/lib/shapely/geometry/multipoint.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/shapely/geometry/multipoint.py (rev 0)
+++ trunk/toolkits/basemap-testing/lib/shapely/geometry/multipoint.py 2007-11-12 16:54:58 UTC (rev 4226)
@@ -0,0 +1,182 @@
+"""
+Multiple points.
+"""
+
+from ctypes import byref, c_double, c_int, c_void_p, cast, POINTER, pointer
+
+from shapely.geos import lgeos
+from shapely.geometry.base import BaseGeometry, GeometrySequence
+from shapely.geometry.point import Point, geos_point_from_py
+
+
+def geos_multipoint_from_py(ob):
+ try:
+ # From array protocol
+ array = ob.__array_interface__
+ assert len(array['shape']) == 2
+ m = array['shape'][0]
+ n = array['shape'][1]
+ assert m >= 1
+ assert n == 2 or n == 3
+
+ # Make pointer to the coordinate array
+ cp = cast(array['data'][0], POINTER(c_double))
+
+ # Array of pointers to sub-geometries
+ subs = (c_void_p * m)()
+
+ for i in xrange(m):
+ geom, ndims = geos_point_from_py(cp[n*i:n*i+2])
+ subs[i] = cast(geom, c_void_p)
+
+ except AttributeError:
+ # Fall back on list
+ m = len(ob)
+ n = len(ob[0])
+ assert n == 2 or n == 3
+
+ # Array of pointers to point geometries
+ subs = (c_void_p * m)()
+
+ # add to coordinate sequence
+ for i in xrange(m):
+ coords = ob[i]
+ geom, ndims = geos_point_from_py(coords)
+ subs[i] = cast(geom, c_void_p)
+
+ return (lgeos.GEOSGeom_createCollection(4, subs, m), n)
+
+
+class MultiPoint(BaseGeometry):
+
+ """A multiple point geometry.
+ """
+
+ def __init__(self, coordinates=None):
+ """Initialize.
+
+ Parameters
+ ----------
+
+ coordinates : sequence or array
+ This may be an object that satisfies the numpy array protocol,
+ providing an M x 2 or M x 3 (with z) array, or it may be a sequence
+ of x, y (,z) coordinate sequences.
+
+ Example
+ -------
+
+ >>> geom = MultiPoint([[0.0, 0.0], [1.0, 2.0]])
+ >>> geom = MultiPoint(array([[0.0, 0.0], [1.0, 2.0]]))
+
+ Each result in a line string from (0.0, 0.0) to (1.0, 2.0).
+ """
+ BaseGeometry.__init__(self)
+
+ if coordinates is None:
+ # allow creation of null lines, to support unpickling
+ pass
+ else:
+ self._geom, self._ndim = geos_multipoint_from_py(coordinates)
+
+
+ @property
+ def __geo_interface__(self):
+ return {
+ 'type': 'MultiPoint',
+ 'coordinates': tuple([g.coords[0] for g in self.geoms])
+ }
+
+ @property
+ def ctypes(self):
+ if not self._ctypes_data:
+ temp = c_double()
+ n = self._ndim
+ m = len(self.geoms)
+ array_type = c_double * (m * n)
+ data = array_type()
+ for i in xrange(m):
+ g = self.geoms[i]._geom
+ cs = lgeos.GEOSGeom_getCoordSeq(g)
+ lgeos.GEOSCoordSeq_getX(cs, 0, byref(temp))
+ data[n*i] = temp.value
+ lgeos.GEOSCoordSeq_getY(cs, 0, byref(temp))
+ data[n*i+1] = temp.value
+ if n == 3: # TODO: use hasz
+ lgeos.GEOSCoordSeq_getZ(cs, 0, byref(temp))
+ data[n*i+2] = temp.value
+ self._ctypes_data = data
+ return self._ctypes_data
+
+ def array_interface(self):
+ """Provide the Numpy array protocol."""
+ ai = self.array_interface_base
+ ai.update({'shape': (len(self.geoms), self._ndim)})
+ return ai
+ __array_interface__ = property(array_interface)
+
+ @property
+ def coords(self):
+ raise NotImplementedError, \
+ "Multipart geometries do not themselves provide coordinate sequences"
+
+ @property
+ def geoms(self):
+ return GeometrySequence(self, Point)
+
+
+class MultiPointAdapter(MultiPoint):
+
+ """Adapts a Python coordinate pair or a numpy array to the multipoint
+ interface.
+ """
+
+ context = None
+
+ def __init__(self, context):
+ self.context = context
+
+ # Override base class __del__
+ def __del__(self):
+ pass
+
+ @property
+ def _ndim(self):
+ try:
+ # From array protocol
+ array = self.context.__array_interface__
+ n = array['shape'][1]
+ assert n == 2 or n == 3
+ return n
+ except AttributeError:
+ # Fall back on list
+ return len(self.context[0])
+
+ @property
+ def _geom(self):
+ """Keeps the GEOS geometry in synch with the context."""
+ return geos_multipoint_from_py(self.context)[0]
+
+ @property
+ def __array_interface__(self):
+ """Provide the Numpy array protocol."""
+ try:
+ return self.context.__array_interface__
+ except AttributeError:
+ return self.array_interface()
+
+
+def asMultiPoint(context):
+ """Factory for MultiPointAdapter instances."""
+ return MultiPointAdapter(context)
+
+
+# Test runner
+def _test():
+ import doctest
+ doctest.testmod()
+
+
+if __name__ == "__main__":
+ _test()
+
Added: trunk/toolkits/basemap-testing/lib/shapely/geometry/multipolygon.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/shapely/geometry/multipolygon.py (rev 0)
+++ trunk/toolkits/basemap-testing/lib/shapely/geometry/multipolygon.py 2007-11-12 16:54:58 UTC (rev 4226)
@@ -0,0 +1,141 @@
+"""
+Multi-part collection of polygons.
+"""
+
+from ctypes import byref, c_double, c_int, c_void_p, cast, POINTER, pointer
+
+from shapely.geos import lgeos
+from shapely.geometry.base import BaseGeometry, GeometrySequence
+from shapely.geometry.polygon import Polygon, geos_polygon_from_py
+
+
+def geos_multipolygon_from_py(ob):
+ """ob must be either a sequence or array of sequences or arrays."""
+ L = len(ob)
+ N = len(ob[0][0][0])
+ assert L >= 1
+ assert N == 2 or N == 3
+
+ subs = (c_void_p * L)()
+ for l in xrange(L):
+ geom, ndims = geos_polygon_from_py(ob[l][0], ob[l][1])
+ subs[l] = cast(geom, c_void_p)
+
+ return (lgeos.GEOSGeom_createCollection(6, subs, L), N)
+
+
+class MultiPolygon(BaseGeometry):
+
+ """a multiple polygon geometry.
+ """
+
+ def __init__(self, polygons=None):
+ """Initialize.
+
+ Parameters
+ ----------
+
+ polygons : sequence
+ A sequence of (shell, holes) tuples where shell is the sequence
+ representation of a linear ring (see linearring.py) and holes is
+ a sequence of such linear rings.
+
+ Example
+ -------
+ >>> geom = MultiPolygon( [
+ ... (
+ ... ((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)),
+ ... [((0.1,0.1), (0.1,0.2), (0.2,0.2), (0.2,0.1))]
+ ... )
+ ... ] )
+ """
+ BaseGeometry.__init__(self)
+
+ if polygons is None:
+ # allow creation of null collections, to support unpickling
+ pass
+ else:
+ self._geom, self._ndim = geos_multipolygon_from_py(polygons)
+
+ @property
+ def __geo_interface__(self):
+ allcoords = []
+ for geom in self.geoms:
+ coords = []
+ coords.append(tuple(geom.exterior.coords))
+ for hole in geom.interiors:
+ coords.append(tuple(hole.coords))
+ allcoords.append(coords)
+ return {
+ 'type': 'MultiPolygon',
+ 'coordinates': allcoords
+ }
+
+ @property
+ def ctypes(self):
+ raise NotImplementedError, \
+ "Multi-part geometries have no ctypes representations"
+
+ @property
+ def __array_interface__(self):
+ """Provide the Numpy array protocol."""
+ raise NotImplementedError, \
+ "Multi-part geometries do not themselves provide the array interface"
+
+ @property
+ def coords(self):
+ raise NotImplementedError, \
+ "Multi-part geometries do not provide a coordinate sequence"
+
+ @property
+ def geoms(self):
+ return GeometrySequence(self, Polygon)
+
+
+class MultiPolygonAdapter(MultiPolygon):
+
+ """Adapts sequences of sequences or numpy arrays to the multipolygon
+ interface.
+ """
+
+ context = None
+
+ def __init__(self, context):
+ self.context = context
+
+ # Override base class __del__
+ def __del__(self):
+ pass
+
+ @property
+ def _ndim(self):
+ try:
+ # From array protocol
+ array = self.context[0][0].__array_interface__
+ n = array['shape'][1]
+ assert n == 2 or n == 3
+ return n
+ except AttributeError:
+ # Fall back on list
+ return len(self.context[0][0][0])
+
+ @property
+ def _geom(self):
+ """Keeps the GEOS geometry in synch with the context."""
+ return geos_multipolygon_from_py(self.context)[0]
+
+
+def asMultiPolygon(context):
+ """Factory for MultiLineStringAdapter instances."""
+ return MultiPolygonAdapter(context)
+
+
+# Test runner
+def _test():
+ import doctest
+ doctest.testmod()
+
+
+if __name__ == "__main__":
+ _test()
+
Added: trunk/toolkits/basemap-testing/lib/shapely/geometry/point.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/shapely/geometry/point.py (rev 0)
+++ trunk/toolkits/basemap-testing/lib/shapely/geometry/point.py 2007-11-12 16:54:58 UTC (rev 4226)
@@ -0,0 +1,242 @@
+"""
+Points.
+"""
+
+from ctypes import string_at, create_string_buffer, \
+ c_char_p, c_double, c_float, c_int, c_uint, c_size_t, c_ubyte, \
+ c_void_p, byref
+from ctypes import cast, POINTER
+
+from shapely.geos import lgeos, DimensionError
+from shapely.geometry.base import BaseGeometry, CoordinateSequence
+
+
+def geos_point_from_py(ob, update_geom=None, update_ndim=0):
+ """Create a GEOS geom from an object that is a coordinate sequence
+ or that provides the array interface.
+
+ Returns the GEOS geometry and the number of its dimensions.
+ """
+ try:
+ # From array protocol
+ array = ob.__array_interface__
+ assert len(array['shape']) == 1
+ n = array['shape'][0]
+ assert n == 2 or n == 3
+
+ cdata = array['data'][0]
+ cp = cast(cdata, POINTER(c_double))
+ dx = c_double(cp[0])
+ dy = c_double(cp[1])
+ dz = None
+ if n == 3:
+ dz = c_double(cp[2])
+ ndim = 3
+ except AttributeError:
+ # Fall back on list
+ coords = ob
+ n = len(coords)
+ dx = c_double(coords[0])
+ dy = c_double(coords[1])
+ dz = None
+ if n == 3:
+ dz = c_double(coords[2])
+
+ if update_geom:
+ cs = lgeos.GEOSGeom_getCoordSeq(update_geom)
+ if n != update_ndim:
+ raise ValueError, \
+ "Wrong coordinate dimensions; this geometry has dimensions: %d" \
+ % update_ndim
+ else:
+ cs = lgeos.GEOSCoordSeq_create(1, n)
+
+ # Because of a bug in the GEOS C API, always set X before Y
+ lgeos.GEOSCoordSeq_setX(cs, 0, dx)
+ lgeos.GEOSCoordSeq_setY(cs, 0, dy)
+ if n == 3:
+ lgeos.GEOSCoordSeq_setZ(cs, 0, dz)
+
+ if update_geom:
+ return None
+ else:
+ return (lgeos.GEOSGeom_createPoint(cs), n)
+
+def update_point_from_py(geom, ob):
+ geos_point_from_py(ob, geom._geom, geom._ndim)
+
+
+class Point(BaseGeometry):
+
+ """A point geometry.
+
+ Attributes
+ ----------
+ x, y, z : float
+ Coordinate values
+
+ Example
+ -------
+ >>> p = Point(1.0, -1.0)
+ >>> str(p)
+ 'POINT (1.0000000000000000 -1.0000000000000000)'
+ >>> p.y = 0.0
+ >>> p.y
+ 0.0
+ >>> p.x
+ 1.0
+ >>> p.array
+ [[1.0, 0.0]]
+ """
+
+ def __init__(self, *args):
+ """This *copies* the given data to a new GEOS geometry.
+
+ Parameters
+ ----------
+
+ There are 2 cases:
+
+ 1) 1 parameter: this must satisfy the numpy array protocol.
+ 2) 2 or more parameters: x, y, z : float
+ Easting, northing, and elevation.
+ """
+ BaseGeometry.__init__(self)
+
+ if len(args) == 0:
+ # allow creation of null points, to support unpickling
+ pass
+ else:
+ if len(args) == 1:
+ self._geom, self._ndim = geos_point_from_py(args[0])
+ else:
+ self._geom, self._ndim = geos_point_from_py(tuple(args))
+
+ # Coordinate getters and setters
+
+ @property
+ def x(self):
+ """Return x coordinate."""
+ cs = lgeos.GEOSGeom_getCoordSeq(self._geom)
+ d = c_double()
+ lgeos.GEOSCoordSeq_getX(cs, 0, byref(d))
+ return d.value
+
+ @property
+ def y(self):
+ """Return y coordinate."""
+ cs = lgeos.GEOSGeom_getCoordSeq(self._geom)
+ d = c_double()
+ lgeos.GEOSCoordSeq_getY(cs, 0, byref(d))
+ return d.value
+
+ @property
+ def z(self):
+ """Return z coordinate."""
+ if self._ndim != 3:
+ raise DimensionError, "This point has no z coordinate."
+ cs = lgeos.GEOSGeom_getCoordSeq(self._geom)
+ d = c_double()
+ lgeos.GEOSCoordSeq_getZ(cs, 0, byref(d))
+ return d.value
+
+ @property
+ def __geo_interface__(self):
+ return {
+ 'type': 'Point',
+ 'coordinates': self.coords[0]
+ }
+
+ @property
+ def ctypes(self):
+ if not self._ctypes_data:
+ array_type = c_double * self._ndim
+ array = array_type()
+ array[0] = self.x
+ array[1] = self.y
+ if self._ndim == 3:
+ array[2] = self.z
+ self._ctypes_data = array
+ return self._ctypes_data
+
+ def array_interface(self):
+ """Provide the Numpy array protocol."""
+ ai = self.array_interface_base
+ ai.update({'shape': (self._ndim,)})
+ return ai
+ __array_interface__ = property(array_interface)
+
+ @property
+ def bounds(self):
+ cs = lgeos.GEOSGeom_getCoordSeq(self._geom)
+ x = c_double()
+ y = c_double()
+ lgeos.GEOSCoordSeq_getX(cs, 0, byref(x))
+ lgeos.GEOSCoordSeq_getY(cs, 0, byref(y))
+ return (x.value, y.value, x.value, y.value)
+
+ # Coordinate access
+
+ def set_coords(self, coordinates):
+ update_point_from_py(self, coordinates)
+
+ coords = property(BaseGeometry.get_coords, set_coords)
+
+
+class PointAdapter(Point):
+
+ """Adapts a Python coordinate pair or a numpy array to the point interface.
+ """
+
+ context = None
+
+ def __init__(self, context):
+ self.context = context
+
+ # Override base class __del__
+ def __del__(self):
+ pass
+
+ @property
+ def _ndim(self):
+ try:
+ # From array protocol
+ array = self.context.__array_interface__
+ n = array['shape'][0]
+ assert n == 2 or n == 3
+ return n
+ except AttributeError:
+ # Fall back on list
+ return len(self.context)
+
+ @property
+ def _geom(self):
+ """Keeps the GEOS geometry in synch with the context."""
+ return geos_point_from_py(self.context)[0]
+
+ # TODO: reimplement x, y, z properties without calling invoking _geom
+
+ @property
+ def __array_interface__(self):
+ """Provide the Numpy array protocol."""
+ try:
+ return self.context.__array_interface__
+ except AttributeError:
+ return self.array_interface()
+
+ coords = property(BaseGeometry.get_coords)
+
+
+def asPoint(context):
+ """Factory for PointAdapter instances."""
+ return PointAdapter(context)
+
+
+# Test runner
+def _test():
+ import doctest
+ doctest.testmod()
+
+if __name__ == "__main__":
+ _test()
+
Added: trunk/toolkits/basemap-testing/lib/shapely/geometry/polygon.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/shapely/geometry/polygon.py (rev 0)
+++ trunk/toolkits/basemap-testing/lib/shapely/geometry/polygon.py 2007-11-12 16:54:58 UTC (rev 4226)
@@ -0,0 +1,414 @@
+"""
+Polygons and their linear ring components.
+"""
+
+from ctypes import byref, c_double, c_int, c_void_p, cast, POINTER, pointer
+
+from shapely.geos import lgeos
+from shapely.geometry.base import BaseGeometry
+from shapely.geometry.linestring import LineString, LineStringAdapter
+
+
+def geos_linearring_from_py(ob, update_geom=None, update_ndim=0):
+ try:
+ # From array protocol
+ array = ob.__array_interface__
+ assert len(array['shape']) == 2
+ m = array['shape'][0]
+ n = array['shape'][1]
+ assert m >= 2
+ assert n == 2 or n == 3
+
+ # Make pointer to the coordinate array
+ try:
+ cp = cast(array['data'][0], POINTER(c_double))
+ except ArgumentError:
+ cp = array['data']
+
+ # Add closing coordinates to sequence?
+ if cp[0] != cp[m*n-n] or cp[1] != cp[m*n-n+1]:
+ M = m + 1
+ else:
+ M = m
+
+ # Create a coordinate sequence
+ if update_geom is not None:
+ cs = lgeos.GEOSGeom_getCoordSeq(update_geom)
+ if n != update_ndim:
+ raise ValueError, \
+ "Wrong coordinate dimensions; this geometry has dimensions: %d" \
+ % update_ndim
+ else:
+ cs = lgeos.GEOSCoordSeq_create(M, n)
+
+ # add to coordinate sequence
+ for i in xrange(m):
+ dx = c_double(cp[n*i])
+ dy = c_double(cp[n*i+1])
+ dz = None
+ if n == 3:
+ dz = c_double(cp[n*i+2])
+
+ # Because of a bug in the GEOS C API,
+ # always set X before Y
+ lgeos.GEOSCoordSeq_setX(cs, i, dx)
+ lgeos.GEOSCoordSeq_setY(cs, i, dy)
+ if n == 3:
+ lgeos.GEOSCoordSeq_setZ(cs, i, dz)
+
+ # Add closing coordinates to sequence?
+ if M > m:
+ dx = c_double(cp[0])
+ dy = c_double(cp[1])
+ dz = None
+ if n == 3:
+ dz = c_double(cp[2])
+
+ # Because of a bug in the GEOS C API,
+ # always set X before Y
+ lgeos.GEOSCoordSeq_setX(cs, M-1, dx)
+ lgeos.GEOSCoordSeq_setY(cs, M-1, dy)
+ if n == 3:
+ lgeos.GEOSCoordSeq_setZ(cs, M-1, dz)
+
+ except AttributeError:
+ # Fall back on list
+ m = len(ob)
+ n = len(ob[0])
+ assert m >= 2
+ assert n == 2 or n == 3
+
+ # Add closing coordinates if not provided
+ if ob[0][0] != ob[-1][0] or ob[0][1] != ob[-1][1]:
+ M = m + 1
+ else:
+ M = m
+
+ # Create a coordinate sequence
+ if update_geom is not None:
+ cs = lgeos.GEOSGeom_getCoordSeq(update_geom)
+ if n != update_ndim:
+ raise ValueError, \
+ "Wrong coordinate dimensions; this geometry has dimensions: %d" \
+ % update_ndim
+ else:
+ cs = lgeos.GEOSCoordSeq_create(M, n)
+
+ # add to coordinate sequence
+ for i in xrange(m):
+ coords = ob[i]
+ dx = c_double(coords[0])
+ dy = c_double(coords[1])
+ dz = None
+ if n == 3:
+ dz = c_double(coords[2])
+
+ # Because of a bug in the GEOS C API,
+ # always set X before Y
+ lgeos.GEOSCoordSeq_setX(cs, i, dx)
+ lgeos.GEOSCoordSeq_setY(cs, i, dy)
+ if n == 3:
+ lgeos.GEOSCoordSeq_setZ(cs, i, dz)
+
+ # Add closing coordinates to sequence?
+ if M > m:
+ coords = ob[0]
+ dx = c_double(coords[0])
+ dy = c_double(coords[1])
+ dz = None
+ if n == 3:
+ dz = c_double(coords[2])
+
+ # Because of a bug in the GEOS C API,
+ # always set X before Y
+ lgeos.GEOSCoordSeq_setX(cs, M-1, dx)
+ lgeos.GEOSCoordSeq_setY(cs, M-1, dy)
+ if n == 3:
+ lgeos.GEOSCoordSeq_setZ(cs, M-1, dz)
+
+ if update_geom is not None:
+ return None
+ else:
+ return (lgeos.GEOSGeom_createLinearRing(cs), n)
+
+def update_linearring_from_py(geom, ob):
+ geos_linearring_from_py(ob, geom._geom, geom._ndim)
+
+
+class LinearRing(LineString):
+
+ """A linear ring.
+ """
+
+ _ndim = 2
+
+ def __init__(self, coordinates=None):
+ """Initialize.
+
+ Parameters
+ ----------
+ coordinates : sequence or array
+ This may be an object that satisfies the numpy array protocol,
+ providing an M x 2 or M x 3 (with z) array, or it may be a sequence
+ of x, y (,z) coordinate sequences.
+
+ Rings are implicitly closed. There is no need to specific a final
+ coordinate pair identical to the first.
+
+ Example
+ -------
+ >>> ring = LinearRing( ((0.,0.), (0.,1.), (1.,1.), (1.,0.)) )
+
+ Produces a 1x1 square.
+ """
+ BaseGeometry.__init__(self)
+
+ if coordinates is None:
+ # allow creation of null lines, to support unpickling
+ pass
+ else:
+ self._geom, self._ndims = geos_linearring_from_py(coordinates)
+
+ @property
+ def __geo_interface__(self):
+ return {
+ 'type': 'LinearRing',
+ 'coordinates': tuple(self.coords)
+ }
+
+ # Coordinate access
+
+ def set_coords(self, coordinates):
+ update_linearring_from_py(self, coordinates)
+
+ coords = property(BaseGeometry.get_coords, set_coords)
+
+
+class LinearRingAdapter(LineStringAdapter):
+
+ @property
+ def _geom(self):
+ """Keeps the GEOS geometry in synch with the context."""
+ return geos_linearring_from_py(self.context)[0]
+
+ @property
+ def __geo_interface__(self):
+ return {
+ 'type': 'LinearRing',
+ 'coordinates': tuple(self.coords)
+ }
+
+ coords = property(BaseGeometry.get_coords)
+
+
+def asLinearRing(context):
+ return LinearRingAdapter(context)
+
+
+class InteriorRingSequence(object):
+
+ _factory = None
+ _geom = None
+ _ndim = None
+ _index = 0
+ _length = 0
+
+ def __init__(self, geom):
+ self._geom = geom._geom
+ self._ndim = geom._ndim
+
+ def __iter__(self):
+ self._index = 0
+ self._length = self.__len__()
+ return self
+
+ def next(self):
+ if self._index < self._length:
+ g = LinearRing()
+ g._owned = True
+ g._geom = lgeos.GEOSGetInteriorRingN(self._geom, self._index)
+ self._index += 1
+ return g
+ else:
+ raise StopIteration
+
+ def __len__(self):
+ return lgeos.GEOSGetNumInteriorRings(self._geom)
+
+ def __getitem__(self, i):
+ M = self.__len__()
+ if i + M < 0 or i >= M:
+ raise IndexError, "index out of range"
+ if i < 0:
+ ii = M + i
+ else:
+ ii = i
+ g = LinearRing()
+ g._owned = True
+ g._geom = lgeos.GEOSGetInteriorRingN(self._geom, ii)
+ return g
+
+ @property
+ def _longest(self):
+ max = 0
+ for g in iter(self):
+ l = len(g.coords)
+ if l > max:
+ max = l
+
+
+def geos_polygon_from_py(shell, holes=None):
+ if shell is not None:
+ geos_shell, ndims = geos_linearring_from_py(shell)
+ ## Polygon geometry takes ownership of the ring
+ #self._exterior._owned = True
+
+ if holes:
+ ob = holes
+ L = len(ob)
+ try:
+ N = len(ob[0][0])
+ except:
+ import pdb; pdb.set_trace()
+ assert L >= 1
+ assert N == 2 or N == 3
+
+ # Array of pointers to ring geometries
+ geos_holes = (c_void_p * L)()
+
+ # add to coordinate sequence
+ for l in xrange(L):
+ geom, ndims = geos_linearring_from_py(ob[l])
+ geos_holes[l] = cast(geom, c_void_p)
+
+ else:
+ geos_holes = POINTER(c_void_p)()
+ L = 0
+
+ return (
+ lgeos.GEOSGeom_createPolygon(
+ c_void_p(geos_shell),
+ geos_holes,
+ L
+ ),
+ ndims
+ )
+
+class Polygon(BaseGeometry):
+
+ """A line string, also known as a polyline.
+ """
+
+ _exterior = None
+ _interiors = []
+ _ndim = 2
+
+ def __init__(self, shell=None, holes=None):
+ """Initialize.
+
+ Parameters
+ ----------
+ exterior : sequence or array
+ This may be an object that satisfies the numpy array protocol,
+ providing an M x 2 or M x 3 (with z) array, or it may be a sequence
+ of x, y (,z) coordinate sequences.
+
+ Example
+ -------
+ >>> coords = ((0., 0.), (0., 1.), (1., 1.), (1., 0.), (0., 0.))
+ >>> polygon = Polygon(coords)
+ """
+ BaseGeometry.__init__(self)
+
+ if shell is not None:
+ self._geom, self._ndims = geos_polygon_from_py(shell, holes)
+
+ @property
+ def exterior(self):
+ if self._exterior is None:
+ # A polygon created from the abstract factory will have a null
+ # _exterior attribute.
+ ring = lgeos.GEOSGetExteriorRing(self._geom)
+ self._exterior = LinearRing()
+ self._exterior._geom = ring
+ self._exterior._owned = True
+ return self._exterior
+
+ @property
+ def interiors(self):
+ return InteriorRingSequence(self)
+
+ @property
+ def ctypes(self):
+ if not self._ctypes_data:
+ self._ctypes_data = self.exterior.ctypes
+ return self._ctypes_data
+
+ @property
+ def __array_interface__(self):
+ raise NotImplementedError, \
+ "A polygon does not itself provide the array interface. Its rings do."
+
+ @property
+ def coords(self):
+ raise NotImplementedError, \
+ "Component rings have coordinate sequences, but the polygon does not"
+
+ @property
+ def __geo_interface__(self):
+ coords = [tuple(self.exterior.coords)]
+ for hole in self.interiors:
+ coords.append(tuple(hole.coords))
+ return {
+ 'type': 'Polygon',
+ 'coordinates': tuple(coords)
+ }
+
+
+class PolygonAdapter(Polygon):
+
+ """Adapts sequences of sequences or numpy arrays to the polygon
+ interface.
+ """
+
+ context = None
+
+ def __init__(self, shell, holes=None):
+ self.shell = shell
+ self.holes = holes
+
+ # Override base class __del__
+ def __del__(self):
+ pass
+
+ @property
+ def _ndim(self):
+ try:
+ # From array protocol
+ array = self.shell.__array_interface__
+ n = array['shape'][1]
+ assert n == 2 or n == 3
+ return n
+ except AttributeError:
+ # Fall back on list
+ return len(self.shell[0])
+
+ @property
+ def _geom(self):
+ """Keeps the GEOS geometry in synch with the context."""
+ return geos_polygon_from_py(self.shell, self.holes)[0]
+
+
+def asPolygon(shell, holes):
+ """Factory for PolygonAdapter instances."""
+ return PolygonAdapter(shell, holes)
+
+
+# Test runner
+def _test():
+ import doctest
+ doctest.testmod()
+
+if __name__ == "__main__":
+ _test()
+
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
