# AI habitat classification with Python: shared data foundation

This folder is the common base for every part of the TechnicalGIS series **"AI habitat classification with
Python"** (series id `ai-habitat-classification`). Every part uses the **same study area, grid, labels,
class scheme and train/test split**, so their accuracies can be compared directly.

| File | What it is |
| --- | --- |
| `habitat_common.py` | The shared module: constants, grid, labels, samples, Sentinel-2 and Sentinel-1 helpers |
| `prepare_data.py` | Downloads the labels, (re)draws the samples, optionally prebuilds composites and time series |
| `samples/habitat_samples.csv` | **Committed** sample design: 9,000 points, 1,000 per class, with split and fold |
| `samples/blocks.csv` | **Committed** 90 spatial blocks: `block_id`, `split`, `fold` (geometry from `hc.load_blocks()`) |
| `requirements.txt` | Pinned environment (Python 3.12) |
| `data/` | Download and composite cache (gitignored; rebuilt on demand) |

## Set-up

```powershell
# from the repository root
py -3.12 -m venv .venv-habitat
.venv-habitat\Scripts\python -m pip install -r tutorial-code\ai-habitat-classification-common\requirements.txt
cd tutorial-code\ai-habitat-classification-common
..\..\.venv-habitat\Scripts\python prepare_data.py                # labels; samples are already committed
..\..\.venv-habitat\Scripts\python prepare_data.py --composites   # + 2023 seasonal S2 and S1 composites (~30-40 s each after the cloud masks)
..\..\.venv-habitat\Scripts\python prepare_data.py --timeseries   # + per-date S2 reflectance at every sample (~4 min)
```

No logins or API keys are needed. Labels come from Natural England's public ArcGIS service, imagery from the
Microsoft Planetary Computer (anonymous, signed URLs). If PROJ or GDAL environment variables from another
install are set (`PROJ_LIB`, `PROJ_DATA`, `GDAL_DATA`), unset them before running.

Use it from a tutorial folder next to this one:

```python
import sys
from pathlib import Path
sys.path.insert(0, str(Path(__file__).resolve().parent.parent / "ai-habitat-classification-common"))
import habitat_common as hc
```

## Study area and grid

| Item | Value |
| --- | --- |
| Name | Dark Peak: Kinder Scout, Bleaklow, Edale, the Hope Valley, Glossop and the Derwent reservoirs |
| CRS | `hc.CRS` = `EPSG:27700` (British National Grid) |
| Bounds | `hc.AOI_BOUNDS` = `(400000, 377500, 425000, 400000)` (xmin, ymin, xmax, ymax), 25 km x 22.5 km |
| Resolution | `hc.RES` = 10 m, north-up, pixel edges on the 10 m BNG lattice |
| Shape | `hc.shape()` = `(2250, 2500)` (rows, cols) |
| Grid objects | `hc.geobox()` (odc-geo GeoBox), `hc.transform()` (affine), `hc.aoi_bbox_wgs84()` |
| Year | `hc.YEAR` = 2023 (matches the Living England 2022-23 epoch) |
| Seed | `hc.SEED` = 42 |

Every raster produced by the helpers is a `(rows, cols)` array aligned to this grid, so `row`/`col` in the
sample table index straight into it. Write rasters with `hc.write_geotiff(path, arr, nodata, dtype)`.

## Labels and class scheme

Labels are the **Living England Habitat Map 2022-23** (Natural England; the latest release, the successor to
Phase 4), from the public FeatureServer
`https://services.arcgis.com/JJzESW51TqeY9uat/ArcGIS/rest/services/Living_England_2022_2023/FeatureServer`
(zone layers 3, 5 and 6). It is a machine-learning product (Sentinel-1/2 plus ancillary data, trained on field
data) with a per-polygon reliability rating, so treat "accuracy" as **agreement with Living England**, not with
the ground. 99.6% of the AOI is labelled; 68.8% is High or Very High reliability.

UKCEH Land Cover Map was not used (download needs a login). The Priority Habitat Inventory is reachable
(`Priority_Habitat_Inventory_North_England_EMP`) but covers priority habitats only, not wall-to-wall.

| id | `key` | Name | Living England `Prmry_H` values | Share of AOI |
| --- | --- | --- | --- | --- |
| 1 | `bog` | Blanket bog | Bog | 16.1% |
| 2 | `heath` | Heath | Dwarf Shrub Heath | 26.9% |
| 3 | `bracken` | Bracken | Bracken | 7.4% |
| 4 | `seminat_grass` | Semi-natural grassland | Unimproved Grassland | 17.1% |
| 5 | `improved_grass` | Improved grassland | Improved and Semi-Improved Grassland; Arable and Horticultural | 16.0% |
| 6 | `broadleaf` | Broadleaf woodland | Broadleaved, Mixed and Yew Woodland; Scrub | 7.1% |
| 7 | `conifer` | Conifer woodland | Coniferous Woodland | 3.1% |
| 8 | `water` | Water | Water | 1.5% |
| 9 | `built_bare` | Built-up and bare | Built-up Areas and Gardens; Bare Ground; Bare Sand | 4.6% |

`Fen, Marsh and Swamp` (about 2 km²) is left unlabelled (0). Semi-natural grassland is acid grassland on the
gritstone moors and calcareous grassland on the limestone south of Castleton; Living England does not split them.
Colours, names and the mapping live in `hc.CLASSES` (list of dicts with `id`, `key`, `name`, `colour`, `le`);
also `hc.CLASS_IDS`, `hc.CLASS_NAMES`, `hc.CLASS_KEYS`, `hc.LE_TO_ID`, and `hc.class_cmap()` for plotting.

| Function | Returns |
| --- | --- |
| `hc.load_labels()` | GeoDataFrame of Living England polygons clipped to the AOI: `class_id`, `class_key`, `le_habitat`, `le_reliability`, `le_source`, ... (downloads once to `data/living_england_aoi.gpkg`) |
| `hc.label_raster()` | uint8 `(2250, 2500)` label map, 0 = unlabelled (cached `data/labels_10m.tif`) |
| `hc.rasterize_labels(gdf, reliable_only=False)` | Same, optionally only High/Very High polygons |

## Sample design and the spatial split (use it unchanged)

`samples/habitat_samples.csv` (sha256 starts `3bfce9bf49fdc9ed`) is committed so every part trains and tests on
identical points. **Do not redraw it.** It is reproducible with `python prepare_data.py --resample`.

- 1,000 "pure" pixels per class: label reliability High or Very High, and all 8 neighbours have the same class.
- The AOI is cut into 90 blocks of 2.5 km (`hc.BLOCK_SIZE`). Samples within 50 m of a block edge are excluded.
- 27 blocks (30%) are **test**, 63 are **train**; train blocks are dealt into 5 folds (0-4) for blocked
  cross-validation. Test samples have `fold = -1`.

| Class | Train | Test |
| --- | --- | --- |
| bog | 743 | 257 |
| heath | 674 | 326 |
| bracken | 715 | 285 |
| seminat_grass | 630 | 370 |
| improved_grass | 718 | 282 |
| broadleaf | 684 | 316 |
| conifer | 880 | 120 |
| water | 862 | 138 |
| built_bare | 638 | 362 |
| **Total** | **6,544** | **2,456** |

| Function | Returns |
| --- | --- |
| `hc.load_samples(split=None)` | DataFrame: `sample_id, x, y, row, col, class_id, class_key, le_habitat, le_reliability, le_source, block_id, split, fold` (`split="train"` or `"test"` filters) |
| `hc.load_blocks()` | GeoDataFrame of blocks: `block_id, split, fold, geometry` |
| `hc.block_id(x, y)` | Block index for BNG coordinates (row-major from the north-west corner) |
| `hc.sample_values(arr, samples)` | Values of a `(rows, cols)` or `(bands, rows, cols)` array at the samples |

Rules for every part: fit and tune only on `split == "train"` (use `fold` for blocked CV and early stopping);
report the headline numbers on `split == "test"` with the confusion matrix, OA, per-class user's and producer's
accuracy and macro F1, as Part 1 does. The sample is class-balanced, so sample OA is not map accuracy.

## Sentinel-2 L2A (Planetary Computer)

| Function | Returns |
| --- | --- |
| `hc.search_s2(start, end, max_cloud=90)` | STAC items over the AOI, one per date and tile, newest processing baseline kept (2023: tile 30UWE only) |
| `hc.load_s2(items, bands=hc.S2_BANDS)` | Lazy dask `xarray.Dataset` on the shared grid: raw uint16 DNs per band plus `SCL`; 20 m bands bilinear-resampled to 10 m |
| `hc.s2_offsets(items)` | Per-date BOA offset (1000 for baseline 04.00 or later). Reflectance = (DN - offset) / 10000 |
| `hc.s2_clear_masks(items)` | `(clear, frac)`: bool DataArray `(time, y, x)` and clear fraction per date. Clear = SCL 4, 5, 6 or 7 and not within 20 m of SCL 3, 8, 9 or 10. Cached `data/s2_clear_2023.npz` |
| `hc.s2_dates(year, min_clear=0.05)` | Series of clear fraction for dates with at least 5% of the AOI clear (50 dates in 2023) |
| `hc.s2_composite(name, months, year, bands, min_clear)` | Median reflectance (float32, 0-1, NaN if never clear) per band plus `n_clear`, as `xarray.Dataset` of `(y, x)` arrays. Cached `data/s2_<year>_<name>.tif` |
| `hc.s2_seasonal_composites()` | `{"winter", "spring", "summer", "autumn"}` -> `s2_composite` (winter = Jan, Feb and Dec of the same year) |
| `hc.s2_sample_timeseries(samples=None)` | Long DataFrame `sample_id, date, clear, B02 ... B12` (float32 reflectance, NaN when not clear) for all 50 dates (450,000 rows, 39% clear). Cached `data/s2_2023_samples_timeseries.parquet` (4 MB) |

Bands: `hc.S2_BANDS` = B02, B03, B04, B05, B06, B07, B08, B8A, B11, B12. Seasons: `hc.SEASONS`. Clear views per
pixel (median): winter 3, spring 5, summer 7, autumn 6; 3.6% of pixels have no clear winter view.

## Sentinel-1 RTC (Planetary Computer)

| Function | Returns |
| --- | --- |
| `hc.search_s1(start, end, orbit_state=None)` | `sentinel-1-rtc` items (gamma0, terrain-flattened); `orbit_state` = `"ascending"` or `"descending"` to restrict |
| `hc.load_s1(items)` | Lazy `xarray.Dataset` with `vv`, `vh` (linear float32) on the shared grid, grouped by solar day |
| `hc.s1_composite(name, months, year, orbit_state=None)` | Median `vv`, `vh` (linear), `vv_db`, `vh_db`, `ratio_db` (VH - VV in dB). Cached `data/s1_<year>_<name>.tif` |

2023 has 27-31 RTC dates per season over the AOI. Speckle is reduced by the temporal median; see the
`sar-speckle-filter-python` tutorial for spatial filters.

## Baseline to beat (Part 1)

Part 1 (`tutorial-code/habitat-mapping-sentinel-2-python/habitat_rules.py`) maps the nine classes with an
8-rule decision list on seasonal Sentinel-2 indices. On the 2,456 test samples:
overall accuracy **0.600** (block-bootstrap 95% CI 0.527-0.664), macro F1 **0.630**, kappa 0.547. Blocked
5-fold CV on the training blocks: 0.684 ± 0.030. Wall-to-wall pixel agreement with Living England: 0.490.

## Licences and attribution

- Living England Habitat Map 2022-23: © Natural England, Open Government Licence v3.0. Contains Ordnance Survey
  data © Crown copyright and database right.
- Sentinel-2 and Sentinel-1: contains modified Copernicus Sentinel data 2023, accessed through the Microsoft
  Planetary Computer.