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.

In [1]:

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__)

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__)

predspot 0.2.0 | geopandas 1.1.4 | pandas 3.0.3

1. Study area: the boundary of Natal¶

get_city_shape geocodes a place name with osmnx (Nominatim) and returns its boundary as a GeoDataFrame in WGS84. The query should be specific enough for the first match to be the municipality.

In [2]:

city = get_city_shape("Natal, RN, Brazil")
city

city = get_city_shape("Natal, RN, Brazil") city

Out[2]:

	geometry	bbox_west	bbox_south	bbox_east	bbox_north	place_id	osm_type	osm_id	lat	lon	class	type	place_rank	importance	addresstype	name	display_name
0	POLYGON ((-35.29124 -5.73221, -35.29121 -5.732...	-35.291237	-5.900211	-35.153103	-5.702738	14119496	relation	301091	-5.805398	-35.20809	boundary	administrative	16	0.583903	municipality	Natal	Natal, Rio Grande do Norte, Northeast Region, ...

In [3]:

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()

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()

geometry: Polygon | vertices: 786
bounds (W, S, E, N): [-35.2912  -5.9002 -35.1531  -5.7027]
area: 166.7 km²

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.

In [4]:

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)

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)

Out[4]:

	tag	t	lon	lat
0	robbery	2019-01-01 08:07:29	-35.262821	-5.826621
1	burglary	2019-01-01 11:04:33	-35.184192	-5.839420
2	assault	2019-01-01 15:45:04	-35.263080	-5.833051
3	robbery	2019-01-01 19:57:59	-35.264664	-5.837595
4	burglary	2019-01-01 20:23:07	-35.256722	-5.819229
5	burglary	2019-01-01 22:25:05	-35.262523	-5.825538
6	burglary	2019-01-01 23:33:14	-35.256481	-5.821109
7	burglary	2019-01-02 00:35:11	-35.192001	-5.805592
8	homicide	2019-01-02 09:55:04	-35.253808	-5.732047
9	homicide	2019-01-02 11:06:07	-35.206620	-5.760610

In [5]:

print(crimes.shape)
crimes.describe(include="all").T

print(crimes.shape) crimes.describe(include="all").T

(12000, 4)

Out[5]:

	count	unique	top	freq	mean	min	25%	50%	75%	max	std
tag	12000	4	burglary	5397	NaN	NaN	NaN	NaN	NaN	NaN	NaN
t	12000	NaN	NaN	NaN	2020-08-04 09:23:35.348416768	2019-01-01 08:07:29	2019-11-18 00:35:28.500000	2020-09-02 18:24:44	2021-04-22 17:33:00.249999872	2021-12-30 22:51:11	NaN
lon	12000.0	NaN	NaN	NaN	-35.23946	-35.290021	-35.263962	-35.252215	-35.222052	-35.153427	0.030905
lat	12000.0	NaN	NaN	NaN	-5.822253	-5.899403	-5.841009	-5.826307	-5.807527	-5.704011	0.039883

In [6]:

hotspots

hotspots

Out[6]:

	sd_km	share	geometry
0	0.7	0.312219	POINT (-35.26454 -5.82894)
1	0.7	0.094299	POINT (-35.18172 -5.8832)
2	0.7	0.051078	POINT (-35.22805 -5.866)
3	0.7	0.114319	POINT (-35.26133 -5.81869)
4	0.7	0.078085	POINT (-35.23374 -5.77675)

In [7]:

crimes["tag"].value_counts().to_frame("events").assign(share=lambda d: (d["events"] / len(crimes)).round(3))

crimes["tag"].value_counts().to_frame("events").assign(share=lambda d: (d["events"] / len(crimes)).round(3))

Out[7]:

	events	share
tag
burglary	5397	0.450
robbery	3564	0.297
assault	2399	0.200
homicide	640	0.053

2.1 Where and when do the events happen?¶

In [8]:

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()

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()

In [9]:

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()

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()

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.

In [10]:

dataset = Dataset(crimes, city)
dataset

dataset = Dataset(crimes, city) dataset

Out[10]:

predspot.Dataset<
  crimes = GeoDataFrame(12000),
    >> {'burglary': 5397, 'robbery': 3564, 'assault': 2399, 'homicide': 640}
  study_area = GeoDataFrame(1),
>

In [11]:

dataset.crimes.head()

dataset.crimes.head()

Out[11]:

	tag	t	lon	lat	geometry
0	robbery	2019-01-01 08:07:29	-35.262821	-5.826621	POINT (-35.26282 -5.82662)
1	burglary	2019-01-01 11:04:33	-35.184192	-5.839420	POINT (-35.18419 -5.83942)
2	assault	2019-01-01 15:45:04	-35.263080	-5.833051	POINT (-35.26308 -5.83305)
3	robbery	2019-01-01 19:57:59	-35.264664	-5.837595	POINT (-35.26466 -5.8376)
4	burglary	2019-01-01 20:23:07	-35.256722	-5.819229	POINT (-35.25672 -5.81923)

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.

In [12]:

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()

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()

638 grid points at 500 m | 207 hexagons of 1 km²

Out[12]:

	lon	lat	geometry
places
55	-35.180730	-5.895619	POINT (-35.18073 -5.89562)
56	-35.176126	-5.895619	POINT (-35.17613 -5.89562)
57	-35.171521	-5.895619	POINT (-35.17152 -5.89562)
58	-35.166917	-5.895619	POINT (-35.16692 -5.89562)
59	-35.162312	-5.895619	POINT (-35.16231 -5.89562)

In [13]:

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()

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()

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.

In [14]:

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)

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)

Series (22968,) | KDE factor: 0.3824

Out[14]:

t           places
2019-01-31  55        63.536374
            56        62.723997
            57        54.013318
            58        40.779398
            59        27.278340
            60        16.455813
            85        71.430944
            86        78.158377
Name: crime_density, dtype: float64

In [15]:

# The same series in wide form: one row per month, one column per place
wide = stseries.unstack("places")
wide.iloc[:6, :8]

# The same series in wide form: one row per month, one column per place wide = stseries.unstack("places") wide.iloc[:6, :8]

Out[15]:

places	55	56	57	58	59	60	85	86
t
2019-01-31	63.536374	62.723997	54.013318	40.779398	27.278340	16.455813	71.430944	78.158377
2019-02-28	50.339803	50.157606	44.957892	37.139687	29.127228	22.190242	54.798466	59.127446
2019-03-31	77.929146	79.196237	71.679799	58.110409	42.398863	27.944718	82.829593	91.021441
2019-04-30	43.955434	42.369459	35.176428	25.057215	15.268911	7.936542	49.364891	53.603443
2019-05-31	57.834267	53.098159	43.560881	32.324654	22.087052	14.227696	68.609979	69.014337
2019-06-30	45.592241	50.773863	50.349845	43.962284	33.557154	22.297426	51.858407	61.938736

In [16]:

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]

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]

events per month recovered by the counts: [320, 312, 290] ...

Out[16]:

places	37	38	39	40	41	50	51	52
t
2019-01-31	0.0	0.0	0.0	1.0	0.0	0.0	0.0	1.0
2019-02-28	0.0	0.0	0.0	1.0	0.0	0.0	0.0	2.0
2019-03-31	0.0	1.0	1.0	0.0	0.0	0.0	1.0	0.0
2019-04-30	0.0	0.0	0.0	1.0	0.0	0.0	0.0	0.0
2019-05-31	0.0	0.0	0.0	1.0	0.0	0.0	0.0	2.0
2019-06-30	0.0	0.0	0.0	0.0	0.0	0.0	0.0	3.0

In [17]:

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()

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()

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).

In [18]:

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()

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()

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.

In [19]:

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)

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)

(19140, 21) | periods: 2019-08-31 -> 2022-01-31

Out[19]:

		ar_1	ar_2	ar_3	seasonal_1	seasonal_2	seasonal_3	seasonal_4	seasonal_5	seasonal_6	trend_1	...	trend_3	trend_4	trend_5	trend_6	diff_1	diff_2	diff_3	diff_4	diff_5	diff_6
t	places
2019-08-31	55	125.066978	33.379474	75.417634	17.792389	-30.120816	-9.378083	-12.053399	28.559599	8.725865	90.150882	...	85.306595	82.496749	79.764318	77.129447	91.687504	-42.038160	5.464522	-54.822791	39.967023	-6.694205
	56	119.722317	43.676995	62.169137	15.704561	-23.082344	-16.896765	-5.064126	32.536963	0.034839	91.472046	...	82.901002	78.376722	73.961567	69.744306	76.045322	-18.492142	-1.109927	-55.373884	41.442033	-10.483457
	57	74.955493	48.107382	44.961841	1.650338	-7.189462	-9.907225	-3.787340	30.708964	-8.703507	69.222787	...	63.221369	60.187650	57.331271	54.778577	26.848111	3.145541	6.223854	-59.937445	39.288169	0.550026
	58	32.885539	37.511566	27.005920	-8.958600	0.552562	-2.519541	-4.530974	23.950183	-4.966147	41.233126	...	39.449199	38.574753	37.873135	37.460754	-4.626026	10.505646	11.468861	-58.689203	30.512404	15.810417
	59	10.548101	19.797947	14.288868	-10.264510	-0.437829	-0.643550	-2.422139	13.554383	4.466334	21.870406	...	22.048987	22.169159	22.412970	22.869190	-9.249846	5.509079	10.369965	-44.347235	15.689751	20.004471
	60	2.208611	7.265361	10.962395	-8.091273	-2.416711	-0.094430	-1.558245	3.652393	12.938206	13.415543	...	13.889126	14.205965	14.612110	15.169423	-5.056750	-3.697034	10.377684	-26.775176	-1.408718	19.260890
	85	151.142469	43.625061	127.863353	19.490251	-48.072805	13.097882	-20.990928	31.129092	5.988812	114.699806	...	112.197090	110.646175	109.190448	107.828026	107.517409	-84.238292	43.927759	-83.758478	65.082126	-23.790069
	86	197.682505	71.207880	119.173254	30.181598	-48.819994	-14.219399	-10.504940	38.803089	5.791861	147.416078	...	135.466737	129.137953	122.845369	116.675448	126.474625	-47.965375	14.792964	-81.530166	67.676212	-37.023901

8 rows × 21 columns

In [20]:

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

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

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.

In [21]:

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())

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())

training rows: 18502
next period to forecast: 2022-01-31

In [22]:

scores = pipeline.evaluate(["r2", "mse"], cv=3)   # time series CV: train on earlier folds, test on the next
scores.loc["mean"] = scores.mean()
scores

scores = pipeline.evaluate(["r2", "mse"], cv=3) # time series CV: train on earlier folds, test on the next scores.loc["mean"] = scores.mean() scores

Out[22]:

	r2	mse
fold 1	0.968710	364.186782
fold 2	0.983170	177.396760
fold 3	0.985974	142.169684
mean	0.979284	227.917742

In [23]:

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)

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)

Out[23]:

	importance
trend_1	0.802854
trend_2	0.112126
trend_3	0.033980
trend_4	0.016472
trend_6	0.015978
seasonal_6	0.007559
diff_6	0.002337
ar_2	0.002062
diff_5	0.002022
ar_1	0.001691
seasonal_2	0.001469
seasonal_4	0.001451

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.

In [24]:

forecasts = pd.concat([pipeline.predict() for _ in range(3)])
forecasts.groupby("t").describe().round(2)

forecasts = pd.concat([pipeline.predict() for _ in range(3)]) forecasts.groupby("t").describe().round(2)

Out[24]:

	crime_density
	count	mean	std	min	25%	50%	75%	max
t
2022-01-31	638.0	67.37	96.78	3.64	24.29	37.09	62.40	638.06
2022-02-28	638.0	66.49	104.55	2.82	17.54	34.48	61.69	711.16
2022-03-31	638.0	65.65	106.42	2.80	16.75	32.32	58.58	715.27

In [25]:

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

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

Top-10 places for January 2022:

Out[25]:

	lon	lat	crime_density
places
502	-35.2636	-5.8267	638.0591
533	-35.2636	-5.8221	629.9240
501	-35.2682	-5.8267	611.3991
471	-35.2636	-5.8313	608.1633
503	-35.2590	-5.8267	602.1349
534	-35.2590	-5.8221	580.6709
532	-35.2682	-5.8221	571.8658
470	-35.2682	-5.8313	539.2236
472	-35.2590	-5.8313	522.6021
564	-35.2636	-5.8175	506.5837

In [26]:

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()

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()

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.

In [27]:

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")

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")

Out[27]:

	share_of_events	nearest_place	distance_km	forecast_rank	percentile
hotspot
0	0.312	502	0.26	1	99.8
3	0.114	564	0.28	10	98.4
1	0.094	148	0.19	44	93.1
4	0.078	849	0.26	58	90.9
2	0.051	231	0.27	96	85.0

In [28]:

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

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

Out[28]:

<matplotlib.legend.Legend at 0x12dbd8e10>

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.

In [29]:

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()

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()

r2 per fold: [0.868 0.916 0.929]

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 for the details of every step.