1 """

     2  This module contains the 'base' GEOSGeometry object -- all GEOS Geometries

     3  inherit from this object.

     4 """

     5 # Python, ctypes and types dependencies.

     6 import re

     7 from ctypes import addressof, byref, c_double, c_size_t

     8 

     9 # super-class for mutable list behavior

    10 from django.contrib.gis.geos.mutable_list import ListMixin

    11 

    12 # GEOS-related dependencies.

    13 from django.contrib.gis.geos.base import GEOSBase, gdal

    14 from django.contrib.gis.geos.coordseq import GEOSCoordSeq

    15 from django.contrib.gis.geos.error import GEOSException, GEOSIndexError

    16 from django.contrib.gis.geos.libgeos import GEOM_PTR, GEOS_PREPARE

    17 from django.contrib.gis.geos.mutable_list import ListMixin

    18 

    19 # All other functions in this module come from the ctypes

    20 # prototypes module -- which handles all interaction with

    21 # the underlying GEOS library.

    22 from django.contrib.gis.geos import prototypes as capi

    23 

    24 # Regular expression for recognizing HEXEWKB and WKT.  A prophylactic measure

    25 # to prevent potentially malicious input from reaching the underlying C

    26 # library.  Not a substitute for good web security programming practices.

    27 hex_regex = re.compile(r'^[0-9A-F]+$', re.I)

    28 wkt_regex = re.compile(r'^(SRID=(?P<srid>\d+);)?(?P<wkt>(POINT|LINESTRING|LINEARRING|POLYGON|MULTIPOINT|MULTILINESTRING|MULTIPOLYGON|GEOMETRYCOLLECTION)[ACEGIMLONPSRUTY\d,\.\-\(\) ]+)$', re.I)

    29 

    30 class GEOSGeometry(GEOSBase, ListMixin):

    31     "A class that, generally, encapsulates a GEOS geometry."

    32 

    33     # Raise GEOSIndexError instead of plain IndexError

    34     # (see ticket #4740 and GEOSIndexError docstring)

    35     _IndexError = GEOSIndexError

    36 

    37     ptr_type = GEOM_PTR

    38 

    39     #### Python 'magic' routines ####

    40     def __init__(self, geo_input, srid=None):

    41         """

    42         The base constructor for GEOS geometry objects, and may take the

    43         following inputs:

    44 

    45          * strings:

    46             - WKT

    47             - HEXEWKB (a PostGIS-specific canonical form)

    48             - GeoJSON (requires GDAL)

    49          * buffer:

    50             - WKB

    51 

    52         The `srid` keyword is used to specify the Source Reference Identifier

    53         (SRID) number for this Geometry.  If not set, the SRID will be None.

    54         """

    55         if isinstance(geo_input, basestring):

    56             if isinstance(geo_input, unicode):

    57                 # Encoding to ASCII, WKT or HEXEWKB doesn't need any more.

    58                 geo_input = geo_input.encode('ascii')

    59 

    60             wkt_m = wkt_regex.match(geo_input)

    61             if wkt_m:

    62                 # Handling WKT input.

    63                 if wkt_m.group('srid'): srid = int(wkt_m.group('srid'))

    64                 g = wkt_r.read(wkt_m.group('wkt'))

    65             elif hex_regex.match(geo_input):

    66                 # Handling HEXEWKB input.

    67                 g = wkb_r.read(geo_input)

    68             elif gdal.GEOJSON and gdal.geometries.json_regex.match(geo_input):

    69                 # Handling GeoJSON input.

    70                 g = wkb_r.read(gdal.OGRGeometry(geo_input).wkb)

    71             else:

    72                 raise ValueError('String or unicode input unrecognized as WKT EWKT, and HEXEWKB.')

    73         elif isinstance(geo_input, GEOM_PTR):

    74             # When the input is a pointer to a geomtry (GEOM_PTR).

    75             g = geo_input

    76         elif isinstance(geo_input, buffer):

    77             # When the input is a buffer (WKB).

    78             g = wkb_r.read(geo_input)

    79         elif isinstance(geo_input, GEOSGeometry):

    80             g = capi.geom_clone(geo_input.ptr)

    81         else:

    82             # Invalid geometry type.

    83             raise TypeError('Improper geometry input type: %s' % str(type(geo_input)))

    84 

    85         if bool(g):

    86             # Setting the pointer object with a valid pointer.

    87             self.ptr = g

    88         else:

    89             raise GEOSException('Could not initialize GEOS Geometry with given input.')

    90 

    91         # Post-initialization setup.

    92         self._post_init(srid)

    93 

    94     def _post_init(self, srid):

    95         "Helper routine for performing post-initialization setup."

    96         # Setting the SRID, if given.

    97         if srid and isinstance(srid, int): self.srid = srid

    98 

    99         # Setting the class type (e.g., Point, Polygon, etc.)

   100         self.__class__ = GEOS_CLASSES[self.geom_typeid]

   101 

   102         # Setting the coordinate sequence for the geometry (will be None on

   103         # geometries that do not have coordinate sequences)

   104         self._set_cs()

   105 

   106     def __del__(self):

   107         """

   108         Destroys this Geometry; in other words, frees the memory used by the

   109         GEOS C++ object.

   110         """

   111         if self._ptr: capi.destroy_geom(self._ptr)

   112 

   113     def __copy__(self):

   114         """

   115         Returns a clone because the copy of a GEOSGeometry may contain an

   116         invalid pointer location if the original is garbage collected.

   117         """

   118         return self.clone()

   119 

   120     def __deepcopy__(self, memodict):

   121         """

   122         The `deepcopy` routine is used by the `Node` class of django.utils.tree;

   123         thus, the protocol routine needs to be implemented to return correct

   124         copies (clones) of these GEOS objects, which use C pointers.

   125         """

   126         return self.clone()

   127 

   128     def __str__(self):

   129         "WKT is used for the string representation."

   130         return self.wkt

   131 

   132     def __repr__(self):

   133         "Short-hand representation because WKT may be very large."

   134         return '<%s object at %s>' % (self.geom_type, hex(addressof(self.ptr)))

   135 

   136     # Pickling support

   137     def __getstate__(self):

   138         # The pickled state is simply a tuple of the WKB (in string form)

   139         # and the SRID.

   140         return str(self.wkb), self.srid

   141 

   142     def __setstate__(self, state):

   143         # Instantiating from the tuple state that was pickled.

   144         wkb, srid = state

   145         ptr = capi.from_wkb(wkb, len(wkb))

   146         if not ptr: raise GEOSException('Invalid Geometry loaded from pickled state.')

   147         self.ptr = ptr

   148         self._post_init(srid)

   149 

   150     # Comparison operators

   151     def __eq__(self, other):

   152         """

   153         Equivalence testing, a Geometry may be compared with another Geometry

   154         or a WKT representation.

   155         """

   156         if isinstance(other, basestring):

   157             return self.wkt == other

   158         elif isinstance(other, GEOSGeometry):

   159             return self.equals_exact(other)

   160         else:

   161             return False

   162 

   163     def __ne__(self, other):

   164         "The not equals operator."

   165         return not (self == other)

   166 

   167     ### Geometry set-like operations ###

   168     # Thanks to Sean Gillies for inspiration:

   169     #  http://lists.gispython.org/pipermail/community/2007-July/001034.html

   170     # g = g1 | g2

   171     def __or__(self, other):

   172         "Returns the union of this Geometry and the other."

   173         return self.union(other)

   174 

   175     # g = g1 & g2

   176     def __and__(self, other):

   177         "Returns the intersection of this Geometry and the other."

   178         return self.intersection(other)

   179 

   180     # g = g1 - g2

   181     def __sub__(self, other):

   182         "Return the difference this Geometry and the other."

   183         return self.difference(other)

   184 

   185     # g = g1 ^ g2

   186     def __xor__(self, other):

   187         "Return the symmetric difference of this Geometry and the other."

   188         return self.sym_difference(other)

   189 

   190     #### Coordinate Sequence Routines ####

   191     @property

   192     def has_cs(self):

   193         "Returns True if this Geometry has a coordinate sequence, False if not."

   194         # Only these geometries are allowed to have coordinate sequences.

   195         if isinstance(self, (Point, LineString, LinearRing)):

   196             return True

   197         else:

   198             return False

   199 

   200     def _set_cs(self):

   201         "Sets the coordinate sequence for this Geometry."

   202         if self.has_cs:

   203             self._cs = GEOSCoordSeq(capi.get_cs(self.ptr), self.hasz)

   204         else:

   205             self._cs = None

   206 

   207     @property

   208     def coord_seq(self):

   209         "Returns a clone of the coordinate sequence for this Geometry."

   210         if self.has_cs:

   211             return self._cs.clone()

   212 

   213     #### Geometry Info ####

   214     @property

   215     def geom_type(self):

   216         "Returns a string representing the Geometry type, e.g. 'Polygon'"

   217         return capi.geos_type(self.ptr)

   218 

   219     @property

   220     def geom_typeid(self):

   221         "Returns an integer representing the Geometry type."

   222         return capi.geos_typeid(self.ptr)

   223 

   224     @property

   225     def num_geom(self):

   226         "Returns the number of geometries in the Geometry."

   227         return capi.get_num_geoms(self.ptr)

   228 

   229     @property

   230     def num_coords(self):

   231         "Returns the number of coordinates in the Geometry."

   232         return capi.get_num_coords(self.ptr)

   233 

   234     @property

   235     def num_points(self):

   236         "Returns the number points, or coordinates, in the Geometry."

   237         return self.num_coords

   238 

   239     @property

   240     def dims(self):

   241         "Returns the dimension of this Geometry (0=point, 1=line, 2=surface)."

   242         return capi.get_dims(self.ptr)

   243 

   244     def normalize(self):

   245         "Converts this Geometry to normal form (or canonical form)."

   246         return capi.geos_normalize(self.ptr)

   247 

   248     #### Unary predicates ####

   249     @property

   250     def empty(self):

   251         """

   252         Returns a boolean indicating whether the set of points in this Geometry

   253         are empty.

   254         """

   255         return capi.geos_isempty(self.ptr)

   256 

   257     @property

   258     def hasz(self):

   259         "Returns whether the geometry has a 3D dimension."

   260         return capi.geos_hasz(self.ptr)

   261 

   262     @property

   263     def ring(self):

   264         "Returns whether or not the geometry is a ring."

   265         return capi.geos_isring(self.ptr)

   266 

   267     @property

   268     def simple(self):

   269         "Returns false if the Geometry not simple."

   270         return capi.geos_issimple(self.ptr)

   271 

   272     @property

   273     def valid(self):

   274         "This property tests the validity of this Geometry."

   275         return capi.geos_isvalid(self.ptr)

   276 

   277     #### Binary predicates. ####

   278     def contains(self, other):

   279         "Returns true if other.within(this) returns true."

   280         return capi.geos_contains(self.ptr, other.ptr)

   281 

   282     def crosses(self, other):

   283         """

   284         Returns true if the DE-9IM intersection matrix for the two Geometries

   285         is T*T****** (for a point and a curve,a point and an area or a line and

   286         an area) 0******** (for two curves).

   287         """

   288         return capi.geos_crosses(self.ptr, other.ptr)

   289 

   290     def disjoint(self, other):

   291         """

   292         Returns true if the DE-9IM intersection matrix for the two Geometries

   293         is FF*FF****.

   294         """

   295         return capi.geos_disjoint(self.ptr, other.ptr)

   296 

   297     def equals(self, other):

   298         """

   299         Returns true if the DE-9IM intersection matrix for the two Geometries

   300         is T*F**FFF*.

   301         """

   302         return capi.geos_equals(self.ptr, other.ptr)

   303 

   304     def equals_exact(self, other, tolerance=0):

   305         """

   306         Returns true if the two Geometries are exactly equal, up to a

   307         specified tolerance.

   308         """

   309         return capi.geos_equalsexact(self.ptr, other.ptr, float(tolerance))

   310 

   311     def intersects(self, other):

   312         "Returns true if disjoint returns false."

   313         return capi.geos_intersects(self.ptr, other.ptr)

   314 

   315     def overlaps(self, other):

   316         """

   317         Returns true if the DE-9IM intersection matrix for the two Geometries

   318         is T*T***T** (for two points or two surfaces) 1*T***T** (for two curves).

   319         """

   320         return capi.geos_overlaps(self.ptr, other.ptr)

   321 

   322     def relate_pattern(self, other, pattern):

   323         """

   324         Returns true if the elements in the DE-9IM intersection matrix for the

   325         two Geometries match the elements in pattern.

   326         """

   327         if not isinstance(pattern, basestring) or len(pattern) > 9:

   328             raise GEOSException('invalid intersection matrix pattern')

   329         return capi.geos_relatepattern(self.ptr, other.ptr, pattern)

   330 

   331     def touches(self, other):

   332         """

   333         Returns true if the DE-9IM intersection matrix for the two Geometries

   334         is FT*******, F**T***** or F***T****.

   335         """

   336         return capi.geos_touches(self.ptr, other.ptr)

   337 

   338     def within(self, other):

   339         """

   340         Returns true if the DE-9IM intersection matrix for the two Geometries

   341         is T*F**F***.

   342         """

   343         return capi.geos_within(self.ptr, other.ptr)

   344 

   345     #### SRID Routines ####

   346     def get_srid(self):

   347         "Gets the SRID for the geometry, returns None if no SRID is set."

   348         s = capi.geos_get_srid(self.ptr)

   349         if s == 0: return None

   350         else: return s

   351 

   352     def set_srid(self, srid):

   353         "Sets the SRID for the geometry."

   354         capi.geos_set_srid(self.ptr, srid)

   355     srid = property(get_srid, set_srid)

   356 

   357     #### Output Routines ####

   358     @property

   359     def ewkt(self):

   360         "Returns the EWKT (WKT + SRID) of the Geometry."

   361         if self.get_srid(): return 'SRID=%s;%s' % (self.srid, self.wkt)

   362         else: return self.wkt

   363 

   364     @property

   365     def wkt(self):

   366         "Returns the WKT (Well-Known Text) of the Geometry."

   367         return wkt_w.write(self)

   368 

   369     @property

   370     def hex(self):

   371         """

   372         Returns the HEX of the Geometry -- please note that the SRID is not

   373         included in this representation, because the GEOS C library uses

   374         -1 by default, even if the SRID is set.

   375         """

   376         # A possible faster, all-python, implementation:

   377         #  str(self.wkb).encode('hex')

   378         return wkb_w.write_hex(self)

   379 

   380     @property

   381     def json(self):

   382         """

   383         Returns GeoJSON representation of this Geometry if GDAL 1.5+

   384         is installed.

   385         """

   386         if gdal.GEOJSON: 

   387             return self.ogr.json

   388         else:

   389             raise GEOSException('GeoJSON output only supported on GDAL 1.5+.')

   390     geojson = json

   391 

   392     @property

   393     def wkb(self):

   394         "Returns the WKB of the Geometry as a buffer."

   395         return wkb_w.write(self)

   396 

   397     @property

   398     def kml(self):

   399         "Returns the KML representation of this Geometry."

   400         gtype = self.geom_type

   401         return '<%s>%s</%s>' % (gtype, self.coord_seq.kml, gtype)

   402 

   403     @property

   404     def prepared(self):

   405         """

   406         Returns a PreparedGeometry corresponding to this geometry -- it is

   407         optimized for the contains, intersects, and covers operations.

   408         """

   409         if GEOS_PREPARE:

   410             return PreparedGeometry(self)

   411         else:

   412             raise GEOSException('GEOS 3.1+ required for prepared geometry support.')

   413 

   414     #### GDAL-specific output routines ####

   415     @property

   416     def ogr(self):

   417         "Returns the OGR Geometry for this Geometry."

   418         if gdal.HAS_GDAL:

   419             if self.srid:

   420                 return gdal.OGRGeometry(self.wkb, self.srid)

   421             else:

   422                 return gdal.OGRGeometry(self.wkb)

   423         else:

   424             raise GEOSException('GDAL required to convert to an OGRGeometry.')

   425 

   426     @property

   427     def srs(self):

   428         "Returns the OSR SpatialReference for SRID of this Geometry."

   429         if gdal.HAS_GDAL:

   430             if self.srid:

   431                 return gdal.SpatialReference(self.srid)

   432             else:

   433                 return None

   434         else:

   435             raise GEOSException('GDAL required to return a SpatialReference object.')

   436 

   437     @property

   438     def crs(self):

   439         "Alias for `srs` property."

   440         return self.srs

   441 

   442     def transform(self, ct, clone=False):

   443         """

   444         Requires GDAL. Transforms the geometry according to the given

   445         transformation object, which may be an integer SRID, and WKT or

   446         PROJ.4 string. By default, the geometry is transformed in-place and

   447         nothing is returned. However if the `clone` keyword is set, then this

   448         geometry will not be modified and a transformed clone will be returned

   449         instead.

   450         """

   451         srid = self.srid

   452         if gdal.HAS_GDAL and srid:

   453             # Creating an OGR Geometry, which is then transformed.

   454             g = gdal.OGRGeometry(self.wkb, srid)

   455             g.transform(ct)

   456             # Getting a new GEOS pointer

   457             ptr = wkb_r.read(g.wkb)

   458             if clone:

   459                 # User wants a cloned transformed geometry returned.

   460                 return GEOSGeometry(ptr, srid=g.srid)

   461             if ptr:

   462                 # Reassigning pointer, and performing post-initialization setup

   463                 # again due to the reassignment.

   464                 capi.destroy_geom(self.ptr)

   465                 self.ptr = ptr

   466                 self._post_init(g.srid)

   467             else:

   468                 raise GEOSException('Transformed WKB was invalid.')

   469 

   470     #### Topology Routines ####

   471     def _topology(self, gptr):

   472         "Helper routine to return Geometry from the given pointer."

   473         return GEOSGeometry(gptr, srid=self.srid)

   474 

   475     @property

   476     def boundary(self):

   477         "Returns the boundary as a newly allocated Geometry object."

   478         return self._topology(capi.geos_boundary(self.ptr))

   479 

   480     def buffer(self, width, quadsegs=8):

   481         """

   482         Returns a geometry that represents all points whose distance from this

   483         Geometry is less than or equal to distance. Calculations are in the

   484         Spatial Reference System of this Geometry. The optional third parameter sets

   485         the number of segment used to approximate a quarter circle (defaults to 8).

   486         (Text from PostGIS documentation at ch. 6.1.3)

   487         """

   488         return self._topology(capi.geos_buffer(self.ptr, width, quadsegs))

   489 

   490     @property

   491     def centroid(self):

   492         """

   493         The centroid is equal to the centroid of the set of component Geometries

   494         of highest dimension (since the lower-dimension geometries contribute zero

   495         "weight" to the centroid).

   496         """

   497         return self._topology(capi.geos_centroid(self.ptr))

   498 

   499     @property

   500     def convex_hull(self):

   501         """

   502         Returns the smallest convex Polygon that contains all the points

   503         in the Geometry.

   504         """

   505         return self._topology(capi.geos_convexhull(self.ptr))

   506 

   507     def difference(self, other):

   508         """

   509         Returns a Geometry representing the points making up this Geometry

   510         that do not make up other.

   511         """

   512         return self._topology(capi.geos_difference(self.ptr, other.ptr))

   513 

   514     @property

   515     def envelope(self):

   516         "Return the envelope for this geometry (a polygon)."

   517         return self._topology(capi.geos_envelope(self.ptr))

   518 

   519     def intersection(self, other):

   520         "Returns a Geometry representing the points shared by this Geometry and other."

   521         return self._topology(capi.geos_intersection(self.ptr, other.ptr))

   522 

   523     @property

   524     def point_on_surface(self):

   525         "Computes an interior point of this Geometry."

   526         return self._topology(capi.geos_pointonsurface(self.ptr))

   527 

   528     def relate(self, other):

   529         "Returns the DE-9IM intersection matrix for this Geometry and the other."

   530         return capi.geos_relate(self.ptr, other.ptr)

   531 

   532     def simplify(self, tolerance=0.0, preserve_topology=False):

   533         """

   534         Returns the Geometry, simplified using the Douglas-Peucker algorithm

   535         to the specified tolerance (higher tolerance => less points).  If no

   536         tolerance provided, defaults to 0.

   537 

   538         By default, this function does not preserve topology - e.g. polygons can

   539         be split, collapse to lines or disappear holes can be created or

   540         disappear, and lines can cross. By specifying preserve_topology=True,

   541         the result will have the same dimension and number of components as the

   542         input. This is significantly slower.

   543         """

   544         if preserve_topology:

   545             return self._topology(capi.geos_preservesimplify(self.ptr, tolerance))

   546         else:

   547             return self._topology(capi.geos_simplify(self.ptr, tolerance))

   548 

   549     def sym_difference(self, other):

   550         """

   551         Returns a set combining the points in this Geometry not in other,

   552         and the points in other not in this Geometry.

   553         """

   554         return self._topology(capi.geos_symdifference(self.ptr, other.ptr))

   555 

   556     def union(self, other):

   557         "Returns a Geometry representing all the points in this Geometry and other."

   558         return self._topology(capi.geos_union(self.ptr, other.ptr))

   559 

   560     #### Other Routines ####

   561     @property

   562     def area(self):

   563         "Returns the area of the Geometry."

   564         return capi.geos_area(self.ptr, byref(c_double()))

   565 

   566     def distance(self, other):

   567         """

   568         Returns the distance between the closest points on this Geometry

   569         and the other. Units will be in those of the coordinate system of

   570         the Geometry.

   571         """

   572         if not isinstance(other, GEOSGeometry):

   573             raise TypeError('distance() works only on other GEOS Geometries.')

   574         return capi.geos_distance(self.ptr, other.ptr, byref(c_double()))

   575 

   576     @property

   577     def extent(self):

   578         """

   579         Returns the extent of this geometry as a 4-tuple, consisting of

   580         (xmin, ymin, xmax, ymax).

   581         """

   582         env = self.envelope

   583         if isinstance(env, Point):

   584             xmin, ymin = env.tuple

   585             xmax, ymax = xmin, ymin

   586         else:

   587             xmin, ymin = env[0][0]

   588             xmax, ymax = env[0][2]

   589         return (xmin, ymin, xmax, ymax)

   590 

   591     @property

   592     def length(self):

   593         """

   594         Returns the length of this Geometry (e.g., 0 for point, or the

   595         circumfrence of a Polygon).

   596         """

   597         return capi.geos_length(self.ptr, byref(c_double()))

   598 

   599     def clone(self):

   600         "Clones this Geometry."

   601         return GEOSGeometry(capi.geom_clone(self.ptr), srid=self.srid)

   602 

   603 # Class mapping dictionary.  Has to be at the end to avoid import

   604 # conflicts with GEOSGeometry.

   605 from django.contrib.gis.geos.linestring import LineString, LinearRing

   606 from django.contrib.gis.geos.point import Point

   607 from django.contrib.gis.geos.polygon import Polygon

   608 from django.contrib.gis.geos.collections import GeometryCollection, MultiPoint, MultiLineString, MultiPolygon

   609 GEOS_CLASSES = {0 : Point,

   610                 1 : LineString,

   611                 2 : LinearRing,

   612                 3 : Polygon,

   613                 4 : MultiPoint,

   614                 5 : MultiLineString,

   615                 6 : MultiPolygon,

   616                 7 : GeometryCollection,

   617                 }

   618 

   619 # Similarly, import the GEOS I/O instances here to avoid conflicts.

   620 from django.contrib.gis.geos.io import wkt_r, wkt_w, wkb_r, wkb_w

   621 

   622 # If supported, import the PreparedGeometry class.

   623 if GEOS_PREPARE:

   624     from django.contrib.gis.geos.prepared import PreparedGeometry