notebooks/bias-adjustment_R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "# Atlas notebooks\n",
    "\n",
    "***\n",
    "\n",
    "> This notebook reproduces and extends parts of the figures and products of the AR6-WGI Atlas. It is part of a notebook collection available at https://github.com/IPCC-WG1/Atlas for reproducibility and reusability purposes. This work is licensed under a [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0).\n",
    ">\n",
    "> ![Creative Commons License >](https://i.creativecommons.org/l/by/4.0/88x31.png)\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(bias-adjustment)=\n",
    "# Bias adjustment of model temperatures for the calculation of extreme indices\n",
    "29/6/2021\n",
    "\n",
    "**A. Casanueva** ([Santander Meteorology Group](http://meteo.unican.es). Dept. Applied Mathematics and Compute Science, University of Cantabria, Santander, Spain)\n",
    "\n",
    "Bias adjustment is used in the IPCC interactive Atlas for the calculation of absolute threshold-based indices, namely the number of days with maximum temperature above 35ºC (TX35) and above 40ºC (TX40), and the number of frost days (minimum temperature below 0ºC). More details about bias adjustment can be found in the Atlas chapter from the IPCC-AR6 (Cross-Chapter Box 10.2: \"Issues in bias adjustment\") and in {cite:authorpar}`TestingBiasAdjustment`.\n",
    "\n",
    "Among the wide variety of methods typically used by the climate community for bias adjustment, this notebook focuses on ISIMIP3 ({cite:authorpar}`GMDTrendpreservingBias`), a parametric quantile mapping alternative which has been designed to robustly adjust biases in all percentiles of a distribution whilst preserving their trends. In particular, we use ISIMIP3 in this example to bias-adjust temperature from the r12i1p1 run of the EC-EARTH global climate model, included in the fifth phase of the Coupled Model Intercomparison Project (CMIP5), over the Iberian Peninsula using as reference the W5E5 observational dataset ({cite}`cucchiWFDE5Bias`). We focus on the boreal summer season (June-July-August) for the calculation of bias-adjusted TX35."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### Loading packages\n",
    "This notebook is based on the R programming language and requires the following {cite}`climate4r` packages:\n",
    "\n",
    " * `loadeR` to load climate model and observational datasets ({cite:authorpar}`RbasedClimate4ROpen`).\n",
    " * `climate4R.UDG` to access harmonized climate data via the [Santander Meteorology Group UDG-TAP](http://meteo.unican.es/udg-tap/).\n",
    " * `transformeR` to manipulate climate data ({cite:authorpar}`RbasedClimate4ROpen`).\n",
    " * `downscaleR` to adjust systematic model biases ({cite:authorpar}`GMDStatisticalDownscaling`).\n",
    " * `climate4R.indices` to compute climate indices ({cite:authorpar}`RbasedClimate4ROpen`).\n",
    " * `visualizeR` to produce graphical representations of the results ({cite:authorpar}`PackageVisualizeCommunicate`).\n",
    " \n",
    "Along with:\n",
    "\n",
    " * `rgdal` to work with spatial data ({cite:authorpar}`CRANPackageRgdal`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: rJava\n",
      "\n",
      "Loading required package: loadeR.java\n",
      "\n",
      "Java version 11x amd64 by Azul Systems, Inc. detected\n",
      "\n",
      "NetCDF Java Library v4.6.0-SNAPSHOT (23 Apr 2015) loaded and ready\n",
      "\n",
      "Loading required package: climate4R.UDG\n",
      "\n",
      "climate4R.UDG version 0.2.4 (2022-06-15) is loaded\n",
      "\n",
      "Please use 'citation(\"climate4R.UDG\")' to cite this package.\n",
      "\n",
      "loadeR version 1.7.1 (2021-07-05) is loaded\n",
      "\n",
      "Please use 'citation(\"loadeR\")' to cite this package.\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "    _______   ____  ___________________  __  ________ \n",
      "   / ___/ /  / /  |/  / __  /_  __/ __/ / / / / __  / \n",
      "  / /  / /  / / /|_/ / /_/ / / / / __/ / /_/ / /_/_/  \n",
      " / /__/ /__/ / /  / / __  / / / / /__ /___  / / \\ \\ \n",
      " \\___/____/_/_/  /_/_/ /_/ /_/  \\___/    /_/\\/   \\_\\ \n",
      " \n",
      "      github.com/SantanderMetGroup/climate4R\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "transformeR version 2.1.5 (2022-05-12) is loaded\n",
      "\n",
      "WARNING: Your current version of transformeR (v2.1.5) is ahead of the master branch version (2.1.4)\n",
      "\n",
      "Development version may have an unexpected behaviour\n",
      "\n",
      "Please see 'citation(\"transformeR\")' to cite this package.\n",
      "\n",
      "downscaleR version 3.3.3 (2021-07-05) is loaded\n",
      "\n",
      "Please use 'citation(\"downscaleR\")' to cite this package.\n",
      "\n",
      "climate4R.indices version 0.2.0 (2021-07-08) is loaded\n",
      "\n",
      "Use 'indexShow()' for an overview of the available climate indices and circIndexShow() for the circulation indices.\n",
      "\n",
      "NOTE: use package climate4R.climdex to calculate ETCCDI indices.\n",
      "\n",
      "\n",
      "Attaching package: ‘climate4R.indices’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:transformeR’:\n",
      "\n",
      "    lambWT\n",
      "\n",
      "\n",
      "visualizeR version 1.6.1 (2021-03-11) is loaded\n",
      "\n",
      "Please see 'citation(\"visualizeR\")' to cite this package.\n",
      "\n",
      "Loading required package: sp\n",
      "\n",
      "rgdal: version: 1.5-16, (SVN revision 1050)\n",
      "Geospatial Data Abstraction Library extensions to R successfully loaded\n",
      "Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29\n",
      "Path to GDAL shared files: /home/phanaur/mambaforge/envs/prueba2/share/gdal\n",
      "GDAL binary built with GEOS: TRUE \n",
      "Loaded PROJ runtime: Rel. 7.2.0, November 1st, 2020, [PJ_VERSION: 720]\n",
      "Path to PROJ shared files: /home/phanaur/mambaforge/envs/prueba2/share/proj\n",
      "PROJ CDN enabled: TRUE\n",
      "Linking to sp version:1.4-5\n",
      "To mute warnings of possible GDAL/OSR exportToProj4() degradation,\n",
      "use options(\"rgdal_show_exportToProj4_warnings\"=\"none\") before loading rgdal.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(loadeR) \n",
    "library(climate4R.UDG) \n",
    "library(transformeR)\n",
    "library(downscaleR)\n",
    "library(climate4R.indices)\n",
    "library(visualizeR)\n",
    "library(rgdal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### Defining the parameters of the experiment\n",
    "We start by setting the general parameters that define the experiment to perform: spatial extent (Iberian Peninsula), season (June-July-August), baseline/historical period (1986-2005) and future period/emission scenario (2041-2060, RCP8.5) of interest. We set likewise the codes which identify the desired datasets in the {cite:authorpar}`UDGTAPHome` , here W5E5 and EC-EARTH (historical and RCP8.5). Check other available, harmonized datasets with `UDG.datasets`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "lons <- c(-9.25, 3.5)  # Iberian Peninsula\n",
    "lats <- c(36, 44)   # Iberian Peninsula\n",
    "\n",
    "season <- 6:8  # June-July-August\n",
    "years.hist <- 1986:2005\n",
    "years.rcp <- 2041:2060\n",
    "\n",
    "dataset.obs <- \"W5E5\" \n",
    "dataset.hist <- \"CMIP5_EC-EARTH_r12i1p1_historical\"\n",
    "dataset.rcp <- \"CMIP5_EC-EARTH_r12i1p1_rcp85\""
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### Loading datasets \n",
    "Next, we will see two alternative ways to load in memory the datasets that are needed to run this experiment. On the one hand, both W5E5 and EC-EARTH can be accessed from the {cite:authorpar}`UDGTAPHome`, a THREDDS-based service that provides access to a wide catalogue of popular climate datasets. On the other hand, for the sake of time, these data have been also stored locally as `NetCDF` files and can be directly loaded from the `auxiliary-material` folder. In both cases, the datasets can be easily accessed via the `loadGridData` function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Online from the User Data Gateway - Thredds Access Portal (UDG-TAP)\n",
    "\n",
    "> This section takes long to execute due to intensive remote data loading. See the next section ([From local files](#From-local-files)) to load these data from a local copy available in the repository and continue with the notebook operations.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Previously defined parameters are used to load daily observed (W5E5) and simulated (EC-EARTH) temperatures from the UDG-TAP. For EC-EARTH, both historical and future (i.e. RCP8.5) simulations are needed. It is important to note that not only mean, but also minimum and maximum temperatures are required by the ISIMIP3 method to adjust maximum and minimum temperature, which are obtained from the correction of the amplitude and skewness of the diurnal temperature cycle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```R\n",
    "\n",
    "load.data <- function (dset, years, var) loadGridData(dataset = dset, var = var,\n",
    "                                                 years = years,\n",
    "                                                 latLim = lats, lonLim = lons,\n",
    "                                                 season = season) \n",
    "# Loading mean temperature\n",
    "y.tas <- load.data(dataset.obs, years.hist, \"tas\")\n",
    "x.tas <- load.data(dataset.hist, years.hist, \"tas\")\n",
    "newdata.tas <- load.data(dataset.rcp, years.rcp, \"tas\")\n",
    "\n",
    "# Loading minimum temperature\n",
    "y.tasmin <- load.data(dataset.obs, years.hist, \"tasmin\")\n",
    "x.tasmin <- load.data(dataset.hist, years.hist, \"tasmin\")\n",
    "newdata.tasmin <- load.data(dataset.rcp, years.rcp, \"tasmin\")\n",
    "\n",
    "# Loading maximum temperature\n",
    "y.tasmax <- load.data(dataset.obs, years.hist, \"tasmax\")\n",
    "x.tasmax <- load.data(dataset.hist, years.hist, \"tasmax\")\n",
    "newdata.tasmax <- load.data(dataset.rcp, years.rcp, \"tasmax\")\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Due to the spatial mismatch between model and observations, we upscale W5E5 from its original resolution (0.5º) to a 1º regular grid similar to that of EC-EARTH (1.125º). See {cite:authorpar}`TestingBiasAdjustment` for a discussion on this procedure. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "```R\n",
    "upscale <- function(grid) redim(upscaleGrid(grid, times = 2, aggr.fun = list(FUN = \"mean\", na.rm = TRUE)), \n",
    "                                drop = TRUE)\n",
    "y.tas <- upscale(y.tas)\n",
    "y.tasmax <- upscale(y.tasmax)\n",
    "y.tasmin <- upscale(y.tasmin)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  From local files\n",
    "\n",
    "Alternatively, these data (both W5E5 and EC-EARTH, for the Iberian Peninsula and the selected time slices) can be directly loaded from the `auxiliary-material` folder by using the `loadGridData` function (skip this section if you have sucessfully loaded your data from the UDG-TAP). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2023-01-16 11:10:48] Defining geo-location parameters\n",
      "\n",
      "[2023-01-16 11:10:48] Defining time selection parameters\n",
      "\n",
      "[2023-01-16 11:10:48] Retrieving data subset ...\n",
      "\n",
      "[2023-01-16 11:10:48] Done\n",
      "\n",
      "[2023-01-16 11:10:52] Defining geo-location parameters\n",
      "\n",
      "[2023-01-16 11:10:52] Defining time selection parameters\n",
      "\n",
      "[2023-01-16 11:10:52] Retrieving data subset ...\n",
      "\n",
      "[2023-01-16 11:10:52] Done\n",
      "\n",
      "[2023-01-16 11:10:56] Defining geo-location parameters\n",
      "\n",
      "[2023-01-16 11:10:56] Defining time selection parameters\n",
      "\n",
      "[2023-01-16 11:10:56] Retrieving data subset ...\n",
      "\n",
      "[2023-01-16 11:10:57] Done\n",
      "\n",
      "[2023-01-16 11:11:01] Defining geo-location parameters\n",
      "\n",
      "[2023-01-16 11:11:01] Defining time selection parameters\n",
      "\n",
      "[2023-01-16 11:11:01] Retrieving data subset ...\n",
      "\n",
      "[2023-01-16 11:11:01] Done\n",
      "\n",
      "[2023-01-16 11:11:05] Defining geo-location parameters\n",
      "\n",
      "[2023-01-16 11:11:05] Defining time selection parameters\n",
      "\n",
      "[2023-01-16 11:11:05] Retrieving data subset ...\n",
      "\n",
      "[2023-01-16 11:11:05] Done\n",
      "\n",
      "[2023-01-16 11:11:09] Defining geo-location parameters\n",
      "\n",
      "[2023-01-16 11:11:09] Defining time selection parameters\n",
      "\n",
      "[2023-01-16 11:11:09] Retrieving data subset ...\n",
      "\n",
      "[2023-01-16 11:11:09] Done\n",
      "\n",
      "[2023-01-16 11:11:13] Defining geo-location parameters\n",
      "\n",
      "[2023-01-16 11:11:13] Defining time selection parameters\n",
      "\n",
      "[2023-01-16 11:11:13] Retrieving data subset ...\n",
      "\n",
      "[2023-01-16 11:11:14] Done\n",
      "\n",
      "[2023-01-16 11:11:18] Defining geo-location parameters\n",
      "\n",
      "[2023-01-16 11:11:18] Defining time selection parameters\n",
      "\n",
      "[2023-01-16 11:11:18] Retrieving data subset ...\n",
      "\n",
      "[2023-01-16 11:11:18] Done\n",
      "\n",
      "[2023-01-16 11:11:22] Defining geo-location parameters\n",
      "\n",
      "[2023-01-16 11:11:22] Defining time selection parameters\n",
      "\n",
      "[2023-01-16 11:11:22] Retrieving data subset ...\n",
      "\n",
      "[2023-01-16 11:11:22] Done\n",
      "\n"
     ]
    }
   ],
   "source": [
    "file.obs <- \"auxiliary-material/W5E5_IP_tas-tasmin-tasmax_1986-2005_JJA.nc\"\n",
    "file.hist <- \"auxiliary-material/CMIP5_EC-EARTH_r12i1p1_historical_IP_tas-tasmin-tasmax_1986-2005_JJA.nc\"\n",
    "file.rcp <- \"auxiliary-material/CMIP5_EC-EARTH_r12i1p1_rcp85_IP_tas-tasmin-tasmax_2041-2060_JJA.nc\"\n",
    "\n",
    "# Loading mean temperature\n",
    "y.tas <- loadGridData(file.obs, var = \"tas\")\n",
    "x.tas <- loadGridData(file.hist, var = \"tas\")\n",
    "newdata.tas <- loadGridData(file.rcp, var = \"tas\")\n",
    "\n",
    "# Loading minimum temperature\n",
    "y.tasmin <- loadGridData(file.obs, var = \"tasmin\")\n",
    "x.tasmin <- loadGridData(file.hist, var = \"tasmin\")\n",
    "newdata.tasmin <- loadGridData(file.rcp, var = \"tasmin\")\n",
    "\n",
    "# Loading maximum temperature\n",
    "y.tasmax <- loadGridData(file.obs, var = \"tasmax\")\n",
    "x.tasmax <- loadGridData(file.hist, var = \"tasmax\")\n",
    "newdata.tasmax <- loadGridData(file.rcp, var = \"tasmax\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### Model biases\n",
    "\n",
    "We continue by exploring EC-EARTH systematic biases compared to W5E5. For this purpose we interpolate the historical EC-EARTH simulations (recall, at 1.125º) to the W5E5's grid (at 1º) with `interpGrid`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2023-01-16 11:11:22] Calculating nearest neighbors...\n",
      "\n",
      "[2023-01-16 11:11:22] Performing nearest interpolation... may take a while\n",
      "\n",
      "[2023-01-16 11:11:22] Done\n",
      "\n",
      "[2023-01-16 11:11:22] Calculating nearest neighbors...\n",
      "\n",
      "[2023-01-16 11:11:22] Performing nearest interpolation... may take a while\n",
      "\n",
      "[2023-01-16 11:11:22] Done\n",
      "\n",
      "[2023-01-16 11:11:22] Calculating nearest neighbors...\n",
      "\n",
      "[2023-01-16 11:11:22] Performing nearest interpolation... may take a while\n",
      "\n",
      "[2023-01-16 11:11:22] Done\n",
      "\n"
     ]
    }
   ],
   "source": [
    "hist.tas <- interpGrid(x.tas, new.coordinates = getGrid(y.tas), method = \"nearest\")\n",
    "hist.tasmin <- interpGrid(x.tasmin, new.coordinates = getGrid(y.tasmin), method = \"nearest\")\n",
    "hist.tasmax <- interpGrid(x.tasmax, new.coordinates = getGrid(y.tasmax), method = \"nearest\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note that nearest neighbour interpolation is internally performed within the `biasCorrection` function (which will be used later) when observations and model grids do not match. Nevertheless, we explicitly do it here for biases calculation. Mean bias of daily mean, maximum and minimum temperature during the historical period are computed with `gridArithmetics`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "bias.tas <- gridArithmetics(climatology(hist.tas), climatology(y.tas), operator = \"-\")\n",
    "bias.tasmax <- gridArithmetics(climatology(hist.tasmax), climatology(y.tasmax), operator = \"-\")\n",
    "bias.tasmin <- gridArithmetics(climatology(hist.tasmin), climatology(y.tasmin), operator = \"-\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are plotted with the `spatialPlot` function from the `visualizeR` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "spatialPlot(suppressMessages(makeMultiGrid(bias.tas, bias.tasmax, bias.tasmin, skip.temporal.check = TRUE)), \n",
    "            backdrop.theme = \"coastline\", \n",
    "            names.attr = c(\"Mean temperature\", \"Maximum Temperature\", \"Minimum Temperature\"),\n",
    "            main = \"EC-EARTH biases (degC)\", layout = c(1, 3), as.table = TRUE, \n",
    "            set.min = -6, set.max = 6, at = seq(-6, 6, 12/20), rev.colors = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, we can also visualize the EC-EARTH bias for TX35 (number of days per year with maximum temperature above 35º) during the historical period. This index can be easily computed with the `indexGrid` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Observed TX35, year by year\n",
    "index.obs <- redim(indexGrid(tx = y.tasmax, index.code = \"TXth\", th = 35, time.resolution = \"year\"), \n",
    "                   drop = TRUE)\n",
    "# Simulated TX35, year by year\n",
    "index.raw <- redim(indexGrid(tx = hist.tasmax, index.code = \"TXth\", th = 35, time.resolution = \"year\"), \n",
    "                   drop = TRUE)\n",
    "\n",
    "tx35.obs <- climatology(index.obs)  # Mean value\n",
    "tx35.raw <- climatology(index.raw)  # Mean value\n",
    "bias.tx35 <- gridArithmetics(tx35.raw, tx35.obs, operator = \"-\")  # Bias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "spatialPlot(bias.tx35, backdrop.theme = \"coastline\", main = \"EC-EARTH bias in TX35 (days/year)\", \n",
    "            set.min = -30, set.max = 30, \n",
    "            at = seq(-30, 30, 60/20),\n",
    "            rev.colors = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### Bias adjustment\n",
    "Three individual adjustments need to be done to get bias-adjusted daily minimum and maximum temperatures with the ISIMIP3 method:\n",
    "\n",
    "1. Bias adjustment of mean temperature\n",
    "2. Bias adjustment of temperature range\n",
    "3. Bias adjustment of temperature skewness\n",
    "\n",
    "In this example, the three steps are next undertaken to correct both historical and future model simulations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by calculating the temperature range and skewness for the observations and (historical and future) simulations using the `gridArithmetics` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Range\n",
    "y.range <- gridArithmetics(y.tasmax, y.tasmin, operator = \"-\")\n",
    "x.range <- gridArithmetics(x.tasmax, x.tasmin, operator = \"-\")\n",
    "newdata.range <- gridArithmetics(newdata.tasmax, newdata.tasmin, operator = \"-\")\n",
    "\n",
    "# Skewness\n",
    "y.skew <- gridArithmetics(gridArithmetics(y.tas, y.tasmin, operator = \"-\"), y.range, operator = \"/\")\n",
    "x.skew <- gridArithmetics(gridArithmetics(x.tas,x.tasmin, operator = \"-\"), x.range, operator = \"/\")\n",
    "newdata.skew <- gridArithmetics(gridArithmetics(newdata.tas, newdata.tasmin, operator = \"-\"), \n",
    "                                newdata.range, operator = \"/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Next, the `biasCorrection` function is used to perform the three steps mentioned above, by indicating the ISIMIP3 method and related arguments as inputs. Other (methodologically simpler) methods are also available in `biasCorrection`; the reader is referred to the function documentation for more details.\n",
    "\n",
    "1.Bias adjustment of daily mean temperature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# List of arguments that have to be passed to the \"biasCorrection\" function when using the ISIMIP3 method\n",
    "isimip3.args = list(lower_bound = NULL,lower_threshold = NULL, upper_bound = NULL, \n",
    "                    upper_threshold = NULL,  randomization_seed =  NULL, \n",
    "                    detrend = array(data = TRUE, dim = 1), rotation_matrices = NULL, \n",
    "                    n_quantiles = 50, distribution = \"normal\", \n",
    "                    trend_preservation = array(data = \"additive\", dim = 1), \n",
    "                    adjust_p_values = array(data = FALSE, dim = 1), if_all_invalid_use = NULL, \n",
    "                    invalid_value_warnings = FALSE)\n",
    "\n",
    "# Adjusting historical simulations\n",
    "bc.tas.hist <- biasCorrection(y = y.tas, x = x.tas, newdata = x.tas, \"precipitation\" = FALSE, \n",
    "                         method = \"isimip3\", isimip3.args = isimip3.args)\n",
    "# Adjusting future simulations\n",
    "bc.tas <- biasCorrection(y = y.tas, x = x.tas, newdata = newdata.tas, \"precipitation\" = FALSE, \n",
    "                         method = \"isimip3\", isimip3.args = isimip3.args)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "2.Bias ajustment of the daily temperature range:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# List of arguments that have to be passed to the \"biasCorrection\" function when using the ISIMIP3 method\n",
    "isimip3.range.args = list(lower_bound = 0, lower_threshold = 0.01, upper_bound = NULL, \n",
    "                            upper_threshold = NULL, randomization_seed = NULL, \n",
    "                            detrend = array(data = FALSE, dim = 1), rotation_matrices =  NULL, \n",
    "                            n_quantiles = 50, distribution = \"rice\", \n",
    "                            trend_preservation = array(data = \"mixed\", dim = 1),\n",
    "                            adjust_p_values = array(data = FALSE, dim = 1), if_all_invalid_use = NULL, \n",
    "                            invalid_value_warnings = FALSE)\n",
    "\n",
    "# Adjusting historical simulations\n",
    "bc.range.hist <- biasCorrection(y = y.range, x = x.range, newdata = x.range, \"precipitation\" = FALSE, \n",
    "                                method = \"isimip3\", isimip3.args = isimip3.range.args)\n",
    "# Adjusting future simulations                           \n",
    "bc.range <- biasCorrection(y = y.range, x = x.range, newdata = newdata.range, \"precipitation\" = FALSE, \n",
    "                           method = \"isimip3\", isimip3.args = isimip3.range.args)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3.Bias adjustment of the daily temperature skewness:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# List of arguments that have to be passed to the \"biasCorrection\" function when using the ISIMIP3 method\n",
    "isimip3.skew.args  =  list(lower_bound =  c(0), lower_threshold =  c(0.0001), upper_bound =  c(1), \n",
    "                           upper_threshold =  c(0.9999), randomization_seed =  NULL,\n",
    "                           detrend =  array(data  =  FALSE, dim = 1), rotation_matrices =  c(NULL), \n",
    "                           n_quantiles = 50, distribution =  c(\"beta\"), \n",
    "                           trend_preservation = array(data = \"bounded\", dim = 1), \n",
    "                           adjust_p_values = array(data  =  FALSE, dim = 1), if_all_invalid_use  =  c(NULL),\n",
    "                           invalid_value_warnings  =  FALSE)\n",
    "\n",
    "# Adjusting historical simulations\n",
    "bc.skew.hist <- biasCorrection(y = y.skew, x = x.skew, newdata = x.skew, \"precipitation\" = FALSE, \n",
    "                          method = \"isimip3\", isimip3.args = isimip3.skew.args)\n",
    "# Adjusting future simulations\n",
    "bc.skew <- biasCorrection(y = y.skew, x = x.skew, newdata = newdata.skew, \"precipitation\" = FALSE, \n",
    "                          method = \"isimip3\", isimip3.args = isimip3.skew.args)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once these three steps are completed, the bias-adjusted daily maximum and minimum temperature can be obtained as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Bias-adjusted historical temperatures\n",
    "bc.tasmin.hist <- gridArithmetics(bc.tas.hist, gridArithmetics(bc.range.hist, bc.skew.hist, operator = \"*\"), \n",
    "                                  operator = \"-\")\n",
    "bc.tasmax.hist <- gridArithmetics(bc.tasmin.hist, bc.range.hist, operator = \"+\")\n",
    "\n",
    "# Bias-adjusted future temperatures\n",
    "bc.tasmin <- gridArithmetics(bc.tas, gridArithmetics(bc.range, bc.skew, operator = \"*\"), \n",
    "                             operator = \"-\")\n",
    "bc.tasmax <- gridArithmetics(bc.tasmin, bc.range, operator = \"+\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### Calculation of extreme temperature indices\n",
    "Finally, we use the `indexGrid` function to compute historical and future TX35 based on the bias-adjusted daily maximum temperature which has been obtained from the application of the ISIMIP3 method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Bias-adjusted TX35 for the historical period, year by year\n",
    "index.hist <- redim(indexGrid(tx = bc.tasmax.hist, index.code = \"TXth\", th = 35, \n",
    "                              time.resolution = \"year\"), drop = TRUE)\n",
    "tx35.hist <- climatology(index.hist)  # Mean value (number of days/year)\n",
    "\n",
    "# Bias-adjusted TX35 for the future period of interest, year by year\n",
    "index.rcp <- redim(indexGrid(tx = bc.tasmax, index.code = \"TXth\", th = 35, \n",
    "                             time.resolution = \"year\"), drop = TRUE)\n",
    "tx35.rcp <- climatology(index.rcp) # Mean value (number of days/year)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "spatialPlot(makeMultiGrid(tx35.hist, tx35.rcp, skip.temporal.check = TRUE), backdrop.theme = \"coastline\", \n",
    "            color.them = \"Reds\", names.attr = c(\"Historical (1986-2005)\",\"Projected (2041-2060, RCP8.5)\"),\n",
    "            main = \"Bias-adjusted TX35 (days/year)\", as.table = TRUE, set.min = 0, set.max = 60, at = seq(0,60,5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### Session info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    },
    "tags": [
     "hide-output"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "sessionInfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### References"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{bibliography}\n",
    ":filter: docname in docnames\n",
    "```\n",
    "\n",
    "<script src=\"https://utteranc.es/client.js\"\n",
    "        repo=\"PhantomAurelia/Atlas\"\n",
    "        issue-term=\"pathname\"\n",
    "        theme=\"preferred-color-scheme\"\n",
    "        crossorigin=\"anonymous\"\n",
    "        async>\n",
    "</script>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.7.12 | packaged by conda-forge | (default, Oct 26 2021, 06:08:53) \n[GCC 9.4.0]"
  },
  "vscode": {
   "interpreter": {
    "hash": "7812f8d43a0b4df1faa1f5f834d677598d68a57a2a8dffe43567739a16087a4c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}