"""Prepare the shared labels and sample design for the "AI habitat classification with Python" series. 1. Downloads Living England 2022-23 habitat polygons for the Dark Peak AOI (Natural England, OGL v3, no login). 2. Maps them to the 9-class scheme in habitat_common.CLASSES and rasterises them onto the shared 10 m grid. 3. Draws a stratified random sample of "pure" pixels from high-reliability polygons and assigns each to a 2.5 km spatial block; blocks are split 70/30 into train/test with a fixed seed, and train blocks into 5 folds for blocked cross-validation. 4. Writes samples/habitat_samples.csv and samples/blocks.csv (committed, so every part uses the same split). Run: python prepare_data.py labels + samples (samples are only redrawn with --resample) python prepare_data.py --composites also build the 2023 seasonal Sentinel-2 and Sentinel-1 composites python prepare_data.py --timeseries also extract per-date Sentinel-2 reflectance at every sample (Part 4) """ import argparse import hashlib import numpy as np import pandas as pd from scipy import ndimage import habitat_common as hc N_PER_CLASS = 1000 # samples per class (fewer if a class has fewer eligible pixels) TEST_FRACTION = 0.3 # share of blocks held out for testing N_FOLDS = 5 def draw_samples(labels_all, labels_rel): rng = np.random.default_rng(hc.SEED) rows, cols = hc.shape() rr, cc = np.mgrid[0:rows, 0:cols] x = hc.AOI_BOUNDS[0] + (cc + 0.5) * hc.RES y = hc.AOI_BOUNDS[3] - (rr + 0.5) * hc.RES # keep pixel centres at least BLOCK_BUFFER from any block edge, so test pixels never touch training pixels off_x = (x - hc.AOI_BOUNDS[0]) % hc.BLOCK_SIZE off_y = (hc.AOI_BOUNDS[3] - y) % hc.BLOCK_SIZE inner = ((off_x > hc.BLOCK_BUFFER) & (off_x < hc.BLOCK_SIZE - hc.BLOCK_BUFFER) & (off_y > hc.BLOCK_BUFFER) & (off_y < hc.BLOCK_SIZE - hc.BLOCK_BUFFER)) out = [] print(f"\n{'class':<20}{'eligible px':>12}{'sampled':>9}") for c in hc.CLASSES: k = c["id"] # "pure" pixels: reliable label k, and all 8 neighbours also have label k in the full map pure = ndimage.binary_erosion((labels_rel == k) & (labels_all == k), structure=np.ones((3, 3)), border_value=0) idx = np.flatnonzero(pure & inner) pick = np.sort(rng.choice(idx, size=min(N_PER_CLASS, idx.size), replace=False)) r, cl = np.unravel_index(pick, (rows, cols)) out.append(pd.DataFrame({"x": x[r, cl], "y": y[r, cl], "row": r, "col": cl, "class_id": k, "class_key": c["key"]})) print(f"{c['name']:<20}{idx.size:>12,}{pick.size:>9,}") return pd.concat(out, ignore_index=True) def assign_blocks(samples): rng = np.random.default_rng(hc.SEED) ncols = int(np.ceil((hc.AOI_BOUNDS[2] - hc.AOI_BOUNDS[0]) / hc.BLOCK_SIZE)) nrows = int(np.ceil((hc.AOI_BOUNDS[3] - hc.AOI_BOUNDS[1]) / hc.BLOCK_SIZE)) blocks = np.arange(ncols * nrows) order = rng.permutation(blocks) n_test = int(round(TEST_FRACTION * blocks.size)) test = set(order[:n_test].tolist()) train_blocks = [b for b in order if b not in test] fold = {b: i % N_FOLDS for i, b in enumerate(train_blocks)} bdf = pd.DataFrame({"block_id": blocks, "split": ["test" if b in test else "train" for b in blocks], "fold": [-1 if b in test else fold[b] for b in blocks]}) samples["block_id"] = hc.block_id(samples["x"], samples["y"]) samples = samples.merge(bdf, on="block_id", how="left") return samples, bdf def main(): ap = argparse.ArgumentParser() ap.add_argument("--resample", action="store_true", help="redraw samples/*.csv even if they exist") ap.add_argument("--composites", action="store_true", help="also build the S2 and S1 seasonal composites") ap.add_argument("--timeseries", action="store_true", help="also extract per-date S2 reflectance at every sample") args = ap.parse_args() print("=== Living England 2022-23 labels ===") gdf = hc.load_labels() print(f"{len(gdf):,} polygons clipped to the AOI {hc.AOI_BOUNDS} ({hc.CRS})") area = gdf.assign(km2=gdf.area / 1e6).groupby(["le_habitat"])["km2"].sum().sort_values(ascending=False) rel = gdf.assign(km2=gdf.area / 1e6).groupby("le_reliability")["km2"].sum() print(area.round(2).to_string()) print("\narea by reliability (km2):\n" + rel.round(1).to_string()) labels_all = hc.label_raster() labels_rel = hc.rasterize_labels(gdf, reliable_only=True) print(f"\nshared grid: {hc.shape()[1]} x {hc.shape()[0]} px at {hc.RES} m; " f"{(labels_all > 0).mean():.1%} labelled, {(labels_rel > 0).mean():.1%} labelled with High/Very High reliability") for c in hc.CLASSES: print(f" {c['id']} {c['name']:<20} {(labels_all == c['id']).mean():6.1%} of the AOI") if hc.SAMPLES_CSV.exists() and not args.resample: samples = hc.load_samples() print(f"\nusing the committed sample design: {hc.SAMPLES_CSV.name} ({len(samples):,} samples)") else: print("\n=== stratified sample of pure, high-reliability pixels ===") samples = draw_samples(labels_all, labels_rel) samples, bdf = assign_blocks(samples) le = gdf[["le_habitat", "le_reliability", "le_source", "geometry"]] import geopandas as gpd pts = gpd.GeoDataFrame(samples, geometry=gpd.points_from_xy(samples.x, samples.y), crs=hc.CRS) pts = gpd.sjoin(pts, le, how="left", predicate="within").drop_duplicates("row col".split()) samples = pd.DataFrame(pts.drop(columns=["geometry", "index_right"])) samples.insert(0, "sample_id", np.arange(len(samples))) cols = ["sample_id", "x", "y", "row", "col", "class_id", "class_key", "le_habitat", "le_reliability", "le_source", "block_id", "split", "fold"] hc.SAMPLES_CSV.parent.mkdir(exist_ok=True) samples[cols].to_csv(hc.SAMPLES_CSV, index=False) bdf.to_csv(hc.BLOCKS_CSV, index=False) samples = hc.load_samples() blocks = pd.read_csv(hc.BLOCKS_CSV) print(f"\n=== spatial split: {hc.BLOCK_SIZE / 1000:g} km blocks, seed {hc.SEED} ===") print(f"{len(blocks)} blocks: {(blocks.split == 'train').sum()} train, {(blocks.split == 'test').sum()} test") tab = pd.crosstab(samples["class_key"], samples["split"]).reindex([c["key"] for c in hc.CLASSES]) tab["total"] = tab.sum(axis=1) print(tab.to_string()) print(f"total {len(samples):,} samples; {samples.block_id.nunique()} blocks contain samples") print(f"sha256(habitat_samples.csv) = {hashlib.sha256(hc.SAMPLES_CSV.read_bytes()).hexdigest()[:16]}") if args.composites: print("\n=== Sentinel-2 seasonal composites ===") frac = hc.s2_dates() print(f"{len(frac)} dates in {hc.YEAR} with >= 5% of the AOI clear") hc.s2_seasonal_composites() print("\n=== Sentinel-1 RTC seasonal composites ===") for s, m in hc.SEASONS.items(): hc.s1_composite(s, m) if args.timeseries: print("\n=== Sentinel-2 per-date reflectance at every sample ===") ts = hc.s2_sample_timeseries() print(f"{len(ts):,} rows: {ts.sample_id.nunique():,} samples x {ts.date.nunique()} dates; " f"{ts.clear.mean():.0%} of observations clear") if __name__ == "__main__": main()