""" This module contains the 'base' GEOSGeometry object -- all GEOS Geometries inherit from this object. """ # Python, ctypes and types dependencies. import re from ctypes import addressof, byref, c_double, c_size_t # super-class for mutable list behavior from django.contrib.gis.geos.mutable_list import ListMixin # GEOS-related dependencies. from django.contrib.gis.geos.base import GEOSBase, gdal from django.contrib.gis.geos.coordseq import GEOSCoordSeq from django.contrib.gis.geos.error import GEOSException, GEOSIndexError from django.contrib.gis.geos.libgeos import GEOM_PTR, GEOS_PREPARE from django.contrib.gis.geos.mutable_list import ListMixin # All other functions in this module come from the ctypes # prototypes module -- which handles all interaction with # the underlying GEOS library. from django.contrib.gis.geos import prototypes as capi # Regular expression for recognizing HEXEWKB and WKT. A prophylactic measure # to prevent potentially malicious input from reaching the underlying C # library. Not a substitute for good web security programming practices. hex_regex = re.compile(r'^[0-9A-F]+$', re.I) wkt_regex = re.compile(r'^(SRID=(?P\d+);)?(?P(POINT|LINESTRING|LINEARRING|POLYGON|MULTIPOINT|MULTILINESTRING|MULTIPOLYGON|GEOMETRYCOLLECTION)[ACEGIMLONPSRUTY\d,\.\-\(\) ]+)$', re.I) class GEOSGeometry(GEOSBase, ListMixin): "A class that, generally, encapsulates a GEOS geometry." # Raise GEOSIndexError instead of plain IndexError # (see ticket #4740 and GEOSIndexError docstring) _IndexError = GEOSIndexError ptr_type = GEOM_PTR #### Python 'magic' routines #### def __init__(self, geo_input, srid=None): """ The base constructor for GEOS geometry objects, and may take the following inputs: * strings: - WKT - HEXEWKB (a PostGIS-specific canonical form) - GeoJSON (requires GDAL) * buffer: - WKB The `srid` keyword is used to specify the Source Reference Identifier (SRID) number for this Geometry. If not set, the SRID will be None. """ if isinstance(geo_input, basestring): if isinstance(geo_input, unicode): # Encoding to ASCII, WKT or HEXEWKB doesn't need any more. geo_input = geo_input.encode('ascii') wkt_m = wkt_regex.match(geo_input) if wkt_m: # Handling WKT input. if wkt_m.group('srid'): srid = int(wkt_m.group('srid')) g = wkt_r.read(wkt_m.group('wkt')) elif hex_regex.match(geo_input): # Handling HEXEWKB input. g = wkb_r.read(geo_input) elif gdal.GEOJSON and gdal.geometries.json_regex.match(geo_input): # Handling GeoJSON input. g = wkb_r.read(gdal.OGRGeometry(geo_input).wkb) else: raise ValueError('String or unicode input unrecognized as WKT EWKT, and HEXEWKB.') elif isinstance(geo_input, GEOM_PTR): # When the input is a pointer to a geomtry (GEOM_PTR). g = geo_input elif isinstance(geo_input, buffer): # When the input is a buffer (WKB). g = wkb_r.read(geo_input) elif isinstance(geo_input, GEOSGeometry): g = capi.geom_clone(geo_input.ptr) else: # Invalid geometry type. raise TypeError('Improper geometry input type: %s' % str(type(geo_input))) if bool(g): # Setting the pointer object with a valid pointer. self.ptr = g else: raise GEOSException('Could not initialize GEOS Geometry with given input.') # Post-initialization setup. self._post_init(srid) super(GEOSGeometry, self).__init__() def _post_init(self, srid): "Helper routine for performing post-initialization setup." # Setting the SRID, if given. if srid and isinstance(srid, int): self.srid = srid # Setting the class type (e.g., Point, Polygon, etc.) self.__class__ = GEOS_CLASSES[self.geom_typeid] # Setting the coordinate sequence for the geometry (will be None on # geometries that do not have coordinate sequences) self._set_cs() def __del__(self): """ Destroys this Geometry; in other words, frees the memory used by the GEOS C++ object. """ if self._ptr: capi.destroy_geom(self._ptr) def __copy__(self): """ Returns a clone because the copy of a GEOSGeometry may contain an invalid pointer location if the original is garbage collected. """ return self.clone() def __deepcopy__(self, memodict): """ The `deepcopy` routine is used by the `Node` class of django.utils.tree; thus, the protocol routine needs to be implemented to return correct copies (clones) of these GEOS objects, which use C pointers. """ return self.clone() def __str__(self): "WKT is used for the string representation." return self.wkt def __repr__(self): "Short-hand representation because WKT may be very large." return '<%s object at %s>' % (self.geom_type, hex(addressof(self.ptr))) # Pickling support def __getstate__(self): # The pickled state is simply a tuple of the WKB (in string form) # and the SRID. return str(self.wkb), self.srid def __setstate__(self, state): # Instantiating from the tuple state that was pickled. wkb, srid = state ptr = capi.from_wkb(wkb, len(wkb)) if not ptr: raise GEOSException('Invalid Geometry loaded from pickled state.') self.ptr = ptr self._post_init(srid) # Comparison operators def __eq__(self, other): """ Equivalence testing, a Geometry may be compared with another Geometry or a WKT representation. """ if isinstance(other, basestring): return self.wkt == other elif isinstance(other, GEOSGeometry): return self.equals_exact(other) else: return False def __ne__(self, other): "The not equals operator." return not (self == other) ### Geometry set-like operations ### # Thanks to Sean Gillies for inspiration: # http://lists.gispython.org/pipermail/community/2007-July/001034.html # g = g1 | g2 def __or__(self, other): "Returns the union of this Geometry and the other." 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) #### Coordinate Sequence Routines #### @property def has_cs(self): "Returns True if this Geometry has a coordinate sequence, False if not." # Only these geometries are allowed to have coordinate sequences. if isinstance(self, (Point, LineString, LinearRing)): return True else: return False def _set_cs(self): "Sets the coordinate sequence for this Geometry." if self.has_cs: self._cs = GEOSCoordSeq(capi.get_cs(self.ptr), self.hasz) else: self._cs = None @property def coord_seq(self): "Returns a clone of the coordinate sequence for this Geometry." if self.has_cs: return self._cs.clone() #### Geometry Info #### @property def geom_type(self): "Returns a string representing the Geometry type, e.g. 'Polygon'" return capi.geos_type(self.ptr) @property def geom_typeid(self): "Returns an integer representing the Geometry type." return capi.geos_typeid(self.ptr) @property def num_geom(self): "Returns the number of geometries in the Geometry." return capi.get_num_geoms(self.ptr) @property def num_coords(self): "Returns the number of coordinates in the Geometry." return capi.get_num_coords(self.ptr) @property def num_points(self): "Returns the number points, or coordinates, in the Geometry." return self.num_coords @property def dims(self): "Returns the dimension of this Geometry (0=point, 1=line, 2=surface)." return capi.get_dims(self.ptr) def normalize(self): "Converts this Geometry to normal form (or canonical form)." return capi.geos_normalize(self.ptr) #### Unary predicates #### @property def empty(self): """ Returns a boolean indicating whether the set of points in this Geometry are empty. """ return capi.geos_isempty(self.ptr) @property def hasz(self): "Returns whether the geometry has a 3D dimension." return capi.geos_hasz(self.ptr) @property def ring(self): "Returns whether or not the geometry is a ring." return capi.geos_isring(self.ptr) @property def simple(self): "Returns false if the Geometry not simple." return capi.geos_issimple(self.ptr) @property def valid(self): "This property tests the validity of this Geometry." return capi.geos_isvalid(self.ptr) #### Binary predicates. #### def contains(self, other): "Returns true if other.within(this) returns true." return capi.geos_contains(self.ptr, other.ptr) def crosses(self, other): """ Returns true if the DE-9IM intersection matrix for the two Geometries is T*T****** (for a point and a curve,a point and an area or a line and an area) 0******** (for two curves). """ return capi.geos_crosses(self.ptr, other.ptr) def disjoint(self, other): """ Returns true if the DE-9IM intersection matrix for the two Geometries is FF*FF****. """ return capi.geos_disjoint(self.ptr, other.ptr) def equals(self, other): """ Returns true if the DE-9IM intersection matrix for the two Geometries is T*F**FFF*. """ return capi.geos_equals(self.ptr, other.ptr) def equals_exact(self, other, tolerance=0): """ Returns true if the two Geometries are exactly equal, up to a specified tolerance. """ return capi.geos_equalsexact(self.ptr, other.ptr, float(tolerance)) def intersects(self, other): "Returns true if disjoint returns false." return capi.geos_intersects(self.ptr, other.ptr) def overlaps(self, other): """ Returns true if the DE-9IM intersection matrix for the two Geometries is T*T***T** (for two points or two surfaces) 1*T***T** (for two curves). """ return capi.geos_overlaps(self.ptr, other.ptr) def relate_pattern(self, other, pattern): """ Returns true if the elements in the DE-9IM intersection matrix for the two Geometries match the elements in pattern. """ if not isinstance(pattern, basestring) or len(pattern) > 9: raise GEOSException('invalid intersection matrix pattern') return capi.geos_relatepattern(self.ptr, other.ptr, pattern) def touches(self, other): """ Returns true if the DE-9IM intersection matrix for the two Geometries is FT*******, F**T***** or F***T****. """ return capi.geos_touches(self.ptr, other.ptr) def within(self, other): """ Returns true if the DE-9IM intersection matrix for the two Geometries is T*F**F***. """ return capi.geos_within(self.ptr, other.ptr) #### SRID Routines #### def get_srid(self): "Gets the SRID for the geometry, returns None if no SRID is set." s = capi.geos_get_srid(self.ptr) if s == 0: return None else: return s def set_srid(self, srid): "Sets the SRID for the geometry." capi.geos_set_srid(self.ptr, srid) srid = property(get_srid, set_srid) #### Output Routines #### @property def ewkt(self): "Returns the EWKT (WKT + SRID) of the Geometry." if self.get_srid(): return 'SRID=%s;%s' % (self.srid, self.wkt) else: return self.wkt @property def wkt(self): "Returns the WKT (Well-Known Text) of the Geometry." return wkt_w.write(self) @property def hex(self): """ Returns the HEX of the Geometry -- please note that the SRID is not included in this representation, because the GEOS C library uses -1 by default, even if the SRID is set. """ # A possible faster, all-python, implementation: # str(self.wkb).encode('hex') return wkb_w.write_hex(self) @property def json(self): """ Returns GeoJSON representation of this Geometry if GDAL 1.5+ is installed. """ if gdal.GEOJSON: return self.ogr.json else: raise GEOSException('GeoJSON output only supported on GDAL 1.5+.') geojson = json @property def wkb(self): "Returns the WKB of the Geometry as a buffer." return wkb_w.write(self) @property def kml(self): "Returns the KML representation of this Geometry." gtype = self.geom_type return '<%s>%s' % (gtype, self.coord_seq.kml, gtype) @property def prepared(self): """ Returns a PreparedGeometry corresponding to this geometry -- it is optimized for the contains, intersects, and covers operations. """ if GEOS_PREPARE: return PreparedGeometry(self) else: raise GEOSException('GEOS 3.1+ required for prepared geometry support.') #### GDAL-specific output routines #### @property def ogr(self): "Returns the OGR Geometry for this Geometry." if gdal.HAS_GDAL: if self.srid: return gdal.OGRGeometry(self.wkb, self.srid) else: return gdal.OGRGeometry(self.wkb) else: raise GEOSException('GDAL required to convert to an OGRGeometry.') @property def srs(self): "Returns the OSR SpatialReference for SRID of this Geometry." if gdal.HAS_GDAL: if self.srid: return gdal.SpatialReference(self.srid) else: return None else: raise GEOSException('GDAL required to return a SpatialReference object.') @property def crs(self): "Alias for `srs` property." return self.srs def transform(self, ct, clone=False): """ Requires GDAL. Transforms the geometry according to the given transformation object, which may be an integer SRID, and WKT or PROJ.4 string. By default, the geometry is transformed in-place and nothing is returned. However if the `clone` keyword is set, then this geometry will not be modified and a transformed clone will be returned instead. """ srid = self.srid if gdal.HAS_GDAL and srid: # Creating an OGR Geometry, which is then transformed. g = gdal.OGRGeometry(self.wkb, srid) g.transform(ct) # Getting a new GEOS pointer ptr = wkb_r.read(g.wkb) if clone: # User wants a cloned transformed geometry returned. return GEOSGeometry(ptr, srid=g.srid) if ptr: # Reassigning pointer, and performing post-initialization setup # again due to the reassignment. capi.destroy_geom(self.ptr) self.ptr = ptr self._post_init(g.srid) else: raise GEOSException('Transformed WKB was invalid.') #### Topology Routines #### def _topology(self, gptr): "Helper routine to return Geometry from the given pointer." return GEOSGeometry(gptr, srid=self.srid) @property def boundary(self): "Returns the boundary as a newly allocated Geometry object." return self._topology(capi.geos_boundary(self.ptr)) def buffer(self, width, quadsegs=8): """ Returns a geometry that represents all points whose distance from this Geometry is less than or equal to distance. Calculations are in the Spatial Reference System of this Geometry. The optional third parameter sets the number of segment used to approximate a quarter circle (defaults to 8). (Text from PostGIS documentation at ch. 6.1.3) """ return self._topology(capi.geos_buffer(self.ptr, width, quadsegs)) @property def centroid(self): """ The centroid is equal to the centroid of the set of component Geometries of highest dimension (since the lower-dimension geometries contribute zero "weight" to the centroid). """ return self._topology(capi.geos_centroid(self.ptr)) @property def convex_hull(self): """ Returns the smallest convex Polygon that contains all the points in the Geometry. """ return self._topology(capi.geos_convexhull(self.ptr)) def difference(self, other): """ Returns a Geometry representing the points making up this Geometry that do not make up other. """ return self._topology(capi.geos_difference(self.ptr, other.ptr)) @property def envelope(self): "Return the envelope for this geometry (a polygon)." return self._topology(capi.geos_envelope(self.ptr)) def intersection(self, other): "Returns a Geometry representing the points shared by this Geometry and other." return self._topology(capi.geos_intersection(self.ptr, other.ptr)) @property def point_on_surface(self): "Computes an interior point of this Geometry." return self._topology(capi.geos_pointonsurface(self.ptr)) def relate(self, other): "Returns the DE-9IM intersection matrix for this Geometry and the other." return capi.geos_relate(self.ptr, other.ptr) def simplify(self, tolerance=0.0, preserve_topology=False): """ Returns the Geometry, simplified using the Douglas-Peucker algorithm to the specified tolerance (higher tolerance => less points). If no tolerance provided, defaults to 0. By default, this function does not preserve topology - e.g. polygons can be split, collapse to lines or disappear holes can be created or disappear, and lines can cross. By specifying preserve_topology=True, the result will have the same dimension and number of components as the input. This is significantly slower. """ if preserve_topology: return self._topology(capi.geos_preservesimplify(self.ptr, tolerance)) else: return self._topology(capi.geos_simplify(self.ptr, tolerance)) def sym_difference(self, other): """ Returns a set combining the points in this Geometry not in other, and the points in other not in this Geometry. """ return self._topology(capi.geos_symdifference(self.ptr, other.ptr)) def union(self, other): "Returns a Geometry representing all the points in this Geometry and other." return self._topology(capi.geos_union(self.ptr, other.ptr)) #### Other Routines #### @property def area(self): "Returns the area of the Geometry." return capi.geos_area(self.ptr, byref(c_double())) def distance(self, other): """ Returns the distance between the closest points on this Geometry and the other. Units will be in those of the coordinate system of the Geometry. """ if not isinstance(other, GEOSGeometry): raise TypeError('distance() works only on other GEOS Geometries.') return capi.geos_distance(self.ptr, other.ptr, byref(c_double())) @property def extent(self): """ Returns the extent of this geometry as a 4-tuple, consisting of (xmin, ymin, xmax, ymax). """ env = self.envelope if isinstance(env, Point): xmin, ymin = env.tuple xmax, ymax = xmin, ymin else: xmin, ymin = env[0][0] xmax, ymax = env[0][2] return (xmin, ymin, xmax, ymax) @property def length(self): """ Returns the length of this Geometry (e.g., 0 for point, or the circumfrence of a Polygon). """ return capi.geos_length(self.ptr, byref(c_double())) def clone(self): "Clones this Geometry." return GEOSGeometry(capi.geom_clone(self.ptr), srid=self.srid) # Class mapping dictionary. Has to be at the end to avoid import # conflicts with GEOSGeometry. from django.contrib.gis.geos.linestring import LineString, LinearRing from django.contrib.gis.geos.point import Point from django.contrib.gis.geos.polygon import Polygon from django.contrib.gis.geos.collections import GeometryCollection, MultiPoint, MultiLineString, MultiPolygon GEOS_CLASSES = {0 : Point, 1 : LineString, 2 : LinearRing, 3 : Polygon, 4 : MultiPoint, 5 : MultiLineString, 6 : MultiPolygon, 7 : GeometryCollection, } # Similarly, import the GEOS I/O instances here to avoid conflicts. from django.contrib.gis.geos.io import wkt_r, wkt_w, wkb_r, wkb_w # If supported, import the PreparedGeometry class. if GEOS_PREPARE: from django.contrib.gis.geos.prepared import PreparedGeometry