diff options
Diffstat (limited to 'webapp/django/contrib/gis/gdal/geometries.py')
| -rw-r--r-- | webapp/django/contrib/gis/gdal/geometries.py | 643 |
1 files changed, 643 insertions, 0 deletions
diff --git a/webapp/django/contrib/gis/gdal/geometries.py b/webapp/django/contrib/gis/gdal/geometries.py new file mode 100644 index 0000000000..ee0ec17559 --- /dev/null +++ b/webapp/django/contrib/gis/gdal/geometries.py @@ -0,0 +1,643 @@ +""" + The OGRGeometry is a wrapper for using the OGR Geometry class + (see http://www.gdal.org/ogr/classOGRGeometry.html). OGRGeometry + may be instantiated when reading geometries from OGR Data Sources + (e.g. SHP files), or when given OGC WKT (a string). + + While the 'full' API is not present yet, the API is "pythonic" unlike + the traditional and "next-generation" OGR Python bindings. One major + advantage OGR Geometries have over their GEOS counterparts is support + for spatial reference systems and their transformation. + + Example: + >>> from django.contrib.gis.gdal import OGRGeometry, OGRGeomType, SpatialReference + >>> wkt1, wkt2 = 'POINT(-90 30)', 'POLYGON((0 0, 5 0, 5 5, 0 5)' + >>> pnt = OGRGeometry(wkt1) + >>> print pnt + POINT (-90 30) + >>> mpnt = OGRGeometry(OGRGeomType('MultiPoint'), SpatialReference('WGS84')) + >>> mpnt.add(wkt1) + >>> mpnt.add(wkt1) + >>> print mpnt + MULTIPOINT (-90 30,-90 30) + >>> print mpnt.srs.name + WGS 84 + >>> print mpnt.srs.proj + +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs + >>> mpnt.transform_to(SpatialReference('NAD27')) + >>> print mpnt.proj + +proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs + >>> print mpnt + MULTIPOINT (-89.999930378602485 29.999797886557641,-89.999930378602485 29.999797886557641) + + The OGRGeomType class is to make it easy to specify an OGR geometry type: + >>> from django.contrib.gis.gdal import OGRGeomType + >>> gt1 = OGRGeomType(3) # Using an integer for the type + >>> gt2 = OGRGeomType('Polygon') # Using a string + >>> gt3 = OGRGeomType('POLYGON') # It's case-insensitive + >>> print gt1 == 3, gt1 == 'Polygon' # Equivalence works w/non-OGRGeomType objects + True +""" +# Python library requisites. +import re, sys +from binascii import a2b_hex +from ctypes import byref, string_at, c_char_p, c_double, c_ubyte, c_void_p +from types import UnicodeType + +# Getting GDAL prerequisites +from django.contrib.gis.gdal.envelope import Envelope, OGREnvelope +from django.contrib.gis.gdal.error import OGRException, OGRIndexError, SRSException +from django.contrib.gis.gdal.geomtype import OGRGeomType +from django.contrib.gis.gdal.srs import SpatialReference, CoordTransform + +# Getting the ctypes prototype functions that interface w/the GDAL C library. +from django.contrib.gis.gdal.prototypes.geom import * +from django.contrib.gis.gdal.prototypes.srs import clone_srs + +# For more information, see the OGR C API source code: +# http://www.gdal.org/ogr/ogr__api_8h.html +# +# The OGR_G_* routines are relevant here. + +# Regular expressions for recognizing HEXEWKB and WKT. +hex_regex = re.compile(r'^[0-9A-F]+$', re.I) +wkt_regex = re.compile(r'^(?P<type>POINT|LINESTRING|LINEARRING|POLYGON|MULTIPOINT|MULTILINESTRING|MULTIPOLYGON|GEOMETRYCOLLECTION)[ACEGIMLONPSRUTY\d,\.\-\(\) ]+$', re.I) +json_regex = re.compile(r'^\{[\s\w,\-\.\"\'\:\[\]]+\}$') + +#### OGRGeometry Class #### +class OGRGeometry(object): + "Generally encapsulates an OGR geometry." + + def __init__(self, geom_input, srs=None): + "Initializes Geometry on either WKT or an OGR pointer as input." + + self._ptr = c_void_p(None) # Initially NULL + str_instance = isinstance(geom_input, basestring) + + # If HEX, unpack input to to a binary buffer. + if str_instance and hex_regex.match(geom_input): + geom_input = buffer(a2b_hex(geom_input.upper())) + str_instance = False + + # Constructing the geometry, + if str_instance: + # Checking if unicode + if isinstance(geom_input, UnicodeType): + # Encoding to ASCII, WKT or HEX doesn't need any more. + geo_input = geo_input.encode('ascii') + + wkt_m = wkt_regex.match(geom_input) + json_m = json_regex.match(geom_input) + if wkt_m: + if wkt_m.group('type').upper() == 'LINEARRING': + # OGR_G_CreateFromWkt doesn't work with LINEARRING WKT. + # See http://trac.osgeo.org/gdal/ticket/1992. + g = create_geom(OGRGeomType(wkt_m.group('type')).num) + import_wkt(g, byref(c_char_p(geom_input))) + else: + g = from_wkt(byref(c_char_p(geom_input)), None, byref(c_void_p())) + elif json_m: + if GEOJSON: + g = from_json(geom_input) + else: + raise NotImplementedError('GeoJSON input only supported on GDAL 1.5+.') + else: + # Seeing if the input is a valid short-hand string + # (e.g., 'Point', 'POLYGON'). + ogr_t = OGRGeomType(geom_input) + g = create_geom(OGRGeomType(geom_input).num) + elif isinstance(geom_input, buffer): + # WKB was passed in + g = from_wkb(str(geom_input), None, byref(c_void_p()), len(geom_input)) + elif isinstance(geom_input, OGRGeomType): + # OGRGeomType was passed in, an empty geometry will be created. + g = create_geom(geom_input.num) + elif isinstance(geom_input, c_void_p): + # OGR pointer (c_void_p) was the input. + g = geom_input + else: + raise OGRException('Invalid input type for OGR Geometry construction: %s' % type(geom_input)) + + # Now checking the Geometry pointer before finishing initialization + # by setting the pointer for the object. + if not g: + raise OGRException('Cannot create OGR Geometry from input: %s' % str(geom_input)) + self._ptr = g + + # Assigning the SpatialReference object to the geometry, if valid. + if bool(srs): self.srs = srs + + # Setting the class depending upon the OGR Geometry Type + self.__class__ = GEO_CLASSES[self.geom_type.num] + + def __del__(self): + "Deletes this Geometry." + if self._ptr: destroy_geom(self._ptr) + + ### Geometry set-like operations ### + # g = g1 | g2 + def __or__(self, other): + "Returns the union of the two geometries." + return self.union(other) + + # g = g1 & g2 + def __and__(self, other): + "Returns the intersection of this Geometry and the other." + return self.intersection(other) + + # g = g1 - g2 + def __sub__(self, other): + "Return the difference this Geometry and the other." + return self.difference(other) + + # g = g1 ^ g2 + def __xor__(self, other): + "Return the symmetric difference of this Geometry and the other." + return self.sym_difference(other) + + def __eq__(self, other): + "Is this Geometry equal to the other?" + return self.equals(other) + + def __ne__(self, other): + "Tests for inequality." + return not self.equals(other) + + def __str__(self): + "WKT is used for the string representation." + return self.wkt + + #### Geometry Properties #### + @property + def dimension(self): + "Returns 0 for points, 1 for lines, and 2 for surfaces." + return get_dims(self._ptr) + + @property + def coord_dim(self): + "Returns the coordinate dimension of the Geometry." + return get_coord_dims(self._ptr) + + @property + def geom_count(self): + "The number of elements in this Geometry." + return get_geom_count(self._ptr) + + @property + def point_count(self): + "Returns the number of Points in this Geometry." + return get_point_count(self._ptr) + + @property + def num_points(self): + "Alias for `point_count` (same name method in GEOS API.)" + return self.point_count + + @property + def num_coords(self): + "Alais for `point_count`." + return self.point_count + + @property + def geom_type(self): + "Returns the Type for this Geometry." + try: + return OGRGeomType(get_geom_type(self._ptr)) + except OGRException: + # VRT datasources return an invalid geometry type + # number, but a valid name -- we'll try that instead. + # See: http://trac.osgeo.org/gdal/ticket/2491 + return OGRGeomType(get_geom_name(self._ptr)) + + @property + def geom_name(self): + "Returns the Name of this Geometry." + return get_geom_name(self._ptr) + + @property + def area(self): + "Returns the area for a LinearRing, Polygon, or MultiPolygon; 0 otherwise." + return get_area(self._ptr) + + @property + def envelope(self): + "Returns the envelope for this Geometry." + # TODO: Fix Envelope() for Point geometries. + return Envelope(get_envelope(self._ptr, byref(OGREnvelope()))) + + @property + def extent(self): + "Returns the envelope as a 4-tuple, instead of as an Envelope object." + return self.envelope.tuple + + #### SpatialReference-related Properties #### + + # The SRS property + def get_srs(self): + "Returns the Spatial Reference for this Geometry." + try: + srs_ptr = get_geom_srs(self._ptr) + return SpatialReference(clone_srs(srs_ptr)) + except SRSException: + return None + + def set_srs(self, srs): + "Sets the SpatialReference for this geometry." + if isinstance(srs, SpatialReference): + srs_ptr = clone_srs(srs._ptr) + elif isinstance(srs, (int, long, basestring)): + sr = SpatialReference(srs) + srs_ptr = clone_srs(sr._ptr) + else: + raise TypeError('Cannot assign spatial reference with object of type: %s' % type(srs)) + assign_srs(self._ptr, srs_ptr) + + srs = property(get_srs, set_srs) + + # The SRID property + def get_srid(self): + if self.srs: return self.srs.srid + else: return None + + def set_srid(self, srid): + if isinstance(srid, (int, long)): + self.srs = srid + else: + raise TypeError('SRID must be set with an integer.') + + srid = property(get_srid, set_srid) + + #### Output Methods #### + @property + def geos(self): + "Returns a GEOSGeometry object from this OGRGeometry." + from django.contrib.gis.geos import GEOSGeometry + return GEOSGeometry(self.wkb, self.srid) + + @property + def gml(self): + "Returns the GML representation of the Geometry." + return to_gml(self._ptr) + + @property + def hex(self): + "Returns the hexadecimal representation of the WKB (a string)." + return str(self.wkb).encode('hex').upper() + #return b2a_hex(self.wkb).upper() + + @property + def json(self): + if GEOJSON: + return to_json(self._ptr) + else: + raise NotImplementedError('GeoJSON output only supported on GDAL 1.5+.') + geojson = json + + @property + def wkb_size(self): + "Returns the size of the WKB buffer." + return get_wkbsize(self._ptr) + + @property + def wkb(self): + "Returns the WKB representation of the Geometry." + if sys.byteorder == 'little': + byteorder = 1 # wkbNDR (from ogr_core.h) + else: + byteorder = 0 # wkbXDR + sz = self.wkb_size + # Creating the unsigned character buffer, and passing it in by reference. + buf = (c_ubyte * sz)() + wkb = to_wkb(self._ptr, byteorder, byref(buf)) + # Returning a buffer of the string at the pointer. + return buffer(string_at(buf, sz)) + + @property + def wkt(self): + "Returns the WKT representation of the Geometry." + return to_wkt(self._ptr, byref(c_char_p())) + + #### Geometry Methods #### + def clone(self): + "Clones this OGR Geometry." + return OGRGeometry(clone_geom(self._ptr), self.srs) + + def close_rings(self): + """ + If there are any rings within this geometry that have not been + closed, this routine will do so by adding the starting point at the + end. + """ + # Closing the open rings. + geom_close_rings(self._ptr) + + def transform(self, coord_trans, clone=False): + """ + Transforms this geometry to a different spatial reference system. + May take a CoordTransform object, a SpatialReference object, string + WKT or PROJ.4, and/or an integer SRID. By default nothing is returned + and the geometry is transformed in-place. However, if the `clone` + keyword is set, then a transformed clone of this geometry will be + returned. + """ + if clone: + klone = self.clone() + klone.transform(coord_trans) + return klone + if isinstance(coord_trans, CoordTransform): + geom_transform(self._ptr, coord_trans._ptr) + elif isinstance(coord_trans, SpatialReference): + geom_transform_to(self._ptr, coord_trans._ptr) + elif isinstance(coord_trans, (int, long, basestring)): + sr = SpatialReference(coord_trans) + geom_transform_to(self._ptr, sr._ptr) + else: + raise TypeError('Transform only accepts CoordTransform, SpatialReference, string, and integer objects.') + + def transform_to(self, srs): + "For backwards-compatibility." + self.transform(srs) + + #### Topology Methods #### + def _topology(self, func, other): + """A generalized function for topology operations, takes a GDAL function and + the other geometry to perform the operation on.""" + if not isinstance(other, OGRGeometry): + raise TypeError('Must use another OGRGeometry object for topology operations!') + + # Returning the output of the given function with the other geometry's + # pointer. + return func(self._ptr, other._ptr) + + def intersects(self, other): + "Returns True if this geometry intersects with the other." + return self._topology(ogr_intersects, other) + + def equals(self, other): + "Returns True if this geometry is equivalent to the other." + return self._topology(ogr_equals, other) + + def disjoint(self, other): + "Returns True if this geometry and the other are spatially disjoint." + return self._topology(ogr_disjoint, other) + + def touches(self, other): + "Returns True if this geometry touches the other." + return self._topology(ogr_touches, other) + + def crosses(self, other): + "Returns True if this geometry crosses the other." + return self._topology(ogr_crosses, other) + + def within(self, other): + "Returns True if this geometry is within the other." + return self._topology(ogr_within, other) + + def contains(self, other): + "Returns True if this geometry contains the other." + return self._topology(ogr_contains, other) + + def overlaps(self, other): + "Returns True if this geometry overlaps the other." + return self._topology(ogr_overlaps, other) + + #### Geometry-generation Methods #### + def _geomgen(self, gen_func, other=None): + "A helper routine for the OGR routines that generate geometries." + if isinstance(other, OGRGeometry): + return OGRGeometry(gen_func(self._ptr, other._ptr), self.srs) + else: + return OGRGeometry(gen_func(self._ptr), self.srs) + + @property + def boundary(self): + "Returns the boundary of this geometry." + return self._geomgen(get_boundary) + + @property + def convex_hull(self): + """ + Returns the smallest convex Polygon that contains all the points in + this Geometry. + """ + return self._geomgen(geom_convex_hull) + + def difference(self, other): + """ + Returns a new geometry consisting of the region which is the difference + of this geometry and the other. + """ + return self._geomgen(geom_diff, other) + + def intersection(self, other): + """ + Returns a new geometry consisting of the region of intersection of this + geometry and the other. + """ + return self._geomgen(geom_intersection, other) + + def sym_difference(self, other): + """ + Returns a new geometry which is the symmetric difference of this + geometry and the other. + """ + return self._geomgen(geom_sym_diff, other) + + def union(self, other): + """ + Returns a new geometry consisting of the region which is the union of + this geometry and the other. + """ + return self._geomgen(geom_union, other) + +# The subclasses for OGR Geometry. +class Point(OGRGeometry): + + @property + def x(self): + "Returns the X coordinate for this Point." + return getx(self._ptr, 0) + + @property + def y(self): + "Returns the Y coordinate for this Point." + return gety(self._ptr, 0) + + @property + def z(self): + "Returns the Z coordinate for this Point." + if self.coord_dim == 3: + return getz(self._ptr, 0) + + @property + def tuple(self): + "Returns the tuple of this point." + if self.coord_dim == 2: + return (self.x, self.y) + elif self.coord_dim == 3: + return (self.x, self.y, self.z) + coords = tuple + +class LineString(OGRGeometry): + + def __getitem__(self, index): + "Returns the Point at the given index." + if index >= 0 and index < self.point_count: + x, y, z = c_double(), c_double(), c_double() + get_point(self._ptr, index, byref(x), byref(y), byref(z)) + dim = self.coord_dim + if dim == 1: + return (x.value,) + elif dim == 2: + return (x.value, y.value) + elif dim == 3: + return (x.value, y.value, z.value) + else: + raise OGRIndexError('index out of range: %s' % str(index)) + + def __iter__(self): + "Iterates over each point in the LineString." + for i in xrange(self.point_count): + yield self[i] + + def __len__(self): + "The length returns the number of points in the LineString." + return self.point_count + + @property + def tuple(self): + "Returns the tuple representation of this LineString." + return tuple([self[i] for i in xrange(len(self))]) + coords = tuple + + def _listarr(self, func): + """ + Internal routine that returns a sequence (list) corresponding with + the given function. + """ + return [func(self._ptr, i) for i in xrange(len(self))] + + @property + def x(self): + "Returns the X coordinates in a list." + return self._listarr(getx) + + @property + def y(self): + "Returns the Y coordinates in a list." + return self._listarr(gety) + + @property + def z(self): + "Returns the Z coordinates in a list." + if self.coord_dim == 3: + return self._listarr(getz) + +# LinearRings are used in Polygons. +class LinearRing(LineString): pass + +class Polygon(OGRGeometry): + + def __len__(self): + "The number of interior rings in this Polygon." + return self.geom_count + + def __iter__(self): + "Iterates through each ring in the Polygon." + for i in xrange(self.geom_count): + yield self[i] + + def __getitem__(self, index): + "Gets the ring at the specified index." + if index < 0 or index >= self.geom_count: + raise OGRIndexError('index out of range: %s' % index) + else: + return OGRGeometry(clone_geom(get_geom_ref(self._ptr, index)), self.srs) + + # Polygon Properties + @property + def shell(self): + "Returns the shell of this Polygon." + return self[0] # First ring is the shell + exterior_ring = shell + + @property + def tuple(self): + "Returns a tuple of LinearRing coordinate tuples." + return tuple([self[i].tuple for i in xrange(self.geom_count)]) + coords = tuple + + @property + def point_count(self): + "The number of Points in this Polygon." + # Summing up the number of points in each ring of the Polygon. + return sum([self[i].point_count for i in xrange(self.geom_count)]) + + @property + def centroid(self): + "Returns the centroid (a Point) of this Polygon." + # The centroid is a Point, create a geometry for this. + p = OGRGeometry(OGRGeomType('Point')) + get_centroid(self._ptr, p._ptr) + return p + +# Geometry Collection base class. +class GeometryCollection(OGRGeometry): + "The Geometry Collection class." + + def __getitem__(self, index): + "Gets the Geometry at the specified index." + if index < 0 or index >= self.geom_count: + raise OGRIndexError('index out of range: %s' % index) + else: + return OGRGeometry(clone_geom(get_geom_ref(self._ptr, index)), self.srs) + + def __iter__(self): + "Iterates over each Geometry." + for i in xrange(self.geom_count): + yield self[i] + + def __len__(self): + "The number of geometries in this Geometry Collection." + return self.geom_count + + def add(self, geom): + "Add the geometry to this Geometry Collection." + if isinstance(geom, OGRGeometry): + if isinstance(geom, self.__class__): + for g in geom: add_geom(self._ptr, g._ptr) + else: + add_geom(self._ptr, geom._ptr) + elif isinstance(geom, basestring): + tmp = OGRGeometry(geom) + add_geom(self._ptr, tmp._ptr) + else: + raise OGRException('Must add an OGRGeometry.') + + @property + def point_count(self): + "The number of Points in this Geometry Collection." + # Summing up the number of points in each geometry in this collection + return sum([self[i].point_count for i in xrange(self.geom_count)]) + + @property + def tuple(self): + "Returns a tuple representation of this Geometry Collection." + return tuple([self[i].tuple for i in xrange(self.geom_count)]) + coords = tuple + +# Multiple Geometry types. +class MultiPoint(GeometryCollection): pass +class MultiLineString(GeometryCollection): pass +class MultiPolygon(GeometryCollection): pass + +# Class mapping dictionary (using the OGRwkbGeometryType as the key) +GEO_CLASSES = {1 : Point, + 2 : LineString, + 3 : Polygon, + 4 : MultiPoint, + 5 : MultiLineString, + 6 : MultiPolygon, + 7 : GeometryCollection, + 101: LinearRing, + } |