notebooks/CSSP_20CRDS_Tutorials/tutorial_3_basic_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# <b>Tutorial 3: Basic data analysis</b>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Objectives:\n",
    "\n",
    "In this session we will learn: \n",
    "1. to calculate and visualise annual and monthly means\n",
    "2. to calculate and visualise seasonal means\n",
    "3. to calculate mean differences (anomalies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "1. [Calculate annual and monthly mean](#annual)\n",
    "2. [Calculate seasonal means](#season)\n",
    "3. [Calculating differences (anomalies)](#percent)\n",
    "4. [Exercises](#exercise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>Prerequisites</b> <br> \n",
    "- Basic programming skills in python<br>\n",
    "- Familiarity with python libraries Iris, Numpy and Matplotlib<br>\n",
    "- Basic understanding of climate data<br>\n",
    "- Tutorials 1 and 2\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Calculating annual and monthly mean<a id='annual'></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 Import libraries <a id='explore'></a>\n",
    "Import the necessary libraries. Current datasets are in zarr format, we need zarr and xarray libraries to access the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import zarr\n",
    "import iris\n",
    "import os\n",
    "from catnip.preparation import extract_rot_cube, add_bounds\n",
    "from scripts.xarray_iris_coord_system import XarrayIrisCoordSystem as xics\n",
    "\n",
    "xi = xics()\n",
    "xr.set_options(\n",
    "    display_style=\"text\"\n",
    ")  # Work around for AML bug that won't display HTML output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Set up authentication for the Azure blob store\n",
    "\n",
    "The data for this course is held online in an Azure Blob Storage Service. To access this we use a SAS (shared access signature).  You should have been given the credentials for this service before the course, but if not please ask your instructor. We use the getpass module here to avoid putting the token into the public domain. Run the cell below and in the box enter your SAS and press return. This will store the password in the variable SAS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import getpass\n",
    "\n",
    "# SAS WITHOUT leading '?'\n",
    "SAS = getpass.getpass()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use the Zarr library to connect to this storage. This is a little like opening a file on a local file system but works without downloading the data. This makes use of the Azure Blob Storage service. The zarr.ABStore method returns a zarr.storage.ABSStore object which we can now use to access the Zarr data in the same way we would use a local file. If you have a Zarr file on a local file system you could skip this step and instead just use the path to the Zarr data below when opening the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "store = zarr.ABSStore(\n",
    "    container=\"metoffice-20cr-ds\",\n",
    "    prefix=\"monthly/\",\n",
    "    account_name=\"metdatasa\",\n",
    "    blob_service_kwargs={\"sas_token\": SAS},\n",
    ")\n",
    "type(store)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Read monthly data\n",
    "A Dataset consists of coordinates and data variables. Let's use the xarray's **open_zarr()** method to read all our zarr data into a dataset object and display it's metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use the open_zarr() method to read in the whole dataset metadata\n",
    "dataset = xr.open_zarr(store)\n",
    "# print out the metadata\n",
    "dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert dataset into iris cubelist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create an empty list to hold the iris cubes\n",
    "cubelist = iris.cube.CubeList([])\n",
    "\n",
    "# use the DataSet.apply() to convert the dataset to Iris Cublelist\n",
    "dataset.apply(lambda da: cubelist.append(xi.to_iris(da)))\n",
    "\n",
    "# print out the cubelist.\n",
    "cubelist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Calculating annual mean, maximum and minimum over an area\n",
    "\n",
    "Here we calculate annual mean, maximum and minimum air_temperature over the Shanghai region from 1981 to 2010. \n",
    "\n",
    "We will first need to extract the required variables, extract the Shanghai region and constrain by time period. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract air_temperature\n",
    "air_temp = cubelist.extract(\"air_temperature\")\n",
    "\n",
    "# extracting maximum air temperature\n",
    "cons = iris.Constraint(\n",
    "    cube_func=lambda c: (\"ukmo__process_flags\" in c.attributes)\n",
    "    and (c.attributes[\"ukmo__process_flags\"][0].split(\" \")[0] == \"Maximum\")\n",
    ")\n",
    "air_temp_max = air_temp.extract_strict(cons)\n",
    "\n",
    "# extracting mainimum air temperature\n",
    "cons = iris.Constraint(\n",
    "    cube_func=lambda c: (\"ukmo__process_flags\" in c.attributes)\n",
    "    and (c.attributes[\"ukmo__process_flags\"][0].split(\" \")[0] == \"Minimum\")\n",
    ")\n",
    "air_temp_min = air_temp.extract_strict(cons)\n",
    "\n",
    "# extracting mean air temperature\n",
    "cons = iris.Constraint(\n",
    "    cube_func=lambda c: (len(c.cell_methods) > 0)\n",
    "    and (c.cell_methods[0].method == \"mean\")\n",
    "    and c.cell_methods[0].intervals[0] == \"1 hour\"\n",
    ")\n",
    "air_temp_mean = air_temp.extract_strict(cons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# defining Shangai region coords\n",
    "min_lat = 29.0\n",
    "max_lat = 32.0\n",
    "min_lon = 118.0\n",
    "max_lon = 123.0\n",
    "\n",
    "\n",
    "# extract data for the the Shanghai region using extract_rot_cube() function\n",
    "max_cube = extract_rot_cube(air_temp_max, min_lat, min_lon, max_lat, max_lon)\n",
    "min_cube = extract_rot_cube(air_temp_min, min_lat, min_lon, max_lat, max_lon)\n",
    "mean_cube = extract_rot_cube(air_temp_mean, min_lat, min_lon, max_lat, max_lon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print out the mean cube\n",
    "mean_cube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define start and end year for our time constraint\n",
    "start_time = 1981\n",
    "end_time = 2010\n",
    "\n",
    "# define the time constraint\n",
    "cons = iris.Constraint(time=lambda cell: start_time <= cell.point.year <= end_time)\n",
    "\n",
    "# load the data into cubes applying the time constraint\n",
    "max_cube = max_cube.extract(cons)\n",
    "min_cube = min_cube.extract(cons)\n",
    "mean_cube = mean_cube.extract(cons)\n",
    "\n",
    "# printing out the mean cube to see the change in the shape of time dimension\n",
    "mean_cube"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b>The <b>CATNIP</b> library function <b>preparation.add_time_coord_cats</b> adds a range of numeric coordinate categorisations to the cube. For more details see the <a href=\"https://metoffice.github.io/CATNIP/tools.html#preparation-module\">documentation</a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have got required cubes. Now we can add categorical coordinates to such as *year*  to the time dimension in our cubes using the CATNIP **preparation.add_time_coord_cats** function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load CATNIP's add_time_coord_cats method\n",
    "from catnip.preparation import add_time_coord_cats\n",
    "\n",
    "# Add other dimension coordinates\n",
    "max_cube = add_time_coord_cats(max_cube)\n",
    "min_cube = add_time_coord_cats(min_cube)\n",
    "mean_cube = add_time_coord_cats(mean_cube)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the *max_cube* and inspect the categorical coordinates that have been added to the time coordinate of our cube."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# printing the max_cube. Note the addtional coordinates under the Auxiliary coordinates\n",
    "max_cube"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that **add_time_coord_cats** has added a few auxiliary coordinates including the *year* coordinate to the *time* dimension.\n",
    "\n",
    "Now we can calculate maximum, minimum and mean values over the *year* coordinate using **aggregated_by** method. This will give us a single value for each year at every grid cell (30 time points in total)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate yearly max, min and mean values\n",
    "yearly_max = max_cube.aggregated_by([\"year\"], iris.analysis.MAX)\n",
    "yearly_min = min_cube.aggregated_by([\"year\"], iris.analysis.MIN)\n",
    "yearly_mean = mean_cube.aggregated_by([\"year\"], iris.analysis.MEAN)\n",
    "\n",
    "# printing the yearly_mean\n",
    "yearly_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the mean, maximum and minimum at every grid point, we can find the maximum, minimum and mean of all grid boxes to give us a timeseries which represents the aveage etc over all grid boxes in our region of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collapse longitude and latitude to get a timeseries\n",
    "yearly_max = yearly_max.collapsed(\n",
    "    [\"grid_longitude\", \"grid_latitude\"], iris.analysis.MAX\n",
    ")\n",
    "yearly_min = yearly_min.collapsed(\n",
    "    [\"grid_longitude\", \"grid_latitude\"], iris.analysis.MIN\n",
    ")\n",
    "yearly_mean = yearly_mean.collapsed(\n",
    "    [\"grid_longitude\", \"grid_latitude\"], iris.analysis.MEAN\n",
    ")\n",
    "\n",
    "# The resulting cube will be time series only\n",
    "yearly_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the *year* coordinate of max cube to see if we have the correct years for our constraint time period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(yearly_max.coord(\"year\").points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5 Calculating mean annual cycle\n",
    "\n",
    "We can calculate the mean annual cycle for precipitation_flux data over the Shanghai region for 1981-2010 (30 years) by averaging together each month(so we average all January data to get the mean for January)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract the precipitation_flux data into an iris cube from the cubelist\n",
    "pflx = cubelist.extract_strict(\"precipitation_flux\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_lat = 29.0\n",
    "max_lat = 32.0\n",
    "min_lon = 118.0\n",
    "max_lon = 123.0\n",
    "\n",
    "# extract data for the the Shanghai region using extract_rot_cube() function\n",
    "pflx_ext = extract_rot_cube(pflx, min_lat, min_lon, max_lat, max_lon)\n",
    "pflx_ext"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next extracting time constraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extracting time constraint\n",
    "start_time = 1981\n",
    "end_time = 2010\n",
    "time_cons = iris.Constraint(time=lambda cell: start_time <= cell.point.year <= end_time)\n",
    "subcube = pflx_ext.extract(time_cons)\n",
    "subcube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove auxiliary coordinates, added by extract_rot_cube but cause unecessary warnings later\n",
    "subcube.remove_coord(\"latitude\")\n",
    "subcube.remove_coord(\"longitude\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we can use the **add_time_coord_cats** method and add categorical coordinates such as *month* to the *time* dimension to our cube."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add the time categorical coordinate to cube\n",
    "subcube = add_time_coord_cats(subcube)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create mean annual cycle\n",
    "monthly_mean = subcube.aggregated_by([\"month_number\"], iris.analysis.MEAN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "monthly_mean = add_bounds(monthly_mean, [\"grid_latitude\", \"grid_longitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can calculate the area weight using **iris.analysis.cartography.area_weights** so that we weight the average to account for the fact that the areas of grid cells are variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import iris.analysis.cartography\n",
    "\n",
    "# calculate the area weight\n",
    "grid_areas = iris.analysis.cartography.area_weights(monthly_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate area averaged monthly mean rainfall\n",
    "monthly_mean = monthly_mean.collapsed(\n",
    "    [\"grid_longitude\", \"grid_latitude\"], iris.analysis.MEAN, weights=grid_areas\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.6 Visualising yearly and monthly means\n",
    "\n",
    "Let's now visualise yearly mean, max, min data for the air temperature and monthly mean data for the precipitation_flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we first need to load libraries for plotting\n",
    "import iris.plot as iplt\n",
    "import iris.quickplot as qplt\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualise yearly max, min and mean data for *air_temperature*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the timeseries for yearly max, min and mean data\n",
    "ax1 = qplt.plot(yearly_max, label=\"Max Temp\")\n",
    "ax2 = qplt.plot(yearly_min, label=\"Min Temp\")\n",
    "ax3 = qplt.plot(yearly_mean, label=\"Mean Temp\")\n",
    "plt.legend(bbox_to_anchor=(1.18, 0.78))\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualise monthly precipitation mean over the thirty years (1980-2010)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qplt.plot(monthly_mean.coord(\"month_number\"), monthly_mean, color=\"seagreen\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "    <b>Task:</b><br><ul>\n",
    "        <li>Calculate and visualise the monthly mean over Tibetan region from 1981 to 2010. Create the monthly mean with month names</li>\n",
    "        <li>Coordinates of Tibetan region: Latitude = [26 36], Longitude = [77 104]</li>\n",
    "    </ul>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# time series plot\n",
    "# write your code here .."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Calculating seasonal means<a id='season'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Calculating seasonal means: djf, mam, jja and son\n",
    "\n",
    "Calculate mean wind speed and wind direction from 1981 to 2010 for different seasons over the entire domain.\n",
    "\n",
    "First we need to calculate wind speed and wind direction. In previous tutorial, we calculated the wind speed using hard coded simple arithmetic operations. In this tutorial, we will use catnip's **windspeed** and **wind_direction** methods.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define constraints for x_wind and y_wind data\n",
    "xcons = iris.Constraint(\n",
    "    cube_func=lambda c: c.standard_name == \"x_wind\"\n",
    "    and (\"pressure\" not in [coord.name() for coord in c.coords()])\n",
    ")\n",
    "ycons = iris.Constraint(\n",
    "    cube_func=lambda c: c.standard_name == \"y_wind\"\n",
    "    and (\"pressure\" not in [coord.name() for coord in c.coords()])\n",
    ")\n",
    "\n",
    "# apply the constraint and load the x_wind and y_wind data\n",
    "u = cubelist.extract_strict(xcons)\n",
    "v = cubelist.extract_strict(ycons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define time constraint\n",
    "start_time = 1981\n",
    "end_time = 2010\n",
    "cons = iris.Constraint(time=lambda cell: start_time <= cell.point.year <= end_time)\n",
    "u = u.extract(cons)\n",
    "v = v.extract(cons)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now import and use catnip's **windspeed** and **wind_direction** methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import catnip methods\n",
    "from catnip.analysis import windspeed\n",
    "from catnip.analysis import wind_direction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate windspeed and wind direction\n",
    "wind_speed_cube = windspeed(u, v)\n",
    "wind_direction_cube = wind_direction(u, v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add coordinates and extract data for different seasons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from catnip.preparation import add_time_coord_cats\n",
    "\n",
    "wind_speed_cube = add_time_coord_cats(wind_speed_cube)\n",
    "wind_direction_cube = add_time_coord_cats(wind_direction_cube)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the windspeed data for all four seasons\n",
    "wndspd_djf = wind_speed_cube.extract(iris.Constraint(season=\"djf\"))\n",
    "wndspd_mam = wind_speed_cube.extract(iris.Constraint(season=\"mam\"))\n",
    "wndspd_jja = wind_speed_cube.extract(iris.Constraint(season=\"jja\"))\n",
    "wndspd_son = wind_speed_cube.extract(iris.Constraint(season=\"son\"))\n",
    "\n",
    "# Extract the wind direction data for the all four season\n",
    "wnddir_djf = wind_direction_cube.extract(iris.Constraint(season=\"djf\"))\n",
    "wnddir_mam = wind_direction_cube.extract(iris.Constraint(season=\"mam\"))\n",
    "wnddir_jja = wind_direction_cube.extract(iris.Constraint(season=\"jja\"))\n",
    "wnddir_son = wind_direction_cube.extract(iris.Constraint(season=\"son\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets print out wnddir_son cube to see its chape\n",
    "wnddir_son"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate seasonal means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the windspeed mean over the seasons\n",
    "wspd_djf_mean = wndspd_djf.aggregated_by([\"season\"], iris.analysis.MEAN)\n",
    "wspd_mam_mean = wndspd_mam.aggregated_by([\"season\"], iris.analysis.MEAN)\n",
    "wspd_jja_mean = wndspd_jja.aggregated_by([\"season\"], iris.analysis.MEAN)\n",
    "wspd_son_mean = wndspd_son.aggregated_by([\"season\"], iris.analysis.MEAN)\n",
    "\n",
    "# calculate the wind direction mean over the seasons\n",
    "wndir_djf_mean = wnddir_djf.aggregated_by([\"season\"], iris.analysis.MEAN)\n",
    "wndir_mam_mean = wnddir_mam.aggregated_by([\"season\"], iris.analysis.MEAN)\n",
    "wndir_jja_mean = wnddir_jja.aggregated_by([\"season\"], iris.analysis.MEAN)\n",
    "wndir_son_mean = wnddir_son.aggregated_by([\"season\"], iris.analysis.MEAN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets print out wndir_son_mean cube to see its chape\n",
    "wndir_son_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now visualise seasonal means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we first need to load libraries for plotting\n",
    "import iris.plot as iplt\n",
    "import iris.quickplot as qplt\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of seasonal cubes to loop through\n",
    "seasonal_cubes = [wspd_djf_mean, wspd_mam_mean, wspd_jja_mean, wspd_son_mean]\n",
    "\n",
    "# here we create a list of mean cubes and we will run for loop over it.\n",
    "\n",
    "# set a figure big enough to hold the subplots\n",
    "plt.figure(figsize=(10, 10))\n",
    "\n",
    "# loop through the seaonal cube list and plot the data\n",
    "# len functions returns the length of the list.\n",
    "for i in range(len(seasonal_cubes)):\n",
    "    plt.subplot(2, 2, i + 1)\n",
    "    # above line will create the 2 rows and 2 cols of subplots and i indicates the position of subplot.\n",
    "    # i starts from 0, that is why we increment 1 in it.\n",
    "    # plot the windspeed at the first timestep\n",
    "    qplt.contourf(seasonal_cubes[i][0, :, :])\n",
    "    # add some coastlines for context\n",
    "    plt.gca().coastlines()\n",
    "    # get the season name from the coordinate\n",
    "    season = seasonal_cubes[i].coord(\"season\").points[0]\n",
    "    # add the name as plot's title\n",
    "    plt.title(\"Season: \" + season)\n",
    "    plt.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot in any other project. Python cartopy library provide us the range of option of ploting in different projection. See the list of option in [cartopy documentation](https://scitools.org.uk/cartopy/docs/latest/crs/projections.html). Lets try the above plot with **PlateCarree** projection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "\n",
    "# list of seasonal cubes to loop through\n",
    "seasonal_cubes = [wspd_djf_mean, wspd_mam_mean, wspd_jja_mean, wspd_son_mean]\n",
    "\n",
    "# set a figure big enough to hold the subplots\n",
    "plt.figure(figsize=(10, 10))\n",
    "\n",
    "# loop through the seaonal cube list and plot the data\n",
    "for i in range(len(seasonal_cubes)):\n",
    "    plt.subplot(2, 2, i + 1, projection=ccrs.PlateCarree())\n",
    "    # plot the windspeed at the first timestep\n",
    "    qplt.contourf(seasonal_cubes[i][0, :, :])\n",
    "    # add some coastlines for context\n",
    "    plt.gca().coastlines()\n",
    "    # get the season name from the coordinate\n",
    "    season = seasonal_cubes[i].coord(\"season\").points[0]\n",
    "    # add the name as plot's title\n",
    "    plt.title(\"Season: \" + season)\n",
    "    plt.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "    <b>Task:</b><br><ul>\n",
    "        <li>Calculate and visualise the seasonal mean of <b>surface temperature</b> over Tibatan region from 1981 to 2010.</li>\n",
    "        <li>Coordinates of Tibetan region: Latitude = [26 36], Longitude = [77 104]</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate seasonal means\n",
    "# write your code here .."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualising seasonal means\n",
    "# write your code here .."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Calculating differences<a id='percent'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 mean surface temperature diffference in winter season (dec, jan, feb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can find the difference of mean surface temperature between the past(1851-1880) and recent(1981-2010) 30 years periods.\n",
    "\n",
    "First, we need to extract out desired data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract air_temperature\n",
    "sft = cubelist.extract_strict(\"surface_temperature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# constraints: two 30 years periods - past and presnet\n",
    "past_cons = iris.Constraint(time=lambda cell: 1851 <= cell.point.year <= 1880)\n",
    "present_cons = iris.Constraint(time=lambda cell: 1981 <= cell.point.year <= 2010)\n",
    "past = sft.extract(past_cons)\n",
    "present = sft.extract(present_cons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load catnip's add_time_coord_cats method\n",
    "from catnip.preparation import add_time_coord_cats\n",
    "\n",
    "# Add other dimension coordinates\n",
    "past = add_time_coord_cats(past)\n",
    "present = add_time_coord_cats(present)\n",
    "\n",
    "# Extract the winter season\n",
    "past_djf = past.extract(iris.Constraint(season=\"djf\"))\n",
    "present_djf = present.extract(iris.Constraint(season=\"djf\"))\n",
    "\n",
    "# extract data for Shanghai region\n",
    "past_djf = extract_rot_cube(past_djf, min_lat, min_lon, max_lat, max_lon)\n",
    "present_djf = extract_rot_cube(present_djf, min_lat, min_lon, max_lat, max_lon)\n",
    "\n",
    "# calculate 30 year mean of winter season\n",
    "past_djf = past_djf.aggregated_by([\"season\"], iris.analysis.MEAN)\n",
    "present_djf = present_djf.aggregated_by([\"season\"], iris.analysis.MEAN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now got our cubes for different climatological periods. We now calcuate the difference by subtracting the past data form present using **iris.analysis.math.subtract** method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "djf_diff = iris.analysis.maths.subtract(present_djf, past_djf)\n",
    "djf_diff.rename(\"surface temperature difference: Winter\")\n",
    "past_djf.rename(\"surface temperature past climate: Winter 1851-1880 \")\n",
    "present_djf.rename(\"surface temperature present climate: Winter 1981-2010 \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add bounds to the cubes\n",
    "past_djf = add_bounds(past_djf, [\"grid_latitude\", \"grid_longitude\"])\n",
    "present_djf = add_bounds(present_djf, [\"grid_latitude\", \"grid_longitude\"])\n",
    "djf_diff = add_bounds(djf_diff, [\"grid_latitude\", \"grid_longitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "</pre>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b> iris.analysis.math provides a range of mathematical and statistical operations. See the <a href=\"https://scitools.org.uk/iris/docs/v1.9.1/iris/iris/analysis/maths.html\">documentation for more information</a>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now visualise the difference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of our djf cubes and diff cube to loop through\n",
    "seasonal_cubes = [present_djf, past_djf, djf_diff]\n",
    "plt.figure(figsize=(15, 10))\n",
    "# loop through the seaonal cube list and plot the data\n",
    "for i in range(len(seasonal_cubes)):\n",
    "    plt.subplot(2, 2, i + 1)\n",
    "    # plot the windspeed at the first timestep\n",
    "    if i == 2:\n",
    "        qplt.pcolormesh(\n",
    "            seasonal_cubes[i][0, :, :], cmap=plt.cm.get_cmap(\"Reds\"), vmin=0, vmax=2\n",
    "        )\n",
    "    else:\n",
    "        qplt.pcolormesh(seasonal_cubes[i][0, :, :], vmin=277.5, vmax=289)\n",
    "\n",
    "    # add some coastlines for context\n",
    "    plt.gca().coastlines()\n",
    "    # get the season name from the coordinate\n",
    "    season = seasonal_cubes[i].coord(\"season\").points[0]\n",
    "    # add the name as plot's title\n",
    "    plt.title(seasonal_cubes[i].name())\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Percentage difference in winter precipitaition "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also calculate the percentage difference. \n",
    "\n",
    "Let's calculate the change in mean precipitation from a past 30 year period (1851-1880) to the most recent 30 years (1981-2010)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract precipitation flux\n",
    "pflx = cubelist.extract_strict(\"precipitation_flux\")\n",
    "\n",
    "# extract the time contraints\n",
    "past_cons = iris.Constraint(time=lambda cell: 1851 <= cell.point.year <= 1880)\n",
    "present_cons = iris.Constraint(time=lambda cell: 1981 <= cell.point.year <= 2010)\n",
    "past = pflx.extract(past_cons)\n",
    "present = pflx.extract(present_cons)\n",
    "\n",
    "# Add other dimension coordinates\n",
    "past = add_time_coord_cats(past)\n",
    "present = add_time_coord_cats(present)\n",
    "\n",
    "# Extract the precipitation data for the winter season\n",
    "past_djf = past.extract(iris.Constraint(season=\"djf\"))\n",
    "present_djf = present.extract(iris.Constraint(season=\"djf\"))\n",
    "\n",
    "# extract data for Shanghai region\n",
    "past_djf = extract_rot_cube(past_djf, min_lat, min_lon, max_lat, max_lon)\n",
    "present_djf = extract_rot_cube(present_djf, min_lat, min_lon, max_lat, max_lon)\n",
    "\n",
    "# calculate the means\n",
    "past_djf = past_djf.aggregated_by([\"season\"], iris.analysis.MEAN)\n",
    "present_djf = present_djf.aggregated_by([\"season\"], iris.analysis.MEAN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now calculate the difference using **subtract** function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "djf_diff = iris.analysis.maths.subtract(present_djf, past_djf)\n",
    "djf_diff.rename(\"precipitation flux difference: Winter\")\n",
    "past_djf.rename(\"precipitation flux past climate: Winter 1851-1880 \")\n",
    "present_djf.rename(\"precipitation flux present climate: Winter 1981-2010 \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add bounds to the cubes\n",
    "past_djf = add_bounds(past_djf, [\"grid_latitude\", \"grid_longitude\"])\n",
    "present_djf = add_bounds(present_djf, [\"grid_latitude\", \"grid_longitude\"])\n",
    "djf_diff = add_bounds(djf_diff, [\"grid_latitude\", \"grid_longitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calcuate the percentage difference, we can use **analysis.maths.multiply** and **iris.analysis.maths.divide** to calculate percentage change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the percentage change\n",
    "pcent_change = iris.analysis.maths.multiply(\n",
    "    iris.analysis.maths.divide(djf_diff, past_djf), 100\n",
    ")\n",
    "\n",
    "# remember to change the title and units to reflect the data processing\n",
    "pcent_change.rename(\"precipitation flux percent difference\")\n",
    "pcent_change.units = \"%\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using iris plot for better colour saturation!\n",
    "import iris.plot as iplt\n",
    "\n",
    "# list of our winter cubes and the diff cube to loop through for plotting\n",
    "seasonal_cubes = [present_djf, past_djf, pcent_change]\n",
    "plt.figure(figsize=(15, 10))\n",
    "# loop through the seaonal cube list and plot the data\n",
    "for i in range(len(seasonal_cubes)):\n",
    "    plt.subplot(2, 2, i + 1)\n",
    "    # plot the windspeed at the first timestep\n",
    "    if i == 2:\n",
    "        qplt.pcolormesh(seasonal_cubes[i][0, :, :], cmap=plt.cm.get_cmap(\"RdBu\"))\n",
    "    else:\n",
    "        qplt.pcolormesh(seasonal_cubes[i][0, :, :])\n",
    "\n",
    "    # add some coastlines for context\n",
    "    plt.gca().coastlines()\n",
    "    # get the season name from the coordinate\n",
    "    season = seasonal_cubes[i].coord(\"season\").points[0]\n",
    "    # add the name as plot's title\n",
    "    plt.title(seasonal_cubes[i].name())\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "    <b>Task:</b><br><ul>\n",
    "        <li>Calculate mean surface temperature difference over Tibetan region from past 30 years (1851-1880) to present 30 years (1981-2010).\n",
    "        <li>Coordinates of Tibetan region: Latitude = [26 36], Longitude = [77 104]</li>\n",
    "    </ul>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Exercises<a id='exercise'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this exercise we will analyse the mean air temperature from past 30 years (1851-1880) to present 30 years (1981-2010), over the Shanghai region, in all four seasons. Visualize past, present and difference in a row."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1: Load monthly data and constraint time and region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excercise 2: Calculate seasonal mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excercise 3: Visualise the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Summary</b><br> \n",
    "    In this session we learned how:<br>\n",
    "    <ul>\n",
    "        <li>to calculate yearly and monthly means</li>\n",
    "        <li>to calculate seasonal means and differences</li>\n",
    "        <li>to visualize the results</li>\n",
    "    </ul>\n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "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": {
  "kernelspec": {
   "display_name": "cssp37",
   "language": "python",
   "name": "cssp37"
  },
  "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}