Upload 16 files

.gitattributes

@@ -33,3 +33,4 @@ saved_model/**/* filter=lfs diff=lfs merge=lfs -text
 *.zip filter=lfs diff=lfs merge=lfs -text
 *.zst filter=lfs diff=lfs merge=lfs -text
 *tfevents* filter=lfs diff=lfs merge=lfs -text
+Task[[:space:]]2[[:space:]]Boundaries[[:space:]]and[[:space:]]Grids.pdf filter=lfs diff=lfs merge=lfs -text

README.md

@@ -0,0 +1,173 @@
+# Task 6: Incomplete Reporting Removal-Experiment Scaffold
+
+This directory contains the simulation scaffold for **Task 6**, designed to evaluate the robustness of crime hotspot identification under incomplete crime reporting / underreporting.
+
+---
+
+## 1. Python Environment Status
+
+The committed `.venv` in this repository was created on another computer with a hardcoded path to:
+```
+C:\Users\padma\Anaconda3\python.exe
+```
+When invoked on this workstation, `.venv\Scripts\python.exe` fails with:
+```
+No Python at '"C:\Users\padma\Anaconda3\python.exe'
+```
+
+### Available & Verified Python Installation:
+- **Executable**: `py` (Windows Python Launcher) -> `Python 3.14.0` at:
+  `C:\Users\mukun\AppData\Local\Python\pythoncore-3.14-64\python.exe`
+- **Installed Geospatial Stack**:
+  - `geopandas 1.2.0`
+  - `shapely 2.2.0`
+  - `pyproj 3.8.0`
+  - `pyogrio 0.13.0`
+  - `pandas 3.0.6`
+  - `numpy 2.5.3`
+
+All commands and tests should be run using `py` (or by recreating a local virtual environment pointing to the available Python installation).
+
+---
+
+## 2. Running the Synthetic Test Suite
+
+The test suite validates all removal mechanics, coordinate math, validations, and edge cases using synthetic grids and points. No real crime data is required or modified.
+
+### Run Tests:
+```powershell
+py test_task6_experiment.py -v
+```
+or via standard `unittest`:
+```powershell
+py -m unittest test_task6_experiment.py -v
+```
+
+### Test Coverage (12 Comprehensive Tests):
+1. `test_zero_removal`: Verifies that a removal probability $p = 0.0$ retains 100% of crimes, produces identical recounts before/after, and retains 100% of original top cells.
+2. `test_complete_removal`: Verifies that $p = 1.0$ removes 100% of crimes, sets all cell counts to 0, and deterministically tie-breaks by ascending `cell_id`.
+3. `test_uneven_removal_changes_ranking`: Verifies that cell-specific removal with differential rates alters the hotspot ranking and reduces the retention fraction.
+4. `test_boundary_assignment`: Verifies boundary coordinate assignment using the lower-left half-open interval rule $[x_{\min}, x_{\min} + 500) \times [y_{\min}, y_{\min} + 500)$, ensuring exact vertices and edges are assigned to exactly one cell.
+5. `test_invalid_rates`: Verifies rejection of negative rates, rates $> 1.0$, `NaN`, `Inf`, non-numeric types, incomplete cell rate mappings, and extraneous cell IDs.
+6. `test_repeatability`: Verifies that identical random seeds produce identical simulations and results, whereas different seeds generate distinct draws.
+7. `test_crs_validation`: Verifies strict enforcement of `EPSG:27700` and rejection of missing CRS or non-27700 layers.
+8. `test_points_outside_city_grid`: Verifies that crime points falling outside the city grid are rejected by default with `ValueError` (or optionally dropped if `allow_outside=True`).
+9. `test_grid_id_validation`: Verifies unique `cell_id` values, correct format `"<x_min>_<y_min>"`, and alignment to 500 m grid multiples.
+10. `test_recount_includes_zero_crime_cells`: Verifies that every cell in the grid is preserved in the recount before and after removal, with count 0 for cells without incidents.
+11. `test_top_cells_selection_and_deterministic_ties`: Verifies $k = \lceil 0.10 \times N \rceil$ calculation and deterministic secondary sorting on `cell_id` (ascending).
+12. `test_partially_clipped_edge_cell_validation`: Verifies that a point located inside the full 500 m unclipped square of an edge cell but outside its clipped city geometry is rejected with `ValueError` when `allow_outside=False`, and is counted and dropped when `allow_outside=True`.
+
+---
+
+## 3. How to Use the Module with Task 3 Leeds Crime Layer
+
+Once Task 3 outputs the city-assigned crime layer (e.g. `data/interim/crime/leeds_crime_EPSG27700.gpkg`), the experiment can be executed either programmatically or via CLI.
+
+### Option A: Command-Line Interface (CLI)
+
+#### 1. Uniform Removal Rate (e.g., 20% underreporting):
+```powershell
+py task6_experiment.py `
+  --crime data/interim/crime/leeds_crime_EPSG27700.gpkg `
+  --grid data/interim/grids/leeds_grid500_EPSG27700.gpkg `
+  --uniform-rate 0.20 `
+  --seed 42 `
+  --output-json outputs/task6_leeds_uniform_0.20.json
+```
+
+#### 2. Cell-Specific Removal Rates:
+The module strictly requires a removal rate for **every grid cell** (e.g., all 2,390 cells in the Leeds 500 m grid). Incomplete mappings or extra cell IDs are rejected by the validator.
+
+You can generate a complete rates JSON or CSV directly from the grid file:
+```python
+import json, geopandas as gpd
+
+grid = gpd.read_file("data/interim/grids/leeds_grid500_EPSG27700.gpkg")
+# Example: 40% removal for edge cells, 15% removal for interior full cells
+cell_rates = {
+    cid: (0.40 if is_edge else 0.15)
+    for cid, is_edge in zip(grid["cell_id"], grid["is_edge"])
+}
+with open("leeds_cell_rates.json", "w") as f:
+    json.dump(cell_rates, f, indent=2)
+```
+
+The resulting `leeds_cell_rates.json` covers all 2,390 cells:
+```json
+{
+  "428500_422500": 0.40,
+  "429000_422500": 0.40,
+  "429500_422500": 0.40,
+  ...
+  "444500_446000": 0.15
+}
+```
+
+Run with the rates file:
+```powershell
+py task6_experiment.py `
+  --crime data/interim/crime/leeds_crime_EPSG27700.gpkg `
+  --grid data/interim/grids/leeds_grid500_EPSG27700.gpkg `
+  --cell-rates leeds_cell_rates.json `
+  --seed 42
+```
+
+---
+
+### Option B: Python API
+
+```python
+import geopandas as gpd
+from task6_experiment import run_removal_experiment
+
+# 1. Load layers
+crime_layer = "data/interim/crime/leeds_crime_EPSG27700.gpkg"
+grid_layer = "data/interim/grids/leeds_grid500_EPSG27700.gpkg"
+
+# 2. Run uniform removal experiment (30% underreporting)
+result_uniform = run_removal_experiment(
+    points=crime_layer,
+    grid=grid_layer,
+    removal=0.30,
+    seed=42,
+    top_fraction=0.10,  # evaluates top ceil(10% of cells)
+)
+
+# 3. View summary
+result_uniform.print_summary()
+print(f"Top 10% hotspot cells retained: {result_uniform.fraction_retained:.2%}")
+print(f"Retained top cells: {result_uniform.retained_top_cells}")
+print(f"Lost top cells: {result_uniform.lost_top_cells}")
+
+# 4. Access full cell recount Series
+counts_before = result_uniform.counts_before  # pd.Series indexed by cell_id
+counts_after = result_uniform.counts_after    # pd.Series indexed by cell_id
+
+# 5. Run cell-specific removal experiment
+grid = gpd.read_file(grid_layer)
+# Example: Vary removal probability by edge vs interior cells
+cell_rates = {
+    cid: (0.40 if is_edge else 0.15)
+    for cid, is_edge in zip(grid["cell_id"], grid["is_edge"])
+}
+
+result_cell = run_removal_experiment(
+    points=crime_layer,
+    grid=grid_layer,
+    removal=cell_rates,
+    seed=42,
+)
+result_cell.print_summary()
+```
+
+---
+
+## 4. Key Implementation Rules & Conventions
+
+- **Working CRS**: `EPSG:27700` (British National Grid, metres). Both crime points and grid must be in `EPSG:27700`.
+- **Lower-Left Assignment Rule**:
+  $$\text{x\_min} = \lfloor x / 500 \rfloor \times 500, \quad \text{y\_min} = \lfloor y / 500 \rfloor \times 500$$
+  $$\text{cell\_id} = \text{"<x\_min>\_<y\_min>"}$$
+  This half-open box rule $[x_{\min}, x_{\min} + 500) \times [y_{\min}, y_{\min} + 500)$ ensures every point is mapped to exactly one cell.
+- **Hotspot Selection**: Top $k = \lceil 0.10 \times N_{\text{cells}} \rceil$.
+- **Tie-Breaking Rule**: When crime counts are tied, cells are ordered deterministically by ascending `cell_id` (alphabetical order of `"<x_min>_<y_min>"`).

Task 2 Boundaries and Grids.pdf

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:55ab7fe83b91be67ea0ce08baec5d16ae52908fbf222d15ee630d9918a3e6c68
+size 840404

__init__.py

@@ -0,0 +1 @@
+"""Reusable project code. Task 2: geomutils, boundaries, grid."""

boundaries.py

@@ -0,0 +1,133 @@
+"""Task 2a: read an ONS Local Authority District (LAD) file and cut out one boundary per city."""
+import re
+from pathlib import Path
+
+import geopandas as gpd
+import pyogrio
+import shapely
+from pyproj import CRS
+
+from src.geomutils import crs_tag, polygonal_part
+
+# ONS boundary resolution codes, as they appear in ONS file and layer names.
+RESOLUTIONS = {
+    "BFC": "full resolution, clipped to the coastline",
+    "BFE": "full resolution, extent of the realm",
+    "BGC": "generalised (20 m), clipped to the coastline",
+    "BSC": "super generalised (200 m), clipped to the coastline",
+    "BUC": "ultra generalised (500 m), clipped to the coastline",
+}
+MONTHS = ("JAN", "FEB", "MAR", "APR", "MAY", "JUN", "JUL", "AUG", "SEP", "OCT", "NOV", "DEC")
+BOUNDARY_COLUMNS = ["city", "lad_code", "lad_name", "area_m2"]
+
+
+def boundary_filename(city, crs):
+    """leeds, EPSG:27700 -> 'leeds_boundary_EPSG27700.gpkg'."""
+    return f"{city}_boundary_{crs_tag(crs)}.gpkg"
+
+
+def describe_lad_file(path):
+    """Print and return what an ONS LAD file contains and which version it is.
+
+    The version comes from the file and layer names (month, year, resolution code
+    such as BFC) and from the code column name (LAD25CD means the 2025 code set).
+    """
+    path = Path(path)
+    layers = [name for name, _ in pyogrio.list_layers(path)]
+    if len(layers) != 1:
+        raise ValueError(f"{path.name}: expected exactly one layer, found {layers}")
+    info = pyogrio.read_info(path, layer=layers[0])
+    columns = list(info["fields"])
+    tokens = set(re.split(r"[^A-Z0-9]+", f"{path.stem}_{layers[0]}".upper()))
+    desc = {
+        "file": path.name,
+        "layer": layers[0],
+        "columns": columns,
+        "crs": info["crs"],
+        "geometry_type": info["geometry_type"],
+        "features": int(info["features"]),
+        "resolution": sorted(tokens & RESOLUTIONS.keys()),
+        "month": sorted(tokens & set(MONTHS)),
+        "year": sorted(t for t in tokens if re.fullmatch(r"(19|20)\d\d", t)),
+        "code_fields": [c for c in columns if re.fullmatch(r"LAD\d\dCD", c)],
+        "name_fields": [c for c in columns if re.fullmatch(r"LAD\d\dNM", c)],
+    }
+    print(f"file:          {desc['file']}")
+    print(f"layers:        {layers}")
+    print(f"columns:       {columns}")
+    print(f"CRS:           {desc['crs']}")
+    print(f"geometry type: {desc['geometry_type']}")
+    print(f"features:      {desc['features']}")
+    print(f"version:       month={desc['month']} year={desc['year']} "
+          f"resolution={desc['resolution']} code column={desc['code_fields']}")
+    return desc
+
+
+def require_bfc(desc):
+    """Stop unless the file is unambiguously the full resolution clipped (BFC) version."""
+    if desc["resolution"] != ["BFC"]:
+        found = ", ".join(f"{r} ({RESOLUTIONS[r]})" for r in desc["resolution"]) or "no resolution code"
+        raise ValueError(f"{desc['file']} is not the BFC version: found {found}")
+
+
+def read_lad_features(path, layer, code_field, codes):
+    """Read only the LAD features whose code is in `codes`, sorted by code."""
+    for code in codes:
+        if not re.fullmatch(r"[EWSN]\d{8}", code):
+            raise ValueError(f"not a valid GSS LAD code: {code!r}")
+    where = f"{code_field} IN ({', '.join(repr(c) for c in codes)})"
+    lad = gpd.read_file(path, layer=layer, where=where, engine="pyogrio")
+    return lad.sort_values(code_field, kind="stable").reset_index(drop=True)
+
+
+def to_target_crs(gdf, target_crs):
+    """Return (gdf in target_crs, reprojected?) using a real to_crs reprojection.
+
+    Never relabels the CRS. A layer with no CRS is an error, because its
+    coordinates cannot be interpreted.
+    """
+    if gdf.crs is None:
+        raise ValueError("boundary layer has no CRS defined")
+    if CRS.from_user_input(gdf.crs).equals(CRS.from_user_input(target_crs)):
+        return gdf, False
+    return gdf.to_crs(target_crs), True
+
+
+def city_boundary(lad, city, lad_code, expected_name, code_field, name_field):
+    """Cut one city out of the LAD layer as a single valid (Multi)Polygon.
+
+    Several features with the same code are dissolved into one. An invalid
+    geometry is repaired with shapely.make_valid and only its polygon parts kept.
+    Returns (GeoDataFrame with BOUNDARY_COLUMNS + geometry, report dict).
+    """
+    rows = lad[lad[code_field] == lad_code]
+    if rows.empty:
+        raise ValueError(f"{city}: LAD code {lad_code} not found in the boundary file")
+    names = sorted(set(rows[name_field]))
+    if names != [expected_name]:
+        raise ValueError(f"{city}: {lad_code} is named {names} in the file, expected '{expected_name}'")
+
+    geom = rows.geometry.union_all() if len(rows) > 1 else rows.geometry.iloc[0]
+    was_valid = bool(shapely.is_valid(geom))
+    invalid_reason = None if was_valid else shapely.is_valid_reason(geom)
+    if not was_valid:
+        geom = shapely.make_valid(geom)
+    geom = polygonal_part(geom)
+    if geom.is_empty or not geom.is_valid:
+        raise ValueError(f"{city}: boundary is empty or still invalid after repair")
+
+    out = gpd.GeoDataFrame(
+        {"city": [city], "lad_code": [lad_code], "lad_name": [names[0]], "area_m2": [geom.area]},
+        geometry=[geom],
+        crs=lad.crs,
+    )
+    report = {
+        "city": city,
+        "n_features": len(rows),
+        "dissolved": len(rows) > 1,
+        "was_valid": was_valid,
+        "invalid_reason": invalid_reason,
+        "repaired": not was_valid,
+        "n_polygons": len(geom.geoms),
+    }
+    return out, report

config.json

@@ -0,0 +1,27 @@
+{
+  "project_year": "2025",
+  "forces": {
+    "west-yorkshire": "Leeds",
+    "south-yorkshire": "Sheffield",
+    "northumbria": "Newcastle upon Tyne",
+    "nottinghamshire": "Nottingham"
+  },
+  "boundary_type": "Local Authority District (BFC)",
+  "primary_crime_category": "Burglary",
+  "crs": "EPSG:27700",
+  "grid_size_m": 500,
+  "data_dir": ".",
+  "osm_junction_tolerance_m": 12,
+  "cities": {
+    "leeds": {"force": "west-yorkshire", "lad_code": "E08000035"},
+    "sheffield": {"force": "south-yorkshire", "lad_code": "E08000039"},
+    "newcastle": {"force": "northumbria", "lad_code": "E08000021"},
+    "nottingham": {"force": "nottinghamshire", "lad_code": "E06000018"}
+  },
+  "boundary_file": "data/raw/boundaries/LAD_DEC_2025_UK_BFC.gpkg",
+  "boundary_code_field": "LAD25CD",
+  "boundary_name_field": "LAD25NM",
+  "boundaries_dir": "data/interim/boundaries",
+  "grids_dir": "data/interim/grids",
+  "qa_dir": "outputs/qa"
+}

conventions.md

@@ -0,0 +1,120 @@
+# Project Conventions
+
+## General Settings
+- **Year**: 2025 (January to December)
+- **Police Forces**: West Yorkshire (Leeds), South Yorkshire (Sheffield), Northumbria (Newcastle), Nottinghamshire (Nottingham)
+- **Boundary Type**: Local Authority District (LAD) using ONS full resolution clipped boundaries (BFC).
+- **Working CRS**: EPSG:27700 (British National Grid), units in metres.
+- **Grid**: 500m x 500m. The origin is strictly snapped down to a multiple of 500 in BNG coordinates, ensuring consistent grid alignment across cities and years.
+- **Grid Origin Coordinates per City**: [To be populated by Task 2]
+- **Data Location**: the project folder, i.e. the folder that holds `config.json` (`data_dir` = `.`). All paths are relative to it: Task 1 crime layers in the project folder itself, raw boundary file in `data/raw/boundaries/`, city boundaries in `data/interim/boundaries/`, grids in `data/interim/grids/`, QA outputs in `outputs/qa/`.
+- **File Naming Rule**: Intermediate outputs are saved as GeoPackages (`.gpkg`) or Parquet. Every spatial layer filename includes the CRS (e.g., `crime_west-yorkshire_EPSG27700.gpkg`).
+
+## Crime Data Rules
+- **Exact Crime Categories Used**: [To be confirmed after Task 1 exploration, placeholder: "Burglary"]
+- **Anti-Social Behaviour**: "Anti-social behaviour" is excluded from all-crime totals but reported separately in the data quality table.
+- **File Filters**: Only merge files ending in `-street.csv`. Outcome and stop-and-search files are explicitly ignored.
+- **Rows Without Coordinates**: Rows with no coordinates are counted and reported in the data quality table, then excluded from all spatial analysis. 2025 counts: West Yorkshire 2,162, South Yorkshire 8,069, Northumbria 0, Nottinghamshire 806.
+
+## OSM Data Rules
+- **Exact OSM Tag Filters Used**: [To be recorded in Task 4 after inspecting raw columns]
+- **Roads Kept**: Kept only roads used by vehicles and general traffic. (Excluded: footway, path, cycleway, steps, track, service roads).
+- **Junction Extraction**: Computed by building a graph from shared vertices of road lines. A true junction is defined as a node with 3 or more connections. Nodes within a **12 metre** tolerance are merged into a single point to handle complex intersections like roundabouts and dual carriageways.
+- **POI Deduplication**: Point POIs that fall completely inside a polygon POI of the exact same category are dropped. Remaining polygons are then converted to a single point using `representative_point()`, ensuring no double counting.
+
+## Task 2: City Boundaries and Analysis Grid (completed 2026-10-09)
+- **Boundary file**: `data/raw/boundaries/LAD_DEC_2025_UK_BFC.gpkg`. This is ONS *Local Authority Districts (December 2025) Boundaries UK BFC*, the full-resolution version clipped to the coastline. It has one layer (`LAD_DEC_2025_UK_BFC`) of 361 MultiPolygon features in EPSG:27700. The code column is `LAD25CD` (2025 code set) and the name column is `LAD25NM`. It was downloaded on 2026-10-09 from the ONS Open Geography Portal (item `92150c7aa60540c5814abe3b26bce6d0`). SHA-256: `29ba41a24d7943db1e8ab50c73506cb486c555f8c0e9bc945b24aece068d3602`. The May 2025 (V2) BFC release has identical areas for all four cities (they differ by less than 1 m²).
+- **City keys and LAD codes** (stored in `config.json` under `cities`):
+
+  | City key | LAD code | LAD name in file | Area (km²) |
+  |---|---|---|---|
+  | leeds | E08000035 | Leeds | 551.70 |
+  | sheffield | E08000039 | Sheffield | 368.00 |
+  | newcastle | E08000021 | Newcastle upon Tyne | 113.45 |
+  | nottingham | E06000018 | Nottingham | 74.64 |
+
+  - Sheffield LAD code is E08000039 (code in use since 2023). The older code E08000019 must not be used.
+  - Newcastle-under-Lyme (E07000195) is in the file and was not selected.
+- **Reprojection**: not needed, because the file is already EPSG:27700.
+- **Geometry repair**: not needed. Each city is one valid feature with one polygon, so no dissolve was needed either. The code would dissolve multiple features and repair invalid geometry with `shapely.make_valid` if a future file needed it.
+- **Boundary outputs**: `data/interim/boundaries/<city>_boundary_EPSG27700.gpkg`, with columns `city, lad_code, lad_name, area_m2, geometry`.
+- **Grid rule** (`src/grid.py`, `make_grid(boundary_geom, cell_size, city, crs=...)`):
+  - **Origin**: snapped down to a multiple of the cell size, using integer arithmetic only:
+    - `x_start = floor(minx / size) * size`, `y_start = floor(miny / size) * size`
+    - The far edges are snapped up with `ceil`.
+  - **Selection**: square cells covering that box are kept if they intersect the boundary. The test uses an STRtree spatial index.
+  - **Clipping**:
+    - Cells fully inside the boundary keep their exact square.
+    - The other cells are clipped to the boundary, and only the polygon parts are kept.
+    - Cells whose clipped area is exactly 0 are dropped. At 500 m no cells were dropped in any city.
+  - **Edge cells and slivers**:
+    - Partial edge cells, including tiny slivers, are kept.
+    - `is_edge = area_share < 0.9999`, where `area_share = area_m2 / size²`.
+  - **Order**: rows run south to north, then west to east.
+- **cell_id rule**:
+  - The format is `"<x_min>_<y_min>"`: the lower-left corner of the full, unclipped cell, in whole metres in EPSG:27700. Example: `413000_422500`.
+  - `x_min` and `y_min` hold the same corner as integers.
+  - IDs are unique within a city and stable across reruns. A grid built from the same boundary and cell size always gives the same cells and IDs.
+  - A cell_id only identifies a cell together with its `cell_size_m`, because a 250 m grid reuses some of the same corner values.
+- **Grid outputs**: `data/interim/grids/<city>_grid<size>_EPSG27700.gpkg`, with columns `cell_id, x_min, y_min, cell_size_m, area_m2, area_share, is_edge, city, is_sliver, geometry` (MultiPolygon).
+- **Grid origin and cell counts per city, 500 m** (this fills in "Grid Origin Coordinates per City" under General Settings):
+
+  | City | Origin (x_start, y_start) | Cells | Full | Edge | Edge cells <5% | <10% | <25% |
+  |---|---|---|---|---|---|---|---|
+  | leeds | (413000, 422500) | 2,390 | 2,034 | 356 | 52 | 79 | 114 |
+  | sheffield | (413000, 378500) | 1,626 | 1,326 | 300 | 39 | 57 | 102 |
+  | newcastle | (414000, 562500) | 526 | 384 | 142 | 22 | 31 | 47 |
+  | nottingham | (450500, 332500) | 366 | 234 | 132 | 22 | 25 | 41 |
+
+  - The sum of cell areas equals the boundary area in every city (difference at most 5e-13 %).
+  - The smallest slivers are 0.16 m² (Leeds `419500_439000`) and 0.11 m² (Sheffield `433000_379500`). Density per area will be extreme for these cells; the later analysis must decide how to handle them.
+- **Sliver rule (PROPOSED TEAM DECISION, not yet agreed)**: Edge cells with area_share below 0.10 are flagged as slivers. They keep their counts, so city totals stay complete. They are excluded from density maps, density rankings and the concentration table. The number of crimes falling in slivers is reported per city.
+  - The grid files carry this flag as the `is_sliver` column (`area_share < 0.10`), added 2026-10-09. Slivers at 500 m: leeds 79, sheffield 57, newcastle 31, nottingham 25 (the same as the "<10%" column above).
+- **Reproducibility**:
+  - Software: Python 3.11.14, geopandas 1.2.0, shapely 2.2.0 (GEOS 3.14.1), pyproj 3.7.2 (PROJ 9.5.1), pyogrio 0.13.0 (GDAL 3.12.4). The full pinned list is in `requirements.txt`.
+  - Each grid has a SHA-256 fingerprint over its rows and geometry, stored in `outputs/qa/task2_grid500_checks.csv`. Two separate runs gave identical fingerprints.
+  - To rebuild: `python task2_boundaries_grids.py`. For the 250 m rerun: `python task2_boundaries_grids.py --cell-size 250`.
+
+### Task 2 reproduced in QGIS
+- **Software**: QGIS 3.44.12 'Solothurn' (LTR) with bundled Python 3.12.13, GDAL 3.13.1, GEOS 3.14.1 and PROJ 9.8.1.
+- **Script**: `task2_qgis.py`. Run it in the QGIS Python Console (Show Editor > Open Script > Run). Set `PROJECT_ROOT` at the top. All settings are read from `config.json`.
+- **Algorithms, in order, per city** (each one a `processing.run()` call to a native QGIS algorithm; the full parameters of every call are in `logs/task2_qgis_log.txt`):
+  1. `native:extractbyattribute`: LAD code from `config.json`, operator `=`.
+  2. `native:checkvalidity`: GEOS method.
+  3. `native:fixgeometries`: only if a geometry is invalid. Not needed for any city.
+  4. `native:dissolve`: one feature per city.
+  5. `native:refactorfields`: boundary schema `city, lad_code, lad_name, area_m2`. `area_m2` uses the planar `area($geometry)`, not `$area`.
+  6. `native:creategrid`: rectangles, 500 m, extent snapped outward to multiples of 500 m.
+  7. `native:clip`: grid clipped by the city boundary.
+  8. `native:extractbyexpression`: `area($geometry) > 0`.
+  9. `native:refactorfields`: the same grid schema as the Python files (names, types and order):
+     - `cell_id = concat(to_int(round("left")), '_', to_int(round("bottom")))`
+     - `x_min`, `y_min`: the same corner as integers
+     - `area_share = area($geometry) / 250000`
+     - `is_edge`: `area_share < 0.9999`
+     - `is_sliver`: `area_share < 0.10`
+  - `native:reprojectlayer` exists in the script but was not needed, because the LAD file is already EPSG:27700.
+- **Comparison with the Python outputs** (`outputs/qa/qgis/task2_qgis_vs_python.csv`; QGIS Python Console run, 2026-10-09 19:27, identical to two earlier headless test runs):
+
+  | City | QGIS cells | Python cells | cell_id sets identical | Max per-cell area difference (m²) | Boundary area difference (m²) | Result |
+  |---|---|---|---|---|---|---|
+  | leeds | 2,390 | 2,390 | yes | 9.6e-08 | 2.8e-05 | match |
+  | sheffield | 1,626 | 1,626 | yes | 4.0e-08 | 1.9e-05 | match |
+  | newcastle | 526 | 526 | yes | 5.2e-08 | 8.0e-07 | match |
+  | nottingham | 366 | 366 | yes | 7.4e-08 | 2.9e-06 | match |
+
+  - Grid origins, edge counts, sliver counts and the `is_edge`/`is_sliver` flags of every cell are identical to Python.
+  - No zero-area cells were dropped.
+  - The Python files were checksum-verified as unchanged.
+- **Where the outputs are**:
+  - Boundaries: `data/interim/qgis/boundaries/<city>_boundary_EPSG27700_qgis.gpkg`
+  - Grids: `data/interim/qgis/grids/<city>_grid500_EPSG27700_qgis.gpkg`
+  - Every intermediate: `data/interim/qgis/steps/<city>/`
+  - QA maps from print layouts: `outputs/qa/qgis/<city>_grid_qgis.png` (200 dpi)
+  - Log: `logs/task2_qgis_log.txt`
+- **Team project**: `QGIS/crimegeo_vgi.qgz`.
+  - Project CRS EPSG:27700, saved with relative paths.
+  - One layer group per city.
+  - OpenStreetMap XYZ basemap at the bottom of the layer list.
+  - Rule-based grid style: Full cell / Edge cell / Sliver, under 10% of a cell.
+  - Four editable print layouts named "<City> grid QA".

geomutils.py

@@ -0,0 +1,53 @@
+"""Small geometry and file helpers shared by the Task 2 modules."""
+import hashlib
+from pathlib import Path
+
+from shapely.geometry import GeometryCollection, MultiPolygon, Polygon
+
+
+def crs_tag(crs):
+    """'EPSG:27700' -> 'EPSG27700', the CRS label used in every spatial filename."""
+    return str(crs).replace(":", "")
+
+
+def polygonal_part(geom):
+    """Keep only the polygon parts of a geometry, always returned as a MultiPolygon.
+
+    Lines and points (for example where a grid cell only touches the boundary)
+    are discarded, so a geometry with no polygon parts becomes an empty MultiPolygon.
+    """
+    return MultiPolygon(_polygons(geom))
+
+
+def _polygons(geom):
+    if geom is None or geom.is_empty:
+        return []
+    if isinstance(geom, Polygon):
+        return [geom]
+    if isinstance(geom, (MultiPolygon, GeometryCollection)):
+        return [p for part in geom.geoms for p in _polygons(part)]
+    return []
+
+
+def write_gpkg(gdf, path):
+    """Write gdf as the only layer (named after the file) of a fresh GeoPackage."""
+    path = Path(path)
+    path.parent.mkdir(parents=True, exist_ok=True)
+    if path.exists():
+        path.unlink()
+    gdf.to_file(path, layer=path.stem, driver="GPKG", engine="pyogrio")
+    return path
+
+
+def fingerprint(gdf, columns):
+    """SHA-256 over the listed columns and the geometry WKB, row by row, in row order.
+
+    Identical layers give identical digests, so a teammate can confirm that a rerun
+    (or a copy of the file) matches the shared output exactly.
+    """
+    h = hashlib.sha256()
+    wkbs = gdf.geometry.to_wkb(hex=False, byte_order=1)
+    for values, wkb in zip(gdf[columns].itertuples(index=False, name=None), wkbs):
+        h.update("|".join(str(v) for v in values).encode())
+        h.update(wkb)
+    return h.hexdigest()

grid.py

@@ -0,0 +1,222 @@
+"""Task 2b: square analysis grids clipped to a city boundary.
+
+Cell corners are integer multiples of the cell size in the projected CRS (metres),
+so grids for different cities, years or runs line up, and the same boundary and
+cell size always give the same cells with the same cell_ids.
+"""
+import math
+
+import geopandas as gpd
+import numpy as np
+import shapely
+from pyproj import CRS
+
+from src.geomutils import crs_tag, polygonal_part
+
+EDGE_SHARE_THRESHOLD = 0.9999       # area_share below this -> is_edge
+SLIVER_SHARE_THRESHOLD = 0.10       # area_share below this -> is_sliver (proposed team decision)
+SLIVER_THRESHOLDS = (0.05, 0.10, 0.25)
+MAX_AREA_DIFF_PCT = 0.01            # sum of cell areas vs boundary area
+GRID_COLUMNS = ["cell_id", "x_min", "y_min", "cell_size_m", "area_m2", "area_share", "is_edge", "city", "is_sliver"]
+
+
+def grid_filename(city, cell_size, crs):
+    """leeds, 500, EPSG:27700 -> 'leeds_grid500_EPSG27700.gpkg'."""
+    return f"{city}_grid{cell_size}_{crs_tag(crs)}.gpkg"
+
+
+def _as_cell_size(cell_size):
+    if isinstance(cell_size, bool) or cell_size <= 0 or int(cell_size) != cell_size:
+        raise ValueError(f"cell_size must be a positive whole number of metres, got {cell_size!r}")
+    return int(cell_size)
+
+
+def grid_extent(bounds, cell_size):
+    """Snap a bounding box outwards to multiples of cell_size, as integers.
+
+    x_start = floor(minx / cell_size) * cell_size, x_end = ceil(maxx / cell_size) * cell_size,
+    and the same for y. Computed as floor(minx) // cell_size * cell_size, which is the
+    same value for a whole-number cell size but uses only integer arithmetic.
+    """
+    cell_size = _as_cell_size(cell_size)
+    minx, miny, maxx, maxy = bounds
+    x_start = math.floor(minx) // cell_size * cell_size
+    y_start = math.floor(miny) // cell_size * cell_size
+    x_end = -(-math.ceil(maxx) // cell_size) * cell_size
+    y_end = -(-math.ceil(maxy) // cell_size) * cell_size
+    return x_start, y_start, x_end, y_end
+
+
+def make_grid(boundary_geom, cell_size, city, *, crs, verbose=True):
+    """Build a square grid over boundary_geom, clipped to it.
+
+    boundary_geom: valid (Multi)Polygon in `crs`, a projected CRS in metres.
+    Returns a GeoDataFrame with GRID_COLUMNS + geometry (MultiPolygon), one row per
+    cell with a clipped area > 0, ordered south to north, then west to east.
+    cell_id, x_min and y_min refer to the lower-left corner of the full, unclipped cell.
+    Grid origin and counts are stored in .attrs.
+    """
+    cell_size = _as_cell_size(cell_size)
+    crs = CRS.from_user_input(crs)
+    if not crs.is_projected or crs.axis_info[0].unit_name not in ("metre", "meter"):
+        raise ValueError(f"make_grid needs a projected CRS in metres, got {crs.to_string()}")
+    if boundary_geom is None or boundary_geom.is_empty or not boundary_geom.is_valid:
+        raise ValueError(f"{city}: boundary geometry is empty or invalid")
+
+    x_start, y_start, x_end, y_end = grid_extent(boundary_geom.bounds, cell_size)
+    xs = np.arange(x_start, x_end, cell_size, dtype=np.int64)
+    ys = np.arange(y_start, y_end, cell_size, dtype=np.int64)
+    x_min, y_min = (a.ravel() for a in np.meshgrid(xs, ys))
+    cells = shapely.box(x_min, y_min, x_min + cell_size, y_min + cell_size)
+
+    # Spatial index: only cells that intersect the boundary. Sorting the hits keeps
+    # the south-to-north, west-to-east order no matter how the tree returns them.
+    shapely.prepare(boundary_geom)
+    hit = np.sort(shapely.STRtree(cells).query(boundary_geom, predicate="intersects"))
+    x_min, y_min, cells = x_min[hit], y_min[hit], cells[hit]
+
+    # Cells fully inside keep their exact square; the rest are clipped to the boundary.
+    inside = shapely.covers(boundary_geom, cells)
+    clipped = cells.copy()
+    clipped[~inside] = shapely.intersection(cells[~inside], boundary_geom)
+    geoms = np.array([polygonal_part(g) for g in clipped], dtype=object)
+    area = shapely.area(geoms)
+
+    keep = area > 0
+    n_dropped = int((~keep).sum())
+    x_min, y_min, geoms, area = x_min[keep], y_min[keep], geoms[keep], area[keep]
+    area_share = area / cell_size**2
+
+    grid = gpd.GeoDataFrame(
+        {
+            "cell_id": [f"{x}_{y}" for x, y in zip(x_min.tolist(), y_min.tolist())],
+            "x_min": x_min,
+            "y_min": y_min,
+            "cell_size_m": np.full(len(area), cell_size, dtype=np.int64),
+            "area_m2": area,
+            "area_share": area_share,
+            "is_edge": area_share < EDGE_SHARE_THRESHOLD,
+            "city": city,
+            "is_sliver": area_share < SLIVER_SHARE_THRESHOLD,
+        },
+        geometry=list(geoms),
+        crs=crs,
+    )
+    grid.attrs.update(
+        x_start=x_start, y_start=y_start, x_end=x_end, y_end=y_end,
+        n_bbox_cells=len(xs) * len(ys), n_intersecting=len(hit),
+        n_inside=int(inside.sum()), n_clipped=int((~inside).sum()), n_dropped_zero_area=n_dropped,
+    )
+    if verbose:
+        print(f"[{city}] {cell_size} m grid: snapped box ({x_start}, {y_start}) to ({x_end}, {y_end}), "
+              f"{len(xs)} x {len(ys)} = {len(xs) * len(ys):,} candidate cells")
+        print(f"[{city}] intersecting the boundary (STRtree): {len(hit):,} "
+              f"(fully inside {int(inside.sum()):,}, clipped {int((~inside).sum()):,})")
+        print(f"[{city}] dropped, clipped area exactly 0 (touch only at a line/point): {n_dropped}")
+        print(f"[{city}] cells kept: {len(grid):,}")
+    return grid
+
+
+def check_grid(grid, boundary_geom):
+    """Run the Task 2 grid checks on a grid built by make_grid (or read back from file).
+
+    Returns (results dict, list of failure messages). An empty list means every check passed.
+    """
+    cs = int(grid["cell_size_m"].iloc[0])
+    x_start, y_start, _, _ = grid_extent(boundary_geom.bounds, cs)
+    x_min = grid["x_min"].to_numpy(np.int64)
+    y_min = grid["y_min"].to_numpy(np.int64)
+    edge_share = grid.loc[grid["is_edge"], "area_share"]
+    cells_area = float(grid["area_m2"].sum())
+    boundary_area = float(boundary_geom.area)
+    b = grid.geometry.bounds
+    tol = 1e-6
+
+    r = {
+        "city": grid["city"].iloc[0],
+        "cell_size_m": cs,
+        "x_start": x_start,
+        "y_start": y_start,
+        "origin_multiple": x_start % cs == 0 and y_start % cs == 0,
+        "origin_matches_cells": int(x_min.min()) == x_start and int(y_min.min()) == y_start,
+        "n_cells": len(grid),
+        "n_full": int((~grid["is_edge"]).sum()),
+        "n_edge": int(grid["is_edge"].sum()),
+        **{f"n_edge_lt_{round(t * 100)}pct": int((edge_share < t).sum()) for t in SLIVER_THRESHOLDS},
+        "n_sliver": int(grid["is_sliver"].sum()),
+        "sliver_flag_ok": bool((grid["is_sliver"] == (grid["area_share"] < SLIVER_SHARE_THRESHOLD)).all()
+                               and not (grid["is_sliver"] & ~grid["is_edge"]).any()),
+        "min_area_share": float(grid["area_share"].min()),
+        "min_area_m2": float(grid["area_m2"].min()),
+        "cells_area_m2": cells_area,
+        "boundary_area_m2": boundary_area,
+        "area_diff_pct": abs(cells_area - boundary_area) / boundary_area * 100,
+        "cell_id_unique": bool(grid["cell_id"].is_unique),
+        "cell_id_matches_corner": bool((grid["cell_id"] == [f"{x}_{y}" for x, y in zip(x_min.tolist(), y_min.tolist())]).all()),
+        "xy_multiple": bool((x_min % cs == 0).all() and (y_min % cs == 0).all()),
+        "geoms_inside_own_square": bool(((b.minx >= x_min - tol) & (b.maxx <= x_min + cs + tol)
+                                         & (b.miny >= y_min - tol) & (b.maxy <= y_min + cs + tol)).all()),
+        "geoms_valid_nonempty": bool(grid.geometry.is_valid.all() and not grid.geometry.is_empty.any()),
+        "crs": grid.crs.to_string() if grid.crs else None,
+    }
+
+    failures = []
+    if not r["origin_multiple"]:
+        failures.append(f"origin ({x_start}, {y_start}) is not a multiple of {cs}")
+    if not r["origin_matches_cells"]:
+        failures.append("lowest cell corner does not equal the snapped origin")
+    if r["area_diff_pct"] >= MAX_AREA_DIFF_PCT:
+        failures.append(f"cell areas differ from boundary area by {r['area_diff_pct']:.6f}% (limit {MAX_AREA_DIFF_PCT}%)")
+    for key in ("cell_id_unique", "cell_id_matches_corner", "xy_multiple", "geoms_inside_own_square", "geoms_valid_nonempty",
+                "sliver_flag_ok"):
+        if not r[key]:
+            failures.append(f"check failed: {key}")
+    if (grid["area_share"] > 1 + 1e-9).any():
+        failures.append("a cell has area_share above 1")
+    return r, failures
+
+
+def plot_grid_qa(grid, boundary, path, title, subtitle):
+    """Save a QA map: full cells, edge cells in a second colour, boundary outline."""
+    import matplotlib
+
+    matplotlib.use("Agg")
+    import matplotlib.pyplot as plt
+    from matplotlib.lines import Line2D
+    from matplotlib.patches import Patch
+    from matplotlib.ticker import FuncFormatter
+
+    surface, ink, muted = "#fcfcfb", "#0b0b0b", "#52514e"
+    full_colour, edge_colour = "#2a78d6", "#eb6834"   # categorical slots 1 and 2
+
+    full = grid[~grid["is_edge"]]
+    edge = grid[grid["is_edge"]]
+    fig, ax = plt.subplots(figsize=(9, 9), dpi=150, facecolor=surface)
+    ax.set_facecolor(surface)
+    full.plot(ax=ax, facecolor=full_colour, edgecolor=surface, linewidth=0.5)
+    edge.plot(ax=ax, facecolor=edge_colour, edgecolor=surface, linewidth=0.5)
+    boundary.boundary.plot(ax=ax, color=ink, linewidth=1.0)
+
+    ax.set_aspect("equal")
+    fmt = FuncFormatter(lambda v, _: f"{v:,.0f}")
+    ax.xaxis.set_major_formatter(fmt)
+    ax.yaxis.set_major_formatter(fmt)
+    ax.tick_params(colors=muted, labelsize=7)
+    for spine in ax.spines.values():
+        spine.set_color("#c9c8c3")
+    ax.set_xlabel("Easting (m)", color=muted, fontsize=8)
+    ax.set_ylabel("Northing (m)", color=muted, fontsize=8)
+
+    handles = [
+        Patch(facecolor=full_colour, label=f"Full cell ({len(full):,})"),
+        Patch(facecolor=edge_colour, label=f"Edge cell, clipped ({len(edge):,})"),
+        Line2D([0], [0], color=ink, linewidth=1.0, label="City boundary"),
+    ]
+    ax.legend(handles=handles, loc="lower left", bbox_to_anchor=(0, 1.0), ncol=3,
+              frameon=False, fontsize=8, labelcolor=ink, borderaxespad=0.4)
+    fig.text(0.02, 0.975, title, fontsize=13, fontweight="bold", color=ink, va="top")
+    fig.text(0.02, 0.945, subtitle, fontsize=8.5, color=muted, va="top")
+    fig.tight_layout(rect=(0, 0, 1, 0.93))
+    fig.savefig(path, facecolor=surface)
+    plt.close(fig)
+    return path

requirements.txt

@@ -0,0 +1,18 @@
+certifi==2026.7.22
+contourpy==1.3.3
+cycler==0.12.1
+fonttools==4.66.1
+geopandas==1.2.0
+kiwisolver==1.5.1
+matplotlib==3.11.2
+numpy==2.4.6
+packaging==26.3
+pandas==3.0.6
+pillow==12.3.0
+pyogrio==0.13.0
+pyparsing==3.3.3
+pyproj==3.7.2
+python-dateutil==2.9.0.post0
+shapely==2.2.0
+six==1.17.0
+tzdata==2026.5

task1_crime.py

@@ -0,0 +1,114 @@
+import os
+import json
+import pandas as pd
+import geopandas as gpd
+
+with open("config.json", "r") as f:
+    config = json.load(f)
+
+data_dir = config["data_dir"]
+crime_dir = os.path.join(data_dir, "Crime_data_uk_cities")
+forces = list(config["forces"].keys())
+crs = config["crs"]
+
+months = [f"2025-{str(m).zfill(2)}" for m in range(1, 13)]
+
+# a. Verify files
+print("--- TASK 1a: Verifying Files ---")
+all_street_files = {force: [] for force in forces}
+other_files = []
+
+for month in months:
+    month_dir = os.path.join(crime_dir, month)
+    if not os.path.exists(month_dir):
+        print(f"Directory missing: {month_dir}")
+        continue
+        
+    for file in os.listdir(month_dir):
+        if file.endswith("-street.csv"):
+            matched_force = None
+            for force in forces:
+                if force in file:
+                    all_street_files[force].append(os.path.join(month_dir, file))
+                    matched_force = force
+                    break
+            if not matched_force:
+                other_files.append(os.path.join(month, file))
+        elif file.endswith(".csv"):
+            other_files.append(os.path.join(month, file))
+
+for force in forces:
+    found_months = []
+    for fpath in all_street_files[force]:
+        basename = os.path.basename(fpath)
+        found_months.append(basename[:7]) 
+    
+    missing = set(months) - set(found_months)
+    if missing:
+        print(f"Force '{force}' is missing months: {sorted(list(missing))}")
+    else:
+        print(f"Force '{force}' has all 12 months.")
+
+print(f"\nFound {len(other_files)} non-street files or files for other forces.")
+if other_files:
+    print(f"Example first 5: {other_files[:5]}")
+
+print("\n--- TASK 1b & 1c: Combine and Count ---")
+crime_dfs = {}
+
+for force in forces:
+    print(f"\nProcessing {force}...")
+    if not all_street_files[force]:
+        print(f"No files found for {force}. Skipping.")
+        continue
+        
+    dfs = []
+    for file in all_street_files[force]:
+        df = pd.read_csv(file)
+        dfs.append(df)
+            
+    combined_df = pd.concat(dfs, ignore_index=True)
+    
+    total_rows = len(combined_df)
+    missing_coords = combined_df['Longitude'].isna() | combined_df['Latitude'].isna()
+    missing_coords_count = missing_coords.sum()
+    
+    print(f"Total rows: {total_rows}")
+    print(f"Rows with missing coordinates: {missing_coords_count} ({missing_coords_count/total_rows*100:.2f}%)")
+    
+    print("\nRows per month:")
+    print(combined_df['Month'].value_counts().sort_index())
+    
+    print("\nRows per crime category:")
+    cats = combined_df['Crime type'].value_counts()
+    print(cats.to_string())
+    
+    asb_count = cats.get("Anti-social behaviour", 0)
+    all_crime_total = total_rows - asb_count
+    print(f"\nTotal excluding ASB: {all_crime_total}")
+    print(f"ASB Total: {asb_count}")
+    
+    crime_dfs[force] = combined_df
+
+print("\n--- TASK 1d: Spatial Conversion & Reprojection ---")
+for force, df in crime_dfs.items():
+    print(f"\nConverting {force} to spatial layer...")
+    df_valid_coords = df.dropna(subset=['Longitude', 'Latitude']).copy()
+    
+    gdf = gpd.GeoDataFrame(
+        df_valid_coords, 
+        geometry=gpd.points_from_xy(df_valid_coords.Longitude, df_valid_coords.Latitude),
+        crs="EPSG:4326"
+    )
+    
+    print("Reprojecting to EPSG:27700...")
+    gdf_27700 = gdf.to_crs("EPSG:27700")
+    
+    print("Sample coordinates (first 3):")
+    for geom in gdf_27700.geometry.head(3):
+        print(f"X: {geom.x:.1f}, Y: {geom.y:.1f}")
+        
+    output_path = os.path.join(data_dir, f"crime_{force}_EPSG27700.gpkg")
+    print(f"Saving to {output_path}...")
+    gdf_27700.to_file(output_path, driver="GPKG")
+    print("Saved.")

task2_boundaries_grids.py

@@ -0,0 +1,225 @@
+"""Task 2: extract the four city boundaries (2a) and build the analysis grids (2b).
+
+Run from the project folder:
+    .venv\\Scripts\\python.exe task2_boundaries_grids.py                  # grid_size_m from config.json
+    .venv\\Scripts\\python.exe task2_boundaries_grids.py --cell-size 250  # e.g. the later 250 m rerun
+
+All settings come from config.json; paths in it are relative to the project folder.
+Any failed check stops the run with an error (no workarounds).
+"""
+import argparse
+import json
+import platform
+import sys
+from datetime import datetime
+from pathlib import Path
+
+import geopandas as gpd
+import pandas as pd
+import pyogrio
+import pyproj
+import shapely
+
+from src.boundaries import (BOUNDARY_COLUMNS, boundary_filename, city_boundary, describe_lad_file,
+                            read_lad_features, require_bfc, to_target_crs)
+from src.geomutils import fingerprint, write_gpkg
+from src.grid import (EDGE_SHARE_THRESHOLD, GRID_COLUMNS, SLIVER_SHARE_THRESHOLD, SLIVER_THRESHOLDS, check_grid,
+                      grid_filename, make_grid, plot_grid_qa)
+
+ROOT = Path(__file__).resolve().parent
+CONFIG = ROOT / "config.json"
+
+# Sanity-check areas from the Task 2 brief (km2); a city must be within 5 percent.
+APPROX_AREA_KM2 = {"leeds": 552, "sheffield": 368, "newcastle": 113, "nottingham": 75}
+AREA_TOLERANCE = 0.05
+
+
+class Tee:
+    """Send printed output to the console and to a run log file."""
+
+    def __init__(self, *streams):
+        self.streams = streams
+
+    def write(self, text):
+        for s in self.streams:
+            if not s.closed:
+                s.write(text)
+
+    def flush(self):
+        for s in self.streams:
+            if not s.closed:
+                s.flush()
+
+
+def fail(message):
+    print(f"\nSTOP: {message}")
+    sys.exit(1)
+
+
+def section(title):
+    print(f"\n{'=' * 78}\n{title}\n{'=' * 78}")
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
+    parser.add_argument("--cell-size", type=int, help="override grid_size_m from config.json")
+    args = parser.parse_args()
+
+    cfg = json.loads(CONFIG.read_text())
+    crs = cfg["crs"]
+    cell_size = args.cell_size or cfg["grid_size_m"]
+    code_field, name_field = cfg["boundary_code_field"], cfg["boundary_name_field"]
+    lad_path = ROOT / cfg["boundary_file"]
+    boundaries_dir = ROOT / cfg["boundaries_dir"]
+    grids_dir = ROOT / cfg["grids_dir"]
+    qa_dir = ROOT / cfg["qa_dir"]
+    qa_dir.mkdir(parents=True, exist_ok=True)
+
+    run_log = open(qa_dir / f"task2_run_grid{cell_size}.txt", "w", encoding="utf-8")
+    sys.stdout = Tee(sys.__stdout__, run_log)
+    print(f"Task 2 run {datetime.now():%Y-%m-%d %H:%M:%S}")
+    print(f"python {platform.python_version()} | geopandas {gpd.__version__} | shapely {shapely.__version__} "
+          f"(GEOS {shapely.geos_version_string}) | pyproj {pyproj.__version__} (PROJ {pyproj.proj_version_str}) "
+          f"| pyogrio {pyogrio.__version__} (GDAL {pyogrio.__gdal_version_string__})")
+    print(f"config: crs={crs} cell_size={cell_size} m (config grid_size_m={cfg['grid_size_m']})")
+
+    # ---------------------------------------------------------------- Step 3a / 3b
+    section("STEP 3a/3b: ONS LAD boundary file")
+    if not lad_path.exists():
+        fail(f"boundary file not found: {lad_path}")
+    desc = describe_lad_file(lad_path)
+    if desc["crs"] is None:
+        fail("boundary file has no CRS defined")
+    if desc["code_fields"] != [code_field] or name_field not in desc["name_fields"]:
+        fail(f"config fields {code_field}/{name_field} do not match file columns {desc['columns']}")
+    try:
+        require_bfc(desc)
+    except ValueError as e:
+        fail(str(e))
+    print(f"-> ONS LAD {' '.join(desc['month'])} {' '.join(desc['year'])}, {code_field} code set; "
+          f"resolution BFC (full resolution, clipped) confirmed")
+
+    # ---------------------------------------------------------------- Step 3c / 3d
+    section("STEP 3c/3d: select cities by LAD code, check CRS")
+    cities = cfg["cities"]
+    all_names = pyogrio.read_dataframe(lad_path, layer=desc["layer"], columns=[code_field, name_field],
+                                       read_geometry=False)
+    lad = read_lad_features(lad_path, desc["layer"], code_field, [c["lad_code"] for c in cities.values()])
+    print(f"features read for the {len(cities)} configured codes: {len(lad)}")
+    for city, spec in cities.items():
+        expected = cfg["forces"][spec["force"]]
+        first_word = expected.split()[0]
+        lookalikes = all_names[all_names[name_field].str.startswith(first_word)
+                               & (all_names[name_field] != expected)]
+        if len(lookalikes):
+            others = ", ".join(f"{n} ({c})" for c, n in lookalikes[[code_field, name_field]].itertuples(index=False))
+            print(f"  {city}: similar names in file NOT selected: {others}")
+
+    lad, reprojected = to_target_crs(lad, crs)
+    print(f"file CRS: {desc['crs']} -> " + (f"reprojected to {crs} with to_crs" if reprojected
+                                            else f"already {crs}, no reprojection needed"))
+
+    # ---------------------------------------------------------------- Step 3c table / 3e / 3f
+    section("STEP 3c/3e/3f: city boundaries (dissolve, validity, save)")
+    boundaries, reports = {}, []
+    for city, spec in cities.items():
+        expected = cfg["forces"][spec["force"]]
+        try:
+            b, rep = city_boundary(lad, city, spec["lad_code"], expected, code_field, name_field)
+        except ValueError as e:
+            fail(str(e))
+        area_km2 = b["area_m2"].iloc[0] / 1e6
+        approx = APPROX_AREA_KM2[city]
+        rep.update(lad_code=spec["lad_code"], lad_name=b["lad_name"].iloc[0], name_ok=b["lad_name"].iloc[0] == expected,
+                   area_km2=area_km2, approx_km2=approx, area_dev_pct=(area_km2 - approx) / approx * 100)
+        if abs(area_km2 - approx) / approx > AREA_TOLERANCE:
+            fail(f"{city}: area {area_km2:.1f} km2 is more than {AREA_TOLERANCE:.0%} from the expected ~{approx} km2")
+        if rep["repaired"]:
+            print(f"  {city}: geometry was INVALID ({rep['invalid_reason']}); repaired with shapely.make_valid")
+        path = write_gpkg(b[BOUNDARY_COLUMNS + ["geometry"]], boundaries_dir / boundary_filename(city, crs))
+        saved = gpd.read_file(path, engine="pyogrio")
+        saved_crs = pyogrio.read_info(path)["crs"]
+        if saved_crs != crs or len(saved) != 1:
+            fail(f"{path.name}: saved CRS {saved_crs}, rows {len(saved)}")
+        rep.update(file=path.name, saved_crs=saved_crs,
+                   fingerprint=fingerprint(saved, BOUNDARY_COLUMNS))
+        boundaries[city] = saved
+        reports.append(rep)
+
+    btable = pd.DataFrame(reports)
+    print(btable[["city", "lad_code", "lad_name", "name_ok", "area_km2", "approx_km2", "area_dev_pct"]]
+          .to_string(index=False, float_format=lambda v: f"{v:,.2f}"))
+    print()
+    print(btable[["city", "n_features", "dissolved", "was_valid", "repaired", "n_polygons", "file", "saved_crs"]]
+          .to_string(index=False))
+    btable.to_csv(qa_dir / "task2_boundaries_summary.csv", index=False)
+
+    # ---------------------------------------------------------------- Step 5 / 6
+    section(f"STEP 5/6: {cell_size} m grids and checks")
+    results = []
+    for city, boundary in boundaries.items():
+        geom = boundary.geometry.iloc[0]   # the boundary exactly as saved and shared
+        grid = make_grid(geom, cell_size, city, crs=crs)
+        rerun = make_grid(geom, cell_size, city, crs=crs, verbose=False)
+        fp = fingerprint(grid, GRID_COLUMNS)
+        deterministic = fp == fingerprint(rerun, GRID_COLUMNS)
+
+        path = write_gpkg(grid, grids_dir / grid_filename(city, cell_size, crs))
+        saved = gpd.read_file(path, engine="pyogrio")
+        saved_crs = pyogrio.read_info(path)["crs"]
+        roundtrip = fingerprint(saved, GRID_COLUMNS) == fp
+
+        r, failures = check_grid(saved, geom)
+        r.update(n_dropped_zero_area=grid.attrs["n_dropped_zero_area"], saved_crs=saved_crs,
+                 deterministic=deterministic, file_roundtrip_identical=roundtrip, file=path.name, fingerprint=fp)
+        if saved_crs != crs:
+            failures.append(f"saved file CRS is {saved_crs}, not {crs}")
+        if not deterministic:
+            failures.append("two runs of make_grid gave different output")
+        if not roundtrip:
+            failures.append("grid read back from file differs from the grid written")
+
+        png = qa_dir / f"{city}_grid{cell_size}_qa.png"
+        name = boundary["lad_name"].iloc[0]
+        plot_grid_qa(saved, boundary, png,
+                     title=f"{name}: {cell_size} m analysis grid ({crs})",
+                     subtitle=(f"{r['n_cells']:,} cells, origin ({r['x_start']:,}, {r['y_start']:,}); "
+                               f"edge = area_share < {EDGE_SHARE_THRESHOLD}; "
+                               f"slivers <5%: {r['n_edge_lt_5pct']}, <10%: {r['n_edge_lt_10pct']}, "
+                               f"<25%: {r['n_edge_lt_25pct']}"))
+        r["qa_png"] = png.name
+        results.append(r)
+
+        print(f"[{city}] a. origin ({r['x_start']}, {r['y_start']}), multiples of {cell_size}: {r['origin_multiple']}, "
+              f"equals lowest cell corner: {r['origin_matches_cells']}")
+        print(f"[{city}] b. cells {r['n_cells']:,} = full {r['n_full']:,} + edge {r['n_edge']:,}")
+        print(f"[{city}] c. edge cells with area_share <5%: {r['n_edge_lt_5pct']}, <10%: {r['n_edge_lt_10pct']}, "
+              f"<25%: {r['n_edge_lt_25pct']} (smallest: {r['min_area_share']:.6f} = {r['min_area_m2']:,.2f} m2)")
+        print(f"[{city}]    is_sliver (area_share < {SLIVER_SHARE_THRESHOLD}): {r['n_sliver']}, "
+              f"flag consistent and only on edge cells: {r['sliver_flag_ok']}")
+        print(f"[{city}] d. cell area sum {r['cells_area_m2']:,.3f} m2 vs boundary {r['boundary_area_m2']:,.3f} m2, "
+              f"diff {r['area_diff_pct']:.2e} %")
+        print(f"[{city}] e. cell_id unique: {r['cell_id_unique']}, matches corner: {r['cell_id_matches_corner']}")
+        print(f"[{city}] f. all x_min/y_min multiples of {cell_size}: {r['xy_multiple']}")
+        print(f"[{city}] g. saved CRS: {saved_crs}")
+        print(f"[{city}] h. QA map: {png.relative_to(ROOT)}")
+        print(f"[{city}]    deterministic rerun identical: {deterministic}; file round-trip identical: {roundtrip}; "
+              f"geometries valid: {r['geoms_valid_nonempty']}; inside own square: {r['geoms_inside_own_square']}")
+        print(f"[{city}]    saved {path.relative_to(ROOT)} ({len(saved):,} rows), fingerprint {fp[:16]}...")
+        if failures:
+            fail(f"{city} grid checks failed: " + "; ".join(failures))
+        print()
+
+    gtable = pd.DataFrame(results)
+    gtable.to_csv(qa_dir / f"task2_grid{cell_size}_checks.csv", index=False)
+    section("SUMMARY")
+    cols = ["city", "x_start", "y_start", "n_cells", "n_full", "n_edge", "n_dropped_zero_area",
+            *[f"n_edge_lt_{round(t * 100)}pct" for t in SLIVER_THRESHOLDS], "n_sliver", "area_diff_pct"]
+    print(gtable[cols].to_string(index=False, float_format=lambda v: f"{v:.2e}"))
+    print("\nALL TASK 2 CHECKS PASSED")
+    sys.stdout = sys.__stdout__
+    run_log.close()
+
+
+if __name__ == "__main__":
+    main()

task2_log.txt

@@ -0,0 +1,144 @@
+TASK 2 LOG: city boundaries and 500 m analysis grids
+Date: 2026-10-09
+Machine: the copied project on a new laptop (the original Task 0/1 work was done on another computer)
+
+WHAT WAS DONE
+1. Environment
+   - Python 3.11.14 (Anaconda). There was no requirements.txt, so I created .venv in the
+     project folder and installed pandas, geopandas, shapely, pyproj, pyogrio, numpy and
+     matplotlib.
+   - Versions: geopandas 1.2.0, shapely 2.2.0 (GEOS 3.14.1), pyproj 3.7.2 (PROJ 9.5.1),
+     pyogrio 0.13.0 (GDAL 3.12.4).
+   - I then wrote requirements.txt (pip freeze) so teammates can install the same versions.
+
+2. Checks on the copied project
+   - config.json line 13: "data_dir" changed from
+       "c:\\Users\\mukun\\Documents\\SEM_3\\6.SI\\MID_project"
+     to ".". Nothing else was changed at that step.
+   - Task 1 outputs load correctly. All four layers are in EPSG:27700, with six-figure
+     coordinates in metres, and each covers 2025-01 to 2025-12.
+   - Row counts in the .gpkg files compared with the expected totals:
+       west-yorkshire   298,385 vs 300,547 (2,162 rows with no coordinates)
+       south-yorkshire  164,140 vs 172,209 (8,069 rows with no coordinates)
+       northumbria      169,796 vs 169,796 (no rows without coordinates)
+       nottinghamshire  120,552 vs 121,358 (806 rows with no coordinates)
+     The expected totals equal the raw 2025 street CSV row counts exactly. Each gap equals
+     the number of rows with no coordinates, which task1_crime.py drops (line 96) before
+     saving. Conclusion: the copy is complete, and the expected totals counted rows before
+     that drop. The user confirmed to continue.
+
+3. Task 2a: boundaries
+   - No LAD boundary file was in the copied folder (checked every folder and all 9 zip
+     archives). config.json also had no LAD codes or boundary path.
+   - With the user's permission I downloaded the ONS "Local Authority Districts (December
+     2025) Boundaries UK BFC" GeoPackage from the ONS Open Geography Portal
+     (item 92150c7aa60540c5814abe3b26bce6d0) into data/raw/boundaries/LAD_DEC_2025_UK_BFC.gpkg.
+       SHA-256: 29ba41a24d7943db1e8ab50c73506cb486c555f8c0e9bc945b24aece068d3602
+       Contents: 1 layer, 361 MultiPolygon features, CRS EPSG:27700.
+       Columns: LAD25CD, LAD25NM, LAD25NMW, BNG_E, BNG_N, LONG, LAT, GlobalID.
+       Version: December 2025, BFC (full resolution, clipped), LAD25 code set.
+   - I chose December 2025 over May 2025 (V2) because it is the latest release and shows
+     boundaries at the end of the study year. Both releases give the same areas for the
+     four cities (difference under 1 m2).
+   - Added to config.json: "cities" (city key -> force and LAD code), "boundary_file",
+     "boundary_code_field", "boundary_name_field", "boundaries_dir", "grids_dir", "qa_dir".
+     Existing keys were not changed.
+   - Cities matched by code. Areas are compared with the brief's approximate values:
+       leeds       E08000035  Leeds                551.70 km2 (brief ~552, -0.05%)
+       sheffield   E08000039  Sheffield            368.00 km2 (brief ~368,  0.00%)
+       newcastle   E08000021  Newcastle upon Tyne  113.45 km2 (brief ~113, +0.40%)
+       nottingham  E06000018  Nottingham            74.64 km2 (brief ~75,  -0.47%)
+     Newcastle-under-Lyme (E07000195) is in the file and was not selected.
+   - Reprojection was not needed (the file is already EPSG:27700).
+   - No dissolve was needed (one feature per city). All four geometries were valid, so no
+     make_valid repair was needed.
+   - Saved data/interim/boundaries/<city>_boundary_EPSG27700.gpkg
+     (columns city, lad_code, lad_name, area_m2, geometry).
+
+4. Task 2b: grids
+   - Code: src/grid.py (make_grid, grid_extent, check_grid, plot_grid_qa),
+     src/boundaries.py, src/geomutils.py, and the runner task2_boundaries_grids.py.
+   - Built 500 m grids (grid_size_m from config.json) and saved
+     data/interim/grids/<city>_grid500_EPSG27700.gpkg.
+
+5. Checks (all passed)
+     city        origin            cells  full  edge  dropped  <5%  <10%  <25%  area diff
+     leeds       (413000, 422500)  2,390  2,034  356       0    52    79   114   4.5e-13 %
+     sheffield   (413000, 378500)  1,626  1,326  300       0    39    57   102   1.1e-13 %
+     newcastle   (414000, 562500)    526    384  142       0    22    31    47   9.2e-14 %
+     nottingham  (450500, 332500)    366    234  132       0    22    25    41   0
+   - Origins and all x_min/y_min values are multiples of 500.
+   - cell_id is unique and equals "<x_min>_<y_min>".
+   - Every clipped cell lies inside its own 500 m square, and all geometries are valid.
+   - The saved files are EPSG:27700.
+   - Running make_grid twice gave identical output, and the files read back identical to
+     what was written. A full second run of the script also gave the same fingerprints:
+       leeds       43f0b4f5addd131a1c6ca37b7cb3227b450748a59371d183e4f754023ba40366
+       sheffield   e2eaab7f9ee18c4485851b0e4dd3cf8ce40741d91efc8af6d925672ed122d061
+       newcastle   0c7f42dbd48911f71014fd3a1e00723fd41f39120439eb344cd6ab1352ac573f
+       nottingham  2c231069a277ad73d8cce6f3a077255a816afa00fae314804a4288f2e127d040
+   - QA maps: outputs/qa/<city>_grid500_qa.png. Full run output:
+     outputs/qa/task2_run_grid500.txt. Tables: outputs/qa/task2_boundaries_summary.csv
+     and outputs/qa/task2_grid500_checks.csv.
+
+6. Records: added a Task 2 section at the end of conventions.md (no other section changed).
+
+THINGS THAT LOOKED SUSPICIOUS OR NEED A DECISION LATER
+- Rows without coordinates: Task 1 silently drops crime rows with no coordinates
+  (South Yorkshire 4.7% of rows, West Yorkshire 0.7%, Nottinghamshire 0.7%). This rule is
+  not written in conventions.md, so the Task 1 owner should add it and put the counts in
+  the data quality table.
+- Crimes placed outside the force area: the first three South Yorkshire records sit at
+  (440514, 350921), which is outside South Yorkshire. Police forces record some crimes
+  outside their own area. Clipping to the city boundary in Task 3 will remove them, but
+  the number removed should be reported.
+- Brief vs reality: the brief said config.json held LAD codes and paths, but it only had
+  data_dir. The boundary file was also missing from the copy. Either the teammate never
+  had one, or it was not copied.
+- Sheffield's LAD code is E08000039 (changed in 2023). Older lists and older boundary files
+  use E08000019, so anyone joining other data by LAD code must use the 2023+ code.
+- Tiny slivers are kept as instructed:
+    Leeds 419500_439000 = 0.16 m2
+    Sheffield 433000_379500 = 0.11 m2
+    23 edge cells under 1% of a cell in both Leeds and Sheffield
+  Densities (crime per km2) for these cells will be extreme. Tasks 3 and later need a rule,
+  for example a minimum area_share for density maps.
+- One Leeds cell (442500_439500) is clipped very slightly (area_share 0.99998) but counts
+  as full under the 0.9999 rule. This is expected behaviour, noted for transparency.
+- conventions.md "Data Location" (General Settings) still shows the old absolute path
+  c:\Users\mukun\... . It was left unchanged because Task 2 must not edit other sections.
+  config.json now uses ".".
+- The project folder also holds 2024 crime data and 2024 OSM extracts. Task 1 used 2025
+  only, which I verified from the month ranges.
+- My script's first full run finished every step and passed every check, then raised a
+  harmless error while closing its own log file (exit code 120). I fixed it and reran:
+  all outputs were identical (same fingerprints) and the run exited cleanly.
+- .venv is inside the OneDrive folder, so OneDrive syncs thousands of package files.
+  Consider excluding .venv from sync. Teammates should build their own from
+  requirements.txt rather than copy it.
+- The raw LAD file (122 MB) is not in the handoff folder because of its size. Teammates can
+  re-download it (see the item ID above) and check the SHA-256, or copy it from
+  data/raw/boundaries/.
+
+ADDENDUM 2026-10-09: is_sliver column and conventions.md edits
+- conventions.md edits:
+  - Data Location now uses the relative project paths from config.json.
+  - Added the "Rows Without Coordinates" rule to the crime data section.
+  - In the Task 2 section:
+    - The Sheffield code sentence now reads "E08000039 ... E08000019 must not be used".
+    - Added the sliver rule, marked PROPOSED TEAM DECISION (not yet agreed).
+    - Added is_sliver to the grid column list.
+- Grid files: added is_sliver = area_share < 0.10. It is the last column before geometry,
+  so all existing columns keep their positions.
+  - Made in make_grid (src/grid.py, SLIVER_SHARE_THRESHOLD = 0.10), so reruns (including
+    the 250 m grid) include it automatically. check_grid also verifies the flag.
+  - Slivers: leeds 79, sheffield 57, newcastle 31, nottingham 25.
+- Checked against the original files: cell counts, total areas, every original column,
+  every cell_id and every geometry (WKB bytes) are identical. The boundary files are
+  unchanged in content.
+- All Task 2 checks passed again. Two separate runs gave identical fingerprints.
+- The fingerprints now include is_sliver, so they replace the ones listed above:
+    leeds       5e6d3d67fd9e7836d1f1cdd2ffaec91b33dd968c217cad3b2e4c4a425a4718bf
+    sheffield   2b56dc075f491a1adf1230d40d5fadb948ca41138fe2d2d849bf33ca5a61a6fb
+    newcastle   95ad4c0669bd220069d8b29a88b95d2e8e145affc162a3a8d183ce2ede1cc68f
+    nottingham  85f36bc798012451f3eaff6dc2895f707600c1575a34d2dc429fb016629b0cd9

task2_qgis.py

@@ -0,0 +1,782 @@
+"""Task 2 reproduced in QGIS: city boundaries (2a) and analysis grids (2b) with native Processing tools.
+
+HOW TO RUN (QGIS 3.44 LTR)
+  1. Start QGIS and create a NEW, empty project (Project > New).
+  2. Plugins > Python Console, then click "Show Editor" in the console toolbar.
+  3. In the editor click "Open Script..." and choose this file (task2_qgis.py).
+  4. Check that PROJECT_ROOT below is the project folder, then click "Run Script".
+  5. Wait for "TASK 2 IN QGIS: FINISHED" in the console (a few minutes; QGIS is busy meanwhile).
+
+Every geoprocessing step is a processing.run() call to a native QGIS algorithm (the same tools as the
+Processing Toolbox). Every call, with its full parameters and output feature counts, is written to
+logs/task2_qgis_log.txt so it can be repeated by hand.
+
+Writes only to:
+  data/interim/qgis/boundaries/<city>_boundary_EPSG27700_qgis.gpkg
+  data/interim/qgis/grids/<city>_grid500_EPSG27700_qgis.gpkg
+  data/interim/qgis/steps/<city>/NN_<step>.gpkg       every intermediate result
+  outputs/qa/qgis/                                    print layout PNGs and task2_qgis_vs_python.csv
+  logs/task2_qgis_log.txt
+  QGIS/crimegeo_vgi.qgz                               team project (relative paths, print layouts)
+The Python Task 2 outputs are only read, and their checksums are verified before and after the run.
+"""
+# ---------------------------------------------------------------------------------------------------
+# EDIT THIS: absolute path of the project folder (the folder that holds config.json).
+PROJECT_ROOT = r"C:\Users\padma\OneDrive\Desktop\SEM 7\SI (Spatial Informatics)\Project\MidPresentation\MID_project"
+# ---------------------------------------------------------------------------------------------------
+
+import csv
+import hashlib
+import json
+import math
+import os
+import re
+import time
+import traceback
+import zipfile
+from datetime import datetime
+
+import processing
+from qgis.core import (Qgis, QgsCoordinateReferenceSystem, QgsExpression, QgsFeatureRequest, QgsFillSymbol,
+                       QgsLayoutExporter, QgsLayoutItemLabel, QgsLayoutItemLegend, QgsLayoutItemMap,
+                       QgsLayoutItemPage, QgsLayoutItemPicture, QgsLayoutItemScaleBar, QgsLayoutPoint,
+                       QgsLayoutSize, QgsPrintLayout, QgsProcessingContext, QgsProject, QgsRasterLayer,
+                       QgsRectangle, QgsRuleBasedRenderer, QgsTextFormat, QgsVectorLayer)
+from qgis.PyQt.QtCore import QCoreApplication, QVariant
+from qgis.PyQt.QtGui import QFont
+
+# Grid rules: the same values as src/grid.py and conventions.md.
+EDGE_SHARE_THRESHOLD = 0.9999
+SLIVER_SHARE_THRESHOLD = 0.10
+SLIVER_REPORT_THRESHOLDS = (0.05, 0.10, 0.25)
+MAX_AREA_DIFF_PCT = 0.01          # sum of cell areas vs boundary area
+MAX_CELL_AREA_DIFF_M2 = 0.01      # QGIS vs Python, per cell
+# Python results at 500 m, from task2_log.txt: (x_start, y_start, cells).
+EXPECTED_PYTHON_500 = {"leeds": (413000, 422500, 2390), "sheffield": (413000, 378500, 1626),
+                       "newcastle": (414000, 562500, 526), "nottingham": (450500, 332500, 366)}
+
+BASEMAP_NAME = "OpenStreetMap"
+BASEMAP_URI = "type=xyz&url=https://tile.openstreetmap.org/{z}/{x}/{y}.png&zmin=0&zmax=19"
+SOURCE_LINE = "Boundary: ONS LAD December 2025 BFC. Basemap: \u00a9 OpenStreetMap contributors."
+NORTH_ARROW_SVG = ":/images/north_arrows/layout_default_north_arrow.svg"
+
+TEXT, BOOL, INT64, DOUBLE = int(QVariant.String), int(QVariant.Bool), int(QVariant.LongLong), int(QVariant.Double)
+TYPE_NAMES = {TEXT: "text", BOOL: "boolean", INT64: "int8", DOUBLE: "double precision"}
+
+
+def project_path(*parts):
+    return os.path.normpath(os.path.join(PROJECT_ROOT, *parts))
+
+
+if not os.path.isfile(project_path("config.json")):
+    raise FileNotFoundError(f"config.json not found: set PROJECT_ROOT at the top of the script (now {PROJECT_ROOT})")
+with open(project_path("config.json"), encoding="utf-8") as _f:
+    CFG = json.load(_f)
+
+CRS = CFG["crs"]
+CS = int(CFG["grid_size_m"])
+TAG = CRS.replace(":", "")
+CODE_FIELD, NAME_FIELD = CFG["boundary_code_field"], CFG["boundary_name_field"]
+LAD_PATH = project_path(*CFG["boundary_file"].split("/"))
+OUT_DIR = project_path("data", "interim", "qgis")
+BND_DIR = os.path.join(OUT_DIR, "boundaries")
+GRID_DIR = os.path.join(OUT_DIR, "grids")
+STEP_DIR = os.path.join(OUT_DIR, "steps")
+QA_DIR = project_path("outputs", "qa", "qgis")
+LOG_PATH = project_path("logs", "task2_qgis_log.txt")
+CSV_PATH = os.path.join(QA_DIR, "task2_qgis_vs_python.csv")
+PROJECT_FILE = project_path("QGIS", "crimegeo_vgi.qgz")
+PY_BND_DIR = project_path(*CFG["boundaries_dir"].split("/"))
+PY_GRID_DIR = project_path(*CFG["grids_dir"].split("/"))
+PROTECTED = ([os.path.join(PY_BND_DIR, f"{c}_boundary_{TAG}.gpkg") for c in CFG["cities"]]
+             + [os.path.join(PY_GRID_DIR, f"{c}_grid{CS}_{TAG}.gpkg") for c in CFG["cities"]]
+             + [project_path("src", "grid.py"), project_path("src", "boundaries.py"),
+                project_path("task2_boundaries_grids.py")])
+
+
+class CheckFailed(Exception):
+    pass
+
+
+class Log:
+    """Print to the console and write the same text to logs/task2_qgis_log.txt."""
+
+    def __init__(self, path):
+        os.makedirs(os.path.dirname(path), exist_ok=True)
+        self.file = open(path, "w", encoding="utf-8")
+
+    def __call__(self, message=""):
+        print(message)
+        self.file.write(message + "\n")
+        self.file.flush()
+
+    def close(self):
+        self.file.close()
+
+
+LOG = None
+CALLS = []
+
+
+# ------------------------------------------------------------------------------------------ helpers
+def boundary_path(city):
+    return os.path.join(BND_DIR, f"{city}_boundary_{TAG}_qgis.gpkg")
+
+
+def grid_path(city):
+    return os.path.join(GRID_DIR, f"{city}_grid{CS}_{TAG}_qgis.gpkg")
+
+
+def step_path(city, name):
+    return os.path.join(STEP_DIR, city, f"{name}.gpkg")
+
+
+def city_outputs(city):
+    steps = ["01_extract_by_code", "02_valid", "02_invalid", "02_validity_errors", "03_fixed", "04_dissolved",
+             "05_grid_full_extent", "06_clipped", "07_area_gt_0", "07_dropped_zero_area"]
+    return [step_path(city, s) for s in steps] + [boundary_path(city), grid_path(city),
+                                                  os.path.join(QA_DIR, f"{city}_grid_qgis.png")]
+
+
+def sha256(path):
+    h = hashlib.sha256()
+    with open(path, "rb") as f:
+        for block in iter(lambda: f.read(1 << 20), b""):
+            h.update(block)
+    return h.hexdigest()
+
+
+def open_layer(path, name="layer"):
+    """Open a vector file read-only (never modifies it)."""
+    options = QgsVectorLayer.LayerOptions()
+    options.forceReadOnly = True
+    options.loadDefaultStyle = False
+    layer = QgsVectorLayer(path, name, "ogr", options)
+    if not layer.isValid():
+        raise RuntimeError(f"cannot open {path}")
+    return layer
+
+
+def feature_count(path):
+    layer = open_layer(path)
+    n = layer.featureCount()
+    del layer
+    return n
+
+
+def make_context():
+    """Processing context: validity is checked explicitly with native:checkvalidity, so no silent filtering."""
+    context = QgsProcessingContext()
+    context.setProject(QgsProject.instance())
+    no_check = (Qgis.InvalidGeometryCheck.NoCheck if hasattr(Qgis, "InvalidGeometryCheck")
+                else QgsFeatureRequest.GeometryNoCheck)
+    context.setInvalidGeometryCheck(no_check)
+    return context
+
+
+def run_alg(city, step, alg_id, params):
+    """processing.run() one native algorithm and log ID, full parameters and output feature counts."""
+    LOG(f"\n[{city}] {step}\n  algorithm: {alg_id}")
+    LOG("  parameters: " + json.dumps(params, indent=2, ensure_ascii=False).replace("\n", "\n  "))
+    result = processing.run(alg_id, params, context=make_context())
+    counts = {}
+    for key, value in result.items():
+        if isinstance(value, str) and value.split("|")[0].lower().endswith(".gpkg"):
+            counts[key] = feature_count(value.split("|")[0])
+        elif key.endswith("_COUNT"):
+            counts[key] = int(value)
+    LOG("  output feature counts: " + ", ".join(f"{k}={v}" for k, v in counts.items()))
+    CALLS.append({"city": city, "step": step, "algorithm": alg_id, "counts": counts})
+    return result, counts
+
+
+def fmap(name, qtype, expression):
+    """One entry of a native:refactorfields FIELDS_MAPPING."""
+    return {"name": name, "type": qtype, "type_name": TYPE_NAMES[qtype], "sub_type": 0, "length": 0,
+            "precision": 0, "expression": expression, "alias": "", "comment": ""}
+
+
+def remove_files(paths, wait_seconds=75):
+    """Delete old outputs before a rerun. Waits for QGIS to release files it still holds open."""
+    pending = [p for p in paths if os.path.exists(p)]
+    deadline = time.time() + wait_seconds
+    while pending:
+        still = []
+        for p in pending:
+            for extra in ("", "-wal", "-shm", "-journal"):
+                if os.path.exists(p + extra):
+                    try:
+                        os.remove(p + extra)
+                    except PermissionError:
+                        still.append(p)
+                        break
+        pending = still
+        if pending:
+            if time.time() > deadline:
+                raise RuntimeError("these old outputs are still open in QGIS and cannot be replaced: "
+                                   f"{pending}. Restart QGIS, create a new project, and run the script again.")
+            QCoreApplication.processEvents()
+            time.sleep(1)
+
+
+def clear_previous_run(project):
+    """Remove layers, groups and layouts that an earlier run of this script added to the open project."""
+    out = os.path.normcase(OUT_DIR)
+    labels = [CFG["forces"][s["force"]] for s in CFG["cities"].values()]
+    for layer in list(project.mapLayers().values()):
+        source = os.path.normcase(os.path.normpath(layer.source().split("|")[0]))
+        if source.startswith(out) or layer.name() == BASEMAP_NAME:
+            project.removeMapLayer(layer.id())
+    root = project.layerTreeRoot()
+    for label in labels:
+        group = root.findGroup(label)
+        if group is not None and not group.children():
+            root.removeChildNode(group)
+        layout = project.layoutManager().layoutByName(f"{label} grid QA")
+        if layout is not None:
+            project.layoutManager().removeLayout(layout)
+    others = [layer.name() for layer in project.mapLayers().values()]
+    if others:
+        LOG(f"WARNING: the open project already contains other layers {others}; they will be saved in "
+            f"{os.path.basename(PROJECT_FILE)} too. Start from Project > New for a clean team project.")
+
+
+def read_grid(path):
+    layer = open_layer(path)
+    cells, duplicates = {}, 0
+    for f in layer.getFeatures():
+        cid = f["cell_id"]
+        if cid in cells:
+            duplicates += 1
+        cells[cid] = {"x_min": int(f["x_min"]), "y_min": int(f["y_min"]), "area_m2": float(f["area_m2"]),
+                      "area_share": float(f["area_share"]), "is_edge": bool(f["is_edge"]),
+                      "is_sliver": bool(f["is_sliver"])}
+    schema = [(fld.name(), fld.typeName()) for fld in layer.fields() if fld.name() != "fid"]
+    info = {"cells": cells, "n": layer.featureCount(), "duplicates": duplicates, "schema": schema,
+            "crs": layer.crs().authid()}
+    del layer
+    return info
+
+
+def read_boundary(path):
+    layer = open_layer(path)
+    feats = list(layer.getFeatures())
+    info = {"n": len(feats), "crs": layer.crs().authid(), "area_attr": float(feats[0]["area_m2"]),
+            "area_geom": feats[0].geometry().area(), "bbox": feats[0].geometry().boundingBox()}
+    del layer
+    return info
+
+
+# ------------------------------------------------------------------------------------- Step 1 (2a)
+def inspect_lad_file():
+    """STEP 1a/1b: load the ONS LAD file, report it, check the CRS (reproject once if needed)."""
+    LOG(f"\n{'=' * 90}\nSTEP 1a/1b: ONS LAD boundary file\n{'=' * 90}")
+    if not os.path.isfile(LAD_PATH):
+        raise CheckFailed(f"boundary file not found: {LAD_PATH}")
+    lad = open_layer(LAD_PATH, "LAD")
+    fields = [f.name() for f in lad.fields()]
+    LOG(f"file:     {LAD_PATH}")
+    LOG(f"CRS:      {lad.crs().authid() or '(none)'}  valid={lad.crs().isValid()}")
+    LOG(f"features: {lad.featureCount()}")
+    LOG(f"fields:   {fields}")
+    code_fields = [f for f in fields if re.fullmatch(r"LAD\d\dCD", f)]
+    name_fields = [f for f in fields if re.fullmatch(r"LAD\d\dNM", f)]
+    LOG(f"LAD code field: {code_fields}, name field: {name_fields} (config: {CODE_FIELD}, {NAME_FIELD})")
+    if code_fields != [CODE_FIELD] or NAME_FIELD not in name_fields:
+        raise CheckFailed("LAD code/name fields in the file do not match config.json")
+    if not lad.crs().isValid():
+        raise CheckFailed("the LAD file has no CRS defined")
+    source = LAD_PATH
+    if lad.crs().authid() != CRS:
+        LOG(f"CRS is {lad.crs().authid()}, not {CRS}: reprojecting with native:reprojectlayer")
+        res, _ = run_alg("all", "STEP 1b reproject LAD file", "native:reprojectlayer", {
+            "INPUT": LAD_PATH, "TARGET_CRS": CRS, "CONVERT_CURVED_GEOMETRIES": False,
+            "OUTPUT": os.path.join(STEP_DIR, "00_lad_reprojected.gpkg")})
+        source = res["OUTPUT"]
+    else:
+        LOG(f"CRS is already {CRS}: no reprojection needed (native:reprojectlayer skipped)")
+    del lad
+    return source
+
+
+def build_boundary(city, spec, lad_source):
+    """STEP 1c-1f for one city. Returns the path of the dissolved layer and the final boundary."""
+    label, code = CFG["forces"][spec["force"]], spec["lad_code"]
+
+    res, n = run_alg(city, "STEP 1c extract by LAD code", "native:extractbyattribute", {
+        "INPUT": lad_source, "FIELD": CODE_FIELD, "OPERATOR": 0, "VALUE": code,
+        "OUTPUT": step_path(city, "01_extract_by_code")})
+    if n["OUTPUT"] != 1:
+        raise CheckFailed(f"{code}: expected exactly 1 feature, got {n['OUTPUT']}")
+    extracted = open_layer(res["OUTPUT"])
+    name = next(extracted.getFeatures())[NAME_FIELD]
+    del extracted
+    if name != label:
+        raise CheckFailed(f"{code} is named '{name}' in the file, expected '{label}'")
+    LOG(f"  name in file: '{name}' = expected '{label}'")
+
+    res_v, nv = run_alg(city, "STEP 1d check validity (GEOS)", "native:checkvalidity", {
+        "INPUT_LAYER": res["OUTPUT"], "METHOD": 2, "IGNORE_RING_SELF_INTERSECTION": False,
+        "VALID_OUTPUT": step_path(city, "02_valid"), "INVALID_OUTPUT": step_path(city, "02_invalid"),
+        "ERROR_OUTPUT": step_path(city, "02_validity_errors")})
+    to_dissolve, repaired = res["OUTPUT"], nv["INVALID_COUNT"] > 0
+    if repaired:
+        LOG(f"  {nv['INVALID_COUNT']} INVALID geometry -> repairing with native:fixgeometries")
+        res_f, _ = run_alg(city, "STEP 1d fix geometries", "native:fixgeometries", {
+            "INPUT": res["OUTPUT"], "METHOD": 0, "OUTPUT": step_path(city, "03_fixed")})
+        to_dissolve = res_f["OUTPUT"]
+    else:
+        LOG("  geometry valid -> native:fixgeometries skipped")
+
+    res_d, nd = run_alg(city, "STEP 1e dissolve to one feature", "native:dissolve", {
+        "INPUT": to_dissolve, "FIELD": [], "SEPARATE_DISJOINT": False, "OUTPUT": step_path(city, "04_dissolved")})
+    if nd["OUTPUT"] != 1:
+        raise CheckFailed(f"dissolve gave {nd['OUTPUT']} features, expected 1")
+
+    run_alg(city, "STEP 1f boundary schema (planar area)", "native:refactorfields", {
+        "INPUT": res_d["OUTPUT"],
+        "FIELDS_MAPPING": [fmap("city", TEXT, QgsExpression.quotedString(city)),
+                           fmap("lad_code", TEXT, QgsExpression.quotedColumnRef(CODE_FIELD)),
+                           fmap("lad_name", TEXT, QgsExpression.quotedColumnRef(NAME_FIELD)),
+                           fmap("area_m2", DOUBLE, "area($geometry)")],
+        "OUTPUT": boundary_path(city)})
+    return res_d["OUTPUT"], repaired
+
+
+# ------------------------------------------------------------------------------------- Step 2 (2b)
+def build_grid(city, dissolved_path, python_origin):
+    """STEP 2a-2f for one city. Returns the snapped origin and the number of zero-area cells dropped."""
+    dissolved = open_layer(dissolved_path)
+    bbox = next(dissolved.getFeatures()).geometry().boundingBox()
+    del dissolved
+    # floor(xmin / CS) * CS etc., computed with integers only so there is no floating point drift.
+    x0 = math.floor(bbox.xMinimum()) // CS * CS
+    y0 = math.floor(bbox.yMinimum()) // CS * CS
+    x1 = -(-math.ceil(bbox.xMaximum()) // CS) * CS
+    y1 = -(-math.ceil(bbox.yMaximum()) // CS) * CS
+    LOG(f"\n[{city}] STEP 2a snap extent: boundary bbox ({bbox.xMinimum():.3f}, {bbox.yMinimum():.3f}) - "
+        f"({bbox.xMaximum():.3f}, {bbox.yMaximum():.3f}) -> ({x0}, {y0}) - ({x1}, {y1})")
+    if (x0, y0) != python_origin:
+        raise CheckFailed(f"snapped origin ({x0}, {y0}) differs from the Python origin {python_origin}")
+    LOG(f"  origin ({x0}, {y0}) = Python origin {python_origin}")
+
+    res_g, ng = run_alg(city, "STEP 2b create grid", "native:creategrid", {
+        "TYPE": 2, "EXTENT": f"{x0},{x1},{y0},{y1} [{CRS}]", "HSPACING": CS, "VSPACING": CS,
+        "HOVERLAY": 0, "VOVERLAY": 0, "CRS": CRS, "OUTPUT": step_path(city, "05_grid_full_extent")})
+    expected = ((x1 - x0) // CS) * ((y1 - y0) // CS)
+    if ng["OUTPUT"] != expected:
+        raise CheckFailed(f"create grid made {ng['OUTPUT']} cells, expected {expected}")
+
+    res_c, _ = run_alg(city, "STEP 2c clip grid to boundary", "native:clip", {
+        "INPUT": res_g["OUTPUT"], "OVERLAY": boundary_path(city), "OUTPUT": step_path(city, "06_clipped")})
+
+    res_e, ne = run_alg(city, "STEP 2d drop zero-area cells", "native:extractbyexpression", {
+        "INPUT": res_c["OUTPUT"], "EXPRESSION": "area($geometry) > 0",
+        "OUTPUT": step_path(city, "07_area_gt_0"), "FAIL_OUTPUT": step_path(city, "07_dropped_zero_area")})
+    LOG(f"  cells dropped (clipped area 0, touch the boundary only at a line or point): {ne['FAIL_OUTPUT']}")
+
+    full = CS * CS
+    run_alg(city, "STEP 2e grid schema (same as the Python grid)", "native:refactorfields", {
+        "INPUT": res_e["OUTPUT"],
+        "FIELDS_MAPPING": [
+            fmap("cell_id", TEXT, "concat(to_int(round(\"left\")), '_', to_int(round(\"bottom\")))"),
+            fmap("x_min", INT64, "to_int(round(\"left\"))"),
+            fmap("y_min", INT64, "to_int(round(\"bottom\"))"),
+            fmap("cell_size_m", INT64, str(CS)),
+            fmap("area_m2", DOUBLE, "area($geometry)"),
+            fmap("area_share", DOUBLE, f"area($geometry) / {full}"),
+            fmap("is_edge", BOOL, f"area($geometry) / {full} < {EDGE_SHARE_THRESHOLD}"),
+            fmap("city", TEXT, QgsExpression.quotedString(city)),
+            fmap("is_sliver", BOOL, f"area($geometry) / {full} < {SLIVER_SHARE_THRESHOLD}")],
+        "OUTPUT": grid_path(city)})
+    return (x0, y0), ne["FAIL_OUTPUT"]
+
+
+# ------------------------------------------------------------------------------- Step 3 and Step 4
+def check_city(city, origin):
+    """STEP 3: checks on the saved QGIS boundary and grid."""
+    b = read_boundary(boundary_path(city))
+    g = read_grid(grid_path(city))
+    cells = list(g["cells"].values())
+    edge = [c for c in cells if c["is_edge"]]
+    area_sum = sum(c["area_m2"] for c in cells)
+    r = {
+        "n_cells": g["n"], "n_full": len(cells) - len(edge), "n_edge": len(edge),
+        **{f"edge_lt_{round(t * 100)}pct": sum(c["area_share"] < t for c in edge) for t in SLIVER_REPORT_THRESHOLDS},
+        "n_sliver": sum(c["is_sliver"] for c in cells),
+        "area_sum": area_sum, "boundary_area": b["area_attr"],
+        "area_diff_pct": abs(area_sum - b["area_attr"]) / b["area_attr"] * 100,
+        "unique": g["duplicates"] == 0 and len(g["cells"]) == g["n"],
+        "multiples": all(c["x_min"] % CS == 0 and c["y_min"] % CS == 0 for c in cells),
+        "id_matches_corner": all(cid == f"{c['x_min']}_{c['y_min']}" for cid, c in g["cells"].items()),
+        "origin_ok": (min(c["x_min"] for c in cells), min(c["y_min"] for c in cells)) == origin,
+        "flags_ok": all(c["is_edge"] == (c["area_share"] < EDGE_SHARE_THRESHOLD)
+                        and c["is_sliver"] == (c["area_share"] < SLIVER_SHARE_THRESHOLD) for c in cells),
+        "crs_boundary": b["crs"], "crs_grid": g["crs"], "boundary_features": b["n"],
+    }
+    LOG(f"\n[{city}] STEP 3 checks")
+    LOG(f"  a. cells {r['n_cells']:,} = full {r['n_full']:,} + edge {r['n_edge']:,}; edge cells with area_share "
+        f"<5%: {r['edge_lt_5pct']}, <10%: {r['edge_lt_10pct']}, <25%: {r['edge_lt_25pct']}; is_sliver: {r['n_sliver']}")
+    LOG(f"  b. sum of area_m2 {area_sum:,.3f} m2 vs boundary {b['area_attr']:,.3f} m2: diff {r['area_diff_pct']:.2e} %")
+    LOG(f"  c. cell_id unique: {r['unique']}; equals '<x_min>_<y_min>': {r['id_matches_corner']}; "
+        f"x_min/y_min multiples of {CS}: {r['multiples']}; lowest corner = origin: {r['origin_ok']}; "
+        f"is_edge/is_sliver flags consistent: {r['flags_ok']}")
+    LOG(f"  d. CRS boundary file: {b['crs']}, grid file: {g['crs']}")
+    problems = []
+    if r["area_diff_pct"] >= MAX_AREA_DIFF_PCT:
+        problems.append(f"area difference {r['area_diff_pct']:.6f}% is not below {MAX_AREA_DIFF_PCT}%")
+    for key in ("unique", "multiples", "id_matches_corner", "origin_ok", "flags_ok"):
+        if not r[key]:
+            problems.append(f"check failed: {key}")
+    if b["crs"] != CRS or g["crs"] != CRS:
+        problems.append(f"saved CRS is boundary {b['crs']}, grid {g['crs']}, expected {CRS}")
+    if b["n"] != 1:
+        problems.append(f"boundary file has {b['n']} features, expected 1")
+    if problems:
+        raise CheckFailed("; ".join(problems))
+    LOG("  all STEP 3 checks passed")
+    return r, g, b
+
+
+def compare_with_python(city, label, g, b, origin):
+    """STEP 4: compare the QGIS grid and boundary with the Python outputs (read-only)."""
+    py_grid_file = os.path.join(PY_GRID_DIR, f"{city}_grid{CS}_{TAG}.gpkg")
+    py_bnd_file = os.path.join(PY_BND_DIR, f"{city}_boundary_{TAG}.gpkg")
+    p, pb = read_grid(py_grid_file), read_boundary(py_bnd_file)
+    q_ids, p_ids = set(g["cells"]), set(p["cells"])
+    common = sorted(q_ids & p_ids)
+    diffs = {cid: abs(g["cells"][cid]["area_m2"] - p["cells"][cid]["area_m2"]) for cid in common}
+    worst = max(diffs, key=diffs.get) if diffs else ""
+    p_cells = list(p["cells"].values())
+    row = {
+        "city": city, "city_name": label,
+        "qgis_cells": g["n"], "python_cells": p["n"],
+        "qgis_edge": sum(c["is_edge"] for c in g["cells"].values()), "python_edge": sum(c["is_edge"] for c in p_cells),
+        "qgis_slivers": sum(c["is_sliver"] for c in g["cells"].values()),
+        "python_slivers": sum(c["is_sliver"] for c in p_cells),
+        "cell_ids_identical": q_ids == p_ids,
+        "n_missing_in_qgis": len(p_ids - q_ids), "n_missing_in_python": len(q_ids - p_ids),
+        "missing_in_qgis": ";".join(sorted(p_ids - q_ids)), "missing_in_python": ";".join(sorted(q_ids - p_ids)),
+        "max_abs_area_diff_m2": diffs[worst] if diffs else float("nan"), "worst_cell_id": worst,
+        "edge_flag_mismatches": sum(g["cells"][c]["is_edge"] != p["cells"][c]["is_edge"] for c in common),
+        "sliver_flag_mismatches": sum(g["cells"][c]["is_sliver"] != p["cells"][c]["is_sliver"] for c in common),
+        "schema_identical": g["schema"] == p["schema"],
+        "qgis_grid_area_m2": sum(c["area_m2"] for c in g["cells"].values()),
+        "python_grid_area_m2": sum(c["area_m2"] for c in p_cells),
+        "qgis_boundary_area_m2": b["area_attr"], "python_boundary_area_m2": pb["area_attr"],
+        "boundary_area_diff_m2": b["area_attr"] - pb["area_attr"],
+        "qgis_origin": f"{origin[0]}_{origin[1]}",
+        "python_origin": f"{min(c['x_min'] for c in p_cells)}_{min(c['y_min'] for c in p_cells)}",
+    }
+    row["match"] = (row["cell_ids_identical"] and row["qgis_cells"] == row["python_cells"]
+                    and row["qgis_edge"] == row["python_edge"] and row["qgis_slivers"] == row["python_slivers"]
+                    and diffs != {} and diffs[worst] < MAX_CELL_AREA_DIFF_M2 and row["edge_flag_mismatches"] == 0
+                    and row["sliver_flag_mismatches"] == 0 and row["schema_identical"]
+                    and row["qgis_origin"] == row["python_origin"])
+    LOG(f"\n[{city}] STEP 4 comparison with the Python outputs (read-only)")
+    LOG(f"  cells QGIS {row['qgis_cells']:,} vs Python {row['python_cells']:,}; edge {row['qgis_edge']} vs "
+        f"{row['python_edge']}; slivers {row['qgis_slivers']} vs {row['python_slivers']}")
+    LOG(f"  cell_id sets identical: {row['cell_ids_identical']} (missing in QGIS: {row['n_missing_in_qgis']}, "
+        f"missing in Python: {row['n_missing_in_python']})")
+    if row["missing_in_qgis"] or row["missing_in_python"]:
+        LOG(f"    missing in QGIS: {row['missing_in_qgis'] or '-'}")
+        LOG(f"    missing in Python: {row['missing_in_python'] or '-'}")
+    LOG(f"  max |area_m2 difference| per cell: {row['max_abs_area_diff_m2']!r} m2 (cell {worst}); "
+        f"flag mismatches: is_edge {row['edge_flag_mismatches']}, is_sliver {row['sliver_flag_mismatches']}")
+    LOG(f"  schema identical (names, types, order): {row['schema_identical']}")
+    if not row["schema_identical"]:
+        LOG(f"    QGIS:   {g['schema']}\n    Python: {p['schema']}")
+    LOG(f"  boundary area QGIS {b['area_attr']:,.6f} vs Python {pb['area_attr']:,.6f} m2: "
+        f"diff {row['boundary_area_diff_m2']:.6f} m2")
+    LOG(f"  -> {'MATCH' if row['match'] else 'MISMATCH'}")
+    if city in EXPECTED_PYTHON_500 and CS == 500 and p["n"] != EXPECTED_PYTHON_500[city][2]:
+        LOG(f"  WARNING: Python grid has {p['n']} cells, task2_log.txt says {EXPECTED_PYTHON_500[city][2]}")
+    return row
+
+
+def python_origin(city):
+    """Grid origin of the existing Python grid (lowest cell corner), cross-checked with task2_log.txt."""
+    p = read_grid(os.path.join(PY_GRID_DIR, f"{city}_grid{CS}_{TAG}.gpkg"))
+    origin = (min(c["x_min"] for c in p["cells"].values()), min(c["y_min"] for c in p["cells"].values()))
+    if CS == 500 and city in EXPECTED_PYTHON_500 and origin != EXPECTED_PYTHON_500[city][:2]:
+        raise CheckFailed(f"Python grid origin {origin} differs from task2_log.txt {EXPECTED_PYTHON_500[city][:2]}")
+    return origin
+
+
+# ------------------------------------------------------------------------------------ Step 5 and 6
+def simple_fill(**props):
+    return QgsFillSymbol.createSimple({k: str(v) for k, v in props.items()})
+
+
+def style_boundary(layer):
+    symbol = simple_fill(style="no", outline_color="#0b0b0b", outline_width=0.8, outline_width_unit="MM")
+    layer.renderer().setSymbol(symbol)
+
+
+GRID_RULES = [  # (legend label, filter expression, fill properties)
+    ("Full cell", 'NOT "is_edge"', dict(style="no", outline_color="#8c8b87", outline_width=0.15)),
+    ("Edge cell", '"is_edge" AND NOT "is_sliver"',
+     dict(style="solid", color="#f7b98e", outline_color="#8c8b87", outline_width=0.15)),
+    ("Sliver, under 10% of a cell", '"is_sliver"',
+     dict(style="solid", color="#d62728", outline_color="#8c8b87", outline_width=0.15)),
+]
+
+
+def style_grid(layer):
+    root = QgsRuleBasedRenderer.Rule(None)
+    for label, expression, props in GRID_RULES:
+        root.appendChild(QgsRuleBasedRenderer.Rule(simple_fill(**props), 0, 0, expression, label))
+    layer.setRenderer(QgsRuleBasedRenderer(root))
+
+
+def rule_counts(layer):
+    return {label: sum(1 for _ in layer.getFeatures(QgsFeatureRequest(QgsExpression(expression))))
+            for label, expression, _ in GRID_RULES}
+
+
+def text_format(size, bold=False):
+    font = QFont("Arial")
+    font.setBold(bold)
+    fmt = QgsTextFormat()
+    fmt.setFont(font)
+    fmt.setSize(size)
+    fmt.setSizeUnit(Qgis.RenderUnit.Points)
+    return fmt
+
+
+def add_label(layout, text, size, x, y, w, h, bold=False):
+    label = QgsLayoutItemLabel(layout)
+    label.setText(text)
+    label.setTextFormat(text_format(size, bold))
+    layout.addLayoutItem(label)
+    label.attemptMove(QgsLayoutPoint(x, y))
+    label.attemptResize(QgsLayoutSize(w, h))
+    return label
+
+
+def build_layout(project, label, boundary_layer, grid_layer, basemap, png_path):
+    """STEP 6: A4 landscape layout with map, title, legend, scale bar (km), north arrow, source line."""
+    manager = project.layoutManager()
+    name = f"{label} grid QA"
+    old = manager.layoutByName(name)
+    if old is not None:
+        manager.removeLayout(old)
+    layout = QgsPrintLayout(project)
+    layout.initializeDefaults()
+    layout.setName(name)
+    layout.pageCollection().page(0).setPageSize("A4", QgsLayoutItemPage.Landscape)
+
+    add_label(layout, f"{label}: {CS} m analysis grid ({CRS})", 18, 10, 7, 277, 11, bold=True)
+
+    map_item = QgsLayoutItemMap(layout)
+    layout.addLayoutItem(map_item)
+    map_item.attemptMove(QgsLayoutPoint(10, 21))
+    map_item.attemptResize(QgsLayoutSize(205, 176))
+    map_item.setCrs(QgsCoordinateReferenceSystem(CRS))
+    map_item.setLayers([boundary_layer, grid_layer] + ([basemap] if basemap else []))
+    map_item.setKeepLayerSet(True)
+    extent = QgsRectangle(boundary_layer.extent())
+    extent.scale(1.06)  # small margin around the city
+    map_item.zoomToExtent(extent)
+    map_item.setFrameEnabled(True)
+
+    legend = QgsLayoutItemLegend(layout)
+    legend.setTitle("Legend")
+    legend.setLinkedMap(map_item)
+    layout.addLayoutItem(legend)
+    legend.setAutoUpdateModel(False)
+    tree = legend.model().rootGroup()
+    tree.removeAllChildren()
+    for layer, title in ((boundary_layer, "City boundary"), (grid_layer, f"Grid cells ({CS} m)")):
+        node = tree.addLayer(layer)
+        node.setCustomProperty("legend/title-label", title)
+        legend.model().refreshLayerLegend(node)
+    legend.updateLegend()
+    legend.attemptMove(QgsLayoutPoint(222, 21))
+
+    arrow = QgsLayoutItemPicture(layout)
+    arrow.setPicturePath(NORTH_ARROW_SVG)
+    layout.addLayoutItem(arrow)
+    arrow.setLinkedMap(map_item)
+    arrow.setNorthMode(QgsLayoutItemPicture.GridNorth)
+    arrow.attemptMove(QgsLayoutPoint(224, 128))
+    arrow.attemptResize(QgsLayoutSize(12, 18))
+
+    scalebar = QgsLayoutItemScaleBar(layout)
+    scalebar.setStyle("Single Box")
+    scalebar.setLinkedMap(map_item)
+    layout.addLayoutItem(scalebar)
+    scalebar.applyDefaultSize(Qgis.DistanceUnit.Kilometers)
+    scalebar.setUnitLabel("km")
+    scalebar.setNumberOfSegmentsLeft(0)
+    scalebar.attemptMove(QgsLayoutPoint(222, 155))
+    scalebar.update()
+
+    add_label(layout, SOURCE_LINE, 8, 10, 200, 277, 6)
+
+    manager.addLayout(layout)
+    settings = QgsLayoutExporter.ImageExportSettings()
+    settings.dpi = 200
+    # QGIS then tries to georeference the image in place; GDAL cannot update a PNG, so it prints a
+    # harmless "ERROR 6: The PNG driver does not support update access" message. The PNG is complete.
+    result = QgsLayoutExporter(layout).exportToImage(png_path, settings)
+    if result != QgsLayoutExporter.Success:
+        raise RuntimeError(f"layout export failed for {label} (code {result})")
+    LOG(f"  layout '{name}' saved in the project; exported {os.path.relpath(png_path, PROJECT_ROOT)} at 200 dpi")
+
+
+def build_project(done, project):
+    """STEP 5 and 6: team project with groups, styles, basemap, layouts; saved with relative paths."""
+    LOG(f"\n{'=' * 90}\nSTEP 5: team QGIS project\n{'=' * 90}")
+    project.setCrs(QgsCoordinateReferenceSystem(CRS))
+    LOG(f"project CRS: {project.crs().authid()}")
+    root = project.layerTreeRoot()
+    layers = {}
+    for city, info in done.items():
+        label = info["label"]
+        group = root.findGroup(label)
+        if group is None:
+            group = root.addGroup(label)
+        boundary = QgsVectorLayer(boundary_path(city), f"{label} boundary (QGIS)", "ogr")
+        grid = QgsVectorLayer(grid_path(city), f"{label} grid {CS} m (QGIS)", "ogr")
+        if not (boundary.isValid() and grid.isValid()):
+            raise RuntimeError(f"cannot load the QGIS outputs for {label}")
+        style_boundary(boundary)
+        style_grid(grid)
+        project.addMapLayer(boundary, False)
+        project.addMapLayer(grid, False)
+        group.addLayer(boundary)
+        group.addLayer(grid)
+        counts = rule_counts(grid)
+        expected = {"Full cell": info["checks"]["n_full"],
+                    "Edge cell": info["checks"]["n_edge"] - info["checks"]["n_sliver"],
+                    "Sliver, under 10% of a cell": info["checks"]["n_sliver"]}
+        LOG(f"  group '{label}': boundary + grid; rule counts {counts} "
+            f"({'match the checks' if counts == expected else 'DO NOT MATCH ' + str(expected)})")
+        if counts != expected:
+            raise CheckFailed(f"{label}: grid style rules do not cover the cells as expected")
+        layers[city] = (boundary, grid)
+
+    basemap = QgsRasterLayer(BASEMAP_URI, BASEMAP_NAME, "wms")
+    if basemap.isValid():
+        project.addMapLayer(basemap, False)
+        root.addLayer(basemap)  # appended last = bottom of the layer list
+        LOG(f"  basemap '{BASEMAP_NAME}' added at the bottom of the layer list")
+    else:
+        basemap = None
+        LOG("  WARNING: OpenStreetMap basemap could not be created; layouts are drawn without it")
+
+    LOG(f"\n{'=' * 90}\nSTEP 6: print layouts and QA maps\n{'=' * 90}")
+    for city, (boundary, grid) in layers.items():
+        build_layout(project, done[city]["label"], boundary, grid, basemap,
+                     os.path.join(QA_DIR, f"{city}_grid_qgis.png"))
+
+    project.setTitle("crimegeo_vgi: Task 2 boundaries and grids")
+    project.setFilePathStorage(Qgis.FilePathType.Relative)
+    if not project.write(PROJECT_FILE):
+        raise RuntimeError(f"could not save the project: {project.error()}")
+    with zipfile.ZipFile(PROJECT_FILE) as z:
+        qgs = z.read([n for n in z.namelist() if n.endswith(".qgs")][0]).decode("utf-8")
+    sources = re.findall(r"<datasource>([^<]*)</datasource>", qgs)
+    local = [s for s in sources if ".gpkg" in s]
+    absolute = [s for s in local if re.match(r"^[A-Za-z]:[\\/]|^/", s)]
+    LOG(f"\nproject saved: {os.path.relpath(PROJECT_FILE, PROJECT_ROOT)} ({len(local)} GeoPackage sources, "
+        f"{len(local) - len(absolute)} relative, {len(absolute)} absolute; "
+        f"{len(project.layoutManager().layouts())} layouts)")
+    if absolute:
+        raise CheckFailed(f"project has absolute data paths: {absolute}")
+    return layers
+
+
+# ------------------------------------------------------------------------------------------- main
+def main():
+    global LOG
+    LOG = Log(LOG_PATH)
+    CALLS.clear()
+    started = datetime.now()
+    project = QgsProject.instance()
+    try:
+        LOG(f"TASK 2 IN QGIS - run started {started:%Y-%m-%d %H:%M:%S}")
+        LOG(f"QGIS {Qgis.version()} | PROJECT_ROOT {PROJECT_ROOT}")
+        LOG(f"config: crs={CRS}, grid_size_m={CS}, boundary_file={CFG['boundary_file']}")
+        LOG("processing context: invalid feature filtering = do not filter "
+            "(validity is checked explicitly with native:checkvalidity)")
+
+        missing = [p for p in PROTECTED if not os.path.isfile(p)]
+        if missing:
+            raise CheckFailed(f"Python Task 2 outputs missing: {missing}")
+        hashes_before = {p: sha256(p) for p in PROTECTED}
+        LOG(f"recorded SHA-256 of {len(PROTECTED)} Python Task 2 files (they are only read)")
+
+        clear_previous_run(project)
+        remove_files([p for c in CFG["cities"] for p in city_outputs(c)]
+                     + [CSV_PATH, os.path.join(STEP_DIR, "00_lad_reprojected.gpkg")])
+        for d in (BND_DIR, GRID_DIR, QA_DIR):
+            os.makedirs(d, exist_ok=True)
+        for c in CFG["cities"]:
+            os.makedirs(os.path.join(STEP_DIR, c), exist_ok=True)
+
+        lad_source = inspect_lad_file()
+
+        done, failures, rows = {}, {}, []
+        for city, spec in CFG["cities"].items():
+            label = CFG["forces"][spec["force"]]
+            LOG(f"\n{'=' * 90}\n{label} ({city}), LAD code {spec['lad_code']}\n{'=' * 90}")
+            try:
+                origin_py = python_origin(city)
+                dissolved, repaired = build_boundary(city, spec, lad_source)
+                origin, dropped = build_grid(city, dissolved, origin_py)
+                checks, g, b = check_city(city, origin)
+                row = compare_with_python(city, label, g, b, origin)
+                row.update(dropped_zero_area=dropped, repaired=repaired)
+                rows.append(row)
+                done[city] = {"label": label, "checks": checks, "origin": origin, "row": row}
+            except Exception as e:  # one city failing must not hide the others
+                failures[city] = f"{type(e).__name__}: {e}"
+                LOG(f"\n[{city}] FAILED: {failures[city]}\n{traceback.format_exc()}")
+
+        if rows:
+            with open(CSV_PATH, "w", newline="", encoding="utf-8") as f:
+                writer = csv.DictWriter(f, fieldnames=list(rows[0].keys()))
+                writer.writeheader()
+                writer.writerows(rows)
+            LOG(f"\ncomparison written: {os.path.relpath(CSV_PATH, PROJECT_ROOT)}")
+
+        if failures:
+            LOG(f"\nSTOPPED: {len(failures)} city/cities failed, so the team project and layouts were NOT built:")
+            for city, message in failures.items():
+                LOG(f"  {city}: {message}")
+        else:
+            build_project(done, project)
+
+        hashes_after = {p: sha256(p) for p in PROTECTED}
+        changed = [os.path.relpath(p, PROJECT_ROOT) for p in PROTECTED if hashes_after[p] != hashes_before[p]]
+        LOG(f"\nPython Task 2 files unchanged: {not changed}" + (f" CHANGED: {changed}" if changed else ""))
+
+        LOG(f"\n{'=' * 90}\nSUMMARY\n{'=' * 90}")
+        LOG(f"{'city':<12}{'QGIS cells':>11}{'Python cells':>13}{'edge':>6}{'slivers':>8}{'ids same':>9}"
+            f"{'max area diff m2':>18}{'boundary diff m2':>18}  result")
+        for r in rows:
+            LOG(f"{r['city']:<12}{r['qgis_cells']:>11,}{r['python_cells']:>13,}{r['qgis_edge']:>6}{r['qgis_slivers']:>8}"
+                f"{'yes' if r['cell_ids_identical'] else 'NO':>9}{r['max_abs_area_diff_m2']:>18.3e}"
+                f"{r['boundary_area_diff_m2']:>18.3e}  {'MATCH' if r['match'] else 'MISMATCH'}")
+        LOG(f"\nalgorithm calls: {len(CALLS)}; IDs used in order (first city): "
+            + " -> ".join(dict.fromkeys(c["algorithm"] for c in CALLS if c["city"] == next(iter(CFG["cities"])))))
+        ok = not failures and not changed and all(r["match"] for r in rows)
+        LOG(f"\nTASK 2 IN QGIS: FINISHED in {(datetime.now() - started).seconds} s - "
+            + ("all checks passed, QGIS matches Python" if ok else "SEE PROBLEMS ABOVE"))
+
+        iface = globals().get("iface")
+        if iface is not None and done:
+            first = next(iter(done))
+            extent = QgsRectangle(read_boundary(boundary_path(first))["bbox"])
+            extent.scale(1.1)
+            iface.mapCanvas().setExtent(extent)
+            iface.mapCanvas().refresh()
+    except Exception as e:
+        LOG(f"\nSTOPPED: {type(e).__name__}: {e}\n{traceback.format_exc()}")
+    finally:
+        LOG.close()
+
+
+main()

task6_experiment.py

@@ -0,0 +1,634 @@
+"""Task 6: Incomplete reporting removal-experiment scaffold for urban crime analysis.
+
+This module provides a reproducible simulation scaffold to evaluate the stability
+of crime hotspot identification under incomplete reporting / underreporting.
+
+Key features:
+1. Accepts a city crime point layer and a 500 m grid, both in EPSG:27700.
+2. Assigns each point to exactly one grid cell using the full square's lower-left
+   coordinates (cell_id: "<x_min>_<y_min>").
+3. Simulates two forms of incomplete reporting:
+   - Uniform removal probability across all crime records.
+   - Supplied removal probability for each cell.
+4. Recounts crimes in every grid cell, explicitly preserving zero-crime cells.
+5. Identifies the top ceil(10% of cells) before and after removal, breaking count ties
+   deterministically by cell_id (ascending), and reports the fraction of original top
+   cells retained.
+6. Validates CRS (EPSG:27700), rates (0.0 to 1.0), grid IDs, and points outside the grid.
+7. Accepts a random seed for full reproducibility.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import math
+import sys
+from dataclasses import asdict, dataclass
+from pathlib import Path
+from typing import Any, Mapping, Sequence
+
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+import pyogrio
+import shapely
+from pyproj import CRS
+
+EXPECTED_CRS = "EPSG:27700"
+DEFAULT_CELL_SIZE = 500
+DEFAULT_TOP_FRACTION = 0.10
+
+
+@dataclass
+class ExperimentResult:
+    """Structured result of a Task 6 removal experiment."""
+
+    fraction_retained: float
+    n_top_cells: int
+    top_cells_before: list[str]
+    top_cells_after: list[str]
+    retained_top_cells: list[str]
+    lost_top_cells: list[str]
+    new_top_cells: list[str]
+    counts_before: pd.Series
+    counts_after: pd.Series
+    initial_crime_count: int
+    retained_crime_count: int
+    removed_crime_count: int
+    removal_type: str
+    seed: int | None
+    top_fraction: float
+    total_grid_cells: int
+
+    def summary(self) -> dict[str, Any]:
+        """Return a serializable dictionary summary of key metrics."""
+        return {
+            "removal_type": self.removal_type,
+            "seed": self.seed,
+            "top_fraction": self.top_fraction,
+            "total_grid_cells": self.total_grid_cells,
+            "n_top_cells": self.n_top_cells,
+            "initial_crime_count": self.initial_crime_count,
+            "retained_crime_count": self.retained_crime_count,
+            "removed_crime_count": self.removed_crime_count,
+            "crime_retention_rate": (
+                self.retained_crime_count / self.initial_crime_count
+                if self.initial_crime_count > 0
+                else 1.0
+            ),
+            "fraction_retained": self.fraction_retained,
+            "n_retained_top_cells": len(self.retained_top_cells),
+            "n_lost_top_cells": len(self.lost_top_cells),
+            "n_new_top_cells": len(self.new_top_cells),
+            "top_cells_before_head": self.top_cells_before[:10],
+            "top_cells_after_head": self.top_cells_after[:10],
+        }
+
+    def print_summary(self) -> None:
+        """Print an interpretable summary to standard output."""
+        s = self.summary()
+        print("\n" + "=" * 65)
+        print("TASK 6 REMOVAL EXPERIMENT SUMMARY")
+        print("=" * 65)
+        print(f"Removal type:                 {s['removal_type']}")
+        print(f"Random seed:                  {s['seed']}")
+        print(f"Total grid cells:             {s['total_grid_cells']:,}")
+        print(f"Top cells evaluated (k):      {s['n_top_cells']} (ceil {s['top_fraction']*100:.1f}%)")
+        print(f"Initial crimes:               {s['initial_crime_count']:,}")
+        print(f"Retained crimes:              {s['retained_crime_count']:,} ({s['crime_retention_rate']*100:.2f}%)")
+        print(f"Removed crimes:               {s['removed_crime_count']:,}")
+        print("-" * 65)
+        print(f"Fraction top cells retained:  {s['fraction_retained']:.4f} ({s['fraction_retained']*100:.2f}%)")
+        print(f"Retained top cells count:     {s['n_retained_top_cells']} / {s['n_top_cells']}")
+        print(f"Lost top cells count:         {s['n_lost_top_cells']}")
+        print(f"New entrant top cells count:  {s['n_new_top_cells']}")
+        print("=" * 65 + "\n")
+
+
+def validate_crs(gdf: gpd.GeoDataFrame, layer_name: str, expected_epsg: int = 27700) -> None:
+    """Validate that the GeoDataFrame has a valid projected CRS matching expected EPSG."""
+    if not isinstance(gdf, gpd.GeoDataFrame):
+        raise TypeError(f"{layer_name} must be a GeoDataFrame, got {type(gdf).__name__}")
+    if gdf.crs is None:
+        raise ValueError(f"{layer_name} has no CRS defined; expected EPSG:{expected_epsg}")
+
+    try:
+        crs_obj = CRS.from_user_input(gdf.crs)
+    except Exception as e:
+        raise ValueError(f"{layer_name} has invalid CRS '{gdf.crs}': {e}") from e
+
+    epsg = crs_obj.to_epsg()
+    if epsg != expected_epsg:
+        raise ValueError(
+            f"{layer_name} CRS must be EPSG:{expected_epsg}, but got '{crs_obj.to_string()}' (EPSG:{epsg})"
+        )
+
+
+def validate_grid(grid_gdf: gpd.GeoDataFrame, cell_size: int = DEFAULT_CELL_SIZE) -> None:
+    """Validate grid structure, unique cell_id values, and format."""
+    validate_crs(grid_gdf, "Grid layer")
+
+    if len(grid_gdf) == 0:
+        raise ValueError("Grid layer is empty (0 cells)")
+
+    if "cell_id" not in grid_gdf.columns:
+        raise ValueError("Grid layer must contain a 'cell_id' column")
+
+    cids = grid_gdf["cell_id"]
+    if cids.isna().any():
+        raise ValueError("Grid contains null/NaN 'cell_id' values")
+
+    if not cids.is_unique:
+        duplicates = cids[cids.duplicated()].unique().tolist()
+        raise ValueError(f"Grid 'cell_id' values must be unique; found duplicate(s): {duplicates[:5]}")
+
+    # Validate cell_id format "<x_min>_<y_min>" and alignment to cell_size
+    malformed = []
+    misaligned = []
+    for cid in cids:
+        parts = str(cid).split("_")
+        if len(parts) != 2:
+            malformed.append(cid)
+            continue
+        try:
+            x_val, y_val = int(parts[0]), int(parts[1])
+            if x_val % cell_size != 0 or y_val % cell_size != 0:
+                misaligned.append(cid)
+        except ValueError:
+            malformed.append(cid)
+
+    if malformed:
+        raise ValueError(f"Found malformed cell_id(s) not in '<x_min>_<y_min>' integer format: {malformed[:5]}")
+    if misaligned:
+        raise ValueError(
+            f"Found cell_id(s) where corner coordinates are not multiples of cell_size {cell_size}: {misaligned[:5]}"
+        )
+
+
+def validate_removal_rate(removal_rate: float) -> float:
+    """Validate that a uniform removal rate is a valid probability in [0.0, 1.0]."""
+    if isinstance(removal_rate, bool) or not isinstance(removal_rate, (int, float)):
+        raise TypeError(f"Removal rate must be a numeric float/int, got {type(removal_rate).__name__}")
+    if math.isnan(removal_rate) or math.isinf(removal_rate):
+        raise ValueError(f"Removal rate must be a finite number, got {removal_rate}")
+    if not (0.0 <= removal_rate <= 1.0):
+        raise ValueError(f"Removal rate must be between 0.0 and 1.0 inclusive, got {removal_rate}")
+    return float(removal_rate)
+
+
+def validate_cell_removal_rates(
+    cell_rates: Mapping[str, float] | pd.Series,
+    grid_cell_ids: Sequence[str] | set[str],
+) -> dict[str, float]:
+    """Validate cell-specific removal rates against the grid's cell IDs."""
+    if isinstance(cell_rates, pd.Series):
+        rates_dict = cell_rates.to_dict()
+    elif isinstance(cell_rates, Mapping):
+        rates_dict = dict(cell_rates)
+    else:
+        raise TypeError(
+            f"cell_rates must be a dict or pandas Series, got {type(cell_rates).__name__}"
+        )
+
+    grid_set = set(grid_cell_ids)
+    rates_set = set(rates_dict.keys())
+
+    # Check for invalid rate values
+    invalid_rates = {}
+    for cid, rate in rates_dict.items():
+        try:
+            rate_val = validate_removal_rate(rate)
+            rates_dict[cid] = rate_val
+        except (ValueError, TypeError) as e:
+            invalid_rates[cid] = (rate, str(e))
+
+    if invalid_rates:
+        sample = list(invalid_rates.items())[:3]
+        raise ValueError(f"Invalid removal rates in cell_rates: {sample}")
+
+    # Check coverage of grid cells
+    missing = grid_set - rates_set
+    if missing:
+        sample_missing = sorted(list(missing))[:5]
+        raise ValueError(
+            f"cell_rates is missing {len(missing)} grid cells (e.g., {sample_missing}). "
+            "Rates must be supplied for all grid cells."
+        )
+
+    # Check for extraneous cell IDs not present in the grid
+    extraneous = rates_set - grid_set
+    if extraneous:
+        sample_extra = sorted(list(extraneous))[:5]
+        raise ValueError(
+            f"cell_rates contains {len(extraneous)} cell IDs not present in the grid: {sample_extra}"
+        )
+
+    return rates_dict
+
+
+def compute_point_cell_id(
+    x: float | np.ndarray | Sequence[float],
+    y: float | np.ndarray | Sequence[float],
+    cell_size: int = DEFAULT_CELL_SIZE,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Assign point coordinates to full square lower-left corners and cell_id.
+
+    Following conventions.md and Task 2 grid rules:
+    x_min = floor(x / cell_size) * cell_size
+    y_min = floor(y / cell_size) * cell_size
+    cell_id = "<x_min>_<y_min>"
+
+    Returns (x_min, y_min, cell_id) as numpy arrays.
+    """
+    x_arr = np.atleast_1d(np.asarray(x, dtype=np.float64))
+    y_arr = np.atleast_1d(np.asarray(y, dtype=np.float64))
+
+    if x_arr.shape != y_arr.shape:
+        raise ValueError(f"x and y coordinate arrays must have matching shape, got {x_arr.shape} vs {y_arr.shape}")
+
+    if np.isnan(x_arr).any() or np.isnan(y_arr).any():
+        raise ValueError("Crime point coordinates contain NaN values")
+
+    # Math floor in integer arithmetic:
+    # Uses floor(coord / cell_size) * cell_size
+    x_min = np.floor(x_arr / cell_size).astype(np.int64) * cell_size
+    y_min = np.floor(y_arr / cell_size).astype(np.int64) * cell_size
+
+    cell_ids = np.array([f"{xi}_{yi}" for xi, yi in zip(x_min, y_min)], dtype=object)
+    return x_min, y_min, cell_ids
+
+
+# Alias for flexible import naming
+create_point_cell_id = compute_point_cell_id
+
+
+def assign_points_to_grid(
+    points_gdf: gpd.GeoDataFrame,
+    grid_gdf: gpd.GeoDataFrame,
+    cell_size: int = DEFAULT_CELL_SIZE,
+    allow_outside: bool = False,
+) -> gpd.GeoDataFrame:
+    """Assign each crime point to exactly one grid cell.
+
+    Validates CRS, computes the full square's lower-left corner, and checks
+    that all points fall inside the city grid cells.
+
+    If allow_outside is False and any point falls outside the grid, raises ValueError.
+    If allow_outside is True, outside points are dropped from the returned GeoDataFrame.
+    """
+    validate_crs(points_gdf, "Crime points layer")
+    validate_grid(grid_gdf, cell_size)
+
+    if len(points_gdf) == 0:
+        assigned = points_gdf.copy()
+        assigned["x_min"] = pd.Series(dtype=np.int64)
+        assigned["y_min"] = pd.Series(dtype=np.int64)
+        assigned["cell_id"] = pd.Series(dtype=str)
+        return assigned
+
+    # Validate geometries are valid Points
+    if points_gdf.geometry.isna().any() or points_gdf.geometry.is_empty.any():
+        raise ValueError("Crime points layer contains null or empty geometries")
+
+    geom_types = points_gdf.geometry.geom_type.unique()
+    if any(t != "Point" for t in geom_types):
+        raise ValueError(f"Crime layer must contain Point geometries only, found {geom_types}")
+
+    x_coords = points_gdf.geometry.x.to_numpy()
+    y_coords = points_gdf.geometry.y.to_numpy()
+    point_geoms = points_gdf.geometry.to_numpy()
+
+    x_min, y_min, cids = compute_point_cell_id(x_coords, y_coords, cell_size)
+
+    grid_geom_map = dict(zip(grid_gdf["cell_id"].astype(str), grid_gdf.geometry))
+
+    # 1. Check if the computed cell_id exists in the grid
+    in_grid_mask = np.array([cid in grid_geom_map for cid in cids], dtype=bool)
+
+    # 2. For points whose cell_id exists in the grid, verify intersection with the clipped cell geometry
+    intersects_mask = np.zeros(len(points_gdf), dtype=bool)
+    if in_grid_mask.any():
+        valid_indices = np.where(in_grid_mask)[0]
+        valid_point_geoms = point_geoms[valid_indices]
+        valid_cell_geoms = [grid_geom_map[cids[i]] for i in valid_indices]
+        intersects_mask[valid_indices] = shapely.intersects(
+            valid_point_geoms,
+            np.array(valid_cell_geoms, dtype=object),
+        )
+
+    # A point is truly valid if it maps to a grid cell AND intersects its clipped geometry
+    valid_mask = in_grid_mask & intersects_mask
+    n_outside = int((~valid_mask).sum())
+
+    if n_outside > 0:
+        n_outside_grid = int((~in_grid_mask).sum())
+        n_outside_clipped = int((in_grid_mask & ~intersects_mask).sum())
+        outside_sample_cids = cids[~valid_mask][:5].tolist()
+        msg = (
+            f"Found {n_outside} crime point(s) outside the city grid "
+            f"({n_outside_grid} outside grid boundaries, {n_outside_clipped} in cell square but outside clipped geometry; "
+            f"sample cell_ids: {outside_sample_cids})."
+        )
+        if not allow_outside:
+            raise ValueError(
+                f"{msg} Set allow_outside=True if you wish to drop points outside the city grid."
+            )
+
+    assigned = points_gdf.copy()
+    assigned["x_min"] = x_min
+    assigned["y_min"] = y_min
+    assigned["cell_id"] = cids.astype(str)
+
+    if n_outside > 0 and allow_outside:
+        assigned = assigned[valid_mask].copy()
+
+    return assigned
+
+
+def recount_crimes_per_cell(
+    grid_gdf: gpd.GeoDataFrame,
+    assigned_points_gdf: gpd.GeoDataFrame,
+) -> pd.Series:
+    """Recount crimes in every grid cell, including cells with zero crimes.
+
+    Returns a pandas Series indexed by cell_id (in the grid's original row order),
+    with integer counts >= 0 for every cell.
+    """
+    grid_cids = grid_gdf["cell_id"].astype(str).tolist()
+    counts = pd.Series(0, index=pd.Index(grid_cids, name="cell_id"), dtype=np.int64)
+
+    if len(assigned_points_gdf) > 0 and "cell_id" in assigned_points_gdf.columns:
+        point_counts = assigned_points_gdf["cell_id"].astype(str).value_counts()
+        counts = counts.add(point_counts, fill_value=0).astype(np.int64)
+
+    return counts
+
+
+def select_top_cells(
+    counts: pd.Series,
+    top_fraction: float = DEFAULT_TOP_FRACTION,
+) -> list[str]:
+    """Select the top ceil(top_fraction * total_cells) cells.
+
+    Ties in crime count are broken deterministically by cell_id in ascending order.
+    Returns a list of the top k cell_id strings.
+    """
+    if len(counts) == 0:
+        return []
+
+    if not (0.0 < top_fraction <= 1.0):
+        raise ValueError(f"top_fraction must be in (0.0, 1.0], got {top_fraction}")
+
+    total_cells = len(counts)
+    k = math.ceil(top_fraction * total_cells)
+    k = min(k, total_cells)
+
+    df = pd.DataFrame({
+        "cell_id": counts.index.astype(str),
+        "count": counts.values,
+    })
+
+    # Sort descending by count, then ascending by cell_id for deterministic tie-breaking
+    sorted_df = df.sort_values(
+        by=["count", "cell_id"],
+        ascending=[False, True],
+        kind="mergesort",  # stable sort
+    )
+
+    return sorted_df.head(k)["cell_id"].tolist()
+
+
+def calculate_retention_fraction(
+    top_cells_before: Sequence[str],
+    top_cells_after: Sequence[str],
+) -> float:
+    """Calculate the fraction of original top cells retained after removal."""
+    if len(top_cells_before) == 0:
+        return 1.0
+
+    orig_set = set(top_cells_before)
+    post_set = set(top_cells_after)
+    retained_set = orig_set.intersection(post_set)
+
+    return len(retained_set) / len(orig_set)
+
+
+def simulate_uniform_removal(
+    points_gdf: gpd.GeoDataFrame,
+    removal_rate: float,
+    seed: int | None = None,
+) -> gpd.GeoDataFrame:
+    """Simulate incomplete reporting with uniform removal probability per crime record."""
+    rate = validate_removal_rate(removal_rate)
+
+    if len(points_gdf) == 0:
+        return points_gdf.copy()
+
+    if rate == 0.0:
+        return points_gdf.copy()
+    if rate == 1.0:
+        return points_gdf.iloc[0:0].copy()
+
+    rng = np.random.default_rng(seed)
+    # Draw uniform random values u in [0, 1)
+    # A record is removed if u < rate; retained if u >= rate
+    u = rng.uniform(0.0, 1.0, size=len(points_gdf))
+    retained_mask = u >= rate
+    return points_gdf[retained_mask].copy()
+
+
+def simulate_cell_removal(
+    points_gdf: gpd.GeoDataFrame,
+    cell_rates: Mapping[str, float] | pd.Series,
+    seed: int | None = None,
+) -> gpd.GeoDataFrame:
+    """Simulate incomplete reporting with cell-specific removal probability."""
+    if len(points_gdf) == 0:
+        return points_gdf.copy()
+
+    if "cell_id" not in points_gdf.columns:
+        raise ValueError("points_gdf must have 'cell_id' column assigned before cell removal")
+
+    # Map each crime point's cell_id to its removal probability
+    point_cids = points_gdf["cell_id"].astype(str)
+    probs = point_cids.map(cell_rates).to_numpy(dtype=np.float64)
+
+    if np.isnan(probs).any():
+        missing_cids = point_cids[np.isnan(probs)].unique().tolist()
+        raise ValueError(f"Some points belong to cells missing from cell_rates: {missing_cids[:5]}")
+
+    rng = np.random.default_rng(seed)
+    u = rng.uniform(0.0, 1.0, size=len(points_gdf))
+    retained_mask = u >= probs
+    return points_gdf[retained_mask].copy()
+
+
+def run_removal_experiment(
+    points: gpd.GeoDataFrame | str | Path,
+    grid: gpd.GeoDataFrame | str | Path,
+    removal: float | Mapping[str, float] | pd.Series,
+    *,
+    seed: int | None = None,
+    cell_size: int = DEFAULT_CELL_SIZE,
+    top_fraction: float = DEFAULT_TOP_FRACTION,
+    allow_outside: bool = False,
+) -> ExperimentResult:
+    """Run a complete Task 6 removal experiment.
+
+    Parameters
+    ----------
+    points : GeoDataFrame or path to GPKG/spatial file
+        City crime point layer in EPSG:27700.
+    grid : GeoDataFrame or path to GPKG/spatial file
+        City grid layer (e.g. 500 m) in EPSG:27700 with 'cell_id'.
+    removal : float or Mapping/Series
+        Either a uniform removal probability (0.0 to 1.0) or
+        a mapping of cell_id -> removal probability for each cell.
+    seed : int or None, optional
+        Random seed for reproducibility.
+    cell_size : int, optional
+        Cell dimension in metres (default 500).
+    top_fraction : float, optional
+        Fraction of top hotspot cells to evaluate (default 0.10).
+    allow_outside : bool, optional
+        Whether to permit and drop points outside the city grid (default False).
+
+    Returns
+    -------
+    ExperimentResult
+        Structured experiment result including counts, top cells, and retention fraction.
+    """
+    # Load layers if paths were provided
+    if isinstance(points, (str, Path)):
+        points_gdf = gpd.read_file(points, engine="pyogrio")
+    else:
+        points_gdf = points
+
+    if isinstance(grid, (str, Path)):
+        grid_gdf = gpd.read_file(grid, engine="pyogrio")
+    else:
+        grid_gdf = grid
+
+    validate_crs(points_gdf, "Crime points layer")
+    validate_grid(grid_gdf, cell_size)
+
+    # 1. Assign points to grid cells
+    assigned_points = assign_points_to_grid(
+        points_gdf,
+        grid_gdf,
+        cell_size=cell_size,
+        allow_outside=allow_outside,
+    )
+
+    # 2. Recount initial crimes per cell (including zero-crime cells)
+    counts_before = recount_crimes_per_cell(grid_gdf, assigned_points)
+
+    # 3. Select top ceil(top_fraction * total_cells) before removal
+    top_before = select_top_cells(counts_before, top_fraction=top_fraction)
+
+    # 4. Simulate removal
+    if isinstance(removal, (int, float)) and not isinstance(removal, bool):
+        removal_type = "uniform"
+        validate_removal_rate(removal)
+        retained_points = simulate_uniform_removal(assigned_points, removal, seed=seed)
+    elif isinstance(removal, (Mapping, pd.Series)):
+        removal_type = "cell_specific"
+        validated_rates = validate_cell_removal_rates(removal, grid_gdf["cell_id"].tolist())
+        retained_points = simulate_cell_removal(assigned_points, validated_rates, seed=seed)
+    else:
+        raise TypeError(
+            f"removal must be float (uniform) or dict/Series (cell-specific), got {type(removal).__name__}"
+        )
+
+    # 5. Recount crimes per cell after removal (including zero-crime cells)
+    counts_after = recount_crimes_per_cell(grid_gdf, retained_points)
+
+    # 6. Select top cells after removal
+    top_after = select_top_cells(counts_after, top_fraction=top_fraction)
+
+    # 7. Compute retention fraction and cell sets
+    retention_fraction = calculate_retention_fraction(top_before, top_after)
+    set_before = set(top_before)
+    set_after = set(top_after)
+
+    retained_top = [cid for cid in top_before if cid in set_after]
+    lost_top = [cid for cid in top_before if cid not in set_after]
+    new_top = [cid for cid in top_after if cid not in set_before]
+
+    return ExperimentResult(
+        fraction_retained=retention_fraction,
+        n_top_cells=len(top_before),
+        top_cells_before=top_before,
+        top_cells_after=top_after,
+        retained_top_cells=retained_top,
+        lost_top_cells=lost_top,
+        new_top_cells=new_top,
+        counts_before=counts_before,
+        counts_after=counts_after,
+        initial_crime_count=len(assigned_points),
+        retained_crime_count=len(retained_points),
+        removed_crime_count=len(assigned_points) - len(retained_points),
+        removal_type=removal_type,
+        seed=seed,
+        top_fraction=top_fraction,
+        total_grid_cells=len(grid_gdf),
+    )
+
+
+def main() -> None:
+    """Command-line interface for running Task 6 removal experiments."""
+    parser = argparse.ArgumentParser(
+        description="Run Task 6 removal experiment on a city crime layer and grid."
+    )
+    parser.add_argument("--crime", type=Path, required=True, help="Path to crime point layer (GPKG)")
+    parser.add_argument("--grid", type=Path, required=True, help="Path to city 500m grid (GPKG)")
+    parser.add_argument("--uniform-rate", type=float, default=None, help="Uniform removal rate in [0.0, 1.0]")
+    parser.add_argument("--cell-rates", type=Path, default=None, help="Path to JSON/CSV with cell_id -> rate")
+    parser.add_argument("--seed", type=int, default=42, help="Random seed for reproducibility (default: 42)")
+    parser.add_argument("--cell-size", type=int, default=DEFAULT_CELL_SIZE, help="Cell size in metres (default: 500)")
+    parser.add_argument("--top-fraction", type=float, default=DEFAULT_TOP_FRACTION, help="Top hotspot fraction (default: 0.10)")
+    parser.add_argument("--allow-outside", action="store_true", help="Drop points outside grid instead of failing")
+    parser.add_argument("--output-json", type=Path, default=None, help="Optional output JSON path for results")
+
+    args = parser.parse_args()
+
+    if args.uniform_rate is None and args.cell_rates is None:
+        parser.error("Must specify either --uniform-rate or --cell-rates")
+    if args.uniform_rate is not None and args.cell_rates is not None:
+        parser.error("Specify either --uniform-rate or --cell-rates, not both")
+
+    if args.uniform_rate is not None:
+        removal = args.uniform_rate
+    else:
+        rates_path = args.cell_rates
+        if rates_path.suffix.lower() == ".json":
+            with open(rates_path, "r", encoding="utf-8") as f:
+                removal = json.load(f)
+        elif rates_path.suffix.lower() == ".csv":
+            df_rates = pd.read_csv(rates_path)
+            removal = dict(zip(df_rates["cell_id"].astype(str), df_rates["removal_rate"]))
+        else:
+            parser.error("Cell rates file must be .json or .csv")
+
+    res = run_removal_experiment(
+        points=args.crime,
+        grid=args.grid,
+        removal=removal,
+        seed=args.seed,
+        cell_size=args.cell_size,
+        top_fraction=args.top_fraction,
+        allow_outside=args.allow_outside,
+    )
+
+    res.print_summary()
+
+    if args.output_json:
+        summary_data = res.summary()
+        with open(args.output_json, "w", encoding="utf-8") as f:
+            json.dump(summary_data, f, indent=2)
+        print(f"Results saved to {args.output_json}")
+
+
+if __name__ == "__main__":
+    main()

test_task6_experiment.py

@@ -0,0 +1,462 @@
+"""Synthetic tests for Task 6 removal experiment scaffold.
+
+Tests all required edge cases and behaviors:
+1. Zero removal (retains 100% of crimes and original top cells).
+2. Complete removal (drops 100% of crimes, recount all 0, deterministic tie-breaking).
+3. Uneven removal changing the hotspot ranking.
+4. Boundary coordinate assignment (lower-left half-open interval rule).
+5. Validation of invalid rates (negative, > 1.0, NaN, missing cell rates).
+6. Repeatability and determinism with random seed.
+7. Validation of CRS (must be EPSG:27700).
+8. Validation of points falling outside the city grid.
+9. Validation of grid IDs (uniqueness, integer multiples of 500m).
+10. Explicit recount preserving zero-crime cells.
+11. Top ceil(10%) selection and deterministic tie-breaking.
+"""
+
+import math
+import unittest
+from typing import Sequence
+
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+import shapely
+from pyproj import CRS
+
+from task6_experiment import (
+    DEFAULT_CELL_SIZE,
+    EXPECTED_CRS,
+    assign_points_to_grid,
+    compute_point_cell_id,
+    recount_crimes_per_cell,
+    run_removal_experiment,
+    select_top_cells,
+    simulate_cell_removal,
+    simulate_uniform_removal,
+    validate_cell_removal_rates,
+    validate_crs,
+    validate_grid,
+    validate_removal_rate,
+)
+
+
+def create_synthetic_grid(
+    x_coords: Sequence[int],
+    y_coords: Sequence[int],
+    cell_size: int = DEFAULT_CELL_SIZE,
+    crs: str = EXPECTED_CRS,
+) -> gpd.GeoDataFrame:
+    """Helper to generate a synthetic square grid conforming to Task 2 structure."""
+    rows = []
+    for x in x_coords:
+        for y in y_coords:
+            cid = f"{x}_{y}"
+            poly = shapely.box(x, y, x + cell_size, y + cell_size)
+            rows.append({
+                "cell_id": cid,
+                "x_min": x,
+                "y_min": y,
+                "cell_size_m": cell_size,
+                "area_m2": float(cell_size * cell_size),
+                "area_share": 1.0,
+                "is_edge": False,
+                "city": "synthetic_city",
+                "is_sliver": False,
+                "geometry": poly,
+            })
+    return gpd.GeoDataFrame(rows, crs=crs)
+
+
+def create_synthetic_points(
+    coords: Sequence[tuple[float, float]],
+    crs: str = EXPECTED_CRS,
+) -> gpd.GeoDataFrame:
+    """Helper to generate a synthetic crime points GeoDataFrame."""
+    geoms = [shapely.Point(x, y) for x, y in coords]
+    return gpd.GeoDataFrame({"geometry": geoms}, crs=crs)
+
+
+class TestTask6Experiment(unittest.TestCase):
+    """Test suite covering Task 6 removal experiment scaffold."""
+
+    def setUp(self) -> None:
+        """Set up standard synthetic 2x5 grid (10 cells, ceil(10% * 10) = 1 top cell)."""
+        # 10 cells in total: xs=[400000, 400500], ys=[500000, 500500, 501000, 501500, 502000]
+        self.xs = [400000, 400500]
+        self.ys = [500000, 500500, 501000, 501500, 502000]
+        self.grid = create_synthetic_grid(self.xs, self.ys)
+        self.total_cells = len(self.grid)  # 10 cells
+
+    def test_zero_removal(self) -> None:
+        """Requirement: Zero removal (p=0.0) retains all crimes and 100% of top cells."""
+        # Place 10 points in cell 400000_500000, 5 points in 400500_500000
+        coords = [(400100.0, 500100.0)] * 10 + [(400600.0, 500100.0)] * 5
+        points = create_synthetic_points(coords)
+
+        res = run_removal_experiment(
+            points=points,
+            grid=self.grid,
+            removal=0.0,
+            seed=42,
+        )
+
+        self.assertEqual(res.initial_crime_count, 15)
+        self.assertEqual(res.retained_crime_count, 15)
+        self.assertEqual(res.removed_crime_count, 0)
+        self.assertEqual(res.fraction_retained, 1.0)
+        self.assertEqual(res.n_top_cells, 1)  # ceil(10% of 10) = 1
+        self.assertEqual(res.top_cells_before, ["400000_500000"])
+        self.assertEqual(res.top_cells_after, ["400000_500000"])
+        self.assertEqual(res.lost_top_cells, [])
+        self.assertTrue((res.counts_before == res.counts_after).all())
+
+    def test_complete_removal(self) -> None:
+        """Requirement: Complete removal (p=1.0) drops all crimes, recounts all 0, tie-breaks deterministically."""
+        coords = [(400100.0, 500100.0)] * 20 + [(400600.0, 502100.0)] * 10
+        points = create_synthetic_points(coords)
+
+        res = run_removal_experiment(
+            points=points,
+            grid=self.grid,
+            removal=1.0,
+            seed=42,
+        )
+
+        self.assertEqual(res.initial_crime_count, 30)
+        self.assertEqual(res.retained_crime_count, 0)
+        self.assertEqual(res.removed_crime_count, 30)
+        # All post-removal counts are zero
+        self.assertEqual(int(res.counts_after.sum()), 0)
+        self.assertTrue((res.counts_after == 0).all())
+
+        # Prior top cell was 400000_500000 (20 crimes)
+        self.assertEqual(res.top_cells_before, ["400000_500000"])
+        # Post removal: all cells tied at 0, tie-broken by cell_id ascending:
+        # First cell_id alphabetically is "400000_500000"
+        self.assertEqual(res.top_cells_after, ["400000_500000"])
+        self.assertEqual(res.fraction_retained, 1.0)
+
+        # Now test where original top cell is NOT the first alphabetically
+        coords_alt = [(400500.0, 502100.0)] * 20  # Cell 400500_502000 has all crimes
+        points_alt = create_synthetic_points(coords_alt)
+        res_alt = run_removal_experiment(
+            points=points_alt,
+            grid=self.grid,
+            removal=1.0,
+            seed=42,
+        )
+        self.assertEqual(res_alt.top_cells_before, ["400500_502000"])
+        # All 0 count, first alphabetically is "400000_500000"
+        self.assertEqual(res_alt.top_cells_after, ["400000_500000"])
+        self.assertEqual(res_alt.fraction_retained, 0.0)
+        self.assertEqual(res_alt.lost_top_cells, ["400500_502000"])
+        self.assertEqual(res_alt.new_top_cells, ["400000_500000"])
+
+    def test_uneven_removal_changes_ranking(self) -> None:
+        """Requirement: Uneven removal changes cell ranking and reduces retention fraction."""
+        # Initial: Cell A (400000_500000) has 50 crimes, Cell B (400500_500000) has 20 crimes.
+        # Top 1 cell is Cell A.
+        coords = [(400100.0, 500100.0)] * 50 + [(400600.0, 500100.0)] * 20
+        points = create_synthetic_points(coords)
+
+        # Cell A gets 95% removal probability, Cell B gets 0% removal probability.
+        # All other cells get 0% removal.
+        rates = {cid: 0.0 for cid in self.grid["cell_id"]}
+        rates["400000_500000"] = 0.95
+        rates["400500_500000"] = 0.0
+
+        res = run_removal_experiment(
+            points=points,
+            grid=self.grid,
+            removal=rates,
+            seed=42,
+        )
+
+        self.assertEqual(res.top_cells_before, ["400000_500000"])
+        # Cell B retained all 20 crimes, Cell A retained ~2-3 crimes
+        self.assertEqual(res.counts_after["400500_500000"], 20)
+        self.assertLess(res.counts_after["400000_500000"], 10)
+        # Cell B is now top cell!
+        self.assertEqual(res.top_cells_after, ["400500_500000"])
+        self.assertEqual(res.fraction_retained, 0.0)
+        self.assertEqual(res.lost_top_cells, ["400000_500000"])
+        self.assertEqual(res.new_top_cells, ["400500_500000"])
+
+    def test_boundary_assignment(self) -> None:
+        """Requirement: Boundary assignment using lower-left half-open interval rule."""
+        # A cell with lower-left (400000, 500000) spans [400000, 400500) x [500000, 500500)
+        # Test exact boundary coordinates:
+        # 1. Exact lower-left vertex (400000.0, 500000.0) -> must map to 400000_500000
+        # 2. Left edge (400000.0, 500250.0) -> must map to 400000_500000
+        # 3. Bottom edge (400250.0, 500000.0) -> must map to 400000_500000
+        # 4. Interior point (400499.99, 500499.99) -> must map to 400000_500000
+        # 5. Right boundary line (400500.0, 500250.0) -> must map to 400500_500000 (next cell east)
+        # 6. Top boundary line (400250.0, 500500.0) -> must map to 400000_500500 (next cell north)
+        # 7. Top-right vertex (400500.0, 500500.0) -> must map to 400500_500500 (diagonal cell)
+
+        test_cases = [
+            ((400000.0, 500000.0), "400000_500000"),
+            ((400000.0, 500250.0), "400000_500000"),
+            ((400250.0, 500000.0), "400000_500000"),
+            ((400499.99, 500499.99), "400000_500000"),
+            ((400500.0, 500250.0), "400500_500000"),
+            ((400250.0, 500500.0), "400000_500500"),
+            ((400500.0, 500500.0), "400500_500500"),
+        ]
+
+        for (x, y), expected_cid in test_cases:
+            x_min, y_min, cid = compute_point_cell_id(x, y, cell_size=500)
+            self.assertEqual(cid[0], expected_cid, f"Point ({x}, {y}) mapped to {cid[0]}, expected {expected_cid}")
+
+        # Check in GeoDataFrame assignment
+        points = create_synthetic_points([pt for pt, _ in test_cases])
+        assigned = assign_points_to_grid(points, self.grid, cell_size=500)
+        self.assertEqual(len(assigned), len(test_cases))
+        for i, (_, expected_cid) in enumerate(test_cases):
+            self.assertEqual(assigned["cell_id"].iloc[i], expected_cid)
+
+    def test_invalid_rates(self) -> None:
+        """Requirement: Validate removal rates (reject negative, >1.0, non-numeric, incomplete)."""
+        # Negative rate
+        with self.assertRaises(ValueError):
+            validate_removal_rate(-0.1)
+
+        # Rate > 1.0
+        with self.assertRaises(ValueError):
+            validate_removal_rate(1.0001)
+
+        # NaN and inf
+        with self.assertRaises(ValueError):
+            validate_removal_rate(float("nan"))
+        with self.assertRaises(ValueError):
+            validate_removal_rate(float("inf"))
+
+        # Non-numeric
+        with self.assertRaises(TypeError):
+            validate_removal_rate("0.5")  # type: ignore
+        with self.assertRaises(TypeError):
+            validate_removal_rate(True)  # bool is not allowed
+
+        # Incomplete cell rates dictionary
+        partial_rates = {"400000_500000": 0.5}  # missing other 9 cells
+        with self.assertRaises(ValueError):
+            validate_cell_removal_rates(partial_rates, self.grid["cell_id"].tolist())
+
+        # Cell rates with invalid value
+        bad_rates = {cid: 0.1 for cid in self.grid["cell_id"]}
+        bad_rates["400000_500000"] = 1.5
+        with self.assertRaises(ValueError):
+            validate_cell_removal_rates(bad_rates, self.grid["cell_id"].tolist())
+
+        # Cell rates with extra non-grid cell_id
+        extra_rates = {cid: 0.1 for cid in self.grid["cell_id"]}
+        extra_rates["999000_999000"] = 0.2
+        with self.assertRaises(ValueError):
+            validate_cell_removal_rates(extra_rates, self.grid["cell_id"].tolist())
+
+    def test_repeatability(self) -> None:
+        """Requirement: Accept random seed and guarantee reproducible results."""
+        # 100 crimes scattered across cells
+        rng = np.random.default_rng(999)
+        xs = rng.uniform(400000, 401000, size=100)
+        ys = rng.uniform(500000, 502500, size=100)
+        points = create_synthetic_points(list(zip(xs, ys)))
+
+        res1 = run_removal_experiment(points=points, grid=self.grid, removal=0.4, seed=12345)
+        res2 = run_removal_experiment(points=points, grid=self.grid, removal=0.4, seed=12345)
+        res3 = run_removal_experiment(points=points, grid=self.grid, removal=0.4, seed=54321)
+
+        # Identical seed -> exactly identical results
+        self.assertEqual(res1.retained_crime_count, res2.retained_crime_count)
+        self.assertEqual(res1.fraction_retained, res2.fraction_retained)
+        self.assertEqual(res1.top_cells_after, res2.top_cells_after)
+        self.assertTrue((res1.counts_after == res2.counts_after).all())
+
+        # Different seed -> distinct random simulation draw
+        # (With 100 points and p=0.4, exact same retained counts vector across different seed is infinitesimal)
+        counts_differ = not (res1.counts_after == res3.counts_after).all()
+        self.assertTrue(counts_differ, "Different seeds should produce different random removal draws")
+
+    def test_crs_validation(self) -> None:
+        """Requirement: Validate CRS is EPSG:27700 on both layers."""
+        coords = [(400100.0, 500100.0)]
+        pts_wgs84 = create_synthetic_points(coords, crs="EPSG:4326")
+        grid_27700 = self.grid
+
+        with self.assertRaises(ValueError) as ctx:
+            run_removal_experiment(points=pts_wgs84, grid=grid_27700, removal=0.1)
+        self.assertIn("EPSG:27700", str(ctx.exception))
+
+        # Test missing CRS on grid
+        grid_no_crs = grid_27700.set_crs(None, allow_override=True)
+        pts_27700 = create_synthetic_points(coords, crs="EPSG:27700")
+
+        with self.assertRaises(ValueError) as ctx2:
+            run_removal_experiment(points=pts_27700, grid=grid_no_crs, removal=0.1)
+        self.assertIn("no CRS", str(ctx2.exception))
+
+    def test_points_outside_city_grid(self) -> None:
+        """Requirement: Validate points falling outside the city grid."""
+        # Cell outside synthetic grid: 450000, 600000
+        outside_coords = [(450100.0, 600100.0), (400100.0, 500100.0)]
+        points = create_synthetic_points(outside_coords)
+
+        # Strict validation (allow_outside=False) must raise ValueError
+        with self.assertRaises(ValueError) as ctx:
+            run_removal_experiment(points=points, grid=self.grid, removal=0.1, allow_outside=False)
+        self.assertIn("outside the city grid", str(ctx.exception))
+
+        # Permissive validation (allow_outside=True) drops the outside point
+        res = run_removal_experiment(points=points, grid=self.grid, removal=0.0, allow_outside=True)
+        self.assertEqual(res.initial_crime_count, 1)  # Only 1 point inside
+
+    def test_grid_id_validation(self) -> None:
+        """Requirement: Validate grid IDs (uniqueness, correct lower-left format)."""
+        # Duplicate cell_id
+        bad_grid = self.grid.copy()
+        bad_grid.loc[1, "cell_id"] = bad_grid.loc[0, "cell_id"]
+        with self.assertRaises(ValueError) as ctx:
+            validate_grid(bad_grid)
+        self.assertIn("unique", str(ctx.exception))
+
+        # Malformed cell_id
+        bad_grid2 = self.grid.copy()
+        bad_grid2.loc[0, "cell_id"] = "cell_one"
+        with self.assertRaises(ValueError) as ctx2:
+            validate_grid(bad_grid2)
+        self.assertIn("malformed", str(ctx2.exception))
+
+        # Misaligned cell_id (not multiple of 500)
+        bad_grid3 = self.grid.copy()
+        bad_grid3.loc[0, "cell_id"] = "400123_500000"
+        with self.assertRaises(ValueError) as ctx3:
+            validate_grid(bad_grid3)
+        self.assertIn("multiples of cell_size", str(ctx3.exception))
+
+    def test_recount_includes_zero_crime_cells(self) -> None:
+        """Requirement: Recount crimes in every grid cell, including cells with zero crimes."""
+        # Only put crimes in 1 of the 10 cells
+        coords = [(400100.0, 500100.0)] * 5
+        points = create_synthetic_points(coords)
+
+        counts = recount_crimes_per_cell(self.grid, assign_points_to_grid(points, self.grid))
+
+        # All 10 cells must be present in the recount series
+        self.assertEqual(len(counts), 10)
+        self.assertEqual(counts["400000_500000"], 5)
+        # Remaining 9 cells must be exactly 0
+        zero_cells = counts[counts == 0]
+        self.assertEqual(len(zero_cells), 9)
+
+    def test_top_cells_selection_and_deterministic_ties(self) -> None:
+        """Requirement: Top ceil(10% of cells) selection and deterministic tie-breaking by cell_id."""
+        # Create grid with 15 cells -> ceil(0.10 * 15) = ceil(1.5) = 2 top cells
+        xs = [400000, 400500, 401000]
+        ys = [500000, 500500, 501000, 501500, 502000]
+        grid15 = create_synthetic_grid(xs, ys)
+        self.assertEqual(len(grid15), 15)
+
+        # All cells have count 0
+        counts_zero = pd.Series(0, index=grid15["cell_id"])
+        top_2_zero = select_top_cells(counts_zero, top_fraction=0.10)
+        self.assertEqual(len(top_2_zero), 2)
+        # Tied at 0, must pick lowest cell_id alphabetically
+        expected_sorted = sorted(grid15["cell_id"])[:2]
+        self.assertEqual(top_2_zero, expected_sorted)
+
+        # Now set cell counts: one clear top cell (count 10), and three tied cells with count 5
+        counts_custom = counts_zero.copy()
+        counts_custom["401000_502000"] = 10  # Highest
+        # Tied at 5:
+        counts_custom["400500_501000"] = 5
+        counts_custom["400000_502000"] = 5
+        counts_custom["400500_500500"] = 5
+
+        # We want top 2: #1 is 401000_502000 (10).
+        # #2 must be chosen from the three tied at 5 by cell_id ascending:
+        # "400000_502000" vs "400500_500500" vs "400500_501000"
+        # Alphabetically, "400000_502000" < "400500_500500" < "400500_501000".
+        # So top 2 must be ["401000_502000", "400000_502000"]
+        top_2_custom = select_top_cells(counts_custom, top_fraction=0.10)
+        self.assertEqual(top_2_custom, ["401000_502000", "400000_502000"])
+
+    def test_partially_clipped_edge_cell_validation(self) -> None:
+        """Requirement: Point in full 500m square but outside clipped city polygon must be rejected or dropped.
+
+        Verifies that:
+        1. A point inside the clipped geometry of an edge cell is accepted.
+        2. A point inside the full 500m square of that same cell, but outside the clipped polygon,
+           is rejected with ValueError when allow_outside=False.
+        3. The outside point is counted and dropped when allow_outside=True.
+        """
+        # Create a synthetic grid with 1 full cell and 1 partially clipped edge cell
+        # Cell 1 (full cell): [400500, 401000] x [500000, 500500]
+        # Cell 2 (edge cell): [400000, 400500] x [500000, 500500], but clipped to y <= 500250 (bottom half)
+        edge_poly = shapely.box(400000.0, 500000.0, 400500.0, 500250.0)
+        full_poly = shapely.box(400500.0, 500000.0, 401000.0, 500500.0)
+
+        grid_df = gpd.GeoDataFrame(
+            [
+                {
+                    "cell_id": "400000_500000",
+                    "x_min": 400000,
+                    "y_min": 500000,
+                    "cell_size_m": 500,
+                    "area_m2": 125000.0,
+                    "area_share": 0.5,
+                    "is_edge": True,
+                    "city": "synthetic_city",
+                    "is_sliver": False,
+                    "geometry": edge_poly,
+                },
+                {
+                    "cell_id": "400500_500000",
+                    "x_min": 400500,
+                    "y_min": 500000,
+                    "cell_size_m": 500,
+                    "area_m2": 250000.0,
+                    "area_share": 1.0,
+                    "is_edge": False,
+                    "city": "synthetic_city",
+                    "is_sliver": False,
+                    "geometry": full_poly,
+                },
+            ],
+            crs=EXPECTED_CRS,
+        )
+
+        # Point inside the clipped polygon of the edge cell
+        pt_inside_edge = (400250.0, 500100.0)
+        # Point in the full 500m square of the edge cell, but in the unclipped area (outside clipped polygon)
+        pt_outside_edge = (400250.0, 500400.0)
+
+        pts = create_synthetic_points([pt_inside_edge, pt_outside_edge])
+
+        # 1. When allow_outside=False: must raise ValueError detailing the outside clipped point
+        with self.assertRaises(ValueError) as ctx:
+            assign_points_to_grid(pts, grid_df, cell_size=500, allow_outside=False)
+        self.assertIn("outside clipped geometry", str(ctx.exception))
+
+        # 2. When allow_outside=True: must drop the outside point and keep the inside point
+        assigned = assign_points_to_grid(pts, grid_df, cell_size=500, allow_outside=True)
+        self.assertEqual(len(assigned), 1)
+        self.assertEqual(assigned["cell_id"].iloc[0], "400000_500000")
+        self.assertEqual(assigned.geometry.iloc[0].x, pt_inside_edge[0])
+        self.assertEqual(assigned.geometry.iloc[0].y, pt_inside_edge[1])
+
+        # 3. Test through run_removal_experiment
+        res = run_removal_experiment(
+            points=pts,
+            grid=grid_df,
+            removal=0.0,
+            allow_outside=True,
+            top_fraction=0.50,
+        )
+        self.assertEqual(res.initial_crime_count, 1)
+        self.assertEqual(res.counts_before["400000_500000"], 1)
+        self.assertEqual(res.counts_before["400500_500000"], 0)
+
+
+if __name__ == "__main__":
+    unittest.main()