In [28]:
import numpy as np
import pysal as ps
import random as rdm
from pysal.contrib.viz import mapping as maps
from matplotlib.collections import LineCollection
In [29]:
shp = ps.open(ps.examples.get_path("taz.shp"))
In [30]:
dbf = ps.open(ps.examples.get_path("taz.dbf"))
              
In [31]:
dbf.header
Out[31]:
['AREA',
 'PERIMETER',
 'CNTY',
 'RSA',
 'AIRDB',
 'TAZ2K',
 'SQ_MILE',
 'ACRE',
 'NEWSEQ',
 'CSA',
 'CSA_NEW',
 'SQMI_TAZ',
 'TAZ_NUM',
 'CountyFIPS']
In [32]:
fig = figure(figsize=(9,9))
base = maps.map_poly_shp(shp)
base.set_linewidth(0.75)
base.set_facecolor('none')
base.set_edgecolor('0.8')
ax = maps.setup_ax([base])
fig.add_axes(ax)
show()

County as unique values

In [33]:
cnty = np.array(dbf.by_col("CNTY"))
In [33]:
 
In [34]:
fig = figure(figsize=(9,9))
base = maps.map_poly_shp(shp)
base.set_linewidth(0.75)
base.set_facecolor('none')
base.set_edgecolor('0.8')
counties = maps.base_choropleth_unique(maps.map_poly_shp(shp), cnty)
counties.set_linewidth(0)
counties.set_alpha(.5)
ax = maps.setup_ax([base, counties])
fig.add_axes(ax)
show()
In [35]:
cents = np.array([poly.centroid for poly in shp])
cents[0]
Out[35]:
array([  601741.78690918,  3939798.30153461])
In [36]:
wrook = ps.rook_from_shapefile("taz.shp")
In [37]:
w = wrook
In [38]:
cents.min(axis=0)
Out[38]:
array([  282150.8269443 ,  3615409.10372805])
In [39]:
def w2line_graph(w, centroids):
    
    
    
    segments = []
    for i in w.id2i:
        origin = cents[i]
        for j in w.neighbors[i]:
            dest = cents[j]
            ij = [i,j]
            ij.sort()
            segments.append([origin, dest])
    #segs = LineCollection(segments)
    
    #maps._add_axes2col(segs, [minx, miny, maxx, maxy])
    return segments    

        
In [40]:
segs = w2line_graph(wrook, cents)
In [41]:
fig = figure(figsize=(9,9))
base = maps.map_poly_shp(shp)
base.set_linewidth(0.75)
base.set_facecolor('none')
base.set_edgecolor('0.8')
segs = LineCollection(segs)
maps._add_axes2col(segs, shp.bbox)
segs.set_linewidth(0.20)
ax = maps.setup_ax([base, segs])
fig.add_axes(ax)
show()

Intersection weights

In [42]:
wb = ps.regime_weights(np.array(dbf.by_col("CNTY")))
In [43]:
wb.n
Out[43]:
4109
In [44]:
wint = ps.weights.Wsets.w_intersection(wb, wrook)
In [45]:
segs = w2line_graph(wint, cents)
In [45]:
 
In [46]:
fig = figure(figsize=(9,9))
base = maps.map_poly_shp(shp)
base.set_linewidth(0.75)
base.set_facecolor('none')
base.set_edgecolor('0.8')
segs = LineCollection(segs)
maps._add_axes2col(segs, shp.bbox)
segs.set_linewidth(0.20)
ax = maps.setup_ax([base, segs])
fig.add_axes(ax)
show()
In [47]:
segments = w2line_graph(wint, cents)
In [48]:
fig = figure(figsize=(9,9))
base = maps.map_poly_shp(shp)
base.set_linewidth(0.75)
base.set_facecolor('none')
base.set_edgecolor('0.8')
counties = maps.base_choropleth_unique(maps.map_poly_shp(shp), cnty)
counties.set_linewidth(0)
counties.set_alpha(.5)
segs = LineCollection(segments)
maps._add_axes2col(segs, shp.bbox)
segs.set_linewidth(0.20)
segs.set_color('0.1')
ax = maps.setup_ax([base, counties, segs])
fig.add_axes(ax)
show()
In [ ]: