notebooks/worksheet3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Worksheet 3: Comparing models and observations for the historical period\n",
    "\n",
    "In the last worksheet we started to undertake basic data analysis and visulisation in python using Iris. This worksheet picks-up where we left off and introduces us to more advanced data analysis.  This worksheet ends with an extended exercise to start you thinking about coding more independently."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>By the end of this worksheet you should be able to:</b><br> \n",
    "- Extract specific regions of interest from large datasets<br>\n",
    "- Apply more advanced statistical operations to multi-annual data. <br>\n",
    "- Produce your own data processing workflow.<br>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "### [3.1: Region extraction](#3.1)\n",
    "### [3.2: Mean annual cycle calculation](#3.2)\n",
    "### [3.3: Comparing models and observations](#3.3)\n",
    "### [3.4: Climatological mean and annual cycle for an ensemble](#3.4)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preamble\n",
    "Run the code preamble below to import the necessary libraries for this worksheet.\n",
    "\n",
    "To run the code, click in the box below and press <kbd>Ctrl</kbd> + <kbd>Enter</kbd>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code preamble - these libraries will be used in this worksheet.\n",
    "# This code block needs to be re-run every time you restart this worksheet!\n",
    "%matplotlib inline\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.dates as mdates\n",
    "import calendar\n",
    "import iris\n",
    "import iris.coord_categorisation\n",
    "import iris.quickplot as qplt\n",
    "import cartopy.crs as ccrs\n",
    "from mpl_toolkits.axes_grid1 import AxesGrid\n",
    "from cartopy.mpl.geoaxes import GeoAxes\n",
    "from utils import copy_s3_files, flush_data\n",
    "\n",
    "# Turn off future warnings in Iris for readability\n",
    "import warnings\n",
    "\n",
    "warnings.simplefilter(action=\"ignore\", category=FutureWarning)\n",
    "\n",
    "# Provide the names of the directories where the netCDF model files are stored\n",
    "DATADIR = \"data_v2/\"\n",
    "DOMAIN = \"EAS-22\"\n",
    "\n",
    "# CHIRPS dir\n",
    "CHIRPSDIR = os.path.join(DATADIR, \"CHIRPS/\")\n",
    "\n",
    "# Directory name where processed data is stored\n",
    "HISTDIR = os.path.join(DATADIR, DOMAIN, \"historical\")\n",
    "\n",
    "# Directory name where climtology data is stored\n",
    "CLIMDIR = os.path.join(DATADIR, DOMAIN, \"climatology\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='3.1'></a>\n",
    "## 3.1 Region extraction \n",
    "\n",
    "If you have an area or region you want to focus on, you can extract data for the region of interest. This example works through how to constrain your cube by latitude and longitude.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a) Extract the area** around Bangkok for each of the three datasets, by specifiying latitude and longitude coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Question:</b> Sort these jumbled values into the correct latitude and longitude extent for Bangkok, then enter them into the code block below so that they can be used to extract CHIRPS precipitation data over the city. \n",
    "    <br> 14.5\n",
    "    <br> 99.5\n",
    "    <br> 100.5\n",
    "    <br> 13.5\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Enter the ordered latitude and longitude coordinates for Bangkok here:\n",
    "# lon=( , ) # longitude (East - West extent)\n",
    "# lat=( , ) # latuitude (South - North extent)\n",
    "lon =\n",
    "lat ="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the code block below to extract the specified lat, lon coordinates from the **CHIRPS** dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load CHIRPS data\n",
    "obs_cube = iris.load_cube(CHIRPSDIR + \"chirps-v2.0.monthly.global.1981_2018.nc\")\n",
    "\n",
    "# Note that the CHIRPS data are on a regular grid (unlike the  model data) so real latitudes and longitudes are\n",
    "# used to define the region around KL (more on this in section 2.6)\n",
    "obs_cube_BK = obs_cube.intersection(latitude=lat, longitude=lon)\n",
    "\n",
    "# save the constrained cube to directory\n",
    "outfile = os.path.join(CHIRPSDIR, \"chirps.mon.1981_2018.BK.nc\")\n",
    "iris.save(obs_cube_BK, outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b> Recall that <b>CORDEX uses a rotated pole grid</b> for these simulations, so coordinates in the rotated pole frame of reference are required to correctly extract Bangkok from the model data.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__b) Extract the area around Bangkok for HadGEM2-ES and MPI-ESM__ driven simulations. To do this we must specify the **rotated latitude and longitude coordinates**. \n",
    "\n",
    "To convert the true Bangkok coordinates onto the rotated pole frame of reference, we need know the latitude and longitude of the pole in our rotated pole coordinate system.  We can find this from the model cube:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the HadGEM2-ES model cube data as need grid information from it\n",
    "infile = os.path.join(HISTDIR, \"hadgem2-es.mon.1986_2005.GERICS-REMO2015.pr.mmday-1.nc\")\n",
    "hadgem2_cube = iris.load_cube(infile)\n",
    "\n",
    "pole_lat = hadgem2_cube.coord_system().grid_north_pole_latitude\n",
    "pole_lon = hadgem2_cube.coord_system().grid_north_pole_longitude\n",
    "print(\"Pole Latitude: \" + pole_lat)\n",
    "print(\"Pole Longitude: \" + pole_lon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iris includes a function to help us convert between coordinates in a rotated and non-rotated coordinate frame.  We can use this to convert our non-rotated Bangkok coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon = np.array([99.5, 100.5])\n",
    "lat = np.array([13.5, 14.5])\n",
    "rotated_lons, rotated_lats = iris.analysis.cartography.rotate_pole(\n",
    "    lon, lat, pole_lon, pole_lat\n",
    ")\n",
    "print(\n",
    "    \"Rotated Longitudes: \"\n",
    "    + round(rotated_lons[0], 2)\n",
    "    + \" , \"\n",
    "    + round(rotated_lons[1], 2)\n",
    ")\n",
    "print(\n",
    "    \"Rotated Latitudes: \"\n",
    "    + round(rotated_lats[0], 2)\n",
    "    + \" , \"\n",
    "    + round(rotated_lats[1], 2)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the code below to constrain the HadGEM2-ES and MPI-ESM data to Bangkok in rotated pole coordinates: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Constrain the cube area over Bangkok (BK).\n",
    "# CORDEX EAS-22 uses a rotated grid, so the co-ordinates required are different to real world coordinates.\n",
    "\n",
    "for gcm in [\"hadgem2-es\", \"mpi-esm-lr\"]:\n",
    "    # Directory name where data are read from and stored to\n",
    "    infile = os.path.join(HISTDIR, gcm + \".mon.1986_2005.GERICS-REMO2015.pr.mmday-1.nc\")\n",
    "\n",
    "    # Load the baseline precipitation data using the BK_constraint - the command below\n",
    "    # loads the data into a cube constrained by the area chosen\n",
    "    data = iris.load_cube(infile)\n",
    "\n",
    "    # All grid cells whose longitudes and latitudes lie within the limits shown will be selected.\n",
    "    data_BK = data.intersection(grid_longitude=rotated_lons, grid_latitude=rotated_lats)\n",
    "\n",
    "    # save the constrained cube\n",
    "    outfile = os.path.join(\n",
    "        CLIMDIR, gcm + \".mon.1986_2005.GERICS-REMO2015.pr.mmday-1.BK.nc\"\n",
    "    )\n",
    "    iris.save(data_BK, outfile)\n",
    "    print(\"Saved: \" + outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Question:</b> How many grid cells are in the newly constrained Bangkok file?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**\n",
    "*Write your answer here...*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='3.2'></a>\n",
    "## 3.2 Mean annual cycle calculation\n",
    "\n",
    "Average over a given area (e.g Bangkok) to calculate multi-annual monthly mean values. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**c) We now calculate monthly mean fields** for 1986-2005 for each of the twelve months for the Bangkok area."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b> To accurately find averages over latitudinal areas, we must account for changes in the \n",
    "    distance between lines of latitude.  These changes make the areas of grid boxes smaller as we move towards the\n",
    "    model pole.  Iris provides a tool to help us do this: <code>iris.analysis.cartography.area_weights()</code>. Find out more in the Iris documentation <a href='https://scitools.org.uk/iris/docs/latest/userguide/cube_statistics.html?highlight=weighting#area-averaging'>here</a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for gcm in [\"hadgem2-es\", \"mpi-esm-lr\"]:\n",
    "    # Set up the path to the data\n",
    "    infile = os.path.join(\n",
    "        CLIMDIR, gcm + \".mon.1986_2005.GERICS-REMO2015.pr.mmday-1.BK.nc\"\n",
    "    )\n",
    "\n",
    "    # Load the data extracted around Bangkok created in previous step\n",
    "    data = iris.load_cube(infile)\n",
    "\n",
    "    # Add monthly coord categorisation to the time dimension coordinate\n",
    "    iris.coord_categorisation.add_month_number(data, \"time\", name=\"month_number\")\n",
    "\n",
    "    # Calculate monthly mean values\n",
    "    monthly_mean = data.aggregated_by([\"month_number\"], iris.analysis.MEAN)\n",
    "\n",
    "    # due to a limitation in Iris remove the 2D 'True' lat and lon coords\n",
    "    monthly_mean.remove_coord(\"longitude\")\n",
    "    monthly_mean.remove_coord(\"latitude\")\n",
    "\n",
    "    # Find latitude weights\n",
    "    monthly_mean.coord(\"grid_longitude\").guess_bounds()\n",
    "    monthly_mean.coord(\"grid_latitude\").guess_bounds()\n",
    "    grid_areas = iris.analysis.cartography.area_weights(monthly_mean)\n",
    "    # Calculate area averaged monthly mean rainfall\n",
    "    monthly_mean = monthly_mean.collapsed(\n",
    "        [\"grid_longitude\", \"grid_latitude\"], iris.analysis.MEAN, weights=grid_areas\n",
    "    )\n",
    "\n",
    "    # Save the area averaged monthly mean data\n",
    "    outfile = os.path.join(CLIMDIR, gcm + \".mon.mean.1986_2005.pr.mmday-1.BK.nc\")\n",
    "    iris.save(monthly_mean, outfile)\n",
    "    print(\"Saved: {}\".format(outfile))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**d) Plot the BK area averaged monthly mean precipitation** amount in mm/day for the HadGEM2-ES and MPI-ESM driven  runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for gcm in [\"hadgem2-es\", \"mpi-esm-lr\"]:\n",
    "    # Load the model cube\n",
    "    inpath = os.path.join(CLIMDIR, gcm + \".mon.mean.1986_2005.pr.mmday-1.BK.nc\")\n",
    "    print(inpath)\n",
    "    cube = iris.load_cube(inpath)\n",
    "    # Quick line plot for each cube\n",
    "    qplt.plot(cube.coord(\"month_number\"), cube)\n",
    "    plt.title(\"BK area averaged \" + gcm + \"\\n monthly average of daily rainfall\")\n",
    "    ax = plt.gca()\n",
    "    ax.xaxis.set_label_text(\"Month Number\")\n",
    "    ax.set_xlim(0.5, 12.5)\n",
    "    ax.set_ylim(0, 16)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Question:</b> Considering the plots of HadGEM2-ES and MPI-ESM-LR data, note down their approximate values in mm day-1 for the following features: \n",
    "\n",
    "- Which month experiences the most rainfall? (also note the amount)\n",
    "- If a second rainfall peak exists in another season, note the month and amount\n",
    "- Which month experiences the least rainfall? (also note the amount)\n",
    "\n",
    "What comparisons can be made about these two datasets?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Answer:</b><br>\n",
    "    \n",
    "<b>HadGEM2-ES</b>\n",
    "- Most rainfall: *month, amount*\n",
    "- Secondary peak of rainfall:  *month, amount*\n",
    "- Least rainfall:  *month, amount*\n",
    "\n",
    "<b>MPI-ESM-LR</b>\n",
    "- Most rainfall: *month, amount*\n",
    "- Secondary peak of rainfall:  *month, amount*\n",
    "- Least rainfall:  *month, amount*\n",
    "\n",
    "<b>How comparable are these datasets? How are they different?</b><br>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**e) Find the multi-annual monthly means (1986-2005) over the BK area for CHIRPS observations** by following the same methodology as above in step d).  As before we've given you the input and output files names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the BK extracted data created in previous step\n",
    "infile = os.path.join(CHIRPSDIR, \"chirps.mon.1986_2005.BK.nc\")\n",
    "\n",
    "\n",
    "# Add monthly coord categorisation to the time dim coordinate\n",
    "\n",
    "\n",
    "# Now calculate monthly means\n",
    "\n",
    "\n",
    "# and change units to mm/day from mm/month\n",
    "\n",
    "\n",
    "# create the area averaged monthly mean of daily rainfall\n",
    "\n",
    "# Save output\n",
    "outfile = os.path.join(CLIMDIR, \"chirps.mon.mean.1986_2005.mmday-1.nc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Question:</b> Plot the CHIRPS observations, compare with the HadGEM2-ES and MPI-ESM-LR driven REMO2015 runs. What are the differences between the observations and models; in which months are the differences greatest?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Complete your plotting script here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Answer:</b><br>\n",
    "*Write your findings here...*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='3.3'></a>\n",
    "## 3.3 Comparing models and observations\n",
    "\n",
    "In Section 2.4 we saw how to plot model output on individual maps, but in order to do a rigorous comparison, **spatial model and observation fields must be on the same grid.**\n",
    "\n",
    "**You should regrid to the coarsest resolution grid.** For the data used here, the observations have the coarsest resolution so we will regrid the model data onto the observation grid.\n",
    "\n",
    "Recall that the REMO2015 model data are on a grid known as a **Rotated Pole** coordinate system. The idea is that the \"real\" north pole in the Arctic is shifted such that the equator relative to our rotated pole then runs through the centre of the regional model domain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**f) Regrid the multiannual OND mean model** fields onto the observations grid.\n",
    "\n",
    "Here we use the `regrid` method to regrid the target cube. Here we will use linear interpolation. First, load in the data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load HadGEM2-ES\n",
    "infile = os.path.join(CLIMDIR, \"hadgem2-es.OND.mean.1986_2005.pr.mmday-1.nc\")\n",
    "hadgem2_model_cube = iris.load_cube(infile)\n",
    "# load MPI-ESM-LR\n",
    "infile = os.path.join(CLIMDIR, \"mpi-esm-lr.OND.mean.1986_2005.pr.mmday-1.nc\")\n",
    "mpi_model_cube = iris.load_cube(infile)\n",
    "# load\n",
    "infile = os.path.join(CLIMDIR, \"chirps.OND.mean.1986_2005.pr.mmday-1.nc\")\n",
    "obs_cube = iris.load_cube(infile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can regrid any data, both the model and observation datasets must have a **specified coordinate system**. \n",
    "\n",
    "In this case, CHIRPS's coordinate system is missing from the original netCDF file, but the observations are on a regular longitude-latitude grid so the correct coordinate system is [WGS84](https://epsg.io/32600/). \n",
    "\n",
    "In the code below, we define the WGS84 coordinate system and then apply it to the x- and y-axes (i.e. longitudes and latitudes) of the observations.  \n",
    "\n",
    "The coordinate system used by the model (the rotated grid) is already defined.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Define WGS84 projection for obs data\n",
    "wgs84 = iris.coord_systems.GeogCS(\n",
    "    semi_major_axis=6378137.0, inverse_flattening=298.257223563\n",
    ")\n",
    "\n",
    "# Apply WGS84 to obs cube\n",
    "obs_cube.coord(axis=\"x\").coord_system = wgs84\n",
    "obs_cube.coord(axis=\"y\").coord_system = wgs84\n",
    "\n",
    "# Print out and compare the two coordinate systems\n",
    "print(obs_cube.coord_system())\n",
    "print(hadgem2_model_cube.coord_system())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next few lines of code **regrid the model data from rotated pole to the regular grid used by the observations**. \n",
    "\n",
    "From the plots created in section 2.4, we know that the area covered by CHIRPS data is larger than the EAS-22 model domain. Hence, the **extrapolation mode** is set to `mask` so that any grid cells on the CHIRPS grid which do not overlap with model grid cells are masked `off`; otherwise, the model data would be interpolated and would produce misleading results.\n",
    "\n",
    "Note also that as CHIRPS is a near global data set, we take a subset of the data over our region of interest.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>Note:</b> When working with data arrays <b>masks</b> can be extremely useful. Masks are an array of boolean values matching the shape of the data arrya, where a particular condition is met. These boolean arrays can then be used to hide or show parts of the dataset, given the value of the mask. You can see more examples of masks <a href='https://jakevdp.github.io/PythonDataScienceHandbook/02.06-boolean-arrays-and-masks.html'>here</a> \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lon = (60.0, 150.0)\n",
    "lat = (10.0, 50.0)\n",
    "\n",
    "obs_cube_sub = obs_cube.intersection(latitude=lat, longitude=lon)\n",
    "\n",
    "# Regrid the climate model data onto CHIRPS grid\n",
    "hadgem2_regrid = hadgem2_model_cube.regrid(\n",
    "    obs_cube_sub, iris.analysis.Nearest(extrapolation_mode=\"mask\")\n",
    ")\n",
    "mpi_regrid = mpi_model_cube.regrid(\n",
    "    obs_cube_sub, iris.analysis.Nearest(extrapolation_mode=\"mask\")\n",
    ")\n",
    "\n",
    "# Save HadGEM2 output\n",
    "outfile = os.path.join(CLIMDIR, \"hadgem2-es.OND.mean.1986_2005.pr.mmday-1.rg.nc\")\n",
    "iris.save(hadgem2_regrid, outfile)\n",
    "print(\"Saved: \" + outfile)\n",
    "\n",
    "# Save MPI-ESM output\n",
    "outfile = os.path.join(CLIMDIR, \"mpi-esm-lr.OND.mean.1986_2005.pr.mmday-1.rg.nc\")\n",
    "iris.save(mpi_regrid, outfile)\n",
    "print(\"Saved: \" + outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Question:</b> Now that the cubes are all on the same grid, what differences do you see? <br> Complete the code template below to help answer this question.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hint:** You need to (i) load the NetCDF files, (ii) plot the CHIRPS and model data again (as in section 2.4) to compare them visually once again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directory name where data is read from\n",
    "indir = os.path.join(DATADIR, \"EAS-22\", \"climatology\")\n",
    "\n",
    "# load HadGEM2-ES downscaled model data\n",
    "\n",
    "# load MPI-ESM-LR downscaled model data\n",
    "\n",
    "# load CHIRPS data and extract region\n",
    "\n",
    "# Do some plotting!\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**\n",
    "\n",
    "_Write your findings here..._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**g) Find the difference between the model and the observation OND multiannual mean** fields and plot maps to view the differences\n",
    "\n",
    "To find differences,  now that the model and observations are on the same grid we can simply subtract the model data from the observations.  \n",
    "\n",
    "There is a subtract function within Iris but it cannot be used here.  The model cubes contain extra coordinates which are not present in the obs cube; Iris requires all coordinates within the cubes to match exactly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sure units are the same\n",
    "obs_cube.units = hadgem2_regrid.units\n",
    "\n",
    "# Make recieving cube\n",
    "hadgem2_obs_diff = obs_cube.copy()\n",
    "mpi_obs_diff = obs_cube.copy()\n",
    "\n",
    "# Replace data with the differences\n",
    "hadgem2_obs_diff.data = hadgem2_regrid.data - obs_cube.data\n",
    "mpi_obs_diff.data = mpi_regrid.data - obs_cube.data\n",
    "\n",
    "# Plotting\n",
    "plt.figure(figsize=(12, 10))\n",
    "plt.subplot(\n",
    "    1, 2, 1\n",
    ")  # Create a new subplot for the first differences, 2 rows, 1 column, 1st plot\n",
    "\n",
    "# Only plot the region containing data. We use the intersection method to do this.\n",
    "qplt.pcolormesh(\n",
    "    hadgem2_obs_diff[0].intersection(longitude=(90, 135), latitude=(-20, 32)),\n",
    "    vmax=10,\n",
    "    vmin=-10,\n",
    "    cmap=plt.get_cmap(\"RdBu\"),\n",
    ")  # Note this is where cube slicing is needed as you can only plot 2-coordinate\n",
    "# dimensions with qplt.contourf, so here we have selected time[0] as there is only\n",
    "# one timestep (the baseline 1986-2005 mean)\n",
    "\n",
    "plt.title(\"HadGEM2-ES - obs\")  # plots a title for the plot\n",
    "ax = plt.gca()  # gca function that returns the current axes\n",
    "ax.coastlines()  # adds coastlines defined by the axes of the plot\n",
    "\n",
    "plt.subplot(\n",
    "    1, 2, 2\n",
    ")  # Create a new subplot for the model data 2 row, 2 columns, 2nd plot\n",
    "qplt.pcolormesh(\n",
    "    mpi_obs_diff[0].intersection(longitude=(90, 135), latitude=(-20, 32)),\n",
    "    vmax=10,\n",
    "    vmin=-10,\n",
    "    cmap=plt.get_cmap(\"RdBu\"),\n",
    ")\n",
    "\n",
    "plt.title(\"MPI-ESM-LR - obs\")  # plots a title for the plot\n",
    "ax = plt.gca()  # gca function that returns the current axes\n",
    "ax.coastlines()  # adds coastlines defined by the axes of the plot\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Questions:</b> We have now created difference plots; to more easily validate the model(s) directly against observations. \n",
    "    \n",
    "- Why has the data over the oceans disappeared?\n",
    "- Which of the two model datasets (HadGEM2-ES and MPI-ESM-LR) has the greatest wet bias? In which region? \n",
    "- Are there any areas in HadGEM2-ES or MPI-ESM-LR with a strong dry bias? Roughly what is the extent of this bias (in mm/day)? \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Answers</b>\n",
    "- The data over the oceans has disappeared because [...] \n",
    "\n",
    "- The greatest wet bias occus in model [...], in region [...] \n",
    "\n",
    "- A strong dry bias exists in model [...] in region [...] (roughly [...] mm/day drier than CHIRPS observations).  \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='3.4'></a>\n",
    "## 3.4 Climatological mean and annual cycle for an ensemble\n",
    "\n",
    "So far data from two GCMs downscaled with the REMO2015 RCM have been analysed. In this section we will add data from an extra RCM (REGCM4) and an extra GCM giving 6 runs in total. We will also make use of the CRU observational data set. \n",
    "\n",
    "The CRU data are a monthly global land-only dataset (1901-present) at 0.5 degree resolution. Nine variables are available, including mean, min, max temperature and precipitation. For further details please see: https://crudata.uea.ac.uk/cru/data/hrg/ \n",
    "\n",
    "This gives us an ensemble of 6 GCM-RCM pairs:\n",
    "\n",
    "* RCM REMO2015 downscaling:\n",
    " * HadGEM2-ES\n",
    " * MPI-ESM-LR\n",
    " * NorESM1-M\n",
    "* RCM REGCM4-7 downscaling:\n",
    " * HadGEM2-ES\n",
    " * MPI-ESM-MR (note medium resolution not low resolution)\n",
    " * NorESM1-M\n",
    "\n",
    "And observational datasets:\n",
    "\n",
    "* CHIRPS\n",
    "* CRU\n",
    "\n",
    "Taking an ensemble approach allows us to account for a range of uncertainty in the model projections.\n",
    "\n",
    "Write a series of scripts to do the following:\n",
    "\n",
    "**h) Calculate the OND seasonal mean and annual cycle** (for the BK area) for 1.5m temperature and precipitation for CRU and CHIRPS observations\n",
    "\n",
    "**i) Calculate OND seasonal-mean and monthly-mean anomalies** for the BK area for the 4 additional models\n",
    "\n",
    "**j) Plot a series of figures** that shows 1) the monthly cycles of temperature and rainfall comparing the 6 models and the observations; and 2) the monthly differences between the models and observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we first need to download CRU and netcdf data\n",
    "copy_s3_files(\"s3://ias-pyprecis/data/CRU/*.nc\", \"data/CRU/\")\n",
    "copy_s3_files(\"s3://ias-pyprecis/data/netcdf/*.nc\", \"data/netcdf/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Here are some useful variables  you might like to use in your scripts\n",
    "\"\"\"\n",
    "# Some helpful data locations\n",
    "DATADIR = \"data_v2\"\n",
    "DOMAIN = \"EAS-22\"\n",
    "CHIRPSDIR = os.path.join(DATADIR, \"CHIRPS\")\n",
    "CRUDIR = os.path.join(DATADIR, \"CRU\")\n",
    "CLIMDIR = os.path.join(DATADIR, DOMAIN, \"climatology\")\n",
    "MODELDIR = os.path.join(DATADIR, \"cordex\", DOMAIN)\n",
    "\n",
    "# Some helpful model variables\n",
    "GCMS = [\"hadgem2-es\", \"mpi-esm-lr\"]\n",
    "RCMS = [\"GERICS-REMO2015\", \"ICTP-RegCM4-7\"]\n",
    "\n",
    "# Bangkok domains...\n",
    "# ... in true lat-lon coordiates\n",
    "lon = (99.5, 100.5)\n",
    "lat = (13.5, 14.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "h) Calculate the OND seasonal-mean and monthly-mean 1.5m temperature and precipitation \n",
    "for the BK area, for CRU and CHIRPS observations\n",
    "\"\"\"\n",
    "# Load CHIRPS data\n",
    "\n",
    "# Load CRU data\n",
    "\n",
    "# Extract BK area\n",
    "\n",
    "# Add OND season catagorisation\n",
    "\n",
    "# Add monthly catagorisation\n",
    "\n",
    "# Extract season\n",
    "\n",
    "# Aggregate cubes\n",
    "\n",
    "# Find BK area average\n",
    "\n",
    "# Check and add cube metadata\n",
    "\n",
    "# Save cubes to CLIMDIR\n",
    "# Remember to use the same naming convention we used earlier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "i) Calculate OND seasonal-mean and monthly-mean anomalies for the BK area \n",
    "for the 4 additional models runs \n",
    "\"\"\"\n",
    "# Load ensemble members\n",
    "# Remember you need to do this for both precipitation AND temperature\n",
    "\n",
    "# Regrid ensemble members onto observational grid\n",
    "# Remember you need to check your model and obs cubes have the appropriate coordinate systems defined\n",
    "\n",
    "# Extract the BK area. Remember you are now working in true lat-lon coordinates!\n",
    "\n",
    "# Find OND and monthly means\n",
    "\n",
    "# Calculate model anomalies\n",
    "# Remember temp anomaly   = model - CRU data\n",
    "#          precip anomaly = model - CHIRPS data\n",
    "\n",
    "# Check cube metadata consistency and save"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "    <b>Question:</b> What difference would it make if we first extracted the BK area and <em>then</em> regrid the models? <br> \n",
    "Which order is best for preserving data integrity?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Answer:</b><br> \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create four figures:\n",
    "    \n",
    "    i) the monthly cycle of temperature (model and observations) \n",
    "    ii) the monthly cycle of rainfall (model and observations)\n",
    "    iii) the monthly temperature anomaly for each model\n",
    "    iv) the monthly precipitation anomaly for each model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Plot 1: The monthly cycle of temperature (model and observations)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Plot 2: The monthly cycle of precipitation (model and observations)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Plot 3: The monthly temperature anomaly for each model\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Plot 4: The monthly precipitation anomaly for each model\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "    <b>Question:</b> How could you summarise the ensemble variability amongst model members in a plot?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Answer:</b><br> \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "    <b>Question:</b> How does the monthly temperature and precipitation anomaly compare to the OND average?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Answer:</b><br> \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Question:</b> What are the advantages and disadvantages of plotting spatial maps of temperature and rainfall variability over Bangkok?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Answer:</b><br> \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "<center><b>This completes worksheet 3.</b></center><br>\n",
    "    You have used Python and Iris to validate the model's climatology (30 years of baseline data: 1986-2005) against observations, in order to analyse the model's capability of representing real-world trends. <br>\n",
    "To do so, you have:\n",
    "    \n",
    "- calculated and plotted multiannual seasonal and monthly mean data,\n",
    "- regridded rotated pole model data to a regular latitude-longitude grid,\n",
    "- created difference plots to compare model data with observations, and \n",
    "- written your own code to do all of these post-processing and analysis tasks with an multi-model ensemble!<br>\n",
    "\n",
    "In worksheet 4, you will consider future model projections and their changes with respect to the baseline period (calculated in this worksheet).\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p><img src=\"img/MO_MASTER_black_mono_for_light_backg_RBG.png\" alt=\"Met Office Logo\" style=\"float: center; height: 100px;\"/></p>\n",
    "<center>© Crown Copyright, Met Office. All rights reserved.</center>\n",
    "<center>This file is part of PyPrecis and is released under the BSD 3-Clause license.</center>\n",
    "<center>See LICENSE in the root of the repository for full licensing details.</center>"
   ]
  }
 ],
 "metadata": {
  "instance_type": "ml.t3.medium",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  },
  "vscode": {
   "interpreter": {
    "hash": "30fc2c9368d5cc76c4fad6f63328738b6cf0a0091c34cff8da6b063602b9d1b0"
   }
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}