{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Flood Early-Warning Model\n",
    "\n",
    "flood_warning.py\n",
    "\n",
    "Capstone project: a next-day flood warning for a river town, from today's\n",
    "rain, upstream rain, river level, soil wetness and dam releases.\n",
    "\n",
    "CBSE Class XII, Artificial Intelligence (843), Part C - Capstone Project.\n",
    "SDG 13 (Climate Action) and SDG 11 (Sustainable Cities and Communities).\n",
    "\n",
    "Data Science Methodology, 10 steps. CLASSIFICATION of a RARE event: floods\n",
    "happen on a few days a year, so accuracy is useless and the project is\n",
    "judged on recall, precision and F1. Two models - Logistic Regression and a\n",
    "Random Forest, both told to weight the rare class - are tested on the two\n",
    "latest monsoons, and checked with leave-one-year-out cross validation.\n",
    "\n",
    "Libraries: pandas, numpy, matplotlib, scikit-learn, joblib."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1 - Business Understanding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"=\" * 64)\n",
    "print(\"A FLOOD EARLY-WARNING MODEL\")\n",
    "print(\"=\" * 64)\n",
    "print(\"Problem : the river crosses its danger level with little warning\")\n",
    "print(\"Goal    : warn a day ahead, so people and cattle can move in time\")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2 - Analytic Approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Analytic approach : predictive - classification of a rare event\")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3 - Data Requirements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import LeaveOneGroupOut\n",
    "from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score, accuracy_score\n",
    "import joblib\n",
    "\n",
    "pd.set_option(\"display.width\", 140)\n",
    "pd.set_option(\"display.max_columns\", 20)\n",
    "DANGER = 104.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4 - Data Collection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"river_gauge.csv\", parse_dates=[\"date\"])\n",
    "print(\"Days collected :\", len(data))\n",
    "print(\"Monsoons       :\", data[\"date\"].dt.year.min(), \"to\", data[\"date\"].dt.year.max())\n",
    "print(data.head(6))\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5 - Data Understanding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 5: Data Understanding ---\")\n",
    "print(\"Missing values:\")\n",
    "print(data.isnull().sum().to_string())\n",
    "print()\n",
    "floods = data[\"flood_next_day\"].sum()\n",
    "print(\"Days followed by a flood :\", floods, \"of\", len(data))\n",
    "print(\"Share of flood days      :\", round(data[\"flood_next_day\"].mean() * 100, 1), \"%\")\n",
    "print(\"Accuracy of a model that NEVER warns :\", round((1 - data[\"flood_next_day\"].mean()) * 100, 1), \"%\")\n",
    "print()\n",
    "\n",
    "# Chart 1: one monsoon, with the danger level.\n",
    "one = data[data[\"date\"].dt.year == 2022]\n",
    "plt.figure(figsize=(10, 4))\n",
    "plt.plot(one[\"date\"], one[\"river_level_m\"], color=\"#4e79a7\")\n",
    "plt.axhline(DANGER, color=\"#e15759\", linestyle=\"--\", label=\"Danger level (104 m)\")\n",
    "plt.title(\"River level through the 2022 monsoon\")\n",
    "plt.ylabel(\"River level (m)\")\n",
    "plt.legend()\n",
    "plt.grid(True, linestyle=\"--\", alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart1_season_2022.png\")\n",
    "plt.close()\n",
    "\n",
    "# Chart 2: how rare floods are.\n",
    "plt.figure(figsize=(6, 4))\n",
    "plt.bar([\"No flood next day\", \"Flood next day\"], data[\"flood_next_day\"].value_counts().sort_index(),\n",
    "        color=[\"#9c9c9c\", \"#e15759\"])\n",
    "plt.title(\"An unbalanced dataset\")\n",
    "plt.ylabel(\"Days\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart2_balance.png\")\n",
    "plt.close()\n",
    "\n",
    "# Chart 3: upstream rain against today's river level.\n",
    "calm = data[data[\"flood_next_day\"] == 0]\n",
    "flood = data[data[\"flood_next_day\"] == 1]\n",
    "plt.figure(figsize=(7, 5))\n",
    "plt.scatter(calm[\"upstream_rain_mm\"], calm[\"river_level_m\"], s=10, alpha=0.4, color=\"#9c9c9c\", label=\"No flood next day\")\n",
    "plt.scatter(flood[\"upstream_rain_mm\"], flood[\"river_level_m\"], s=30, color=\"#e15759\", label=\"Flood next day\")\n",
    "plt.axhline(DANGER, color=\"#e15759\", linestyle=\"--\", linewidth=0.8)\n",
    "plt.title(\"Upstream rain and today's river level\")\n",
    "plt.xlabel(\"Upstream rain today (mm)\")\n",
    "plt.ylabel(\"River level today (m)\")\n",
    "plt.legend()\n",
    "plt.grid(True, linestyle=\"--\", alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart3_upstream.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6 - Data Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 6: Data Preparation ---\")\n",
    "before = len(data)\n",
    "data = data.drop_duplicates(subset=\"date\").sort_values(\"date\").reset_index(drop=True)\n",
    "print(\"Duplicate days removed      :\", before - len(data))\n",
    "data[\"rain_mm\"] = pd.to_numeric(data[\"rain_mm\"])\n",
    "blank = data[\"rain_mm\"].isnull().sum()\n",
    "# A blank local rain reading is filled with that day's upstream rain, the\n",
    "# nearest gauge - better than a seasonal average on a stormy day.\n",
    "data[\"rain_mm\"] = data[\"rain_mm\"].fillna(data[\"upstream_rain_mm\"])\n",
    "print(\"Blank rain readings filled  :\", blank, \"(from the upstream gauge)\")\n",
    "\n",
    "# Feature engineering: rain over the last three days, and the gap between\n",
    "# today's level and the danger mark.\n",
    "data[\"year\"] = data[\"date\"].dt.year\n",
    "data[\"rain_3day_mm\"] = data.groupby(\"year\")[\"rain_mm\"].transform(lambda r: r.rolling(3, min_periods=1).sum())\n",
    "data[\"upstream_3day_mm\"] = data.groupby(\"year\")[\"upstream_rain_mm\"].transform(lambda r: r.rolling(3, min_periods=1).sum())\n",
    "data[\"metres_below_danger\"] = DANGER - data[\"river_level_m\"]\n",
    "features = [\"rain_mm\", \"upstream_rain_mm\", \"rain_3day_mm\", \"upstream_3day_mm\", \"river_level_m\",\n",
    "            \"rise_24h_m\", \"metres_below_danger\", \"soil_wetness\", \"dam_release_cumecs\"]\n",
    "print(\"Features used               :\", len(features))\n",
    "print()\n",
    "\n",
    "# Chart 4: three-day upstream rain for flood and calm days.\n",
    "plt.figure(figsize=(6.5, 4.3))\n",
    "plt.boxplot([data[data[\"flood_next_day\"] == 0][\"upstream_3day_mm\"], data[data[\"flood_next_day\"] == 1][\"upstream_3day_mm\"]],\n",
    "            labels=[\"No flood next day\", \"Flood next day\"])\n",
    "plt.title(\"Upstream rain over the last three days\")\n",
    "plt.ylabel(\"mm\")\n",
    "plt.grid(True, axis=\"y\", linestyle=\"--\", alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart4_upstream_3day.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 7 - Modelling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 7: Modelling ---\")\n",
    "# A warning system is used on future monsoons, so the test is the two\n",
    "# latest monsoons, 2023 and 2024 - never mixed into training.\n",
    "train = data[data[\"year\"] <= 2022]\n",
    "test = data[data[\"year\"] >= 2023]\n",
    "X_train, y_train = train[features], train[\"flood_next_day\"]\n",
    "X_test, y_test = test[features], test[\"flood_next_day\"]\n",
    "print(\"Training days (2015-2022) :\", len(train), \" flood days :\", y_train.sum())\n",
    "print(\"Testing days (2023-2024)  :\", len(test), \" flood days :\", y_test.sum())\n",
    "\n",
    "# class_weight=\"balanced\" tells each model that a missed flood day counts\n",
    "# as much as all the calm days together - otherwise it learns to say \"no\".\n",
    "logistic = make_pipeline(StandardScaler(), LogisticRegression(class_weight=\"balanced\", max_iter=2000))\n",
    "forest = RandomForestClassifier(n_estimators=300, min_samples_leaf=3, class_weight=\"balanced\", random_state=0)\n",
    "logistic.fit(X_train, y_train)\n",
    "forest.fit(X_train, y_train)\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 8 - Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 8: Evaluation ---\")\n",
    "results = {}\n",
    "for name, model in [[\"Logistic Regression\", logistic], [\"Random Forest\", forest]]:\n",
    "    pred = model.predict(X_test)\n",
    "    tn, fp, fn, tp = confusion_matrix(y_test, pred).ravel()\n",
    "    results[name] = pred\n",
    "    print(name)\n",
    "    print(\"  \" + name + \" TP, FN, FP, TN :\", tp, fn, fp, tn)\n",
    "    print(\"  \" + name + \" accuracy  :\", round(accuracy_score(y_test, pred), 3))\n",
    "    print(\"  \" + name + \" precision :\", round(precision_score(y_test, pred, zero_division=0), 3))\n",
    "    print(\"  \" + name + \" recall    :\", round(recall_score(y_test, pred), 3))\n",
    "    print(\"  \" + name + \" F1        :\", round(f1_score(y_test, pred), 3))\n",
    "print()\n",
    "\n",
    "# Leave-one-year-out cross validation: each monsoon in turn is held out and\n",
    "# predicted by a model trained on all the other monsoons. Days from the same\n",
    "# season never sit on both sides of the split.\n",
    "groups = data[\"year\"]\n",
    "logo = LeaveOneGroupOut()\n",
    "pooled = {\"Logistic Regression\": np.zeros(len(data), dtype=int), \"Random Forest\": np.zeros(len(data), dtype=int)}\n",
    "year_recall = {\"Logistic Regression\": [], \"Random Forest\": []}\n",
    "years = []\n",
    "year_floods = []\n",
    "for tr, te in logo.split(data[features], data[\"flood_next_day\"], groups):\n",
    "    years.append(int(groups.iloc[te[0]]))\n",
    "    yt = data[\"flood_next_day\"].iloc[te]\n",
    "    year_floods.append(int(yt.sum()))\n",
    "    for name in pooled:\n",
    "        if name == \"Logistic Regression\":\n",
    "            m = make_pipeline(StandardScaler(), LogisticRegression(class_weight=\"balanced\", max_iter=2000))\n",
    "        else:\n",
    "            m = RandomForestClassifier(n_estimators=300, min_samples_leaf=3, class_weight=\"balanced\", random_state=0)\n",
    "        m.fit(data[features].iloc[tr], data[\"flood_next_day\"].iloc[tr])\n",
    "        p = m.predict(data[features].iloc[te])\n",
    "        pooled[name][te] = p\n",
    "        # Recall is undefined for a monsoon with no floods: nothing to catch.\n",
    "        year_recall[name].append(recall_score(yt, p) if yt.sum() > 0 else np.nan)\n",
    "print(\"Leave-one-year-out recall for each held-out monsoon:\")\n",
    "print(\"  Year  Floods  Logistic  Forest\")\n",
    "for i in range(len(years)):\n",
    "    if year_floods[i] == 0:\n",
    "        print(\"  \", years[i], \"    0     no floods to catch\")\n",
    "    else:\n",
    "        print(\"  \", years[i], str(year_floods[i]).rjust(5), \"   \", round(year_recall[\"Logistic Regression\"][i], 2),\n",
    "              \"    \", round(year_recall[\"Random Forest\"][i], 2))\n",
    "# Pooled: every day's held-out prediction, from the model that never saw its monsoon.\n",
    "y_all = data[\"flood_next_day\"]\n",
    "for name in pooled:\n",
    "    print(name, \"- pooled recall    :\", round(recall_score(y_all, pooled[name]), 3))\n",
    "    print(name, \"- pooled precision :\", round(precision_score(y_all, pooled[name]), 3))\n",
    "    print(name, \"- pooled F1        :\", round(f1_score(y_all, pooled[name]), 3))\n",
    "print()\n",
    "\n",
    "# Chart 5: confusion matrices on 2023-2024.\n",
    "plt.figure(figsize=(10, 4.3))\n",
    "for k, name in enumerate([\"Logistic Regression\", \"Random Forest\"]):\n",
    "    cm = confusion_matrix(y_test, results[name], labels=[1, 0])\n",
    "    plt.subplot(1, 2, k + 1)\n",
    "    plt.imshow(cm, cmap=\"Reds\")\n",
    "    labels = [[\"TP\", \"FN\"], [\"FP\", \"TN\"]]\n",
    "    for r in range(2):\n",
    "        for c in range(2):\n",
    "            plt.text(c, r, labels[r][c] + \"\\n\" + str(cm[r][c]), ha=\"center\", va=\"center\", fontsize=12,\n",
    "                     color=\"white\" if cm[r][c] > cm.max() / 2 else \"black\")\n",
    "    plt.xticks([0, 1], [\"Warned\", \"Not warned\"])\n",
    "    plt.yticks([0, 1], [\"Flood\", \"No flood\"])\n",
    "    plt.title(name)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart5_confusion.png\")\n",
    "plt.close()\n",
    "\n",
    "# Chart 6: recall in each held-out year (years with no floods left blank).\n",
    "x = np.arange(len(years))\n",
    "plt.figure(figsize=(9, 4))\n",
    "plt.bar(x - 0.2, np.nan_to_num(year_recall[\"Logistic Regression\"]), width=0.4, color=\"#4e79a7\", label=\"Logistic Regression\")\n",
    "plt.bar(x + 0.2, np.nan_to_num(year_recall[\"Random Forest\"]), width=0.4, color=\"#59a14f\", label=\"Random Forest\")\n",
    "for i in range(len(years)):\n",
    "    if year_floods[i] == 0:\n",
    "        plt.text(x[i], 0.04, \"no\\nfloods\", ha=\"center\", fontsize=8)\n",
    "plt.xticks(x, [str(y) + \"\\n(\" + str(f) + \")\" for y, f in zip(years, year_floods)])\n",
    "plt.ylim(0, 1.05)\n",
    "plt.ylabel(\"Recall (floods warned of)\")\n",
    "plt.xlabel(\"Held-out monsoon (flood days)\")\n",
    "plt.title(\"Leave-one-year-out cross validation\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart6_by_year.png\")\n",
    "plt.close()\n",
    "\n",
    "# The threshold: Logistic Regression's probability, and how many floods\n",
    "# are caught - and false alarms raised - at each level.\n",
    "prob = logistic.predict_proba(X_test)[:, 1]\n",
    "thresholds = np.round(np.arange(0.1, 0.91, 0.1), 1)\n",
    "recalls = []\n",
    "alarms = []\n",
    "for t in thresholds:\n",
    "    warn = prob >= t\n",
    "    recalls.append(recall_score(y_test, warn))\n",
    "    alarms.append(int((warn & (y_test == 0)).sum()))\n",
    "chosen = 0.5\n",
    "print(\"Threshold used                     :\", chosen)\n",
    "print(\"False alarms in two monsoons at 0.5 :\", alarms[list(thresholds).index(chosen)])\n",
    "print()\n",
    "\n",
    "# Chart 7: floods caught and false alarms at each threshold.\n",
    "fig, ax1 = plt.subplots(figsize=(8, 4.3))\n",
    "ax1.plot(thresholds, recalls, marker=\"o\", color=\"#e15759\", label=\"Floods warned of (recall)\")\n",
    "ax1.set_ylabel(\"Recall\")\n",
    "ax1.set_ylim(0, 1.05)\n",
    "ax1.set_xlabel(\"Probability needed to issue a warning\")\n",
    "ax2 = ax1.twinx()\n",
    "ax2.bar(thresholds, alarms, width=0.05, alpha=0.35, color=\"#4e79a7\", label=\"False alarms\")\n",
    "ax2.set_ylabel(\"False alarms in 2023-2024\")\n",
    "ax1.set_title(\"Warnings caught against false alarms\")\n",
    "fig.legend(loc=\"upper right\", bbox_to_anchor=(0.88, 0.88))\n",
    "fig.tight_layout()\n",
    "plt.savefig(\"chart7_threshold.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 9 - Deployment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 9: Deployment ---\")\n",
    "joblib.dump(logistic, \"flood_model.pkl\")\n",
    "model = joblib.load(\"flood_model.pkl\")\n",
    "today = pd.DataFrame({\n",
    "    \"rain_mm\": [12.0, 85.0], \"upstream_rain_mm\": [20.0, 140.0], \"rain_3day_mm\": [30.0, 190.0],\n",
    "    \"upstream_3day_mm\": [45.0, 260.0], \"river_level_m\": [98.4, 102.9], \"rise_24h_m\": [0.1, 1.4],\n",
    "    \"metres_below_danger\": [5.6, 1.1], \"soil_wetness\": [0.45, 0.95], \"dam_release_cumecs\": [0, 900],\n",
    "})\n",
    "p_today = model.predict_proba(today)[:, 1]\n",
    "for i, label in enumerate([\"An ordinary monsoon day\", \"After three days of heavy rain\"]):\n",
    "    action = \"ISSUE FLOOD WARNING for tomorrow\" if p_today[i] >= chosen else \"no warning\"\n",
    "    print(label.ljust(32), \"probability\", round(p_today[i], 3), \"->\", action)\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 10 - Feedback"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 10: Feedback ---\")\n",
    "# The 2025 monsoon arrives. Each evening the model issues (or does not\n",
    "# issue) a warning for the next day; afterwards the warnings are checked.\n",
    "s25 = pd.read_csv(\"monsoon_2025.csv\", parse_dates=[\"date\"])\n",
    "s25[\"year\"] = 2025\n",
    "s25[\"rain_3day_mm\"] = s25[\"rain_mm\"].rolling(3, min_periods=1).sum()\n",
    "s25[\"upstream_3day_mm\"] = s25[\"upstream_rain_mm\"].rolling(3, min_periods=1).sum()\n",
    "s25[\"metres_below_danger\"] = DANGER - s25[\"river_level_m\"]\n",
    "warn25 = model.predict_proba(s25[features])[:, 1] >= chosen\n",
    "flood25 = s25[\"flood_next_day\"] == 1\n",
    "print(\"Flood days in the 2025 monsoon      :\", int(flood25.sum()))\n",
    "print(\"Of these, warned the evening before :\", int((warn25 & flood25).sum()))\n",
    "print(\"False alarms in 2025                :\", int((warn25 & ~flood25).sum()))\n",
    "print(\"Recall in 2025                      :\", round(recall_score(flood25, warn25), 3))\n",
    "\n",
    "# Chart 8: the 2025 monsoon, warnings and floods.\n",
    "plt.figure(figsize=(10, 4))\n",
    "plt.plot(s25[\"date\"], s25[\"river_level_m\"], color=\"#4e79a7\", label=\"River level\")\n",
    "plt.axhline(DANGER, color=\"#e15759\", linestyle=\"--\", label=\"Danger level\")\n",
    "warned_days = s25[\"date\"][warn25]\n",
    "plt.scatter(warned_days, [DANGER + 1.2] * len(warned_days), marker=\"v\", color=\"#f28e2b\", s=40, label=\"Warning issued for next day\")\n",
    "plt.title(\"Feedback: the 2025 monsoon\")\n",
    "plt.ylabel(\"River level (m)\")\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.grid(True, linestyle=\"--\", alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart8_monsoon_2025.png\")\n",
    "plt.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}