{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1 - Features: from 40 years of Landsat to one row per building pixel\n",
    "\n",
    "This notebook turns the two raw parquet files into two feature tables.\n",
    "\n",
    "It does **two** things at once:\n",
    "\n",
    "1. **Reproduces the organisers' 60 features exactly** (Notebook 3), so every\n",
    "   number we report stays comparable with the published Random-Forest\n",
    "   baseline of 0.618 +- 0.004 macro F1 on Madrid.\n",
    "2. **Adds an extended set** on top: proper surface reflectance instead of raw\n",
    "   digital numbers, correct spectral indices, construction-event detection,\n",
    "   per-decade means, trends, and label-free neighbourhood context.\n",
    "\n",
    "Nothing here touches `y_amsterdam`. Amsterdam's *features* are used (for\n",
    "per-city statistics), its *labels* are not.\n",
    "\n",
    "Runs on CPU in a few minutes - no GPU needed. Do this once, locally, and the\n",
    "GPU notebooks just load the parquet files it writes."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "# Make src/ importable whether we are local, on Kaggle, or on Colab.\n",
    "for cand in [Path.cwd() / \"src\", Path.cwd().parent / \"src\",\n",
    "             Path(\"/kaggle/working/src\"), Path(\"/content/hack4dev/src\")]:\n",
    "    if cand.exists():\n",
    "        sys.path.insert(0, str(cand))\n",
    "        break\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import config as C\n",
    "import features as F\n",
    "import retrieval as R\n",
    "\n",
    "pd.set_option(\"display.width\", 200)\n",
    "pd.set_option(\"display.max_columns\", 60)\n",
    "print(\"platform:\", C.PLATFORM)\n",
    "print(\"looking for raw parquet under:\", C.RAW_DIR)"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 Find the raw files\n",
    "\n",
    "Expected: `madrid_train.parquet` (about 3.17 M rows) and\n",
    "`amsterdam_data.parquet` (about 1.06 M rows). On Kaggle they will be inside\n",
    "`/kaggle/input/<dataset-name>/`; the search below handles that."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "def find_parquet(pattern: str) -> Path:\n",
    "    roots = [C.RAW_DIR, C.DATA_DIR, C.PROJECT_ROOT, Path(\"/kaggle/input\")]\n",
    "    for root in roots:\n",
    "        if not root.exists():\n",
    "            continue\n",
    "        hits = sorted(root.rglob(pattern))\n",
    "        if hits:\n",
    "            return hits[0]\n",
    "    raise FileNotFoundError(\n",
    "        f\"{pattern} not found. Put the two parquet files in {C.RAW_DIR} \"\n",
    "        \"(local) or attach them as a Kaggle dataset.\")\n",
    "\n",
    "MADRID_PARQUET = find_parquet(\"madrid_train.parquet\")\n",
    "AMS_PARQUET = find_parquet(\"amsterdam_data.parquet\")\n",
    "print(\"Madrid   :\", MADRID_PARQUET)\n",
    "print(\"Amsterdam:\", AMS_PARQUET)"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 What is actually in the raw file\n",
    "\n",
    "One row per *(pixel, year)*. Up to three summer observations per year are\n",
    "stored side by side as `Blue_1/Blue_2/Blue_3`, each with a `qa_valid_N`\n",
    "cloud flag. Band values are raw Landsat Collection-2 integers:\n",
    "\n",
    "$$\\text{surface reflectance} = \\text{DN} \\times 2.75\\times10^{-5} - 0.2$$"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "peek = pd.read_parquet(MADRID_PARQUET)\n",
    "print(f\"{len(peek):,} rows x {peek.shape[1]} columns\")\n",
    "print(\"\\ncolumns:\")\n",
    "print(list(peek.columns))\n",
    "print(\"\\nyears:\", peek[\"year\"].min(), \"-\", peek[\"year\"].max(),\n",
    "      f\"({peek['year'].nunique()} distinct)\")\n",
    "print(\"unique pixels:\", f\"{peek.groupby(C.PID).ngroups:,}\")\n",
    "print(\"\\nfirst row:\")\n",
    "print(peek.iloc[0].to_string())"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# Reflectance sanity check on one band, and the QA-flag hit rate.\n",
    "dn = peek[\"Red_1\"].dropna()\n",
    "sr = dn * C.SR_SCALE + C.SR_OFFSET\n",
    "print(\"Red_1 digital numbers :\", dn.quantile([0, .01, .5, .99, 1]).round(0).to_dict())\n",
    "print(\"Red_1 as reflectance  :\", sr.quantile([0, .01, .5, .99, 1]).round(3).to_dict())\n",
    "for s in (1, 2, 3):\n",
    "    col = f\"qa_valid_{s}\"\n",
    "    if col in peek.columns:\n",
    "        print(f\"{col}: {peek[col].mean():.1%} valid\")\n",
    "del peek"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 Build both feature sets\n",
    "\n",
    "The heavy lifting lives in `src/features.py`. Per city it:\n",
    "\n",
    "* picks the first QA-valid observation per band per year, drops haze\n",
    "  (`Blue > 15000`);\n",
    "* records `has_early_data` / `has_late_data` on the *observed* years;\n",
    "* fills the annual series (linear inside, nearest at the edges);\n",
    "* pivots to one `(pixels x years)` matrix per band, then computes every\n",
    "  statistic with numpy - a whole city in about a minute."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "madrid = F.build_city(MADRID_PARQUET, city_hint=\"Madrid\")"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "amsterdam = F.build_city(AMS_PARQUET, city_hint=\"Amsterdam\")"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.4 Does this match the organisers' own output?\n",
    "\n",
    "If `preprocessed_data.pkl` from the organisers' Notebook 3 is around, compare\n",
    "the 60 official features against it. Their arrays are standardised with\n",
    "Madrid's `StandardScaler`, so we compare **correlation per feature**, which is\n",
    "invariant to that. Anything below 0.999 means our pipeline has drifted."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "pkl = None\n",
    "for root in [C.DATA_DIR, C.PROJECT_ROOT, Path(\"/kaggle/input\")]:\n",
    "    if root.exists():\n",
    "        hits = sorted(root.rglob(\"preprocessed_data.pkl\"))\n",
    "        if hits:\n",
    "            pkl = hits[0]\n",
    "            break\n",
    "\n",
    "if pkl is None:\n",
    "    print(\"preprocessed_data.pkl not found - skipping the cross-check.\")\n",
    "else:\n",
    "    import pickle\n",
    "    with open(pkl, \"rb\") as fh:\n",
    "        ref = pickle.load(fh)\n",
    "    print(\"reference file:\", pkl)\n",
    "    ours = madrid.set_index([\"px_key\", \"py_key\"])\n",
    "    ref_ids = pd.MultiIndex.from_arrays(\n",
    "        [ref[\"pixel_ids_madrid\"][:, 0], ref[\"pixel_ids_madrid\"][:, 1]],\n",
    "        names=[\"px_key\", \"py_key\"])\n",
    "    common = ref_ids.intersection(ours.index)\n",
    "    print(f\"pixels: ours {len(ours):,}, theirs {len(ref_ids):,}, \"\n",
    "          f\"in common {len(common):,}\")\n",
    "    ref_df = pd.DataFrame(ref[\"X_madrid\"], index=ref_ids,\n",
    "                          columns=ref[\"feature_names\"]).loc[common]\n",
    "    our_df = ours.loc[common, ref[\"feature_names\"]]\n",
    "    corr = {c: float(np.corrcoef(ref_df[c], our_df[c])[0, 1])\n",
    "            for c in ref[\"feature_names\"]\n",
    "            if ref_df[c].std() > 0 and our_df[c].std() > 0}\n",
    "    bad = {k: v for k, v in corr.items() if v < 0.999}\n",
    "    print(f\"features compared: {len(corr)}  |  min correlation: \"\n",
    "          f\"{min(corr.values()):.5f}\")\n",
    "    print(\"dead features in their file (zero variance):\",\n",
    "          [c for c in ref[\"feature_names\"] if ref_df[c].std() == 0])\n",
    "    print(\"features below 0.999:\", bad if bad else \"none - pipeline matches\")\n",
    "    # Labels must agree exactly.\n",
    "    y_ref = pd.Series(ref[\"y_madrid\"], index=ref_ids).loc[common]\n",
    "    print(\"label agreement:\",\n",
    "          f\"{(y_ref.to_numpy() == our_df.index.map(ours['age_class']).to_numpy()).mean():.4%}\")"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.5 Class balance, and the domain gap\n",
    "\n",
    "Two facts drive every later design decision:\n",
    "\n",
    "* the classes are unbalanced, differently in each city;\n",
    "* the same feature has a different *level* in each city, by much more than\n",
    "  the gap between classes - which is why a Madrid model collapses to ~0.34 on\n",
    "  Amsterdam."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "for name, df in [(\"Madrid\", madrid), (\"Amsterdam\", amsterdam)]:\n",
    "    n = len(df)\n",
    "    counts = df[\"age_class\"].value_counts().sort_index()\n",
    "    print(f\"\\n{name}: {n:,} labelled pixels\")\n",
    "    for c in C.CLASS_LABELS[name]:\n",
    "        k = int(counts.get(c, 0))\n",
    "        print(f\"  class {c} ({C.CLASS_LABELS[name][c]:>11}): {k:6,}  {k / n:6.1%}\")"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# How far apart are the two cities, feature by feature?\n",
    "# Expressed in Madrid standard deviations, so it is directly comparable.\n",
    "key = [\"sr_brightness_mean\", \"sr_Red_mean\", \"sr_NIR_mean\", \"sr_SWIR1_mean\",\n",
    "       \"sr_NDVI_mean\", \"sr_NDBI_mean\", \"sr_UI_mean\", \"sr_MNDWI_mean\",\n",
    "       \"sr_BSI_mean\", \"sr_brightness_std\", \"sr_brightness_late_minus_early\",\n",
    "       \"step_brightness_mag\", \"step_brightness_njumps\", \"coverage\"]\n",
    "key = [k for k in key if k in madrid.columns and k in amsterdam.columns]\n",
    "gap = pd.DataFrame({\n",
    "    \"madrid_mean\": [madrid[k].mean() for k in key],\n",
    "    \"amsterdam_mean\": [amsterdam[k].mean() for k in key],\n",
    "}, index=key)\n",
    "gap[\"madrid_sd\"] = [madrid[k].std() for k in key]\n",
    "gap[\"gap_in_madrid_sd\"] = (gap.amsterdam_mean - gap.madrid_mean) / gap.madrid_sd\n",
    "print(gap.round(3).sort_values(\"gap_in_madrid_sd\").to_string())"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A feature whose `gap_in_madrid_sd` is near zero is one a Madrid-trained model\n",
    "can still read in Amsterdam. A feature at +/-2 is telling the model which\n",
    "*city* it is looking at, not how old the building is."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# Within-class spread versus between-city gap, for the top offenders.\n",
    "fig, axes = plt.subplots(2, 4, figsize=(18, 7))\n",
    "for ax, k in zip(axes.ravel(), key[:8]):\n",
    "    for name, df, col in [(\"Madrid\", madrid, \"tab:red\"),\n",
    "                          (\"Amsterdam\", amsterdam, \"tab:blue\")]:\n",
    "        ax.hist(df[k], bins=60, density=True, alpha=0.45, label=name, color=col)\n",
    "    ax.set_title(k, fontsize=9)\n",
    "    ax.tick_params(labelsize=7)\n",
    "axes.ravel()[0].legend(fontsize=8)\n",
    "plt.suptitle(\"The same feature in the two cities\", fontweight=\"bold\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.6 Does the construction-event detector actually work?\n",
    "\n",
    "A genuine test of the extended features. For classes 3 and 4 the building went\n",
    "up *during* the satellite record, so `step_brightness_year` should land near\n",
    "`weighted_mean_year`. For classes 1 and 2 it was already standing in 1984, so\n",
    "the detector should find nothing meaningful. If that pattern does not show up,\n",
    "the detector is noise and should be dropped."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "for name, df in [(\"Madrid\", madrid), (\"Amsterdam\", amsterdam)]:\n",
    "    print(f\"\\n{name}\")\n",
    "    for sig in (\"brightness\", \"SWIR1\", \"NDVI\"):\n",
    "        sub = df[df.age_class.isin([3, 4])]\n",
    "        err = (sub[f\"step_{sig}_year\"] - sub[\"weighted_mean_year\"]).abs()\n",
    "        r = np.corrcoef(sub[f\"step_{sig}_year\"], sub[\"weighted_mean_year\"])[0, 1]\n",
    "        print(f\"  step_{sig}_year vs true year (classes 3+4): r={r:+.3f}  \"\n",
    "              f\"median |error| = {err.median():.1f} years  \"\n",
    "              f\"within 5 y: {(err <= 5).mean():.1%}\")\n",
    "    print(\"  mean |step magnitude| by class:\",\n",
    "          df.groupby(\"age_class\")[\"step_brightness_mag\"].apply(\n",
    "              lambda s: round(float(s.abs().mean()), 4)).to_dict())"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))\n",
    "for ax, (name, df) in zip(axes, [(\"Madrid\", madrid), (\"Amsterdam\", amsterdam)]):\n",
    "    sub = df[df.age_class.isin([3, 4])].sample(min(4000, len(df)), random_state=0)\n",
    "    ax.scatter(sub[\"weighted_mean_year\"], sub[\"step_brightness_year\"], s=3,\n",
    "               alpha=0.25)\n",
    "    lims = [1984, 2024]\n",
    "    ax.plot(lims, lims, \"r--\", lw=1)\n",
    "    ax.set(xlabel=\"true construction year\", ylabel=\"detected brightness step year\",\n",
    "           title=f\"{name}: classes 3 and 4\", xlim=lims, ylim=lims)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.7 Neighbourhood context - free information, no labels\n",
    "\n",
    "Buildings go up a street at a time, so a pixel's neighbours carry real signal\n",
    "about its era. These features average a few key quantities over the 8 nearest\n",
    "pixels (within 10 pixel units = 300 m) and also store the *difference* between\n",
    "the pixel and its surroundings. No labels are involved, so this is available\n",
    "in Amsterdam at n = 0."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "madrid = R.add_neighbourhood_features(madrid)\n",
    "amsterdam = R.add_neighbourhood_features(amsterdam)\n",
    "nbr_cols = [c for c in madrid.columns if c.startswith(\"nbr_\")]\n",
    "print(f\"added {len(nbr_cols)} neighbourhood features\")\n",
    "print(madrid[[\"age_class\"] + nbr_cols[:6]].groupby(\"age_class\").mean().round(4).to_string())"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.8 Save\n",
    "\n",
    "One parquet per city. Every later notebook starts from these."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "C.PROC_DIR.mkdir(parents=True, exist_ok=True)\n",
    "madrid.to_parquet(C.PROC_DIR / \"madrid_features.parquet\", index=False)\n",
    "amsterdam.to_parquet(C.PROC_DIR / \"amsterdam_features.parquet\", index=False)\n",
    "for n, df in [(\"madrid\", madrid), (\"amsterdam\", amsterdam)]:\n",
    "    p = C.PROC_DIR / f\"{n}_features.parquet\"\n",
    "    print(f\"{p}  {p.stat().st_size / 1e6:.1f} MB  {df.shape}\")\n",
    "print(\"\\nfeature groups:\",\n",
    "      f\"{len(C.FEATURE_COLS_60)} official,\",\n",
    "      f\"{sum(1 for c in madrid.columns if c.startswith('sr_'))} reflectance-based,\",\n",
    "      f\"{sum(1 for c in madrid.columns if c.startswith('step_'))} change-event,\",\n",
    "      f\"{len(nbr_cols)} neighbourhood\")"
   ],
   "execution_count": null,
   "outputs": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  },
  "accelerator": "GPU",
  "colab": {
   "provenance": [],
   "gpuType": "L4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}