{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0acd499d",
   "metadata": {},
   "source": [
    "# Python for Climate Science\n",
    "\n",
    "An introduction for climate scientists. Follow along by running each cell as we go.\n",
    "\n",
    "**What we will cover**\n",
    "\n",
    "1. Python basics: variables, lists, loops, functions, modules\n",
    "2. Reading station data with **pandas** (`Sampledata.csv`)\n",
    "3. Plotting time series with **matplotlib**\n",
    "4. Reading gridded data with **xarray** (`fnl_africa_202302-202303.nc`)\n",
    "5. Spatial plots: maps, regions, multi-panel figures\n",
    "6. Saving your results, the Climate Data Store, and useful links\n",
    "\n",
    "**Files you need in the same folder as this notebook**\n",
    "\n",
    "- `Sampledata.csv`: daily rainfall (mm) for three sites, 2005 to 2024\n",
    "- `fnl_africa_202302-202303.nc`: NCEP FNL analysis cut to Africa, Feb to Mar 2023 (a subset of the 1.3 GB global file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6af70b17",
   "metadata": {},
   "source": [
    "## How to use this notebook\n",
    "\n",
    "A Jupyter notebook is made of **cells**. A cell is either *Markdown* (text, like this one) or *code*.\n",
    "\n",
    "- Click a cell and press **Shift + Enter** to run it and move to the next cell.\n",
    "- Code cells show a number in `[ ]` on the left once they have run. `[*]` means it is still running.\n",
    "- The result of the **last line** of a code cell is displayed automatically. Use `print()` to show anything else.\n",
    "- If things go wrong, use the menu **Kernel → Restart** and run the cells again from the top.\n",
    "\n",
    "Try it now: run the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f599ff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Hello, climate scientists!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94bf83ba",
   "metadata": {},
   "source": [
    "## Part 1: Python basics\n",
    "\n",
    "Python is a high-level programming language known for its simplicity and readability. It was created by Guido van Rossum and first released in 1991.\n",
    "\n",
    "It is free and open source, and it has a large ecosystem of scientific libraries. That is why it has become the most common language in climate science.\n",
    "\n",
    "### Variables\n",
    "\n",
    "A variable is a name that stores a value. You create one with the assignment operator `=`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc97b6ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "months = 10\n",
    "months"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a144d07",
   "metadata": {},
   "source": [
    "Rules for variable names:\n",
    "\n",
    "- They cannot start with a number. They can start with a letter or an underscore.\n",
    "- They are case sensitive: `Temp` and `temp` are different variables.\n",
    "- They can only contain letters, numbers and underscores."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49968685",
   "metadata": {},
   "source": [
    "### Types\n",
    "\n",
    "Every value has a type. The most common ones are:\n",
    "\n",
    "| Type | Example | Used for |\n",
    "|---|---|---|\n",
    "| `int` | `10` | whole numbers |\n",
    "| `float` | `10.12` | decimal numbers |\n",
    "| `str` | `\"Hello World\"` | text |\n",
    "| `bool` | `True`, `False` | yes/no |\n",
    "| `list` | `[1, 2, 4, 5]` | an ordered collection you can change |\n",
    "| `tuple` | `(1, 2, 4, 5)` | an ordered collection you cannot change |\n",
    "\n",
    "The `type()` function tells you the type of a value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f294bc45",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_years = 30                  # int\n",
    "temp_mean = 16.4              # float\n",
    "station = \"Cape Town\"         # str\n",
    "is_raining = False            # bool\n",
    "\n",
    "print(type(n_years), type(temp_mean), type(station), type(is_raining))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc4943bc",
   "metadata": {},
   "source": [
    "### Arithmetic\n",
    "\n",
    "Addition `+`, subtraction `-`, multiplication `*`, division `/`, exponentiation `**`, remainder `%`, integer division `//`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8a3104e",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_kelvin = 298.15\n",
    "temp_celsius = temp_kelvin - 273.15\n",
    "print(\"Temperature in °C:\", temp_celsius)\n",
    "\n",
    "# Converting a precipitation rate from kg m-2 s-1 to mm/day\n",
    "prate = 0.00005\n",
    "mm_per_day = prate * 86400     # seconds in a day\n",
    "print(\"Rainfall in mm/day:\", mm_per_day)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9675bd00",
   "metadata": {},
   "source": [
    "### Strings\n",
    "\n",
    "Strings are text. You can join them with `+`, and you can build text that includes values using an **f-string**: `f\"...\"` with variables inside `{}`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbff98bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "site = \"Site A\"\n",
    "rain = 12.5\n",
    "print(site + \" had rain today\")\n",
    "print(f\"{site} recorded {rain} mm of rain\")\n",
    "print(f\"{rain:.1f} mm\")     # :.1f means 1 decimal place"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9354b9d4",
   "metadata": {},
   "source": [
    "### Lists\n",
    "\n",
    "A list holds many values in order. Python counts from **0**, so the first element is `[0]`.\n",
    "\n",
    "Negative indices count from the end: `[-1]` is the last element. A *slice* `[1:3]` gives elements 1 and 2 (the end is not included)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74dab151",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_day = [18, 20, 24, 16]\n",
    "\n",
    "print(temp_day[0])       # first element\n",
    "print(temp_day[1])       # second element\n",
    "print(temp_day[-1])      # last element\n",
    "print(temp_day[1:3])     # a slice\n",
    "print(len(temp_day))     # how many elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cddf1bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_day.append(22)      # add a value to the end\n",
    "temp_day[1] = 21         # change the second value\n",
    "temp_day"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f503556f",
   "metadata": {},
   "source": [
    "A **tuple** looks similar but uses round brackets and **cannot be changed**. Trying to change it gives an error. Errors are normal, read the last line of the message."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7be41bea",
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = (-33.9, 18.4)     # latitude, longitude of Cape Town\n",
    "print(coords[0])\n",
    "\n",
    "# Uncomment the next line to see the error\n",
    "# coords[0] = -30"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eaaf6f7",
   "metadata": {},
   "source": [
    "### Control flow: `if`, `elif`, `else`\n",
    "\n",
    "Let's imagine a task: a hot day warning system. We are given the temperature for the day and the past mean for that day. If the temperature is 4 degrees or more above the mean, print a warning.\n",
    "\n",
    "Notice the **indentation**: Python uses the spaces at the start of a line to know which code belongs inside the `if`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33f647d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_mean = 16\n",
    "temp_day = 20\n",
    "\n",
    "if temp_day >= temp_mean + 4:\n",
    "    print(\"HOT - Issue a warning\")\n",
    "else:\n",
    "    print(\"No warning necessary\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2e276dd",
   "metadata": {},
   "source": [
    "### Loops: `for`\n",
    "\n",
    "What if we have more than one temperature? A `for` loop runs the same code once for each item in a list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9627e447",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_mean = 16\n",
    "temp_day = [18, 20, 24, 16]\n",
    "\n",
    "for t in temp_day:\n",
    "    if t >= temp_mean + 4:\n",
    "        print(t, \"HOT - Issue a warning\")\n",
    "    else:\n",
    "        print(t, \"No warning necessary\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99ab66ee",
   "metadata": {},
   "source": [
    "### Functions\n",
    "\n",
    "A function is a piece of code written to carry out a specific task. It may take inputs (arguments) and it may `return` a value.\n",
    "\n",
    "Python already has many built-in functions, for example `print()`, `len()`, `max()`. You can also define your own with `def`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfc0f85e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def warn_temp(temp_day, temp_mean):\n",
    "    if temp_day >= temp_mean + 4:\n",
    "        return \"HOT - Issue a warning\"\n",
    "    else:\n",
    "        return \"No warning necessary\"\n",
    "\n",
    "\n",
    "temp_mean = 16\n",
    "temp_day = [18, 20, 24, 16]\n",
    "\n",
    "for t in temp_day:\n",
    "    print(t, warn_temp(t, temp_mean))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae518542",
   "metadata": {},
   "source": [
    "### Modules\n",
    "\n",
    "A module is a file of Python code that someone else has written: functions, classes and variables you can reuse. You load one with `import`.\n",
    "\n",
    "```python\n",
    "import module                      # use as module.function()\n",
    "import module as md                # shorter name, use as md.function()\n",
    "from module import function1       # use as function1()\n",
    "```\n",
    "\n",
    "The modules we use today:\n",
    "\n",
    "| Module | Usual short name | What it does |\n",
    "|---|---|---|\n",
    "| `numpy` | `np` | fast arrays and maths |\n",
    "| `pandas` | `pd` | tables (rows and columns), CSV files |\n",
    "| `xarray` | `xr` | labelled multi-dimensional arrays, NetCDF files |\n",
    "| `matplotlib.pyplot` | `plt` | plotting |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee42faab",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "temps = np.array([18, 20, 24, 16])\n",
    "print(\"mean:\", temps.mean())\n",
    "print(\"max:\", temps.max())\n",
    "print(\"in Kelvin:\", temps + 273.15)      # arithmetic on the whole array at once"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f9037e4",
   "metadata": {},
   "source": [
    "### Exercise: Python basics\n",
    "\n",
    "1. Declare an integer variable and display it.\n",
    "2. Declare a float variable and display it.\n",
    "3. Declare a list variable with 5 temperatures and display the second element.\n",
    "4. Write a function `to_celsius(kelvin)` that converts a temperature from Kelvin to °C and test it with 300.\n",
    "5. Use a `for` loop to print the temperature in °C for every value in the list `[280.5, 295.0, 310.2]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "094b8380",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e90289bb",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>Click to show a solution</b></summary>\n",
    "\n",
    "```python\n",
    "n_stations = 12\n",
    "print(n_stations)\n",
    "\n",
    "rainfall = 4.5\n",
    "print(rainfall)\n",
    "\n",
    "temps = [17, 21, 25, 19, 23]\n",
    "print(temps[1])\n",
    "\n",
    "def to_celsius(kelvin):\n",
    "    return kelvin - 273.15\n",
    "\n",
    "print(to_celsius(300))\n",
    "\n",
    "for k in [280.5, 295.0, 310.2]:\n",
    "    print(f\"{k} K = {to_celsius(k):.1f} °C\")\n",
    "```\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35f857a5",
   "metadata": {},
   "source": [
    "## Part 2: Station data with pandas\n",
    "\n",
    "**pandas** works with tables: rows and columns, like a spreadsheet. The table is called a **DataFrame**, and a single column is a **Series**.\n",
    "\n",
    "We start with `Sampledata.csv`, daily rainfall for three sites from 2005 to 2024. The file uses `;` to separate columns, so we tell pandas that with `sep=\";\"`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0344033",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we import our libraries\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e014151",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"Sampledata.csv\", sep=\";\")\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1b15d32",
   "metadata": {},
   "source": [
    "### Looking at the data\n",
    "\n",
    "A DataFrame has some useful methods and attributes for a first look. Methods need `()`, attributes do not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "194d8ab9",
   "metadata": {},
   "outputs": [],
   "source": [
    "data.head()          # first 5 rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b637e7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "data.tail(3)         # last 3 rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64337599",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data.shape)    # (rows, columns)\n",
    "print(data.columns)  # column names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ad87973",
   "metadata": {},
   "outputs": [],
   "source": [
    "data.describe()      # summary statistics for every numeric column"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b38e0bc2",
   "metadata": {},
   "source": [
    "### Selecting a column\n",
    "\n",
    "Use square brackets with the column name. The result is a Series, and a Series has statistics built in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70a3d303",
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"Site_A_Mediterranean_mm\"].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10c0d762",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Site A max:\", data[\"Site_A_Mediterranean_mm\"].max())\n",
    "print(\"Site B max:\", data[\"Site_B_TropicalMonsoon_mm\"].max())\n",
    "print(\"Site C max:\", data[\"Site_C_TemperateMaritime_mm\"].max())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d503b4ea",
   "metadata": {},
   "source": [
    "### Working with dates\n",
    "\n",
    "The date is split over three columns (`year`, `month`, `day`). Combine them into one proper date column with `pd.to_datetime`. Dates let pandas and matplotlib understand time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ccdaded",
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"Date\"] = pd.to_datetime(data[[\"year\", \"month\", \"day\"]])\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d1f80cf",
   "metadata": {},
   "source": [
    "### Grouping: climatology and annual totals\n",
    "\n",
    "`groupby` splits the table into groups and applies a statistic to each group. Grouping by `month` and taking the mean gives the **mean annual cycle**. Grouping by `year` and taking the sum gives **annual totals**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49542bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "sites = [\"Site_A_Mediterranean_mm\", \"Site_B_TropicalMonsoon_mm\", \"Site_C_TemperateMaritime_mm\"]\n",
    "\n",
    "monthly_mean = data.groupby(\"month\")[sites].mean()\n",
    "monthly_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3162de7",
   "metadata": {},
   "outputs": [],
   "source": [
    "annual = data.groupby(\"year\")[sites].sum()\n",
    "annual.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cb3ed57",
   "metadata": {},
   "source": [
    "### Selecting rows with a condition\n",
    "\n",
    "A comparison on a column gives `True`/`False` for every row. Put that inside `[]` to keep only the `True` rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04cda6d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "wet_days = data[data[\"Site_A_Mediterranean_mm\"] > 20]\n",
    "print(\"Number of days with more than 20 mm at Site A:\", len(wet_days))\n",
    "wet_days.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd8e8473",
   "metadata": {},
   "source": [
    "### Exercise: pandas\n",
    "\n",
    "1. What is the mean daily rainfall at Site B?\n",
    "2. How many days had **no rain at all** (0.0 mm) at Site A?\n",
    "3. Calculate the mean rainfall per **year** for Site C (not the sum).\n",
    "4. Which year had the highest annual total at Site B? Hint: `annual[\"Site_B_TropicalMonsoon_mm\"].idxmax()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed027933",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edbbfaee",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>Click to show a solution</b></summary>\n",
    "\n",
    "```python\n",
    "print(data[\"Site_B_TropicalMonsoon_mm\"].mean())\n",
    "\n",
    "dry = data[data[\"Site_A_Mediterranean_mm\"] == 0.0]\n",
    "print(len(dry))\n",
    "\n",
    "print(data.groupby(\"year\")[\"Site_C_TemperateMaritime_mm\"].mean())\n",
    "\n",
    "print(annual[\"Site_B_TropicalMonsoon_mm\"].idxmax())\n",
    "```\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10c097e2",
   "metadata": {},
   "source": [
    "## Part 3: Time series plots with matplotlib\n",
    "\n",
    "The basic recipe for every plot:\n",
    "\n",
    "```python\n",
    "plt.figure(figsize=(width, height))   # 1. make a figure\n",
    "plt.plot(x, y, label=\"...\")            # 2. draw the data\n",
    "plt.xlabel(\"...\")                      # 3. label the axes\n",
    "plt.ylabel(\"...\")\n",
    "plt.title(\"...\")\n",
    "plt.legend()                           # 4. legend (if more than one line)\n",
    "plt.show()                             # 5. show it\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3db350f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting daily rainfall at one site\n",
    "plt.figure(figsize=(12, 4))\n",
    "\n",
    "plt.plot(data[\"Date\"], data[\"Site_A_Mediterranean_mm\"])\n",
    "\n",
    "plt.xlabel(\"Date\")\n",
    "plt.ylabel(\"Rainfall (mm)\")\n",
    "plt.title(\"Daily rainfall - Site A\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96c7fde1",
   "metadata": {},
   "source": [
    "### Several lines on one plot\n",
    "\n",
    "Call `plt.plot` once per line and give each a `label`. Then `plt.legend()` shows which line is which."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f080e5fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "\n",
    "plt.plot(monthly_mean.index, monthly_mean[\"Site_A_Mediterranean_mm\"], marker=\"o\", label=\"Site A (Mediterranean)\")\n",
    "plt.plot(monthly_mean.index, monthly_mean[\"Site_B_TropicalMonsoon_mm\"], marker=\"o\", label=\"Site B (Tropical monsoon)\")\n",
    "plt.plot(monthly_mean.index, monthly_mean[\"Site_C_TemperateMaritime_mm\"], marker=\"o\", label=\"Site C (Temperate maritime)\")\n",
    "\n",
    "plt.xlabel(\"Month\")\n",
    "plt.ylabel(\"Mean daily rainfall (mm/day)\")\n",
    "plt.title(\"Mean annual cycle of rainfall, 2005–2024\")\n",
    "plt.xticks(range(1, 13))\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39eb3071",
   "metadata": {},
   "source": [
    "### Anomalies\n",
    "\n",
    "An anomaly is the difference from the long-term mean. Subtracting the mean from the annual totals shows which years were wetter or drier than usual. `plt.axhline(0)` draws the zero line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40cfee80",
   "metadata": {},
   "outputs": [],
   "source": [
    "anomaly = annual - annual.mean()\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "\n",
    "plt.plot(anomaly.index, anomaly[\"Site_A_Mediterranean_mm\"], marker=\"o\", label=\"Site A\")\n",
    "plt.plot(anomaly.index, anomaly[\"Site_B_TropicalMonsoon_mm\"], marker=\"o\", label=\"Site B\")\n",
    "plt.plot(anomaly.index, anomaly[\"Site_C_TemperateMaritime_mm\"], marker=\"o\", label=\"Site C\")\n",
    "\n",
    "plt.axhline(0, color=\"black\", linestyle=\"--\")\n",
    "\n",
    "plt.xlabel(\"Year\")\n",
    "plt.ylabel(\"Rainfall anomaly (mm)\")\n",
    "plt.title(\"Annual rainfall anomalies, 2005–2024\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "006c6c14",
   "metadata": {},
   "source": [
    "### Bar charts and saving a figure\n",
    "\n",
    "`plt.bar` draws bars. `plt.savefig(\"name.png\")` saves the figure to a file. Call it **before** `plt.show()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7642e90e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7, 5))\n",
    "\n",
    "plt.bar([\"Site A\", \"Site B\", \"Site C\"], annual.mean(), yerr=annual.std(), capsize=5)\n",
    "\n",
    "plt.ylabel(\"Annual rainfall (mm)\")\n",
    "plt.title(\"Mean annual rainfall ± standard deviation\")\n",
    "\n",
    "plt.savefig(\"mean_annual_rainfall.png\", dpi=150, bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a122aaa0",
   "metadata": {},
   "source": [
    "### Exercise: time series plots\n",
    "\n",
    "1. Plot the annual totals (`annual`) for all three sites on one figure, with labels, a legend and a grid.\n",
    "2. Plot the daily rainfall for Site B for the year **2015 only**. Hint: `data[data[\"year\"] == 2015]`.\n",
    "3. Save the second plot as `site_b_2015.png`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b94ea3b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8125a97f",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>Click to show a solution</b></summary>\n",
    "\n",
    "```python\n",
    "plt.figure(figsize=(10, 5))\n",
    "for site in sites:\n",
    "    plt.plot(annual.index, annual[site], marker=\"o\", label=site)\n",
    "plt.xlabel(\"Year\")\n",
    "plt.ylabel(\"Annual rainfall (mm)\")\n",
    "plt.title(\"Annual rainfall totals\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()\n",
    "\n",
    "data_2015 = data[data[\"year\"] == 2015]\n",
    "\n",
    "plt.figure(figsize=(12, 4))\n",
    "plt.plot(data_2015[\"Date\"], data_2015[\"Site_B_TropicalMonsoon_mm\"])\n",
    "plt.xlabel(\"Date\")\n",
    "plt.ylabel(\"Rainfall (mm)\")\n",
    "plt.title(\"Daily rainfall - Site B, 2015\")\n",
    "plt.savefig(\"site_b_2015.png\", dpi=150, bbox_inches=\"tight\")\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "531dc156",
   "metadata": {},
   "source": [
    "## Part 4: Gridded data with xarray\n",
    "\n",
    "Climate data is usually **gridded**: a value at every latitude, longitude, (level) and time. That is more than two dimensions, so a table does not fit. The standard file format is **NetCDF** (`.nc`), and the standard Python library for it is **xarray**.\n",
    "\n",
    "A NetCDF file is *self-describing*: it stores the data, the coordinates (time, lat, lon, level) and metadata (units, long names) together, following the **CF conventions**.\n",
    "\n",
    "xarray has two main objects:\n",
    "\n",
    "- **Dataset**: the whole file, a collection of variables sharing coordinates\n",
    "- **DataArray**: one variable, with its dimensions, coordinates and attributes\n",
    "\n",
    "Our file is the NCEP **FNL** (Final) global analysis for February and March 2023, cut down to Africa so it downloads quickly. The original global file is 1.3 GB and has the same structure, just more of it.\n",
    "\n",
    "| | |\n",
    "|---|---|\n",
    "| Grid | 1° × 1°, 86 latitudes × 91 longitudes, 40°N to 45°S and 25°W to 65°E |\n",
    "| Time | every 6 hours, 236 time steps |\n",
    "| Levels | 3 pressure levels: 500, 850, 1000 hPa (stored in Pa) |\n",
    "| Variables | `t` temperature (K), `gh` geopotential height (gpm), `u_2`/`v_2` wind (m/s), `prate` precipitation rate (kg m⁻² s⁻¹) |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb6c34e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1856d771",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reading in our netcdf file\n",
    "ds = xr.open_dataset(\"fnl_africa_202302-202303.nc\")\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edeb24af",
   "metadata": {},
   "source": [
    "Click the small icons in the output above to expand the coordinates, data variables and attributes.\n",
    "\n",
    "Notice this cell ran instantly, even though the file holds tens of millions of values. `open_dataset` is *lazy*: it reads the structure and metadata only. Data values are loaded when you actually use them. The full 1.3 GB global file opens just as fast. That is why the habit is to **select** what you need before computing.\n",
    "\n",
    "### Picking a variable\n",
    "\n",
    "Use `ds[\"t\"]` or `ds.t` to get one variable as a DataArray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8db9370d",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = ds[\"t\"]\n",
    "t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3407681d",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(t.dims)                   # the dimensions\n",
    "print(t.shape)                  # their sizes\n",
    "print(t.attrs[\"units\"])         # metadata (attributes)\n",
    "print(t.attrs[\"long_name\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4798e76f",
   "metadata": {},
   "source": [
    "### Selecting by label: `.sel`\n",
    "\n",
    "`.sel` selects using coordinate **values** (a real latitude, a real date). `.isel` selects by **position** (0, 1, 2, ...), like a list.\n",
    "\n",
    "- `method=\"nearest\"` finds the closest grid point when your value is not exactly on the grid.\n",
    "- `slice(start, stop)` selects a range. Unlike Python lists, **both ends are included**.\n",
    "- Latitude in this file runs from **40 down to −45** (north to south), so a latitude slice goes from the larger value to the smaller: `slice(0, -40)`.\n",
    "- Longitude runs from **−25 to 65**. Many files store longitude as 0 to 359 instead, so always check `ds.lon` before selecting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1febb3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Surface (1000 hPa) temperature, first time step\n",
    "t_surface_first = t.sel(plev=100000).isel(time=0)\n",
    "t_surface_first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5a8241e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The grid point nearest Cape Town, at 1000 hPa, for all times\n",
    "t_cape_town = t.sel(plev=100000, lat=-33.9, lon=18.4, method=\"nearest\")\n",
    "t_cape_town"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8579491",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A region: southern Africa, all times, 1000 hPa\n",
    "t_sa = t.sel(plev=100000, lat=slice(0, -40), lon=slice(10, 52))\n",
    "print(t_sa.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76a37a71",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A range of dates\n",
    "t_feb = t.sel(time=slice(\"2023-02-01\", \"2023-02-28\"))\n",
    "print(t_feb.time.values[0], \"to\", t_feb.time.values[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "641dd9f5",
   "metadata": {},
   "source": [
    "### Converting units\n",
    "\n",
    "Arithmetic works on the whole array at once. Temperature is in Kelvin and precipitation rate is in kg m⁻² s⁻¹ (which is mm/s). Let's convert them, and keep the metadata up to date so our plots label themselves correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fbcce4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_c = t - 273.15\n",
    "t_c.attrs[\"units\"] = \"°C\"\n",
    "t_c.attrs[\"long_name\"] = \"Temperature\"\n",
    "\n",
    "pr = ds[\"prate\"] * 86400          # kg m-2 s-1 -> mm/day\n",
    "pr.attrs[\"units\"] = \"mm/day\"\n",
    "pr.attrs[\"long_name\"] = \"Precipitation\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98215754",
   "metadata": {},
   "source": [
    "### Reductions: `.mean()`, `.max()`, `.sum()` along a dimension\n",
    "\n",
    "Pass `dim=` to say which dimension to collapse. Averaging over `time` gives a map. Averaging over `lat` and `lon` gives a time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5c5b9af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A map: mean 1000 hPa temperature over the two months\n",
    "t_c_mean_map = t_c.sel(plev=100000).mean(dim=\"time\")\n",
    "t_c_mean_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c6cda9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A time series: area-average rainfall over southern Africa\n",
    "pr_sa = pr.sel(lat=slice(0, -40), lon=slice(10, 52))\n",
    "pr_sa_series = pr_sa.mean(dim=(\"lat\", \"lon\"))\n",
    "pr_sa_series"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df836a5f",
   "metadata": {},
   "source": [
    "### Time operations: `resample` and `groupby`\n",
    "\n",
    "Our data is 6-hourly. `resample(time=\"1D\")` turns it into daily values. `groupby(\"time.month\")` groups all time steps by calendar month. Both work like pandas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19a31da0",
   "metadata": {},
   "outputs": [],
   "source": [
    "pr_sa_daily = pr_sa_series.resample(time=\"1D\").mean()\n",
    "pr_sa_daily"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7995faf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monthly mean maps for the region: February and March\n",
    "pr_sa_monthly = pr_sa.groupby(\"time.month\").mean(dim=\"time\")\n",
    "pr_sa_monthly"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3aabefa3",
   "metadata": {},
   "source": [
    "### Exercise: xarray\n",
    "\n",
    "1. Select the 500 hPa geopotential height (`ds[\"gh\"]`, `plev=50000`) on 15 February 2023 at 12:00. Hint: `time=\"2023-02-15T12\"`.\n",
    "2. Calculate the **maximum** 1000 hPa temperature over time at the grid point nearest Nairobi (lat −1.3, lon 36.8). Print it in °C.\n",
    "3. Calculate the mean 850 hPa wind speed over the two months for southern Africa. Wind speed is `sqrt(u² + v²)`, use `np.sqrt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d7e5ba0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a4641aa",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>Click to show a solution</b></summary>\n",
    "\n",
    "```python\n",
    "gh_500 = ds[\"gh\"].sel(plev=50000, time=\"2023-02-15T12\")\n",
    "print(gh_500)\n",
    "\n",
    "t_nairobi = t_c.sel(plev=100000, lat=-1.3, lon=36.8, method=\"nearest\")\n",
    "print(float(t_nairobi.max()))\n",
    "\n",
    "u = ds[\"u_2\"].sel(plev=85000, lat=slice(0, -40), lon=slice(10, 52))\n",
    "v = ds[\"v_2\"].sel(plev=85000, lat=slice(0, -40), lon=slice(10, 52))\n",
    "speed = np.sqrt(u**2 + v**2)\n",
    "speed_mean = speed.mean(dim=\"time\")\n",
    "print(speed_mean)\n",
    "```\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7a9f2bf",
   "metadata": {},
   "source": [
    "## Part 5: Spatial plots\n",
    "\n",
    "A DataArray has a built-in `.plot()` that uses matplotlib. With **two** dimensions left (lat, lon) it draws a map. With **one** dimension left (time) it draws a line. The axis labels and colour bar come from the metadata, which is why we kept the `units` attribute up to date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1cdce70",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 7))\n",
    "t_c_mean_map.plot()\n",
    "plt.title(\"Mean 1000 hPa temperature, Feb–Mar 2023\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f4b720e",
   "metadata": {},
   "source": [
    "### Choosing a colour map\n",
    "\n",
    "Use `cmap=` to set the colours. Match the colour map to the data:\n",
    "\n",
    "- **Sequential** (one hue, light to dark) for magnitudes: `\"viridis\"`, `\"Blues\"`, `\"YlGnBu\"`, `\"Oranges\"`\n",
    "- **Diverging** (two hues around a neutral middle) for anomalies and differences: `\"RdBu_r\"`, `\"BrBG\"`, `\"coolwarm\"`. Use `center=0` or `vmin`/`vmax` so the middle is zero.\n",
    "- Avoid `\"jet\"` / rainbow colour maps: they invent features that are not in the data and are hard to read for colour-blind readers.\n",
    "\n",
    "`vmin` and `vmax` fix the colour range, useful when comparing panels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcd410f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "pr_mean_map = pr.mean(dim=\"time\")\n",
    "\n",
    "plt.figure(figsize=(9, 7))\n",
    "pr_mean_map.plot(cmap=\"YlGnBu\", vmin=0, vmax=15)\n",
    "plt.title(\"Mean precipitation, Feb–Mar 2023\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f77b77a",
   "metadata": {},
   "source": [
    "### A region\n",
    "\n",
    "Selecting a region first makes the plot faster and clearer. Here is March minus February rainfall over southern Africa, a difference, so we use a diverging colour map centred on zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c04295e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "pr_diff = pr_sa_monthly.sel(month=3) - pr_sa_monthly.sel(month=2)\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "pr_diff.plot(cmap=\"BrBG\", center=0)\n",
    "plt.title(\"Precipitation difference: March − February 2023 (mm/day)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbe8dd70",
   "metadata": {},
   "source": [
    "### Adding coastlines with cartopy\n",
    "\n",
    "Plain lat/lon axes are fine for a quick look, but a real map needs coastlines. **cartopy** adds map projections and coastlines to matplotlib. It is an optional extra: if it is not installed, this cell prints a message and moves on.\n",
    "\n",
    "Install it with `conda install -c conda-forge cartopy` (or `pip install cartopy`). The first run downloads the coastline data, so it needs internet access."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80e9b815",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import cartopy.crs as ccrs\n",
    "    import cartopy.feature as cfeature\n",
    "\n",
    "    fig = plt.figure(figsize=(9, 7))\n",
    "    ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "\n",
    "    t_c_mean_map.sel(lat=slice(0, -40), lon=slice(10, 52)).plot(\n",
    "        ax=ax, transform=ccrs.PlateCarree(), cmap=\"RdYlBu_r\",\n",
    "        cbar_kwargs={\"label\": \"Temperature (°C)\"},\n",
    "    )\n",
    "    ax.coastlines()\n",
    "    ax.add_feature(cfeature.BORDERS, linewidth=0.5)\n",
    "    ax.gridlines(draw_labels=True)\n",
    "    ax.set_title(\"Mean 1000 hPa temperature, Feb–Mar 2023\")\n",
    "    plt.show()\n",
    "\n",
    "except ImportError:\n",
    "    print(\"cartopy is not installed. Install it with: conda install -c conda-forge cartopy\")\n",
    "except Exception as e:\n",
    "    print(\"cartopy could not draw the map (it may need internet access for coastlines):\", e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "653f649a",
   "metadata": {},
   "source": [
    "### Several panels in one figure\n",
    "\n",
    "`plt.subplots(rows, cols)` makes a grid of axes. Pass `ax=` to `.plot()` to say which panel to draw in. Using the same `vmin`/`vmax` in every panel makes them comparable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a14c882",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
    "\n",
    "pr_sa_monthly.sel(month=2).plot(ax=axes[0], cmap=\"YlGnBu\", vmin=0, vmax=15)\n",
    "axes[0].set_title(\"February 2023\")\n",
    "\n",
    "pr_sa_monthly.sel(month=3).plot(ax=axes[1], cmap=\"YlGnBu\", vmin=0, vmax=15)\n",
    "axes[1].set_title(\"March 2023\")\n",
    "\n",
    "fig.suptitle(\"Mean precipitation (mm/day), southern Africa\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38a9eb5c",
   "metadata": {},
   "source": [
    "Look at the band of heavy rain running from the Mozambique Channel inland over Mozambique and Malawi in the March panel. That is the track of **Tropical Cyclone Freddy**, which made landfall in Mozambique on 24 February and again on 11 March 2023. Two months of reanalysis, four lines of code, and the event is visible."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6061cc52",
   "metadata": {},
   "source": [
    "### Contours and wind vectors\n",
    "\n",
    "`.plot.contour()` draws contour lines, the classic way to show geopotential height. `plt.quiver` draws wind arrows. We draw every 3rd grid point (`[::3]`) so the arrows do not overlap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e51b25dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 500 hPa geopotential height and 850 hPa winds on 15 February 2023, 12:00\n",
    "when = \"2023-02-15T12\"\n",
    "region = dict(lat=slice(0, -40), lon=slice(10, 52))\n",
    "\n",
    "gh500 = ds[\"gh\"].sel(plev=50000, time=when, **region)\n",
    "u850 = ds[\"u_2\"].sel(plev=85000, time=when, **region)\n",
    "v850 = ds[\"v_2\"].sel(plev=85000, time=when, **region)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 7))\n",
    "\n",
    "cs = gh500.plot.contour(ax=ax, levels=15, colors=\"black\", linewidths=0.8)\n",
    "ax.clabel(cs, fontsize=8, fmt=\"%.0f\")\n",
    "\n",
    "ax.quiver(u850.lon[::3], u850.lat[::3], u850[::3, ::3], v850[::3, ::3], color=\"tab:blue\")\n",
    "\n",
    "ax.set_title(f\"500 hPa geopotential height (gpm) and 850 hPa wind, {when}\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fe5d50d",
   "metadata": {},
   "source": [
    "### Time series from gridded data\n",
    "\n",
    "Once you reduce the data to one dimension (time), `.plot()` draws a line. Here is the daily area-averaged rainfall over southern Africa, and the temperature at the grid point nearest Cape Town."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d6d8ede",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 1, figsize=(12, 7), sharex=True)\n",
    "\n",
    "pr_sa_daily.plot(ax=axes[0])\n",
    "axes[0].set_title(\"Daily mean precipitation, southern Africa (0 to 40°S, 10 to 52°E)\")\n",
    "axes[0].set_ylabel(\"mm/day\")\n",
    "axes[0].grid(True)\n",
    "\n",
    "(t_cape_town - 273.15).plot(ax=axes[1])\n",
    "axes[1].set_title(\"6-hourly 1000 hPa temperature near Cape Town\")\n",
    "axes[1].set_ylabel(\"°C\")\n",
    "axes[1].grid(True)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffc59c84",
   "metadata": {},
   "source": [
    "### Exercise: spatial plots\n",
    "\n",
    "1. Plot a map of the mean 500 hPa geopotential height (`ds[\"gh\"]`, `plev=50000`) over the whole domain for February 2023.\n",
    "2. Plot the 1000 hPa temperature **anomaly** on 15 February 2023 12:00 relative to the two-month mean (`t_c_mean_map`), over southern Africa. Use a diverging colour map centred on zero.\n",
    "3. Make a 1 × 2 figure showing mean precipitation for East Africa (lat 15 to −15, lon 25 to 55) in February and March with the same colour scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28aea1ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fd5e1c0",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>Click to show a solution</b></summary>\n",
    "\n",
    "```python\n",
    "gh_feb = ds[\"gh\"].sel(plev=50000, time=slice(\"2023-02-01\", \"2023-02-28\")).mean(dim=\"time\")\n",
    "plt.figure(figsize=(9, 7))\n",
    "gh_feb.plot(cmap=\"viridis\")\n",
    "plt.title(\"Mean 500 hPa geopotential height, February 2023\")\n",
    "plt.show()\n",
    "\n",
    "anom = t_c.sel(plev=100000, time=\"2023-02-15T12\") - t_c_mean_map\n",
    "plt.figure(figsize=(8, 6))\n",
    "anom.sel(lat=slice(0, -40), lon=slice(10, 52)).plot(cmap=\"RdBu_r\", center=0)\n",
    "plt.title(\"1000 hPa temperature anomaly, 2023-02-15 12:00 (°C)\")\n",
    "plt.show()\n",
    "\n",
    "pr_ea = pr.sel(lat=slice(15, -15), lon=slice(25, 55)).groupby(\"time.month\").mean(dim=\"time\")\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
    "pr_ea.sel(month=2).plot(ax=axes[0], cmap=\"YlGnBu\", vmin=0, vmax=15)\n",
    "axes[0].set_title(\"February 2023\")\n",
    "pr_ea.sel(month=3).plot(ax=axes[1], cmap=\"YlGnBu\", vmin=0, vmax=15)\n",
    "axes[1].set_title(\"March 2023\")\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abe02fa0",
   "metadata": {},
   "source": [
    "## Part 6: Saving results\n",
    "\n",
    "- A DataArray or Dataset: `.to_netcdf(\"file.nc\")`\n",
    "- A DataFrame: `.to_csv(\"file.csv\")`\n",
    "- A figure: `plt.savefig(\"file.png\")` before `plt.show()`\n",
    "\n",
    "Give saved variables a name and keep the attributes, so the file you write is as self-describing as the one you read."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8993136c",
   "metadata": {},
   "outputs": [],
   "source": [
    "pr_sa_monthly.name = \"pr\"\n",
    "pr_sa_monthly.to_netcdf(\"pr_southern_africa_monthly_2023.nc\")\n",
    "\n",
    "# Convert the daily series to pandas and save as CSV\n",
    "pr_sa_daily.to_dataframe(name=\"pr_mm_day\").to_csv(\"pr_southern_africa_daily_2023.csv\")\n",
    "\n",
    "# Read the NetCDF back to check\n",
    "xr.open_dataset(\"pr_southern_africa_monthly_2023.nc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "396ea902",
   "metadata": {},
   "source": [
    "## Part 7: The Climate Data Store (CDS)\n",
    "\n",
    "The [Copernicus Climate Data Store](https://cds.climate.copernicus.eu/) holds ERA5 reanalysis, seasonal forecasts, CMIP6 projections and many more datasets. You can download data by hand from the website, or from Python with the `cdsapi` package.\n",
    "\n",
    "**One-time setup**\n",
    "\n",
    "1. Create an account at https://cds.climate.copernicus.eu/ and log in.\n",
    "2. Open your profile page and copy your **Personal Access Token**.\n",
    "3. Create a file called `.cdsapirc` in your home folder containing:\n",
    "\n",
    "```\n",
    "url: https://cds.climate.copernicus.eu/api\n",
    "key: YOUR-PERSONAL-ACCESS-TOKEN\n",
    "```\n",
    "\n",
    "4. Accept the licence of each dataset you want to use, on that dataset's download page.\n",
    "5. Install the client: `pip install cdsapi`\n",
    "\n",
    "**Requesting data**\n",
    "\n",
    "The easiest way to build a request is to use the download form on the dataset page and click *Show API request*. It produces code like the cell below. The cell is not run here because it needs your token and an internet connection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5517057e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example request: ERA5 2 m temperature, daily at 12:00 for February 2023, southern Africa\n",
    "# Remove the triple quotes to run it once you have set up your .cdsapirc file.\n",
    "\n",
    "'''\n",
    "import cdsapi\n",
    "\n",
    "client = cdsapi.Client()\n",
    "\n",
    "client.retrieve(\n",
    "    \"reanalysis-era5-single-levels\",\n",
    "    {\n",
    "        \"product_type\": \"reanalysis\",\n",
    "        \"variable\": \"2m_temperature\",\n",
    "        \"year\": \"2023\",\n",
    "        \"month\": \"02\",\n",
    "        \"day\": [f\"{d:02d}\" for d in range(1, 29)],\n",
    "        \"time\": \"12:00\",\n",
    "        \"area\": [0, 10, -40, 52],      # North, West, South, East\n",
    "        \"format\": \"netcdf\",\n",
    "    },\n",
    "    \"era5_t2m_feb2023.nc\",\n",
    ")\n",
    "\n",
    "era5 = xr.open_dataset(\"era5_t2m_feb2023.nc\")\n",
    "era5\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ea645c3",
   "metadata": {},
   "source": [
    "### Working in JupyterLab on your own computer\n",
    "\n",
    "Today you ran this notebook in **JupyterLab** on your own machine. To set it up again later, or on another computer:\n",
    "\n",
    "1. Install [Miniconda](https://www.anaconda.com/download/success) (or Anaconda).\n",
    "2. Open a terminal (Anaconda Prompt on Windows) and create an environment with everything we used:\n",
    "\n",
    "```\n",
    "conda create -n climate -c conda-forge python=3.12 xarray netcdf4 pandas matplotlib cartopy jupyterlab\n",
    "conda activate climate\n",
    "```\n",
    "\n",
    "3. Go to the folder with your notebook and data, then start JupyterLab. It opens in your browser:\n",
    "\n",
    "```\n",
    "cd path/to/your/folder\n",
    "jupyter lab\n",
    "```\n",
    "\n",
    "4. On the left is the **file browser**. Double-click a notebook to open it. The kernel name is in the top right corner: it should say `Python 3`.\n",
    "5. Save with **Ctrl/Cmd + S**. Export a finished notebook with **File → Save and Export Notebook As → HTML** to share it with someone who does not have Python.\n",
    "6. When you are done: **File → Shut Down**, then close the browser tab.\n",
    "\n",
    "No install possible? [Google Colab](https://colab.research.google.com/) runs notebooks in the browser for free. Upload the notebook and data files, and run `!pip install xarray netcdf4 cartopy` in the first cell.\n",
    "\n",
    "## Useful links\n",
    "\n",
    "**Documentation and tutorials**\n",
    "\n",
    "- xarray: https://docs.xarray.dev/ (start with the *User Guide* and the *Tutorial*)\n",
    "- xarray tutorial (hands-on): https://tutorial.xarray.dev/\n",
    "- pandas: https://pandas.pydata.org/docs/user_guide/\n",
    "- matplotlib gallery, copy from the examples: https://matplotlib.org/stable/gallery/\n",
    "- cartopy: https://scitools.org.uk/cartopy/docs/latest/\n",
    "- Project Pythia, Python for the geosciences: https://foundations.projectpythia.org/\n",
    "- Software Carpentry, programming for scientists: https://software-carpentry.org/lessons/\n",
    "\n",
    "**Data**\n",
    "\n",
    "- Copernicus Climate Data Store: https://cds.climate.copernicus.eu/\n",
    "- cdsapi documentation: https://cds.climate.copernicus.eu/how-to-api\n",
    "- NCEP FNL analysis (today's file): https://rda.ucar.edu/datasets/d083003/\n",
    "- CRU TS station-based gridded data: https://crudata.uea.ac.uk/cru/data/hrg/\n",
    "- Climate Explorer (KNMI): https://climexp.knmi.nl/\n",
    "\n",
    "**Getting help**\n",
    "\n",
    "- Read the error message from the bottom up. The last line says what went wrong.\n",
    "- `help(xr.open_dataset)` or `xr.open_dataset?` in a cell shows the documentation.\n",
    "- Press **Tab** after a dot (`ds.`) to see what is available.\n",
    "- Search the exact error text online. Stack Overflow usually has an answer."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
