from ctypes import string_at, byref, c_int, c_size_t, c_char_p, c_double
    #, create_string_buffer \
    #, c_char_p, c_double, c_float, c_int, c_uint, c_size_t, c_ubyte \
    #, c_void_p, byref

from shapely.geos import lgeos
from shapely.predicates import BinaryPredicate, UnaryPredicate
from shapely.topology import BinaryTopologicalOp, UnaryTopologicalOp


def geom_factory(g):
    ob = BaseGeometry()
    geom_type = string_at(lgeos.GEOSGeomType(g))
    mod = __import__(
        'shapely.geometry', 
        globals(), 
        locals(), 
        [geom_type],
        )
    ob.__class__ = getattr(mod, geom_type)
    ob._geom = g
    return ob


class BaseGeometry(object):
    
    """Defines methods common to all geometries.
    """

    _geom = None
    _ctypes_data = None
    _ndim = None
    _crs = None

    def __init__(self):
        self._geom = lgeos.GEOSGeomFromWKT(c_char_p('GEOMETRYCOLLECTION EMPTY'))

    def __del__(self):
        if self._geom is not None:
            lgeos.GEOSGeom_destroy(self._geom)

    def __str__(self):
        return self.to_wkt()

    def __reduce__(self):
        return (self.__class__, (), self.to_wkb())

    def __setstate__(self, state):
        self._geom = lgeos.GEOSGeomFromWKB_buf(
                        c_char_p(state), 
                        c_size_t(len(state))
                        )

    # Array and ctypes interfaces

    @property
    def array(self):
        """Return a GeoJSON coordinate array.
        
        To be overridden by extension classes."""
        raise NotImplemented

    @property
    def ctypes(self):
        """Return a ctypes representation.
        
        To be overridden by extension classes."""
        raise NotImplemented

    @property
    def __array_interface__(self):
        """Provide the Numpy array protocol."""
        return {
            'version': 3,
            'shape': (self._ndim,),
            'typestr': '>f8',
            'data': self.ctypes
            }

    # Python feature protocol

    @property
    def type(self):
        return self.geometryType()

    @property
    def coordinates(self):
        return self.array


    def geometryType(self):
        """Returns a string representing the geometry type, e.g. 'Polygon'."""
        return string_at(lgeos.GEOSGeomType(self._geom))

    def to_wkb(self):
        """Returns a WKT string representation of the geometry."""
        size = c_int()
        bytes = lgeos.GEOSGeomToWKB_buf(self._geom, byref(size))
        return string_at(bytes, size.value)

    def to_wkt(self):
        """Returns a WKT string representation of the geometry."""
        return string_at(lgeos.GEOSGeomToWKT(self._geom))

    # Properties
    geom_type = property(geometryType)
    wkt = property(to_wkt)
    wkb = property(to_wkb)

    # Topology operations
   
    envelope = UnaryTopologicalOp(lgeos.GEOSEnvelope, geom_factory)
    intersection = BinaryTopologicalOp(lgeos.GEOSIntersection, geom_factory)
    
    def buffer(self, distance, quadsegs=16):
        return geom_factory(
            lgeos.GEOSBuffer(self._geom, c_double(distance), c_int(quadsegs))
            )

    convex_hull = UnaryTopologicalOp(lgeos.GEOSConvexHull, geom_factory)
    difference = BinaryTopologicalOp(lgeos.GEOSDifference, geom_factory)
    symmetric_difference = BinaryTopologicalOp(lgeos.GEOSSymDifference, 
        geom_factory)
    boundary = UnaryTopologicalOp(lgeos.GEOSBoundary, geom_factory)
    union = BinaryTopologicalOp(lgeos.GEOSUnion, geom_factory)
    centroid = UnaryTopologicalOp(lgeos.GEOSGetCentroid, geom_factory)

    def relate(self, other):
        func = lgeos.GEOSRelate
        func.restype = c_char_p
        return lgeos.GEOSRelate(self._geom, other._geom)

    # Binary predicates

    # Relate Pattern (TODO?)
    disjoint = BinaryPredicate(lgeos.GEOSDisjoint)
    touches = BinaryPredicate(lgeos.GEOSTouches)
    intersects = BinaryPredicate(lgeos.GEOSIntersects)
    crosses = BinaryPredicate(lgeos.GEOSCrosses)
    within = BinaryPredicate(lgeos.GEOSWithin)
    contains = BinaryPredicate(lgeos.GEOSContains)
    overlaps = BinaryPredicate(lgeos.GEOSOverlaps)
    equals = BinaryPredicate(lgeos.GEOSEquals)

    # Unary predicates

    is_empty = UnaryPredicate(lgeos.GEOSisEmpty)
    is_valid = UnaryPredicate(lgeos.GEOSisValid)
    is_simple = UnaryPredicate(lgeos.GEOSisSimple)
    is_ring = UnaryPredicate(lgeos.GEOSisRing)
    has_z = UnaryPredicate(lgeos.GEOSHasZ)


@@ -1,12 +1,21 @@
+%22%22%22%0A%22%22%22%0A%0A
 from ctypes 
@@ -78,132 +78,8 @@
 uble
-%0A    #, create_string_buffer %5C%0A    #, c_char_p, c_double, c_float, c_int, c_uint, c_size_t, c_ubyte %5C%0A    #, c_void_p, byref
 %0A%0Afr
@@ -241,16 +241,84 @@
 calOp%0A%0A%0A
+# Abstract geometry factory for use with topological methods below%0A%0A
 def geom
@@ -403,16 +403,70 @@
 ype(g))%0A
+    # TODO: check cost of dynamic import by profiling%0A
     mod 
@@ -2828,20 +2828,91 @@
 rations%0A
+    #%0A    # These use descriptors to reduce the amount of boilerplate.%0A
    %0A
-
     enve
@@ -3052,180 +3052,8 @@
 ry)%0A
-    %0A    def buffer(self, distance, quadsegs=16):%0A        return geom_factory(%0A            lgeos.GEOSBuffer(self._geom, c_double(distance), c_int(quadsegs))%0A            )%0A%0A
     
@@ -3267,16 +3267,55 @@
 rence, %0A
+                                       
         
@@ -3324,24 +3324,24 @@
 om_factory)%0A
-
     boundary
@@ -3531,16 +3531,294 @@
 ctory)%0A%0A
+    # Buffer has a unique distance argument, so not a descriptor%0A    def buffer(self, distance, quadsegs=16):%0A        return geom_factory(%0A            lgeos.GEOSBuffer(self._geom, c_double(distance), c_int(quadsegs))%0A            )%0A%0A    # Relate has a unique string return value%0A
     def 
@@ -3979,24 +3979,95 @@
 y predicates
+%0A    #%0A    # These use descriptors to reduce the amount of boilerplate.
 %0A%0A    # Rela
@@ -4486,16 +4486,16 @@
 quals)%0A%0A
-
     # Un
@@ -4508,16 +4508,87 @@
 edicates
+%0A    #%0A    # These use descriptors to reduce the amount of boilerplate.
 %0A%0A    is

Improve readibility