Paper Code pyscan library

Accepted as a full paper at ACM SIGSPATIAL 2026 — 58 of 248 submissions. Joint work with Drew McClelland (Taptap Send), Michael Matheny (Meta), and Jeff M. Phillips.

The problem. Public-health and crime maps almost never come with addresses — you get one total for a whole county or zip code, and some counties are the size of a small state. The fastest algorithms for spotting a hotspot need actual points on the map, so the standard shortcut squashes each county down to a single dot at its center. That dot erases the county’s real shape and size — and a disease outbreak that a health department finds late, or not at all, is measured in lives, not in decimals of statistical power.

How I solved it. Our fix is almost embarrassingly simple: instead of one dot at the center, scatter 20–50 points across the county’s actual shape and share its case count among them. The map gets its geography back, and the fast algorithms never notice the difference. The proof it works: given nothing but county totals of Valley Fever cases, the method redraws California’s San Joaquin Valley endemic region — a region fixed by soil ecology — far more faithfully than the centroid shortcut ever could, while scanning the entire continental-U.S. map — 3,711 region polygons spanning 3,108 counties — in 0.33 seconds, roughly 3,000× faster than FlexScan’s 1,109 s (scan step only; FlexScan additionally runs a significance test the other methods do not).

Spatial scan statistics are a core tool for anomaly detection in geospatial data — locating regions where a measured quantity (disease cases, crime, and so on) is significantly elevated relative to a baseline. The most efficient scan algorithms operate on point data, but real-world data is usually aggregated into predefined regions such as census tracts, zip codes, or counties. The standard workaround, used by widely adopted tools like SaTScan, collapses each region to its centroid — convenient, but it discards the region’s spatial extent and substantially reduces statistical power.

This work proposes a simple, scalable alternative: replace each region with 20–50 points sampled uniformly from its geometry, spreading the region’s baseline and measured values evenly across them (Geom k). It preserves the region’s spatial structure while staying fully compatible with fast point-based scan algorithms, and pairs with pyscan’s C++ backend and adaptive gridding so that even a 50× increase in points adds little runtime.

A convergence analysis explains why so few samples per region suffice, and — perhaps surprisingly — that as a map is divided into more regions, fewer sample points per region are needed; empirically the recovery error follows 1/√k. Across six datasets (NYC zip codes and the counties of Arkansas, Utah, California, Georgia, and the continental U.S.), the method recovers planted anomalies at far smaller effect sizes than the centroid baseline. On the continental U.S. the scan takes 0.33 s, against 114 s for an exact area-based method — which is additionally handed the target size a practitioner never has — and 1,109 s for FlexScan, whose figure also includes a Monte-Carlo significance test the others do not run. All timings are the scan step alone.

On a real public-health dataset — county-level Valley Fever incidence in California — it recovers the known San Joaquin Valley endemic region far more accurately than the centroid approach: Point Jaccard distance 0.448 against 0.717 — within 0.036 of the 0.412 floor that no axis-aligned rectangle can beat on that target.

We recommend this sampling-based conversion as the default way to apply point-based spatial scan statistics to region-aggregated data. The companion repository contains the Python experiments and figure-rendering scripts that reproduce every figure and runtime table in the paper, built on pyscan’s C++ backend.