Python pysal.open函数代码示例

    def __init__(self, path, time_col):
        shp = pysal.open(path + '.shp')
        dbf = pysal.open(path + '.dbf')

        # extract the spatial coordinates from the shapefile
        x = []
        y = []
        n = 0
        for i in shp:
            count = 0
            for j in i:
                if count == 0:
                    x.append(j)
                elif count == 1:
                    y.append(j)
                count += 1
            n += 1

        self.n = n
        x = np.array(x)
        y = np.array(y)
        self.x = np.reshape(x, (n, 1))
        self.y = np.reshape(y, (n, 1))
        self.space = np.hstack((self.x, self.y))

        # extract the temporal information from the database
        t = np.array(dbf.by_col(time_col))
        line = np.ones((n, 1))
        self.t = np.reshape(t, (n, 1))
        self.time = np.hstack((self.t, line))

        # close open objects
        dbf.close()
        shp.close()

def dl_merge(outname="tracts",sumlevel="140"):
    os.chdir('/tmp')
    outshp = pysal.open(outname+'.shp','w')
    outdbf = pysal.open(outname+'.dbf','w')
    for st in config.STATE_FIPS:
        fname = name_template%(st)#,sumlevel)
        if DEBUG: print fname
        url = urllib.urlopen(base_url+fname)
        dat = url.read()
        if not os.path.exists(fname):
            with open(fname,'wb') as o:
                o.write(dat)
            os.system('unzip '+fname)
        shp = pysal.open(fname.replace('.zip','.shp'),'r')
        for x in shp:
            outshp.write(x)
        dbf = pysal.open(fname.replace('.zip','.dbf'),'r')
        outdbf.header = dbf.header
        outdbf.field_spec = dbf.field_spec
        for row in dbf:
            outdbf.write(row)
        os.remove(fname)
        os.remove(fname.replace('.zip','.shp'))
        os.remove(fname.replace('.zip','.shx'))
        os.remove(fname.replace('.zip','.dbf'))
        os.remove(fname.replace('.zip','.shp.xml'))
        os.remove(fname.replace('.zip','.prj'))
    outshp.close()
    outdbf.close()

def test_d3viz():
    import pysal
    shp = pysal.open(pysal.examples.get_path('NAT.shp'),'r')
    dbf = pysal.open(pysal.examples.get_path('NAT.dbf'),'r')    

    #shp = pysal.open('/Users/xun/github/PySAL-Viz/test_data/man_road.shp','r')
    #dbf = pysal.open('/Users/xun/github/PySAL-Viz/test_data/man_road.dbf','r')
    import d3viz
    d3viz.setup()
    
    d3viz.init_map(shp)
    d3viz.show_map(shp)
   
    d3viz.scatter_plot(shp, field_x="HR90", field_y="HC90") 
    w = pysal.queen_from_shapefile(pysal.examples.get_path('NAT.shp'))
    d3viz.moran_scatter_plot(shp, dbf, "HR90", w)
    
    d3viz.quantile_map(shp, "HR90", 5)
    d3viz.quantile_map(shp, "HR90", 5, basemap="leaflet")
   
    import numpy as np
    y = np.array(dbf.by_col["HR90"])
    lm = pysal.Moran_Local(y, w)
    
    d3viz.lisa_map(shp, "HR90", lm)

    def test_chi2(self):
        import pysal

        f = pysal.open(pysal.examples.get_path("usjoin.csv"))
        pci = np.array([f.by_col[str(y)] for y in range(1929, 2010)])
        pci = pci.transpose()
        rpci = pci / (pci.mean(axis=0))
        w = pysal.open(pysal.examples.get_path("states48.gal")).read()
        w.transform = "r"
        sm = pysal.Spatial_Markov(rpci, w, fixed=True, k=5)
        chi = np.matrix(
            [
                [4.06139105e01, 6.32961385e-04, 1.60000000e01],
                [5.55485793e01, 2.88879565e-06, 1.60000000e01],
                [1.77772638e01, 3.37100315e-01, 1.60000000e01],
                [4.00925436e01, 7.54729084e-04, 1.60000000e01],
                [4.68588786e01, 7.16364084e-05, 1.60000000e01],
            ]
        ).getA()
        obs = np.matrix(sm.chi2).getA()
        np.testing.assert_array_almost_equal(obs, chi)
        obs = np.matrix(
            [
                [4.61209613e02, 0.00000000e00, 4.00000000e00],
                [1.48140694e02, 0.00000000e00, 4.00000000e00],
                [6.33129261e01, 5.83089133e-13, 4.00000000e00],
                [7.22778509e01, 7.54951657e-15, 4.00000000e00],
                [2.32659201e02, 0.00000000e00, 4.00000000e00],
            ]
        )
        np.testing.assert_array_almost_equal(obs.getA(), np.matrix(sm.shtest).getA())

    def __init__(self, shapefile, idvariable=None, attribute=False):
        self.points = {}
        self.npoints = 0

        if idvariable:
            ids = get_ids(shapefile, idvariable)
        else:
            ids = None

        pts = ps.open(shapefile)

        # Get attributes if requested
        if attribute == True:
            dbname = os.path.splitext(shapefile)[0] + '.dbf'
            db = ps.open(dbname)
        else:
            db = None

        for i, pt in enumerate(pts):
            if ids and db:
                self.points[ids[i]] = {'coordinates':pt, 'properties':db[i]}
            elif ids and not db:
                self.points[ids[i]] = {'coordinates':pt, 'properties':None}
            elif not ids and db:
                self.points[i] = {'coordinates':pt, 'properties':db[i]}
            else:
                self.points[i] = {'coordinates':pt, 'properties':None}

        pts.close()
        if db:
            db.close()
        self.npoints = len(self.points.keys())

 def db(self, headerOnly=False):
     if self.data['config']['other_missingValueCheck']:
         pysal.MISSINGVALUE = self.data['config']['other_missingValue']
     if 'fname' in self.data:
         fileType = self.data['fname'].rsplit('.')[-1].lower()
         self.fileType = fileType
         if fileType == 'csv':
             if headerOnly:
                 f = pysal.open(self.data['fname'], 'rU')
                 db = {}
                 db['header'] = f.header
                 f.close()
                 return db
             else:
                 return pysal.open(self.data['fname'], 'rU')
         elif fileType == 'dbf':
             db = pysal.open(self.data['fname'], 'r')
             header = []
             # grab only the numeric fields.
             if headerOnly:
                 for field, spec in zip(db.header, db.field_spec):
                     typ = spec[0].lower()
                     if typ in ['d', 'n', 'f']:
                         header.append(field)
                 return {'header': header}
             else:
                 return db
         else:
             print "Unknown File Type"
             return False
     else:
         return None

def test2():
    #d6cd52286e5d3e9e08a5a42489180df3.shp
    import pysal
    shp = pysal.open('../www/tmp/d6cd52286e5d3e9e08a5a42489180df3.shp')
    dbf = pysal.open('../www/tmp/d6cd52286e5d3e9e08a5a42489180df3.dbf')        
    
    w = pysal.queen_from_shapefile('../www/tmp/d6cd52286e5d3e9e08a5a42489180df3.shp')
    #d3viz.moran_scatter_plot(shp, dbf, "dog_cnt", w)    
    import numpy as np
    y = np.array(dbf.by_col["dog_cnt"])
    lm = pysal.Moran_Local(y, w)    
    
    for i,j in enumerate(lm.q):
        if lm.p_sim[i] >= 0.05:
            print 0
        else:
            print j
        
    import d3viz
    d3viz.setup()
    
    d3viz.show_map(shp)
    d3viz.scatter_plot(shp, ["dog_cnt","home_cnt"])
    
    d3viz.lisa_map(shp, "dog_cnt", lm)

def map_line_shp(shp_link, which='all'):
    '''
    Create a map object from a line shapefile
    ...

    Arguments
    ---------

    shp_link        : str
                      Path to shapefile
    which           : str/list

    Returns
    -------

    map             : PatchCollection
                      Map object with the polygons from the shapefile

    '''
    shp = ps.open(shp_link)
    if which == 'all':
        db = ps.open(shp_link.replace('.shp', '.dbf'))
        n = len(db.by_col(db.header[0]))
        db.close()
        which = [True] * n
    patches = []
    for inwhich, shape in zip(which, shp):
        if inwhich:
            for xy in shape.parts:
                patches.append(xy)
    lc = LineCollection(patches)
    _ = _add_axes2col(lc, shp.bbox)
    return lc

def test():
    # Test
    shp = pysal.open(pysal.examples.get_path('NAT.shp'),'r')
    dbf = pysal.open(pysal.examples.get_path('NAT.dbf'),'r')
    
    show_map(shp)
    
    ids = get_selected(shp)
    print ids
    
    w = pysal.rook_from_shapefile(pysal.examples.get_path('NAT.shp'))
    moran_scatter_plot(shp, dbf, "HR90", w)
    
    scatter_plot(shp, ["HR90", "PS90"])
    scatter_plot_matrix(shp, ["HR90", "PS90"])
    
    quantile_map(shp, dbf, "HC60", 5, basemap="leaflet_map")
    
    
    select_ids = [i for i,v in enumerate(dbf.by_col["HC60"]) if v < 20.0]
    select(shp, ids=select_ids)
    
    
    quantile_map(shp, dbf, "HC60", 5)
    
    
    lisa_map(shp, dbf, "HC60", w)

 def test_sids(self):
     w = pysal.open(pysal.examples.get_path("sids2.gal")).read()
     f = pysal.open(pysal.examples.get_path("sids2.dbf"))
     SIDR = np.array(f.by_col("SIDR74"))
     mi = pysal.Moran(SIDR, w, two_tailed=False)
     np.testing.assert_allclose(mi.I, 0.24772519320480135, atol=ATOL, rtol=RTOL)
     self.assertAlmostEquals(mi.p_norm,  5.7916539074498452e-05)

def map_point_shp(shp_link, which='all'):
    '''
    Create a map object from a point shapefile
    ...

    Arguments
    ---------

    shp_link        : str
                      Path to shapefile
    which           : str/list

    Returns
    -------

    map             : PatchCollection
                      Map object with the points from the shapefile

    '''
    shp = ps.open(shp_link)
    if which == 'all':
        db = ps.open(shp_link.replace('.shp', '.dbf'))
        n = len(db.by_col(db.header[0]))
        db.close()
        which = [True] * n
    pts = []
    for inwhich, pt in zip(which, shp):
        if inwhich:
                pts.append(pt)
    pts = np.array(pts)
    sc = plt.scatter(pts[:, 0], pts[:, 1])
    _ = _add_axes2col(sc, shp.bbox)
    return sc

def getHexShp(codeBook,outFile):
    """Create the hexagon shapeFile
    Do NOT specify file extension in outFile"""
    inCOD = open(codeBook)
    inCOD = inCOD.readlines()
    header = inCOD.pop(0)
    dims,topo,cols,rows,neigh = header.split()
    cols=int(cols)
    rows=int(rows)
    fieldSpec = [('N',20,10)]*int(dims)
    fieldHead = ['plane%d'%i for i in range(int(dims))]
    #hexpts=hexPts(x,y) #Orig WRONG!!!
    hexpts=hexPts(cols,rows)
    polys=hexPoly(hexpts)# polygons
    out = pysal.open(outFile+'.shp','w')
    outd = pysal.open(outFile+'.dbf','w')
    outd.header = fieldHead
    outd.field_spec = fieldSpec
    for line in inCOD:
        outd.write(map(float,line.split()))
    outd.close()
    for bbox,poly in polys:
        poly = map(pysal.cg.shapes.Point,poly)
        poly = pysal.cg.shapes.Polygon(poly)
        out.write(poly)
    out.close()
    return 'Shapefile created'

def map_poly_shp(shp_link, which='all'):
    '''
    Create a map object from a shapefile
    ...

    Arguments
    ---------

    shp_link        : str
                      Path to shapefile
    which           : str/list

    Returns
    -------

    map             : PatchCollection
                      Map object with the polygons from the shapefile

    '''
    shp = ps.open(shp_link)
    if which == 'all':
        db = ps.open(shp_link.replace('.shp', '.dbf'))
        n = len(db.by_col(db.header[0]))
        db.close()
        which = [True] * n
    patches = []
    for inwhich, shape in zip(which, shp):
        if inwhich:
            for ring in shape.parts:
                xy = np.array(ring)
                patches.append(xy)
    return PolyCollection(patches)

 def setUp(self):
     f = pysal.open(pysal.examples.get_path("sids2.dbf"))
     varnames = ['SIDR74', 'SIDR79', 'NWR74', 'NWR79']
     self.names = varnames
     vars = [np.array(f.by_col[var]) for var in varnames]
     self.vars = vars
     self.w = pysal.open(pysal.examples.get_path("sids2.gal")).read()

def _writeShp():
    oShp = ps.open(out_shp, 'w')
    shp = ps.open(pt_shp)
    oDbf = ps.open(out_shp[:-3]+'dbf', 'w')
    dbf = ps.open(pt_shp[:-3]+'dbf')
    oDbf.header = dbf.header
    col_name = 'in_poly'
    col_spec = ('C', 14, 0)
    if polyID_col:
        col_name = polyID_col
        #db = ps.open(poly_shp[:-3]+'dbf')
        #col_spec = db.field_spec[db.header.index(polyID_col)]
    oDbf.header.append(col_name)
    oDbf.field_spec = dbf.field_spec
    oDbf.field_spec.append(col_spec)
    for poly, rec, i in zip(shp, dbf, correspondences):
        oShp.write(poly)
        rec.append(i)
        oDbf.write(rec)
    shp.close()
    oShp.close()
    dbf.close()
    oDbf.close()
    prj = open(pt_shp[:-3]+'prj').read()
    oPrj = open(out_shp[:-3]+'prj', 'w')
    oPrj.write(prj); oPrj.close()
    t4 = time.time()
    print '\t', t4-t3, ' seconds to write shapefile'
    print 'Shapefile written to %s'%out_shp