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 from types import UnicodeType

     9 

    10 # GEOS-related dependencies.

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

    12 from django.contrib.gis.geos.error import GEOSException

    13 from django.contrib.gis.geos.libgeos import GEOM_PTR

    14 

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

    16 # prototypes module -- which handles all interaction with

    17 # the underlying GEOS library.

    18 from django.contrib.gis.geos.prototypes import * 

    19 

    20 # Trying to import GDAL libraries, if available.  Have to place in

    21 # try/except since this package may be used outside GeoDjango.

    22 try:

    23     from django.contrib.gis.gdal import OGRGeometry, SpatialReference, GEOJSON

    24     HAS_GDAL = True

    25 except:

    26     HAS_GDAL, GEOJSON = False, False

    27 

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

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

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

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

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

    33 json_regex = re.compile(r'^\{.+\}$')

    34 

    35 class GEOSGeometry(object):

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

    37 

    38     # Initially, the geometry pointer is NULL

    39     _ptr = None

    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          * string: WKT

    48          * string: HEXEWKB (a PostGIS-specific canonical form)

    49          * buffer: WKB

    50         

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

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

    53         """ 

    54         if isinstance(geo_input, basestring):

    55             if isinstance(geo_input, UnicodeType):

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

    57                 geo_input = geo_input.encode('ascii')

    58                             

    59             wkt_m = wkt_regex.match(geo_input)

    60             if wkt_m:

    61                 # Handling WKT input.

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

    63                 g = from_wkt(wkt_m.group('wkt'))

    64             elif hex_regex.match(geo_input):

    65                 # Handling HEXEWKB input.

    66                 g = from_hex(geo_input, len(geo_input))

    67             elif GEOJSON and json_regex.match(geo_input):

    68                 # Handling GeoJSON input.

    69                 wkb_input = str(OGRGeometry(geo_input).wkb)

    70                 g = from_wkb(wkb_input, len(wkb_input))

    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             wkb_input = str(geo_input)

    79             g = from_wkb(wkb_input, len(wkb_input))

    80         else:

    81             # Invalid geometry type.

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

    83 

    84         if bool(g):

    85             # Setting the pointer object with a valid pointer.

    86             self._ptr = g

    87         else:

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

    89 

    90         # Post-initialization setup.

    91         self._post_init(srid)

    92 

    93     def _post_init(self, srid):

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

    95         # Setting the SRID, if given.

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

    97         

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

    99         self.__class__ = GEOS_CLASSES[self.geom_typeid]

   100 

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

   102         # geometries that do not have coordinate sequences)

   103         self._set_cs()

   104 

   105     @property

   106     def ptr(self):

   107         """

   108         Property for controlling access to the GEOS geometry pointer.  Using

   109         this raises an exception when the pointer is NULL, thus preventing

   110         the C library from attempting to access an invalid memory location.

   111         """

   112         if self._ptr: 

   113             return self._ptr

   114         else:

   115             raise GEOSException('NULL GEOS pointer encountered; was this geometry modified?')

   116 

   117     def __del__(self):

   118         """

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

   120         GEOS C++ object.

   121         """

   122         if self._ptr: destroy_geom(self._ptr)

   123 

   124     def __copy__(self):

   125         """

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

   127         invalid pointer location if the original is garbage collected.

   128         """

   129         return self.clone()

   130 

   131     def __deepcopy__(self, memodict):

   132         """

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

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

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

   136         """

   137         return self.clone()

   138 

   139     def __str__(self):

   140         "WKT is used for the string representation."

   141         return self.wkt

   142 

   143     def __repr__(self):

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

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

   146 

   147     # Pickling support

   148     def __getstate__(self):

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

   150         # and the SRID.

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

   152 

   153     def __setstate__(self, state):

   154         # Instantiating from the tuple state that was pickled.

   155         wkb, srid = state

   156         ptr = from_wkb(wkb, len(wkb))

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

   158         self._ptr = ptr

   159         self._post_init(srid)

   160 

   161     # Comparison operators

   162     def __eq__(self, other):

   163         """

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

   165         or a WKT representation.

   166         """

   167         if isinstance(other, basestring):

   168             return self.wkt == other

   169         elif isinstance(other, GEOSGeometry):

   170             return self.equals_exact(other)

   171         else:

   172             return False

   173 

   174     def __ne__(self, other):

   175         "The not equals operator."

   176         return not (self == other)

   177 

   178     ### Geometry set-like operations ###

   179     # Thanks to Sean Gillies for inspiration:

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

   181     # g = g1 | g2

   182     def __or__(self, other):

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

   184         return self.union(other)

   185 

   186     # g = g1 & g2

   187     def __and__(self, other):

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

   189         return self.intersection(other)

   190 

   191     # g = g1 - g2

   192     def __sub__(self, other):

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

   194         return self.difference(other)

   195 

   196     # g = g1 ^ g2

   197     def __xor__(self, other):

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

   199         return self.sym_difference(other)

   200 

   201     #### Coordinate Sequence Routines ####

   202     @property

   203     def has_cs(self):

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

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

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

   207             return True

   208         else:

   209             return False

   210 

   211     def _set_cs(self):

   212         "Sets the coordinate sequence for this Geometry."

   213         if self.has_cs:

   214             self._cs = GEOSCoordSeq(get_cs(self.ptr), self.hasz)

   215         else:

   216             self._cs = None

   217 

   218     @property

   219     def coord_seq(self):

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

   221         if self.has_cs:

   222             return self._cs.clone()

   223 

   224     #### Geometry Info ####

   225     @property

   226     def geom_type(self):

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

   228         return geos_type(self.ptr)

   229 

   230     @property

   231     def geom_typeid(self):

   232         "Returns an integer representing the Geometry type."

   233         return geos_typeid(self.ptr)

   234 

   235     @property

   236     def num_geom(self):

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

   238         return get_num_geoms(self.ptr)

   239 

   240     @property

   241     def num_coords(self):

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

   243         return get_num_coords(self.ptr)

   244 

   245     @property

   246     def num_points(self):

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

   248         return self.num_coords

   249 

   250     @property

   251     def dims(self):

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

   253         return get_dims(self.ptr)

   254 

   255     def normalize(self):

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

   257         return geos_normalize(self.ptr)

   258 

   259     #### Unary predicates ####

   260     @property

   261     def empty(self):

   262         """

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

   264         are empty.

   265         """

   266         return geos_isempty(self.ptr)

   267 

   268     @property

   269     def hasz(self):

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

   271         return geos_hasz(self.ptr)

   272 

   273     @property

   274     def ring(self):

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

   276         return geos_isring(self.ptr)

   277 

   278     @property

   279     def simple(self):

   280         "Returns false if the Geometry not simple."

   281         return geos_issimple(self.ptr)

   282 

   283     @property

   284     def valid(self):

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

   286         return geos_isvalid(self.ptr)

   287 

   288     #### Binary predicates. ####

   289     def contains(self, other):

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

   291         return 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 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 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 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 geos_equalsexact(self.ptr, other.ptr, float(tolerance))

   321 

   322     def intersects(self, other):

   323         "Returns true if disjoint returns false."

   324         return 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 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, str) or len(pattern) > 9:

   339             raise GEOSException('invalid intersection matrix pattern')

   340         return 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 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 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 = 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         geos_set_srid(self.ptr, srid)

   366     srid = property(get_srid, set_srid)

   367 

   368     #### Output Routines ####

   369     @property

   370     def ewkt(self):

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

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

   373         else: return self.wkt

   374 

   375     @property

   376     def wkt(self):

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

   378         return to_wkt(self.ptr)

   379 

   380     @property

   381     def hex(self):

   382         """

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

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

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

   386         """

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

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

   389         return to_hex(self.ptr, byref(c_size_t()))

   390 

   391     @property

   392     def json(self):

   393         """

   394         Returns GeoJSON representation of this Geometry if GDAL 1.5+ 

   395         is installed.

   396         """

   397         if GEOJSON: return self.ogr.json

   398     geojson = json

   399 

   400     @property

   401     def wkb(self):

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

   403         bin = to_wkb(self.ptr, byref(c_size_t()))

   404         return buffer(bin)

   405 

   406     @property

   407     def kml(self):

   408         "Returns the KML representation of this Geometry."

   409         gtype = self.geom_type

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

   411 

   412     #### GDAL-specific output routines ####

   413     @property

   414     def ogr(self):

   415         "Returns the OGR Geometry for this Geometry."

   416         if HAS_GDAL:

   417             if self.srid:

   418                 return OGRGeometry(self.wkb, self.srid)

   419             else:

   420                 return OGRGeometry(self.wkb)

   421         else:

   422             return None

   423 

   424     @property

   425     def srs(self):

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

   427         if HAS_GDAL and self.srid:

   428             return SpatialReference(self.srid)

   429         else:

   430             return None

   431 

   432     @property

   433     def crs(self):

   434         "Alias for `srs` property."

   435         return self.srs

   436 

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

   438         """

   439         Requires GDAL. Transforms the geometry according to the given 

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

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

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

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

   444         instead.

   445         """

   446         srid = self.srid

   447         if HAS_GDAL and srid:

   448             g = OGRGeometry(self.wkb, srid)

   449             g.transform(ct)

   450             wkb = str(g.wkb)

   451             ptr = from_wkb(wkb, len(wkb))

   452             if clone: 

   453                 # User wants a cloned transformed geometry returned.

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

   455             if ptr:

   456                 # Reassigning pointer, and performing post-initialization setup

   457                 # again due to the reassignment.

   458                 destroy_geom(self.ptr)

   459                 self._ptr = ptr

   460                 self._post_init(g.srid)

   461             else:

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

   463 

   464     #### Topology Routines ####

   465     def _topology(self, gptr):

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

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

   468 

   469     @property

   470     def boundary(self):

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

   472         return self._topology(geos_boundary(self.ptr))

   473 

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

   475         """

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

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

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

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

   480         (Text from PostGIS documentation at ch. 6.1.3)

   481         """

   482         return self._topology(geos_buffer(self.ptr, width, quadsegs))

   483 

   484     @property

   485     def centroid(self):

   486         """

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

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

   489         "weight" to the centroid).

   490         """

   491         return self._topology(geos_centroid(self.ptr))

   492 

   493     @property

   494     def convex_hull(self):

   495         """

   496         Returns the smallest convex Polygon that contains all the points 

   497         in the Geometry.

   498         """

   499         return self._topology(geos_convexhull(self.ptr))

   500 

   501     def difference(self, other):

   502         """

   503         Returns a Geometry representing the points making up this Geometry

   504         that do not make up other.

   505         """

   506         return self._topology(geos_difference(self.ptr, other.ptr))

   507 

   508     @property

   509     def envelope(self):

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

   511         return self._topology(geos_envelope(self.ptr))

   512 

   513     def intersection(self, other):

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

   515         return self._topology(geos_intersection(self.ptr, other.ptr))

   516 

   517     @property

   518     def point_on_surface(self):

   519         "Computes an interior point of this Geometry."

   520         return self._topology(geos_pointonsurface(self.ptr))

   521 

   522     def relate(self, other):

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

   524         return geos_relate(self.ptr, other.ptr)

   525 

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

   527         """

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

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

   530         tolerance provided, defaults to 0.

   531 

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

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

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

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

   536         input. This is significantly slower.         

   537         """

   538         if preserve_topology:

   539             return self._topology(geos_preservesimplify(self.ptr, tolerance))

   540         else:

   541             return self._topology(geos_simplify(self.ptr, tolerance))

   542 

   543     def sym_difference(self, other):

   544         """

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

   546         and the points in other not in this Geometry.

   547         """

   548         return self._topology(geos_symdifference(self.ptr, other.ptr))

   549 

   550     def union(self, other):

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

   552         return self._topology(geos_union(self.ptr, other.ptr))

   553 

   554     #### Other Routines ####

   555     @property

   556     def area(self):

   557         "Returns the area of the Geometry."

   558         return geos_area(self.ptr, byref(c_double()))

   559 

   560     def distance(self, other):

   561         """

   562         Returns the distance between the closest points on this Geometry

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

   564         the Geometry.

   565         """

   566         if not isinstance(other, GEOSGeometry): 

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

   568         return geos_distance(self.ptr, other.ptr, byref(c_double()))

   569 

   570     @property

   571     def extent(self):

   572         """

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

   574         (xmin, ymin, xmax, ymax).

   575         """

   576         env = self.envelope

   577         if isinstance(env, Point):

   578             xmin, ymin = env.tuple

   579             xmax, ymax = xmin, ymin

   580         else:

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

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

   583         return (xmin, ymin, xmax, ymax)

   584 

   585     @property

   586     def length(self):

   587         """

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

   589         circumfrence of a Polygon).

   590         """

   591         return geos_length(self.ptr, byref(c_double()))

   592     

   593     def clone(self):

   594         "Clones this Geometry."

   595         return GEOSGeometry(geom_clone(self.ptr), srid=self.srid)

   596 

   597 # Class mapping dictionary

   598 from django.contrib.gis.geos.geometries import Point, Polygon, LineString, LinearRing

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

   600 GEOS_CLASSES = {0 : Point,

   601                 1 : LineString,

   602                 2 : LinearRing,

   603                 3 : Polygon,

   604                 4 : MultiPoint,

   605                 5 : MultiLineString,

   606                 6 : MultiPolygon,

   607                 7 : GeometryCollection,

   608                 }