Django

Code

Changeset 5637

Show
Ignore:
Timestamp:
07/09/07 01:16:06 (1 year ago)
Author:
jbronn
Message:

gis: GEOS library has been re-factored to, changes include:

(1) allow write-access to Geometries (for future lazy-geometries)
(2) Point, LineString?, LinearRing? have their own constructors (e.g., p = Point(5, 17))
(3) improved ctypes memory management (as part of writability enhancements)
(4) improved tests and comments
(5) numpy support

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • django/branches/gis/django/contrib/gis/geos/GEOSGeometry.py

    r5587 r5637  
    1 # Copyright (c) 2007, Justin Bronn 
    2 # All rights reserved. 
    3 # 
    4 # Redistribution and use in source and binary forms, with or without modification, 
    5 # are permitted provided that the following conditions are met: 
    6 # 
    7 #     1. Redistributions of source code must retain the above copyright notice,  
    8 #        this list of conditions and the following disclaimer. 
    9 #     
    10 #     2. Redistributions in binary form must reproduce the above copyright  
    11 #        notice, this list of conditions and the following disclaimer in the 
    12 #        documentation and/or other materials provided with the distribution. 
    13 # 
    14 #     3. Neither the name of GEOSGeometry nor the names of its contributors may be used 
    15 #        to endorse or promote products derived from this software without 
    16 #        specific prior written permission. 
    17 # 
    18 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 
    19 # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
    20 # WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
    21 # DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR 
    22 # ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 
    23 # (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 
    24 # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 
    25 # ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
    26 # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
    27 # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
    28 # 
    29  
    301# Trying not to pollute the namespace. 
    312from ctypes import \ 
    32      CDLL, CFUNCTYPE, byref, string_at, create_string_buffer, \ 
     3     byref, string_at, create_string_buffer, pointer, \ 
    334     c_char_p, c_double, c_float, c_int, c_uint, c_size_t 
    34 import os, sys 
    35 from types import StringType, IntType, FloatType 
    36  
    37 """ 
    38   The goal of this module is to be a ctypes wrapper around the GEOS library 
    39   that will work on both *NIX and Windows systems.  Specifically, this uses 
    40   the GEOS C api. 
    41  
    42   I have several motivations for doing this: 
    43     (1) The GEOS SWIG wrapper is no longer maintained, and requires the 
    44         installation of SWIG. 
    45     (2) The PCL implementation is over 2K+ lines of C and would make 
    46         PCL a requisite package for the GeoDjango application stack. 
    47     (3) Windows and Mac compatibility becomes substantially easier, and does not 
    48         require the additional compilation of PCL or GEOS and SWIG -- all that 
    49         is needed is a Win32 or Mac compiled GEOS C library (dll or dylib) 
    50         in a location that Python can read (e.g. 'C:\Python25'). 
    51  
    52   In summary, I wanted to wrap GEOS in a more maintainable and portable way using 
    53    only Python and the excellent ctypes library (now standard in Python 2.5). 
    54  
    55   In the spirit of loose coupling, this library does not require Django or 
    56    GeoDjango.  Only the GEOS C library and ctypes are needed for the platform 
    57    of your choice. 
    58  
    59   For more information about GEOS: 
    60     http://geos.refractions.net 
    61    
    62   For more info about PCL and the discontinuation of the Python GEOS 
    63    library see Sean Gillies' writeup (and subsequent update) at: 
    64      http://zcologia.com/news/150/geometries-for-python/ 
    65      http://zcologia.com/news/429/geometries-for-python-update/ 
    66 """ 
    67  
    68 # Setting the appropriate name for the GEOS-C library, depending on which 
    69 # OS and POSIX platform we're running. 
    70 if os.name == 'nt': 
    71     # Windows NT library 
    72     lib_name = 'libgeos_c-1.dll' 
    73 elif os.name == 'posix': 
    74     platform = os.uname()[0] # Using os.uname() 
    75     if platform in ('Linux', 'SunOS'): 
    76         # Linux or Solaris shared library 
    77         lib_name = 'libgeos_c.so' 
    78     elif platform == 'Darwin': 
    79         # Mac OSX Shared Library (Thanks Matt!) 
    80         lib_name = 'libgeos_c.dylib' 
    81     else: 
    82         raise GEOSException, 'Unknown POSIX platform "%s"' % platform 
    83 else: 
    84     raise GEOSException, 'Unsupported OS "%s"' % os.name 
    85  
    86 # The GEOSException class 
    87 class GEOSException(Exception): pass 
    88 class GEOSGeometryIndexError(GEOSException, KeyError): 
    89     """This exception is raised when an invalid index is encountered, and has 
    90     the 'silent_variable_feature' attribute set to true.  This ensures that 
    91     django's templates proceed to use the next lookup type gracefully when 
    92     an Exception is raised.  Fixes ticket #4740. 
    93     """ 
    94     silent_variable_failure = True 
    95  
    96 # Getting the GEOS C library.  The C interface (CDLL) is used for 
    97 #  both *NIX and Windows. 
    98 # See the GEOS C API source code for more details on the library function calls: 
    99 #  http://geos.refractions.net/ro/doxygen_docs/html/geos__c_8h-source.html 
    100 lgeos = CDLL(lib_name) 
    101  
    102 # The notice and error handler C function callback definitions. 
    103 #  Supposed to mimic the GEOS message handler (C below): 
    104 #  "typedef void (*GEOSMessageHandler)(const char *fmt, ...);" 
    105 NOTICEFUNC = CFUNCTYPE(None, c_char_p, c_char_p) 
    106 def notice_h(fmt, list): 
    107     sys.stdout.write((fmt + '\n') % list) 
    108 notice_h = NOTICEFUNC(notice_h) 
    109  
    110 ERRORFUNC = CFUNCTYPE(None, c_char_p, c_char_p) 
    111 def error_h(fmt, list): 
    112     if list: 
    113         err_msg = fmt % list 
    114     else: 
    115         err_msg = fmt 
    116     sys.stderr.write(err_msg) 
    117 error_h = ERRORFUNC(error_h) 
    118  
    119 # The initGEOS routine should be called first, however, that routine takes 
    120 #  the notice and error functions as parameters.  Here is the C code that 
    121 #  we need to wrap: 
    122 #  "extern void GEOS_DLL initGEOS(GEOSMessageHandler notice_function, GEOSMessageHandler error_function);" 
    123 lgeos.initGEOS(notice_h, error_h) 
     5from types import StringType, IntType, FloatType, TupleType, ListType 
     6 
     7# Getting GEOS-related dependencies. 
     8from django.contrib.gis.geos.libgeos import lgeos, GEOSPointer, HAS_NUMPY 
     9from django.contrib.gis.geos.GEOSError import GEOSException 
     10from django.contrib.gis.geos.GEOSCoordSeq import GEOSCoordSeq, create_cs 
     11 
     12if HAS_NUMPY: 
     13    from numpy import ndarray, array 
    12414 
    12515class GEOSGeometry(object): 
    12616    "A class that, generally, encapsulates a GEOS geometry." 
    12717     
    128     _g = 0  # Initially NULL 
    129  
    13018    #### Python 'magic' routines #### 
    131     def __init__(self, input, input_type='wkt'): 
    132         "Takes an input and the type of the input for initialization." 
    133  
    134         if isinstance(input, StringType): 
     19    def __init__(self, geo_input, input_type='wkt', child=False): 
     20        """The constructor for GEOS geometry objects.  May take the following 
     21        strings as inputs, WKT ("wkt"), HEXEWKB ("hex", PostGIS-specific canonical form). 
     22 
     23        When a hex string is to be used, the `input_type` keyword should be set with 'hex'. 
     24 
     25        The `child` keyword is for internal use only, and indicates to the garbage collector 
     26          not to delete this geometry if it was spawned from a parent (e.g., the exterior 
     27          ring from a polygon). 
     28        """ 
     29 
     30        # Initially, setting the pointer to NULL 
     31        self._ptr = GEOSPointer(0) 
     32 
     33        if isinstance(geo_input, StringType): 
    13534            if input_type == 'wkt': 
    13635                # If the geometry is in WKT form 
    137                 g = lgeos.GEOSGeomFromWKT(c_char_p(input)) 
     36                g = lgeos.GEOSGeomFromWKT(c_char_p(geo_input)) 
    13837            elif input_type == 'hex': 
    13938                # If the geometry is in HEX form. 
    140                 sz = c_size_t(len(input)) 
    141                 buf = create_string_buffer(input) 
     39                sz = c_size_t(len(geo_input)) 
     40                buf = create_string_buffer(geo_input) 
    14241                g = lgeos.GEOSGeomFromHEX_buf(buf, sz) 
    14342            else: 
    144                 raise GEOSException, 'GEOS input geometry type "%s" not supported!' % input_type 
    145         elif isinstance(input, IntType): 
    146             # When the input is a C pointer (Python integer) 
    147             g = input 
     43                raise TypeError, 'GEOS input geometry type "%s" not supported.' % input_type 
     44        elif isinstance(geo_input, (IntType, GEOSPointer)): 
     45            # When the input is either a raw pointer value (an integer), or a GEOSPointer object. 
     46            g = geo_input 
    14847        else: 
    14948            # Invalid geometry type. 
    150             raise GEOSException, 'Improper geometry input type: %s' % str(type(input)) 
    151  
    152         # If the geometry pointer is NULL (0), raise an exception. 
    153         if not g: 
    154             raise GEOSException, 'Invalid %s given!' % input_type.upper() 
    155         else: 
    156             self._g = g 
    157  
    158         # Setting the class type (e.g. 'Point', 'Polygon', etc.) 
     49            raise TypeError, 'Improper geometry input type: %s' % str(type(geo_input)) 
     50 
     51        if bool(g): 
     52            # If we have a GEOSPointer object, just set the '_ptr' attribute with g 
     53            if isinstance(g, GEOSPointer): self._ptr = g 
     54            else: self._ptr.set(g) # Otherwise, set the address 
     55        else: 
     56            raise GEOSException, 'Could not initialize GEOS Geometry with given input.' 
     57 
     58        # Setting the 'child' flag -- when the object is labeled with this flag 
     59        #  it will not be destroyed by __del__().  This is used for child geometries from 
     60        #  parent geometries (e.g., LinearRings from a Polygon, Points from a MultiPoint, etc.). 
     61        self._child = child 
     62 
     63        # Setting the class type (e.g., 'Point', 'Polygon', etc.) 
    15964        self.__class__ = GEOS_CLASSES[self.geom_type] 
    16065 
     66        # Extra setup needed for Geometries that may be parents. 
     67        if isinstance(self, GeometryCollection): self._geoms = {} 
     68        if isinstance(self, Polygon): self._rings = {} 
     69 
    16170    def __del__(self): 
    162         "This cleans up the memory allocated for the geometry." 
    163         if self._g: lgeos.GEOSGeom_destroy(self._g) 
     71        "Destroys this geometry -- only if the pointer is valid and this is not a child geometry." 
     72        #print 'Deleting %s (child=%s, valid=%s)' % (self.geom_type, self._child, self._ptr.valid) 
     73        if self._ptr.valid and not self._child: lgeos.GEOSGeom_destroy(self._ptr()) 
    16474 
    16575    def __str__(self): 
    …… 
    17181        return self.equals(other) 
    17282 
     83    #### Coordinate Sequence Routines #### 
     84    def _cache_cs(self): 
     85        "Caches the coordinate sequence for this Geometry." 
     86        if not hasattr(self, '_cs'): 
     87            # Only these geometries are allowed to have coordinate sequences. 
     88            if self.geom_type in ('LineString', 'LinearRing', 'Point'): 
     89                self._cs = GEOSCoordSeq(GEOSPointer(lgeos.GEOSGeom_getCoordSeq(self._ptr())), self.hasz) 
     90            else: 
     91                self._cs = None 
     92 
     93    @property 
     94    def coord_seq(self): 
     95        "Returns the coordinate sequence for the geometry." 
     96        # Getting the coordinate sequence for the geometry 
     97        self._cache_cs() 
     98 
     99        # Returning a GEOSCoordSeq wrapped around the pointer. 
     100        return self._cs 
     101 
    173102    #### Geometry Info #### 
    174103    @property 
    175104    def geom_type(self): 
    176105        "Returns a string representing the geometry type, e.g. 'Polygon'" 
    177         return string_at(lgeos.GEOSGeomType(self._g)) 
     106        return string_at(lgeos.GEOSGeomType(self._ptr())) 
    178107 
    179108    @property 
    180109    def geom_typeid(self): 
    181110        "Returns an integer representing the geometry type." 
    182         return lgeos.GEOSGeomTypeId(self._g) 
     111        return lgeos.GEOSGeomTypeId(self._ptr()) 
    183112 
    184113    @property 
    185114    def num_geom(self): 
    186115        "Returns the number of geometries in the geometry." 
    187         n = lgeos.GEOSGetNumGeometries(self._g) 
    188         if n == -1: raise GEOSException, 'Error getting number of geometries!' 
     116        n = lgeos.GEOSGetNumGeometries(self._ptr()) 
     117        if n == -1: raise GEOSException, 'Error getting number of geometries.' 
    189118        else: return n 
    190119 
    …… 
    192121    def num_coords(self): 
    193122        "Returns the number of coordinates in the geometry." 
    194         n = lgeos.GEOSGetNumCoordinates(self._g) 
    195         if n == -1: raise GEOSException, 'Error getting number of coordinates!' 
     123        n = lgeos.GEOSGetNumCoordinates(self._ptr()) 
     124        if n == -1: raise GEOSException, 'Error getting number of coordinates.' 
    196125        else: return n 
    197126 
    …… 
    204133    def dims(self): 
    205134        "Returns the dimension of this Geometry (0=point, 1=line, 2=surface)." 
    206         return lgeos.GEOSGeom_getDimensions(self._g) 
    207  
    208     @property 
    209     def coord_seq(self): 
    210         "Returns the coordinate sequence for the geometry." 
    211  
    212         # Only these geometries can return coordinate sequences 
    213         if self.geom_type not in ['LineString', 'LinearRing', 'Point']: 
    214             return None 
    215  
    216         # Getting the coordinate sequence for the geometry 
    217         cs = lgeos.GEOSGeom_getCoordSeq(self._g) 
    218  
    219         # Cloning the coordinate sequence (if the original is returned, 
    220         #  and it is garbage-collected we will get a segmentation fault!) 
    221         clone = lgeos.GEOSCoordSeq_clone(cs) 
    222         return GEOSCoordSeq(clone, z=self.hasz) 
     135        return lgeos.GEOSGeom_getDimensions(self._ptr()) 
    223136 
    224137    def normalize(self): 
    225138        "Converts this Geometry to normal form (or canonical form)." 
    226         status = lgeos.GEOSNormalize(self._g) 
     139        status = lgeos.GEOSNormalize(self._ptr()) 
    227140        if status == -1: raise GEOSException, 'failed to normalize geometry' 
    228141 
    229     def _predicate(self, val): 
    230         "Checks the result, 2 for exception, 1 on true, 0 on false." 
    231         if val == 0: 
    232             return False 
    233         elif val == 1: 
    234             return True 
    235         else: 
    236             raise GEOSException, 'Predicate exception occurred!' 
    237  
    238     ### Unary predicates ### 
     142    def _unary_predicate(self, func): 
     143        "Returns the result, or raises an exception for the given unary predicate function." 
     144        val = func(self._ptr()) 
     145        if val == 0: return False 
     146        elif val == 1: return True 
     147        else: raise GEOSException, '%s: exception occurred.' % func.__name__ 
     148 
     149    def _binary_predicate(self, func, other, *args): 
     150        "Returns the result, or raises an exception for the given binary predicate function." 
     151        if not isinstance(other, GEOSGeometry): 
     152            raise TypeError, 'Binary predicate operation ("%s") requires another GEOSGeometry instance.' % func.__name__ 
     153        val = func(self._ptr(), other._ptr(), *args) 
     154        if val == 0: return False 
     155        elif val == 1: return True 
     156        else: raise GEOSException, '%s: exception occurred.' % func.__name__ 
     157 
     158    #### Unary predicates #### 
    239159    @property 
    240160    def empty(self): 
    241161        "Returns a boolean indicating whether the set of points in this geometry are empty." 
    242         return self._predicate(lgeos.GEOSisEmpty(self._g)) 
     162        return self._unary_predicate(lgeos.GEOSisEmpty) 
    243163 
    244164    @property 
    245165    def valid(self): 
    246166        "This property tests the validity of this geometry." 
    247         return self._predicate(lgeos.GEOSisValid(self._g)) 
     167        return self._unary_predicate(lgeos.GEOSisValid) 
    248168 
    249169    @property 
    250170    def simple(self): 
    251171        "Returns false if the Geometry not simple." 
    252         return self._predicate(lgeos.GEOSisSimple(self._g)) 
     172        return self._unary_predicate(lgeos.GEOSisSimple) 
    253173 
    254174    @property 
    255175    def ring(self): 
    256176        "Returns whether or not the geometry is a ring." 
    257         return self._predicate(lgeos.GEOSisRing(self._g)) 
     177        return self._unary_predicate(lgeos.GEOSisRing) 
    258178 
    259179    @property 
    260180    def hasz(self): 
    261181        "Returns whether the geometry has a 3D dimension." 
    262         return self._predicate(lgeos.GEOSHasZ(self._g)) 
     182        return self._unary_predicate(lgeos.GEOSHasZ) 
    263183 
    264184    #### Binary predicates. #### 
    …… 
    268188        if len(pattern) > 9: 
    269189            raise GEOSException, 'invalid intersection matrix pattern' 
    270         return self._predicate(lgeos.GEOSRelatePattern(self._g, other._g, c_char_p(pattern))) 
     190        return self._binary_predicate(lgeos.GEOSRelatePattern, other, c_char_p(pattern)) 
    271191 
    272192    def disjoint(self, other): 
    273193        "Returns true if the DE-9IM intersection matrix for the two Geometrys is FF*FF****." 
    274         return self._predicate(lgeos.GEOSDisjoint(self._g, other._g)) 
     194        return self._binary_predicate(lgeos.GEOSDisjoint, other) 
    275195 
    276196    def touches(self, other): 
    277197        "Returns true if the DE-9IM intersection matrix for the two Geometrys is FT*******, F**T***** or F***T****." 
    278         return self._predicate(lgeos.GEOSTouches(self._g, other._g)) 
     198        return self._binary_predicate(lgeos.GEOSTouches, other) 
    279199 
    280200    def intersects(self, other): 
    281201        "Returns true if disjoint returns false." 
    282         return self._predicate(lgeos.GEOSIntersects(self._g, other._g)) 
     202        return self._binary_predicate(lgeos.GEOSIntersects, other) 
    283203 
    284204    def crosses(self, other): 
    285205        """Returns true if the DE-9IM intersection matrix for the two Geometrys is T*T****** (for a point and a curve, 
    286206        a point and an area or a line and an area) 0******** (for two curves).""" 
    287         return self._predicate(lgeos.GEOSCrosses(self._g, other._g)) 
     207        return self._binary_predicate(lgeos.GEOSCrosses, other) 
    288208 
    289209    def within(self, other): 
    290210        "Returns true if the DE-9IM intersection matrix for the two Geometrys is T*F**F***." 
    291         return self._predicate(lgeos.GEOSWithin(self._g, other._g)) 
     211        return self._binary_predicate(lgeos.GEOSWithin, other) 
    292212 
    293213    def contains(self, other): 
    294214        "Returns true if other.within(this) returns true." 
    295         return self._predicate(lgeos.GEOSContains(self._g, other._g)) 
     215        return self._binary_predicate(lgeos.GEOSContains, other) 
    296216 
    297217    def overlaps(self, other): 
    298218        """Returns true if the DE-9IM intersection matrix for the two Geometrys is T*T***T** (for two points 
    299219        or two surfaces) 1*T***T** (for two curves).""" 
    300         return self._predicate(lgeos.GEOSOverlaps(self._g, other._g)) 
     220        return self._binary_predicate(lgeos.GEOSOverlaps, other) 
    301221 
    302222    def equals(self, other): 
    303223        "Returns true if the DE-9IM intersection matrix for the two Geometrys is T*F**FFF*." 
    304         return self._predicate(lgeos.GEOSEquals(self._g, other._g)) 
     224        return self._binary_predicate(lgeos.GEOSEquals, other) 
    305225 
    306226    def equals_exact(self, other, tolerance=0): 
    307227        "Returns true if the two Geometrys are exactly equal, up to a specified tolerance." 
    308228        tol = c_double(tolerance) 
    309         return self._predicate(lgeos.GEOSEqualsExact(self._g, other._g, tol)) 
     229        return self._binary_predicate(lgeos.GEOSEqualsExact, other, tol) 
    310230 
    311231    #### SRID Routines #### 
    …… 
    313233    def srid(self): 
    314234        "Gets the SRID for the geometry, returns None if no SRID is set." 
    315         s = lgeos.GEOSGetSRID(self._g) 
     235        s = lgeos.GEOSGetSRID(self._ptr()) 
    316236        if s == 0: 
    317237            return None 
    …… 
    321241    def set_srid(self, srid): 
    322242        "Sets the SRID for the geometry." 
    323         lgeos.GEOSSetSRID(self._g, c_int(srid)) 
     243        lgeos.GEOSSetSRID(self._ptr(), c_int(srid)) 
    324244     
    325245    #### Output Routines #### 
    …… 
    327247    def wkt(self): 
    328248        "Returns the WKT of the Geometry." 
    329         return string_at(lgeos.GEOSGeomToWKT(self._g)) 
     249        return string_at(lgeos.GEOSGeomToWKT(self._ptr())) 
    330250 
    331251    @property 
    …… 
    333253        "Returns the WKBHEX of the Geometry." 
    334254        sz = c_size_t() 
    335         h = lgeos.GEOSGeomToHEX_buf(self._g, byref(sz)) 
     255        h = lgeos.GEOSGeomToHEX_buf(self._ptr(), byref(sz)) 
    336256        return string_at(h, sz.value) 
    337257 
    338258    #### Topology Routines #### 
     259    def _unary_topology(self, func, *args): 
     260        "Returns a GEOSGeometry for the given unary (only takes one geomtry as a paramter) topological operation." 
     261        return GEOSGeometry(func(self._ptr(), *args)) 
     262 
     263    def _binary_topology(self, func, other, *args): 
     264        "Returns a GEOSGeometry for the given binary (takes two geometries as parameters) topological operation." 
     265        return GEOSGeometry(func(self._ptr(), other._ptr(), *args)) 
     266 
    339267    def buffer(self, width, quadsegs=8): 
    340268        """Returns a geometry that represents all points whose distance from this 
    …… 
    344272        (Text from PostGIS documentation at ch. 6.1.3) 
    345273        """ 
    346         if not isinstance(width, FloatType): 
     274        if not isinstance(width, (FloatType, IntType)): 
    347275            raise TypeError, 'width parameter must be a float' 
    348276        if not isinstance(quadsegs, IntType): 
    349277            raise TypeError, 'quadsegs parameter must be an integer' 
    350         b = lgeos.GEOSBuffer(self._g, c_double(width), c_int(quadsegs)) 
    351         return GEOSGeometry(b) 
     278        return self._unary_topology(lgeos.GEOSBuffer, c_double(width), c_int(quadsegs)) 
    352279 
    353280    @property 
    354281    def envelope(self): 
    355282        "Return the envelope for this geometry (a polygon)." 
    356         return GEOSGeometry(lgeos.GEOSEnvelope(self._g)) 
     283        return self._unary_topology(lgeos.GEOSEnvelope) 
    357284 
    358285    @property 
    …… 
    361288        of highest dimension (since the lower-dimension geometries contribute zero 
    362289        "weight" to the centroid).""" 
    363         return GEOSGeometry(lgeos.GEOSGetCentroid(self._g)) 
     290        return self._unary_topology(lgeos.GEOSGetCentroid) 
    364291 
    365292    @property 
    366293    def boundary(self): 
    367294        "Returns the boundary as a newly allocated Geometry object." 
    368         return GEOSGeometry(lgeos.GEOSBoundary(self._g)) 
     295        return self._unary_topology(lgeos.GEOSBoundary) 
    369296 
    370297    @property 
    371298    def convex_hull(self): 
    372299        "Returns the smallest convex Polygon that contains all the points in the Geometry." 
    373         return GEOSGeometry(lgeos.GEOSConvexHull(self._g)) 
     300        return self._unary_topology(lgeos.GEOSConvexHull) 
    374301 
    375302    @property 
    376303    def point_on_surface(self): 
    377304        "Computes an interior point of this Geometry." 
    378         return GEOSGeometry(lgeos.GEOSPointOnSurface(self._g)) 
     305        return self._unary_topology(lgeos.GEOSPointOnSurface) 
    379306 
    380307    def relate(self, other): 
    381308        "Returns the DE-9IM intersection matrix for this geometry and the other." 
    382         return string_at(lgeos.GEOSRelate(self._g, other._g)) 
     309        return string_at(lgeos.GEOSRelate(self._ptr(), other._ptr())) 
    383310 
    384311    def difference(self, other): 
    385312        """Returns a Geometry representing the points making up this Geometry 
    386313        that do not make up other.""" 
    387         return GEOSGeometry(lgeos.GEOSDifference(self._g, other._g)) 
     314        return self._binary_topology(lgeos.GEOSDifference, other) 
    388315 
    389316    def sym_difference(self, other): 
    390317        """Returns a set combining the points in this Geometry not in other, 
    391318        and the points in other not in this Geometry.""" 
    392         return GEOSGeometry(lgeos.GEOSSymDifference(self._g, other._g)) 
     319        return self._binary_topology(lgeos.GEOSSymDifference, other) 
    393320 
    394321    def intersection(self, other): 
    395322        "Returns a Geometry representing the points shared by this Geometry and other." 
    396         return GEOSGeometry(lgeos.GEOSIntersection(self._g, other._g)) 
     323        return self._binary_topology(lgeos.GEOSIntersection, other) 
    397324 
    398325    def union(self, other): 
    399326        "Returns a Geometry representing all the points in this Geometry and other." 
    400         return GEOSGeometry(lgeos.GEOSUnion(self._g, other._g)) 
     327        return self._binary_topology(lgeos.GEOSUnion, other) 
    401328 
    402329    #### Other Routines #### 
    …… 
    405332        "Returns the area of the Geometry." 
    406333        a = c_double() 
    407         status = lgeos.GEOSArea(self._g, byref(a)) 
    408         if not status: 
    409             return None 
    410         else: 
    411             return a.value 
    412  
    413 class GEOSCoordSeq(object): 
    414     "The internal representation of a list of coordinates inside a Geometry." 
    415  
    416     def __init__(self, ptr, z=False): 
    417         "Initializes from a GEOS pointer." 
    418         self._cs = ptr 
    419         self._z = z 
    420  
    421     def __del__(self): 
    422         "Destroys the reference to this coordinate sequence." 
    423         if self._cs: lgeos.GEOSCoordSeq_destroy(self._cs) 
    424  
    425     def __iter__(self): 
    426         "Iterates over each point in the coordinate sequence." 
    427         for i in xrange(self.size): 
    428             yield self.__getitem__(i) 
    429  
    430     def __len__(self): 
    431         "Returns the number of points in the coordinate sequence." 
    432         return int(self.size) 
    433  
    434     def __str__(self): 
    435         "The string representation of the coordinate sequence." 
    436         return str(self.tuple) 
    437  
    438     def _checkindex(self, index): 
    439         "Checks the index." 
    440         sz = self.size 
    441         if (sz < 1) or (index < 0) or (index >= sz): 
    442             raise GEOSGeometryIndexError, 'invalid GEOS Geometry index: %s' % str(index) 
    443  
    444     def _checkdim(self, dim): 
    445         "Checks the given dimension." 
    446         if dim < 0 or dim > 2: 
    447             raise GEOSException, 'invalid ordinate dimension "%d"' % dim 
     334        status = lgeos.GEOSArea(self._ptr(), byref(a)) 
     335        if not status: return None 
     336        else: return a.value 
     337 
     338    def clone(self): 
     339        "Clones this Geometry." 
     340        return GEOSGeometry(lgeos.GEOSGeom_clone(self._ptr())) 
     341     
     342class Point(GEOSGeometry): 
     343 
     344    def __init__(self, x, y=None, z=None): 
     345        """The Point object may be initialized with either a tuple, or individual 
     346        parameters.  For example: 
     347          >>> p = Point((5, 23)) # 2D point, passed in as a tuple 
     348          >>> p = Point(5, 23, 8) # 3D point, passed in with individual parameters 
     349        """ 
    448350         
    449     def __getitem__(self, index): 
    450         "Can use the index [] operator to get coordinate sequence at an index." 
    451         coords = [self.getX(index), self.getY(index)] 
    452         if self.dims == 3 and self._z: 
    453             coords.append(self.getZ(index)) 
    454         return tuple(coords) 
    455  
    456     def __setitem__(self, index, value): 
    457         "Can use the index [] operator to set coordinate sequence at an index." 
    458         if self.dims == 3 and self._z: 
    459             n_args = 3 
    460             set_3d = True 
    461         else: 
    462             n_args = 2 
    463             set_3d = False 
    464         if len(value) != n_args: 
    465             raise GEOSException, 'Improper value given!' 
    466         self.setX(index, value[0]) 
    467         self.setY(index, value[1]) 
    468         if set_3d: self.setZ(index, value[2]) 
    469          
    470     # Getting and setting the X coordinate for the given index. 
    471     def getX(self, index): 
    472         return self.getOrdinate(0, index) 
    473  
    474     def setX(self, index, value): 
    475         self.setOrdinate(0, index, value) 
    476  
    477     # Getting and setting the Y coordinate for the given index. 
    478     def getY(self, index): 
    479         return self.getOrdinate(1, index) 
    480  
    481     def setY(self, index, value): 
    482         self.setOrdinate(1, index, value) 
    483  
    484     # Getting and setting the Z coordinate for the given index 
    485     def getZ(self, index): 
    486         return self.getOrdinate(2, index) 
    487  
    488     def setZ(self, index, value): 
    489         self.setOrdinate(2, index, value) 
    490  
    491     def getOrdinate(self, dimension, index): 
    492         "Gets the value for the given dimension and index." 
    493         self._checkindex(index) 
    494         self._checkdim(dimension) 
    495  
    496         # Wrapping the dimension, index 
    497         dim = c_uint(dimension) 
    498         idx = c_uint(index) 
    499  
    500         # 'd' is the value of the point, passed in by reference 
    501         d = c_double() 
    502         status = lgeos.GEOSCoordSeq_getOrdinate(self._cs, idx, dim, byref(d)) 
    503         if status == 0: 
    504             raise GEOSException, 'Could not get the ordinate for (dim, index): (%d, %d)' % (dimension, index) 
    505         return d.value 
    506  
    507     def setOrdinate(self, dimension, index, value): 
    508         "Sets the value for the given dimension and index." 
    509         self._checkindex(index) 
    510         self._checkdim(dimension) 
    511  
    512         # Wrapping the dimension, index 
    513         dim = c_uint(dimension) 
    514         idx = c_uint(index) 
    515  
    516         # Setting the ordinate 
    517         status = lgeos.GEOSCoordSeq_setOrdinate(self._cs, idx, dim, c_double(value)) 
    518         if status == 0: 
    519             raise GEOSException, 'Could not set the ordinate for (dim, index): (%d, %d)' % (dimension, index) 
    520  
    521     ### Dimensions ### 
    522     @property 
    523     def size(self): 
    524         "Returns the size of this coordinate sequence." 
    525         n = c_uint(0) 
    526         status = lgeos.GEOSCoordSeq_getSize(self._cs, byref(n)) 
    527         if status == 0: 
    528             raise GEOSException, 'Could not get CoordSeq size!' 
    529         return n.value 
    530  
    531     @property 
    532     def dims(self): 
    533         "Returns the dimensions of this coordinate sequence." 
    534         n = c_uint(0) 
    535         status = lgeos.GEOSCoordSeq_getDimensions(self._cs, byref(n)) 
    536         if status == 0: 
    537             raise GEOSException, 'Could not get CoordSeq dimensoins!' 
    538         return n.value 
    539  
    540     @property 
    541     def hasz(self): 
    542         "Inherits this from the parent geometry." 
    543         return self._z 
    544  
    545     ### Other Methods ### 
    546     @property 
    547     def tuple(self): 
    548         "Returns a tuple version of this coordinate sequence." 
    549         n = self.size 
    550         if n == 1: 
    551             return self.__getitem__(0) 
    552         else: 
    553             return tuple(self.__getitem__(i) for i in xrange(n)) 
    554              
    555 # Factory coordinate sequence Function 
    556 def createCoordSeq(size, dims): 
    557         return GEOSCoordSeq(lgeos.GEOSCoordSeq_create(c_uint(size), c_uint(dims))) 
    558  
    559 class Point(GEOSGeometry): 
    560  
    561     def _cache_cs(self): 
    562         "Caches the coordinate sequence." 
    563         if not hasattr(self, '_cs'): self._cs = self.coord_seq         
     351        if isinstance(x, (TupleType, ListType)): 
     352            # Here a tuple or list was passed in under the ``x`` parameter. 
     353            ndim = len(x) 
     354            if ndim < 2 or ndim > 3: 
     355                raise TypeError, 'Invalid sequence parameter: %s' % str(x) 
     356            coords = x 
     357        elif isinstance(x, (IntType, FloatType)) and isinstance(y, (IntType, FloatType)): 
     358            # Here X, Y, and (optionally) Z were passed in individually as parameters. 
     359            if isinstance(z, (IntType, FloatType)): 
     360                ndim = 3 
     361                coords = [x, y, z] 
     362            else: 
     363                ndim = 2 
     364                coords = [x, y] 
     365        else: 
     366            raise TypeError, 'Invalid parameters given for Point initialization.' 
     367 
     368        # Creating the coordinate sequence 
     369        cs = create_cs(c_uint(1), c_uint(ndim)) 
     370 
     371        # Setting the X 
     372        status = lgeos.GEOSCoordSeq_setX(cs, c_uint(0), c_double(coords[0])) 
     373        if not status: raise GEOSException, 'Could not set X during Point initialization.' 
     374 
     375        # Setting the Y 
     376        status = lgeos.GEOSCoordSeq_setY(cs, c_uint(0), c_double(coords[1])) 
     377        if not status: raise GEOSException, 'Could not set Y during Point initialization.' 
     378 
     379        # Setting the Z 
     380        if ndim == 3: 
     381            status = lgeos.GEOSCoordSeq_setZ(cs, c_uint(0), c_double(coords[2])) 
     382 
     383        # Initializing from the geometry, and getting a Python object 
     384        super(Point, self).__init__(lgeos.GEOSGeom_createPoint(cs)) 
    564385 
    565386    def _getOrdinate(self, dim, idx): 
    …… 
    568389        return self._cs.getOrdinate(dim, idx) 
    569390 
    570     @property 
    571     def x(self): 
     391    def _setOrdinate(self, dim, idx, value): 
     392        "The coordinate sequence setOrdinate() wrapper." 
     393        self._cache_cs() 
     394        self._cs.setOrdinate(dim, idx, value) 
     395 
     396    def get_x(self): 
    572397        "Returns the X component of the Point." 
    573398        return self._getOrdinate(0, 0) 
    574399 
    575     @property 
    576     def y(self): 
     400    def set_x(self, value): 
     401        "Sets the X component of the Point." 
     402        self._setOrdinate(0, 0, value) 
     403 
     404    def get_y(self): 
    577405        "Returns the Y component of the Point." 
    578406        return self._getOrdinate(1, 0) 
    579407 
    580     @property 
    581     def z(self): 
     408    def set_y(self, value): 
     409        "Sets the Y component of the Point." 
     410        self._setOrdinate(1, 0, value) 
     411 
     412    def get_z(self): 
    582413        "Returns the Z component of the Point." 
    583414        if self.hasz: 
    …… 
    586417            return None 
    587418 
     419    def set_z(self, value): 
     420        "Sets the Z component of the Point." 
     421        if self.hasz: 
     422            self._setOrdinate(2, 0, value) 
     423        else: 
     424            raise GEOSException, 'Cannot set Z on 2D Point.' 
     425     
     426    # X, Y, Z properties 
     427    x = property(get_x, set_x) 
     428    y = property(get_y, set_y) 
     429    z = property(get_z, set_z) 
     430 
    588431    @property 
    589432    def tuple(self): 
    …… 
    594437class LineString(GEOSGeometry): 
    595438 
    596     def _cache_cs(self): 
    597         "Caches the coordinate sequence." 
    598         if not hasattr(self, '_cs'): self._cs = self.coord_seq  
     439    #### Python 'magic' routines #### 
     440    def __init__(self, coords, ring=False): 
     441        """Initializes on the given sequence, may take lists, tuples, or NumPy arrays 
     442        of X,Y pairs.""" 
     443 
     444        if isinstance(coords, (TupleType, ListType)): 
     445            ncoords = len(coords) 
     446            first = True 
     447            for coord in coords: 
     448                if not isinstance(coord, (TupleType, ListType)): 
     449                    raise TypeError, 'each coordinate should be a sequence (list or tuple)' 
     450                if first: 
     451                    ndim = len(coord) 
     452                    self._checkdim(ndim) 
     453                    first = False 
     454                else: 
     455                    if len(coord) != ndim: raise TypeError, 'Dimension mismatch.' 
     456            numpy_coords = False 
     457        elif HAS_NUMPY and isinstance(coords, ndarray): 
     458            shape = coords.shape 
     459            if len(shape) != 2: raise TypeError, 'Too many dimensions.' 
     460            self._checkdim(shape[1]) 
     461            ncoords = shape[0] 
     462            ndim = shape[1] 
     463            numpy_coords = True 
     464        else: 
     465            raise TypeError, 'Invalid initialization input for LineStrings.' 
     466 
     467        # Creating the coordinate sequence 
     468        cs = GEOSCoordSeq(GEOSPointer(create_cs(c_uint(ncoords), c_uint(ndim)))) 
     469 
     470        # Setting each point in the coordinate sequence 
     471        for i in xrange(ncoords): 
     472            if numpy_coords: cs[i] = coords[i,:] 
     473            else: cs[i] = coords[i]         
     474 
     475        # Getting the initialization function 
     476        if ring: 
     477            func = lgeos.GEOSGeom_createLinearRing 
     478        else: 
     479            func = lgeos.GEOSGeom_createLineString 
     480        
     481        # Calling the base geometry initialization with the returned pointer from the function. 
     482        super(LineString, self).__init__(func(cs._ptr())) 
    599483 
    600484    def __getitem__(self, index): 
    601485        "Gets the point at the specified index." 
    602486        self._cache_cs() 
    603         if index < 0 or index >= self._cs.size: 
    604             raise GEOSGeometryIndexError, 'invalid GEOS Geometry index: %s' % str(index) 
    605         else: 
    606             return self._cs[index] 
     487        return self._cs[index] 
     488 
     489    def __setitem__(self, index, value): 
     490        "Sets the point at the specified index, e.g., line_str[0] = (1, 2)." 
     491        self._cache_cs() 
     492        self._cs[index] = value 
    607493 
    608494    def __iter__(self): 
    609495        "Allows iteration over this LineString." 
    610         self._cache_cs() 
    611         for i in xrange(self._cs.size): 
     496        for i in xrange(self.__len__()): 
    612497            yield self.__getitem__(index) 
    613498 
    …… 
    617502        return len(self._cs) 
    618503 
     504    def _checkdim(self, dim): 
     505        if dim not in (2, 3): raise TypeError, 'Dimension mismatch.' 
     506 
     507    #### Sequence Properties #### 
    619508    @property 
    620509    def tuple(self): 
    …… 
    623512        return self._cs.tuple 
    624513 
     514    def _listarr(self, func): 
     515        """Internal routine that returns a sequence (list) corresponding with 
     516        the given function.  Will return a numpy array if possible.""" 
     517        lst = [func(i) for i in xrange(self.__len__())] # constructing the list, using the function 
     518        if HAS_NUMPY: return array(lst) # ARRRR! 
     519        else: return lst 
     520 
     521    @property 
     522    def array(self): 
     523        "Returns a numpy array for the LineString." 
     524        self._cache_cs() 
     525        return self._listarr(self._cs.__getitem__) 
     526 
     527    @property 
     528    def x(self): 
     529        "Returns a list or numpy array of the X variable." 
     530        self._cache_cs() 
     531        return self._listarr(self._cs.getX) 
     532     
     533    @property 
     534    def y(self): 
     535        "Returns a list or numpy array of the Y variable." 
     536        self._cache_cs() 
     537        return self._listarr(self._cs.getY) 
     538 
     539    @property 
     540    def z(self): 
     541        "Returns a list or numpy array of the Z variable." 
     542        self._cache_cs() 
     543        if not self.hasz: return None 
     544        else: return self._listarr(self._cs.getZ) 
     545 
    625546# LinearRings are LineStrings used within Polygons. 
    626 class LinearRing(LineString): pass 
     547class LinearRing(LineString): 
     548    def __init__(self, coords): 
     549        "Overriding the initialization function to set the ring keyword." 
     550        super(LinearRing, self).__init__(coords, ring=True) 
    627551 
    628552class Polygon(GEOSGeometry): 
    629553 
     554    def __del__(self): 
     555        "Override the GEOSGeometry delete routine to safely take care of any spawned rings." 
     556        # Nullifying the pointers to internal rings, preventing any attempted future access 
     557        for k in self._rings: self._rings[k].nullify() 
     558        super(Polygon, self).__del__() # Calling the parent __del__() method. 
     559     
    630560    def __getitem__(self, index): 
    631561        """Returns the ring at the specified index.  The first index, 0, will always 
    …… 
    650580 
    651581    def get_interior_ring(self, ring_i): 
    652         "Gets the interior ring at the specified index." 
     582        """Gets the interior ring at the specified index, 
     583        0 is for the first interior ring, not the exterior ring.""" 
    653584 
    654585        # Making sure the ring index is within range 
    655         if ring_i >= self.num_interior_rings: 
    656             raise GEOSException, 'Invalid ring index.' 
    657          
    658         # Getting a clone of the ring geometry at the given ring index. 
    659         r = lgeos.GEOSGeom_clone(lgeos.GEOSGetInteriorRingN(self._g, c_int(ring_i))) 
    660         return GEOSGeometry(r, 'geos') 
     586        if ring_i < 0 or ring_i >= self.num_interior_rings: 
     587            raise IndexError, 'ring index out of range' 
     588 
     589        # Placing the ring in internal rings dictionary. 
     590&n