<jupyter_start><jupyter_text>---------------**If any part of this notebook is used in your research, please cite with the reference found in** **[README.md](https://github.com/pysal/spaghettibibtex-citation).**---------------- Spatial network analysis Demonstrating network representation and cluster detection**Author: <NAME>** ******This notebook is an advanced walk-through for:**1. Exploring the attributes of network objects and point patterns2. Understanding the difference between a network and its graph-theoretic representation3. Performing spatial network analysis<jupyter_code>%load_ext watermark
%watermark
import esda
import libpysal
import matplotlib
import matplotlib.pyplot as plt
import matplotlib_scalebar
from matplotlib_scalebar.scalebar import ScaleBar
import numpy
import spaghetti
import warnings
%matplotlib inline
%watermark -w
%watermark -iv
try:
    from IPython.display import set_matplotlib_formats
    set_matplotlib_formats("retina")
except ImportError:
    pass<jupyter_output><empty_output><jupyter_text>---------------- 1. Instantiating a `spaghetti.Network` object Instantiate the network from a `.shp` file<jupyter_code>ntw = spaghetti.Network(in_data=libpysal.examples.get_path("streets.shp"))<jupyter_output><empty_output><jupyter_text>------------------------------ 2. Allocating observations (snapping points) to a network:**A network is composed of a single topological representation of network elements (arcs and vertices) to which point patterns may be snapped.**<jupyter_code># Crimes
ntw.snapobservations(
    libpysal.examples.get_path("crimes.shp"), "crimes", attribute=True
)
# Schools
ntw.snapobservations(
    libpysal.examples.get_path("schools.shp"), "schools", attribute=False
)
ntw.pointpatterns<jupyter_output><empty_output><jupyter_text>**Attributes for every point pattern**1. ``dist_snapped`` dict keyed by point id with the value as snapped distance from observation to network arc<jupyter_code>ntw.pointpatterns["crimes"].dist_snapped[0]<jupyter_output><empty_output><jupyter_text>2. ``dist_to_vertex`` dict keyed by pointid with the value being a dict in the form         {node: distance to vertex, node: distance to vertex}<jupyter_code>ntw.pointpatterns["crimes"].dist_to_vertex[0]<jupyter_output><empty_output><jupyter_text>3. ``npoints`` point observations in set<jupyter_code>ntw.pointpatterns["crimes"].npoints<jupyter_output><empty_output><jupyter_text>4. ``obs_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), ... }<jupyter_code>ntw.pointpatterns["crimes"].obs_to_arc[(161, 162)]<jupyter_output><empty_output><jupyter_text>5. ``obs_to_vertex`` list of incident network vertices to snapped observation points<jupyter_code>ntw.pointpatterns["crimes"].obs_to_vertex[0]<jupyter_output><empty_output><jupyter_text>6. ``points`` geojson like representation of the point pattern.  Includes properties if read with attributes=``True``<jupyter_code>ntw.pointpatterns["crimes"].points[0]<jupyter_output><empty_output><jupyter_text>7. ``snapped_coordinates`` dict keyed by pointid with the value being (x-coord, y-coord)<jupyter_code>ntw.pointpatterns["crimes"].snapped_coordinates[0]<jupyter_output><empty_output><jupyter_text>**Counts per link (arc or edge) are important, but should not be precomputed since we have different representations of the network (spatial and graph currently).  (Relatively) Uniform segmentation still needs to be done.**<jupyter_code>counts = ntw.count_per_link(ntw.pointpatterns["crimes"].obs_to_arc, graph=False)
orig_counts_per_link = sum(list(counts.values())) / float(len(counts.keys()))
orig_counts_per_link<jupyter_output><empty_output><jupyter_text>----------------------------- 3. Network segmentation Split network arcs into segments of 200 ([US feet is this case](https://github.com/pysal/libpysal/blob/master/libpysal/examples/geodanet/streets.prj))<jupyter_code>n200 = ntw.split_arcs(200.0)<jupyter_output><empty_output><jupyter_text>Calculate observations per link on the segmented network<jupyter_code>counts = n200.count_per_link(n200.pointpatterns["crimes"].obs_to_arc, graph=False)
segm_counts_per_link = sum(counts.values()) / float(len(counts.keys()))
segm_counts_per_link<jupyter_output><empty_output><jupyter_text>--------------------- 4. Extract `geopandas.GeoDataFrame` objects of the vertices and arcs<jupyter_code># 'full' unsegmented network
vertices_df, arcs_df = spaghetti.element_as_gdf(
    ntw, vertices=ntw.vertex_coords, arcs=ntw.arcs
)
# network segmented at 200-meter increments
vertices200_df, arcs200_df = spaghetti.element_as_gdf(
    n200, vertices=n200.vertex_coords, arcs=n200.arcs
)<jupyter_output><empty_output><jupyter_text>---------------------- 5. Visualization demonstration* shapefile derived, unsegmented network with vertices in a larger, blue, semi-opaque form* distance segmented network with small, red, fully opaque vertices<jupyter_code>base = arcs_df.plot(color="k", alpha=0.25, figsize=(12, 12))
vertices_df.plot(ax=base, color="b", markersize=300, alpha=0.25)
vertices200_df.plot(ax=base, color="r", markersize=25, alpha=1.0)
# add scale bar
scalebar = ScaleBar(1, units="m", location="lower left")
base.add_artist(scalebar);<jupyter_output><empty_output><jupyter_text>--------------------------- 6. Spatial Analysis with [pysal/esda](https://esda.readthedocs.io/en/latest/) 6.1 Moran's *I* using the digitized network<jupyter_code># Compute the counts
counts = ntw.count_per_link(ntw.pointpatterns["crimes"].obs_to_arc, graph=False)

# Binary Adjacency
w = ntw.contiguityweights(graph=False)

# Build the y vector
arcs = w.neighbors.keys()
y = numpy.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))<jupyter_output>['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']<jupyter_text>6.2 Moran's *I* using the graph representation to generate the *W** Note that we have to regenerate the counts per arc, since the graph will have less edges.<jupyter_code># Compute the counts
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 = numpy.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))<jupyter_output>['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']<jupyter_text>6.3 Moran's *I* using the segmented network and intensities instead of counts<jupyter_code># Compute the counts
counts = n200.count_per_link(n200.pointpatterns["crimes"].obs_to_arc, graph=False)

# Binary Adjacency
w = n200.contiguityweights(graph=False)

# Build the y vector and convert from raw counts to intensities
arcs = w.neighbors.keys()
y = numpy.zeros(len(n200.arcs))

for i, a in enumerate(arcs):
    if a in counts.keys():
        length = n200.arc_lengths[a]
        y[i] = float(counts[a]) / float(length)

# Moran's I
res = esda.moran.Moran(y, w, permutations=99)

print(dir(res))<jupyter_output>['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']<jupyter_text>---------------------------- 7. Simulate a point pattern on the network* Need to supply a count of the number of points and a distirbution (default is uniform).  Generally, this will not be called by the user, since the simulation will be used for Monte Carlo permutation.<jupyter_code>npts = ntw.pointpatterns["crimes"].npoints
sim = ntw.simulate_observations(npts)
print(dir(sim))<jupyter_output>['__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__', 'dist_to_vertex', 'npoints', 'obs_to_arc', 'obs_to_vertex', 'points', 'snapped_coordinates']<jupyter_text>--------------------- 8. Network-constrained observation cluster testing *K*-function<jupyter_code>kres = ntw.NetworkK(ntw.pointpatterns["crimes"], permutations=99)
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()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16768440438111853, 'fraction_numerical': 0.011631288165164292, 'mean_word_length': 6.364739471805853, 'pattern_counts': {'":': 0, '<': 78, '<?xml version=': 0, '>': 78, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '330688', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3856, 'n_tokens_neox': 3669, 'n_words': 828, 'path': 'notebooks/network-analysis.ipynb', 'repo_name': 'alex-treebeard/spaghetti', 'size': '10328'}