nkde

Network kernel density estimation along street networks, in Python.

Crashes, crimes, shops and bus stops happen on streets. A planar (2D) kernel density spreads every event over a disc, so it puts density on the street behind a block that is 100 m away as the crow flies but 300 m away on foot or by car. It also counts every parallel street that the disc touches. Network kernel density estimation (NKDE) spreads each event along the network using shortest-path or branch distances, so the result is a density per metre of street.

nkde provides:

• the three NKDE estimators described by Okabe & Sugihara (2012) and implemented in the R package spNetwork (Gelb 2021): simple, discontinuous and continuous
• lixelisation (edges split into lixels of length L) and event snapping to the nearest edge, with offsets
• five kernels: quartic, triangular, Epanechnikov, Gaussian and uniform; a fixed or adaptive (Abramson) bandwidth; event weights
• TNKDE: a spatio-temporal density on a time grid
• planar KDE for comparison, the network K-function with Monte Carlo envelopes, and hot-spot ranking and clustering of lixels
• GeoJSON in and out, a dependency-free lon/lat to UTM helper, optional geopandas/networkx bridges, and a command line tool

The only required dependency is numpy. scipy (optional) speeds up the planar KDE.

Installation

pip install nkde              # numpy only
pip install "nkde[fast]"      # + scipy (planar KDE through a KD-tree)
pip install "nkde[geo]"       # + scipy, shapely, geopandas, networkx

Python 3.9 to 3.13.

Quick start

Command line

nkde run roads.geojson events.geojson --bw 300 --lixel 10 --method continuous -o density.geojson

Inputs are GeoJSON LineString/MultiLineString roads and Point/MultiPoint events in projected metres. With --lonlat, WGS84 inputs are projected to their UTM zone and the output is written back in lon/lat. Other options:

nkde kfunction roads.geojson events.geojson --tmax 1000 --step 50 --nsim 99 -o k.csv writes the network K-function with its simulation envelope.

The output is one LineString per lixel with properties lixel, edge, source_id, length, density (per metre, integrates to 1 over the network) and intensity (events per metre, integrates to the number of events or the total weight). TNKDE output has density_t<j> / intensity_t<j> columns, and the times are stored in the collection's nkde member.

Python

import nkde

net = nkde.network_from_geojson("roads.geojson")             # projected metres
events = nkde.events_from_geojson("crashes.geojson", net, max_distance=25)

res = nkde.nkde(net, events, bw=300, lixel_length=10, method="continuous")
res.to_geojson("density.geojson")

top = nkde.rank_hotspots(res, top=20)                         # highest-density lixels
clusters = nkde.hotspot_clusters(res.lixels, res.density >= 0.002)

Lon/lat data:

lines, _ = nkde.read_lines("roads_wgs84.geojson")
import numpy as np
allc = np.vstack(lines)
fwd, inv, zone = nkde.utm_transformers(allc[:, 0], allc[:, 1])  # zone.epsg -> 326xx
net = nkde.network_from_geojson("roads_wgs84.geojson", transform=fwd)
events = nkde.events_from_geojson("events_wgs84.geojson", net, transform=fwd)
res = nkde.nkde(net, events, 300, lixel_length=10)
res.to_geojson("density_wgs84.geojson", inverse_transform=inv)

Spatio-temporal density, comparison with a planar KDE, and K-function:

ev = net.snap(xy, times=days)                                 # times in any unit
tres = nkde.nkde(net, ev, 300, lixel_length=10, times=range(0, 365, 7), time_bw=30)
tres.density.shape                                            # (n_lixels, n_times)

planar = nkde.planar_kde(res.lixels, events, 300)             # per square metre

k = nkde.network_k_function(net, events, distances=range(50, 1001, 50), n_simulations=99)
k.clustered()                                                 # observed K above the envelope

examples/quickstart.py builds a small network and writes GeoJSON outputs. examples/benchmark.py reproduces the timings below.

Methods

Lixels and events

Each edge of length l is cut into ceil(l / L) equal lixels, so no lixel is longer than L. The density is evaluated at lixel mid points. Each event is snapped to the nearest point on the nearest edge and is stored as (edge, offset), where the offset is the distance in metres from the start of the edge. Junctions are the end points of lines and the vertices shared by several lines. split_at_intersections=True (--split-crossings) also adds junctions where lines cross without a shared vertex.

Kernels

All kernels have support [-h, h] and unit mass in 1D (u = d / h):

The planar KDE uses the 2D versions of these kernels, which have unit mass over the disc of radius h. Note that spNetwork's Gaussian kernel uses sigma = h without truncation, so Gaussian results differ from spNetwork by construction.

The three estimators

Take an event e and a point x on the network.

• simple: f(x) = k(d(e, x)), where d is the shortest-path distance. It is computed by a multi-source Dijkstra from the two ends of the event's edge, stopped at h. At a junction of degree n the full kernel continues on all n - 1 other edges, so one event has a mass greater than 1 near junctions (2.3 to 3.6 on a city grid, see below). Dead ends lose mass. This is the classical "NKDE with network distance". It can over-weight dense networks, but its ranking of lixels is often still useful.
• discontinuous (Okabe, Satoh & Sugihara 2009, equal split): the kernel spreads from the event in both directions. At a node of degree n each branch splits into n - 1 branches with weight w / (n - 1). Each branch stops when its distance reaches h. Mass is preserved exactly, except at dead ends. The density jumps at junctions.
• continuous (Sugihara, Satoh & Okabe 2010; Okabe & Sugihara 2012): at a node of degree n a branch continues on each other edge with weight 2w/n and is reflected back on its own edge with weight (2/n - 1) w. The density is continuous at nodes, unbiased, and mass is preserved exactly. A dead end (n = 1) reflects the whole branch, so mass is preserved even at dead ends.

As in spNetwork, the discontinuous and continuous kernels follow every walk shorter than h, not only shortest paths, so a branch can go around a block and come back. The number of walks grows exponentially with h / block length. nkde merges branches that enter the same directed edge at the same distance, because their future contributions are identical. On regular street grids this keeps the cost polynomial. max_depth (number of junctions crossed) and min_weight truncate the enumeration if needed. Both default to the exact estimator.

Implementation: each event produces a list of traversals (edge, entry end, start distance, weight) with a priority queue ordered by distance (the bounded Dijkstra for simple). The traversals of thousands of events are then expanded into (lixel, value) pairs in one vectorised numpy pass and summed with bincount. Only lixels within h are touched. Events are processed in chunks (chunk_size), which bounds memory use.

Adaptive bandwidth

With adaptive=True, a fixed-bandwidth pilot f~ is evaluated at the lixel of each event, and h_i = h0 * sqrt(g / f~(x_i)), where g is the geometric mean of the pilot values (Abramson 1982). h_i is capped at trim_bw (default 3 h0), as with spNetwork's trim_bw.

TNKDE

f(x, t) = (1 / W) sum_i w_i k_net(x; e_i) k_time(t - t_i), where both kernels have unit mass. It is evaluated at the lixels for each requested time. NKDEResult.integral() integrates over space and time.

K-function and hot spots

The global network K-function (Okabe & Yamada 2001) is K(t) = L / (n (n - 1)) * #{ordered pairs i != j : d(i, j) <= t}, where L is the network length. The envelope comes from n_simulations patterns of n uniform random points on the network (CSR on the network). Pointwise min/max or quantiles are available. rank_hotspots orders lixels by density and flags those above a quantile. monte_carlo_pvalues reruns the estimator on CSR patterns to give a p-value per lixel. hotspot_clusters labels connected groups of hot lixels.

Validation

The test suite (110 tests) checks:

• Unit mass of every kernel in 1D and 2D, by numerical integration.

• Hand-checkable cases. On a single line every estimator equals k(|x - e|) exactly. On a 3-arm star the value on another arm is k, k/2 and 2k/3 for the simple, discontinuous and continuous kernels. Between the event and the centre the continuous kernel is k(d) - k(d')/3 because of the reflected branch. At a dead end the continuous kernel adds the full mirror image.

• Agreement where the methods must coincide. With no junction within h, or with only degree-2 junctions, all three methods agree to 1e-15.

• Symmetry. The density at lixel b from an event at a equals the density at a from an event at b, for all three methods on a jittered grid. A centred event on a regular grid gives a mirror-symmetric density.

• Continuity. The continuous kernel has the same value on every edge next to a node. The discontinuous kernel does not.

• Mass conservation, integral over lixels with the midpoint rule and 1 m lixels. The table below uses 50 random events on a jittered 6 x 6 grid (no dead ends; h up to 4 block lengths, so branches loop around blocks):

  method	h = 60 m	h = 230 m	h = 420 m
  simple	1.000 to 1.98	1.29 to 3.20	1.89 to 5.04
  discontinuous	1 +/- 1e-6	1 +/- 1e-6	1 +/- 1e-6
  continuous	1 +/- 1e-6	1 +/- 1e-6	1 +/- 1e-6

  On a star with two dead ends within reach, the mass is 1.065 for simple, 0.895 for discontinuous and 1.000000 for continuous. Per kernel (continuous, h = 230 m) the mass error is 1e-8 for quartic, 8e-7 for triangular, 2e-6 for Epanechnikov, 4e-5 for Gaussian and 2e-3 for uniform. The last two errors come from the midpoint rule at the kernel's jump at h, not from the estimator.

• TNKDE integrates to the number of events over space and time. The K-function is checked against hand counts, and clustered patterns lie above the CSR envelope. UTM coordinates match PROJ within 1 mm, and snapping matches a brute-force search.

Cross-check against spNetwork. No numeric comparison with published spNetwork outputs was made, because no reference values were available that could be reproduced reliably. The estimators follow the definitions in Okabe & Sugihara (2012) and in the spNetwork documentation (Gelb 2021): equal split by n - 1, continuous factors 2/n and 2/n - 1, and enumeration of all walks within the bandwidth. The hand-checkable cases above test those definitions directly. The known difference is the Gaussian kernel convention (see Kernels). spNetwork may also place lixel sample points or handle the last, shorter lixel of each edge differently.

Performance

These timings come from python examples/benchmark.py on a laptop (AMD Ryzen 5 4600H, Python 3.12, numpy 2.5, single thread). The network is a jittered 80 x 80 street grid with 100 m blocks: 6,400 nodes, 12,640 edges and 1,272 km of street, cut into 133,462 lixels of 10 m. There are 10,000 random events, the bandwidth is 300 m and the kernel is quartic.

API overview

References

• Abramson, I. S. (1982). On bandwidth variation in kernel estimates: a square root law. The Annals of Statistics, 10(4), 1217-1223.
• Gelb, J. (2021). spNetwork: a package for network kernel density estimation. The R Journal, 13(2), 528-544.
• Okabe, A., & Sugihara, K. (2012). Spatial Analysis Along Networks: Statistical and Computational Methods. Wiley.
• Okabe, A., Satoh, T., & Sugihara, K. (2009). A kernel density estimation method for networks, its computational method and a GIS-based tool. International Journal of Geographical Information Science, 23(1), 7-32.
• Okabe, A., & Yamada, I. (2001). The K-function method on a network and its computational implementation. Geographical Analysis, 33(3), 271-290.
• Sugihara, K., Satoh, T., & Okabe, A. (2010). Simple and unbiased kernel function for network analysis. Proceedings of the 10th International Symposium on Communications and Information Technologies (ISCIT), 827-832.
• Xie, Z., & Yan, J. (2008). Kernel density estimation of traffic accidents in a network space. Computers, Environment and Urban Systems, 32(5), 396-406.

Citation

See CITATION.cff.

License

MIT (c) 2026 Yusuf Eminoğlu