diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 343b82a..a88ddbe 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -37,7 +37,7 @@ jobs: python-version: ${{ matrix.python-version }} cache: pip - name: Install - run: python -m pip install -e ".[dev,contour]" + run: python -m pip install -e ".[dev,osm,contour]" - name: Test run: pytest --cov --cov-report=term-missing diff --git a/CHANGELOG.md b/CHANGELOG.md index 5b45ac8..2c94242 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,11 @@ 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.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 + area: Gaussian hotspots plus uniform background, with trend, annual cycle, + day-of-week and hour-of-day patterns; reproducible with `seed`. - `QuadratCount` mapping (event counts per cell) as a first-class alternative to `KDE`, usable with hexagonal (`create_gridhexagonal`) and square (`create_gridsquares`) grids inside `PredictionPipeline`. diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 5a46c84..b510773 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -8,7 +8,7 @@ Thanks for your interest in Predspot! Issues and pull requests are welcome. git clone https://github.com/adaj/predspot.git cd predspot python -m venv .venv && source .venv/bin/activate -pip install -e ".[dev,contour]" +pip install -e ".[dev,osm,contour]" ``` ## Checks @@ -17,6 +17,7 @@ pip install -e ".[dev,contour]" ruff check src tests # lint ruff format src tests # format pytest # tests (~10 s) +PREDSPOT_NETWORK_TESTS=1 pytest tests/test_load_study_area.py # also query OpenStreetMap ``` CI runs the same checks on Python 3.10 to 3.13 for every pull request. diff --git a/README.md b/README.md index 81c7214..d0965e0 100644 --- a/README.md +++ b/README.md @@ -67,6 +67,23 @@ from predspot.crime_mapping import QuadratCount, create_gridhexagonal mapping = QuadratCount(tfreq='W', grid=create_gridhexagonal(study_area_gdf, resolution=1)) ``` +### Study area from OpenStreetMap and synthetic data + +You do not need real data to try Predspot. Fetch a city boundary from +OpenStreetMap (`pip install "predspot[osm]"`) and generate synthetic events +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 + +study_area = load_study_area("Natal, Rio Grande do Norte, Brazil") +crimes = generate_crimes(study_area, n_events=5000, n_hotspots=4, + start="2019-01-01", end="2020-12-31", seed=0) +dataset = Dataset(crimes, study_area) +dataset.plot() +``` + Or run the default pipeline in one call: ```python @@ -80,7 +97,7 @@ print(pipeline.evaluate('r2', cv=3)) ## Development ⚡ -Predspot has four main modules: +Predspot has five main modules: `dataset_preparation`: Module for preparing and managing crime datasets and study areas. @@ -90,12 +107,15 @@ Predspot has four main modules: `ml_modelling`: Module that implements the prediction pipeline and model evaluation. +`synthetic`: Module that generates synthetic crime events (hotspots + temporal patterns) inside any study area. + ### Installation steps 🛠️ Predspot requires Python 3.10 or newer. ```bash pip install predspot # from PyPI +pip install "predspot[osm]" # + study areas from OpenStreetMap (osmnx) pip install "predspot[contour]" # + GeoJSON contour export (geojsoncontour) ``` @@ -104,7 +124,7 @@ From source, for development: ```bash git clone https://github.com/adaj/predspot.git cd predspot -pip install -e ".[dev,contour]" +pip install -e ".[dev,osm,contour]" ``` Core dependencies (installed automatically): pandas, geopandas, shapely, diff --git a/src/predspot/__init__.py b/src/predspot/__init__.py index 2dae115..f85d774 100644 --- a/src/predspot/__init__.py +++ b/src/predspot/__init__.py @@ -14,6 +14,7 @@ dataset_preparation, feature_engineering, ml_modelling, + synthetic, utilities, ) from predspot.crime_mapping import ( @@ -22,9 +23,11 @@ create_gridhexagonal, create_gridpoints, create_gridsquares, + load_study_area, ) from predspot.dataset_preparation import Dataset from predspot.ml_modelling import PredictionPipeline +from predspot.synthetic import generate_crimes from predspot.utilities import PandasFeatureUnion __version__ = "0.2.0" @@ -38,6 +41,9 @@ "create_gridpoints", "create_gridhexagonal", "create_gridsquares", + "load_study_area", + "generate_crimes", + "synthetic", "crime_mapping", "dataset_preparation", "feature_engineering", diff --git a/src/predspot/crime_mapping.py b/src/predspot/crime_mapping.py index 941538d..78d16f6 100644 --- a/src/predspot/crime_mapping.py +++ b/src/predspot/crime_mapping.py @@ -14,6 +14,9 @@ Both produce the same output format, a :class:`pandas.Series` named ``crime_density`` indexed by ``(t, places)``, so they are interchangeable inside :class:`predspot.ml_modelling.PredictionPipeline`. + +The study area itself can be fetched from OpenStreetMap with +:func:`load_study_area` (requires the optional ``osmnx`` dependency). """ __author__ = "Adelson Araujo" @@ -83,6 +86,54 @@ def tfreq_offset(tfreq): return pd.offsets.Day(1) +def load_study_area(place, crs=WGS84, which_result=None): + """ + Fetch the boundary polygon of a place from OpenStreetMap. + + Uses `osmnx `_ (optional dependency: + ``pip install predspot[osm]``) to geocode ``place`` with Nominatim and + return its administrative boundary, ready to be used as the + ``study_area`` of :class:`predspot.Dataset` or as the ``bbox`` of the + ``create_grid*`` functions. + + Args: + place (str or list): Name of the place as Nominatim understands it, + e.g. ``"Natal, Rio Grande do Norte, Brazil"``. A list of names + returns one row per place. + crs: CRS of the returned GeoDataFrame (default WGS84). + which_result (int, optional): Forwarded to + ``osmnx.geocode_to_gdf`` to pick a specific Nominatim match when + the first one is not the boundary you want. + + Returns: + GeoDataFrame: One row per place with ``name``, ``display_name``, + ``osm_type``, ``osm_id`` and a (Multi)Polygon ``geometry``. + + Raises: + ImportError: If ``osmnx`` is not installed. + ValueError: If Nominatim returns a point instead of a boundary + polygon for the query. + """ + try: + import osmnx as ox + except ImportError as exc: + raise ImportError( + "load_study_area requires the optional dependency `osmnx`: pip install predspot[osm]" + ) from exc + logger.debug("Geocoding study area %r with osmnx", place) + gdf = ox.geocode_to_gdf(place, which_result=which_result) + if not gdf.geom_type.isin(["Polygon", "MultiPolygon"]).all(): + bad = gdf.loc[~gdf.geom_type.isin(["Polygon", "MultiPolygon"]), "display_name"] + raise ValueError( + "OpenStreetMap returned a non-polygon geometry for " + f"{bad.tolist()}. Try a more specific query (e.g. add the state " + "and country) or a different `which_result`." + ) + columns = [c for c in ("name", "display_name", "osm_type", "osm_id") if c in gdf.columns] + gdf = gdf[columns + ["geometry"]].reset_index(drop=True) + return gdf.to_crs(crs) + + def _check_bbox(bbox): if not isinstance(bbox, gpd.GeoDataFrame): raise TypeError("bbox must be a geopandas GeoDataFrame.") diff --git a/src/predspot/pipeline.py b/src/predspot/pipeline.py index 6b51933..d7354d6 100644 --- a/src/predspot/pipeline.py +++ b/src/predspot/pipeline.py @@ -25,6 +25,7 @@ from sklearn.preprocessing import QuantileTransformer from predspot import crime_mapping, dataset_preparation, feature_engineering, ml_modelling +from predspot.synthetic import generate_crimes from predspot.utilities import PandasFeatureUnion logger = logging.getLogger(__name__) @@ -35,7 +36,10 @@ def generate_testdata(n_points, start_time, end_time, bounds=DEFAULT_BOUNDS, seed=None): """ - Generate uniformly random crime events inside a rectangular study area. + Generate synthetic crime events inside a rectangular study area. + + A thin wrapper around :func:`predspot.synthetic.generate_crimes` with + three hotspots and the default temporal patterns. Args: n_points (int): Number of events. @@ -48,27 +52,13 @@ def generate_testdata(n_points, start_time, end_time, bounds=DEFAULT_BOUNDS, see tuple: ``(crimes, study_area)`` — a DataFrame with ``tag``, ``t``, ``lon``, ``lat`` and a one-row GeoDataFrame with the study area. """ - rng = np.random.default_rng(seed) west, south, east, north = bounds study_area = gpd.GeoDataFrame( {"name": ["study_area"]}, geometry=[box(west, south, east, north)], crs="EPSG:4326" ) - start, end = pd.Timestamp(start_time), pd.Timestamp(end_time) - seconds = rng.integers(0, int((end - start).total_seconds()), n_points) - tags = rng.choice( - ["burglary", "assault", "drugs", "homicide"], - size=n_points, - p=np.array([1000, 100, 10, 1]) / 1111, - ) - crimes = pd.DataFrame( - { - "tag": tags, - "t": start + pd.to_timedelta(seconds, unit="s"), - "lon": rng.uniform(west, east, n_points), - "lat": rng.uniform(south, north, n_points), - } + crimes = generate_crimes( + study_area, n_events=n_points, start=start_time, end=end_time, seed=seed ) - logger.debug("Generated %d synthetic events", n_points) return crimes, study_area diff --git a/src/predspot/synthetic.py b/src/predspot/synthetic.py new file mode 100644 index 0000000..bea3b9a --- /dev/null +++ b/src/predspot/synthetic.py @@ -0,0 +1,340 @@ +""" +Synthetic Data Module +===================== + +Generates realistic-looking synthetic crime events inside a study area, so +that Predspot can be tried, demonstrated and tested without real data. + +The generator is a simple inhomogeneous space-time point process: + +* **Space** — a mixture of ``n_hotspots`` Gaussian hotspots (centres drawn + uniformly inside the study area) plus a uniform background. The fraction + of events that belong to hotspots is ``hotspot_share``. +* **Time** — an intensity built from a linear trend, an annual cycle, a + day-of-week profile and an hour-of-day profile. Timestamps are drawn by + thinning uniform candidates, which is exact for a fixed number of events. + +Example:: + + from predspot.synthetic import generate_crimes + crimes = generate_crimes(study_area, n_events=5000, seed=0) + dataset = Dataset(crimes, study_area) +""" + +__author__ = "Adelson Araujo" + +import logging +import math + +import geopandas as gpd +import numpy as np +import pandas as pd +import shapely + +from predspot.crime_mapping import KM_PER_DEG_LAT, KM_PER_DEG_LON, WGS84, _check_bbox + +logger = logging.getLogger(__name__) + +DEFAULT_TAGS = {"burglary": 0.45, "robbery": 0.30, "assault": 0.20, "homicide": 0.05} + +# Relative intensity Monday..Sunday (normalised to mean 1 internally). +DEFAULT_WEEKLY_PROFILE = (0.90, 0.85, 0.90, 0.95, 1.10, 1.30, 1.00) + +# Relative intensity by hour of day, 0..23: quiet early morning, busy evening. +DEFAULT_HOURLY_PROFILE = ( + 0.7, 0.5, 0.4, 0.3, 0.25, 0.3, 0.45, 0.7, 0.9, 1.0, 1.05, 1.1, + 1.1, 1.05, 1.1, 1.15, 1.25, 1.4, 1.55, 1.7, 1.75, 1.6, 1.3, 1.0, +) # fmt: skip + + +def _study_polygon(study_area): + """Return the union of the study area geometries as a WGS84 shapely geometry.""" + _check_bbox(study_area) + return study_area.to_crs(WGS84).geometry.union_all() + + +def sample_points_in_polygon(polygon, n, rng, max_iterations=1000): + """ + Draw ``n`` points uniformly inside a polygon by rejection sampling. + + Args: + polygon (shapely.Geometry): Polygon in WGS84. + n (int): Number of points. + rng (numpy.random.Generator): Random generator. + max_iterations (int): Safety cap on rejection rounds. + + Returns: + tuple: ``(lon, lat)`` arrays of length ``n``. + """ + minx, miny, maxx, maxy = polygon.bounds + lon = np.empty(0) + lat = np.empty(0) + fill = polygon.area / ((maxx - minx) * (maxy - miny)) if polygon.area > 0 else 1.0 + batch = int(max(n / max(fill, 1e-3) * 1.2, 64)) + for _ in range(max_iterations): + cx = rng.uniform(minx, maxx, batch) + cy = rng.uniform(miny, maxy, batch) + inside = shapely.contains_xy(polygon, cx, cy) + lon = np.concatenate([lon, cx[inside]]) + lat = np.concatenate([lat, cy[inside]]) + if len(lon) >= n: + return lon[:n], lat[:n] + raise RuntimeError("Could not sample enough points inside the study area.") + + +def sample_points_around(polygon, centers, sd_km, n_per_center, rng, max_iterations=1000): + """ + Draw points from Gaussian clouds around hotspot centres, kept inside the polygon. + + Args: + polygon (shapely.Geometry): Study area in WGS84. + centers (numpy.ndarray): ``(k, 2)`` array of ``(lon, lat)`` centres. + sd_km (float or array): Standard deviation of each cloud in km. + n_per_center (array): Number of points to draw per centre. + rng (numpy.random.Generator): Random generator. + max_iterations (int): Safety cap on rejection rounds. + + Returns: + tuple: ``(lon, lat)`` arrays. + """ + sd_km = np.broadcast_to(np.asarray(sd_km, dtype=float), (len(centers),)) + lons, lats = [], [] + for (c_lon, c_lat), sd, n in zip(centers, sd_km, n_per_center, strict=True): + if n == 0: + continue + sd_lat = sd / KM_PER_DEG_LAT + sd_lon = sd / (KM_PER_DEG_LON * max(math.cos(math.radians(c_lat)), 0.05)) + lon = np.empty(0) + lat = np.empty(0) + remaining = int(n) + for _ in range(max_iterations): + cx = rng.normal(c_lon, sd_lon, remaining * 2) + cy = rng.normal(c_lat, sd_lat, remaining * 2) + inside = shapely.contains_xy(polygon, cx, cy) + lon = np.concatenate([lon, cx[inside]]) + lat = np.concatenate([lat, cy[inside]]) + if len(lon) >= n: + break + remaining = int(n - len(lon)) + else: + raise RuntimeError("Could not sample enough hotspot points inside the study area.") + lons.append(lon[:n]) + lats.append(lat[:n]) + if not lons: + return np.empty(0), np.empty(0) + return np.concatenate(lons), np.concatenate(lats) + + +def temporal_intensity( + timestamps, + start, + end, + trend=0.0, + annual_amplitude=0.0, + annual_peak_month=1, + weekly_profile=None, + hourly_profile=None, +): + """ + Relative event intensity at each timestamp (mean around 1). + + Args: + timestamps (DatetimeIndex): Times to evaluate. + start, end (Timestamp): Bounds of the simulation, used for the trend. + trend (float): Relative change of the intensity from ``start`` to + ``end`` (``0.5`` means +50% at the end, ``-0.3`` means -30%). + annual_amplitude (float): Amplitude of the annual cosine (0-1). + annual_peak_month (int): Month (1-12) where the annual cycle peaks. + weekly_profile (sequence): 7 relative weights, Monday to Sunday. + hourly_profile (sequence): 24 relative weights, hour 0 to 23. + + Returns: + numpy.ndarray: Intensity values, one per timestamp. + """ + timestamps = pd.DatetimeIndex(timestamps) + span = max((end - start).total_seconds(), 1.0) + frac = (timestamps - start).total_seconds() / span + lam = np.clip(1.0 + trend * frac, 0.0, None) + if annual_amplitude: + peak_doy = (annual_peak_month - 1) * 365.25 / 12 + 15 + lam = lam * ( + 1 + annual_amplitude * np.cos(2 * np.pi * (timestamps.dayofyear - peak_doy) / 365.25) + ) + if weekly_profile is not None: + weekly = np.asarray(weekly_profile, dtype=float) + if weekly.shape != (7,): + raise ValueError("weekly_profile must have 7 values (Monday..Sunday).") + lam = lam * (weekly / weekly.mean())[timestamps.dayofweek] + if hourly_profile is not None: + hourly = np.asarray(hourly_profile, dtype=float) + if hourly.shape != (24,): + raise ValueError("hourly_profile must have 24 values (hour 0..23).") + lam = lam * (hourly / hourly.mean())[timestamps.hour] + return np.asarray(lam) + + +def sample_timestamps(n, start, end, rng, max_iterations=1000, **intensity_kwargs): + """ + Draw ``n`` timestamps from an inhomogeneous process by thinning. + + Args: + n (int): Number of timestamps. + start, end (Timestamp): Simulation bounds. + rng (numpy.random.Generator): Random generator. + max_iterations (int): Safety cap on thinning rounds. + **intensity_kwargs: Forwarded to :func:`temporal_intensity`. + + Returns: + pandas.DatetimeIndex: ``n`` timestamps, unsorted. + """ + span = (end - start).total_seconds() + if span <= 0: + raise ValueError("end must be after start.") + # Upper bound of the intensity: product of the maxima of each factor. + lam_max = max(1.0, 1.0 + intensity_kwargs.get("trend", 0.0)) + lam_max *= 1 + abs(intensity_kwargs.get("annual_amplitude", 0.0)) + for key, size in (("weekly_profile", 7), ("hourly_profile", 24)): + profile = intensity_kwargs.get(key) + if profile is not None: + profile = np.asarray(profile, dtype=float) + if profile.shape != (size,): + raise ValueError(f"{key} must have {size} values.") + lam_max *= profile.max() / profile.mean() + accepted = [] + total = 0 + batch = int(max(n * lam_max * 1.2, 64)) + for _ in range(max_iterations): + candidates = start + pd.to_timedelta(rng.uniform(0, span, batch), unit="s") + lam = temporal_intensity(candidates, start, end, **intensity_kwargs) + keep = rng.uniform(0, lam_max, batch) < lam + accepted.append(candidates[keep]) + total += int(keep.sum()) + if total >= n: + break + else: + raise RuntimeError("Could not sample enough timestamps.") + stamps = accepted[0].append(accepted[1:]) if len(accepted) > 1 else accepted[0] + return stamps[:n].floor("s") + + +def generate_crimes( + study_area, + n_events=5000, + start="2019-01-01", + end="2020-12-31", + n_hotspots=3, + hotspot_share=0.7, + hotspot_sd_km=0.5, + tags=None, + trend=0.0, + annual_amplitude=0.2, + annual_peak_month=1, + weekly_profile=DEFAULT_WEEKLY_PROFILE, + hourly_profile=DEFAULT_HOURLY_PROFILE, + seed=None, + return_hotspots=False, +): + """ + Generate synthetic crime events inside a study area. + + Args: + study_area (GeoDataFrame): Boundary of the study area (any CRS). See + :func:`predspot.crime_mapping.load_study_area` to fetch one from + OpenStreetMap. + n_events (int): Number of events to generate. + start (str or Timestamp): First possible timestamp. + end (str or Timestamp): Last possible timestamp. + n_hotspots (int): Number of Gaussian hotspots (0 for a uniform map). + hotspot_share (float): Fraction of events that belong to hotspots; + the rest is uniform background (0-1). + hotspot_sd_km (float or sequence): Standard deviation of the hotspot + clouds in km (one value or one per hotspot). + tags (dict or sequence): Crime types. A dict maps type to relative + weight; a sequence gives equal weights. Defaults to + :data:`DEFAULT_TAGS`. + trend (float): Relative change of the event rate from ``start`` to + ``end`` (``0.5`` = +50%). + annual_amplitude (float): Amplitude of the annual cycle (0 disables). + annual_peak_month (int): Month (1-12) where the annual cycle peaks. + weekly_profile (sequence or None): 7 weights Monday..Sunday + (``None`` disables the weekly pattern). + hourly_profile (sequence or None): 24 weights, hour 0..23 + (``None`` disables the hour-of-day pattern). + seed (int, optional): Seed for reproducibility. + return_hotspots (bool): Also return the hotspot centres. + + Returns: + pandas.DataFrame: Events with ``tag``, ``t``, ``lon``, ``lat`` columns + sorted by time, ready for :class:`predspot.Dataset`. If + ``return_hotspots`` is True, a tuple ``(crimes, hotspots)`` where + ``hotspots`` is a GeoDataFrame with the centre, ``sd_km`` and + ``share`` of each hotspot. + """ + if n_events <= 0: + raise ValueError("n_events must be positive.") + if not 0 <= hotspot_share <= 1: + raise ValueError("hotspot_share must be between 0 and 1.") + if n_hotspots < 0: + raise ValueError("n_hotspots must be >= 0.") + rng = np.random.default_rng(seed) + polygon = _study_polygon(study_area) + start, end = pd.Timestamp(start), pd.Timestamp(end) + + # --- tags --------------------------------------------------------------- + if tags is None: + tags = DEFAULT_TAGS + if isinstance(tags, dict): + names, weights = list(tags.keys()), np.asarray(list(tags.values()), dtype=float) + else: + names, weights = list(tags), np.ones(len(tags)) + if len(names) == 0 or (weights < 0).any() or weights.sum() == 0: + raise ValueError("tags must contain at least one type with a positive weight.") + tag_values = rng.choice(names, size=n_events, p=weights / weights.sum()) + + # --- space -------------------------------------------------------------- + if n_hotspots == 0: + n_hot = 0 + else: + n_hot = int(round(n_events * hotspot_share)) + n_bg = n_events - n_hot + hot_lon, hot_lat = sample_points_in_polygon(polygon, n_hotspots, rng) + centers = np.column_stack([hot_lon, hot_lat]) if n_hotspots else np.empty((0, 2)) + shares = rng.dirichlet(np.full(n_hotspots, 2.0)) if n_hotspots else np.empty(0) + per_center = rng.multinomial(n_hot, shares) if n_hot else np.zeros(n_hotspots, dtype=int) + sd_km = np.broadcast_to(np.asarray(hotspot_sd_km, dtype=float), (n_hotspots,)) + lon_h, lat_h = sample_points_around(polygon, centers, sd_km, per_center, rng) + lon_b, lat_b = ( + sample_points_in_polygon(polygon, n_bg, rng) if n_bg else (np.empty(0), np.empty(0)) + ) + lon = np.concatenate([lon_h, lon_b]) + lat = np.concatenate([lat_h, lat_b]) + order = rng.permutation(n_events) + lon, lat = lon[order], lat[order] + + # --- time --------------------------------------------------------------- + timestamps = sample_timestamps( + n_events, + start, + end, + rng, + trend=trend, + annual_amplitude=annual_amplitude, + annual_peak_month=annual_peak_month, + weekly_profile=weekly_profile, + hourly_profile=hourly_profile, + ) + + crimes = ( + pd.DataFrame({"tag": tag_values, "t": timestamps, "lon": lon, "lat": lat}) + .sort_values("t") + .reset_index(drop=True) + ) + logger.debug("Generated %d synthetic events with %d hotspots", n_events, n_hotspots) + if not return_hotspots: + return crimes + hotspots = gpd.GeoDataFrame( + {"sd_km": sd_km, "share": shares * hotspot_share}, + geometry=gpd.points_from_xy(centers[:, 0], centers[:, 1]), + crs=WGS84, + ) + return crimes, hotspots diff --git a/tests/test_load_study_area.py b/tests/test_load_study_area.py new file mode 100644 index 0000000..70efc14 --- /dev/null +++ b/tests/test_load_study_area.py @@ -0,0 +1,72 @@ +import os +import sys +import types + +import geopandas as gpd +import pytest +from shapely.geometry import Point, box + +from predspot.crime_mapping import load_study_area + + +def _fake_osmnx(monkeypatch, geometry): + """Install a fake `osmnx` module whose geocode_to_gdf returns `geometry`.""" + calls = {} + + def geocode_to_gdf(query, which_result=None): + calls["query"], calls["which_result"] = query, which_result + return gpd.GeoDataFrame( + { + "osm_type": ["relation"], + "osm_id": [1], + "name": ["Natal"], + "display_name": ["Natal, Rio Grande do Norte, Brasil"], + "place_rank": [16], + }, + geometry=[geometry], + crs="EPSG:4326", + ) + + fake = types.ModuleType("osmnx") + fake.geocode_to_gdf = geocode_to_gdf + monkeypatch.setitem(sys.modules, "osmnx", fake) + return calls + + +def test_load_study_area_returns_polygon(monkeypatch): + calls = _fake_osmnx(monkeypatch, box(-35.3, -5.9, -35.2, -5.8)) + area = load_study_area("Natal, Brazil") + assert calls == {"query": "Natal, Brazil", "which_result": None} + assert list(area.columns) == ["name", "display_name", "osm_type", "osm_id", "geometry"] + assert area.crs.to_epsg() == 4326 + assert area.geometry.iloc[0].geom_type == "Polygon" + + +def test_load_study_area_reprojects(monkeypatch): + _fake_osmnx(monkeypatch, box(-35.3, -5.9, -35.2, -5.8)) + area = load_study_area("Natal, Brazil", crs="EPSG:31985") + assert area.crs.to_epsg() == 31985 + + +def test_load_study_area_rejects_points(monkeypatch): + _fake_osmnx(monkeypatch, Point(-35.2, -5.8)) + with pytest.raises(ValueError, match="non-polygon"): + load_study_area("Somewhere") + + +def test_load_study_area_without_osmnx(monkeypatch): + monkeypatch.setitem(sys.modules, "osmnx", None) + with pytest.raises(ImportError, match="predspot\\[osm\\]"): + load_study_area("Natal, Brazil") + + +@pytest.mark.skipif( + not os.environ.get("PREDSPOT_NETWORK_TESTS"), + reason="set PREDSPOT_NETWORK_TESTS=1 to query OpenStreetMap", +) +def test_load_study_area_network(): + pytest.importorskip("osmnx") + area = load_study_area("Natal, Rio Grande do Norte, Brazil") + 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 diff --git a/tests/test_synthetic.py b/tests/test_synthetic.py new file mode 100644 index 0000000..592f0a8 --- /dev/null +++ b/tests/test_synthetic.py @@ -0,0 +1,124 @@ +import geopandas as gpd +import numpy as np +import pandas as pd +import pytest +import shapely +from shapely.geometry import Polygon + +from predspot import Dataset, synthetic +from predspot.synthetic import generate_crimes + + +@pytest.fixture(scope="module") +def irregular_area(): + poly = Polygon( + [(-35.30, -5.90), (-35.18, -5.92), (-35.16, -5.85), (-35.20, -5.82), + (-35.19, -5.78), (-35.28, -5.77), (-35.32, -5.84)] + ) # fmt: skip + return gpd.GeoDataFrame(geometry=[poly], crs="EPSG:4326") + + +def test_generate_crimes_basic(irregular_area): + crimes, hotspots = generate_crimes( + irregular_area, n_events=3000, n_hotspots=3, seed=0, return_hotspots=True + ) + assert list(crimes.columns) == ["tag", "t", "lon", "lat"] + assert len(crimes) == 3000 + assert crimes["t"].is_monotonic_increasing + assert crimes["t"].between("2019-01-01", "2020-12-31").all() + polygon = irregular_area.geometry.iloc[0] + assert shapely.contains_xy(polygon, crimes["lon"].values, crimes["lat"].values).all() + assert len(hotspots) == 3 + assert hotspots.geometry.within(polygon).all() + assert np.isclose(hotspots["share"].sum(), 0.7) + # it plugs directly into Dataset + assert Dataset(crimes, irregular_area).shape["crimes"][0] == 3000 + + +def test_generate_crimes_is_reproducible(study_area): + a = generate_crimes(study_area, n_events=500, seed=123) + b = generate_crimes(study_area, n_events=500, seed=123) + pd.testing.assert_frame_equal(a, b) + c = generate_crimes(study_area, n_events=500, seed=124) + assert not a.equals(c) + + +def test_hotspots_concentrate_events(study_area): + crimes, hotspots = generate_crimes( + study_area, n_events=4000, n_hotspots=1, hotspot_share=0.8, hotspot_sd_km=0.3, + seed=1, return_hotspots=True, + ) # fmt: skip + center = hotspots.geometry.iloc[0] + d_lon = (crimes["lon"] - center.x) * 111.32 * np.cos(np.radians(center.y)) + d_lat = (crimes["lat"] - center.y) * 110.57 + within_1km = np.hypot(d_lon, d_lat) < 1.0 + # ~80% hotspot events within ~3 sd, plus a little background in a 100 km2 box + assert 0.7 < within_1km.mean() < 0.9 + + +def test_uniform_when_no_hotspots(study_area): + crimes = generate_crimes(study_area, n_events=4000, n_hotspots=0, seed=2) + west, south, east, north = study_area.total_bounds + # split the box into 4 quadrants; each should hold ~25% of the events + q = ((crimes["lon"] > (west + east) / 2).astype(int) * 2 + + (crimes["lat"] > (south + north) / 2).astype(int)) # fmt: skip + shares = q.value_counts(normalize=True) + assert len(shares) == 4 and (shares.between(0.2, 0.3)).all() + + +def test_temporal_patterns(study_area): + crimes = generate_crimes( + study_area, n_events=20000, trend=1.0, annual_amplitude=0.5, annual_peak_month=7, + weekly_profile=(1, 1, 1, 1, 1, 1, 6), hourly_profile=None, seed=3, + ) # fmt: skip + t = crimes["t"] + # weekly: Sunday gets 6/12 = 50% of the events + assert 0.45 < (t.dt.dayofweek == 6).mean() < 0.55 + # no hourly profile: hours roughly uniform + hours = t.dt.hour.value_counts(normalize=True) + assert hours.max() < 0.07 + # trend: the last year has more events than the first year + per_year = t.dt.year.value_counts() + assert per_year[2020] > 1.3 * per_year[2019] + # annual cycle peaks in July + assert t.dt.month.value_counts().idxmax() in (6, 7, 8) + + +def test_tags_weights(study_area): + crimes = generate_crimes(study_area, n_events=5000, tags={"theft": 3, "fraud": 1}, seed=4) + shares = crimes["tag"].value_counts(normalize=True) + assert set(shares.index) == {"theft", "fraud"} + assert 0.7 < shares["theft"] < 0.8 + listed = generate_crimes(study_area, n_events=200, tags=["a", "b", "c"], seed=4) + assert set(listed["tag"]) == {"a", "b", "c"} + + +def test_projected_study_area(study_area): + projected = study_area.to_crs(study_area.estimate_utm_crs()) + crimes = generate_crimes(projected, n_events=300, seed=5) + west, south, east, north = study_area.total_bounds + assert crimes["lon"].between(west, east).all() + assert crimes["lat"].between(south, north).all() + + +def test_validation(study_area): + with pytest.raises(ValueError): + generate_crimes(study_area, n_events=0) + with pytest.raises(ValueError): + generate_crimes(study_area, n_events=10, hotspot_share=1.5) + with pytest.raises(ValueError): + generate_crimes(study_area, n_events=10, tags={}) + with pytest.raises(ValueError, match="7 values"): + generate_crimes(study_area, n_events=10, weekly_profile=(1, 2)) + with pytest.raises(ValueError, match="24 values"): + generate_crimes(study_area, n_events=10, hourly_profile=(1, 2)) + with pytest.raises(ValueError, match="end must be after"): + generate_crimes(study_area, n_events=10, start="2020-01-01", end="2019-01-01") + + +def test_temporal_intensity_shape(): + start, end = pd.Timestamp("2019-01-01"), pd.Timestamp("2019-12-31") + stamps = pd.date_range(start, end, freq="D") + lam = synthetic.temporal_intensity(stamps, start, end, trend=1.0) + assert lam.shape == (len(stamps),) + assert np.isclose(lam[0], 1.0) and np.isclose(lam[-1], 2.0)