pysal.spaghetti
¶import os
last_modified = None
if os.name == "posix":
last_modified = !stat -f\
"# This notebook was last updated: %Sm"\
Network_Usage.ipynb
elif os.name == "nt":
last_modified = !for %a in (Network_Usage.ipynb)\
do echo # This notebook was last updated: %~ta
if last_modified:
get_ipython().set_next_input(last_modified[-1])
# This notebook was last updated: Nov 12 12:29:49 2019
# pysal submodule imports
from libpysal import examples
import spaghetti as spgh
import esda
import numpy as np
import matplotlib.pyplot as plt
import time
%matplotlib inline
__author__ = "James Gaboardi <jgaboardi@gmail.com>"
ntw = spgh.Network(in_data=examples.get_path('streets.shp'))
# Crimes
ntw.snapobservations(
examples.get_path('crimes.shp'), 'crimes', attribute=True
)
# Schools
ntw.snapobservations(
examples.get_path('schools.shp'), 'schools', attribute=False
)
ntw.pointpatterns
{'crimes': <spaghetti.network.PointPattern at 0x11c758fd0>, 'schools': <spaghetti.network.PointPattern at 0x11c733c88>}
dist_snapped
dict keyed by pointid with the value as snapped distance from observation to network arcdist_to_vertex
dict keyed by pointid with the value being a dict in the form
{node: distance to vertex, node: distance to vertex}npoints
point observations in setobs_to_arc
dict keyed by arc with the value being a dict in the form
{pointID:(x-coord, y-coord), pointID:(x-coord, y-coord), ... }obs_to_vertex
list of incident network vertices to snapped observation pointspoints
geojson like representation of the point pattern. Includes properties if read with attributes=Truesnapped_coordinates
dict keyed by pointid with the value being (x-coord, y-coord)counts = ntw.count_per_link(
ntw.pointpatterns['crimes'].obs_to_arc, graph=False
)
sum(list(counts.values())) / float(len(counts.keys()))
2.682242990654206
n200 = ntw.split_arcs(200.0)
counts = n200.count_per_link(
n200.pointpatterns['crimes'].obs_to_arc, graph=False
)
sum(counts.values()) / float(len(counts.keys()))
2.0354609929078014
geopandas.GeoDataFrame
objects of the vertices and arcs¶# 'full' unsegmented network
vertices_df, arcs_df = spgh.element_as_gdf(
ntw, vertices=ntw.vertex_coords, arcs=ntw.arcs
)
# network segmented at 200-meter increments
vertices200_df, arcs200_df = spgh.element_as_gdf(
n200, vertices=n200.vertex_coords, arcs=n200.arcs
)
base = arcs_df.plot(color='k', alpha=.25, figsize=(12,12))
vertices_df.plot(ax=base, color='b', markersize=300, alpha=.25)
arcs200_df.plot(ax=base, color='k', alpha=.25)
vertices200_df.plot(ax=base, color='r', markersize=25, alpha=1.)
<matplotlib.axes._subplots.AxesSubplot at 0x11ca5bda0>
# Binary Adjacency
w = ntw.contiguityweights(graph=False)
# Build the y vector
arcs = w.neighbors.keys()
y = np.zeros(len(arcs))
for i, a in enumerate(arcs):
if a in counts.keys():
y[i] = counts[a]
# Moran's I
res = esda.moran.Moran(y, w, permutations=99)
print(dir(res))
['EI', 'EI_sim', 'I', 'VI_norm', 'VI_rand', 'VI_sim', '_Moran__calc', '_Moran__moments', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_statistic', 'by_col', 'n', 'p_norm', 'p_rand', 'p_sim', 'p_z_sim', 'permutations', 'seI_norm', 'seI_rand', 'seI_sim', 'sim', 'w', 'y', 'z', 'z2ss', 'z_norm', 'z_rand', 'z_sim']
counts = ntw.count_per_link(ntw.pointpatterns['crimes'].obs_to_arc, graph=True)
# Binary Adjacency
w = ntw.contiguityweights(graph=True)
# Build the y vector
edges = w.neighbors.keys()
y = np.zeros(len(edges))
for i, e in enumerate(edges):
if e in counts.keys():
y[i] = counts[e]
# Moran's I
res = esda.moran.Moran(y, w, permutations=99)
print(dir(res))
['EI', 'EI_sim', 'I', 'VI_norm', 'VI_rand', 'VI_sim', '_Moran__calc', '_Moran__moments', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_statistic', 'by_col', 'n', 'p_norm', 'p_rand', 'p_sim', 'p_z_sim', 'permutations', 'seI_norm', 'seI_rand', 'seI_sim', 'sim', 'w', 'y', 'z', 'z2ss', 'z_norm', 'z_rand', 'z_sim']
# Binary Adjacency
w = n200.contiguityweights(graph=False)
# Compute the counts
counts = n200.count_per_link(
n200.pointpatterns['crimes'].obs_to_arc, graph=False
)
# Build the y vector and convert from raw counts to intensities
arcs = w.neighbors.keys()
y = np.zeros(len(arcs))
for i, a in enumerate(edges):
if a in counts.keys():
length = n200.arc_lengths[a]
y[i] = counts[a] / length
# Moran's I
res = esda.moran.Moran(y, w, permutations=99)
print(dir(res))
['EI', 'EI_sim', 'I', 'VI_norm', 'VI_rand', 'VI_sim', '_Moran__calc', '_Moran__moments', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_statistic', 'by_col', 'n', 'p_norm', 'p_rand', 'p_sim', 'p_z_sim', 'permutations', 'seI_norm', 'seI_rand', 'seI_sim', 'sim', 'w', 'y', 'z', 'z2ss', 'z_norm', 'z_rand', 'z_sim']
/Users/jgaboardi/miniconda3/envs/py3_spgh_dev/lib/python3.6/site-packages/esda/moran.py:201: RuntimeWarning: invalid value encountered in double_scalars k = k_num / k_den /Users/jgaboardi/miniconda3/envs/py3_spgh_dev/lib/python3.6/site-packages/esda/moran.py:212: RuntimeWarning: invalid value encountered in double_scalars return self.n / self.w.s0 * inum / self.z2ss /Users/jgaboardi/miniconda3/envs/py3_spgh_dev/lib/python3.6/site-packages/esda/moran.py:160: RuntimeWarning: invalid value encountered in greater_equal above = sim >= self.I /Users/jgaboardi/miniconda3/envs/py3_spgh_dev/lib/python3.6/site-packages/esda/moran.py:176: RuntimeWarning: invalid value encountered in true_divide self.z /= sy
t1 = time.time()
n0 = ntw.allneighbordistances(ntw.pointpatterns['crimes'])
print(time.time()-t1)
0.5213339328765869
t1 = time.time()
n1 = n200.allneighbordistances(n200.pointpatterns['crimes'])
print(time.time()-t1)
1.7624461650848389
t1 = time.time()
n0 = ntw.allneighbordistances(ntw.pointpatterns['crimes'])
print(time.time()-t1)
0.17091989517211914
t1 = time.time()
n1 = n200.allneighbordistances(n200.pointpatterns['crimes'])
print(time.time()-t1)
0.19501996040344238
npts = ntw.pointpatterns['crimes'].npoints
sim = ntw.simulate_observations(npts)
sim
<spaghetti.network.SimulatedPointPattern at 0x11ee5bf98>
fres = ntw.NetworkF(ntw.pointpatterns['crimes'],
permutations=99)
plt.figure(figsize=(8,8))
plt.plot(fres.xaxis, fres.observed, 'b-', linewidth=1.5, label='Observed')
plt.plot(fres.xaxis, fres.upperenvelope, 'r--', label='Upper')
plt.plot(fres.xaxis, fres.lowerenvelope, 'k--', label='Lower')
plt.legend(loc='best', fontsize='x-large')
plt.title('Network F Function', fontsize='xx-large')
plt.show()
gres = ntw.NetworkG(ntw.pointpatterns['crimes'], permutations=99)
plt.figure(figsize=(8,8))
plt.plot(gres.xaxis, gres.observed, 'b-', linewidth=1.5, label='Observed')
plt.plot(gres.xaxis, gres.upperenvelope, 'r--', label='Upper')
plt.plot(gres.xaxis, gres.lowerenvelope, 'k--', label='Lower')
plt.legend(loc='best', fontsize='x-large')
plt.title('Network G Function', fontsize='xx-large')
plt.show()
kres = ntw.NetworkK(ntw.pointpatterns['crimes'], permutations=99)
/Users/jgaboardi/spaghetti/spaghetti/analysis.py:401: RuntimeWarning: invalid value encountered in less_equal y[i] = len(nearest[nearest <= r])
plt.figure(figsize=(8,8))
plt.plot(kres.xaxis, kres.observed, 'b-', linewidth=1.5, label='Observed')
plt.plot(kres.xaxis, kres.upperenvelope, 'r--', label='Upper')
plt.plot(kres.xaxis, kres.lowerenvelope, 'k--', label='Lower')
plt.legend(loc='best', fontsize='x-large')
plt.title('Network K Function', fontsize='xx-large')
plt.show()