diff --git a/django/contrib/gis/geos/GEOSCoordSeq.py b/django/contrib/gis/geos/GEOSCoordSeq.py
new file mode 100644
index 0000000000..faacccc471
--- /dev/null
+++ b/django/contrib/gis/geos/GEOSCoordSeq.py
@@ -0,0 +1,162 @@
+from django.contrib.gis.geos.libgeos import lgeos
+from django.contrib.gis.geos.GEOSError import GEOSException, GEOSGeometryIndexError
+from ctypes import c_double, c_int, c_uint, byref
+
+"""
+  This module houses the GEOSCoordSeq object, and is used internally
+  by GEOSGeometry to house the actual coordinates of the Point,
+  LineString, and LinearRing geometries.
+"""
+
+class GEOSCoordSeq(object):
+    "The internal representation of a list of coordinates inside a Geometry."
+
+    #### Python 'magic' routines ####
+    def __init__(self, ptr, z=False):
+        "Initializes from a GEOS pointer."
+        self._ptr = ptr
+        self._z = z
+
+    def __iter__(self):
+        "Iterates over each point in the coordinate sequence."
+        for i in xrange(self.size):
+            yield self.__getitem__(i)
+
+    def __len__(self):
+        "Returns the number of points in the coordinate sequence."
+        return int(self.size)
+
+    def __str__(self):
+        "The string representation of the coordinate sequence."
+        return str(self.tuple)
+
+    def __getitem__(self, index):
+        "Can use the index [] operator to get coordinate sequence at an index."
+        coords = [self.getX(index), self.getY(index)]
+        if self.dims == 3 and self._z:
+            coords.append(self.getZ(index))
+        return tuple(coords)
+
+    def __setitem__(self, index, value):
+        "Can use the index [] operator to set coordinate sequence at an index."
+        if self.dims == 3 and self._z:
+            n_args = 3
+            set_3d = True
+        else:
+            n_args = 2
+            set_3d = False
+        if len(value) != n_args:
+            raise TypeError, 'Dimension of value does not match.'
+        self.setX(index, value[0])
+        self.setY(index, value[1])
+        if set_3d: self.setZ(index, value[2])
+
+    #### Internal Routines ####
+    def _checkindex(self, index):
+        "Checks the index."
+        sz = self.size
+        if (sz < 1) or (index < 0) or (index >= sz):
+            raise GEOSGeometryIndexError, 'invalid GEOS Geometry index: %s' % str(index)
+
+    def _checkdim(self, dim):
+        "Checks the given dimension."
+        if dim < 0 or dim > 2:
+            raise GEOSException, 'invalid ordinate dimension "%d"' % dim
+
+    #### Ordinate getting and setting routines ####
+    def getOrdinate(self, dimension, index):
+        "Gets the value for the given dimension and index."
+        self._checkindex(index)
+        self._checkdim(dimension)
+
+        # Wrapping the dimension, index
+        dim = c_uint(dimension)
+        idx = c_uint(index)
+
+        # 'd' is the value of the point, passed in by reference
+        d = c_double()
+        status = lgeos.GEOSCoordSeq_getOrdinate(self._ptr(), idx, dim, byref(d))
+        if status == 0:
+            raise GEOSException, 'could not retrieve %th ordinate at index: %s' % (str(dimension), str(index))
+        return d.value
+
+    def setOrdinate(self, dimension, index, value):
+        "Sets the value for the given dimension and index."
+        self._checkindex(index)
+        self._checkdim(dimension)
+
+        # Wrapping the dimension, index
+        dim = c_uint(dimension)
+        idx = c_uint(index)
+
+        # Setting the ordinate
+        status = lgeos.GEOSCoordSeq_setOrdinate(self._ptr(), idx, dim, c_double(value))
+        if status == 0:
+            raise GEOSException, 'Could not set the ordinate for (dim, index): (%d, %d)' % (dimension, index)
+
+    def getX(self, index):
+        "Get the X value at the index."
+        return self.getOrdinate(0, index)
+
+    def setX(self, index, value):
+        "Set X with the value at the given index."
+        self.setOrdinate(0, index, value)
+
+    def getY(self, index):
+        "Get the Y value at the given index."
+        return self.getOrdinate(1, index)
+
+    def setY(self, index, value):
+        "Set Y with the value at the given index."
+        self.setOrdinate(1, index, value)
+
+    def getZ(self, index):
+        "Get Z with the value at the given index."
+        return self.getOrdinate(2, index)
+
+    def setZ(self, index, value):
+        "Set Z with the value at the given index."
+        self.setOrdinate(2, index, value)
+
+    ### Dimensions ###
+    @property
+    def size(self):
+        "Returns the size of this coordinate sequence."
+        n = c_uint(0)
+        status = lgeos.GEOSCoordSeq_getSize(self._ptr(), byref(n))
+        if status == 0:
+            raise GEOSException, 'Could not get CoordSeq size.'
+        return n.value
+
+    @property
+    def dims(self):
+        "Returns the dimensions of this coordinate sequence."
+        n = c_uint(0)
+        status = lgeos.GEOSCoordSeq_getDimensions(self._ptr(), byref(n))
+        if status == 0:
+            raise GEOSException, 'Could not get CoordSeq dimensions.'
+        return n.value
+
+    @property
+    def hasz(self):
+        "Inherits this from the parent geometry."
+        return self._z
+
+    ### Other Methods ###
+    @property
+    def clone(self):
+        "Clones this coordinate sequence."
+        pass
+
+    @property
+    def tuple(self):
+        "Returns a tuple version of this coordinate sequence."
+        n = self.size
+        if n == 1:
+            return self.__getitem__(0)
+        else:
+            return tuple(self.__getitem__(i) for i in xrange(n))
+
+# ctypes prototype for the Coordinate Sequence creation factory
+create_cs = lgeos.GEOSCoordSeq_create
+create_cs.argtypes = [c_uint, c_uint]
diff --git a/django/contrib/gis/geos/GEOSError.py b/django/contrib/gis/geos/GEOSError.py
new file mode 100644
index 0000000000..219a71929c
--- /dev/null
+++ b/django/contrib/gis/geos/GEOSError.py
@@ -0,0 +1,15 @@
+
+class GEOSException(Exception):
+    "The base GEOS exception, indicates a GEOS-related error."
+    pass
+
+class GEOSGeometryIndexError(GEOSException, KeyError):
+    """This exception is raised when an invalid index is encountered, and has
+    the 'silent_variable_feature' attribute set to true.  This ensures that
+    django's templates proceed to use the next lookup type gracefully when
+    an Exception is raised.  Fixes ticket #4740.
+    """
+    # "If, during the method lookup, a method raises an exception, the exception
+    #  will be propagated, unless the exception has an attribute silent_variable_failure
+    #  whose value is True." -- django template docs.
+    silent_variable_failure = True
diff --git a/django/contrib/gis/geos/GEOSGeometry.py b/django/contrib/gis/geos/GEOSGeometry.py
index d3ad0b5e2a..8ee2b566b6 100644
--- a/django/contrib/gis/geos/GEOSGeometry.py
+++ b/django/contrib/gis/geos/GEOSGeometry.py
@@ -1,166 +1,76 @@
-# Copyright (c) 2007, Justin Bronn
-# All rights reserved.
-#
-# Redistribution and use in source and binary forms, with or without modification,
-# are permitted provided that the following conditions are met:
-#
-#     1. Redistributions of source code must retain the above copyright notice, 
-#        this list of conditions and the following disclaimer.
-#    
-#     2. Redistributions in binary form must reproduce the above copyright 
-#        notice, this list of conditions and the following disclaimer in the
-#        documentation and/or other materials provided with the distribution.
-#
-#     3. Neither the name of GEOSGeometry nor the names of its contributors may be used
-#        to endorse or promote products derived from this software without
-#        specific prior written permission.
-#
-# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
-# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
-# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
-# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
-# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
-# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
-# ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
-# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-#
-
 # Trying not to pollute the namespace.
 from ctypes import \
-     CDLL, CFUNCTYPE, byref, string_at, create_string_buffer, \
+     byref, string_at, create_string_buffer, pointer, \
      c_char_p, c_double, c_float, c_int, c_uint, c_size_t
-import os, sys
-from types import StringType, IntType, FloatType
+from types import StringType, IntType, FloatType, TupleType, ListType
 
-"""
-  The goal of this module is to be a ctypes wrapper around the GEOS library
-  that will work on both *NIX and Windows systems.  Specifically, this uses
-  the GEOS C api.
+# Getting GEOS-related dependencies.
+from django.contrib.gis.geos.libgeos import lgeos, GEOSPointer, HAS_NUMPY
+from django.contrib.gis.geos.GEOSError import GEOSException
+from django.contrib.gis.geos.GEOSCoordSeq import GEOSCoordSeq, create_cs
 
-  I have several motivations for doing this:
-    (1) The GEOS SWIG wrapper is no longer maintained, and requires the
-        installation of SWIG.
-    (2) The PCL implementation is over 2K+ lines of C and would make
-        PCL a requisite package for the GeoDjango application stack.
-    (3) Windows and Mac compatibility becomes substantially easier, and does not
-        require the additional compilation of PCL or GEOS and SWIG -- all that
-        is needed is a Win32 or Mac compiled GEOS C library (dll or dylib)
-        in a location that Python can read (e.g. 'C:\Python25').
-
-  In summary, I wanted to wrap GEOS in a more maintainable and portable way using
-   only Python and the excellent ctypes library (now standard in Python 2.5).
-
-  In the spirit of loose coupling, this library does not require Django or
-   GeoDjango.  Only the GEOS C library and ctypes are needed for the platform
-   of your choice.
-
-  For more information about GEOS:
-    http://geos.refractions.net
-  
-  For more info about PCL and the discontinuation of the Python GEOS
-   library see Sean Gillies' writeup (and subsequent update) at:
-     http://zcologia.com/news/150/geometries-for-python/
-     http://zcologia.com/news/429/geometries-for-python-update/
-"""
-
-# Setting the appropriate name for the GEOS-C library, depending on which
-# OS and POSIX platform we're running.
-if os.name == 'nt':
-    # Windows NT library
-    lib_name = 'libgeos_c-1.dll'
-elif os.name == 'posix':
-    platform = os.uname()[0] # Using os.uname()
-    if platform in ('Linux', 'SunOS'):
-        # Linux or Solaris shared library
-        lib_name = 'libgeos_c.so'
-    elif platform == 'Darwin':
-        # Mac OSX Shared Library (Thanks Matt!)
-        lib_name = 'libgeos_c.dylib'
-    else:
-        raise GEOSException, 'Unknown POSIX platform "%s"' % platform
-else:
-    raise GEOSException, 'Unsupported OS "%s"' % os.name
-
-# The GEOSException class
-class GEOSException(Exception): pass
-class GEOSGeometryIndexError(GEOSException, KeyError):
-    """This exception is raised when an invalid index is encountered, and has
-    the 'silent_variable_feature' attribute set to true.  This ensures that
-    django's templates proceed to use the next lookup type gracefully when
-    an Exception is raised.  Fixes ticket #4740.
-    """
-    silent_variable_failure = True
-
-# Getting the GEOS C library.  The C interface (CDLL) is used for
-#  both *NIX and Windows.
-# See the GEOS C API source code for more details on the library function calls:
-#  http://geos.refractions.net/ro/doxygen_docs/html/geos__c_8h-source.html
-lgeos = CDLL(lib_name)
-
-# The notice and error handler C function callback definitions.
-#  Supposed to mimic the GEOS message handler (C below):
-#  "typedef void (*GEOSMessageHandler)(const char *fmt, ...);"
-NOTICEFUNC = CFUNCTYPE(None, c_char_p, c_char_p)
-def notice_h(fmt, list):
-    sys.stdout.write((fmt + '\n') % list)
-notice_h = NOTICEFUNC(notice_h)
-
-ERRORFUNC = CFUNCTYPE(None, c_char_p, c_char_p)
-def error_h(fmt, list):
-    if list:
-        err_msg = fmt % list
-    else:
-        err_msg = fmt
-    sys.stderr.write(err_msg)
-error_h = ERRORFUNC(error_h)
-
-# The initGEOS routine should be called first, however, that routine takes
-#  the notice and error functions as parameters.  Here is the C code that
-#  we need to wrap:
-#  "extern void GEOS_DLL initGEOS(GEOSMessageHandler notice_function, GEOSMessageHandler error_function);"
-lgeos.initGEOS(notice_h, error_h)
+if HAS_NUMPY:
+    from numpy import ndarray, array
 
 class GEOSGeometry(object):
     "A class that, generally, encapsulates a GEOS geometry."
     
-    _g = 0  # Initially NULL
-
     #### Python 'magic' routines ####
-    def __init__(self, input, input_type='wkt'):
-        "Takes an input and the type of the input for initialization."
+    def __init__(self, geo_input, input_type='wkt', child=False):
+        """The constructor for GEOS geometry objects.  May take the following
+        strings as inputs, WKT ("wkt"), HEXEWKB ("hex", PostGIS-specific canonical form).
 
-        if isinstance(input, StringType):
+        When a hex string is to be used, the `input_type` keyword should be set with 'hex'.
+
+        The `child` keyword is for internal use only, and indicates to the garbage collector
+          not to delete this geometry if it was spawned from a parent (e.g., the exterior
+          ring from a polygon).
+        """
+
+        # Initially, setting the pointer to NULL
+        self._ptr = GEOSPointer(0)
+
+        if isinstance(geo_input, StringType):
             if input_type == 'wkt':
                 # If the geometry is in WKT form
-                g = lgeos.GEOSGeomFromWKT(c_char_p(input))
+                g = lgeos.GEOSGeomFromWKT(c_char_p(geo_input))
             elif input_type == 'hex':
                 # If the geometry is in HEX form.
-                sz = c_size_t(len(input))
-                buf = create_string_buffer(input)
+                sz = c_size_t(len(geo_input))
+                buf = create_string_buffer(geo_input)
                 g = lgeos.GEOSGeomFromHEX_buf(buf, sz)
             else:
-                raise GEOSException, 'GEOS input geometry type "%s" not supported!' % input_type
-        elif isinstance(input, IntType):
-            # When the input is a C pointer (Python integer)
-            g = input
+                raise TypeError, 'GEOS input geometry type "%s" not supported.' % input_type
+        elif isinstance(geo_input, (IntType, GEOSPointer)):
+            # When the input is either a raw pointer value (an integer), or a GEOSPointer object.
+            g = geo_input
         else:
             # Invalid geometry type.
-            raise GEOSException, 'Improper geometry input type: %s' % str(type(input))
+            raise TypeError, 'Improper geometry input type: %s' % str(type(geo_input))
 
-        # If the geometry pointer is NULL (0), raise an exception.
-        if not g:
-            raise GEOSException, 'Invalid %s given!' % input_type.upper()
+        if bool(g):
+            # If we have a GEOSPointer object, just set the '_ptr' attribute with g
+            if isinstance(g, GEOSPointer): self._ptr = g
+            else: self._ptr.set(g) # Otherwise, set the address
         else:
-            self._g = g
+            raise GEOSException, 'Could not initialize GEOS Geometry with given input.'
 
-        # Setting the class type (e.g. 'Point', 'Polygon', etc.)
+        # Setting the 'child' flag -- when the object is labeled with this flag
+        #  it will not be destroyed by __del__().  This is used for child geometries from
+        #  parent geometries (e.g., LinearRings from a Polygon, Points from a MultiPoint, etc.).
+        self._child = child
+
+        # Setting the class type (e.g., 'Point', 'Polygon', etc.)
         self.__class__ = GEOS_CLASSES[self.geom_type]
 
+        # Extra setup needed for Geometries that may be parents.
+        if isinstance(self, GeometryCollection): self._geoms = {}
+        if isinstance(self, Polygon): self._rings = {}
+
     def __del__(self):
-        "This cleans up the memory allocated for the geometry."
-        if self._g: lgeos.GEOSGeom_destroy(self._g)
+        "Destroys this geometry -- only if the pointer is valid and this is not a child geometry."
+        #print 'Deleting %s (child=%s, valid=%s)' % (self.geom_type, self._child, self._ptr.valid)
+        if self._ptr.valid and not self._child: lgeos.GEOSGeom_destroy(self._ptr())
 
     def __str__(self):
         "WKT is used for the string representation."
@@ -170,29 +80,48 @@ class GEOSGeometry(object):
         "Equivalence testing."
         return self.equals(other)
 
+    #### Coordinate Sequence Routines ####
+    def _cache_cs(self):
+        "Caches the coordinate sequence for this Geometry."
+        if not hasattr(self, '_cs'):
+            # Only these geometries are allowed to have coordinate sequences.
+            if self.geom_type in ('LineString', 'LinearRing', 'Point'):
+                self._cs = GEOSCoordSeq(GEOSPointer(lgeos.GEOSGeom_getCoordSeq(self._ptr())), self.hasz)
+            else:
+                self._cs = None
+
+    @property
+    def coord_seq(self):
+        "Returns the coordinate sequence for the geometry."
+        # Getting the coordinate sequence for the geometry
+        self._cache_cs()
+
+        # Returning a GEOSCoordSeq wrapped around the pointer.
+        return self._cs
+
     #### Geometry Info ####
     @property
     def geom_type(self):
         "Returns a string representing the geometry type, e.g. 'Polygon'"
-        return string_at(lgeos.GEOSGeomType(self._g))
+        return string_at(lgeos.GEOSGeomType(self._ptr()))
 
     @property
     def geom_typeid(self):
         "Returns an integer representing the geometry type."
-        return lgeos.GEOSGeomTypeId(self._g)
+        return lgeos.GEOSGeomTypeId(self._ptr())
 
     @property
     def num_geom(self):
         "Returns the number of geometries in the geometry."
-        n = lgeos.GEOSGetNumGeometries(self._g)
-        if n == -1: raise GEOSException, 'Error getting number of geometries!'
+        n = lgeos.GEOSGetNumGeometries(self._ptr())
+        if n == -1: raise GEOSException, 'Error getting number of geometries.'
         else: return n
 
     @property
     def num_coords(self):
         "Returns the number of coordinates in the geometry."
-        n = lgeos.GEOSGetNumCoordinates(self._g)
-        if n == -1: raise GEOSException, 'Error getting number of coordinates!'
+        n = lgeos.GEOSGetNumCoordinates(self._ptr())
+        if n == -1: raise GEOSException, 'Error getting number of coordinates.'
         else: return n
 
     @property
@@ -203,63 +132,54 @@ class GEOSGeometry(object):
     @property
     def dims(self):
         "Returns the dimension of this Geometry (0=point, 1=line, 2=surface)."
-        return lgeos.GEOSGeom_getDimensions(self._g)
-
-    @property
-    def coord_seq(self):
-        "Returns the coordinate sequence for the geometry."
-
-        # Only these geometries can return coordinate sequences
-        if self.geom_type not in ['LineString', 'LinearRing', 'Point']:
-            return None
-
-        # Getting the coordinate sequence for the geometry
-        cs = lgeos.GEOSGeom_getCoordSeq(self._g)
-
-        # Cloning the coordinate sequence (if the original is returned,
-        #  and it is garbage-collected we will get a segmentation fault!)
-        clone = lgeos.GEOSCoordSeq_clone(cs)
-        return GEOSCoordSeq(clone, z=self.hasz)
+        return lgeos.GEOSGeom_getDimensions(self._ptr())
 
     def normalize(self):
         "Converts this Geometry to normal form (or canonical form)."
-        status = lgeos.GEOSNormalize(self._g)
+        status = lgeos.GEOSNormalize(self._ptr())
         if status == -1: raise GEOSException, 'failed to normalize geometry'
 
-    def _predicate(self, val):
-        "Checks the result, 2 for exception, 1 on true, 0 on false."
-        if val == 0:
-            return False
-        elif val == 1:
-            return True
-        else:
-            raise GEOSException, 'Predicate exception occurred!'
+    def _unary_predicate(self, func):
+        "Returns the result, or raises an exception for the given unary predicate function."
+        val = func(self._ptr())
+        if val == 0: return False
+        elif val == 1: return True
+        else: raise GEOSException, '%s: exception occurred.' % func.__name__
 
-    ### Unary predicates ###
+    def _binary_predicate(self, func, other, *args):
+        "Returns the result, or raises an exception for the given binary predicate function."
+        if not isinstance(other, GEOSGeometry):
+            raise TypeError, 'Binary predicate operation ("%s") requires another GEOSGeometry instance.' % func.__name__
+        val = func(self._ptr(), other._ptr(), *args)
+        if val == 0: return False
+        elif val == 1: return True
+        else: raise GEOSException, '%s: exception occurred.' % func.__name__
+
+    #### Unary predicates ####
     @property
     def empty(self):
         "Returns a boolean indicating whether the set of points in this geometry are empty."
-        return self._predicate(lgeos.GEOSisEmpty(self._g))
+        return self._unary_predicate(lgeos.GEOSisEmpty)
 
     @property
     def valid(self):
         "This property tests the validity of this geometry."
-        return self._predicate(lgeos.GEOSisValid(self._g))
+        return self._unary_predicate(lgeos.GEOSisValid)
 
     @property
     def simple(self):
         "Returns false if the Geometry not simple."
-        return self._predicate(lgeos.GEOSisSimple(self._g))
+        return self._unary_predicate(lgeos.GEOSisSimple)
 
     @property
     def ring(self):
         "Returns whether or not the geometry is a ring."
-        return self._predicate(lgeos.GEOSisRing(self._g))
+        return self._unary_predicate(lgeos.GEOSisRing)
 
     @property
     def hasz(self):
         "Returns whether the geometry has a 3D dimension."
-        return self._predicate(lgeos.GEOSHasZ(self._g))
+        return self._unary_predicate(lgeos.GEOSHasZ)
 
     #### Binary predicates. ####
     def relate_pattern(self, other, pattern):
@@ -267,52 +187,52 @@ class GEOSGeometry(object):
         the two Geometrys match the elements in pattern."""
         if len(pattern) > 9:
             raise GEOSException, 'invalid intersection matrix pattern'
-        return self._predicate(lgeos.GEOSRelatePattern(self._g, other._g, c_char_p(pattern)))
+        return self._binary_predicate(lgeos.GEOSRelatePattern, other, c_char_p(pattern))
 
     def disjoint(self, other):
         "Returns true if the DE-9IM intersection matrix for the two Geometrys is FF*FF****."
-        return self._predicate(lgeos.GEOSDisjoint(self._g, other._g))
+        return self._binary_predicate(lgeos.GEOSDisjoint, other)
 
     def touches(self, other):
         "Returns true if the DE-9IM intersection matrix for the two Geometrys is FT*******, F**T***** or F***T****."
-        return self._predicate(lgeos.GEOSTouches(self._g, other._g))
+        return self._binary_predicate(lgeos.GEOSTouches, other)
 
     def intersects(self, other):
         "Returns true if disjoint returns false."
-        return self._predicate(lgeos.GEOSIntersects(self._g, other._g))
+        return self._binary_predicate(lgeos.GEOSIntersects, other)
 
     def crosses(self, other):
         """Returns true if the DE-9IM intersection matrix for the two Geometrys is T*T****** (for a point and a curve,
         a point and an area or a line and an area) 0******** (for two curves)."""
-        return self._predicate(lgeos.GEOSCrosses(self._g, other._g))
+        return self._binary_predicate(lgeos.GEOSCrosses, other)
 
     def within(self, other):
         "Returns true if the DE-9IM intersection matrix for the two Geometrys is T*F**F***."
-        return self._predicate(lgeos.GEOSWithin(self._g, other._g))
+        return self._binary_predicate(lgeos.GEOSWithin, other)
 
     def contains(self, other):
         "Returns true if other.within(this) returns true."
-        return self._predicate(lgeos.GEOSContains(self._g, other._g))
+        return self._binary_predicate(lgeos.GEOSContains, other)
 
     def overlaps(self, other):
         """Returns true if the DE-9IM intersection matrix for the two Geometrys is T*T***T** (for two points
         or two surfaces) 1*T***T** (for two curves)."""
-        return self._predicate(lgeos.GEOSOverlaps(self._g, other._g))
+        return self._binary_predicate(lgeos.GEOSOverlaps, other)
 
     def equals(self, other):
         "Returns true if the DE-9IM intersection matrix for the two Geometrys is T*F**FFF*."
-        return self._predicate(lgeos.GEOSEquals(self._g, other._g))
+        return self._binary_predicate(lgeos.GEOSEquals, other)
 
     def equals_exact(self, other, tolerance=0):
         "Returns true if the two Geometrys are exactly equal, up to a specified tolerance."
         tol = c_double(tolerance)
-        return self._predicate(lgeos.GEOSEqualsExact(self._g, other._g, tol))
+        return self._binary_predicate(lgeos.GEOSEqualsExact, other, tol)
 
     #### SRID Routines ####
     @property
     def srid(self):
         "Gets the SRID for the geometry, returns None if no SRID is set."
-        s = lgeos.GEOSGetSRID(self._g)
+        s = lgeos.GEOSGetSRID(self._ptr())
         if s == 0:
             return None
         else:
@@ -320,22 +240,30 @@ class GEOSGeometry(object):
 
     def set_srid(self, srid):
         "Sets the SRID for the geometry."
-        lgeos.GEOSSetSRID(self._g, c_int(srid))
+        lgeos.GEOSSetSRID(self._ptr(), c_int(srid))
     
     #### Output Routines ####
     @property
     def wkt(self):
         "Returns the WKT of the Geometry."
-        return string_at(lgeos.GEOSGeomToWKT(self._g))
+        return string_at(lgeos.GEOSGeomToWKT(self._ptr()))
 
     @property
     def hex(self):
         "Returns the WKBHEX of the Geometry."
         sz = c_size_t()
-        h = lgeos.GEOSGeomToHEX_buf(self._g, byref(sz))
+        h = lgeos.GEOSGeomToHEX_buf(self._ptr(), byref(sz))
         return string_at(h, sz.value)
 
     #### Topology Routines ####
+    def _unary_topology(self, func, *args):
+        "Returns a GEOSGeometry for the given unary (only takes one geomtry as a paramter) topological operation."
+        return GEOSGeometry(func(self._ptr(), *args))
+
+    def _binary_topology(self, func, other, *args):
+        "Returns a GEOSGeometry for the given binary (takes two geometries as parameters) topological operation."
+        return GEOSGeometry(func(self._ptr(), other._ptr(), *args))
+
     def buffer(self, width, quadsegs=8):
         """Returns a geometry that represents all points whose distance from this
         Geometry is less than or equal to distance. Calculations are in the
@@ -343,248 +271,163 @@ class GEOSGeometry(object):
         the number of segment used to approximate a quarter circle (defaults to 8).
         (Text from PostGIS documentation at ch. 6.1.3)
         """
-        if not isinstance(width, FloatType):
+        if not isinstance(width, (FloatType, IntType)):
             raise TypeError, 'width parameter must be a float'
         if not isinstance(quadsegs, IntType):
             raise TypeError, 'quadsegs parameter must be an integer'
-        b = lgeos.GEOSBuffer(self._g, c_double(width), c_int(quadsegs))
-        return GEOSGeometry(b)
+        return self._unary_topology(lgeos.GEOSBuffer, c_double(width), c_int(quadsegs))
 
     @property
     def envelope(self):
         "Return the envelope for this geometry (a polygon)."
-        return GEOSGeometry(lgeos.GEOSEnvelope(self._g))
+        return self._unary_topology(lgeos.GEOSEnvelope)
 
     @property
     def centroid(self):
         """The centroid is equal to the centroid of the set of component Geometrys
         of highest dimension (since the lower-dimension geometries contribute zero
         "weight" to the centroid)."""
-        return GEOSGeometry(lgeos.GEOSGetCentroid(self._g))
+        return self._unary_topology(lgeos.GEOSGetCentroid)
 
     @property
     def boundary(self):
         "Returns the boundary as a newly allocated Geometry object."
-        return GEOSGeometry(lgeos.GEOSBoundary(self._g))
+        return self._unary_topology(lgeos.GEOSBoundary)
 
     @property
     def convex_hull(self):
         "Returns the smallest convex Polygon that contains all the points in the Geometry."
-        return GEOSGeometry(lgeos.GEOSConvexHull(self._g))
+        return self._unary_topology(lgeos.GEOSConvexHull)
 
     @property
     def point_on_surface(self):
         "Computes an interior point of this Geometry."
-        return GEOSGeometry(lgeos.GEOSPointOnSurface(self._g))
+        return self._unary_topology(lgeos.GEOSPointOnSurface)
 
     def relate(self, other):
         "Returns the DE-9IM intersection matrix for this geometry and the other."
-        return string_at(lgeos.GEOSRelate(self._g, other._g))
+        return string_at(lgeos.GEOSRelate(self._ptr(), other._ptr()))
 
     def difference(self, other):
         """Returns a Geometry representing the points making up this Geometry
         that do not make up other."""
-        return GEOSGeometry(lgeos.GEOSDifference(self._g, other._g))
+        return self._binary_topology(lgeos.GEOSDifference, other)
 
     def sym_difference(self, other):
         """Returns a set combining the points in this Geometry not in other,
         and the points in other not in this Geometry."""
-        return GEOSGeometry(lgeos.GEOSSymDifference(self._g, other._g))
+        return self._binary_topology(lgeos.GEOSSymDifference, other)
 
     def intersection(self, other):
         "Returns a Geometry representing the points shared by this Geometry and other."
-        return GEOSGeometry(lgeos.GEOSIntersection(self._g, other._g))
+        return self._binary_topology(lgeos.GEOSIntersection, other)
 
     def union(self, other):
         "Returns a Geometry representing all the points in this Geometry and other."
-        return GEOSGeometry(lgeos.GEOSUnion(self._g, other._g))
+        return self._binary_topology(lgeos.GEOSUnion, other)
 
     #### Other Routines ####
     @property
     def area(self):
         "Returns the area of the Geometry."
         a = c_double()
-        status = lgeos.GEOSArea(self._g, byref(a))
-        if not status:
-            return None
-        else:
-            return a.value
-
-class GEOSCoordSeq(object):
-    "The internal representation of a list of coordinates inside a Geometry."
-
-    def __init__(self, ptr, z=False):
-        "Initializes from a GEOS pointer."
-        self._cs = ptr
-        self._z = z
-
-    def __del__(self):
-        "Destroys the reference to this coordinate sequence."
-        if self._cs: lgeos.GEOSCoordSeq_destroy(self._cs)
-
-    def __iter__(self):
-        "Iterates over each point in the coordinate sequence."
-        for i in xrange(self.size):
-            yield self.__getitem__(i)
-
-    def __len__(self):
-        "Returns the number of points in the coordinate sequence."
-        return int(self.size)
-
-    def __str__(self):
-        "The string representation of the coordinate sequence."
-        return str(self.tuple)
-
-    def _checkindex(self, index):
-        "Checks the index."
-        sz = self.size
-        if (sz < 1) or (index < 0) or (index >= sz):
-            raise GEOSGeometryIndexError, 'invalid GEOS Geometry index: %s' % str(index)
-
-    def _checkdim(self, dim):
-        "Checks the given dimension."
-        if dim < 0 or dim > 2:
-            raise GEOSException, 'invalid ordinate dimension "%d"' % dim
-        
-    def __getitem__(self, index):
-        "Can use the index [] operator to get coordinate sequence at an index."
-        coords = [self.getX(index), self.getY(index)]
-        if self.dims == 3 and self._z:
-            coords.append(self.getZ(index))
-        return tuple(coords)
-
-    def __setitem__(self, index, value):
-        "Can use the index [] operator to set coordinate sequence at an index."
-        if self.dims == 3 and self._z:
-            n_args = 3
-            set_3d = True
-        else:
-            n_args = 2
-            set_3d = False
-        if len(value) != n_args:
-            raise GEOSException, 'Improper value given!'
-        self.setX(index, value[0])
-        self.setY(index, value[1])
-        if set_3d: self.setZ(index, value[2])
-        
-    # Getting and setting the X coordinate for the given index.
-    def getX(self, index):
-        return self.getOrdinate(0, index)
-
-    def setX(self, index, value):
-        self.setOrdinate(0, index, value)
-
-    # Getting and setting the Y coordinate for the given index.
-    def getY(self, index):
-        return self.getOrdinate(1, index)
-
-    def setY(self, index, value):
-        self.setOrdinate(1, index, value)
-
-    # Getting and setting the Z coordinate for the given index
-    def getZ(self, index):
-        return self.getOrdinate(2, index)
-
-    def setZ(self, index, value):
-        self.setOrdinate(2, index, value)
-
-    def getOrdinate(self, dimension, index):
-        "Gets the value for the given dimension and index."
-        self._checkindex(index)
-        self._checkdim(dimension)
-
-        # Wrapping the dimension, index
-        dim = c_uint(dimension)
-        idx = c_uint(index)
-
-        # 'd' is the value of the point, passed in by reference
-        d = c_double()
-        status = lgeos.GEOSCoordSeq_getOrdinate(self._cs, idx, dim, byref(d))
-        if status == 0:
-            raise GEOSException, 'Could not get the ordinate for (dim, index): (%d, %d)' % (dimension, index)
-        return d.value
-
-    def setOrdinate(self, dimension, index, value):
-        "Sets the value for the given dimension and index."
-        self._checkindex(index)
-        self._checkdim(dimension)
-
-        # Wrapping the dimension, index
-        dim = c_uint(dimension)
-        idx = c_uint(index)
-
-        # Setting the ordinate
-        status = lgeos.GEOSCoordSeq_setOrdinate(self._cs, idx, dim, c_double(value))
-        if status == 0:
-            raise GEOSException, 'Could not set the ordinate for (dim, index): (%d, %d)' % (dimension, index)
-
-    ### Dimensions ###
-    @property
-    def size(self):
-        "Returns the size of this coordinate sequence."
-        n = c_uint(0)
-        status = lgeos.GEOSCoordSeq_getSize(self._cs, byref(n))
-        if status == 0:
-            raise GEOSException, 'Could not get CoordSeq size!'
-        return n.value
-
-    @property
-    def dims(self):
-        "Returns the dimensions of this coordinate sequence."
-        n = c_uint(0)
-        status = lgeos.GEOSCoordSeq_getDimensions(self._cs, byref(n))
-        if status == 0:
-            raise GEOSException, 'Could not get CoordSeq dimensoins!'
-        return n.value
-
-    @property
-    def hasz(self):
-        "Inherits this from the parent geometry."
-        return self._z
-
-    ### Other Methods ###
-    @property
-    def tuple(self):
-        "Returns a tuple version of this coordinate sequence."
-        n = self.size
-        if n == 1:
-            return self.__getitem__(0)
-        else:
-            return tuple(self.__getitem__(i) for i in xrange(n))
-            
-# Factory coordinate sequence Function
-def createCoordSeq(size, dims):
-        return GEOSCoordSeq(lgeos.GEOSCoordSeq_create(c_uint(size), c_uint(dims)))
+        status = lgeos.GEOSArea(self._ptr(), byref(a))
+        if not status: return None
+        else: return a.value
 
+    def clone(self):
+        "Clones this Geometry."
+        return GEOSGeometry(lgeos.GEOSGeom_clone(self._ptr()))
+    
 class Point(GEOSGeometry):
 
-    def _cache_cs(self):
-        "Caches the coordinate sequence."
-        if not hasattr(self, '_cs'): self._cs = self.coord_seq        
+    def __init__(self, x, y=None, z=None):
+        """The Point object may be initialized with either a tuple, or individual
+        parameters.  For example:
+          >>> p = Point((5, 23)) # 2D point, passed in as a tuple
+          >>> p = Point(5, 23, 8) # 3D point, passed in with individual parameters
+        """
+        
+        if isinstance(x, (TupleType, ListType)):
+            # Here a tuple or list was passed in under the ``x`` parameter.
+            ndim = len(x)
+            if ndim < 2 or ndim > 3:
+                raise TypeError, 'Invalid sequence parameter: %s' % str(x)
+            coords = x
+        elif isinstance(x, (IntType, FloatType)) and isinstance(y, (IntType, FloatType)):
+            # Here X, Y, and (optionally) Z were passed in individually as parameters.
+            if isinstance(z, (IntType, FloatType)):
+                ndim = 3
+                coords = [x, y, z]
+            else:
+                ndim = 2
+                coords = [x, y]
+        else:
+            raise TypeError, 'Invalid parameters given for Point initialization.'
+
+        # Creating the coordinate sequence
+        cs = create_cs(c_uint(1), c_uint(ndim))
+
+        # Setting the X
+        status = lgeos.GEOSCoordSeq_setX(cs, c_uint(0), c_double(coords[0]))
+        if not status: raise GEOSException, 'Could not set X during Point initialization.'
+
+        # Setting the Y
+        status = lgeos.GEOSCoordSeq_setY(cs, c_uint(0), c_double(coords[1]))
+        if not status: raise GEOSException, 'Could not set Y during Point initialization.'
+
+        # Setting the Z
+        if ndim == 3:
+            status = lgeos.GEOSCoordSeq_setZ(cs, c_uint(0), c_double(coords[2]))
+
+        # Initializing from the geometry, and getting a Python object
+        super(Point, self).__init__(lgeos.GEOSGeom_createPoint(cs))
 
     def _getOrdinate(self, dim, idx):
         "The coordinate sequence getOrdinate() wrapper."
         self._cache_cs()
         return self._cs.getOrdinate(dim, idx)
 
-    @property
-    def x(self):
+    def _setOrdinate(self, dim, idx, value):
+        "The coordinate sequence setOrdinate() wrapper."
+        self._cache_cs()
+        self._cs.setOrdinate(dim, idx, value)
+
+    def get_x(self):
         "Returns the X component of the Point."
         return self._getOrdinate(0, 0)
 
-    @property
-    def y(self):
+    def set_x(self, value):
+        "Sets the X component of the Point."
+        self._setOrdinate(0, 0, value)
+
+    def get_y(self):
         "Returns the Y component of the Point."
         return self._getOrdinate(1, 0)
 
-    @property
-    def z(self):
+    def set_y(self, value):
+        "Sets the Y component of the Point."
+        self._setOrdinate(1, 0, value)
+
+    def get_z(self):
         "Returns the Z component of the Point."
         if self.hasz:
             return self._getOrdinate(2, 0)
         else:
             return None
 
+    def set_z(self, value):
+        "Sets the Z component of the Point."
+        if self.hasz:
+            self._setOrdinate(2, 0, value)
+        else:
+            raise GEOSException, 'Cannot set Z on 2D Point.'
+    
+    # X, Y, Z properties
+    x = property(get_x, set_x)
+    y = property(get_y, set_y)
+    z = property(get_z, set_z)
+
     @property
     def tuple(self):
         "Returns a tuple of the point."
@@ -593,22 +436,64 @@ class Point(GEOSGeometry):
 
 class LineString(GEOSGeometry):
 
-    def _cache_cs(self):
-        "Caches the coordinate sequence."
-        if not hasattr(self, '_cs'): self._cs = self.coord_seq 
+    #### Python 'magic' routines ####
+    def __init__(self, coords, ring=False):
+        """Initializes on the given sequence, may take lists, tuples, or NumPy arrays
+        of X,Y pairs."""
+
+        if isinstance(coords, (TupleType, ListType)):
+            ncoords = len(coords)
+            first = True
+            for coord in coords:
+                if not isinstance(coord, (TupleType, ListType)):
+                    raise TypeError, 'each coordinate should be a sequence (list or tuple)'
+                if first:
+                    ndim = len(coord)
+                    self._checkdim(ndim)
+                    first = False
+                else:
+                    if len(coord) != ndim: raise TypeError, 'Dimension mismatch.'
+            numpy_coords = False
+        elif HAS_NUMPY and isinstance(coords, ndarray):
+            shape = coords.shape
+            if len(shape) != 2: raise TypeError, 'Too many dimensions.'
+            self._checkdim(shape[1])
+            ncoords = shape[0]
+            ndim = shape[1]
+            numpy_coords = True
+        else:
+            raise TypeError, 'Invalid initialization input for LineStrings.'
+
+        # Creating the coordinate sequence
+        cs = GEOSCoordSeq(GEOSPointer(create_cs(c_uint(ncoords), c_uint(ndim))))
+
+        # Setting each point in the coordinate sequence
+        for i in xrange(ncoords):
+            if numpy_coords: cs[i] = coords[i,:]
+            else: cs[i] = coords[i]        
+
+        # Getting the initialization function
+        if ring:
+            func = lgeos.GEOSGeom_createLinearRing
+        else:
+            func = lgeos.GEOSGeom_createLineString
+       
+        # Calling the base geometry initialization with the returned pointer from the function.
+        super(LineString, self).__init__(func(cs._ptr()))
 
     def __getitem__(self, index):
         "Gets the point at the specified index."
         self._cache_cs()
-        if index < 0 or index >= self._cs.size:
-            raise GEOSGeometryIndexError, 'invalid GEOS Geometry index: %s' % str(index)
-        else:
-            return self._cs[index]
+        return self._cs[index]
+
+    def __setitem__(self, index, value):
+        "Sets the point at the specified index, e.g., line_str[0] = (1, 2)."
+        self._cache_cs()
+        self._cs[index] = value
 
     def __iter__(self):
         "Allows iteration over this LineString."
-        self._cache_cs()
-        for i in xrange(self._cs.size):
+        for i in xrange(self.__len__()):
             yield self.__getitem__(index)
 
     def __len__(self):
@@ -616,17 +501,62 @@ class LineString(GEOSGeometry):
         self._cache_cs()
         return len(self._cs)
 
+    def _checkdim(self, dim):
+        if dim not in (2, 3): raise TypeError, 'Dimension mismatch.'
+
+    #### Sequence Properties ####
     @property
     def tuple(self):
         "Returns a tuple version of the geometry from the coordinate sequence."
         self._cache_cs()
         return self._cs.tuple
 
+    def _listarr(self, func):
+        """Internal routine that returns a sequence (list) corresponding with
+        the given function.  Will return a numpy array if possible."""
+        lst = [func(i) for i in xrange(self.__len__())] # constructing the list, using the function
+        if HAS_NUMPY: return array(lst) # ARRRR!
+        else: return lst
+
+    @property
+    def array(self):
+        "Returns a numpy array for the LineString."
+        self._cache_cs()
+        return self._listarr(self._cs.__getitem__)
+
+    @property
+    def x(self):
+        "Returns a list or numpy array of the X variable."
+        self._cache_cs()
+        return self._listarr(self._cs.getX)
+    
+    @property
+    def y(self):
+        "Returns a list or numpy array of the Y variable."
+        self._cache_cs()
+        return self._listarr(self._cs.getY)
+
+    @property
+    def z(self):
+        "Returns a list or numpy array of the Z variable."
+        self._cache_cs()
+        if not self.hasz: return None
+        else: return self._listarr(self._cs.getZ)
+
 # LinearRings are LineStrings used within Polygons.
-class LinearRing(LineString): pass
+class LinearRing(LineString):
+    def __init__(self, coords):
+        "Overriding the initialization function to set the ring keyword."
+        super(LinearRing, self).__init__(coords, ring=True)
 
 class Polygon(GEOSGeometry):
 
+    def __del__(self):
+        "Override the GEOSGeometry delete routine to safely take care of any spawned rings."
+        # Nullifying the pointers to internal rings, preventing any attempted future access
+        for k in self._rings: self._rings[k].nullify()
+        super(Polygon, self).__del__() # Calling the parent __del__() method.
+    
     def __getitem__(self, index):
         """Returns the ring at the specified index.  The first index, 0, will always
         return the exterior ring.  Indices > 0 will return the interior ring."""
@@ -649,15 +579,20 @@ class Polygon(GEOSGeometry):
         return self.num_interior_rings + 1
 
     def get_interior_ring(self, ring_i):
-        "Gets the interior ring at the specified index."
+        """Gets the interior ring at the specified index,
+        0 is for the first interior ring, not the exterior ring."""
 
         # Making sure the ring index is within range
-        if ring_i >= self.num_interior_rings:
-            raise GEOSException, 'Invalid ring index.'
-        
-        # Getting a clone of the ring geometry at the given ring index.
-        r = lgeos.GEOSGeom_clone(lgeos.GEOSGetInteriorRingN(self._g, c_int(ring_i)))
-        return GEOSGeometry(r, 'geos')
+        if ring_i < 0 or ring_i >= self.num_interior_rings:
+            raise IndexError, 'ring index out of range'
+
+        # Placing the ring in internal rings dictionary.
+        idx = ring_i+1 # the index for the polygon is +1 because of the exterior ring
+        if not idx in self._rings:
+            self._rings[idx] = GEOSPointer(lgeos.GEOSGetInteriorRingN(self._ptr(), c_int(ring_i)))
+
+        # Returning the ring at the given index.
+        return GEOSGeometry(self._rings[idx], child=True)
                                                         
     #### Polygon Properties ####
     @property
@@ -665,19 +600,18 @@ class Polygon(GEOSGeometry):
         "Returns the number of interior rings."
 
         # Getting the number of rings
-        n = lgeos.GEOSGetNumInteriorRings(self._g)
+        n = lgeos.GEOSGetNumInteriorRings(self._ptr())
 
         # -1 indicates an exception occurred
-        if n == -1: raise GEOSException, 'Error getting the number of interior rings!'
+        if n == -1: raise GEOSException, 'Error getting the number of interior rings.'
         else: return n
 
     @property
     def exterior_ring(self):
         "Gets the exterior ring of the Polygon."
-
-        # Getting a clone of the ring geometry
-        r = lgeos.GEOSGeom_clone(lgeos.GEOSGetExteriorRing(self._g))
-        return GEOSGeometry(r, 'geos')
+        # Returns exterior ring 
+        self._rings[0] = GEOSPointer(lgeos.GEOSGetExteriorRing((self._ptr())))
+        return GEOSGeometry(self._rings[0], child=True)
 
     @property
     def shell(self):
@@ -691,26 +625,37 @@ class Polygon(GEOSGeometry):
 
 class GeometryCollection(GEOSGeometry):
 
-    def _checkindex(self, index):
-        "Checks the given geometry index."
-        if index < 0 or index >= self.num_geom:
-            raise GEOSGeometryIndexError, 'invalid GEOS Geometry index: %s' % str(index)
-
-    def __iter__(self):
-        "For iteration on the multiple geometries."
-        for i in xrange(self.num_geom):
-            yield self.__getitem__(i)
+    def __del__(self):
+        "Override the GEOSGeometry delete routine to safely take care of any spawned geometries."
+        # Nullifying the pointers to internal geometries, preventing any attempted future access
+        for k in self._geoms: self._geoms[k].nullify()
+        super(GeometryCollection, self).__del__() # Calling the parent __del__() method.
 
     def __getitem__(self, index):
         "For indexing on the multiple geometries."
         self._checkindex(index)
+
+        # Setting an entry in the _geoms dictionary for the requested geometry.
+        if not index in self._geoms:
+            self._geoms[index] = GEOSPointer(lgeos.GEOSGetGeometryN(self._ptr(), c_int(index)))
+
         # Cloning the GEOS Geometry first, before returning it.
-        return GEOSGeometry(lgeos.GEOSGeom_clone(lgeos.GEOSGetGeometryN(self._g, c_int(index))))
+        return GEOSGeometry(self._geoms[index], child=True)
+
+    def __iter__(self):
+        "For iteration on the multiple geometries."
+        for i in xrange(self.__len__()):
+            yield self.__getitem__(i)
 
     def __len__(self):
         "Returns the number of geometries in this collection."
         return self.num_geom
 
+    def _checkindex(self, index):
+        "Checks the given geometry index."
+        if index < 0 or index >= self.num_geom:
+            raise GEOSGeometryIndexError, 'invalid GEOS Geometry index: %s' % str(index)
+
 # MultiPoint, MultiLineString, and MultiPolygon class definitions.
 class MultiPoint(GeometryCollection): pass
 class MultiLineString(GeometryCollection): pass
diff --git a/django/contrib/gis/geos/LICENSE b/django/contrib/gis/geos/LICENSE
new file mode 100644
index 0000000000..84cf485d00
--- /dev/null
+++ b/django/contrib/gis/geos/LICENSE
@@ -0,0 +1,27 @@
+Copyright (c) 2007, Justin Bronn
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+
+    1. Redistributions of source code must retain the above copyright notice, 
+       this list of conditions and the following disclaimer.
+   
+    2. Redistributions in binary form must reproduce the above copyright 
+       notice, this list of conditions and the following disclaimer in the
+       documentation and/or other materials provided with the distribution.
+
+    3. Neither the name of GEOSGeometry nor the names of its contributors may be used
+       to endorse or promote products derived from this software without
+       specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/django/contrib/gis/geos/__init__.py b/django/contrib/gis/geos/__init__.py
index 3ced1ad602..d6ba06ca7a 100644
--- a/django/contrib/gis/geos/__init__.py
+++ b/django/contrib/gis/geos/__init__.py
@@ -1,4 +1,36 @@
-from GEOSGeometry import GEOSGeometry, GEOSException
+"""
+  The goal of this module is to be a ctypes wrapper around the GEOS library
+  that will work on both *NIX and Windows systems.  Specifically, this uses
+  the GEOS C api.
+
+  I have several motivations for doing this:
+    (1) The GEOS SWIG wrapper is no longer maintained, and requires the
+        installation of SWIG.
+    (2) The PCL implementation is over 2K+ lines of C and would make
+        PCL a requisite package for the GeoDjango application stack.
+    (3) Windows and Mac compatibility becomes substantially easier, and does not
+        require the additional compilation of PCL or GEOS and SWIG -- all that
+        is needed is a Win32 or Mac compiled GEOS C library (dll or dylib)
+        in a location that Python can read (e.g. 'C:\Python25').
+
+  In summary, I wanted to wrap GEOS in a more maintainable and portable way using
+   only Python and the excellent ctypes library (now standard in Python 2.5).
+
+  In the spirit of loose coupling, this library does not require Django or
+   GeoDjango.  Only the GEOS C library and ctypes are needed for the platform
+   of your choice.
+
+  For more information about GEOS:
+    http://geos.refractions.net
+  
+  For more info about PCL and the discontinuation of the Python GEOS
+   library see Sean Gillies' writeup (and subsequent update) at:
+     http://zcologia.com/news/150/geometries-for-python/
+     http://zcologia.com/news/429/geometries-for-python-update/
+"""
+
+from GEOSGeometry import GEOSGeometry, Point, LineString, LinearRing, HAS_NUMPY
+from GEOSError import GEOSException
 
 def hex_to_wkt(hex):
     "Converts HEXEWKB into WKT."
diff --git a/django/contrib/gis/geos/libgeos.py b/django/contrib/gis/geos/libgeos.py
new file mode 100644
index 0000000000..17c8898cf6
--- /dev/null
+++ b/django/contrib/gis/geos/libgeos.py
@@ -0,0 +1,106 @@
+"""
+  This module houses the ctypes initialization procedures, as well
+  as the notice and error handler function callbacks (get called
+  when an error occurs in GEOS).
+"""
+
+from django.contrib.gis.geos.GEOSError import GEOSException
+from ctypes import \
+     c_char_p, c_int, c_size_t, c_ubyte, pointer, addressof, \
+     CDLL, CFUNCTYPE, POINTER, Structure
+import os, sys
+
+# NumPy supported?
+try:
+    from numpy import array, ndarray
+    HAS_NUMPY = True
+except ImportError:
+    HAS_NUMPY = False
+
+# Setting the appropriate name for the GEOS-C library, depending on which
+# OS and POSIX platform we're running.
+if os.name == 'nt':
+    # Windows NT library
+    lib_name = 'libgeos_c-1.dll'
+elif os.name == 'posix':
+    platform = os.uname()[0] # Using os.uname()
+    if platform in ('Linux', 'SunOS'):
+        # Linux or Solaris shared library
+        lib_name = 'libgeos_c.so'
+    elif platform == 'Darwin':
+        # Mac OSX Shared Library (Thanks Matt!)
+        lib_name = 'libgeos_c.dylib'
+    else:
+        raise GEOSException, 'Unknown POSIX platform "%s"' % platform
+else:
+    raise GEOSException, 'Unsupported OS "%s"' % os.name
+
+# Getting the GEOS C library.  The C interface (CDLL) is used for
+#  both *NIX and Windows.
+# See the GEOS C API source code for more details on the library function calls:
+#  http://geos.refractions.net/ro/doxygen_docs/html/geos__c_8h-source.html
+lgeos = CDLL(lib_name)
+
+# The notice and error handler C function callback definitions.
+#  Supposed to mimic the GEOS message handler (C below):
+#  "typedef void (*GEOSMessageHandler)(const char *fmt, ...);"
+NOTICEFUNC = CFUNCTYPE(None, c_char_p, c_char_p)
+def notice_h(fmt, list, output_h=sys.stdout):
+    output_h.write('GEOS_NOTICE: %s\n' % (fmt % list))
+notice_h = NOTICEFUNC(notice_h)
+
+ERRORFUNC = CFUNCTYPE(None, c_char_p, c_char_p)
+def error_h(fmt, lst, output_h=sys.stderr):
+    try:
+        err_msg = fmt % lst
+    except:
+        err_msg = fmt
+    output_h.write('GEOS_ERROR: %s\n' % err_msg)
+error_h = ERRORFUNC(error_h)
+
+# The initGEOS routine should be called first, however, that routine takes
+#  the notice and error functions as parameters.  Here is the C code that
+#  we need to wrap:
+#  "extern void GEOS_DLL initGEOS(GEOSMessageHandler notice_function, GEOSMessageHandler error_function);"
+lgeos.initGEOS(notice_h, error_h)
+
+class GEOSPointer(object):
+    """The GEOSPointer provides a layer of abstraction to accessing the values returned by
+    GEOS geometry creation routines."""
+
+    ### Python 'magic' routines ###
+    def __init__(self, ptr):
+        "Given a ctypes pointer(c_int)"
+        if isinstance(ptr, int):
+            self._ptr = pointer(c_int(ptr))
+        else:
+            raise TypeError, 'GEOSPointer object must initialize with an integer.'
+        
+    def __call__(self):
+        """If the pointer is NULL, then an exception will be raised, otherwise the
+        address value (a GEOSGeom_ptr) will be returned."""
+        if self.valid: return self.address
+        else: raise GEOSException, 'GEOS pointer no longer valid (was the parent geometry deleted?)'
+
+    ### GEOSPointer Properties ###
+    @property
+    def address(self):
+        "Returns the address of the GEOSPointer (represented as an integer)."
+        return self._ptr.contents.value
+
+    @property
+    def valid(self):
+        "Tests whether this GEOSPointer is valid."
+        if bool(self.address): return True
+        else: return False
+    
+    ### GEOSPointer Methods ###
+    def set(self, address):
+        "Sets this pointer with the new address (represented as an integer)"
+        if not isinstance(address, int):
+            raise TypeError, 'GEOSPointer must be set with an address (an integer).'
+        self._ptr.contents = c_int(address)
+
+    def nullify(self):
+        "Nullify this geometry pointer (set the address to 0)."
+        self.set(0)
diff --git a/django/contrib/gis/tests/geometries.py b/django/contrib/gis/tests/geometries.py
index ae28d81929..40b73e4737 100644
--- a/django/contrib/gis/tests/geometries.py
+++ b/django/contrib/gis/tests/geometries.py
@@ -90,7 +90,11 @@ multipoints = (TestGeom('MULTIPOINT(10 10, 20 20 )', n_p=2, points=((10., 10.),
 
 # LineStrings
 linestrings = (TestGeom('LINESTRING (60 180, 120 100, 180 180)', n_p=3, centroid=(120, 140), tup=((60, 180), (120, 100), (180, 180))),
-               TestGeom('LINESTRING (0 0, 5 5, 10 5, 10 10)', n_p=4, centroid=(6.1611652351681556, 4.6966991411008934), tup=((0, 0), (5, 5), (10, 5), (10, 10)),)
+               TestGeom('LINESTRING (0 0, 5 5, 10 5, 10 10)', n_p=4, centroid=(6.1611652351681556, 4.6966991411008934), tup=((0, 0), (5, 5), (10, 5), (10, 10)),),
+               )
+
+# Linear Rings
+linearrings = (TestGeom('LINEARRING (649899.3065171393100172 4176512.3807915160432458, 649902.7294133581453934 4176512.7834989596158266, 649906.5550170192727819 4176514.3942507002502680, 649910.5820134161040187 4176516.0050024418160319, 649914.4076170771149918 4176518.0184616246260703, 649917.2264131171396002 4176519.4278986593708396, 649920.0452871860470623 4176521.6427505780011415, 649922.0587463703704998 4176522.8507948759943247, 649924.2735982896992937 4176524.4616246484220028, 649926.2870574744883925 4176525.4683542405255139, 649927.8978092158213258 4176526.8777912775985897, 649929.3072462501004338 4176528.0858355751261115, 649930.1126611357321963 4176529.4952726080082357, 649927.4951798024121672 4176506.9444361114874482, 649899.3065171393100172 4176512.3807915160432458)', n_p=15),
                )
 
 # MultiLineStrings
@@ -128,5 +132,12 @@ relate_geoms = ( (TestGeom('MULTIPOINT(80 70, 20 20, 200 170, 140 120)'),
                   TestGeom('MULTILINESTRING((90 20, 170 100, 170 140), (130 140, 130 60, 90 20, 20 90, 90 20))'),
                   'FF10F0102', True,),
                  )
-#((TestGeom('POLYGON ((545 317, 617 379, 581 321, 545 317))')),
-# (TestGeom('POLYGON ((484 290, 558 359, 543 309, 484 290))'))),
+
+buffer_geoms = ( (TestGeom('POINT(0 0)'),
+                  TestGeom('POLYGON ((5 0,4.903926402016153 -0.97545161008064,4.619397662556435 -1.913417161825447,4.157348061512728 -2.777851165098009,3.53553390593274 -3.535533905932735,2.777851165098015 -4.157348061512724,1.913417161825454 -4.619397662556431,0.975451610080648 -4.903926402016151,0.000000000000008 -5.0,-0.975451610080632 -4.903926402016154,-1.913417161825439 -4.619397662556437,-2.777851165098002 -4.157348061512732,-3.53553390593273 -3.535533905932746,-4.157348061512719 -2.777851165098022,-4.619397662556429 -1.913417161825462,-4.903926402016149 -0.975451610080656,-5.0 -0.000000000000016,-4.903926402016156 0.975451610080624,-4.619397662556441 1.913417161825432,-4.157348061512737 2.777851165097995,-3.535533905932752 3.535533905932723,-2.777851165098029 4.157348061512714,-1.913417161825468 4.619397662556426,-0.975451610080661 4.903926402016149,-0.000000000000019 5.0,0.975451610080624 4.903926402016156,1.913417161825434 4.61939766255644,2.777851165097998 4.157348061512735,3.535533905932727 3.535533905932748,4.157348061512719 2.777851165098022,4.619397662556429 1.91341716182546,4.90392640201615 0.975451610080652,5 0))'),
+                  5.0, 8),
+                 (TestGeom('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'),
+                  TestGeom('POLYGON ((-2 0,-2 10,-1.961570560806461 10.390180644032258,-1.847759065022573 10.765366864730179,-1.662939224605091 11.111140466039204,-1.414213562373095 11.414213562373096,-1.111140466039204 11.662939224605092,-0.765366864730179 11.847759065022574,-0.390180644032256 11.961570560806461,0 12,10 12,10.390180644032256 11.961570560806461,10.765366864730179 11.847759065022574,11.111140466039204 11.66293922460509,11.414213562373096 11.414213562373096,11.66293922460509 11.111140466039204,11.847759065022574 10.765366864730179,11.961570560806461 10.390180644032256,12 10,12 0,11.961570560806461 -0.390180644032256,11.847759065022574 -0.76536686473018,11.66293922460509 -1.111140466039204,11.414213562373096 -1.414213562373095,11.111140466039204 -1.66293922460509,10.765366864730179 -1.847759065022573,10.390180644032256 -1.961570560806461,10 -2,0.0 -2.0,-0.390180644032255 -1.961570560806461,-0.765366864730177 -1.847759065022575,-1.1111404660392 -1.662939224605093,-1.41421356237309 -1.4142135623731,-1.662939224605086 -1.111140466039211,-1.84775906502257 -0.765366864730189,-1.961570560806459 -0.390180644032268,-2 0))'),
+                  2.0, 8),
+                 )
+
diff --git a/django/contrib/gis/tests/test_geos.py b/django/contrib/gis/tests/test_geos.py
index e2955a8a66..b356cac208 100644
--- a/django/contrib/gis/tests/test_geos.py
+++ b/django/contrib/gis/tests/test_geos.py
@@ -1,39 +1,39 @@
 import unittest
-from copy import copy
-from django.contrib.gis.geos import GEOSGeometry, GEOSException
+from django.contrib.gis.geos import GEOSGeometry, GEOSException, Point, LineString, LinearRing, HAS_NUMPY
 from geometries import *
 
-class GeosTest2(unittest.TestCase):
+if HAS_NUMPY:
+    from numpy import array
 
-    def test0100_wkt(self):
+class GEOSTest(unittest.TestCase):
+
+    def test01a_wkt(self):
         "Testing WKT output."
         for g in wkt_out:
             geom = GEOSGeometry(g.wkt)
             self.assertEqual(g.ewkt, geom.wkt)
 
-    def test0101_hex(self):
+    def test01b_hex(self):
         "Testing HEX output."
         for g in hex_wkt:
             geom = GEOSGeometry(g.wkt)
             self.assertEqual(g.hex, geom.hex)
 
-    def test0102_errors(self):
+    def test01c_errors(self):
         "Testing the Error handlers."
-
-        print "\nBEGIN - expecting ParseError; safe to ignore.\n"
+        print "\nBEGIN - expecting GEOS_ERROR; safe to ignore.\n"
         for err in errors:
             if err.hex:
                 self.assertRaises(GEOSException, GEOSGeometry, err.wkt, 'hex')
             else:
                 self.assertRaises(GEOSException, GEOSGeometry, err.wkt)
-        print "\nEND - expecting ParseError; safe to ignore.\n"
-
+        print "\nEND - expecting GEOS_ERROR; safe to ignore.\n"
                 
-    def test02_points(self):
+    def test02a_points(self):
         "Testing Point objects."
         prev = GEOSGeometry('POINT(0 0)')
-
         for p in points:
+            # Creating the point from the WKT
             pnt = GEOSGeometry(p.wkt)
             self.assertEqual(pnt.geom_type, 'Point')
             self.assertEqual(pnt.geom_typeid, 0)
@@ -41,22 +41,39 @@ class GeosTest2(unittest.TestCase):
             self.assertEqual(p.y, pnt.y)
             self.assertEqual(True, pnt == GEOSGeometry(p.wkt))
             self.assertEqual(False, pnt == prev)
-            prev = pnt
 
+            # Making sure that the point's X, Y components are what we expect
             self.assertAlmostEqual(p.x, pnt.tuple[0], 9)
             self.assertAlmostEqual(p.y, pnt.tuple[1], 9)
+
+            # Testing the third dimension, and getting the tuple arguments
             if hasattr(p, 'z'):
                 self.assertEqual(True, pnt.hasz)
                 self.assertEqual(p.z, pnt.z)
                 self.assertEqual(p.z, pnt.tuple[2], 9)
+                tup_args = (p.x, p.y, p.z)
             else:
                 self.assertEqual(False, pnt.hasz)
                 self.assertEqual(None, pnt.z)
+                tup_args = (p.x, p.y)
 
-                
+            # Centroid operation on point should be point itself
             self.assertEqual(p.centroid, pnt.centroid.tuple)
 
-    def test02_multipoints(self):
+            # Now testing the different constructors
+            pnt2 = Point(tup_args)  # e.g., Point((1, 2))
+            pnt3 = Point(*tup_args) # e.g., Point(1, 2)
+            self.assertEqual(True, pnt == pnt2)
+            self.assertEqual(True, pnt == pnt3)
+
+            # Now testing setting the x and y
+            pnt.y = 3.14
+            pnt.x = 2.71
+            self.assertEqual(3.14, pnt.y)
+            self.assertEqual(2.71, pnt.x)
+            prev = pnt # setting the previous geometry
+
+    def test02b_multipoints(self):
         "Testing MultiPoint objects."
         for mp in multipoints:
             mpnt = GEOSGeometry(mp.wkt)
@@ -73,81 +90,38 @@ class GeosTest2(unittest.TestCase):
                 self.assertEqual(p.geom_typeid, 0)
                 self.assertEqual(p.empty, False)
                 self.assertEqual(p.valid, True)
-    
-    def test03_polygons(self):
-        "Testing Polygon objects."
 
-        prev = GEOSGeometry('POINT(0 0)')
-
-        for p in polygons:
-            poly = GEOSGeometry(p.wkt)
-            self.assertEqual(poly.geom_type, 'Polygon')
-            self.assertEqual(poly.geom_typeid, 3)
-            self.assertEqual(poly.empty, False)
-            self.assertEqual(poly.ring, False)
-            self.assertEqual(p.n_i, poly.num_interior_rings)
-            self.assertEqual(p.n_i + 1, len(poly)) # Testing __len__
-            self.assertEqual(p.n_p, poly.num_points)
-
-            # Area & Centroid
-            self.assertAlmostEqual(p.area, poly.area, 9)
-            self.assertAlmostEqual(p.centroid[0], poly.centroid.tuple[0], 9)
-            self.assertAlmostEqual(p.centroid[1], poly.centroid.tuple[1], 9)
-
-            # Testing the geometry equivalence
-            self.assertEqual(True, poly == GEOSGeometry(p.wkt))
-            self.assertEqual(False, poly == prev)
-            prev = poly
-
-            # Testing the exterior ring
-            ring = poly.exterior_ring
-            self.assertEqual(ring.geom_type, 'LinearRing')
-            self.assertEqual(ring.geom_typeid, 2)
-            if p.ext_ring_cs:
-                self.assertEqual(p.ext_ring_cs, ring.tuple)
-                self.assertEqual(p.ext_ring_cs, poly[0].tuple) # Testing __getitem__
-
-    def test03_multipolygons(self):
-        "Testing MultiPolygon objects."
-
-        prev = GEOSGeometry('POINT (0 0)')
-
-        for mp in multipolygons:
-            mpoly = GEOSGeometry(mp.wkt)
-            self.assertEqual(mpoly.geom_type, 'MultiPolygon')
-            self.assertEqual(mpoly.geom_typeid, 6)
-            self.assertEqual(mp.valid, mpoly.valid)
-
-            if mp.valid:
-                self.assertEqual(mp.num_geom, mpoly.num_geom)
-                self.assertEqual(mp.n_p, mpoly.num_coords)
-                self.assertEqual(mp.num_geom, len(mpoly))
-                for p in mpoly:
-                    self.assertEqual(p.geom_type, 'Polygon')
-                    self.assertEqual(p.geom_typeid, 3)
-                    self.assertEqual(p.valid, True)
-
-    def test04_linestring(self):
+    def test03a_linestring(self):
         "Testing LineString objects."
-
         prev = GEOSGeometry('POINT(0 0)')
-
         for l in linestrings:
             ls = GEOSGeometry(l.wkt)
             self.assertEqual(ls.geom_type, 'LineString')
             self.assertEqual(ls.geom_typeid, 1)
             self.assertEqual(ls.empty, False)
             self.assertEqual(ls.ring, False)
-            self.assertEqual(l.centroid, ls.centroid.tuple)
+            if hasattr(l, 'centroid'):
+                self.assertEqual(l.centroid, ls.centroid.tuple)
+            if hasattr(l, 'tup'):
+                self.assertEqual(l.tup, ls.tuple)
+                
             self.assertEqual(True, ls == GEOSGeometry(l.wkt))
             self.assertEqual(False, ls == prev)
+
             prev = ls
 
-    def test04_multilinestring(self):
+            # Creating a LineString from a tuple, list, and numpy array
+            ls2 = LineString(ls.tuple)
+            self.assertEqual(ls, ls2)
+            ls3 = LineString([list(tup) for tup in ls.tuple])
+            self.assertEqual(ls, ls3)
+            if HAS_NUMPY:
+                ls4 = LineString(array(ls.tuple))
+                self.assertEqual(ls, ls4)
+
+    def test03b_multilinestring(self):
         "Testing MultiLineString objects."
-
         prev = GEOSGeometry('POINT(0 0)')
-
         for l in multilinestrings:
             ml = GEOSGeometry(l.wkt)
             self.assertEqual(ml.geom_type, 'MultiLineString')
@@ -165,14 +139,152 @@ class GeosTest2(unittest.TestCase):
                 self.assertEqual(ls.geom_typeid, 1)
                 self.assertEqual(ls.empty, False)
 
-    #def test05_linearring(self):
-    #    "Testing LinearRing objects."
-    #    pass
+    def test04a_linearring(self):
+        "Testing LinearRing objects."
+        for rr in linearrings:
+            lr = GEOSGeometry(rr.wkt)
+            self.assertEqual(lr.geom_type, 'LinearRing')
+            self.assertEqual(lr.geom_typeid, 2)
+            self.assertEqual(rr.n_p, len(lr))
+            self.assertEqual(True, lr.valid)
+            self.assertEqual(False, lr.empty)
 
+            # Creating a LinearRing from a tuple, list, and numpy array
+            lr2 = LinearRing(lr.tuple)
+            self.assertEqual(lr, lr2)
+            lr3 = LinearRing([list(tup) for tup in lr.tuple])
+            self.assertEqual(lr, lr3)
+            if HAS_NUMPY:
+                lr4 = LineString(array(lr.tuple))
+                self.assertEqual(lr, lr4)
+    
+    def test05a_polygons(self):
+        "Testing Polygon objects."
+        prev = GEOSGeometry('POINT(0 0)')
+        for p in polygons:
+            # Creating the Polygon, testing its properties.
+            poly = GEOSGeometry(p.wkt)
+            self.assertEqual(poly.geom_type, 'Polygon')
+            self.assertEqual(poly.geom_typeid, 3)
+            self.assertEqual(poly.empty, False)
+            self.assertEqual(poly.ring, False)
+            self.assertEqual(p.n_i, poly.num_interior_rings)
+            self.assertEqual(p.n_i + 1, len(poly)) # Testing __len__
+            self.assertEqual(p.n_p, poly.num_points)
+
+            # Area & Centroid
+            self.assertAlmostEqual(p.area, poly.area, 9)
+            self.assertAlmostEqual(p.centroid[0], poly.centroid.tuple[0], 9)
+            self.assertAlmostEqual(p.centroid[1], poly.centroid.tuple[1], 9)
+
+            # Testing the geometry equivalence
+            self.assertEqual(True, poly == GEOSGeometry(p.wkt))
+            self.assertEqual(False, poly == prev) # Should not be equal to previous geometry
+
+            # Testing the exterior ring
+            ring = poly.exterior_ring
+            self.assertEqual(ring.geom_type, 'LinearRing')
+            self.assertEqual(ring.geom_typeid, 2)
+            if p.ext_ring_cs:
+                self.assertEqual(p.ext_ring_cs, ring.tuple)
+                self.assertEqual(p.ext_ring_cs, poly[0].tuple) # Testing __getitem__
+
+            # Testing __iter__
+            for r in poly:
+                self.assertEqual(ring.geom_type, 'LinearRing')
+                self.assertEqual(ring.geom_typeid, 2)
+
+            # Setting the second point of the first ring (which should set the
+            #  first point of the polygon).
+            prev = poly.clone() # Using clone() to get a copy of the current polygon
+            self.assertEqual(True, poly == prev) # They clone should be equal to the first
+            newval = (poly[0][1][0] + 5.0, poly[0][1][1] + 5.0) # really testing __getitem__ ([ring][point][tuple])
+            try:
+                poly[0][1] = ('cannot assign with', 'string values')
+            except TypeError:
+                pass
+            poly[0][1] = newval # setting the second point in the polygon with the newvalue (based on the old)
+            self.assertEqual(newval, poly[0][1]) # The point in the polygon should be the
+            self.assertEqual(False, poly == prev) # Even different from the clone we just made
+            
+    def test05b_multipolygons(self):
+        "Testing MultiPolygon objects."
+        print "\nBEGIN - expecting GEOS_NOTICE; safe to ignore.\n"
+        prev = GEOSGeometry('POINT (0 0)')
+        for mp in multipolygons:
+            mpoly = GEOSGeometry(mp.wkt)
+            self.assertEqual(mpoly.geom_type, 'MultiPolygon')
+            self.assertEqual(mpoly.geom_typeid, 6)
+            self.assertEqual(mp.valid, mpoly.valid)
+
+            if mp.valid:
+                self.assertEqual(mp.num_geom, mpoly.num_geom)
+                self.assertEqual(mp.n_p, mpoly.num_coords)
+                self.assertEqual(mp.num_geom, len(mpoly))
+                for p in mpoly:
+                    self.assertEqual(p.geom_type, 'Polygon')
+                    self.assertEqual(p.geom_typeid, 3)
+                    self.assertEqual(p.valid, True)
+        print "\nEND - expecting GEOS_NOTICE; safe to ignore.\n"  
+
+    def test06_memory_hijinks(self):
+        "Testing Geometry __del__() in different scenarios"
+        #### Memory issues with rings and polygons
+
+        # These tests are needed to ensure sanity with writable geometries.
+
+        # Getting a polygon with interior rings, and pulling out the interior rings
+        poly = GEOSGeometry(polygons[1].wkt)
+        ring1 = poly[0]
+        ring2 = poly[1]
+
+        # These deletes should be 'harmless' since they are done on child geometries
+        del ring1 
+        del ring2
+        ring1 = poly[0]
+        ring2 = poly[1]
+
+        # Deleting the polygon
+        del poly
+
+        # Ensuring that trying to access the deleted memory (by getting the string
+        #  representation of the ring of a deleted polygon) raises a GEOSException
+        #  instead of something worse..
+        self.assertRaises(GEOSException, str, ring1)
+        self.assertRaises(GEOSException, str, ring2)
+
+        #### Memory issues with geometries from Geometry Collections
+        mp = GEOSGeometry('MULTIPOINT(85 715, 235 1400, 4620 1711)')
+        
+        # Getting the points
+        pts = [p for p in mp]
+
+        # More 'harmless' child geometry deletes
+        for p in pts: del p
+
+        # Cloning for comparisons
+        clones = [p.clone() for p in pts]
+
+        for i in xrange(len(clones)):
+            # Testing equivalence before & after modification
+            self.assertEqual(True, pts[i] == clones[i]) # before
+            pts[i].x = 3.14159
+            pts[i].y = 2.71828
+            self.assertEqual(False, pts[i] == clones[i]) # after
+            self.assertEqual(3.14159, mp[i].x) # parent x,y should be modified
+            self.assertEqual(2.71828, mp[i].y)
+
+        # Should raise GEOSException when trying to get geometries from the multipoint
+        #  after it has been deleted.
+        del mp
+        for p in pts:
+            self.assertRaises(GEOSException, str, p)
+        
     def test08_coord_seq(self):
         "Testing Coordinate Sequence objects."
         for p in polygons:
             if p.ext_ring_cs:
+                # Constructing the polygon and getting the coordinate sequence
                 poly = GEOSGeometry(p.wkt)
                 cs = poly.exterior_ring.coord_seq
 
@@ -181,15 +293,19 @@ class GeosTest2(unittest.TestCase):
 
                 # Checks __getitem__ and __setitem__
                 for i in xrange(len(p.ext_ring_cs)):
-                    c1 = p.ext_ring_cs[i]
-                    c2 = cs[i]
+                    c1 = p.ext_ring_cs[i] # Expected value
+                    c2 = cs[i] # Value from coordseq
                     self.assertEqual(c1, c2)
+
+                    # Constructing the test value to set the coordinate sequence with
                     if len(c1) == 2: tset = (5, 23)
                     else: tset = (5, 23, 8)
                     cs[i] = tset
+                    
+                    # Making sure every set point matches what we expect
                     for j in range(len(tset)):
                         cs[i] = tset
-                        self.assertEqual(tset, cs[i])
+                        self.assertEqual(tset[j], cs[i][j])
 
     def test09_relate_pattern(self):
         "Testing relate() and relate_pattern()."
@@ -208,7 +324,6 @@ class GeosTest2(unittest.TestCase):
 
     def test10_intersection(self):
         "Testing intersects() and intersection()."
-
         for i in xrange(len(topology_geoms)):
             g_tup = topology_geoms[i]
             a = GEOSGeometry(g_tup[0].wkt)
@@ -239,10 +354,36 @@ class GeosTest2(unittest.TestCase):
             d2 = a.difference(b)
             self.assertEqual(d1, d2)
 
+    def test13_buffer(self):
+        "Testing buffer()."
+        for i in xrange(len(buffer_geoms)):
+            g_tup = buffer_geoms[i]
+            g = GEOSGeometry(g_tup[0].wkt)
+
+            # The buffer we expect
+            exp_buf = GEOSGeometry(g_tup[1].wkt)
+
+            # Can't use a floating-point for the number of quadsegs.
+            self.assertRaises(TypeError, g.buffer, g_tup[2], float(g_tup[3]))
+
+            # Constructing our buffer
+            buf = g.buffer(g_tup[2], g_tup[3])
+            self.assertEqual(exp_buf.num_coords, buf.num_coords)
+            self.assertEqual(len(exp_buf), len(buf))
+
+            # Now assuring that each point in the buffer is almost equal
+            for j in xrange(len(exp_buf)):
+                exp_ring = exp_buf[j]
+                buf_ring = buf[j]
+                self.assertEqual(len(exp_ring), len(buf_ring))
+                for k in xrange(len(exp_ring)):
+                    # Asserting the X, Y of each point are almost equal (due to floating point imprecision)
+                    self.assertAlmostEqual(exp_ring[k][0], buf_ring[k][0], 9)
+                    self.assertAlmostEqual(exp_ring[k][1], buf_ring[k][1], 9)
 
 def suite():
     s = unittest.TestSuite()
-    s.addTest(unittest.makeSuite(GeosTest2))
+    s.addTest(unittest.makeSuite(GEOSTest))
     return s
 
 def run(verbosity=2):