1 """

     2  The OGRGeometry is a wrapper for using the OGR Geometry class

     3  (see http://www.gdal.org/ogr/classOGRGeometry.html).  OGRGeometry

     4  may be instantiated when reading geometries from OGR Data Sources

     5  (e.g. SHP files), or when given OGC WKT (a string).

     6 

     7  While the 'full' API is not present yet, the API is "pythonic" unlike

     8  the traditional and "next-generation" OGR Python bindings.  One major

     9  advantage OGR Geometries have over their GEOS counterparts is support

    10  for spatial reference systems and their transformation.

    11 

    12  Example:

    13   >>> from django.contrib.gis.gdal import OGRGeometry, OGRGeomType, SpatialReference

    14   >>> wkt1, wkt2 = 'POINT(-90 30)', 'POLYGON((0 0, 5 0, 5 5, 0 5)'

    15   >>> pnt = OGRGeometry(wkt1)

    16   >>> print pnt

    17   POINT (-90 30)

    18   >>> mpnt = OGRGeometry(OGRGeomType('MultiPoint'), SpatialReference('WGS84'))

    19   >>> mpnt.add(wkt1)

    20   >>> mpnt.add(wkt1)

    21   >>> print mpnt

    22   MULTIPOINT (-90 30,-90 30)

    23   >>> print mpnt.srs.name

    24   WGS 84

    25   >>> print mpnt.srs.proj

    26   +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

    27   >>> mpnt.transform_to(SpatialReference('NAD27'))

    28   >>> print mpnt.proj

    29   +proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs

    30   >>> print mpnt

    31   MULTIPOINT (-89.999930378602485 29.999797886557641,-89.999930378602485 29.999797886557641)

    32   

    33   The OGRGeomType class is to make it easy to specify an OGR geometry type:

    34   >>> from django.contrib.gis.gdal import OGRGeomType

    35   >>> gt1 = OGRGeomType(3) # Using an integer for the type

    36   >>> gt2 = OGRGeomType('Polygon') # Using a string

    37   >>> gt3 = OGRGeomType('POLYGON') # It's case-insensitive

    38   >>> print gt1 == 3, gt1 == 'Polygon' # Equivalence works w/non-OGRGeomType objects

    39   True

    40 """

    41 # Python library requisites.

    42 import re, sys

    43 from binascii import a2b_hex

    44 from ctypes import byref, string_at, c_char_p, c_double, c_ubyte, c_void_p

    45 from types import UnicodeType

    46 

    47 # Getting GDAL prerequisites

    48 from django.contrib.gis.gdal.envelope import Envelope, OGREnvelope

    49 from django.contrib.gis.gdal.error import OGRException, OGRIndexError, SRSException

    50 from django.contrib.gis.gdal.geomtype import OGRGeomType

    51 from django.contrib.gis.gdal.srs import SpatialReference, CoordTransform

    52 

    53 # Getting the ctypes prototype functions that interface w/the GDAL C library.

    54 from django.contrib.gis.gdal.prototypes.geom import *

    55 from django.contrib.gis.gdal.prototypes.srs import clone_srs

    56 

    57 # For more information, see the OGR C API source code:

    58 #  http://www.gdal.org/ogr/ogr__api_8h.html

    59 #

    60 # The OGR_G_* routines are relevant here.

    61 

    62 # Regular expressions for recognizing HEXEWKB and WKT.

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

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

    65 json_regex = re.compile(r'^\{[\s\w,\-\.\"\'\:\[\]]+\}$')

    66 

    67 #### OGRGeometry Class ####

    68 class OGRGeometry(object):

    69     "Generally encapsulates an OGR geometry."

    70 

    71     def __init__(self, geom_input, srs=None):

    72         "Initializes Geometry on either WKT or an OGR pointer as input."

    73 

    74         self._ptr = c_void_p(None) # Initially NULL

    75         str_instance = isinstance(geom_input, basestring)

    76 

    77         # If HEX, unpack input to to a binary buffer.

    78         if str_instance and hex_regex.match(geom_input):

    79             geom_input = buffer(a2b_hex(geom_input.upper()))

    80             str_instance = False

    81 

    82         # Constructing the geometry, 

    83         if str_instance:

    84             # Checking if unicode

    85             if isinstance(geom_input, UnicodeType):

    86                 # Encoding to ASCII, WKT or HEX doesn't need any more.

    87                 geo_input = geo_input.encode('ascii')

    88 

    89             wkt_m = wkt_regex.match(geom_input)

    90             json_m = json_regex.match(geom_input)

    91             if wkt_m:

    92                 if wkt_m.group('type').upper() == 'LINEARRING':

    93                     # OGR_G_CreateFromWkt doesn't work with LINEARRING WKT.

    94                     #  See http://trac.osgeo.org/gdal/ticket/1992.

    95                     g = create_geom(OGRGeomType(wkt_m.group('type')).num)

    96                     import_wkt(g, byref(c_char_p(geom_input)))

    97                 else:

    98                     g = from_wkt(byref(c_char_p(geom_input)), None, byref(c_void_p()))

    99             elif json_m:

   100                 if GEOJSON:

   101                     g = from_json(geom_input)

   102                 else:

   103                     raise NotImplementedError('GeoJSON input only supported on GDAL 1.5+.')

   104             else:

   105                 # Seeing if the input is a valid short-hand string

   106                 # (e.g., 'Point', 'POLYGON').

   107                 ogr_t = OGRGeomType(geom_input)

   108                 g = create_geom(OGRGeomType(geom_input).num)

   109         elif isinstance(geom_input, buffer):

   110             # WKB was passed in

   111             g = from_wkb(str(geom_input), None, byref(c_void_p()), len(geom_input))

   112         elif isinstance(geom_input, OGRGeomType):

   113             # OGRGeomType was passed in, an empty geometry will be created.

   114             g = create_geom(geom_input.num)

   115         elif isinstance(geom_input, c_void_p):

   116             # OGR pointer (c_void_p) was the input.

   117             g = geom_input

   118         else:

   119             raise OGRException('Invalid input type for OGR Geometry construction: %s' % type(geom_input))

   120 

   121         # Now checking the Geometry pointer before finishing initialization

   122         # by setting the pointer for the object.

   123         if not g:

   124             raise OGRException('Cannot create OGR Geometry from input: %s' % str(geom_input))

   125         self._ptr = g

   126 

   127         # Assigning the SpatialReference object to the geometry, if valid.

   128         if bool(srs): self.srs = srs

   129 

   130         # Setting the class depending upon the OGR Geometry Type

   131         self.__class__ = GEO_CLASSES[self.geom_type.num]

   132 

   133     def __del__(self):

   134         "Deletes this Geometry."

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

   136 

   137     ### Geometry set-like operations ###

   138     # g = g1 | g2

   139     def __or__(self, other):

   140         "Returns the union of the two geometries."

   141         return self.union(other)

   142 

   143     # g = g1 & g2

   144     def __and__(self, other):

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

   146         return self.intersection(other)

   147 

   148     # g = g1 - g2

   149     def __sub__(self, other):

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

   151         return self.difference(other)

   152 

   153     # g = g1 ^ g2

   154     def __xor__(self, other):

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

   156         return self.sym_difference(other)

   157 

   158     def __eq__(self, other):

   159         "Is this Geometry equal to the other?"

   160         return self.equals(other)

   161 

   162     def __ne__(self, other):

   163         "Tests for inequality."

   164         return not self.equals(other)

   165 

   166     def __str__(self):

   167         "WKT is used for the string representation."

   168         return self.wkt

   169 

   170     #### Geometry Properties ####

   171     @property

   172     def dimension(self):

   173         "Returns 0 for points, 1 for lines, and 2 for surfaces."

   174         return get_dims(self._ptr)

   175 

   176     @property

   177     def coord_dim(self):

   178         "Returns the coordinate dimension of the Geometry."

   179         return get_coord_dims(self._ptr)

   180 

   181     @property

   182     def geom_count(self):

   183         "The number of elements in this Geometry."

   184         return get_geom_count(self._ptr)

   185 

   186     @property

   187     def point_count(self):

   188         "Returns the number of Points in this Geometry."

   189         return get_point_count(self._ptr)

   190 

   191     @property

   192     def num_points(self):

   193         "Alias for `point_count` (same name method in GEOS API.)"

   194         return self.point_count

   195 

   196     @property

   197     def num_coords(self):

   198         "Alais for `point_count`."

   199         return self.point_count

   200 

   201     @property

   202     def geom_type(self):

   203         "Returns the Type for this Geometry."

   204         try:

   205             return OGRGeomType(get_geom_type(self._ptr))

   206         except OGRException:

   207             # VRT datasources return an invalid geometry type

   208             # number, but a valid name -- we'll try that instead.

   209             # See: http://trac.osgeo.org/gdal/ticket/2491

   210             return OGRGeomType(get_geom_name(self._ptr))

   211 

   212     @property

   213     def geom_name(self):

   214         "Returns the Name of this Geometry."

   215         return get_geom_name(self._ptr)

   216 

   217     @property

   218     def area(self):

   219         "Returns the area for a LinearRing, Polygon, or MultiPolygon; 0 otherwise."

   220         return get_area(self._ptr)

   221 

   222     @property

   223     def envelope(self):

   224         "Returns the envelope for this Geometry."

   225         # TODO: Fix Envelope() for Point geometries.

   226         return Envelope(get_envelope(self._ptr, byref(OGREnvelope())))

   227 

   228     @property

   229     def extent(self):

   230         "Returns the envelope as a 4-tuple, instead of as an Envelope object."

   231         return self.envelope.tuple

   232 

   233     #### SpatialReference-related Properties ####

   234     

   235     # The SRS property

   236     def get_srs(self):

   237         "Returns the Spatial Reference for this Geometry."

   238         try:

   239             srs_ptr = get_geom_srs(self._ptr)

   240             return SpatialReference(clone_srs(srs_ptr))

   241         except SRSException:

   242             return None

   243 

   244     def set_srs(self, srs):

   245         "Sets the SpatialReference for this geometry."

   246         if isinstance(srs, SpatialReference):

   247             srs_ptr = clone_srs(srs._ptr)

   248         elif isinstance(srs, (int, long, basestring)):

   249             sr = SpatialReference(srs)

   250             srs_ptr = clone_srs(sr._ptr)

   251         else:

   252             raise TypeError('Cannot assign spatial reference with object of type: %s' % type(srs))

   253         assign_srs(self._ptr, srs_ptr)

   254 

   255     srs = property(get_srs, set_srs)

   256 

   257     # The SRID property

   258     def get_srid(self):

   259         if self.srs: return self.srs.srid

   260         else: return None

   261 

   262     def set_srid(self, srid):

   263         if isinstance(srid, (int, long)):

   264             self.srs = srid

   265         else:

   266             raise TypeError('SRID must be set with an integer.')

   267 

   268     srid = property(get_srid, set_srid)

   269 

   270     #### Output Methods ####

   271     @property

   272     def geos(self):

   273         "Returns a GEOSGeometry object from this OGRGeometry."

   274         from django.contrib.gis.geos import GEOSGeometry

   275         return GEOSGeometry(self.wkb, self.srid)

   276 

   277     @property

   278     def gml(self):

   279         "Returns the GML representation of the Geometry."

   280         return to_gml(self._ptr)

   281 

   282     @property

   283     def hex(self):

   284         "Returns the hexadecimal representation of the WKB (a string)."

   285         return str(self.wkb).encode('hex').upper()

   286         #return b2a_hex(self.wkb).upper()

   287 

   288     @property

   289     def json(self):

   290         if GEOJSON: 

   291             return to_json(self._ptr)

   292         else:

   293             raise NotImplementedError('GeoJSON output only supported on GDAL 1.5+.')

   294     geojson = json

   295 

   296     @property

   297     def wkb_size(self):

   298         "Returns the size of the WKB buffer."

   299         return get_wkbsize(self._ptr)

   300 

   301     @property

   302     def wkb(self):

   303         "Returns the WKB representation of the Geometry."

   304         if sys.byteorder == 'little':

   305             byteorder = 1 # wkbNDR (from ogr_core.h)

   306         else:

   307             byteorder = 0 # wkbXDR

   308         sz = self.wkb_size

   309         # Creating the unsigned character buffer, and passing it in by reference.

   310         buf = (c_ubyte * sz)()

   311         wkb = to_wkb(self._ptr, byteorder, byref(buf))

   312         # Returning a buffer of the string at the pointer.

   313         return buffer(string_at(buf, sz))

   314 

   315     @property

   316     def wkt(self):

   317         "Returns the WKT representation of the Geometry."

   318         return to_wkt(self._ptr, byref(c_char_p()))

   319     

   320     #### Geometry Methods ####

   321     def clone(self):

   322         "Clones this OGR Geometry."

   323         return OGRGeometry(clone_geom(self._ptr), self.srs)

   324 

   325     def close_rings(self):

   326         """

   327         If there are any rings within this geometry that have not been

   328         closed, this routine will do so by adding the starting point at the

   329         end.

   330         """

   331         # Closing the open rings.

   332         geom_close_rings(self._ptr)

   333 

   334     def transform(self, coord_trans, clone=False):

   335         """

   336         Transforms this geometry to a different spatial reference system.

   337         May take a CoordTransform object, a SpatialReference object, string

   338         WKT or PROJ.4, and/or an integer SRID.  By default nothing is returned

   339         and the geometry is transformed in-place.  However, if the `clone`

   340         keyword is set, then a transformed clone of this geometry will be

   341         returned.

   342         """

   343         if clone:

   344             klone = self.clone()

   345             klone.transform(coord_trans)

   346             return klone

   347         if isinstance(coord_trans, CoordTransform):

   348             geom_transform(self._ptr, coord_trans._ptr)

   349         elif isinstance(coord_trans, SpatialReference):

   350             geom_transform_to(self._ptr, coord_trans._ptr)

   351         elif isinstance(coord_trans, (int, long, basestring)):

   352             sr = SpatialReference(coord_trans)

   353             geom_transform_to(self._ptr, sr._ptr)

   354         else:

   355             raise TypeError('Transform only accepts CoordTransform, SpatialReference, string, and integer objects.')

   356 

   357     def transform_to(self, srs):

   358         "For backwards-compatibility."

   359         self.transform(srs)

   360 

   361     #### Topology Methods ####

   362     def _topology(self, func, other):

   363         """A generalized function for topology operations, takes a GDAL function and

   364         the other geometry to perform the operation on."""

   365         if not isinstance(other, OGRGeometry):

   366             raise TypeError('Must use another OGRGeometry object for topology operations!')

   367 

   368         # Returning the output of the given function with the other geometry's

   369         # pointer.

   370         return func(self._ptr, other._ptr)

   371 

   372     def intersects(self, other):

   373         "Returns True if this geometry intersects with the other."

   374         return self._topology(ogr_intersects, other)

   375     

   376     def equals(self, other):

   377         "Returns True if this geometry is equivalent to the other."

   378         return self._topology(ogr_equals, other)

   379 

   380     def disjoint(self, other):

   381         "Returns True if this geometry and the other are spatially disjoint."

   382         return self._topology(ogr_disjoint, other)

   383 

   384     def touches(self, other):

   385         "Returns True if this geometry touches the other."

   386         return self._topology(ogr_touches, other)

   387 

   388     def crosses(self, other):

   389         "Returns True if this geometry crosses the other."

   390         return self._topology(ogr_crosses, other)

   391 

   392     def within(self, other):

   393         "Returns True if this geometry is within the other."

   394         return self._topology(ogr_within, other)

   395 

   396     def contains(self, other):

   397         "Returns True if this geometry contains the other."

   398         return self._topology(ogr_contains, other)

   399 

   400     def overlaps(self, other):

   401         "Returns True if this geometry overlaps the other."

   402         return self._topology(ogr_overlaps, other)

   403 

   404     #### Geometry-generation Methods ####

   405     def _geomgen(self, gen_func, other=None):

   406         "A helper routine for the OGR routines that generate geometries."

   407         if isinstance(other, OGRGeometry):

   408             return OGRGeometry(gen_func(self._ptr, other._ptr), self.srs)

   409         else:

   410             return OGRGeometry(gen_func(self._ptr), self.srs)

   411 

   412     @property

   413     def boundary(self):

   414         "Returns the boundary of this geometry."

   415         return self._geomgen(get_boundary)

   416 

   417     @property

   418     def convex_hull(self):

   419         """

   420         Returns the smallest convex Polygon that contains all the points in 

   421         this Geometry.

   422         """

   423         return self._geomgen(geom_convex_hull)

   424 

   425     def difference(self, other):

   426         """

   427         Returns a new geometry consisting of the region which is the difference

   428         of this geometry and the other.

   429         """

   430         return self._geomgen(geom_diff, other)

   431 

   432     def intersection(self, other):

   433         """

   434         Returns a new geometry consisting of the region of intersection of this

   435         geometry and the other.

   436         """

   437         return self._geomgen(geom_intersection, other)

   438 

   439     def sym_difference(self, other):

   440         """                                                                                                                                                

   441         Returns a new geometry which is the symmetric difference of this

   442         geometry and the other.

   443         """

   444         return self._geomgen(geom_sym_diff, other)

   445 

   446     def union(self, other):

   447         """

   448         Returns a new geometry consisting of the region which is the union of

   449         this geometry and the other.

   450         """

   451         return self._geomgen(geom_union, other)

   452 

   453 # The subclasses for OGR Geometry.

   454 class Point(OGRGeometry):

   455 

   456     @property

   457     def x(self):

   458         "Returns the X coordinate for this Point."

   459         return getx(self._ptr, 0)

   460 

   461     @property

   462     def y(self):

   463         "Returns the Y coordinate for this Point."

   464         return gety(self._ptr, 0)

   465 

   466     @property

   467     def z(self):

   468         "Returns the Z coordinate for this Point."

   469         if self.coord_dim == 3:

   470             return getz(self._ptr, 0)

   471 

   472     @property

   473     def tuple(self):

   474         "Returns the tuple of this point."

   475         if self.coord_dim == 2:

   476             return (self.x, self.y)

   477         elif self.coord_dim == 3:

   478             return (self.x, self.y, self.z)

   479     coords = tuple

   480 

   481 class LineString(OGRGeometry):

   482 

   483     def __getitem__(self, index):

   484         "Returns the Point at the given index."

   485         if index >= 0 and index < self.point_count:

   486             x, y, z = c_double(), c_double(), c_double()

   487             get_point(self._ptr, index, byref(x), byref(y), byref(z))

   488             dim = self.coord_dim

   489             if dim == 1:

   490                 return (x.value,)

   491             elif dim == 2:

   492                 return (x.value, y.value)

   493             elif dim == 3:

   494                 return (x.value, y.value, z.value)

   495         else:

   496             raise OGRIndexError('index out of range: %s' % str(index))

   497 

   498     def __iter__(self):

   499         "Iterates over each point in the LineString."

   500         for i in xrange(self.point_count):

   501             yield self[i]

   502 

   503     def __len__(self):

   504         "The length returns the number of points in the LineString."

   505         return self.point_count

   506 

   507     @property

   508     def tuple(self):

   509         "Returns the tuple representation of this LineString."

   510         return tuple([self[i] for i in xrange(len(self))])

   511     coords = tuple

   512 

   513     def _listarr(self, func):

   514         """

   515         Internal routine that returns a sequence (list) corresponding with

   516         the given function.

   517         """

   518         return [func(self._ptr, i) for i in xrange(len(self))]

   519 

   520     @property

   521     def x(self):

   522         "Returns the X coordinates in a list."

   523         return self._listarr(getx)

   524 

   525     @property

   526     def y(self):

   527         "Returns the Y coordinates in a list."

   528         return self._listarr(gety)

   529     

   530     @property

   531     def z(self):

   532         "Returns the Z coordinates in a list."

   533         if self.coord_dim == 3:

   534             return self._listarr(getz)

   535 

   536 # LinearRings are used in Polygons.

   537 class LinearRing(LineString): pass

   538 

   539 class Polygon(OGRGeometry):

   540 

   541     def __len__(self):

   542         "The number of interior rings in this Polygon."

   543         return self.geom_count

   544 

   545     def __iter__(self):

   546         "Iterates through each ring in the Polygon."

   547         for i in xrange(self.geom_count):

   548             yield self[i]

   549 

   550     def __getitem__(self, index):

   551         "Gets the ring at the specified index."

   552         if index < 0 or index >= self.geom_count:

   553             raise OGRIndexError('index out of range: %s' % index)

   554         else:

   555             return OGRGeometry(clone_geom(get_geom_ref(self._ptr, index)), self.srs)

   556 

   557     # Polygon Properties

   558     @property

   559     def shell(self):

   560         "Returns the shell of this Polygon."

   561         return self[0] # First ring is the shell

   562     exterior_ring = shell

   563 

   564     @property

   565     def tuple(self):

   566         "Returns a tuple of LinearRing coordinate tuples."

   567         return tuple([self[i].tuple for i in xrange(self.geom_count)])

   568     coords = tuple

   569 

   570     @property

   571     def point_count(self):

   572         "The number of Points in this Polygon."

   573         # Summing up the number of points in each ring of the Polygon.

   574         return sum([self[i].point_count for i in xrange(self.geom_count)])

   575 

   576     @property

   577     def centroid(self):

   578         "Returns the centroid (a Point) of this Polygon."

   579         # The centroid is a Point, create a geometry for this.

   580         p = OGRGeometry(OGRGeomType('Point'))

   581         get_centroid(self._ptr, p._ptr)

   582         return p

   583 

   584 # Geometry Collection base class.

   585 class GeometryCollection(OGRGeometry):

   586     "The Geometry Collection class."

   587 

   588     def __getitem__(self, index):

   589         "Gets the Geometry at the specified index."

   590         if index < 0 or index >= self.geom_count:

   591             raise OGRIndexError('index out of range: %s' % index)

   592         else:

   593             return OGRGeometry(clone_geom(get_geom_ref(self._ptr, index)), self.srs)

   594         

   595     def __iter__(self):

   596         "Iterates over each Geometry."

   597         for i in xrange(self.geom_count):

   598             yield self[i]

   599 

   600     def __len__(self):

   601         "The number of geometries in this Geometry Collection."

   602         return self.geom_count

   603 

   604     def add(self, geom):

   605         "Add the geometry to this Geometry Collection."

   606         if isinstance(geom, OGRGeometry):

   607             if isinstance(geom, self.__class__):

   608                 for g in geom: add_geom(self._ptr, g._ptr)

   609             else:

   610                 add_geom(self._ptr, geom._ptr)

   611         elif isinstance(geom, basestring):

   612             tmp = OGRGeometry(geom)

   613             add_geom(self._ptr, tmp._ptr)

   614         else:

   615             raise OGRException('Must add an OGRGeometry.')

   616 

   617     @property

   618     def point_count(self):

   619         "The number of Points in this Geometry Collection."

   620         # Summing up the number of points in each geometry in this collection

   621         return sum([self[i].point_count for i in xrange(self.geom_count)])

   622 

   623     @property

   624     def tuple(self):

   625         "Returns a tuple representation of this Geometry Collection."

   626         return tuple([self[i].tuple for i in xrange(self.geom_count)])

   627     coords = tuple

   628 

   629 # Multiple Geometry types.

   630 class MultiPoint(GeometryCollection): pass

   631 class MultiLineString(GeometryCollection): pass

   632 class MultiPolygon(GeometryCollection): pass

   633 

   634 # Class mapping dictionary (using the OGRwkbGeometryType as the key)

   635 GEO_CLASSES = {1 : Point,

   636                2 : LineString,

   637                3 : Polygon,

   638                4 : MultiPoint,

   639                5 : MultiLineString,

   640                6 : MultiPolygon,

   641                7 : GeometryCollection,

   642                101: LinearRing, 

   643                }