This page was generated from notebooks/network-spatial-dependence.ipynb. Interactive online version: Binder badge

If any part of this notebook is used in your research, please cite with the reference found in

Network-constrained spatial dependence

Demonstrating cluster detection along networks with the Global Auto K function

Author: James D. Gaboardi

This notebook is an advanced walk-through for:

  1. Understanding the global auto K function with an elementary geometric object

  2. Basic examples with synthetic data

  3. Empirical examples

%config InlineBackend.figure_format = "retina"
%load_ext watermark
Last updated: 2022-04-15T11:24:32.839694-04:00

Python implementation: CPython
Python version       : 3.10.4
IPython version      : 8.2.0

Compiler    : Clang 12.0.1
OS          : Darwin
Release     : 21.4.0
Machine     : x86_64
Processor   : i386
CPU cores   : 8
Architecture: 64bit

import geopandas
import libpysal
import matplotlib
import matplotlib_scalebar
from matplotlib_scalebar.scalebar import ScaleBar
import numpy
import spaghetti

%matplotlib inline
%watermark -w
%watermark -iv
/Users/the-gaboardi/spaghetti/spaghetti/ FutureWarning: The next major release of pysal/spaghetti (2.0.0) will drop support for all ```` geometries. This change is a first step in refactoring ``spaghetti`` that is expected to result in dramatically reduced runtimes for network instantiation and operations. Users currently requiring network and point pattern input as ```` geometries should prepare for this simply by converting to ``shapely`` geometries.
  warnings.warn(f"{dep_msg}", FutureWarning)
Watermark: 2.3.0

matplotlib_scalebar: 0.8.0
spaghetti          : 1.6.5+44.g17b88a0.dirty
libpysal           : 4.6.2
geopandas          : 0.10.2
matplotlib         : 3.5.1
numpy              : 1.21.6
json               : 2.0.9

The K function considers all pairwise distances of nearest neighbors to determine the existence of clustering, or lack thereof, over a delineated range of distances. For further description see O’Sullivan and Unwin (2010) and Okabe and Sugihara (2012).

  • D. O’Sullivan and D. J. Unwin. Point Pattern Analysis, chapter 5, pages 121–156. John Wiley & Sons, Ltd, 2010. doi:10.1002/9780470549094.ch5.

  • Atsuyki Okabe and Kokichi Sugihara. Network K Function Methods, chapter 6, pages 119–136. John Wiley & Sons, Ltd, 2012. doi:10.1002/9781119967101.ch6.

1. A demonstration of clustering

Results plotting helper function

def plot_k(k, _arcs, df1, df2, obs, scale=True, wr=[1, 1.2], size=(14, 7)):
    """Plot a Global Auto K-function and spatial context."""
    def function_plot(f, ax):
        """Plot a Global Auto K-function."""
        ax.plot(k.xaxis, k.observed, "b-", linewidth=1.5, label="Observed")
        ax.plot(k.xaxis, k.upperenvelope, "r--", label="Upper")
        ax.plot(k.xaxis, k.lowerenvelope, "k--", label="Lower")
        ax.legend(loc="best", fontsize="x-large")
        title_text = "Global Auto $K$ Function: %s\n" % obs
        title_text += "%s steps, %s permutations," % (k.nsteps, k.permutations)
        title_text += " %s distribution" % k.distribution
        f.suptitle(title_text, fontsize=25, y=1.1)
        ax.set_xlabel("Distance $(r)$", fontsize="x-large")
        ax.set_ylabel("$K(r)$", fontsize="x-large")

    def spatial_plot(ax):
        """Plot spatial context."""
        base = _arcs.plot(ax=ax, color="k", alpha=0.25)
        df1.plot(ax=base, color="g", markersize=30, alpha=0.25)
        df2.plot(ax=base, color="g", marker="x", markersize=100, alpha=0.5)
        carto_elements(base, scale)

    sub_args = {"gridspec_kw":{"width_ratios": wr}, "figsize":size}
    fig, arr = matplotlib.pyplot.subplots(1, 2, **sub_args)
    function_plot(fig, arr[0])

def carto_elements(b, scale):
    """Add/adjust cartographic elements."""
    if scale:
        kw = {"units":"ft", "dimension":"imperial-length", "fixed_value":1000}
        b.add_artist(ScaleBar(1, **kw))
    b.set(xticklabels=[], xticks=[], yticklabels=[], yticks=[]);

Equilateral triangle

def equilateral_triangle(x1, y1, x2, mids=True):
    """Return an equilateral triangle and its side midpoints."""
    x3 = (x1+x2)/2.
    y3 = numpy.sqrt((x1-x2)**2 - (x3-x1)**2) + y1
    p1, p2, p3 = (x1, y1), (x2, y1), (x3, y3)
    eqitri =[p1, p2, p3, p1])
    if mids:
        eqvs = eqitri.vertices[:-1]
        eqimps, vcount = [], len(eqvs),
        for i in range(vcount):
            for j in range(i+1, vcount):
                (_x1, _y1), (_x2, _y2) = eqvs[i], eqvs[j]
                mp =, (_y1+_y2)/2.))
    return eqitri, eqimps
eqtri_sides, eqtri_midps = equilateral_triangle(0., 0., 6., 1)
ntw = spaghetti.Network(eqtri_sides)
ntw.snapobservations(eqtri_midps, "eqtri_midps")
vertices_df, arcs_df = spaghetti.element_as_gdf(
    ntw, vertices=ntw.vertex_coords, arcs=ntw.arcs
eqv = spaghetti.element_as_gdf(ntw, pp_name="eqtri_midps")
eqv_snapped = spaghetti.element_as_gdf(ntw, pp_name="eqtri_midps", snapped=True)
id geometry comp_label
0 0 POINT (3.00000 0.00000) 0
1 1 POINT (1.50000 2.59808) 0
2 2 POINT (4.50000 2.59808) 0
kres = ntw.GlobalAutoK(
plot_k(kres, arcs_df, eqv, eqv_snapped, "eqtri_mps", wr=[1, 1.8], scale=False)


  • This example demonstrates a complete lack of clustering with a strong indication of dispersion when approaching 5 units of distance.

2. Synthetic examples

Regular lattice — distinguishing visual clustering from statistical clustering

bounds = (0,0,3,3)
h, v = 2, 2
lattice = spaghetti.regular_lattice(bounds, h, nv=v, exterior=True)
ntw = spaghetti.Network(in_data=lattice)

Network arc midpoints: statistical clustering

midpoints = []
for chain in lattice:
    (v1x, v1y), (v2x, v2y) = chain.vertices
    mp =, (v1y+v2y)/2.))
ntw.snapobservations(midpoints, "midpoints")

All observations on two network arcs: visual clustering

npts = len(midpoints) * 2
xs = [0.0] * npts + [2.0] * npts
ys = list(numpy.linspace(0.4,0.6, npts)) + list(numpy.linspace(2.1,2.9, npts))
pclusters = [ for xy in zip(xs,ys)]
ntw.snapobservations(pclusters, "pclusters")
vertices_df, arcs_df = spaghetti.element_as_gdf(ntw, vertices=True, arcs=True)
midpoints = spaghetti.element_as_gdf(ntw, pp_name="midpoints", snapped=True)
pclusters = spaghetti.element_as_gdf(ntw, pp_name="pclusters", snapped=True)

Visual clustering

kres = ntw.GlobalAutoK(ntw.pointpatterns["pclusters"], nsteps=100, permutations=100)
plot_k(kres, arcs_df, pclusters, pclusters, "pclusters", wr=[1, 1.8], scale=False)


  • This example exhibits a high degree of clustering within 1 unit of distance followed by a complete lack of clustering, then a strong indication of clustering around 3.5 units of distance and above. Both colloquilly and statistically, this pattern is clustered.

Statistical clustering

kres = ntw.GlobalAutoK(ntw.pointpatterns["midpoints"], nsteps=100, permutations=100)
plot_k(kres, arcs_df, midpoints, midpoints, "midpoints", wr=[1, 1.8], scale=False)


  • This example exhibits no clustering within 1 unit of distance followed by large increases in clustering at each 1-unit increment. After 3 units of distance, this pattern is highly clustered. Statistically speaking, this pattern is clustered, but not colloquilly.

3. Empircal examples

Instantiate the network from a .shp file

ntw = spaghetti.Network(in_data=libpysal.examples.get_path("streets.shp"))
vertices_df, arcs_df = spaghetti.element_as_gdf(
    ntw, vertices=ntw.vertex_coords, arcs=ntw.arcs

Associate the network with point patterns

for pp_name in ["crimes", "schools"]:
    pp_shp = libpysal.examples.get_path("%s.shp" % pp_name)
    ntw.snapobservations(pp_shp, pp_name, attribute=True)
{'crimes': < at 0x10f748ac0>,
 'schools': < at 0x10f7496c0>}

Empircal — schools

schools = spaghetti.element_as_gdf(ntw, pp_name="schools")
schools_snapped = spaghetti.element_as_gdf(ntw, pp_name="schools", snapped=True)
kres = ntw.GlobalAutoK(
plot_k(kres, arcs_df, schools, schools_snapped, "schools")


  • Schools exhibit no clustering until roughly 1,000 feet then display more clustering up to approximately 3,000 feet, followed by high clustering up to 6,000 feet.

Empircal — crimes

crimes = spaghetti.element_as_gdf(ntw, pp_name="crimes")
crimes_snapped = spaghetti.element_as_gdf(ntw, pp_name="crimes", snapped=True)
kres = ntw.GlobalAutoK(
plot_k(kres, arcs_df, crimes, crimes_snapped, "crimes")


  • There is strong evidence to suggest crimes are clustered at all distances along the network, but more so within several feet of each other and at 10,000 feet.