How to use the pysal.cg function in pysal

def internal_or_external(polys,filament, vertices):
    #Modification of Serge's code to find poly in poly for line in poly
    #pl = ps.cg.PolygonLocator(polys) #Trying to use this for polyline in polygon
    
    #Spatial Index of Filaments and minimal cycles (polygons)
    polyline = ps.cg.Chain([ps.cg.Point(vertices[pnt]) for pnt in filament])
    polyline_mbr =  polyline.bounding_box
    pl = ps.cg.PolygonLocator(polys)
    overlaps = pl.overlapping(polyline_mbr)
    
    #For the overlapping MBRs check to see if the polyline is internal or external to the min cycle
    for k in range(len(overlaps)):
        s = sum(overlaps[k].contains_point(v) for v in polyline.vertices)
        if s == len(polyline.vertices):
            #Internal
            return True
        else:
            return False

if numadj == 0:
                nodes, node_coord, primitives, vertices, ext_edges = extractisolated(nodes,node_coord,start_key,primitives, vertices, ext_edges)
                sorted_nodes.pop(0)
            elif numadj == 1:
                sorted_nodes, edges, nodes, node_coord, primitives, vertices, ext_edges = extractfilament(start_key, adj_nodes(start_key, edges)[0],nodes, node_coord, sorted_nodes, edges,primitives,cycle_edge, vertices, ext_edges)
            else:
                sorted_nodes, edges, nodes, node_coord, primitives, minimal_cycles,cycle_edge, vertices, ext_edges = extract_primitives(start_key,sorted_nodes, edges, nodes, node_coord, primitives, minimal_cycles,cycle_edge, vertices, ext_edges)

        #4. Remove holes from the graph
        if remove_holes == True:
            polys = []
            for cycle in minimal_cycles:
                polys.append(ps.cg.Polygon([ps.cg.Point(vertices[pnt]) for pnt in cycle]))

            pl = ps.cg.PolygonLocator(polys)

            # find all overlapping polygon mbrs
            overlaps ={}
            nump = len(minimal_cycles)
            for i in range(nump):
                overlaps[i] = pl.overlapping(polys[i].bounding_box)

            # for overlapping mbrs (left,right) check if right polygon is contained in left
            holes = []
            for k in overlaps:
                for  pc in overlaps[k]:
                    s = sum( [polys[k].contains_point(v) for v in pc.vertices])
                    if s == len(pc.vertices):
                        # print k, pc
                        holes.append((k,pc))

def net_global_stats(G, boundary=None, detour=True):
    v = no_nodes(G)
    e = no_edges(G)
    L = tot_net_length(G)
    p = no_components(G)
    u = e - v + p # cyclomatic number
    alpha = u*1.0/(2*v - 5)
    beta = e*1.0/v
    emax = 3*(v-2)
    gamma = e*1.0/emax
    eta = L*1.0/e
    net_den = None
    if boundary:
        s = pysal.open(boundary)
        if s.type != pysal.cg.shapes.Polygon: 
            raise ValueError, 'File is not of type POLYGON'
        net_den = s.next().area
    net_dist, eucl_dist = 0.0, 0.0
    det = None
    if detour:
        nodes = G.keys()
        for n in nodes:
            net_D = dijkstras(G, n)
            net_dist += sum(net_D.values())
            eucl_D = [ math.sqrt((n[0] - m[0])**2 + (n[1] - m[1])**2) for m in nodes]
            eucl_dist += sum(eucl_D)
        net_dist /= 2.0
        eucl_dist /= 2.0
        if net_dist > 0.0:
            det = eucl_dist*1.0/net_dist
    return v, e, L, p, u, alpha, beta, emax, gamma, eta, net_den, det

dist_to_node:   dict
                        Dict with point ids as keys and values as dicts with keys for
                        node ids and values as distances from point to node.
        """

        obs_to_edge = {}
        dist_to_node = {}

        pointpattern.snapped_coordinates = {}
        segments = []
        s2e = {}
        for edge in self.edges:
            head = self.node_coords[edge[0]]
            tail = self.node_coords[edge[1]]
            segments.append(ps.cg.Chain([head,tail]))
            s2e[(head,tail)] = edge
            

        points = {}
        p2id = {}
        for pointIdx, point in pointpattern.points.iteritems(): 
            points[pointIdx] = point['coordinates']

        snapped = util.snapPointsOnSegments(points, segments)

        for pointIdx, snapInfo in snapped.iteritems():
            x,y = snapInfo[1].tolist()
            edge = s2e[tuple(snapInfo[0])]
            if edge not in obs_to_edge:
                obs_to_edge[edge] = {}
            obs_to_edge[edge][pointIdx] = (x,y)

Examples
    --------
    >>> points=[(10, 10), (20, 10), (40, 10), (15, 20), (30, 20), (30, 30)]
    >>> wcheck = pysal.W({0: [1, 3], 1: [0, 3, ], 2: [], 3: [1, 0], 4: [5], 5: [4]})
    WARNING: there is one disconnected observation (no neighbors)
    Island id:  [2]
    >>> w=threshold_binaryW_from_array(points,threshold=11.2)
    WARNING: there is one disconnected observation (no neighbors)
    Island id:  [2]
    >>> pysal.weights.util.neighbor_equality(w, wcheck)
    True

    >>>
    """
    if radius is not None:
        array = pysal.cg.KDTree(array, distance_metric='Arc', radius=radius)
    return DistanceBand(array, threshold=threshold, p=p)

** Note **
                        each element is a segment represented as a chain with 
                        *one head and one tail node*, in other words one link only.
              
    Returns
    ------- 
    p2s:        dict
                key:    point id (see points in arguments)
         
                value:  a 2-tuple: ((head, tail), point)
                        where (head, tail) is the target segment, and point is the snapped
                        location on the segment
    """
    
    # Put segments in an Rtree.
    rt = ps.cg.Rtree()
    SMALL = 0.01
    node2segs = {}
    
    for segment in segments:
        head,tail = segment.vertices
        x0,y0 = head
        x1,y1 = tail
        if (x0,y0) not in node2segs:
            node2segs[(x0,y0)] = []
        if (x1,y1) not in node2segs:
            node2segs[(x1,y1)] = []
        node2segs[(x0,y0)].append(segment)
        node2segs[(x1,y1)].append(segment)
        x0,y0,x1,y1 =  segment.bounding_box
        x0 -= SMALL
        y0 -= SMALL

start_key = sorted_nodes[0][0]
            numadj = len(adj_nodes(start_key, edges))

            if numadj == 0:
                nodes, node_coord, primitives, vertices, ext_edges = extractisolated(nodes,node_coord,start_key,primitives, vertices, ext_edges)
                sorted_nodes.pop(0)
            elif numadj == 1:
                sorted_nodes, edges, nodes, node_coord, primitives, vertices, ext_edges = extractfilament(start_key, adj_nodes(start_key, edges)[0],nodes, node_coord, sorted_nodes, edges,primitives,cycle_edge, vertices, ext_edges)
            else:
                sorted_nodes, edges, nodes, node_coord, primitives, minimal_cycles,cycle_edge, vertices, ext_edges = extract_primitives(start_key,sorted_nodes, edges, nodes, node_coord, primitives, minimal_cycles,cycle_edge, vertices, ext_edges)

        #4. Remove holes from the graph
        if remove_holes == True:
            polys = []
            for cycle in minimal_cycles:
                polys.append(ps.cg.Polygon([ps.cg.Point(vertices[pnt]) for pnt in cycle]))

            pl = ps.cg.PolygonLocator(polys)

            # find all overlapping polygon mbrs
            overlaps ={}
            nump = len(minimal_cycles)
            for i in range(nump):
                overlaps[i] = pl.overlapping(polys[i].bounding_box)

            # for overlapping mbrs (left,right) check if right polygon is contained in left
            holes = []
            for k in overlaps:
                for  pc in overlaps[k]:
                    s = sum( [polys[k].contains_point(v) for v in pc.vertices])
                    if s == len(pc.vertices):
                        # print k, pc

#edge = (rec['FNODE'],rec['TNODE'])
        #TODO: Calculate location along edge and distance to edge"
        #return edge
        return x2,y2
        
    
if __name__=='__main__':
    import random
    net = SpatialNetwork('beth_network.shp',ARC_DISTANCE)

    n = 1000
    minX,minY,maxX,maxY = net.shp.bbox
    xRange = maxX-minX
    yRange = maxY-minY
    qpts = [(random.random(), random.random()) for i in xrange(n)]
    qpts = [pysal.cg.Point((minX+(xRange*x),minY+(yRange*y))) for x,y in qpts]
    o = pysal.open('random_qpts.shp','w')
    for p in qpts:
        o.write(p)
    o.close()
    o = pysal.open('random_qpts_snapped.shp','w')
    for qpt in qpts:
        spt = net.snap(qpt)
        o.write(pysal.cg.Chain([qpt,spt]))
    o.close()

Arguments
    ---------
    tolerance -- float -- snapping tolerance in meters
    arc -- bool -- If true, Ard Distance will be used instead of Euclidean

    Returns
    -------
    generator -- each element is a new pysal.cg.Chain with corrected vertices.
    """
    kmtol = tolerance/1000.

    data = numpy.concatenate([rec.vertices for rec in shp])
    
    if arc:
        kd = pysal.cg.KDTree(data,distance_metric="Arc",radius = pysal.cg.sphere.RADIUS_EARTH_KM)
    else:
        kd = pysal.cg.KDTree(data)
    q = kd.query_ball_tree(kd,kmtol)
    ### Next three lines assert that snappings are mutual... if 1 snaps to 8, 8 must snap to 1.
    for r,a in enumerate(q):
        for o in a:
            assert a==q[o]
    ### non-mutual snapping can happen.
    ### consider the three points, A (-1,0), B (0,0), C (1,0) and a snapping tolerance of 1.
    ### A-> B
    ### B-> A,C
    ### C-> B
    ### For now, try lowering adjusting the tolerance to avoid this.

    data2 = numpy.empty_like(data)
    for i,r in enumerate(q):