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 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 

    46 # Getting GDAL prerequisites

    47 from django.contrib.gis.gdal.base import GDALBase

    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.libgdal import GEOJSON, GDAL_VERSION

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

    53 

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

    55 from django.contrib.gis.gdal.prototypes import geom as capi, srs as srs_api

    56 

    57 # For recognizing geometry input.

    58 from django.contrib.gis.geometry.regex import hex_regex, wkt_regex, json_regex

    59 

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

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

    62 #

    63 # The OGR_G_* routines are relevant here.

    64 

    65 #### OGRGeometry Class ####

    66 class OGRGeometry(GDALBase):

    67     "Generally encapsulates an OGR geometry."

    68 

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

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

    71 

    72         str_instance = isinstance(geom_input, basestring)

    73 

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

    75         if str_instance and hex_regex.match(geom_input):

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

    77             str_instance = False

    78 

    79         # Constructing the geometry,

    80         if str_instance:

    81             # Checking if unicode

    82             if isinstance(geom_input, unicode):

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

    84                 geom_input = geom_input.encode('ascii')

    85 

    86             wkt_m = wkt_regex.match(geom_input)

    87             json_m = json_regex.match(geom_input)

    88             if wkt_m:

    89                 if wkt_m.group('srid'):

    90                     # If there's EWKT, set the SRS w/value of the SRID.

    91                     srs = int(wkt_m.group('srid'))

    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 = capi.create_geom(OGRGeomType(wkt_m.group('type')).num)

    96                     capi.import_wkt(g, byref(c_char_p(wkt_m.group('wkt'))))

    97                 else:

    98                     g = capi.from_wkt(byref(c_char_p(wkt_m.group('wkt'))), None, byref(c_void_p()))

    99             elif json_m:

   100                 if GEOJSON:

   101                     g = capi.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 = capi.create_geom(OGRGeomType(geom_input).num)

   109         elif isinstance(geom_input, buffer):

   110             # WKB was passed in

   111             g = capi.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 = capi.create_geom(geom_input.num)

   115         elif isinstance(geom_input, self.ptr_type):

   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: capi.destroy_geom(self._ptr)

   136 

   137     # Pickle routines

   138     def __getstate__(self):

   139         srs = self.srs

   140         if srs:

   141             srs = srs.wkt

   142         else:

   143             srs = None

   144         return str(self.wkb), srs

   145 

   146     def __setstate__(self, state):

   147         wkb, srs = state

   148         ptr = capi.from_wkb(wkb, None, byref(c_void_p()), len(wkb))

   149         if not ptr: raise OGRException('Invalid OGRGeometry loaded from pickled state.')

   150         self.ptr = ptr

   151         self.srs = srs

   152 

   153     @classmethod

   154     def from_bbox(cls, bbox):

   155         "Constructs a Polygon from a bounding box (4-tuple)."

   156         x0, y0, x1, y1 = bbox

   157         return OGRGeometry( 'POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))' %  (

   158                 x0, y0, x0, y1, x1, y1, x1, y0, x0, y0) )

   159 

   160     ### Geometry set-like operations ###

   161     # g = g1 | g2

   162     def __or__(self, other):

   163         "Returns the union of the two geometries."

   164         return self.union(other)

   165 

   166     # g = g1 & g2

   167     def __and__(self, other):

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

   169         return self.intersection(other)

   170 

   171     # g = g1 - g2

   172     def __sub__(self, other):

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

   174         return self.difference(other)

   175 

   176     # g = g1 ^ g2

   177     def __xor__(self, other):

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

   179         return self.sym_difference(other)

   180 

   181     def __eq__(self, other):

   182         "Is this Geometry equal to the other?"

   183         if isinstance(other, OGRGeometry):

   184             return self.equals(other)

   185         else:

   186             return False

   187 

   188     def __ne__(self, other):

   189         "Tests for inequality."

   190         return not (self == other)

   191 

   192     def __str__(self):

   193         "WKT is used for the string representation."

   194         return self.wkt

   195 

   196     #### Geometry Properties ####

   197     @property

   198     def dimension(self):

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

   200         return capi.get_dims(self.ptr)

   201 

   202     def _get_coord_dim(self):

   203         "Returns the coordinate dimension of the Geometry."

   204         if isinstance(self, GeometryCollection) and GDAL_VERSION < (1, 5, 2):

   205             # On GDAL versions prior to 1.5.2, there exists a bug in which

   206             # the coordinate dimension of geometry collections is always 2:

   207             #   http://trac.osgeo.org/gdal/ticket/2334

   208             # Here we workaround by returning the coordinate dimension of the

   209             # first geometry in the collection instead.

   210             if len(self):

   211                 return capi.get_coord_dim(capi.get_geom_ref(self.ptr, 0))

   212         return capi.get_coord_dim(self.ptr)

   213 

   214     def _set_coord_dim(self, dim):

   215         "Sets the coordinate dimension of this Geometry."

   216         if not dim in (2, 3):

   217             raise ValueError('Geometry dimension must be either 2 or 3')

   218         capi.set_coord_dim(self.ptr, dim)

   219 

   220     coord_dim = property(_get_coord_dim, _set_coord_dim)

   221 

   222     @property

   223     def geom_count(self):

   224         "The number of elements in this Geometry."

   225         return capi.get_geom_count(self.ptr)

   226 

   227     @property

   228     def point_count(self):

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

   230         return capi.get_point_count(self.ptr)

   231 

   232     @property

   233     def num_points(self):

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

   235         return self.point_count

   236 

   237     @property

   238     def num_coords(self):

   239         "Alais for `point_count`."

   240         return self.point_count

   241 

   242     @property

   243     def geom_type(self):

   244         "Returns the Type for this Geometry."

   245         return OGRGeomType(capi.get_geom_type(self.ptr))

   246 

   247     @property

   248     def geom_name(self):

   249         "Returns the Name of this Geometry."

   250         return capi.get_geom_name(self.ptr)

   251 

   252     @property

   253     def area(self):

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

   255         return capi.get_area(self.ptr)

   256 

   257     @property

   258     def envelope(self):

   259         "Returns the envelope for this Geometry."

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

   261         return Envelope(capi.get_envelope(self.ptr, byref(OGREnvelope())))

   262 

   263     @property

   264     def extent(self):

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

   266         return self.envelope.tuple

   267 

   268     #### SpatialReference-related Properties ####

   269 

   270     # The SRS property

   271     def _get_srs(self):

   272         "Returns the Spatial Reference for this Geometry."

   273         try:

   274             srs_ptr = capi.get_geom_srs(self.ptr)

   275             return SpatialReference(srs_api.clone_srs(srs_ptr))

   276         except SRSException:

   277             return None

   278 

   279     def _set_srs(self, srs):

   280         "Sets the SpatialReference for this geometry."

   281         # Do not have to clone the `SpatialReference` object pointer because

   282         # when it is assigned to this `OGRGeometry` it's internal OGR

   283         # reference count is incremented, and will likewise be released

   284         # (decremented) when this geometry's destructor is called.

   285         if isinstance(srs, SpatialReference):

   286             srs_ptr = srs.ptr

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

   288             sr = SpatialReference(srs)

   289             srs_ptr = sr.ptr

   290         else:

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

   292         capi.assign_srs(self.ptr, srs_ptr)

   293 

   294     srs = property(_get_srs, _set_srs)

   295 

   296     # The SRID property

   297     def _get_srid(self):

   298         srs = self.srs

   299         if srs: return srs.srid

   300         return None

   301 

   302     def _set_srid(self, srid):

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

   304             self.srs = srid

   305         else:

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

   307 

   308     srid = property(_get_srid, _set_srid)

   309 

   310     #### Output Methods ####

   311     @property

   312     def geos(self):

   313         "Returns a GEOSGeometry object from this OGRGeometry."

   314         from django.contrib.gis.geos import GEOSGeometry

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

   316 

   317     @property

   318     def gml(self):

   319         "Returns the GML representation of the Geometry."

   320         return capi.to_gml(self.ptr)

   321 

   322     @property

   323     def hex(self):

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

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

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

   327 

   328     @property

   329     def json(self):

   330         """

   331         Returns the GeoJSON representation of this Geometry (requires

   332         GDAL 1.5+).

   333         """

   334         if GEOJSON:

   335             return capi.to_json(self.ptr)

   336         else:

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

   338     geojson = json

   339 

   340     @property

   341     def kml(self):

   342         "Returns the KML representation of the Geometry."

   343         if GEOJSON:

   344             return capi.to_kml(self.ptr, None)

   345         else:

   346             raise NotImplementedError('KML output only supported on GDAL 1.5+.')

   347 

   348     @property

   349     def wkb_size(self):

   350         "Returns the size of the WKB buffer."

   351         return capi.get_wkbsize(self.ptr)

   352 

   353     @property

   354     def wkb(self):

   355         "Returns the WKB representation of the Geometry."

   356         if sys.byteorder == 'little':

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

   358         else:

   359             byteorder = 0 # wkbXDR

   360         sz = self.wkb_size

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

   362         buf = (c_ubyte * sz)()

   363         wkb = capi.to_wkb(self.ptr, byteorder, byref(buf))

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

   365         return buffer(string_at(buf, sz))

   366 

   367     @property

   368     def wkt(self):

   369         "Returns the WKT representation of the Geometry."

   370         return capi.to_wkt(self.ptr, byref(c_char_p()))

   371 

   372     @property

   373     def ewkt(self):

   374         "Returns the EWKT representation of the Geometry."

   375         srs = self.srs

   376         if srs and srs.srid:

   377             return 'SRID=%s;%s' % (srs.srid, self.wkt)

   378         else:

   379             return self.wkt

   380 

   381     #### Geometry Methods ####

   382     def clone(self):

   383         "Clones this OGR Geometry."

   384         return OGRGeometry(capi.clone_geom(self.ptr), self.srs)

   385 

   386     def close_rings(self):

   387         """

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

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

   390         end.

   391         """

   392         # Closing the open rings.

   393         capi.geom_close_rings(self.ptr)

   394 

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

   396         """

   397         Transforms this geometry to a different spatial reference system.

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

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

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

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

   402         returned.

   403         """

   404         if clone:

   405             klone = self.clone()

   406             klone.transform(coord_trans)

   407             return klone

   408 

   409         # Have to get the coordinate dimension of the original geometry

   410         # so it can be used to reset the transformed geometry's dimension

   411         # afterwards.  This is done because of GDAL bug (in versions prior

   412         # to 1.7) that turns geometries 3D after transformation, see:

   413         #  http://trac.osgeo.org/gdal/changeset/17792

   414         if GDAL_VERSION < (1, 7):

   415             orig_dim = self.coord_dim

   416 

   417         # Depending on the input type, use the appropriate OGR routine

   418         # to perform the transformation.

   419         if isinstance(coord_trans, CoordTransform):

   420             capi.geom_transform(self.ptr, coord_trans.ptr)

   421         elif isinstance(coord_trans, SpatialReference):

   422             capi.geom_transform_to(self.ptr, coord_trans.ptr)

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

   424             sr = SpatialReference(coord_trans)

   425             capi.geom_transform_to(self.ptr, sr.ptr)

   426         else:

   427             raise TypeError('Transform only accepts CoordTransform, '

   428                             'SpatialReference, string, and integer objects.')

   429 

   430         # Setting with original dimension, see comment above.

   431         if GDAL_VERSION < (1, 7):

   432             if isinstance(self, GeometryCollection):

   433                 # With geometry collections have to set dimension on

   434                 # each internal geometry reference, as the collection

   435                 # dimension isn't affected.

   436                 for i in xrange(len(self)):

   437                     internal_ptr = capi.get_geom_ref(self.ptr, i)

   438                     if orig_dim != capi.get_coord_dim(internal_ptr):

   439                         capi.set_coord_dim(internal_ptr, orig_dim)

   440             else:

   441                 if self.coord_dim != orig_dim:

   442                     self.coord_dim = orig_dim

   443 

   444     def transform_to(self, srs):

   445         "For backwards-compatibility."

   446         self.transform(srs)

   447 

   448     #### Topology Methods ####

   449     def _topology(self, func, other):

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

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

   452         if not isinstance(other, OGRGeometry):

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

   454 

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

   456         # pointer.

   457         return func(self.ptr, other.ptr)

   458 

   459     def intersects(self, other):

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

   461         return self._topology(capi.ogr_intersects, other)

   462 

   463     def equals(self, other):

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

   465         return self._topology(capi.ogr_equals, other)

   466 

   467     def disjoint(self, other):

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

   469         return self._topology(capi.ogr_disjoint, other)

   470 

   471     def touches(self, other):

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

   473         return self._topology(capi.ogr_touches, other)

   474 

   475     def crosses(self, other):

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

   477         return self._topology(capi.ogr_crosses, other)

   478 

   479     def within(self, other):

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

   481         return self._topology(capi.ogr_within, other)

   482 

   483     def contains(self, other):

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

   485         return self._topology(capi.ogr_contains, other)

   486 

   487     def overlaps(self, other):

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

   489         return self._topology(capi.ogr_overlaps, other)

   490 

   491     #### Geometry-generation Methods ####

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

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

   494         if isinstance(other, OGRGeometry):

   495             return OGRGeometry(gen_func(self.ptr, other.ptr), self.srs)

   496         else:

   497             return OGRGeometry(gen_func(self.ptr), self.srs)

   498 

   499     @property

   500     def boundary(self):

   501         "Returns the boundary of this geometry."

   502         return self._geomgen(capi.get_boundary)

   503 

   504     @property

   505     def convex_hull(self):

   506         """

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

   508         this Geometry.

   509         """

   510         return self._geomgen(capi.geom_convex_hull)

   511 

   512     def difference(self, other):

   513         """

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

   515         of this geometry and the other.

   516         """

   517         return self._geomgen(capi.geom_diff, other)

   518 

   519     def intersection(self, other):

   520         """

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

   522         geometry and the other.

   523         """

   524         return self._geomgen(capi.geom_intersection, other)

   525 

   526     def sym_difference(self, other):

   527         """

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

   529         geometry and the other.

   530         """

   531         return self._geomgen(capi.geom_sym_diff, other)

   532 

   533     def union(self, other):

   534         """

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

   536         this geometry and the other.

   537         """

   538         return self._geomgen(capi.geom_union, other)

   539 

   540 # The subclasses for OGR Geometry.

   541 class Point(OGRGeometry):

   542 

   543     @property

   544     def x(self):

   545         "Returns the X coordinate for this Point."

   546         return capi.getx(self.ptr, 0)

   547 

   548     @property

   549     def y(self):

   550         "Returns the Y coordinate for this Point."

   551         return capi.gety(self.ptr, 0)

   552 

   553     @property

   554     def z(self):

   555         "Returns the Z coordinate for this Point."

   556         if self.coord_dim == 3:

   557             return capi.getz(self.ptr, 0)

   558 

   559     @property

   560     def tuple(self):

   561         "Returns the tuple of this point."

   562         if self.coord_dim == 2:

   563             return (self.x, self.y)

   564         elif self.coord_dim == 3:

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

   566     coords = tuple

   567 

   568 class LineString(OGRGeometry):

   569 

   570     def __getitem__(self, index):

   571         "Returns the Point at the given index."

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

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

   574             capi.get_point(self.ptr, index, byref(x), byref(y), byref(z))

   575             dim = self.coord_dim

   576             if dim == 1:

   577                 return (x.value,)

   578             elif dim == 2:

   579                 return (x.value, y.value)

   580             elif dim == 3:

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

   582         else:

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

   584 

   585     def __iter__(self):

   586         "Iterates over each point in the LineString."

   587         for i in xrange(self.point_count):

   588             yield self[i]

   589 

   590     def __len__(self):

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

   592         return self.point_count

   593 

   594     @property

   595     def tuple(self):

   596         "Returns the tuple representation of this LineString."

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

   598     coords = tuple

   599 

   600     def _listarr(self, func):

   601         """

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

   603         the given function.

   604         """

   605         return [func(self.ptr, i) for i in xrange(len(self))]

   606 

   607     @property

   608     def x(self):

   609         "Returns the X coordinates in a list."

   610         return self._listarr(capi.getx)

   611 

   612     @property

   613     def y(self):

   614         "Returns the Y coordinates in a list."

   615         return self._listarr(capi.gety)

   616 

   617     @property

   618     def z(self):

   619         "Returns the Z coordinates in a list."

   620         if self.coord_dim == 3:

   621             return self._listarr(capi.getz)

   622 

   623 # LinearRings are used in Polygons.

   624 class LinearRing(LineString): pass

   625 

   626 class Polygon(OGRGeometry):

   627 

   628     def __len__(self):

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

   630         return self.geom_count

   631 

   632     def __iter__(self):

   633         "Iterates through each ring in the Polygon."

   634         for i in xrange(self.geom_count):

   635             yield self[i]

   636 

   637     def __getitem__(self, index):

   638         "Gets the ring at the specified index."

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

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

   641         else:

   642             return OGRGeometry(capi.clone_geom(capi.get_geom_ref(self.ptr, index)), self.srs)

   643 

   644     # Polygon Properties

   645     @property

   646     def shell(self):

   647         "Returns the shell of this Polygon."

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

   649     exterior_ring = shell

   650 

   651     @property

   652     def tuple(self):

   653         "Returns a tuple of LinearRing coordinate tuples."

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

   655     coords = tuple

   656 

   657     @property

   658     def point_count(self):

   659         "The number of Points in this Polygon."

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

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

   662 

   663     @property

   664     def centroid(self):

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

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

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

   668         capi.get_centroid(self.ptr, p.ptr)

   669         return p

   670 

   671 # Geometry Collection base class.

   672 class GeometryCollection(OGRGeometry):

   673     "The Geometry Collection class."

   674 

   675     def __getitem__(self, index):

   676         "Gets the Geometry at the specified index."

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

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

   679         else:

   680             return OGRGeometry(capi.clone_geom(capi.get_geom_ref(self.ptr, index)), self.srs)

   681 

   682     def __iter__(self):

   683         "Iterates over each Geometry."

   684         for i in xrange(self.geom_count):

   685             yield self[i]

   686 

   687     def __len__(self):

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

   689         return self.geom_count

   690 

   691     def add(self, geom):

   692         "Add the geometry to this Geometry Collection."

   693         if isinstance(geom, OGRGeometry):

   694             if isinstance(geom, self.__class__):

   695                 for g in geom: capi.add_geom(self.ptr, g.ptr)

   696             else:

   697                 capi.add_geom(self.ptr, geom.ptr)

   698         elif isinstance(geom, basestring):

   699             tmp = OGRGeometry(geom)

   700             capi.add_geom(self.ptr, tmp.ptr)

   701         else:

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

   703 

   704     @property

   705     def point_count(self):

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

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

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

   709 

   710     @property

   711     def tuple(self):

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

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

   714     coords = tuple

   715 

   716 # Multiple Geometry types.

   717 class MultiPoint(GeometryCollection): pass

   718 class MultiLineString(GeometryCollection): pass

   719 class MultiPolygon(GeometryCollection): pass

   720 

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

   722 GEO_CLASSES = {1 : Point,

   723                2 : LineString,

   724                3 : Polygon,

   725                4 : MultiPoint,

   726                5 : MultiLineString,

   727                6 : MultiPolygon,

   728                7 : GeometryCollection,

   729                101: LinearRing,

   730                1 + OGRGeomType.wkb25bit : Point,

   731                2 + OGRGeomType.wkb25bit : LineString,

   732                3 + OGRGeomType.wkb25bit : Polygon,

   733                4 + OGRGeomType.wkb25bit : MultiPoint,

   734                5 + OGRGeomType.wkb25bit : MultiLineString,

   735                6 + OGRGeomType.wkb25bit : MultiPolygon,

   736                7 + OGRGeomType.wkb25bit : GeometryCollection,

   737                }