{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fa5c2941",
   "metadata": {},
   "source": [
    "# Pandas — Lab (solutions)\n",
    "\n",
    "**AI for Hydraulics · Part 1 · Week 2**\n",
    "\n",
    "In this lab you will use pandas on **real data**: ten years (2015–2024) of daily discharge measured by the U.S. Geological Survey (USGS) at eight gauging stations, from the humid east of the United States to the semi-arid south-west and the Mediterranean climate of California.\n",
    "\n",
    "**How to work**\n",
    "- Run the cells from top to bottom (`Shift + Enter`).\n",
    "- Replace every `...` with your code.\n",
    "- Each exercise ends with a ✅ **Check** giving the values you should obtain.\n",
    "- The two data files are read directly from the course repository: you need an internet connection.\n",
    "\n",
    "**Units.** USGS publishes values in US units. You will convert them:\n",
    "\n",
    "| Quantity | USGS unit | SI conversion |\n",
    "|---|---|---|\n",
    "| Discharge | ft³/s (cfs) | 1 ft³/s = 0.0283168 m³/s |\n",
    "| Drainage area | mi² | 1 mi² = 2.589988 km² |\n",
    "| Altitude | ft | 1 ft = 0.3048 m |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ee34c81",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "print(pd.__version__)\n",
    "\n",
    "DATA_URL = \"https://raw.githubusercontent.com/rafikoulebsir-usthb/ai-hydraulic-course/main/data/\"\n",
    "SITES_CSV = DATA_URL + \"usgs_sites.csv\"\n",
    "DAILY_CSV = DATA_URL + \"usgs_daily_discharge.csv\"\n",
    "\n",
    "CFS_TO_M3S = 0.0283168\n",
    "MI2_TO_KM2 = 2.589988\n",
    "FT_TO_M = 0.3048"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1200ce02",
   "metadata": {},
   "source": [
    "## Exercise 1 — The station table\n",
    "\n",
    "**a)** Read `SITES_CSV` into a DataFrame `sites`. Make sure the column `site_no` is read as **text**. Display the table and its `info()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7f8e97e",
   "metadata": {},
   "outputs": [],
   "source": [
    "sites = pd.read_csv(SITES_CSV, dtype={\"site_no\": str})\n",
    "sites.info()\n",
    "sites"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ae03f3e",
   "metadata": {},
   "source": [
    "**b)** What happens if you forget the `dtype` option? Read the file again without it into `sites_bad` and print the first three values of `site_no`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f1e0698",
   "metadata": {},
   "outputs": [],
   "source": [
    "sites_bad = pd.read_csv(SITES_CSV)\n",
    "print(sites_bad[\"site_no\"].head(3).tolist())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e318e19",
   "metadata": {},
   "source": [
    "**c)** Add two columns: `area_km2` (drainage area in km²) and `altitude_m` (gauge altitude in m). Then set `site_no` as the index of `sites`.\n",
    "\n",
    "**d)** Using `sites`, find:\n",
    "- the name of the station with the **largest** and the **smallest** drainage area;\n",
    "- the stations located **above 500 m** (names and altitudes);\n",
    "- the number of stations per state.\n",
    "\n",
    "✅ **Check:** `sites` has 8 rows; without `dtype` the codes become `[1646500, 7010000, 8158000]` (leading zeros lost). Largest basin: Mississippi River at St. Louis (≈ 1 805 222 km²); smallest: Arroyo Seco near Pasadena (≈ 41.4 km²). 3 stations are above 500 m, the highest being the Rio Grande at Embudo (≈ 1764.5 m). Arizona and California have 2 stations each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1175a67b",
   "metadata": {},
   "outputs": [],
   "source": [
    "sites[\"area_km2\"] = sites[\"drainage_area_mi2\"] * MI2_TO_KM2\n",
    "sites[\"altitude_m\"] = sites[\"altitude_ft\"] * FT_TO_M\n",
    "sites = sites.set_index(\"site_no\")\n",
    "\n",
    "largest = sites.loc[sites[\"area_km2\"].idxmax(), \"station_name\"]\n",
    "smallest = sites.loc[sites[\"area_km2\"].idxmin(), \"station_name\"]\n",
    "above_500 = sites.loc[sites[\"altitude_m\"] > 500, [\"station_name\", \"altitude_m\"]].sort_values(\"altitude_m\", ascending=False)\n",
    "per_state = sites[\"state\"].value_counts()\n",
    "\n",
    "print(\"Largest :\", largest, f\"({sites['area_km2'].max():,.0f} km²)\")\n",
    "print(\"Smallest:\", smallest, f\"({sites['area_km2'].min():.1f} km²)\")\n",
    "print(above_500.round(1))\n",
    "print(per_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea4882c9",
   "metadata": {},
   "source": [
    "## Exercise 2 — The daily discharge record\n",
    "\n",
    "**a)** Read `DAILY_CSV` into `daily`, with `site_no` as text and `date` as a date. Print its shape, its dtypes and the first and last dates of the record."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92f2282f",
   "metadata": {},
   "outputs": [],
   "source": [
    "daily = pd.read_csv(DAILY_CSV, dtype={\"site_no\": str}, parse_dates=[\"date\"])\n",
    "print(daily.shape)\n",
    "print(daily.dtypes)\n",
    "print(daily[\"date\"].min(), \"→\", daily[\"date\"].max())\n",
    "daily.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb525621",
   "metadata": {},
   "source": [
    "**b)** Add a column `Q` with the discharge in m³/s, and a column `year`.\n",
    "\n",
    "**c)** Ten years from 2015-01-01 to 2024-12-31 contain 3653 days. Count the number of records of each station. Which station has a missing day?\n",
    "\n",
    "**d)** Find **which date** is missing. Hint: build the full calendar with `pd.date_range(\"2015-01-01\", \"2024-12-31\")` and use its `.difference(...)` method with the dates of that station.\n",
    "\n",
    "✅ **Check:** `daily` has shape `(29223, 4)`. Station `08158000` has 3652 records; the missing date is 2017-01-26."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83fc0a64",
   "metadata": {},
   "outputs": [],
   "source": [
    "daily[\"Q\"] = daily[\"discharge_cfs\"] * CFS_TO_M3S\n",
    "daily[\"year\"] = daily[\"date\"].dt.year\n",
    "\n",
    "counts = daily[\"site_no\"].value_counts().sort_index()\n",
    "print(counts)\n",
    "\n",
    "short_site = counts.idxmin()\n",
    "calendar = pd.date_range(\"2015-01-01\", \"2024-12-31\")\n",
    "missing = calendar.difference(daily.loc[daily[\"site_no\"] == short_site, \"date\"])\n",
    "print(short_site, missing)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a09f82e",
   "metadata": {},
   "source": [
    "## Exercise 3 — Data quality flags\n",
    "\n",
    "Every USGS value can carry **qualifiers** in the column `qualifier`: `ESTIMATED` (the value was estimated, e.g. because the gauge was damaged or the channel was frozen), `REVISED`, `EQUIP` (equipment malfunction)… Several qualifiers are separated by `;`. An empty cell means a normal measured value.\n",
    "\n",
    "**a)** How many values of `qualifier` are missing (`NaN`)? Replace them by the empty string `\"\"`, then display the `value_counts()` of the column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dba2f230",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_no_qualifier = daily[\"qualifier\"].isna().sum()\n",
    "print(n_no_qualifier)\n",
    "daily[\"qualifier\"] = daily[\"qualifier\"].fillna(\"\")\n",
    "daily[\"qualifier\"].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbb668f1",
   "metadata": {},
   "source": [
    "**b)** Create a boolean column `estimated`, `True` when the qualifier **contains** `ESTIMATED` (alone or combined with another qualifier). How many estimated values are there in total? Compute the **percentage** of estimated values for each station.\n",
    "\n",
    "**c)** We decide not to trust estimated values. Create a column `Q_clean`, equal to `Q` but with `NaN` where `estimated` is `True` (use `.loc`, not chained indexing). For the Potomac (`01646500`), compare the mean of `Q` and the mean of `Q_clean`, and the number of values used by each (`count()`).\n",
    "\n",
    "✅ **Check:** 28 564 values without qualifier, 378 estimated values. The highest share is at station `11098000` (3.18 %). Potomac: mean of `Q` ≈ 321.3 m³/s on 3653 values, mean of `Q_clean` ≈ 321.6 m³/s on only 3576 values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2991a69e",
   "metadata": {},
   "outputs": [],
   "source": [
    "daily[\"estimated\"] = daily[\"qualifier\"].str.contains(\"ESTIMATED\")\n",
    "n_estimated = daily[\"estimated\"].sum()\n",
    "pct_estimated = (daily.groupby(\"site_no\")[\"estimated\"].mean() * 100).round(2)\n",
    "print(n_estimated)\n",
    "print(pct_estimated)\n",
    "\n",
    "daily[\"Q_clean\"] = daily[\"Q\"]\n",
    "daily.loc[daily[\"estimated\"], \"Q_clean\"] = np.nan\n",
    "\n",
    "potomac = daily[daily[\"site_no\"] == \"01646500\"]\n",
    "print(potomac[[\"Q\", \"Q_clean\"]].agg([\"mean\", \"count\"]).round(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e64349e",
   "metadata": {},
   "source": [
    "## Exercise 4 — Selecting floods and dry days\n",
    "\n",
    "**a)** Select the record of the **American River at Fair Oaks** (`11446500`) into `american`. Display the 5 largest daily discharges with their dates.\n",
    "\n",
    "**b)** In how many days of 2023 did the American River exceed 500 m³/s?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7adb9707",
   "metadata": {},
   "outputs": [],
   "source": [
    "american = daily[daily[\"site_no\"] == \"11446500\"]\n",
    "top5 = american.nlargest(5, \"Q\")[[\"date\", \"Q\"]]\n",
    "print(top5.round(1))\n",
    "\n",
    "n_above_500 = ((american[\"year\"] == 2023) & (american[\"Q\"] > 500)).sum()\n",
    "print(n_above_500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cb01dd8",
   "metadata": {},
   "source": [
    "**c)** The **Arroyo Seco** (`11098000`) drains a small mountain catchment above Pasadena (California). How many days have **zero flow**? How are they distributed by year?\n",
    "\n",
    "**d)** The Mississippi reached its highest discharge of the period on 2019-06-08. Display the discharge of **all** stations on that day, with the station names (hint: filter by date, then `set_index(\"site_no\")` and add the `station_name` column from `sites`).\n",
    "\n",
    "✅ **Check:** largest American River flow ≈ 2333.3 m³/s on 2017-02-10. 18 days above 500 m³/s in 2023. 114 zero-flow days at Arroyo Seco, in years 2015 (3), 2017 (2), 2021 (52), 2022 (44) and 2023 (13). On 2019-06-08 the Mississippi carried ≈ 26 362.94 m³/s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b64f01f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "arroyo = daily[daily[\"site_no\"] == \"11098000\"]\n",
    "n_zero = (arroyo[\"Q\"] == 0).sum()\n",
    "zero_by_year = arroyo.loc[arroyo[\"Q\"] == 0, \"year\"].value_counts().sort_index()\n",
    "print(n_zero)\n",
    "print(zero_by_year)\n",
    "\n",
    "flood_day = daily[daily[\"date\"] == \"2019-06-08\"].set_index(\"site_no\")[[\"Q\"]]\n",
    "flood_day[\"station_name\"] = sites[\"station_name\"]\n",
    "flood_day.round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b91c782",
   "metadata": {},
   "source": [
    "## Exercise 5 — Statistics per station with `groupby`\n",
    "\n",
    "**a)** Using `groupby` and named aggregation, build a table `stats` (one row per station) with the columns `mean`, `median`, `min`, `max` and `std` of `Q`.\n",
    "\n",
    "**b)** Add the **coefficient of variation** `cv = std / mean` and the ratio `max / mean`. Sort the stations by `cv`. Which river is the most **irregular**, and which the most **regular**? Can you explain the most regular one? (Hint: look up *Glen Canyon Dam*.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3bc4dcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats = daily.groupby(\"site_no\").agg(mean=(\"Q\", \"mean\"),\n",
    "                                     median=(\"Q\", \"median\"),\n",
    "                                     min=(\"Q\", \"min\"),\n",
    "                                     max=(\"Q\", \"max\"),\n",
    "                                     std=(\"Q\", \"std\"))\n",
    "stats[\"cv\"] = stats[\"std\"] / stats[\"mean\"]\n",
    "stats[\"max_over_mean\"] = stats[\"max\"] / stats[\"mean\"]\n",
    "stats.sort_values(\"cv\").round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80defbf8",
   "metadata": {},
   "source": [
    "**c)** Compute the **annual mean discharge** of each station as a station × year table (`groupby` + `unstack`). For each station, find the wettest year (`idxmax(axis=1)`).\n",
    "\n",
    "**d)** Compute the **mean monthly discharge** (all years together) as a station × month table: this is the *flow regime* of each river. For each station, find the month with the highest mean flow. Relate the results to the climate: snowmelt (Rio Grande), Mediterranean winter rains (California).\n",
    "\n",
    "✅ **Check:** most irregular: Verde River (cv ≈ 3.50, max/mean ≈ 102.8); most regular: Colorado River at Lees Ferry (cv ≈ 0.23). Wettest years: Potomac 2018, Mississippi 2019, Colorado TX 2018, Rio Grande 2019, Colorado AZ 2016, Verde 2023, Arroyo Seco 2023, American 2017. Months of maximum flow: Potomac, Mississippi and Rio Grande 5 (May), Colorado TX 6, Colorado AZ 7, Verde and Arroyo Seco 3, American 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "807cc423",
   "metadata": {},
   "outputs": [],
   "source": [
    "annual_mean = daily.groupby([\"site_no\", \"year\"])[\"Q\"].mean().unstack()\n",
    "wettest_year = annual_mean.idxmax(axis=1)\n",
    "print(annual_mean.round(1))\n",
    "print(wettest_year)\n",
    "\n",
    "daily[\"month\"] = daily[\"date\"].dt.month\n",
    "regime = daily.groupby([\"site_no\", \"month\"])[\"Q\"].mean().unstack()\n",
    "peak_month = regime.idxmax(axis=1)\n",
    "print(regime.round(1))\n",
    "print(peak_month)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1da3318",
   "metadata": {},
   "source": [
    "## Exercise 6 — Comparing catchments: `merge`\n",
    "\n",
    "To compare rivers of very different sizes, hydrologists divide the discharge by the drainage area.\n",
    "\n",
    "**a)** Merge `stats` with the columns `station_name`, `state` and `area_km2` of `sites` into a table `catch` (both have the station code as index: use `left_index=True, right_index=True`, or `join`). Check that you still have 8 rows.\n",
    "\n",
    "**b)** Compute:\n",
    "- the **specific discharge** `q_spec` = mean discharge / area, in L/s/km²;\n",
    "- the **runoff depth** `runoff_mm` = mean discharge × (365.25 × 86400 s) / area, in mm/year;\n",
    "- the **specific peak** `peak_spec` = maximum discharge / area, in m³/s/km².\n",
    "\n",
    "**c)** Sort by `runoff_mm`. Which catchment produces the most runoff per km²? The least? Which has the highest specific peak?\n",
    "\n",
    "💬 **Discussion:** small catchments usually have the most violent floods per km², yet the tiny Arroyo Seco (41 km²) is only second. Why? (Hint: these are **daily means**. How long does a flood last on a 41 km² mountain catchment?)\n",
    "\n",
    "**d)** Save `catch` (rounded to 3 decimals) to a file `catchment_summary.csv`, then read it back to check it.\n",
    "\n",
    "✅ **Check:** most runoff: American River (q ≈ 21.99 L/s/km², ≈ 693.8 mm/year); least: Colorado River at Austin (≈ 8.0 mm/year). Highest specific peak: American River (≈ 0.477 m³/s/km²), then Arroyo Seco (≈ 0.249)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eafff864",
   "metadata": {},
   "outputs": [],
   "source": [
    "catch = stats.merge(sites[[\"station_name\", \"state\", \"area_km2\"]], left_index=True, right_index=True)\n",
    "print(len(catch))\n",
    "\n",
    "catch[\"q_spec\"] = catch[\"mean\"] / catch[\"area_km2\"] * 1000\n",
    "catch[\"runoff_mm\"] = catch[\"mean\"] * 365.25 * 86400 / (catch[\"area_km2\"] * 1e6) * 1000\n",
    "catch[\"peak_spec\"] = catch[\"max\"] / catch[\"area_km2\"]\n",
    "\n",
    "catch.sort_values(\"runoff_mm\", ascending=False)[[\"station_name\", \"area_km2\", \"q_spec\", \"runoff_mm\", \"peak_spec\"]].round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fa46a2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "catch.round(3).to_csv(\"catchment_summary.csv\")\n",
    "pd.read_csv(\"catchment_summary.csv\", dtype={\"site_no\": str}).head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46377ab7",
   "metadata": {},
   "source": [
    "## Bonus — Do neighbouring rivers flood together?\n",
    "\n",
    "**a)** Use `pivot` (or `pivot_table`) to build a **wide** table `wide` with one row per date and one column per station (values: `Q`).\n",
    "\n",
    "**b)** Compute the correlation matrix between stations with `wide.corr()`. Which pair of **different** stations has the highest correlation? (Hint: use `.stack()` and keep only the pairs where the first code is smaller than the second.)\n",
    "\n",
    "💬 **Discussion:** the Mississippi and the Rio Grande are 1500 km apart. Do they really flood together, or do they only share the same **seasonal** cycle (look at the regimes of Exercise 5d)? Which of the top pairs is more likely to share the same storms?\n",
    "\n",
    "✅ **Check:** `wide` has shape `(3653, 8)`. Most correlated pair: Mississippi – Rio Grande (r ≈ 0.446), then Verde – Arroyo Seco (r ≈ 0.427)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d3c10ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "wide = daily.pivot(index=\"date\", columns=\"site_no\", values=\"Q\")\n",
    "print(wide.shape)\n",
    "\n",
    "corr = wide.corr()\n",
    "pairs = corr.stack()\n",
    "pairs = pairs[pairs.index.get_level_values(0) < pairs.index.get_level_values(1)]\n",
    "print(pairs.sort_values(ascending=False).head(3).round(3))\n",
    "corr.round(2)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}