Metadata-Version: 2.4
Name: nkde
Version: 0.1.0
Summary: Network kernel density estimation (NKDE / TNKDE) along street networks: simple, discontinuous and continuous kernels, network K-function and hot spots.
Author-email: Yusuf Eminoğlu <yusufeminoglu@gmail.com>
License-Expression: MIT
Project-URL: Homepage, https://github.com/YusufEminoglu/nkde_sdk
Project-URL: Repository, https://github.com/YusufEminoglu/nkde_sdk
Project-URL: Issues, https://github.com/YusufEminoglu/nkde_sdk/issues
Project-URL: Changelog, https://github.com/YusufEminoglu/nkde_sdk/blob/main/CHANGELOG.md
Keywords: kernel-density,nkde,tnkde,network-kde,spatial-analysis,street-network,road-safety,crash-hotspots,crime-analysis,retail,lixel,k-function,hotspot,gis
Classifier: Development Status :: 4 - Beta
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Operating System :: OS Independent
Classifier: Topic :: Scientific/Engineering :: GIS
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Developers
Classifier: Typing :: Typed
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.20.0
Provides-Extra: fast
Requires-Dist: scipy>=1.7.0; extra == "fast"
Provides-Extra: geo
Requires-Dist: scipy>=1.7.0; extra == "geo"
Requires-Dist: shapely>=2.0.0; extra == "geo"
Requires-Dist: geopandas>=0.12.0; extra == "geo"
Requires-Dist: networkx>=2.8.0; extra == "geo"
Provides-Extra: dev
Requires-Dist: pytest>=7.0.0; extra == "dev"
Requires-Dist: pytest-cov>=4.0.0; extra == "dev"
Requires-Dist: ruff>=0.1.0; extra == "dev"
Requires-Dist: mypy>=1.0.0; extra == "dev"
Dynamic: license-file

# 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

```bash
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

```bash
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:

| option | meaning |
|---|---|
| `--kernel quartic\|triangular\|epanechnikov\|gaussian\|uniform` | kernel (default quartic) |
| `--method simple\|discontinuous\|continuous` | estimator (default continuous) |
| `--adaptive`, `--trim-bw M` | Abramson adaptive bandwidth, capped at `M` (default 3 x bw) |
| `--weight-field F` | event weights, for example casualties |
| `--time-field F --times 0:365:7 --time-bw 30` | TNKDE (numbers or ISO dates; ISO dates become days) |
| `--max-snap M` | drop events farther than `M` metres from any road |
| `--split-crossings` | add junctions where roads cross without a shared vertex |
| `--max-depth N` | limit branch enumeration to N junctions |
| `--top N` | print the N highest-density lixels |

`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

```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:

```python
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:

```python
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`):

| kernel | `k(d)` |
|---|---|
| quartic | `15/16 (1 - u^2)^2 / h` |
| triangular | `(1 - |u|) / h` |
| Epanechnikov | `3/4 (1 - u^2) / h` |
| Gaussian | `exp(-u^2 / (2 (1/3)^2))`, truncated at `h` and renormalised (sigma = h/3) |
| uniform | `1 / (2h)` |

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.

| step | time | mass per event |
|---|---|---|
| snap 10,000 points | 1.3 s | |
| simple NKDE | 1.5 s | 3.57 |
| discontinuous NKDE | 2.6 s | 1.000000 |
| continuous NKDE | 5.9 s | 1.000000 |
| continuous, adaptive (pilot + final) | 12.2 s | |
| planar KDE at lixel mid points (scipy) | 1.1 s | |
| K-function, 1,000 events, 20 distances, 19 simulations | 3.6 s | |

## API overview

| object | purpose |
|---|---|
| `Network.from_lines`, `Network.from_edges`, `network_from_geojson` | build the network |
| `Network.lixelize(L)` -> `Lixels` | lixels (`edge`, `start`, `end`, `mid`, `length`, `midpoints_xy()`, `locate()`) |
| `Network.snap(xy, times=, weights=, max_distance=)` -> `Events` | snap points; `events_on_edges`, `random_events` |
| `nkde(network_or_lixels, events, bw, ...)` -> `NKDEResult` | NKDE / TNKDE (`intensity`, `density`, `bandwidths`, `integral()`, `to_geojson()`) |
| `event_kernel`, `abramson_bandwidths` | one event's kernel; adaptive bandwidths |
| `planar_kde(targets, events, bw)` | Euclidean KDE |
| `network_k_function`, `pairwise_network_distances` | K-function and envelopes |
| `rank_hotspots`, `monte_carlo_pvalues`, `hotspot_clusters` | hot spots |
| `read_lines`, `read_points`, `events_from_geojson`, `result_to_geojson` | GeoJSON I/O |
| `lonlat_to_utm`, `utm_to_lonlat`, `utm_transformers`, `UTMZone` | projection |
| `nkde.interop` | `network_from_geodataframe`, `events_from_geodataframe`, `result_to_geodataframe`, `network_from_networkx`, `network_to_networkx` |

## 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
