diff --git a/.gitignore b/.gitignore
index 73b2e17..7591d61 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,3 +12,5 @@ dist/
htmlcov/
.ruff_cache/
site/
+examples/cache/
+.cache/
diff --git a/CHANGELOG.md b/CHANGELOG.md
index cff9573..e400e6e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,14 @@ Revival release: the code base now targets Python 3.10+ with current
versions of pandas (>= 2.2), GeoPandas (>= 1.0), scikit-learn and statsmodels.
### Added
+- `crime_mapping.get_city_shape("Natal, RN, Brazil")` — thin wrapper around
+ `osmnx.geocode_to_gdf` returning the raw city shape.
+- `examples/natal.ipynb`: an executed end-to-end walkthrough on Natal with
+ synthetic data, also rendered in the documentation.
+- `PredictionPipeline.evaluate` accepts a list of scorings and returns a
+ DataFrame (one CV pass for all metrics).
+- README section explaining the framework (thesis, Chapter 3, Figures 7-12);
+ the README is now the documentation home page.
- `crime_mapping.load_study_area("City, Country")` fetches a study area
boundary from OpenStreetMap via `osmnx` (`pip install predspot[osm]`).
- `synthetic.generate_crimes` generates synthetic events inside any study
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index e578cda..06ae5fb 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -34,7 +34,15 @@ mkdocs serve # live preview at http://127.0.0.1:8000
mkdocs build --strict # what CI runs
```
-Pushing to `master` deploys the site to GitHub Pages automatically.
+Pushing to `master` deploys the site to GitHub Pages automatically. The home
+page is generated from `README.md` (see `docs/hooks/readme.py`), and the
+example notebook is rendered from `examples/natal.ipynb`; regenerate and
+re-execute it with:
+
+```bash
+python examples/build_natal_notebook.py
+jupyter nbconvert --to notebook --execute --inplace examples/natal.ipynb
+```
## Releasing
diff --git a/README.md b/README.md
index 39cae86..a341d08 100644
--- a/README.md
+++ b/README.md
@@ -1,33 +1,44 @@
# Predspot
[](https://github.com/adaj/predspot/actions/workflows/ci.yml)
+[](https://adaj.github.io/predspot/)
[](https://pypi.org/project/predspot/)
[](https://pypi.org/project/predspot/)
-[](LICENSE)
+[](https://github.com/adaj/predspot/blob/master/LICENSE)
## Overview 📖
-Predspot is a Python library for spatio-temporal crime prediction and hotspot detection. It combines machine learning techniques with spatial analysis to help predict and visualize crime patterns across time and space.
+Predspot is a Python library for spatio-temporal crime prediction and hotspot
+detection. It turns a table of georeferenced, timestamped crime events into a
+grid of *places* and a sequence of *periods*, builds time series features for
+every place and trains a scikit-learn model to forecast where the next period's
+hotspots will be.
Key features:
-- Spatial and temporal crime mapping
-- Feature engineering for time series data
-- Machine learning-based prediction pipeline
-- Crime hotspot detection using Kernel Density Estimation
-- Visualization tools for crime patterns
-## Status 🚧
+- Spatio-temporal crime mapping: kernel density estimation (KDE) on point grids,
+ or event counts on hexagonal and square grids, at daily, weekly or monthly resolution
+- Time series feature engineering: lagged autoregressive, difference, seasonal and
+ trend features (STL decomposition)
+- A prediction pipeline that accepts any scikit-learn regressor, with time series
+ cross-validation and recursive multi-step forecasts
+- Study areas fetched from OpenStreetMap by name
+- A synthetic crime generator (hotspots, trend, annual, weekly and hourly patterns)
+ to try everything without real data
-Predspot started as part of a master's thesis (2018-2019) and is being revived
-and modernised. The code base now targets Python 3.10+ with current versions of
-pandas (>= 2.2), GeoPandas (>= 1.0), scikit-learn and statsmodels, and is
-covered by a test suite. It remains research software: use it as a reference
-implementation and adapt it to your own data.
+Full documentation, with a quickstart, a user guide and the API reference, lives at
+**https://adaj.github.io/predspot/**. A complete, executed walkthrough for the city
+of Natal is in [`examples/natal.ipynb`](https://github.com/adaj/predspot/blob/master/examples/natal.ipynb).
## How to use? 🚀
-Full documentation, with a quickstart, a user guide and the API reference, lives at
-**https://adaj.github.io/predspot/**.
+Install from PyPI (Python 3.10 or newer):
+
+```bash
+pip install predspot # core
+pip install "predspot[osm]" # + study areas from OpenStreetMap (osmnx)
+pip install "predspot[contour]" # + GeoJSON contour export (geojsoncontour)
+```
Basic usage example:
@@ -35,15 +46,14 @@ Basic usage example:
from predspot import Dataset, PredictionPipeline
from predspot.crime_mapping import KDE, create_gridpoints
from predspot.feature_engineering import Seasonality, Trend, Diff
-
from predspot.utilities import PandasFeatureUnion
from sklearn.ensemble import RandomForestRegressor
-# Load and prepare data: crimes_df needs `tag`, `t`, `lon`, `lat` columns and
+# crimes_df needs `tag`, `t`, `lon`, `lat` columns;
# study_area_gdf is a GeoDataFrame with the boundary of the study area
dataset = Dataset(crimes_df, study_area_gdf)
-# Create prediction pipeline (monthly KDE on a 1 km point grid)
+# Monthly KDE on a 1 km point grid, lag features, random forest
pipeline = PredictionPipeline(
mapping=KDE(tfreq='M', grid=create_gridpoints(study_area_gdf, resolution=1)),
fextraction=PandasFeatureUnion([
@@ -76,12 +86,12 @@ with spatial hotspots and realistic temporal patterns (trend, annual cycle,
day-of-week and hour-of-day profiles):
```python
-from predspot import Dataset, load_study_area, generate_crimes
+from predspot import Dataset, get_city_shape, generate_crimes
-study_area = load_study_area("Natal, Rio Grande do Norte, Brazil")
-crimes = generate_crimes(study_area, n_events=5000, n_hotspots=4,
+city = get_city_shape("Natal, RN, Brazil")
+crimes = generate_crimes(city, n_events=5000, n_hotspots=4,
start="2019-01-01", end="2020-12-31", seed=0)
-dataset = Dataset(crimes, study_area)
+dataset = Dataset(crimes, city)
dataset.plot()
```
@@ -95,72 +105,186 @@ predictions, pipeline = run_prediction_pipeline(crimes, study_area, grid_resolut
print(pipeline.evaluate('r2', cv=3))
```
+### Input data format 📊
-## Development ⚡
+The crime data should be a pandas DataFrame with the following required columns:
-Predspot has five main modules:
+- `tag`: crime type
+- `t`: timestamp
+- `lon`: longitude (WGS84 degrees)
+- `lat`: latitude (WGS84 degrees)
+
+The study area should be a GeoDataFrame (with a CRS) defining the boundaries
+of interest.
+
+## The Predspot framework 🧭
+
+Predspot implements the framework described in Chapter 3 of the master's
+thesis [*Predspot: predicting crime hotspots with machine learning*](https://repositorio.ufrn.br/server/api/core/bitstreams/3655b8e1-2f32-4ce9-af9c-0e6b64d7af84/content)
+(Araújo Jr., 2019). The framework is split into two phases, mirroring the
+training and prediction steps of a machine learning system: **model
+selection**, where a model is trained, evaluated and saved, and **prediction
+service**, where that model is used operationally, period after period. The
+figures below are reproduced from the thesis.
+
+### Model selection
+
+**1. Dataset preparation** *(Figure 7)* — everything starts from three inputs:
+a crime database, the city shape and, optionally, auxiliary Points of Interest
+(PoI) from OpenStreetMap.
+
+
+
+- Crime records must carry at least latitude, longitude, timestamp and crime type.
+- The city shape acts as a spatial filter: events and PoI outside the boundary are dropped,
+ along with duplicates and "default" locations assigned to badly registered events.
+- Crimes are split into **crime scenarios** (one per crime type, and possibly
+ per day/night period) that are modelled separately, since aggregating different
+ crime types blurs their distinct spatial patterns.
+- In the package: `Dataset` (validation, WGS84 points) and `get_city_shape` /
+ `load_study_area` (city boundary from OpenStreetMap).
+
+**2. Feature ingest** *(Figures 8, 9 and 10)* — the step that assembles the
+feature matrix `X` and the target `y`. Before it starts, two units of analysis
+are chosen: the **spatial unit** (a grid derived from the crime mapping method)
+and the **temporal unit** (daily, weekly or monthly samples — finer units give
+sparser, harder-to-predict series).
+
+
+
+- **Spatio-temporal aggregation**: the mapping method (KDE on a point grid, or
+ counts on polygonal cells) turns the events of each period into one value per
+ grid cell, producing the time series `C_ij` of cell *i* at period *j*.
+- **Temporal feature extraction** *(Figure 9)*: each cell's series is decomposed
+ with STL into **trend** (`T`) and **seasonal** (`S`) components, and
+ **differenced** (`D`); the features are the *k* most recent lags of each
+ component, and the target is the series value one period ahead.
+
+
+
+- **Spatial aggregation of PoI** (optional): the density of each PoI category
+ around each cell (`G`) describes places geographically, complementing the
+ temporal features.
+- **Join**: temporal and geographic features are joined into one row per
+ `(place, time)` pair — the layout shown in Figure 10.
+
+
+
+- In the package: `create_gridpoints` / `create_gridhexagonal` /
+ `create_gridsquares` (grid), `KDE` / `QuadratCount` (spatio-temporal
+ aggregation), `Trend`, `Seasonality`, `Diff`, `AR` and `PandasFeatureUnion`
+ (temporal feature extraction and join). PoI features are not part of the
+ current package release.
+
+**3. Machine learning modelling** *(Figure 11)* — the feature set feeds a
+model-agnostic training loop.
+
+
+
+- **Feature selection** first: noisy lags and PoI layers are filtered by a
+ learning-based selector (an embedded, tree-based method in the thesis).
+- **Several supervised algorithms** are trained and **tuned**, rather than
+ betting on a single one; the thesis compared random forests and gradient boosting.
+- **Evaluation** uses time series K-fold cross-validation (train on the first
+ *k* folds, test on fold *k*+1), so no information from the future leaks into
+ training. The best model per crime scenario is the **golden model**, saved
+ together with its selected features.
+- In the package: `PredictionPipeline` (with `FeatureScaling`,
+ `FeatureSelection` and `Model` wrappers and `evaluate()` for the time series CV).
+
+### Prediction service
+
+**4. Prediction pipeline** *(Figure 12)* — the model selection steps are
+tailored for operation.
+
+
+
+- For each new period, only the most recent events are loaded (enough to
+ compute the *k* lags), filtered and split into scenarios as before.
+- The spatio-temporal aggregation and temporal feature extraction are repeated;
+ PoI features are reused, since they do not change over time.
+- The previously selected features are kept and the golden model returns the
+ crime incidence level `y_i` of every place *i* one period ahead — a
+ prediction layer ready to be mapped. The process repeats every period.
+- In the package: `PredictionPipeline.predict()` appends each forecast to the
+ series and advances one period, so repeated calls walk forward in time.
+
+**5. Web service** — the thesis also outlines how to wrap the pipeline in a
+decoupled service (with a file *volume* for data, models and predictions, and
+an ETL controller triggering the pipeline each period) that serves predictions
+as GeoJSON to existing GIS tools. That layer is outside the scope of this
+package, which covers the model selection phase and the prediction pipeline.
-`dataset_preparation`: Module for preparing and managing crime datasets and study areas.
+## Resources 📚
-`crime_mapping`: Module for spatial and temporal crime mapping: point, hexagonal and square grids, KDE-based density surfaces and per-cell counts (`QuadratCount`).
+- Master's thesis (full description of the framework, evaluation on Natal and
+ Boston, feature importance analysis): Araújo Jr., A. (2019).
+ [*Predspot: Predicting Crime Hotspots with Machine Learning*](https://repositorio.ufrn.br/server/api/core/bitstreams/3655b8e1-2f32-4ce9-af9c-0e6b64d7af84/content).
+ M.Sc. dissertation, PPgSC/UFRN, Natal, Brazil.
+- The earlier version of the framework: Araújo et al. (2018), *Towards a crime
+ hotspot detection framework for patrol planning* (HPCC/SmartCity/DSS).
+- Methods worth reading about: kernel density estimation for hotspot mapping
+ (Chainey, Tompson & Uhlig, 2008), STL time series decomposition (Cleveland
+ et al., 1990) and time series cross-validation (Bergmeir, Hyndman & Koo, 2018).
-`feature_engineering`: Module for time series feature engineering, including seasonality, trend, and difference features.
+## Cite us
-`ml_modelling`: Module that implements the prediction pipeline and model evaluation.
+If you use Predspot in your research, please cite us:
-`synthetic`: Module that generates synthetic crime events (hotspots + temporal patterns) inside any study area.
+APA:
+```
+Araújo Jr., A. (2019). Predspot: Predicting Crime Hotspots with Machine Learning. Master's dissertation, UFRN (Universidade Federal do Rio Grande do Norte), Natal, Brazil.
-### Installation steps 🛠️
+Araújo, A., Cacho, N., Bezerra, L., Vieira, C., & Borges, J. (2018, June). Towards a crime hotspot detection framework for patrol planning. In 2018 IEEE 20th International Conference on High Performance Computing and Communications; IEEE 16th International Conference on Smart City; IEEE 4th International Conference on Data Science and Systems (HPCC/SmartCity/DSS) (pp. 1256-1263). IEEE.
+```
-Predspot requires Python 3.10 or newer.
+or bibtex:
+```
+@mastersthesis{araujo2019predspot,
+ title={Predspot: Predicting crime hotspots with machine learning},
+ author={Araujo, Adelson},
+ year={2019},
+ school={Universidade Federal do Rio Grande do Norte},
+ url={https://repositorio.ufrn.br/server/api/core/bitstreams/3655b8e1-2f32-4ce9-af9c-0e6b64d7af84/content}
+}
-```bash
-pip install predspot # from PyPI
-pip install "predspot[osm]" # + study areas from OpenStreetMap (osmnx)
-pip install "predspot[contour]" # + GeoJSON contour export (geojsoncontour)
+@inproceedings{araujo2018towards,
+ title={Towards a crime hotspot detection framework for patrol planning},
+ author={Ara{\'u}jo, Adelson and Cacho, N{\'a}dia and Bezerra, Lucas and Vieira, Carlos and Borges, Jo{\~a}o},
+ booktitle={2018 IEEE 20th International Conference on High Performance Computing and Communications; IEEE 16th International Conference on Smart City; IEEE 4th International Conference on Data Science and Systems (HPCC/SmartCity/DSS)},
+ pages={1256--1263},
+ year={2018},
+ organization={IEEE}
+}
```
+## Development ⚡
+
+Predspot has five main modules:
+
+- `dataset_preparation`: preparing and managing crime datasets and study areas.
+- `crime_mapping`: spatial and temporal crime mapping — point, hexagonal and
+ square grids, KDE density surfaces, per-cell counts (`QuadratCount`) and
+ study areas from OpenStreetMap.
+- `feature_engineering`: time series feature engineering (seasonality, trend,
+ difference and autoregressive lags).
+- `ml_modelling`: the prediction pipeline and model evaluation.
+- `synthetic`: synthetic crime events (hotspots + temporal patterns) inside any study area.
+
From source, for development:
```bash
git clone https://github.com/adaj/predspot.git
cd predspot
pip install -e ".[dev,osm,contour]"
-```
-
-Core dependencies (installed automatically): pandas, geopandas, shapely,
-numpy, scipy, scikit-learn, statsmodels and matplotlib.
-
-### Tests 🧪
-
-```bash
ruff check src tests # lint
pytest # ~10 s
```
-See [CONTRIBUTING.md](CONTRIBUTING.md) for the release process and
-[CHANGELOG.md](CHANGELOG.md) for what changed between versions.
-
-### Input Data Format 📊
-
-The crime data should be a pandas DataFrame with the following required columns:
-- `tag`: Crime type
-- `t`: Timestamp
-- `lon`: Longitude
-- `lat`: Latitude
-
-The study area should be a GeoDataFrame defining the geographical boundaries of interest.
-
-## Resources 📚
-
-For more information on the methods used in Predspot, please search more about these methods:
-- Kernel Density Estimation for crime hotspot detection
-- Time series decomposition for feature engineering
-- Spatio-temporal crime prediction techniques
-
-## License 📜
-
-BSD 3-Clause
+See [CONTRIBUTING.md](https://github.com/adaj/predspot/blob/master/CONTRIBUTING.md)
+for the release process and
+[CHANGELOG.md](https://github.com/adaj/predspot/blob/master/CHANGELOG.md) for
+what changed between versions.
## Contributing 💡
@@ -173,35 +297,14 @@ Guidelines for contributing:
4. Push to the branch
5. Create a new Pull Request
+## License 📜
-## Cite us
-
-If you use Predspot in your research, please cite us:
-
-APA:
-```
-Araujo, A., & Cacho, N. (2019). Predspot: Predicting crime hotspots with machine learning. Master’s thesis, UFRN (Universidade Federal do Rio Grande do Norte), Natal, Brazil.
-
-Araújo, A., Cacho, N., Bezerra, L., Vieira, C., & Borges, J. (2018, June). Towards a crime hotspot detection framework for patrol planning. In 2018 IEEE 20th International Conference on High Performance Computing and Communications; IEEE 16th International Conference on Smart City; IEEE 4th International Conference on Data Science and Systems (HPCC/SmartCity/DSS) (pp. 1256-1263). IEEE.
-```
-
-or bibtex:
-```
-@article{araujo2019predspot,
- title={Predspot: Predicting crime hotspots with machine learning},
- author={Araujo, Adelson},
- year={2019},
- school={Universidade Federal do Rio Grande do Norte}
-}
-
-@inproceedings{araujo2018towards,
- title={Towards a crime hotspot detection framework for patrol planning},
- author={Ara{\'u}jo, Adelson and Cacho, N{\'a}dia and Bezerra, Lucas and Vieira, Carlos and Borges, Jo{\~a}o},
- booktitle={2018 IEEE 20th International Conference on High Performance Computing and Communications; IEEE 16th International Conference on Smart City; IEEE 4th International Conference on Data Science and Systems (HPCC/SmartCity/DSS)},
- pages={1256--1263},
- year={2018},
- organization={IEEE}
-}
-```
+BSD 3-Clause
+## Status 🚧
+Predspot started as part of a master's thesis (2018-2019) and was revived in
+2026. The code base now targets Python 3.10+ with current versions of pandas
+(>= 2.2), GeoPandas (>= 1.0), scikit-learn and statsmodels, and is covered by a
+test suite and continuous integration. It remains research software: use it as
+a reference implementation and adapt it to your own data.
diff --git a/docs/assets/thesis/fig07-dataset-preparation.png b/docs/assets/thesis/fig07-dataset-preparation.png
new file mode 100644
index 0000000..8e9b663
Binary files /dev/null and b/docs/assets/thesis/fig07-dataset-preparation.png differ
diff --git a/docs/assets/thesis/fig08-feature-ingest.png b/docs/assets/thesis/fig08-feature-ingest.png
new file mode 100644
index 0000000..f5e4e4f
Binary files /dev/null and b/docs/assets/thesis/fig08-feature-ingest.png differ
diff --git a/docs/assets/thesis/fig09-temporal-feature-extraction.png b/docs/assets/thesis/fig09-temporal-feature-extraction.png
new file mode 100644
index 0000000..0c12518
Binary files /dev/null and b/docs/assets/thesis/fig09-temporal-feature-extraction.png differ
diff --git a/docs/assets/thesis/fig10-feature-set-example.png b/docs/assets/thesis/fig10-feature-set-example.png
new file mode 100644
index 0000000..0cb40d4
Binary files /dev/null and b/docs/assets/thesis/fig10-feature-set-example.png differ
diff --git a/docs/assets/thesis/fig11-ml-modelling.png b/docs/assets/thesis/fig11-ml-modelling.png
new file mode 100644
index 0000000..f2c1f9c
Binary files /dev/null and b/docs/assets/thesis/fig11-ml-modelling.png differ
diff --git a/docs/assets/thesis/fig12-prediction-pipeline.png b/docs/assets/thesis/fig12-prediction-pipeline.png
new file mode 100644
index 0000000..d26ec40
Binary files /dev/null and b/docs/assets/thesis/fig12-prediction-pipeline.png differ
diff --git a/docs/examples/natal.ipynb b/docs/examples/natal.ipynb
new file mode 120000
index 0000000..a40a4d7
--- /dev/null
+++ b/docs/examples/natal.ipynb
@@ -0,0 +1 @@
+../../examples/natal.ipynb
\ No newline at end of file
diff --git a/docs/guide/modelling.md b/docs/guide/modelling.md
index 639be8f..c3fa9c3 100644
--- a/docs/guide/modelling.md
+++ b/docs/guide/modelling.md
@@ -32,7 +32,7 @@ estimator sees them as independent samples.
## Evaluating
```python
-pipeline.evaluate("r2", cv=5) # or "mse"
+pipeline.evaluate("r2", cv=5) # or "mse", or ["r2", "mse"] for a DataFrame
```
Periods are ordered and split with scikit-learn's `TimeSeriesSplit`, so each
diff --git a/docs/hooks/readme.py b/docs/hooks/readme.py
new file mode 100644
index 0000000..516e311
--- /dev/null
+++ b/docs/hooks/readme.py
@@ -0,0 +1,31 @@
+"""MkDocs hook: the documentation home page is the repository README.
+
+The README is written for GitHub (paths relative to the repository root); this
+hook rewrites those paths so the same content renders on the site.
+"""
+
+import re
+from pathlib import Path
+
+README = Path(__file__).resolve().parents[2] / "README.md"
+
+REPO_BLOB = "https://github.com/adaj/predspot/blob/master/"
+LINKS = {
+ REPO_BLOB + "CONTRIBUTING.md": "contributing.md",
+ REPO_BLOB + "CHANGELOG.md": "changelog.md",
+ REPO_BLOB + "examples/natal.ipynb": "examples/natal.ipynb",
+ "https://adaj.github.io/predspot/": "index.md",
+}
+
+
+def on_page_markdown(markdown, page, config, files):
+ if page.file.src_path != "index.md":
+ return markdown
+ text = README.read_text(encoding="utf-8")
+ text = text.replace('src="docs/assets/', 'src="assets/')
+ text = text.replace("](docs/assets/", "](assets/")
+ for url, target in LINKS.items():
+ text = text.replace(f"]({url})", f"]({target})")
+ # Keep the page's own front matter (title etc.), if any.
+ front = re.match(r"^---\n.*?\n---\n", markdown, flags=re.S)
+ return (front.group(0) if front else "") + text
diff --git a/docs/index.md b/docs/index.md
index d9ba291..c4f5189 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -1,60 +1,4 @@
-# Predspot
-
-**Predicting crime hotspots with machine learning.**
-
-Predspot is a Python library for spatio-temporal crime prediction. It turns a
-table of georeferenced, timestamped crime events into a grid of *places* and a
-sequence of *periods*, builds time series features for every place and trains a
-scikit-learn model to forecast where the next period's hotspots will be.
-
-
- { width="900" }
- Left: forecast density for the next month on a 500 m point grid.
- Right: observed monthly density and the forecast at the hottest grid point.
-
-
-## What it does
-
-- **Spatio-temporal mapping** — kernel density estimation (KDE) on a grid of
- points, or event counts on hexagonal / square grids, at daily, weekly or
- monthly resolution.
-- **Feature engineering** — lagged autoregressive, difference, seasonal and
- trend features (STL decomposition) for every place.
-- **Prediction pipeline** — any scikit-learn regressor (or `Pipeline` with
- scaling and feature selection) trained to forecast the next period, with
- time series cross-validation.
-- **Study areas from OpenStreetMap** — fetch a city boundary from a name.
-- **Synthetic data** — generate realistic events (hotspots, trend, annual,
- weekly and hourly patterns) inside any study area to try things out.
-
-## In a nutshell
-
-```python
-from predspot import Dataset, load_study_area, generate_crimes
-from predspot.pipeline import build_default_pipeline
-
-study_area = load_study_area("Natal, Rio Grande do Norte, Brazil") # (1)!
-crimes = generate_crimes(study_area, n_events=5000, seed=0) # (2)!
-
-dataset = Dataset(crimes, study_area)
-pipeline = build_default_pipeline(study_area, tfreq="M", grid_resolution=1)
-pipeline.fit(dataset)
-
-print(pipeline.evaluate("r2", cv=3)) # time series cross-validation
-forecast = pipeline.predict() # density per grid point, next month
-```
-
-1. Needs `pip install "predspot[osm]"`. Any GeoDataFrame with a boundary works too.
-2. Replace with your own DataFrame with `tag`, `t`, `lon`, `lat` columns.
-
-Head to the [installation](getting-started/installation.md) and
-[quickstart](getting-started/quickstart.md) pages, or read
-[how Predspot works](guide/concepts.md).
-
-## About
-
-Predspot was created by [Adelson Araujo](https://github.com/adaj) as part of
-his master's thesis at the Universidade Federal do Rio Grande do Norte (UFRN),
-Brazil, and revived in 2026 for current versions of Python and its scientific
-stack. It is research software released under the BSD-3-Clause license; see
-[Citing Predspot](citing.md) if you use it in your work.
+---
+title: Home
+---
+
diff --git a/examples/build_natal_notebook.py b/examples/build_natal_notebook.py
new file mode 100644
index 0000000..b6e8863
--- /dev/null
+++ b/examples/build_natal_notebook.py
@@ -0,0 +1,515 @@
+"""Build examples/natal.ipynb from source cells.
+
+The notebook itself is the deliverable; this script only exists so that the
+notebook can be regenerated (and re-executed) reproducibly:
+
+ python examples/build_natal_notebook.py
+ jupyter nbconvert --to notebook --execute --inplace examples/natal.ipynb
+"""
+
+from pathlib import Path
+
+import nbformat as nbf
+
+cells = []
+
+
+def md(text):
+ cells.append(nbf.v4.new_markdown_cell(text.strip()))
+
+
+def code(text):
+ cells.append(nbf.v4.new_code_cell(text.strip()))
+
+
+md("""
+# Predspot end to end: crime hotspots for Natal, Brazil
+
+This notebook walks through the whole Predspot workflow on the city of Natal
+(Rio Grande do Norte, Brazil), using **synthetic** crime events so that it
+runs anywhere without confidential police data:
+
+1. fetch the city boundary from OpenStreetMap;
+2. generate synthetic crime events with spatial hotspots and temporal patterns;
+3. build the `Dataset`, the spatial grids and the spatio-temporal series (KDE and counts);
+4. extract time series features;
+5. fit, evaluate and inspect a `PredictionPipeline`;
+6. forecast the next months and check how well the true hotspots are recovered.
+
+Every intermediate object is displayed so you can see exactly what flows
+between the steps. Requirements: `pip install "predspot[osm]" matplotlib`.
+""")
+
+code("""
+import warnings
+
+import geopandas as gpd
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+import predspot
+from predspot import Dataset, PredictionPipeline, PandasFeatureUnion, generate_crimes, get_city_shape
+from predspot.crime_mapping import KDE, QuadratCount, create_gridhexagonal, create_gridpoints
+from predspot.feature_engineering import AR, Diff, Seasonality, Trend
+
+warnings.filterwarnings("ignore", category=UserWarning)
+pd.set_option("display.width", 120)
+pd.set_option("display.max_columns", 20)
+plt.rcParams["figure.dpi"] = 100
+
+print("predspot", predspot.__version__, "| geopandas", gpd.__version__, "| pandas", pd.__version__)
+""")
+
+md("""
+## 1. Study area: the boundary of Natal
+
+`get_city_shape` geocodes a place name with [osmnx](https://osmnx.readthedocs.io)
+(Nominatim) and returns its boundary as a GeoDataFrame in WGS84. The query
+should be specific enough for the first match to be the municipality.
+""")
+
+code("""
+city = get_city_shape("Natal, RN, Brazil")
+city
+""")
+
+code("""
+polygon = city.geometry.iloc[0]
+print("geometry:", polygon.geom_type, "| vertices:", len(polygon.exterior.coords))
+print("bounds (W, S, E, N):", np.round(city.total_bounds, 4))
+area_km2 = city.to_crs(city.estimate_utm_crs()).area.iloc[0] / 1e6
+print(f"area: {area_km2:.1f} km²")
+
+ax = city.plot(figsize=(5, 6), color="#f2f2f2", edgecolor="black")
+ax.set_title("Natal, RN (OpenStreetMap boundary)")
+ax.set_axis_off()
+""")
+
+md("""
+## 2. Synthetic crime events
+
+`generate_crimes` draws events from a space-time point process inside the
+polygon: a mixture of Gaussian **hotspots** plus a uniform background in space,
+and a temporal intensity with a linear **trend**, an **annual** cycle and
+**day-of-week** / **hour-of-day** profiles. With `return_hotspots=True` we also
+get the hotspot centres, which lets us check later whether the model finds them.
+""")
+
+code("""
+crimes, hotspots = generate_crimes(
+ city,
+ n_events=12_000,
+ start="2019-01-01",
+ end="2021-12-31",
+ n_hotspots=5,
+ hotspot_share=0.65,
+ hotspot_sd_km=0.7,
+ trend=0.4, # +40% events from start to end
+ annual_amplitude=0.25, # busier around the peak month...
+ annual_peak_month=12, # ...December
+ seed=2019,
+ return_hotspots=True,
+)
+crimes.head(10)
+""")
+
+code("""
+print(crimes.shape)
+crimes.describe(include="all").T
+""")
+
+code("""
+hotspots
+""")
+
+code("""
+crimes["tag"].value_counts().to_frame("events").assign(share=lambda d: (d["events"] / len(crimes)).round(3))
+""")
+
+md("""
+### 2.1 Where and when do the events happen?
+""")
+
+code("""
+fig, ax = plt.subplots(figsize=(7, 8))
+city.plot(ax=ax, color="#f7f7f7", edgecolor="black")
+sample = crimes.sample(4000, random_state=0)
+ax.scatter(sample["lon"], sample["lat"], s=3, alpha=0.35, color="#1f4e79", label="events (sample)")
+hotspots.plot(ax=ax, color="#d62728", marker="*", markersize=150, zorder=5, label="hotspot centres")
+ax.legend(loc="lower left")
+ax.set_title("Synthetic crime events in Natal")
+ax.set_axis_off()
+""")
+
+code("""
+monthly = crimes.set_index("t").resample("ME").size().rename("events")
+fig, axes = plt.subplots(1, 3, figsize=(15, 3.6))
+monthly.plot(ax=axes[0], marker="o")
+axes[0].set_title("Events per month (trend + annual cycle)")
+axes[0].set_xlabel("")
+crimes["t"].dt.day_name().value_counts().reindex(
+ ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"]
+).plot.bar(ax=axes[1], color="#1f4e79")
+axes[1].set_title("Events per weekday")
+crimes["t"].dt.hour.value_counts().sort_index().plot.bar(ax=axes[2], color="#1f4e79", width=0.9)
+axes[2].set_title("Events per hour of day")
+fig.tight_layout()
+""")
+
+md("""
+## 3. Dataset, grids and the spatio-temporal series
+
+`Dataset` validates the events (columns `tag`, `t`, `lon`, `lat`) and turns
+them into a GeoDataFrame of points in WGS84, together with the study area.
+""")
+
+code("""
+dataset = Dataset(crimes, city)
+dataset
+""")
+
+code("""
+dataset.crimes.head()
+""")
+
+md("""
+### 3.1 Grids
+
+Predspot discretises the city into **places**. Two kinds of grid are built
+below: a grid of **points** every 500 m (used by `KDE`) and a grid of
+**hexagons** of 1 km² (used by `QuadratCount`). Only cells intersecting the
+city are kept; each grid has `lon`/`lat` columns with the cell centroid and an
+index named `places`.
+""")
+
+code("""
+points = create_gridpoints(city, resolution=0.5)
+hexes = create_gridhexagonal(city, resolution=1.0)
+print(f"{len(points)} grid points at 500 m | {len(hexes)} hexagons of 1 km²")
+points.head()
+""")
+
+code("""
+fig, axes = plt.subplots(1, 2, figsize=(12, 7))
+city.boundary.plot(ax=axes[0], color="black")
+points.plot(ax=axes[0], markersize=4, color="#1f4e79")
+axes[0].set_title(f"Point grid, 500 m ({len(points)} places)")
+city.boundary.plot(ax=axes[1], color="black")
+hexes.plot(ax=axes[1], facecolor="none", edgecolor="#1f4e79", linewidth=0.6)
+axes[1].set_title(f"Hexagonal grid, 1 km² ({len(hexes)} places)")
+for ax in axes:
+ ax.set_axis_off()
+""")
+
+md("""
+### 3.2 Spatio-temporal mapping
+
+A mapping assigns one value to every `(period, place)` pair. `KDE` fits a
+Gaussian kernel density estimate to the events of each period and evaluates it
+at the grid points; `QuadratCount` counts the events inside each polygon. Both
+return the same object: a `pandas.Series` named `crime_density` with a
+`(t, places)` MultiIndex — the **spatio-temporal series**.
+""")
+
+code("""
+kde = KDE(tfreq="M", grid=points, bandwidth="silverman")
+stseries = kde.fit_transform(dataset.crimes)
+print(type(stseries).__name__, stseries.shape, "| KDE factor:", round(kde.factor, 4))
+stseries.head(8)
+""")
+
+code("""
+# The same series in wide form: one row per month, one column per place
+wide = stseries.unstack("places")
+wide.iloc[:6, :8]
+""")
+
+code("""
+counts = QuadratCount(tfreq="M", grid=hexes).fit_transform(dataset.crimes)
+counts_wide = counts.unstack("places")
+print("events per month recovered by the counts:", counts_wide.sum(axis=1).astype(int).head(3).tolist(), "...")
+counts_wide.iloc[:6, :8]
+""")
+
+code("""
+month = pd.Timestamp("2021-06-30")
+fig, axes = plt.subplots(1, 2, figsize=(13, 6.5))
+
+g = points.copy()
+g["density"] = stseries.xs(month, level="t").reindex(points.index).values
+g.plot(ax=axes[0], column="density", cmap="viridis", markersize=12, marker="s", legend=True,
+ legend_kwds={"shrink": 0.6, "label": "KDE density"})
+city.boundary.plot(ax=axes[0], color="black", linewidth=1)
+hotspots.plot(ax=axes[0], color="#d62728", marker="*", markersize=150, zorder=5)
+axes[0].set_title(f"KDE on the point grid - {month:%B %Y}")
+
+h = hexes.copy()
+h["events"] = counts.xs(month, level="t").reindex(hexes.index).values
+h.plot(ax=axes[1], column="events", cmap="viridis", edgecolor="white", linewidth=0.3, legend=True,
+ legend_kwds={"shrink": 0.6, "label": "events in the cell"})
+city.boundary.plot(ax=axes[1], color="black", linewidth=1)
+hotspots.plot(ax=axes[1], color="#d62728", marker="*", markersize=150, zorder=5)
+axes[1].set_title(f"QuadratCount on hexagons - {month:%B %Y}")
+for ax in axes:
+ ax.set_axis_off()
+""")
+
+md("""
+### 3.3 The KDE bandwidth
+
+`bandwidth="silverman"` estimates the kernel width from the spread of the
+events in the first period and keeps it fixed, which makes densities comparable
+over time but tends to over-smooth a whole city. A numeric `bandwidth` is used
+directly as the KDE factor (a multiple of the city-wide spread of the events):
+smaller values give sharper maps. Compare three settings for the same month —
+this is the trade-off discussed in the thesis (Chapter 2, Figure 3).
+""")
+
+code("""
+fig, axes = plt.subplots(1, 3, figsize=(16, 5.5))
+for ax, bw in zip(axes, ["silverman", 0.2, 0.08]):
+ m = KDE(tfreq="M", grid=points, bandwidth=bw)
+ st = m.fit_transform(dataset.crimes)
+ g = points.copy()
+ g["density"] = st.xs(month, level="t").reindex(points.index).values
+ g.plot(ax=ax, column="density", cmap="viridis", markersize=10, marker="s")
+ city.boundary.plot(ax=ax, color="black", linewidth=1)
+ hotspots.plot(ax=ax, color="#d62728", marker="*", markersize=120, zorder=5)
+ ax.set_title(f"bandwidth={bw!r} (factor {m.factor:.3f})")
+ ax.set_axis_off()
+""")
+
+md("""
+## 4. Time series features
+
+From here on we use `bandwidth=0.2`, a middle ground between the
+over-smoothed Silverman estimate and a noisy small kernel.
+
+Each place has a monthly series. The feature classes transform it (STL trend,
+STL seasonal component, first difference or the raw series) and build `lags`
+lagged columns. `PandasFeatureUnion` aligns them on the `(t, places)` index and
+drops the warm-up rows. Note the extra row for the month **after** the last
+observed one — that is the row the model will forecast.
+""")
+
+code("""
+BANDWIDTH = 0.2
+stseries = KDE(tfreq="M", grid=points, bandwidth=BANDWIDTH).fit_transform(dataset.crimes)
+
+LAGS = 6
+fextraction = PandasFeatureUnion([
+ ("ar", AR(lags=3)),
+ ("seasonal", Seasonality(lags=LAGS)),
+ ("trend", Trend(lags=LAGS)),
+ ("diff", Diff(lags=LAGS)),
+])
+X = fextraction.fit_transform(stseries)
+print(X.shape, "| periods:", X.index.get_level_values("t").min().date(), "->", X.index.get_level_values("t").max().date())
+X.head(8)
+""")
+
+code("""
+# What the decomposition looks like for the busiest place
+busiest = stseries.groupby("places").mean().idxmax()
+ts = stseries.xs(busiest, level="places")
+from statsmodels.tsa.seasonal import STL
+res = STL(ts, period=LAGS).fit()
+fig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True)
+ts.plot(ax=axes[0], title=f"place {busiest}: KDE density per month")
+res.trend.plot(ax=axes[1], title="STL trend")
+res.seasonal.plot(ax=axes[2], title=f"STL seasonal (period={LAGS})")
+ts.diff().plot(ax=axes[3], title="first difference")
+for ax in axes:
+ ax.set_xlabel("")
+fig.tight_layout()
+""")
+
+md("""
+## 5. Fit and evaluate a prediction pipeline
+
+The pipeline chains the mapping, the feature extraction and a scikit-learn
+estimator. Here the estimator is itself a `Pipeline`: quantile scaling,
+recursive feature elimination and a random forest, each wrapped so that
+DataFrames (and the `(t, places)` index) survive every step.
+""")
+
+code("""
+from sklearn.ensemble import RandomForestRegressor
+from sklearn.feature_selection import RFE
+from sklearn.pipeline import Pipeline
+from sklearn.preprocessing import QuantileTransformer
+
+from predspot.feature_engineering import FeatureScaling
+from predspot.ml_modelling import FeatureSelection, Model
+
+pipeline = PredictionPipeline(
+ mapping=KDE(tfreq="M", grid=points, bandwidth=BANDWIDTH),
+ fextraction=PandasFeatureUnion([
+ ("ar", AR(lags=3)),
+ ("seasonal", Seasonality(lags=LAGS)),
+ ("trend", Trend(lags=LAGS)),
+ ("diff", Diff(lags=LAGS)),
+ ]),
+ estimator=Pipeline([
+ ("scaling", FeatureScaling(QuantileTransformer(n_quantiles=20, output_distribution="uniform"))),
+ ("selection", FeatureSelection(RFE(RandomForestRegressor(n_estimators=20, random_state=0, n_jobs=-1),
+ n_features_to_select=12, step=3))),
+ ("model", Model(RandomForestRegressor(n_estimators=100, min_samples_leaf=2, random_state=0, n_jobs=-1))),
+ ]),
+ random_state=0,
+)
+pipeline.fit(dataset)
+print("training rows:", pipeline.features.loc[: stseries.index.get_level_values("t").max()].shape[0])
+print("next period to forecast:", pipeline.next_time.date())
+""")
+
+code("""
+scores = pipeline.evaluate(["r2", "mse"], cv=3) # time series CV: train on earlier folds, test on the next
+scores.loc["mean"] = scores.mean()
+scores
+""")
+
+code("""
+fi = pipeline.feature_importances
+ax = fi.sort_values("importance").plot.barh(figsize=(7, 5), legend=False, color="#1f4e79")
+ax.set_title("Selected features and their importance")
+fi.head(12)
+""")
+
+md("""
+## 6. Forecast the next months
+
+`predict()` returns the density of every place for the period after the last
+observed one. Each call appends its forecast to the series, recomputes the
+features and moves the horizon one period forward, so calling it three times
+gives a three-month recursive forecast.
+""")
+
+code("""
+forecasts = pd.concat([pipeline.predict() for _ in range(3)])
+forecasts.groupby("t").describe().round(2)
+""")
+
+code("""
+first_month = forecasts.index.get_level_values("t").min()
+fc = points.join(forecasts.xs(first_month, level="t"))
+top10 = fc.nlargest(10, "crime_density")[["lon", "lat", "crime_density"]].round(4)
+print(f"Top-10 places for {first_month:%B %Y}:")
+top10
+""")
+
+code("""
+fig, ax = plt.subplots(figsize=(7, 8))
+fc.plot(ax=ax, column="crime_density", cmap="magma", markersize=14, marker="s", legend=True,
+ legend_kwds={"shrink": 0.6, "label": "forecast density"})
+city.boundary.plot(ax=ax, color="black", linewidth=1)
+hotspots.plot(ax=ax, color="#00e5ff", marker="*", markersize=170, zorder=5, label="true hotspot centres")
+top = fc.nlargest(20, "crime_density")
+ax.scatter(top["lon"], top["lat"], s=110, facecolors="none", edgecolors="#00e5ff", linewidths=1.6, label="top-20 forecast")
+ax.legend(loc="lower left")
+ax.set_title(f"Forecast hotspots for {first_month:%B %Y}")
+ax.set_axis_off()
+""")
+
+md("""
+### 6.1 Does the forecast find the true hotspots?
+
+Since the data are synthetic we know where the hotspots are. For each true
+hotspot centre we look up the nearest grid point and check how it ranks in the
+forecast (1 = hottest place of the city). Hotspots holding a larger share of the
+events should rank near the top; small hotspots on a 500 m grid compete with the
+many grid points that surround the biggest one.
+""")
+
+code("""
+def distance_km(lon1, lat1, lon2, lat2):
+ dx = (lon1 - lon2) * 111.32 * np.cos(np.radians((lat1 + lat2) / 2))
+ dy = (lat1 - lat2) * 110.57
+ return np.hypot(dx, dy)
+
+ranked = fc.sort_values("crime_density", ascending=False)
+ranked["rank"] = np.arange(1, len(ranked) + 1)
+rows = []
+for i, h in hotspots.iterrows():
+ d = distance_km(ranked["lon"].values, ranked["lat"].values, h.geometry.x, h.geometry.y)
+ nearest = ranked.iloc[int(d.argmin())]
+ rows.append({
+ "hotspot": i,
+ "share_of_events": round(h["share"], 3),
+ "nearest_place": int(nearest.name),
+ "distance_km": round(float(d.min()), 2),
+ "forecast_rank": int(nearest["rank"]),
+ "percentile": round(100 * (1 - nearest["rank"] / len(ranked)), 1),
+ })
+pd.DataFrame(rows).set_index("hotspot").sort_values("forecast_rank")
+""")
+
+code("""
+# The full history + forecast of the top forecast place
+place = fc["crime_density"].idxmax()
+series = pipeline.stseries.xs(place, level="places")
+observed = series.loc[: stseries.index.get_level_values("t").max()]
+predicted = series.loc[forecasts.index.get_level_values("t").min():]
+ax = observed.plot(figsize=(10, 3.5), marker="o", label="observed (KDE)")
+predicted.plot(ax=ax, marker="*", markersize=12, linestyle="--", color="#d62728", label="forecast")
+ax.set_title(f"Place {place}: monthly density and 3-month recursive forecast")
+ax.set_xlabel("")
+ax.legend()
+""")
+
+md("""
+## 7. The same pipeline with counts on hexagons
+
+`QuadratCount` is a drop-in replacement for `KDE`: swap the mapping and the
+grid, keep everything else.
+""")
+
+code("""
+hex_pipeline = PredictionPipeline(
+ mapping=QuadratCount(tfreq="M", grid=hexes),
+ fextraction=PandasFeatureUnion([
+ ("ar", AR(lags=3)),
+ ("seasonal", Seasonality(lags=LAGS)),
+ ("trend", Trend(lags=LAGS)),
+ ]),
+ estimator=RandomForestRegressor(n_estimators=100, min_samples_leaf=2, random_state=0, n_jobs=-1),
+ random_state=0,
+).fit(dataset)
+print("r2 per fold:", np.round(hex_pipeline.evaluate("r2", cv=3), 3))
+hex_fc = hexes.join(hex_pipeline.predict().droplevel("t"))
+
+fig, ax = plt.subplots(figsize=(7, 8))
+hex_fc.plot(ax=ax, column="crime_density", cmap="magma", edgecolor="white", linewidth=0.3, legend=True,
+ legend_kwds={"shrink": 0.6, "label": "forecast events"})
+city.boundary.plot(ax=ax, color="black", linewidth=1)
+hotspots.plot(ax=ax, color="#00e5ff", marker="*", markersize=170, zorder=5)
+ax.set_title(f"Forecast events per hexagon for {hex_pipeline.stseries.index.get_level_values('t').max():%B %Y}")
+ax.set_axis_off()
+""")
+
+md("""
+## Wrapping up
+
+- `get_city_shape` gave us the study area; `generate_crimes` produced events with known hotspots.
+- `Dataset` → grid → `KDE` / `QuadratCount` turned the events into a spatio-temporal series.
+- Lagged STL trend/seasonal, difference and autoregressive features fed a scikit-learn pipeline
+ wrapped by `PredictionPipeline`, evaluated with time series cross-validation.
+- `predict()` forecast the following months; the grid points nearest to the true hotspot
+ centres rank at the top of the forecast, in proportion to each hotspot's share of events.
+
+To run this on real data, replace the synthetic `crimes` DataFrame with your own
+`tag, t, lon, lat` table and model each crime type (`tag`) separately, as the
+framework recommends. See the [documentation](https://adaj.github.io/predspot/)
+for the details of every step.
+""")
+
+nb = nbf.v4.new_notebook()
+nb["cells"] = cells
+nb["metadata"] = {
+ "kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"},
+ "language_info": {"name": "python"},
+}
+out = Path(__file__).with_name("natal.ipynb")
+nbf.write(nb, out)
+print(f"wrote {out} with {len(cells)} cells")
diff --git a/examples/natal.ipynb b/examples/natal.ipynb
new file mode 100644
index 0000000..e995e61
--- /dev/null
+++ b/examples/natal.ipynb
@@ -0,0 +1,728 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "7e4e38dc",
+ "metadata": {},
+ "source": [
+ "# Predspot end to end: crime hotspots for Natal, Brazil\n",
+ "\n",
+ "This notebook walks through the whole Predspot workflow on the city of Natal\n",
+ "(Rio Grande do Norte, Brazil), using **synthetic** crime events so that it\n",
+ "runs anywhere without confidential police data:\n",
+ "\n",
+ "1. fetch the city boundary from OpenStreetMap;\n",
+ "2. generate synthetic crime events with spatial hotspots and temporal patterns;\n",
+ "3. build the `Dataset`, the spatial grids and the spatio-temporal series (KDE and counts);\n",
+ "4. extract time series features;\n",
+ "5. fit, evaluate and inspect a `PredictionPipeline`;\n",
+ "6. forecast the next months and check how well the true hotspots are recovered.\n",
+ "\n",
+ "Every intermediate object is displayed so you can see exactly what flows\n",
+ "between the steps. Requirements: `pip install \"predspot[osm]\" matplotlib`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1dac71a9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import warnings\n",
+ "\n",
+ "import geopandas as gpd\n",
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "\n",
+ "import predspot\n",
+ "from predspot import Dataset, PredictionPipeline, PandasFeatureUnion, generate_crimes, get_city_shape\n",
+ "from predspot.crime_mapping import KDE, QuadratCount, create_gridhexagonal, create_gridpoints\n",
+ "from predspot.feature_engineering import AR, Diff, Seasonality, Trend\n",
+ "\n",
+ "warnings.filterwarnings(\"ignore\", category=UserWarning)\n",
+ "pd.set_option(\"display.width\", 120)\n",
+ "pd.set_option(\"display.max_columns\", 20)\n",
+ "plt.rcParams[\"figure.dpi\"] = 100\n",
+ "\n",
+ "print(\"predspot\", predspot.__version__, \"| geopandas\", gpd.__version__, \"| pandas\", pd.__version__)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "cb9de6ba",
+ "metadata": {},
+ "source": [
+ "## 1. Study area: the boundary of Natal\n",
+ "\n",
+ "`get_city_shape` geocodes a place name with [osmnx](https://osmnx.readthedocs.io)\n",
+ "(Nominatim) and returns its boundary as a GeoDataFrame in WGS84. The query\n",
+ "should be specific enough for the first match to be the municipality."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b0fbc889",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "city = get_city_shape(\"Natal, RN, Brazil\")\n",
+ "city"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "88ba713c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "polygon = city.geometry.iloc[0]\n",
+ "print(\"geometry:\", polygon.geom_type, \"| vertices:\", len(polygon.exterior.coords))\n",
+ "print(\"bounds (W, S, E, N):\", np.round(city.total_bounds, 4))\n",
+ "area_km2 = city.to_crs(city.estimate_utm_crs()).area.iloc[0] / 1e6\n",
+ "print(f\"area: {area_km2:.1f} km²\")\n",
+ "\n",
+ "ax = city.plot(figsize=(5, 6), color=\"#f2f2f2\", edgecolor=\"black\")\n",
+ "ax.set_title(\"Natal, RN (OpenStreetMap boundary)\")\n",
+ "ax.set_axis_off()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e1f4ca42",
+ "metadata": {},
+ "source": [
+ "## 2. Synthetic crime events\n",
+ "\n",
+ "`generate_crimes` draws events from a space-time point process inside the\n",
+ "polygon: a mixture of Gaussian **hotspots** plus a uniform background in space,\n",
+ "and a temporal intensity with a linear **trend**, an **annual** cycle and\n",
+ "**day-of-week** / **hour-of-day** profiles. With `return_hotspots=True` we also\n",
+ "get the hotspot centres, which lets us check later whether the model finds them."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "caf16e7f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "crimes, hotspots = generate_crimes(\n",
+ " city,\n",
+ " n_events=12_000,\n",
+ " start=\"2019-01-01\",\n",
+ " end=\"2021-12-31\",\n",
+ " n_hotspots=5,\n",
+ " hotspot_share=0.65,\n",
+ " hotspot_sd_km=0.7,\n",
+ " trend=0.4, # +40% events from start to end\n",
+ " annual_amplitude=0.25, # busier around the peak month...\n",
+ " annual_peak_month=12, # ...December\n",
+ " seed=2019,\n",
+ " return_hotspots=True,\n",
+ ")\n",
+ "crimes.head(10)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "256277d3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(crimes.shape)\n",
+ "crimes.describe(include=\"all\").T"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "eb00b456",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "hotspots"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "999d61a3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "crimes[\"tag\"].value_counts().to_frame(\"events\").assign(share=lambda d: (d[\"events\"] / len(crimes)).round(3))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "bd6b0177",
+ "metadata": {},
+ "source": [
+ "### 2.1 Where and when do the events happen?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4f4153c8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig, ax = plt.subplots(figsize=(7, 8))\n",
+ "city.plot(ax=ax, color=\"#f7f7f7\", edgecolor=\"black\")\n",
+ "sample = crimes.sample(4000, random_state=0)\n",
+ "ax.scatter(sample[\"lon\"], sample[\"lat\"], s=3, alpha=0.35, color=\"#1f4e79\", label=\"events (sample)\")\n",
+ "hotspots.plot(ax=ax, color=\"#d62728\", marker=\"*\", markersize=150, zorder=5, label=\"hotspot centres\")\n",
+ "ax.legend(loc=\"lower left\")\n",
+ "ax.set_title(\"Synthetic crime events in Natal\")\n",
+ "ax.set_axis_off()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b15e0701",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "monthly = crimes.set_index(\"t\").resample(\"ME\").size().rename(\"events\")\n",
+ "fig, axes = plt.subplots(1, 3, figsize=(15, 3.6))\n",
+ "monthly.plot(ax=axes[0], marker=\"o\")\n",
+ "axes[0].set_title(\"Events per month (trend + annual cycle)\")\n",
+ "axes[0].set_xlabel(\"\")\n",
+ "crimes[\"t\"].dt.day_name().value_counts().reindex(\n",
+ " [\"Monday\", \"Tuesday\", \"Wednesday\", \"Thursday\", \"Friday\", \"Saturday\", \"Sunday\"]\n",
+ ").plot.bar(ax=axes[1], color=\"#1f4e79\")\n",
+ "axes[1].set_title(\"Events per weekday\")\n",
+ "crimes[\"t\"].dt.hour.value_counts().sort_index().plot.bar(ax=axes[2], color=\"#1f4e79\", width=0.9)\n",
+ "axes[2].set_title(\"Events per hour of day\")\n",
+ "fig.tight_layout()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b0340fd2",
+ "metadata": {},
+ "source": [
+ "## 3. Dataset, grids and the spatio-temporal series\n",
+ "\n",
+ "`Dataset` validates the events (columns `tag`, `t`, `lon`, `lat`) and turns\n",
+ "them into a GeoDataFrame of points in WGS84, together with the study area."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8a810174",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "dataset = Dataset(crimes, city)\n",
+ "dataset"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d915daaf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "dataset.crimes.head()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5552a527",
+ "metadata": {},
+ "source": [
+ "### 3.1 Grids\n",
+ "\n",
+ "Predspot discretises the city into **places**. Two kinds of grid are built\n",
+ "below: a grid of **points** every 500 m (used by `KDE`) and a grid of\n",
+ "**hexagons** of 1 km² (used by `QuadratCount`). Only cells intersecting the\n",
+ "city are kept; each grid has `lon`/`lat` columns with the cell centroid and an\n",
+ "index named `places`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "534bcc4d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "points = create_gridpoints(city, resolution=0.5)\n",
+ "hexes = create_gridhexagonal(city, resolution=1.0)\n",
+ "print(f\"{len(points)} grid points at 500 m | {len(hexes)} hexagons of 1 km²\")\n",
+ "points.head()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "77bc64fc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig, axes = plt.subplots(1, 2, figsize=(12, 7))\n",
+ "city.boundary.plot(ax=axes[0], color=\"black\")\n",
+ "points.plot(ax=axes[0], markersize=4, color=\"#1f4e79\")\n",
+ "axes[0].set_title(f\"Point grid, 500 m ({len(points)} places)\")\n",
+ "city.boundary.plot(ax=axes[1], color=\"black\")\n",
+ "hexes.plot(ax=axes[1], facecolor=\"none\", edgecolor=\"#1f4e79\", linewidth=0.6)\n",
+ "axes[1].set_title(f\"Hexagonal grid, 1 km² ({len(hexes)} places)\")\n",
+ "for ax in axes:\n",
+ " ax.set_axis_off()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e96fee0a",
+ "metadata": {},
+ "source": [
+ "### 3.2 Spatio-temporal mapping\n",
+ "\n",
+ "A mapping assigns one value to every `(period, place)` pair. `KDE` fits a\n",
+ "Gaussian kernel density estimate to the events of each period and evaluates it\n",
+ "at the grid points; `QuadratCount` counts the events inside each polygon. Both\n",
+ "return the same object: a `pandas.Series` named `crime_density` with a\n",
+ "`(t, places)` MultiIndex — the **spatio-temporal series**."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "110f6eb3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "kde = KDE(tfreq=\"M\", grid=points, bandwidth=\"silverman\")\n",
+ "stseries = kde.fit_transform(dataset.crimes)\n",
+ "print(type(stseries).__name__, stseries.shape, \"| KDE factor:\", round(kde.factor, 4))\n",
+ "stseries.head(8)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d9cf26da",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# The same series in wide form: one row per month, one column per place\n",
+ "wide = stseries.unstack(\"places\")\n",
+ "wide.iloc[:6, :8]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b92ac677",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "counts = QuadratCount(tfreq=\"M\", grid=hexes).fit_transform(dataset.crimes)\n",
+ "counts_wide = counts.unstack(\"places\")\n",
+ "print(\"events per month recovered by the counts:\", counts_wide.sum(axis=1).astype(int).head(3).tolist(), \"...\")\n",
+ "counts_wide.iloc[:6, :8]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ef9896a8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "month = pd.Timestamp(\"2021-06-30\")\n",
+ "fig, axes = plt.subplots(1, 2, figsize=(13, 6.5))\n",
+ "\n",
+ "g = points.copy()\n",
+ "g[\"density\"] = stseries.xs(month, level=\"t\").reindex(points.index).values\n",
+ "g.plot(ax=axes[0], column=\"density\", cmap=\"viridis\", markersize=12, marker=\"s\", legend=True,\n",
+ " legend_kwds={\"shrink\": 0.6, \"label\": \"KDE density\"})\n",
+ "city.boundary.plot(ax=axes[0], color=\"black\", linewidth=1)\n",
+ "hotspots.plot(ax=axes[0], color=\"#d62728\", marker=\"*\", markersize=150, zorder=5)\n",
+ "axes[0].set_title(f\"KDE on the point grid - {month:%B %Y}\")\n",
+ "\n",
+ "h = hexes.copy()\n",
+ "h[\"events\"] = counts.xs(month, level=\"t\").reindex(hexes.index).values\n",
+ "h.plot(ax=axes[1], column=\"events\", cmap=\"viridis\", edgecolor=\"white\", linewidth=0.3, legend=True,\n",
+ " legend_kwds={\"shrink\": 0.6, \"label\": \"events in the cell\"})\n",
+ "city.boundary.plot(ax=axes[1], color=\"black\", linewidth=1)\n",
+ "hotspots.plot(ax=axes[1], color=\"#d62728\", marker=\"*\", markersize=150, zorder=5)\n",
+ "axes[1].set_title(f\"QuadratCount on hexagons - {month:%B %Y}\")\n",
+ "for ax in axes:\n",
+ " ax.set_axis_off()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c0a2c16f",
+ "metadata": {},
+ "source": [
+ "### 3.3 The KDE bandwidth\n",
+ "\n",
+ "`bandwidth=\"silverman\"` estimates the kernel width from the spread of the\n",
+ "events in the first period and keeps it fixed, which makes densities comparable\n",
+ "over time but tends to over-smooth a whole city. A numeric `bandwidth` is used\n",
+ "directly as the KDE factor (a multiple of the city-wide spread of the events):\n",
+ "smaller values give sharper maps. Compare three settings for the same month —\n",
+ "this is the trade-off discussed in the thesis (Chapter 2, Figure 3)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "859b64f4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig, axes = plt.subplots(1, 3, figsize=(16, 5.5))\n",
+ "for ax, bw in zip(axes, [\"silverman\", 0.2, 0.08]):\n",
+ " m = KDE(tfreq=\"M\", grid=points, bandwidth=bw)\n",
+ " st = m.fit_transform(dataset.crimes)\n",
+ " g = points.copy()\n",
+ " g[\"density\"] = st.xs(month, level=\"t\").reindex(points.index).values\n",
+ " g.plot(ax=ax, column=\"density\", cmap=\"viridis\", markersize=10, marker=\"s\")\n",
+ " city.boundary.plot(ax=ax, color=\"black\", linewidth=1)\n",
+ " hotspots.plot(ax=ax, color=\"#d62728\", marker=\"*\", markersize=120, zorder=5)\n",
+ " ax.set_title(f\"bandwidth={bw!r} (factor {m.factor:.3f})\")\n",
+ " ax.set_axis_off()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "47e84097",
+ "metadata": {},
+ "source": [
+ "## 4. Time series features\n",
+ "\n",
+ "From here on we use `bandwidth=0.2`, a middle ground between the\n",
+ "over-smoothed Silverman estimate and a noisy small kernel.\n",
+ "\n",
+ "Each place has a monthly series. The feature classes transform it (STL trend,\n",
+ "STL seasonal component, first difference or the raw series) and build `lags`\n",
+ "lagged columns. `PandasFeatureUnion` aligns them on the `(t, places)` index and\n",
+ "drops the warm-up rows. Note the extra row for the month **after** the last\n",
+ "observed one — that is the row the model will forecast."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1d901dc3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "BANDWIDTH = 0.2\n",
+ "stseries = KDE(tfreq=\"M\", grid=points, bandwidth=BANDWIDTH).fit_transform(dataset.crimes)\n",
+ "\n",
+ "LAGS = 6\n",
+ "fextraction = PandasFeatureUnion([\n",
+ " (\"ar\", AR(lags=3)),\n",
+ " (\"seasonal\", Seasonality(lags=LAGS)),\n",
+ " (\"trend\", Trend(lags=LAGS)),\n",
+ " (\"diff\", Diff(lags=LAGS)),\n",
+ "])\n",
+ "X = fextraction.fit_transform(stseries)\n",
+ "print(X.shape, \"| periods:\", X.index.get_level_values(\"t\").min().date(), \"->\", X.index.get_level_values(\"t\").max().date())\n",
+ "X.head(8)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "66dd5428",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# What the decomposition looks like for the busiest place\n",
+ "busiest = stseries.groupby(\"places\").mean().idxmax()\n",
+ "ts = stseries.xs(busiest, level=\"places\")\n",
+ "from statsmodels.tsa.seasonal import STL\n",
+ "res = STL(ts, period=LAGS).fit()\n",
+ "fig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True)\n",
+ "ts.plot(ax=axes[0], title=f\"place {busiest}: KDE density per month\")\n",
+ "res.trend.plot(ax=axes[1], title=\"STL trend\")\n",
+ "res.seasonal.plot(ax=axes[2], title=f\"STL seasonal (period={LAGS})\")\n",
+ "ts.diff().plot(ax=axes[3], title=\"first difference\")\n",
+ "for ax in axes:\n",
+ " ax.set_xlabel(\"\")\n",
+ "fig.tight_layout()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0d3cca51",
+ "metadata": {},
+ "source": [
+ "## 5. Fit and evaluate a prediction pipeline\n",
+ "\n",
+ "The pipeline chains the mapping, the feature extraction and a scikit-learn\n",
+ "estimator. Here the estimator is itself a `Pipeline`: quantile scaling,\n",
+ "recursive feature elimination and a random forest, each wrapped so that\n",
+ "DataFrames (and the `(t, places)` index) survive every step."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9a613809",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sklearn.ensemble import RandomForestRegressor\n",
+ "from sklearn.feature_selection import RFE\n",
+ "from sklearn.pipeline import Pipeline\n",
+ "from sklearn.preprocessing import QuantileTransformer\n",
+ "\n",
+ "from predspot.feature_engineering import FeatureScaling\n",
+ "from predspot.ml_modelling import FeatureSelection, Model\n",
+ "\n",
+ "pipeline = PredictionPipeline(\n",
+ " mapping=KDE(tfreq=\"M\", grid=points, bandwidth=BANDWIDTH),\n",
+ " fextraction=PandasFeatureUnion([\n",
+ " (\"ar\", AR(lags=3)),\n",
+ " (\"seasonal\", Seasonality(lags=LAGS)),\n",
+ " (\"trend\", Trend(lags=LAGS)),\n",
+ " (\"diff\", Diff(lags=LAGS)),\n",
+ " ]),\n",
+ " estimator=Pipeline([\n",
+ " (\"scaling\", FeatureScaling(QuantileTransformer(n_quantiles=20, output_distribution=\"uniform\"))),\n",
+ " (\"selection\", FeatureSelection(RFE(RandomForestRegressor(n_estimators=20, random_state=0, n_jobs=-1),\n",
+ " n_features_to_select=12, step=3))),\n",
+ " (\"model\", Model(RandomForestRegressor(n_estimators=100, min_samples_leaf=2, random_state=0, n_jobs=-1))),\n",
+ " ]),\n",
+ " random_state=0,\n",
+ ")\n",
+ "pipeline.fit(dataset)\n",
+ "print(\"training rows:\", pipeline.features.loc[: stseries.index.get_level_values(\"t\").max()].shape[0])\n",
+ "print(\"next period to forecast:\", pipeline.next_time.date())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "eee77d05",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "scores = pipeline.evaluate([\"r2\", \"mse\"], cv=3) # time series CV: train on earlier folds, test on the next\n",
+ "scores.loc[\"mean\"] = scores.mean()\n",
+ "scores"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8f797ad0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fi = pipeline.feature_importances\n",
+ "ax = fi.sort_values(\"importance\").plot.barh(figsize=(7, 5), legend=False, color=\"#1f4e79\")\n",
+ "ax.set_title(\"Selected features and their importance\")\n",
+ "fi.head(12)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "47c7fc5c",
+ "metadata": {},
+ "source": [
+ "## 6. Forecast the next months\n",
+ "\n",
+ "`predict()` returns the density of every place for the period after the last\n",
+ "observed one. Each call appends its forecast to the series, recomputes the\n",
+ "features and moves the horizon one period forward, so calling it three times\n",
+ "gives a three-month recursive forecast."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "18edefe0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "forecasts = pd.concat([pipeline.predict() for _ in range(3)])\n",
+ "forecasts.groupby(\"t\").describe().round(2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "af03b229",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "first_month = forecasts.index.get_level_values(\"t\").min()\n",
+ "fc = points.join(forecasts.xs(first_month, level=\"t\"))\n",
+ "top10 = fc.nlargest(10, \"crime_density\")[[\"lon\", \"lat\", \"crime_density\"]].round(4)\n",
+ "print(f\"Top-10 places for {first_month:%B %Y}:\")\n",
+ "top10"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2bc5e223",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig, ax = plt.subplots(figsize=(7, 8))\n",
+ "fc.plot(ax=ax, column=\"crime_density\", cmap=\"magma\", markersize=14, marker=\"s\", legend=True,\n",
+ " legend_kwds={\"shrink\": 0.6, \"label\": \"forecast density\"})\n",
+ "city.boundary.plot(ax=ax, color=\"black\", linewidth=1)\n",
+ "hotspots.plot(ax=ax, color=\"#00e5ff\", marker=\"*\", markersize=170, zorder=5, label=\"true hotspot centres\")\n",
+ "top = fc.nlargest(20, \"crime_density\")\n",
+ "ax.scatter(top[\"lon\"], top[\"lat\"], s=110, facecolors=\"none\", edgecolors=\"#00e5ff\", linewidths=1.6, label=\"top-20 forecast\")\n",
+ "ax.legend(loc=\"lower left\")\n",
+ "ax.set_title(f\"Forecast hotspots for {first_month:%B %Y}\")\n",
+ "ax.set_axis_off()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "822a11ea",
+ "metadata": {},
+ "source": [
+ "### 6.1 Does the forecast find the true hotspots?\n",
+ "\n",
+ "Since the data are synthetic we know where the hotspots are. For each true\n",
+ "hotspot centre we look up the nearest grid point and check how it ranks in the\n",
+ "forecast (1 = hottest place of the city). Hotspots holding a larger share of the\n",
+ "events should rank near the top; small hotspots on a 500 m grid compete with the\n",
+ "many grid points that surround the biggest one."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cb82e3a9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def distance_km(lon1, lat1, lon2, lat2):\n",
+ " dx = (lon1 - lon2) * 111.32 * np.cos(np.radians((lat1 + lat2) / 2))\n",
+ " dy = (lat1 - lat2) * 110.57\n",
+ " return np.hypot(dx, dy)\n",
+ "\n",
+ "ranked = fc.sort_values(\"crime_density\", ascending=False)\n",
+ "ranked[\"rank\"] = np.arange(1, len(ranked) + 1)\n",
+ "rows = []\n",
+ "for i, h in hotspots.iterrows():\n",
+ " d = distance_km(ranked[\"lon\"].values, ranked[\"lat\"].values, h.geometry.x, h.geometry.y)\n",
+ " nearest = ranked.iloc[int(d.argmin())]\n",
+ " rows.append({\n",
+ " \"hotspot\": i,\n",
+ " \"share_of_events\": round(h[\"share\"], 3),\n",
+ " \"nearest_place\": int(nearest.name),\n",
+ " \"distance_km\": round(float(d.min()), 2),\n",
+ " \"forecast_rank\": int(nearest[\"rank\"]),\n",
+ " \"percentile\": round(100 * (1 - nearest[\"rank\"] / len(ranked)), 1),\n",
+ " })\n",
+ "pd.DataFrame(rows).set_index(\"hotspot\").sort_values(\"forecast_rank\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c3350c32",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# The full history + forecast of the top forecast place\n",
+ "place = fc[\"crime_density\"].idxmax()\n",
+ "series = pipeline.stseries.xs(place, level=\"places\")\n",
+ "observed = series.loc[: stseries.index.get_level_values(\"t\").max()]\n",
+ "predicted = series.loc[forecasts.index.get_level_values(\"t\").min():]\n",
+ "ax = observed.plot(figsize=(10, 3.5), marker=\"o\", label=\"observed (KDE)\")\n",
+ "predicted.plot(ax=ax, marker=\"*\", markersize=12, linestyle=\"--\", color=\"#d62728\", label=\"forecast\")\n",
+ "ax.set_title(f\"Place {place}: monthly density and 3-month recursive forecast\")\n",
+ "ax.set_xlabel(\"\")\n",
+ "ax.legend()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f01c22b2",
+ "metadata": {},
+ "source": [
+ "## 7. The same pipeline with counts on hexagons\n",
+ "\n",
+ "`QuadratCount` is a drop-in replacement for `KDE`: swap the mapping and the\n",
+ "grid, keep everything else."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8c8e650f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "hex_pipeline = PredictionPipeline(\n",
+ " mapping=QuadratCount(tfreq=\"M\", grid=hexes),\n",
+ " fextraction=PandasFeatureUnion([\n",
+ " (\"ar\", AR(lags=3)),\n",
+ " (\"seasonal\", Seasonality(lags=LAGS)),\n",
+ " (\"trend\", Trend(lags=LAGS)),\n",
+ " ]),\n",
+ " estimator=RandomForestRegressor(n_estimators=100, min_samples_leaf=2, random_state=0, n_jobs=-1),\n",
+ " random_state=0,\n",
+ ").fit(dataset)\n",
+ "print(\"r2 per fold:\", np.round(hex_pipeline.evaluate(\"r2\", cv=3), 3))\n",
+ "hex_fc = hexes.join(hex_pipeline.predict().droplevel(\"t\"))\n",
+ "\n",
+ "fig, ax = plt.subplots(figsize=(7, 8))\n",
+ "hex_fc.plot(ax=ax, column=\"crime_density\", cmap=\"magma\", edgecolor=\"white\", linewidth=0.3, legend=True,\n",
+ " legend_kwds={\"shrink\": 0.6, \"label\": \"forecast events\"})\n",
+ "city.boundary.plot(ax=ax, color=\"black\", linewidth=1)\n",
+ "hotspots.plot(ax=ax, color=\"#00e5ff\", marker=\"*\", markersize=170, zorder=5)\n",
+ "ax.set_title(f\"Forecast events per hexagon for {hex_pipeline.stseries.index.get_level_values('t').max():%B %Y}\")\n",
+ "ax.set_axis_off()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "08c02abe",
+ "metadata": {},
+ "source": [
+ "## Wrapping up\n",
+ "\n",
+ "- `get_city_shape` gave us the study area; `generate_crimes` produced events with known hotspots.\n",
+ "- `Dataset` → grid → `KDE` / `QuadratCount` turned the events into a spatio-temporal series.\n",
+ "- Lagged STL trend/seasonal, difference and autoregressive features fed a scikit-learn pipeline\n",
+ " wrapped by `PredictionPipeline`, evaluated with time series cross-validation.\n",
+ "- `predict()` forecast the following months; the grid points nearest to the true hotspot\n",
+ " centres rank at the top of the forecast, in proportion to each hotspot's share of events.\n",
+ "\n",
+ "To run this on real data, replace the synthetic `crimes` DataFrame with your own\n",
+ "`tag, t, lon, lat` table and model each crime type (`tag`) separately, as the\n",
+ "framework recommends. See the [documentation](https://adaj.github.io/predspot/)\n",
+ "for the details of every step."
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "name": "python"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/mkdocs.yml b/mkdocs.yml
index 78309a3..d91c589 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -48,6 +48,8 @@ nav:
- Grids and spatio-temporal mapping: guide/mapping.md
- Feature engineering: guide/features.md
- Prediction and evaluation: guide/modelling.md
+ - Examples:
+ - Natal, Brazil (end to end): examples/natal.ipynb
- API reference:
- predspot: api/predspot.md
- dataset_preparation: api/dataset_preparation.md
@@ -62,8 +64,15 @@ nav:
- Contributing: contributing.md
- Citing Predspot: citing.md
+hooks:
+ - docs/hooks/readme.py
+
plugins:
- search
+ - mkdocs-jupyter:
+ execute: false
+ include_source: true
+ include_requirejs: true
- autorefs
- mkdocstrings:
handlers:
diff --git a/pyproject.toml b/pyproject.toml
index 5c61551..89b35b0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -49,6 +49,7 @@ docs = [
"mkdocs-material>=9.5,<10",
"mkdocstrings[python]>=0.26",
"mkdocs-autorefs>=1.2",
+ "mkdocs-jupyter>=0.25",
]
dev = [
"pytest>=8",
@@ -88,6 +89,7 @@ branch = true
[tool.ruff]
line-length = 100
+exclude = ["examples"] # tutorial code: readability over line length
target-version = "py310"
src = ["src", "tests"]
diff --git a/src/predspot/__init__.py b/src/predspot/__init__.py
index f85d774..5d2ca3d 100644
--- a/src/predspot/__init__.py
+++ b/src/predspot/__init__.py
@@ -23,6 +23,7 @@
create_gridhexagonal,
create_gridpoints,
create_gridsquares,
+ get_city_shape,
load_study_area,
)
from predspot.dataset_preparation import Dataset
@@ -42,6 +43,7 @@
"create_gridhexagonal",
"create_gridsquares",
"load_study_area",
+ "get_city_shape",
"generate_crimes",
"synthetic",
"crime_mapping",
diff --git a/src/predspot/crime_mapping.py b/src/predspot/crime_mapping.py
index 710d7c7..f2ccfa9 100644
--- a/src/predspot/crime_mapping.py
+++ b/src/predspot/crime_mapping.py
@@ -138,6 +138,36 @@ def load_study_area(place, crs=WGS84, which_result=None):
return gdf.to_crs(crs)
+def get_city_shape(place_query):
+ """
+ Fetch the shape (polygon) of a city or region from OpenStreetMap.
+
+ A thin wrapper around ``osmnx.geocode_to_gdf`` that returns the raw
+ Nominatim result. Prefer [`load_study_area`][predspot.crime_mapping.load_study_area]
+ when you want the result validated (polygon geometry, tidy columns).
+ Requires the optional ``osmnx`` dependency (``pip install predspot[osm]``).
+
+ Args:
+ place_query (str): Name of the place, in a format Nominatim accepts,
+ e.g. ``"Natal, RN, Brazil"`` or ``"Rio Grande do Norte, Brazil"``.
+
+ Returns:
+ GeoDataFrame: One row (or more, if the query is ambiguous) with the
+ place geometry (Polygon/MultiPolygon) in WGS84.
+
+ Example:
+ >>> city = get_city_shape("Natal, RN, Brazil")
+ >>> city.geometry.iloc[0] # shapely Polygon/MultiPolygon
+ """
+ try:
+ import osmnx as ox
+ except ImportError as exc:
+ raise ImportError(
+ "get_city_shape requires the optional dependency `osmnx`: pip install predspot[osm]"
+ ) from exc
+ return ox.geocode_to_gdf(place_query)
+
+
def _check_bbox(bbox):
if not isinstance(bbox, gpd.GeoDataFrame):
raise TypeError("bbox must be a geopandas GeoDataFrame.")
diff --git a/src/predspot/ml_modelling.py b/src/predspot/ml_modelling.py
index ec7e4e1..32c6e6c 100644
--- a/src/predspot/ml_modelling.py
+++ b/src/predspot/ml_modelling.py
@@ -227,16 +227,19 @@ def evaluate(self, scoring="r2", cv=5):
refitted on the original data afterwards.
Args:
- scoring (str): ``'r2'`` or ``'mse'``.
+ scoring (str or list): ``'r2'``, ``'mse'`` or a list of them.
cv (int): Number of folds (must be lower than the number of periods).
Returns:
- list: One score per fold.
+ list or pandas.DataFrame: One score per fold; with a list of
+ scorings, a DataFrame with one column per scoring and one row
+ per fold.
"""
self._check_fitted()
- if scoring not in SCORERS:
+ scorings = [scoring] if isinstance(scoring, str) else list(scoring)
+ unknown = [s for s in scorings if s not in SCORERS]
+ if unknown or not scorings:
raise ValueError('invalid scoring. Try "r2" or "mse".')
- scorer = SCORERS[scoring]
timestamps = (
self._X.index.get_level_values("t")
.unique()
@@ -245,7 +248,7 @@ def evaluate(self, scoring="r2", cv=5):
)
if not isinstance(cv, int) or cv >= len(timestamps):
raise ValueError("cv must be an integer lower than the number of periods.")
- scores = []
+ scores = {name: [] for name in scorings}
for train_t, test_t in TimeSeriesSplit(cv).split(timestamps):
X_train = self._X.loc[idx[timestamps[train_t], :], :].sample(
frac=1, random_state=self.random_state
@@ -255,7 +258,10 @@ def evaluate(self, scoring="r2", cv=5):
y_test = self._stseries.loc[X_test.index]
self.estimator.fit(X_train, y_train)
y_pred = self.estimator.predict(X_test)
- scores.append(scorer(y_test, y_pred))
- logger.debug("%s-fold CV %s scores: %s", cv, scoring, scores)
+ for name in scorings:
+ scores[name].append(SCORERS[name](y_test, y_pred))
+ logger.debug("%s-fold CV scores: %s", cv, scores)
self.fit(self._dataset) # back to normal
- return scores
+ if isinstance(scoring, str):
+ return scores[scoring]
+ return pd.DataFrame(scores, index=[f"fold {i + 1}" for i in range(cv)])
diff --git a/tests/test_load_study_area.py b/tests/test_load_study_area.py
index 70efc14..92b7640 100644
--- a/tests/test_load_study_area.py
+++ b/tests/test_load_study_area.py
@@ -6,7 +6,7 @@
import pytest
from shapely.geometry import Point, box
-from predspot.crime_mapping import load_study_area
+from predspot.crime_mapping import get_city_shape, load_study_area
def _fake_osmnx(monkeypatch, geometry):
@@ -70,3 +70,17 @@ def test_load_study_area_network():
assert len(area) == 1
assert area.geometry.iloc[0].geom_type in ("Polygon", "MultiPolygon")
assert area.geometry.iloc[0].contains(Point(-35.2094, -5.7945)) # city centre
+
+
+def test_get_city_shape_returns_raw_geocode(monkeypatch):
+ calls = _fake_osmnx(monkeypatch, box(-35.3, -5.9, -35.2, -5.8))
+ city = get_city_shape("Natal, RN, Brazil")
+ assert calls["query"] == "Natal, RN, Brazil"
+ assert city.geometry.iloc[0].geom_type == "Polygon"
+ assert "place_rank" in city.columns # raw Nominatim columns are kept
+
+
+def test_get_city_shape_without_osmnx(monkeypatch):
+ monkeypatch.setitem(sys.modules, "osmnx", None)
+ with pytest.raises(ImportError, match="predspot\\[osm\\]"):
+ get_city_shape("Natal, RN, Brazil")
diff --git a/tests/test_ml_modelling.py b/tests/test_ml_modelling.py
index 0a1479f..d6cac36 100644
--- a/tests/test_ml_modelling.py
+++ b/tests/test_ml_modelling.py
@@ -73,8 +73,12 @@ def test_evaluate(dataset, study_area):
assert len(scores) == 3
mse = pipe.evaluate("mse", cv=3)
assert all(s >= 0 for s in mse)
+ both = pipe.evaluate(["r2", "mse"], cv=3)
+ assert list(both.columns) == ["r2", "mse"] and len(both) == 3
with pytest.raises(ValueError):
pipe.evaluate("mae")
+ with pytest.raises(ValueError):
+ pipe.evaluate([])
with pytest.raises(ValueError):
pipe.evaluate("r2", cv=100)
# evaluate refits on the full data: predictions still work afterwards
