/django/contrib/gis/geos/geometry.py

  1"""
  2 This module contains the 'base' GEOSGeometry object -- all GEOS Geometries
  3 inherit from this object.
  4"""
  5# Python, ctypes and types dependencies.
  6import re
  7import warnings
  8from ctypes import addressof, byref, c_double, c_size_t
  9
 10# super-class for mutable list behavior
 11from django.contrib.gis.geos.mutable_list import ListMixin
 12
 13# GEOS-related dependencies.
 14from django.contrib.gis.geos.base import GEOSBase, gdal
 15from django.contrib.gis.geos.coordseq import GEOSCoordSeq
 16from django.contrib.gis.geos.error import GEOSException, GEOSIndexError
 17from django.contrib.gis.geos.libgeos import GEOM_PTR, GEOS_PREPARE
 18from django.contrib.gis.geos.mutable_list import ListMixin
 19
 20# All other functions in this module come from the ctypes
 21# prototypes module -- which handles all interaction with
 22# the underlying GEOS library.
 23from django.contrib.gis.geos import prototypes as capi
 24
 25# These functions provide access to a thread-local instance
 26# of their corresponding GEOS I/O class.
 27from django.contrib.gis.geos.prototypes.io import wkt_r, wkt_w, wkb_r, wkb_w, ewkb_w, ewkb_w3d
 28
 29# For recognizing geometry input.
 30from django.contrib.gis.geometry.regex import hex_regex, wkt_regex, json_regex
 31
 32class GEOSGeometry(GEOSBase, ListMixin):
 33    "A class that, generally, encapsulates a GEOS geometry."
 34
 35    # Raise GEOSIndexError instead of plain IndexError
 36    # (see ticket #4740 and GEOSIndexError docstring)
 37    _IndexError = GEOSIndexError
 38
 39    ptr_type = GEOM_PTR
 40
 41    #### Python 'magic' routines ####
 42    def __init__(self, geo_input, srid=None):
 43        """
 44        The base constructor for GEOS geometry objects, and may take the
 45        following inputs:
 46
 47         * strings:
 48            - WKT
 49            - HEXEWKB (a PostGIS-specific canonical form)
 50            - GeoJSON (requires GDAL)
 51         * buffer:
 52            - WKB
 53
 54        The `srid` keyword is used to specify the Source Reference Identifier
 55        (SRID) number for this Geometry.  If not set, the SRID will be None.
 56        """
 57        if isinstance(geo_input, basestring):
 58            if isinstance(geo_input, unicode):
 59                # Encoding to ASCII, WKT or HEXEWKB doesn't need any more.
 60                geo_input = geo_input.encode('ascii')
 61
 62            wkt_m = wkt_regex.match(geo_input)
 63            if wkt_m:
 64                # Handling WKT input.
 65                if wkt_m.group('srid'): srid = int(wkt_m.group('srid'))
 66                g = wkt_r().read(wkt_m.group('wkt'))
 67            elif hex_regex.match(geo_input):
 68                # Handling HEXEWKB input.
 69                g = wkb_r().read(geo_input)
 70            elif gdal.GEOJSON and json_regex.match(geo_input):
 71                # Handling GeoJSON input.
 72                g = wkb_r().read(gdal.OGRGeometry(geo_input).wkb)
 73            else:
 74                raise ValueError('String or unicode input unrecognized as WKT EWKT, and HEXEWKB.')
 75        elif isinstance(geo_input, GEOM_PTR):
 76            # When the input is a pointer to a geomtry (GEOM_PTR).
 77            g = geo_input
 78        elif isinstance(geo_input, buffer):
 79            # When the input is a buffer (WKB).
 80            g = wkb_r().read(geo_input)
 81        elif isinstance(geo_input, GEOSGeometry):
 82            g = capi.geom_clone(geo_input.ptr)
 83        else:
 84            # Invalid geometry type.
 85            raise TypeError('Improper geometry input type: %s' % str(type(geo_input)))
 86
 87        if bool(g):
 88            # Setting the pointer object with a valid pointer.
 89            self.ptr = g
 90        else:
 91            raise GEOSException('Could not initialize GEOS Geometry with given input.')
 92
 93        # Post-initialization setup.
 94        self._post_init(srid)
 95
 96    def _post_init(self, srid):
 97        "Helper routine for performing post-initialization setup."
 98        # Setting the SRID, if given.
 99        if srid and isinstance(srid, int): self.srid = srid
100
101        # Setting the class type (e.g., Point, Polygon, etc.)
102        self.__class__ = GEOS_CLASSES[self.geom_typeid]
103
104        # Setting the coordinate sequence for the geometry (will be None on
105        # geometries that do not have coordinate sequences)
106        self._set_cs()
107
108    def __del__(self):
109        """
110        Destroys this Geometry; in other words, frees the memory used by the
111        GEOS C++ object.
112        """
113        if self._ptr: capi.destroy_geom(self._ptr)
114
115    def __copy__(self):
116        """
117        Returns a clone because the copy of a GEOSGeometry may contain an
118        invalid pointer location if the original is garbage collected.
119        """
120        return self.clone()
121
122    def __deepcopy__(self, memodict):
123        """
124        The `deepcopy` routine is used by the `Node` class of django.utils.tree;
125        thus, the protocol routine needs to be implemented to return correct
126        copies (clones) of these GEOS objects, which use C pointers.
127        """
128        return self.clone()
129
130    def __str__(self):
131        "WKT is used for the string representation."
132        return self.wkt
133
134    def __repr__(self):
135        "Short-hand representation because WKT may be very large."
136        return '<%s object at %s>' % (self.geom_type, hex(addressof(self.ptr)))
137
138    # Pickling support
139    def __getstate__(self):
140        # The pickled state is simply a tuple of the WKB (in string form)
141        # and the SRID.
142        return str(self.wkb), self.srid
143
144    def __setstate__(self, state):
145        # Instantiating from the tuple state that was pickled.
146        wkb, srid = state
147        ptr = wkb_r().read(buffer(wkb))
148        if not ptr: raise GEOSException('Invalid Geometry loaded from pickled state.')
149        self.ptr = ptr
150        self._post_init(srid)
151
152    # Comparison operators
153    def __eq__(self, other):
154        """
155        Equivalence testing, a Geometry may be compared with another Geometry
156        or a WKT representation.
157        """
158        if isinstance(other, basestring):
159            return self.wkt == other
160        elif isinstance(other, GEOSGeometry):
161            return self.equals_exact(other)
162        else:
163            return False
164
165    def __ne__(self, other):
166        "The not equals operator."
167        return not (self == other)
168
169    ### Geometry set-like operations ###
170    # Thanks to Sean Gillies for inspiration:
171    #  http://lists.gispython.org/pipermail/community/2007-July/001034.html
172    # g = g1 | g2
173    def __or__(self, other):
174        "Returns the union of this Geometry and the other."
175        return self.union(other)
176
177    # g = g1 & g2
178    def __and__(self, other):
179        "Returns the intersection of this Geometry and the other."
180        return self.intersection(other)
181
182    # g = g1 - g2
183    def __sub__(self, other):
184        "Return the difference this Geometry and the other."
185        return self.difference(other)
186
187    # g = g1 ^ g2
188    def __xor__(self, other):
189        "Return the symmetric difference of this Geometry and the other."
190        return self.sym_difference(other)
191
192    #### Coordinate Sequence Routines ####
193    @property
194    def has_cs(self):
195        "Returns True if this Geometry has a coordinate sequence, False if not."
196        # Only these geometries are allowed to have coordinate sequences.
197        if isinstance(self, (Point, LineString, LinearRing)):
198            return True
199        else:
200            return False
201
202    def _set_cs(self):
203        "Sets the coordinate sequence for this Geometry."
204        if self.has_cs:
205            self._cs = GEOSCoordSeq(capi.get_cs(self.ptr), self.hasz)
206        else:
207            self._cs = None
208
209    @property
210    def coord_seq(self):
211        "Returns a clone of the coordinate sequence for this Geometry."
212        if self.has_cs:
213            return self._cs.clone()
214
215    #### Geometry Info ####
216    @property
217    def geom_type(self):
218        "Returns a string representing the Geometry type, e.g. 'Polygon'"
219        return capi.geos_type(self.ptr)
220
221    @property
222    def geom_typeid(self):
223        "Returns an integer representing the Geometry type."
224        return capi.geos_typeid(self.ptr)
225
226    @property
227    def num_geom(self):
228        "Returns the number of geometries in the Geometry."
229        return capi.get_num_geoms(self.ptr)
230
231    @property
232    def num_coords(self):
233        "Returns the number of coordinates in the Geometry."
234        return capi.get_num_coords(self.ptr)
235
236    @property
237    def num_points(self):
238        "Returns the number points, or coordinates, in the Geometry."
239        return self.num_coords
240
241    @property
242    def dims(self):
243        "Returns the dimension of this Geometry (0=point, 1=line, 2=surface)."
244        return capi.get_dims(self.ptr)
245
246    def normalize(self):
247        "Converts this Geometry to normal form (or canonical form)."
248        return capi.geos_normalize(self.ptr)
249
250    #### Unary predicates ####
251    @property
252    def empty(self):
253        """
254        Returns a boolean indicating whether the set of points in this Geometry
255        are empty.
256        """
257        return capi.geos_isempty(self.ptr)
258
259    @property
260    def hasz(self):
261        "Returns whether the geometry has a 3D dimension."
262        return capi.geos_hasz(self.ptr)
263
264    @property
265    def ring(self):
266        "Returns whether or not the geometry is a ring."
267        return capi.geos_isring(self.ptr)
268
269    @property
270    def simple(self):
271        "Returns false if the Geometry not simple."
272        return capi.geos_issimple(self.ptr)
273
274    @property
275    def valid(self):
276        "This property tests the validity of this Geometry."
277        return capi.geos_isvalid(self.ptr)
278
279    @property
280    def valid_reason(self):
281        """
282        Returns a string containing the reason for any invalidity.
283        """
284        if not GEOS_PREPARE:
285            raise GEOSException('Upgrade GEOS to 3.1 to get validity reason.')
286        return capi.geos_isvalidreason(self.ptr)
287
288    #### Binary predicates. ####
289    def contains(self, other):
290        "Returns true if other.within(this) returns true."
291        return capi.geos_contains(self.ptr, other.ptr)
292
293    def crosses(self, other):
294        """
295        Returns true if the DE-9IM intersection matrix for the two Geometries
296        is T*T****** (for a point and a curve,a point and an area or a line and
297        an area) 0******** (for two curves).
298        """
299        return capi.geos_crosses(self.ptr, other.ptr)
300
301    def disjoint(self, other):
302        """
303        Returns true if the DE-9IM intersection matrix for the two Geometries
304        is FF*FF****.
305        """
306        return capi.geos_disjoint(self.ptr, other.ptr)
307
308    def equals(self, other):
309        """
310        Returns true if the DE-9IM intersection matrix for the two Geometries
311        is T*F**FFF*.
312        """
313        return capi.geos_equals(self.ptr, other.ptr)
314
315    def equals_exact(self, other, tolerance=0):
316        """
317        Returns true if the two Geometries are exactly equal, up to a
318        specified tolerance.
319        """
320        return capi.geos_equalsexact(self.ptr, other.ptr, float(tolerance))
321
322    def intersects(self, other):
323        "Returns true if disjoint returns false."
324        return capi.geos_intersects(self.ptr, other.ptr)
325
326    def overlaps(self, other):
327        """
328        Returns true if the DE-9IM intersection matrix for the two Geometries
329        is T*T***T** (for two points or two surfaces) 1*T***T** (for two curves).
330        """
331        return capi.geos_overlaps(self.ptr, other.ptr)
332
333    def relate_pattern(self, other, pattern):
334        """
335        Returns true if the elements in the DE-9IM intersection matrix for the
336        two Geometries match the elements in pattern.
337        """
338        if not isinstance(pattern, basestring) or len(pattern) > 9:
339            raise GEOSException('invalid intersection matrix pattern')
340        return capi.geos_relatepattern(self.ptr, other.ptr, pattern)
341
342    def touches(self, other):
343        """
344        Returns true if the DE-9IM intersection matrix for the two Geometries
345        is FT*******, F**T***** or F***T****.
346        """
347        return capi.geos_touches(self.ptr, other.ptr)
348
349    def within(self, other):
350        """
351        Returns true if the DE-9IM intersection matrix for the two Geometries
352        is T*F**F***.
353        """
354        return capi.geos_within(self.ptr, other.ptr)
355
356    #### SRID Routines ####
357    def get_srid(self):
358        "Gets the SRID for the geometry, returns None if no SRID is set."
359        s = capi.geos_get_srid(self.ptr)
360        if s == 0: return None
361        else: return s
362
363    def set_srid(self, srid):
364        "Sets the SRID for the geometry."
365        capi.geos_set_srid(self.ptr, srid)
366    srid = property(get_srid, set_srid)
367
368    #### Output Routines ####
369    @property
370    def ewkt(self):
371        """
372        Returns the EWKT (WKT + SRID) of the Geometry.  Note that Z values
373        are *not* included in this representation because GEOS does not yet
374        support serializing them.
375        """
376        if self.get_srid(): return 'SRID=%s;%s' % (self.srid, self.wkt)
377        else: return self.wkt
378
379    @property
380    def wkt(self):
381        "Returns the WKT (Well-Known Text) representation of this Geometry."
382        return wkt_w().write(self)
383
384    @property
385    def hex(self):
386        """
387        Returns the WKB of this Geometry in hexadecimal form.  Please note
388        that the SRID and Z values are not included in this representation
389        because it is not a part of the OGC specification (use the `hexewkb`
390        property instead).
391        """
392        # A possible faster, all-python, implementation:
393        #  str(self.wkb).encode('hex')
394        return wkb_w().write_hex(self)
395
396    @property
397    def hexewkb(self):
398        """
399        Returns the EWKB of this Geometry in hexadecimal form.  This is an
400        extension of the WKB specification that includes SRID and Z values
401        that are a part of this geometry.
402        """
403        if self.hasz:
404            if not GEOS_PREPARE:
405                # See: http://trac.osgeo.org/geos/ticket/216
406                raise GEOSException('Upgrade GEOS to 3.1 to get valid 3D HEXEWKB.')
407            return ewkb_w3d().write_hex(self)
408        else:
409            return ewkb_w().write_hex(self)
410
411    @property
412    def json(self):
413        """
414        Returns GeoJSON representation of this Geometry if GDAL 1.5+
415        is installed.
416        """
417        if gdal.GEOJSON:
418            return self.ogr.json
419        else:
420            raise GEOSException('GeoJSON output only supported on GDAL 1.5+.')
421    geojson = json
422
423    @property
424    def wkb(self):
425        """
426        Returns the WKB (Well-Known Binary) representation of this Geometry
427        as a Python buffer.  SRID and Z values are not included, use the
428        `ewkb` property instead.
429        """
430        return wkb_w().write(self)
431
432    @property
433    def ewkb(self):
434        """
435        Return the EWKB representation of this Geometry as a Python buffer.
436        This is an extension of the WKB specification that includes any SRID
437        and Z values that are a part of this geometry.
438        """
439        if self.hasz:
440            if not GEOS_PREPARE:
441                # See: http://trac.osgeo.org/geos/ticket/216
442                raise GEOSException('Upgrade GEOS to 3.1 to get valid 3D EWKB.')
443            return ewkb_w3d().write(self)
444        else:
445            return ewkb_w().write(self)
446
447    @property
448    def kml(self):
449        "Returns the KML representation of this Geometry."
450        gtype = self.geom_type
451        return '<%s>%s</%s>' % (gtype, self.coord_seq.kml, gtype)
452
453    @property
454    def prepared(self):
455        """
456        Returns a PreparedGeometry corresponding to this geometry -- it is
457        optimized for the contains, intersects, and covers operations.
458        """
459        if GEOS_PREPARE:
460            return PreparedGeometry(self)
461        else:
462            raise GEOSException('GEOS 3.1+ required for prepared geometry support.')
463
464    #### GDAL-specific output routines ####
465    @property
466    def ogr(self):
467        "Returns the OGR Geometry for this Geometry."
468        if gdal.HAS_GDAL:
469            if self.srid:
470                return gdal.OGRGeometry(self.wkb, self.srid)
471            else:
472                return gdal.OGRGeometry(self.wkb)
473        else:
474            raise GEOSException('GDAL required to convert to an OGRGeometry.')
475
476    @property
477    def srs(self):
478        "Returns the OSR SpatialReference for SRID of this Geometry."
479        if gdal.HAS_GDAL:
480            if self.srid:
481                return gdal.SpatialReference(self.srid)
482            else:
483                return None
484        else:
485            raise GEOSException('GDAL required to return a SpatialReference object.')
486
487    @property
488    def crs(self):
489        "Alias for `srs` property."
490        return self.srs
491
492    def transform(self, ct, clone=False):
493        """
494        Requires GDAL. Transforms the geometry according to the given
495        transformation object, which may be an integer SRID, and WKT or
496        PROJ.4 string. By default, the geometry is transformed in-place and
497        nothing is returned. However if the `clone` keyword is set, then this
498        geometry will not be modified and a transformed clone will be returned
499        instead.
500        """
501        srid = self.srid
502
503        if ct == srid:
504            # short-circuit where source & dest SRIDs match
505            if clone:
506                return self.clone()
507            else:
508                return
509
510        if (srid is None) or (srid < 0):
511            warnings.warn("Calling transform() with no SRID set does no transformation!",
512                          stacklevel=2)
513            warnings.warn("Calling transform() with no SRID will raise GEOSException in v1.5",
514                          FutureWarning, stacklevel=2)
515            return
516
517        if not gdal.HAS_GDAL:
518            raise GEOSException("GDAL library is not available to transform() geometry.")
519
520        # Creating an OGR Geometry, which is then transformed.
521        g = gdal.OGRGeometry(self.wkb, srid)
522        g.transform(ct)
523        # Getting a new GEOS pointer
524        ptr = wkb_r().read(g.wkb)
525        if clone:
526            # User wants a cloned transformed geometry returned.
527            return GEOSGeometry(ptr, srid=g.srid)
528        if ptr:
529            # Reassigning pointer, and performing post-initialization setup
530            # again due to the reassignment.
531            capi.destroy_geom(self.ptr)
532            self.ptr = ptr
533            self._post_init(g.srid)
534        else:
535            raise GEOSException('Transformed WKB was invalid.')
536
537    #### Topology Routines ####
538    def _topology(self, gptr):
539        "Helper routine to return Geometry from the given pointer."
540        return GEOSGeometry(gptr, srid=self.srid)
541
542    @property
543    def boundary(self):
544        "Returns the boundary as a newly allocated Geometry object."
545        return self._topology(capi.geos_boundary(self.ptr))
546
547    def buffer(self, width, quadsegs=8):
548        """
549        Returns a geometry that represents all points whose distance from this
550        Geometry is less than or equal to distance. Calculations are in the
551        Spatial Reference System of this Geometry. The optional third parameter sets
552        the number of segment used to approximate a quarter circle (defaults to 8).
553        (Text from PostGIS documentation at ch. 6.1.3)
554        """
555        return self._topology(capi.geos_buffer(self.ptr, width, quadsegs))
556
557    @property
558    def centroid(self):
559        """
560        The centroid is equal to the centroid of the set of component Geometries
561        of highest dimension (since the lower-dimension geometries contribute zero
562        "weight" to the centroid).
563        """
564        return self._topology(capi.geos_centroid(self.ptr))
565
566    @property
567    def convex_hull(self):
568        """
569        Returns the smallest convex Polygon that contains all the points
570        in the Geometry.
571        """
572        return self._topology(capi.geos_convexhull(self.ptr))
573
574    def difference(self, other):
575        """
576        Returns a Geometry representing the points making up this Geometry
577        that do not make up other.
578        """
579        return self._topology(capi.geos_difference(self.ptr, other.ptr))
580
581    @property
582    def envelope(self):
583        "Return the envelope for this geometry (a polygon)."
584        return self._topology(capi.geos_envelope(self.ptr))
585
586    def intersection(self, other):
587        "Returns a Geometry representing the points shared by this Geometry and other."
588        return self._topology(capi.geos_intersection(self.ptr, other.ptr))
589
590    @property
591    def point_on_surface(self):
592        "Computes an interior point of this Geometry."
593        return self._topology(capi.geos_pointonsurface(self.ptr))
594
595    def relate(self, other):
596        "Returns the DE-9IM intersection matrix for this Geometry and the other."
597        return capi.geos_relate(self.ptr, other.ptr)
598
599    def simplify(self, tolerance=0.0, preserve_topology=False):
600        """
601        Returns the Geometry, simplified using the Douglas-Peucker algorithm
602        to the specified tolerance (higher tolerance => less points).  If no
603        tolerance provided, defaults to 0.
604
605        By default, this function does not preserve topology - e.g. polygons can
606        be split, collapse to lines or disappear holes can be created or
607        disappear, and lines can cross. By specifying preserve_topology=True,
608        the result will have the same dimension and number of components as the
609        input. This is significantly slower.
610        """
611        if preserve_topology:
612            return self._topology(capi.geos_preservesimplify(self.ptr, tolerance))
613        else:
614            return self._topology(capi.geos_simplify(self.ptr, tolerance))
615
616    def sym_difference(self, other):
617        """
618        Returns a set combining the points in this Geometry not in other,
619        and the points in other not in this Geometry.
620        """
621        return self._topology(capi.geos_symdifference(self.ptr, other.ptr))
622
623    def union(self, other):
624        "Returns a Geometry representing all the points in this Geometry and other."
625        return self._topology(capi.geos_union(self.ptr, other.ptr))
626
627    #### Other Routines ####
628    @property
629    def area(self):
630        "Returns the area of the Geometry."
631        return capi.geos_area(self.ptr, byref(c_double()))
632
633    def distance(self, other):
634        """
635        Returns the distance between the closest points on this Geometry
636        and the other. Units will be in those of the coordinate system of
637        the Geometry.
638        """
639        if not isinstance(other, GEOSGeometry):
640            raise TypeError('distance() works only on other GEOS Geometries.')
641        return capi.geos_distance(self.ptr, other.ptr, byref(c_double()))
642
643    @property
644    def extent(self):
645        """
646        Returns the extent of this geometry as a 4-tuple, consisting of
647        (xmin, ymin, xmax, ymax).
648        """
649        env = self.envelope
650        if isinstance(env, Point):
651            xmin, ymin = env.tuple
652            xmax, ymax = xmin, ymin
653        else:
654            xmin, ymin = env[0][0]
655            xmax, ymax = env[0][2]
656        return (xmin, ymin, xmax, ymax)
657
658    @property
659    def length(self):
660        """
661        Returns the length of this Geometry (e.g., 0 for point, or the
662        circumfrence of a Polygon).
663        """
664        return capi.geos_length(self.ptr, byref(c_double()))
665
666    def clone(self):
667        "Clones this Geometry."
668        return GEOSGeometry(capi.geom_clone(self.ptr), srid=self.srid)
669
670# Class mapping dictionary.  Has to be at the end to avoid import
671# conflicts with GEOSGeometry.
672from django.contrib.gis.geos.linestring import LineString, LinearRing
673from django.contrib.gis.geos.point import Point
674from django.contrib.gis.geos.polygon import Polygon
675from django.contrib.gis.geos.collections import GeometryCollection, MultiPoint, MultiLineString, MultiPolygon
676GEOS_CLASSES = {0 : Point,
677                1 : LineString,
678                2 : LinearRing,
679                3 : Polygon,
680                4 : MultiPoint,
681                5 : MultiLineString,
682                6 : MultiPolygon,
683                7 : GeometryCollection,
684                }
685
686# If supported, import the PreparedGeometry class.
687if GEOS_PREPARE:
688    from django.contrib.gis.geos.prepared import PreparedGeometry