2022_Bano_GMD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Downscaling Multi-Model Climate Projection Ensembles with Deep Learning (DeepESD): Contribution to CORDEX EUR-44\n",
    "### *Geoscientific Model Development*\n",
    "### J. Baño-Medina, R. Manzanas, E. Cimadevilla, J. Fernández, J. González-Abad, A.S. Cofiño, and J.M. Gutiérrez\n",
    "\n",
    "GitHub repository at https://github.com/SantanderMetGroup/DeepDownscaling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook reproduces the results presented in **Downscaling Multi-Model Climate Projection Ensembles with Deep Learning (DeepESD): Contribution to CORDEX EUR-44***, submitted to *Geoscientific Model Development* by *J. Baño-Medina, R. Manzanas, E. Cimadevilla, J. Fernández, J. González-Abad, A.S. Cofiño and J.M. Gutiérrez*. \n",
    "This paper presents *DeepESD*, the first dataset of high-resolution (0.5º) climate change projections (up to 2100) of daily precipitation and temperature over Europe obtained with deep learning techniques (in particular convolutional neural networks) from an ensemble of eight global climate models from the Coupled Model Intercomparison Project version 5 (CMIP5). \n",
    "\n",
    "**Note:** The technical specifications of the machine used to run the code presented herein can be found at the end of the notebook. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Preparing the R environment and working directories\n",
    "This notebook is written in the free programming language `R` (version 3.6.1) and builds on [`climate4R`](https://github.com/SantanderMetGroup/climate4R) (hereafter C4R), a suite of `R` packages developed by the [Santander Met Group](http://meteo.unican.es) for transparent climate data access, post processing (including bias correction and downscaling) and visualization. For details on climate4R (C4R hereafter), the interested reader is referred to [Iturbide et al. 2019](https://www.sciencedirect.com/science/article/pii/S1364815218303049?via%3Dihub).\n",
    "\n",
    "In particular, the following C4R libraries are used along the notebook: [`loadeR`](https://github.com/SantanderMetGroup/loadeR) and [`loadeR.2nc`](https://github.com/SantanderMetGroup/loadeR.2nc) (data loading), [`transformeR`](https://github.com/SantanderMetGroup/transformeR) (data manipulation), [`downscaleR`](https://github.com/SantanderMetGroup/downscaleR) and [`downscaleR.keras`](https://github.com/SantanderMetGroup/downscaleR.keras) (downscaling with neural networks) and [`visualizeR`](https://github.com/SantanderMetGroup/visualizeR) (visualization). To install them you may use the devtools package (e.g., ``devtools::install_github(\"SantanderMetGroup/downscaleR.keras@v1.0.0\")`` to install `downscaleR.keras`). Alternatively, you may directly install the entire C4R framework through ``conda``, C4R version 1.5.0, following the instructions provided at the end of this [page](https://github.com/SantanderMetGroup/climate4R). The latter option is highly recommended. Note that even with C4R v1.5.0 installed via conda, you still need to upgrade libraries `climate4R.UDG` and `VALUE` with the `devtools` package by typing: ``devtools::install_github(c(\"SantanderMetGroup/climate4R.UDG@v0.2.2\",\"SantanderMetGroup/VALUE@v2.2.2\"))``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(java.parameters = \"-Xmx8g\")  # expanding Java memory\n",
    "\n",
    "# C4R libraries\n",
    "library(loadeR) # version v1.7.0\n",
    "library(loadeR.2nc) # version v0.1.1\n",
    "library(transformeR) # version v2.1.0\n",
    "library(downscaleR) # version v3.3.2\n",
    "library(visualizeR) # version v1.6.0\n",
    "library(downscaleR.keras) # version v1.0.0 that build on Keras version 2.3.0 and tensorflow version 2.2.0 \n",
    "library(climate4R.value) # version v0.0.2 and relies on VALUE version v2.2.2\n",
    "library(climate4R.UDG) # version v0.2.2\n",
    "\n",
    "# Other useful libraries\n",
    "library(magrittr) # to operate with '%>%' or '%<>%'\n",
    "\n",
    "# For visualization purposes\n",
    "library(RColorBrewer)\n",
    "library(gridExtra)\n",
    "library(ggplot2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The predictions and models generated along the notebook are saved in a well-organized set of directories. Please use the `dir.create` function to create two new folders (*Data* and *models*) in your working directory. Within each of these folders, create subsequently two more subfolders, named *temperature* and *precip*. Finally, we create also the *figures* directory directly in the working directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Uncomment to create directories\n",
    "# dir.create(\"./Data\")\n",
    "# dir.create(\"./Data/precip\")\n",
    "# dir.create(\"./Data/temperature\")\n",
    "# dir.create(\"./models\")\n",
    "# dir.create(\"./models/precip\")\n",
    "# dir.create(\"./models/temperature\")\n",
    "# dir.create(\"./figures\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to load into our `R` environment all the data we are going to work with, which can be freely accessed through the [Climate Data Service](http://meteo.unican.es/cds) developed by the [Santander Met Group](http://meteo.unican.es) (non registered users need to register first [here](http://meteo.unican.es/udg-tap/signup)). Use the `loginUDG` function to log into the service with your own credentials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "loginUDG(username = \"youruser\", password = \"yourpassword\") # login into the Santander CDS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following Table lists the [Santander Climate Data Service (CDS)](http://www.meteo.unican.es/en/cds) endpoint of the datasets used in this study (except for E-OBS whose observational records can be found in their [website](https://www.ecad.eu/download/ensembles/download.php)). For a wide variety of datasets, C4R uses labels to point to these endpoints. The available labels can be displayed by typing `UDG.datasets()` into an R terminal. Throughout the notebook we lean on the corresponding labels to load the data into our R environment.\n",
    "\n",
    "| Dataset  | CDS endpoint: |\t\n",
    "|---|---|\n",
    "|ERA-Interim|https://data.meteo.unican.es/tds5/catalog/catalogs/interim/interim_DM_predictors.html?dataset=interim/daily/interim20_daily.ncml|\n",
    "|CMIP5|https://data.meteo.unican.es/tds5/catalog/catalogs/cmip5/cmip5Datasets.html|\n",
    "|CORDEX|https://data.meteo.unican.es/thredds/catalog/devel/c3s34d/catalog.html|\n",
    "|DeepESD|https://data.meteo.unican.es/thredds/catalog/esgcet/catalog.html|\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following block of code allows for loading the ERA-Interim predictor variables, which are needed to train our neural networks, for the period 1979-2005 by using the `loadGridData` function. Subsequently, the `makeMultiGrid` creates a unique C4R object containing all this information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predictor variables considered (see -*-)\n",
    "vars  <- c(\"psl\",\"z@500\",\"z@700\",\"z@850\", \n",
    "          \"hus@500\",\"hus@700\",\"hus@850\",\n",
    "          \"ta@500\",\"ta@700\",\"ta@850\",\n",
    "          \"ua@500\",\"ua@700\",\"ua@850\",\n",
    "          \"va@500\",\"va@700\",\"va@850\")\n",
    "# We loop over the variables and then use  makeMultiGrid, to bind the variables in a single C4R object\n",
    "x <- lapply(vars, function(z) {\n",
    "  loadGridData(dataset = \"ECMWF_ERA-Interim-ESD\",  # \"ECMWF_ERA-Interim-ESD\" is the label that identifies the dataset of interest in the Santander CDS\n",
    "               var = z,\n",
    "               lonLim = c(-8,34),  # domain of interest for the predictors\n",
    "               latLim = c(34,76),  # domain of interest for the predictors\n",
    "               years = 1979:2005)\n",
    "}) %>% makeMultiGrid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As predictands we use temperature and precipitation from E-OBS, which can be obtained as netCDF files [here](https://www.ecad.eu/download/ensembles/download.php). Once downloaded, these data can be imported in `R` with the `loadGridData` function. Subsequently, we upscale these E-OBS fields from their native 0.25º to the 0.5º regular grid our projections are delivered by using the `interpGrid` from `transformeR`. In the cell below, we illustrate how to load both precipitation and temperature with `loadGridData`, however, note that you should only load one at a time, or change the object name `y` to e.g., `y_rr` or `y_tg`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid05 = list(\"x\" = c(-9.75,30.25),\"y\" = c(34.25,74.25))  # boundaries of our projections' domain\n",
    "attr(grid05,\"resX\") <- attr(grid05,\"resY\") <- 0.5  # target spatial resolution for our projections\n",
    "\n",
    "## Please load only one predictand variable (either temperature or precipitation) \n",
    "## or give a different name to each variable (e.g. 'y_rr' and 'y_tg') ----------------------------------\n",
    "\n",
    "# To load E-OBS precipitation (previously downloaded as netCDF file)\n",
    "y <- loadGridData(dataset = \"rr_ens_mean_0.25deg_reg_v20.0e.nc\",\n",
    "                  var = \"rr\",\n",
    "                  lonLim = c(-10,30),\n",
    "                  latLim = c(34,74), \n",
    "                  years = 1979:2005) %>% interpGrid(new.coordinates = grid05, method = \"bilinear\")\n",
    "\n",
    "# To load E-OBS temperature (previously downloaded as netCDF file)\n",
    "y <- loadGridData(dataset = \"tg_ens_mean_0.25deg_reg_v20.0e.nc\",\n",
    "                  var = \"tg\",\n",
    "                  lonLim = c(-10,30),\n",
    "                  latLim = c(34,74), \n",
    "                  years = 1979:2005) %>% interpGrid(new.coordinates = grid05, method = \"bilinear\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We recommend the user to save the predictor `x` and predictand `y` data into `.rda` objects since these loading steps can be quite time-consuming."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Convolutional Neural Networks (CNNs)**\n",
    "\n",
    "To build *DeepESD* we rely on the convolutional neural networks (CNN) presented in [Baño-Medina et al. 2020](https://gmd.copernicus.org/articles/13/2109/2020/); in particular, on the CNN -*- models, which were found to provide robust results for precipitation (temperature) both in ''perfect-prognosis'' conditions but also in the GCM space. The cell below shows how to build these CNN models based on `Keras`, and save them in a custom function called `modelCNN`. Note that precipitation and temperature CNN models are different, so please (un)comment the needed lines depending on your particular target variable of interest. \n",
    "\n",
    "**Note:** We refer the reader to [Baño-Medina et al. 2020](https://gmd.copernicus.org/articles/13/2109/2020/) for further details about the exact configuration of the CNNs used herein."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Please select one ----------------------------------------------------------------------------------------------------------------\n",
    "## ----------------------------------------------------------------------------------------------------------------------------------\n",
    "## See https://gmd.copernicus.org/articles/13/2109/2020 for technical details\n",
    "## ----------------------------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "# Precipitation model\n",
    "modelCNN <- function(inp) {\n",
    "    # Input layer\n",
    "    inputs <- layer_input(shape = dim(inp$x.global)[2:4])\n",
    "    # Hidden layers\n",
    "    l1 = layer_conv_2d(inputs,filters = 50, kernel_size = c(3,3), activation = 'relu', padding = \"same\")\n",
    "    l2 = layer_conv_2d(l1,filters = 25, kernel_size = c(3,3), activation = 'relu', padding = \"same\")\n",
    "    l3 = layer_conv_2d(l2,filters = 1, kernel_size = c(3,3), activation = 'relu', padding = \"same\")\n",
    "    l4 = layer_flatten(l3)\n",
    "    # Output layer\n",
    "    l51 = layer_dense(l4,units = dim(inp$y$Data)[2], activation = 'sigmoid') \n",
    "    l52 = layer_dense(l4,units = dim(inp$y$Data)[2], activation = 'linear') \n",
    "    l53 = layer_dense(l4,units = dim(inp$y$Data)[2], activation = 'linear') \n",
    "    outputs <- layer_concatenate(list(l51,l52,l53))      \n",
    "    model <- keras_model(inputs = inputs, outputs = outputs) \n",
    "}\n",
    "\n",
    "## ----------------------------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "# Temperature model\n",
    "modelCNN <- function(inp) {\n",
    "    # Input layer    \n",
    "    inputs <- layer_input(shape = dim(inp$x.global)[2:4])\n",
    "    # Hidden layers    \n",
    "    l1 = layer_conv_2d(inputs,filters = 50, kernel_size = c(3,3), activation = 'relu', padding = \"valid\")\n",
    "    l2 = layer_conv_2d(l1,filters = 25, kernel_size = c(3,3), activation = 'relu', padding = \"valid\")\n",
    "    l3 = layer_conv_2d(l2,filters = 10, kernel_size = c(3,3), activation = 'relu', padding = \"valid\") \n",
    "    l4 = layer_flatten(l3)\n",
    "    # Output layer    \n",
    "    l51 = layer_dense(l4,units = dim(inp$y$Data)[2], activation = 'linear') \n",
    "    l52 = layer_dense(l4,units = dim(inp$y$Data)[2], activation = 'linear') \n",
    "    outputs <- layer_concatenate(list(l51,l52))      \n",
    "    model <- keras_model(inputs = inputs, outputs = outputs) \n",
    "}\n",
    "\n",
    "## ----------------------------------------------------------------------------------------------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. DeepESD\n",
    "In this section we 1) load the predictor variables of interest from the 8 GCM considered from the Santander CDS, 2) harmonize and standardize these predictor fields, 3) save these processed fields in `rda` objects to avoid repeating these steps in the future, 4) build the CNN models based on ERA-Interim predictors and E-OBS predictands and 5) apply these models to the GCM predictor variables to obtain the final high-resolution (downscaled at 0.5º) projections up to 2100.\n",
    "\n",
    "### 2.1 Preparing the predictor datasets\n",
    "The `dh` and `df` objects below contain the labels that identify the 8 GCMs considered in this work in the Santander CDS, for the historical and RCP.8.5 scenario, respectively. These lables are used when calling the `loadGridData` function for data loading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Use UDG.datasets() to obtain the labels of the desired GCMs\n",
    "\n",
    "# Historical scenario\n",
    "dh <- c(\"CMIP5-subset_CanESM2_r1i1p1_historical\",\n",
    "       \"CMIP5-subset_CNRM-CM5_r1i1p1_historical\",\n",
    "       \"CMIP5-subset_MPI-ESM-MR_r1i1p1_historical\",\n",
    "       \"CMIP5-subset_MPI-ESM-LR_r1i1p1_historical\",\n",
    "       \"CMIP5-subset_NorESM1-M_r1i1p1_historical\", \n",
    "       \"CMIP5-subset_GFDL-ESM2M_r1i1p1_historical\",\n",
    "       \"CMIP5-subset_EC-EARTH_r12i1p1_historical\",\n",
    "       \"CMIP5-subset_IPSL-CM5A-MR_r1i1p1_historical\")\n",
    "\n",
    "# RCP8.5 scenario\n",
    "df <- c(\"CMIP5-subset_CanESM2_r1i1p1_rcp85\",\n",
    "       \"CMIP5-subset_CNRM-CM5_r1i1p1_rcp85\", \n",
    "       \"CMIP5-subset_MPI-ESM-MR_r1i1p1_rcp85\",\n",
    "       \"CMIP5-subset_MPI-ESM-LR_r1i1p1_rcp85\",\n",
    "       \"CMIP5-subset_NorESM1-M_r1i1p1_rcp85\",\n",
    "       \"CMIP5-subset_GFDL-ESM2M_r1i1p1_rcp85\",\n",
    "       \"CMIP5-subset_EC-EARTH_r12i1p1_rcp85\",\n",
    "       \"CMIP5-subset_IPSL-CM5A-MR_r1i1p1_rcp85\")\n",
    "\n",
    "# Labels used to identify the 8 GCMs\n",
    "dName <- c(\"CanESM2\",\"CNRM-CM5\",\"MPI-ESM-MR\",\"MPI-ESM-LR\",\"NorESM1\",\"GFDL\",\"EC-Earth\",\"IPSL\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following loop allows us to load the predictors from the above GCMs over our target domain for the reference (1975-2005) and future (early-future: 2006-2040, mid-future: 2041-2070, far-future: 2071-2100) periods of interest. Note that the historical (RCP8.5) scenario is used for the reference (future) periods. Note also that, within the loop, all the GCMs are interpolated to the spatial resolution of the ERA-Interim predictors which were used to fit the CNNs. Once loaded, the GCM predictors are saved as `.rda` files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lapply(c(\"h\",\"ef\",\"mf\",\"ff\"), FUN = function(sc) {\n",
    "    if (sc == \"h\")  years <- 1975:2005  # historical period of interest\n",
    "    if (sc == \"ef\") years <- 2006:2040  # early-future\n",
    "    if (sc == \"mf\") years <- 2041:2070  # mid-future\n",
    "    if (sc == \"ff\") years <- 2071:2100  # far-future\n",
    "    if (sc == \"h\"){d <- dh} else {d <- df}\n",
    "    \n",
    "    # We loop over the GCMs\n",
    "    lapply(1:length(d), FUN = function(zz) {    \n",
    "    x <- lapply(vars, function(z) loadGridData(dataset = d[zz],var = z,\n",
    "                                               lonLim = c(-8,34),latLim = c(34,76),\n",
    "                                               years = years) %>% interpGrid(new.coordinates = getGrid(x))) %>% makeMultiGrid()      \n",
    "    \n",
    "    # Since the IPSL contains NA values in the 850hPa level at certain gridpoints, we replace these NA values with the numeric values of their closest neighbours.\n",
    "    if (dName[zz] == \"IPSL\") {\n",
    "      ind850 <- grepl(\"@850\",x$Variable$varName,fixed = TRUE) %>% which()\n",
    "      indGP <- apply(x$Data[ind850[1],1,,,], MARGIN = c(2,3), anyNA) %>% which(arr.ind = TRUE)\n",
    "      for (i in 1:nrow(indGP)) {\n",
    "        indTime <- is.na(x$Data[ind850[1],1,,indGP[i,1],indGP[i,2]]) %>% which()\n",
    "        x$Data[ind850,1,indTime,indGP[i,1],indGP[i,2]] <- x$Data[ind850,1,indTime,indGP[i,1],indGP[i,2]-1]\n",
    "      }\n",
    "    }\n",
    "                \n",
    "    # We save the predictor fields as .rda files            \n",
    "    if (sc == \"h\") {\n",
    "        xh <- x\n",
    "        save(xh, file = paste0(\"./Data/xh_\",dName[zz],\".rda\"))\n",
    "        rm(x,xh)\n",
    "    } else {\n",
    "        xf <- x\n",
    "        save(xf, file = paste0(\"./Data/x\",sc,\"_\",dName[zz],\".rda\"))\n",
    "        rm(x,xf)\n",
    "    }            \n",
    "  })  \n",
    "})  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following loop allows us to harmonize and standardize the GCM predictors loaded in the previous step to assure they reasonable resemble the ERA-Interim variables used to train the CNN models (note this is one of the key assumptions that are done in ''perfect-prognosis'' downscaling). For this harmonization+standardization step, which is different depending on the particular scenario of interest (the reader is referred again to [Baño-Medina et al. 2020](https://gmd.copernicus.org/articles/13/2109/2020/) for details about this process), the `scaleGrid` function from `transformeR` is used. The so-processed GCM predictors, wich will be used as inputs to the CNN models, are saved as `rda` files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We loop over the GCMs\n",
    "lapply(1:length(dName), FUN = function(zz) {\n",
    "  load(paste0(\"./Data/xh_\",dName[zz],\".rda\")) \n",
    "  # We loop over the temporal periods  \n",
    "  lapply(c(\"h\",\"ef\",\"mf\",\"ff\"), FUN = function(sc) {\n",
    "      \n",
    "    # We harmonize and standardize the historical scenario  \n",
    "    if (sc == \"h\") {\n",
    "      xh <- scaleGrid(xh, base = subsetGrid(xh, years = 1979:2005), ref = x, type = \"standardize\", spatial.frame = \"gridbox\", time.frame = \"monthly\") # harmonization      \n",
    "      xn <- scaleGrid(xh, base = subsetGrid(xh, years = 1979:2005), type = \"standardize\")  # standardization\n",
    "        \n",
    "    # We harmonize and standardize the RCP8.5 scenario    \n",
    "    } else {\n",
    "      load(paste0(\"./Data/x\",sc,\"_\",dName[zz],\".rda\"))  \n",
    "      xf <- scaleGrid(xf, base = subsetGrid(xh, years = 1979:2005), ref = x, type = \"standardize\", spatial.frame = \"gridbox\", time.frame = \"monthly\") # harmonization    \n",
    "      xh <- scaleGrid(xh, base = subsetGrid(xh, years = 1979:2005), ref = x, type = \"standardize\", spatial.frame = \"gridbox\", time.frame = \"monthly\") # harmonization     \n",
    "      xn <- xf %>% scaleGrid(base = subsetGrid(xh, years = 1979:2005), type = \"standardize\")  # standardization                     \n",
    "    }\n",
    "      \n",
    "    # We save the standardized predictor fields as `rda` objects  \n",
    "    save(xn, file = paste0(\"./Data/xn\",sc,\"_\",dName[zz],\".rda\"))  \n",
    "  })\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Precipitation downscaling\n",
    "This section shows how to fit the CNN model which links the large-scale predictors from ERA-Interim with the high-resolution E-OBS precipitation at surface. The steps to take would be the following:\n",
    "- Prepare the predictor and predictand tensors with the `prepareData.keras` function from `downscaleR.keras`.\n",
    "- Standardize the ERA-Interim predictors with the `scaleGrid` function from `transformeR`. \n",
    "- For a better fit of the Gamma distribution, 0.99 is substracted from observed precipitation and negative values are ignored (note that this step implies that rainy days are defined as those receiving 1 or more mm of precipitation). To do this, the `gridArithmetics` and `binaryGrid` functions from `transformeR` are used.\n",
    "- Train the CNN model encapsuled in the `modelCNN` (which has been previously defined) with the `downscaleTrain.keras` function from `downscaleR.keras`. To optimize the negative log-likelihood of the Bernoulli-Gamma distribution, we employ the custom loss function `bernouilliGammaLoss` from `downscaleR.keras`. The network is fitted using the adam optimizer and a learning rate of 1e-4. Early-stopping with a patience of 30 epochs is applied and the best model (epoch) is saved in the working directory as a `.h5` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: Running this cell takes about 1 hour \n",
    "\n",
    "# Preparing predictor and predictand data for downscaling with downscaleR.keras\n",
    "xyT <- prepareData.keras(x = scaleGrid(x,type = \"standardize\"), y = binaryGrid(gridArithmetics(y,0.99,operator = \"-\"),condition = \"GE\",threshold = 0,partial = TRUE),\n",
    "                         first.connection = \"conv\",last.connection = \"dense\",channels = \"last\")  \n",
    "\n",
    "# Training the CNN model to downscale precipitation\n",
    "downscaleTrain.keras(obj = xyT,model = modelCNN(xyT),clear.session = TRUE,\n",
    "                     compile.args = list(\"loss\" = bernouilliGammaLoss(last.connection = \"dense\"),\"optimizer\" = optimizer_adam(lr = 0.0001)),\n",
    "                     fit.args = list(\"batch_size\" = 100,\"epochs\" = 10000,\"validation_split\" = 0.1,\"verbose\" = 1,\n",
    "                                     \"callbacks\" = list(callback_early_stopping(patience = 30),callback_model_checkpoint(filepath='./models/precip/CNN1.h5',monitor='val_loss', save_best_only=TRUE))))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the model is trained, we use it to predict in both the train (training period using ERA-Interim variables) and the GCM spaces. As per the former, we are interested in the estimation of the parameter `p` (probability of rain), since it is needed to later adjust the frequency of rain in the high-resolution projections obtained from the GCM (see the manuscript for details). To compute `p` in the train period the following is done:\n",
    "- Prepare the predictors which will serve as inputs for the CNN model with the `prepareNewData.keras` function. Subsequently, use them to predict in the train set with the `downscalePredict.keras` function. The `model` argument indicates the path where the CNN model was previously stored, and `C4Rtemplate` is a C4R object used as template for the predictions which provides de proper metadata. Since `downscalePredict.keras` outputs 3 parameters (the probability of rain, `p`, and the logarithmic of the shape and scale parameters of the Gamma distribution, `log_alpha` and `log_beta`), the `subsetGrid` is applied in order to keep only `p`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preparing predictor data\n",
    "xyt <- prepareNewData.keras(scaleGrid(x,type = \"standardize\"), xyT) \n",
    "pred_ocu_train <- downscalePredict.keras(newdata = xyt,C4R.template = y,clear.session = TRUE,loss = \"bernouilliGammaLoss\",\n",
    "                                         model = list(\"filepath\" = './models/precip/CNN1.h5',\"custom_objects\" = c(\"custom_loss\" = bernouilliGammaLoss(last.connection = \"dense\")))) %>% subsetGrid(var = \"p\")  \n",
    "rm(xyt) # to save memory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, the trained CNN model is used to generate the high-resolution projections building on the 8 GCMs considered in this work. To do so, we perform a loop over the distinct GCMs in which the corresponding predictors (which had been previously saved) are loaded and conveniently transformed using the `prepareNewData.keras` function. Finally, the `log_alpha`, `log_beta` and `p` parameters, which are obtained with the `downscalePredict.keras` function are saved in the `pred` object.\n",
    "- On the one hand, `log_alpha` and `log_beta` are used to obtain the rainfall amount with the `computeRainfall` function from `downscaleR.keras`. The argument `simulate` allows us for specifying if either a stochastic or a deterministic outcome is wanted. The argument `bias` is used to re-center the Gamma distribution to 1mm/day. \n",
    "- On the other hand, we use the `p` parameter to derive the binary event occurrence/non occurrence through the `bynaryGrid` function.\n",
    "- Finally, both series (binary and continuous) are multiplied to produce the complete precipitation time series.\n",
    "\n",
    "In the following block of code we compute the (deterministic) frequency and (stochastic) amount of rain. The generated projections are saved in `.nc` format with the `grid2nc` function.\n",
    "\n",
    "***Note:*** If a purely deterministic, or a purely stochastic version of the projections is wanted, please uncomment the corresponding lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: Running this cell takes about 4 hours\n",
    "\n",
    "# We loop over the GCMs\n",
    "lapply(1:length(dName), FUN = function(zz) {\n",
    "    \n",
    "  # We loop over the temporal periods  \n",
    "  lapply(c(\"ef\",\"mf\"), FUN = function(sc) {\n",
    "    load(paste0(\"./Data/xn2\",sc,\"_\",dName[zz],\".rda\"))\n",
    "    xyt <- prepareNewData.keras(xn,xyT)  \n",
    "    pred <- downscalePredict.keras(newdata = xyt,C4R.template = y,clear.session = TRUE,loss = \"bernouilliGammaLoss\",\n",
    "                                   model = list(\"filepath\" = './models/precip/CNN1.h5',\"custom_objects\" = c(\"custom_loss\" = bernouilliGammaLoss(last.connection = \"dense\")))) \n",
    "    \n",
    "    ## Frequency (deterministic) and amount of rain (stochastic) ------------------------------------------\n",
    "    pred2 <- computeRainfall(log_alpha = subsetGrid(pred,var = \"log_alpha\"),log_beta = subsetGrid(pred,var = \"log_beta\"),bias = 1,simulate = TRUE) %>% gridArithmetics(binaryGrid(subsetGrid(pred,var = \"p\"),ref.obs = binaryGrid(y,threshold = 1, condition = \"GE\"),ref.pred = pred_ocu_train))  \n",
    "    grid2nc(pred2,NetCDFOutFile = paste0(\"./Data/precip/CNN2_\",sc,\"_\",dName[zz],\".nc4\"))   \n",
    "  })    \n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Temperature downscaling\n",
    "This section shows how to train the CNN model which links the large-scale predictors from ERA-Interim with the high-resolution E-OBS temperature at surface. As for precipitation, the steps to take would be the following:\n",
    "\n",
    "- Prepare the predictor and predictand tensors with the `prepareData.keras` function from `downscaleR.keras`.\n",
    "- Standardize the ERA-Interim predictors with the `scaleGrid` function from `transformeR`. \n",
    "- Train the CNN model encapsuled in the `modelCNN` (which has been previously defined) with the `downscaleTrain.keras` function from `downscaleR.keras`. To optimize the negative log-likelihood of the Gaussian distribution, we use the custom loss function `GaussianLoss` from `downscaleR.keras`. Whe network is fitted using the adam optimizer and a learning rate of 1e-4. Early-stopping with a patience of 30 epochs is applied and the best model (epoch) is saved in the working directory as a `.h5` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: Running this cell takes about 1 hour\n",
    "\n",
    "# Preparing predictor and predictand data for downscaling with downscaleR.keras\n",
    "xyT <- prepareData.keras(x = scaleGrid(x,type = \"standardize\"),y = y,\n",
    "                         first.connection = \"conv\",last.connection = \"dense\",channels = \"last\")  \n",
    "\n",
    "# Training the CNN model to downscale temperature\n",
    "downscaleTrain.keras(obj = xyT,model = modelCNN(xyT),clear.session = TRUE,\n",
    "                     compile.args = list(\"loss\" = gaussianLoss(last.connection = \"dense\"),\"optimizer\" = optimizer_adam(lr = 0.0001)),\n",
    "                     fit.args = list(\"batch_size\" = 100,\"epochs\" = 10000,\"validation_split\" = 0.1,\"verbose\" = 1,\n",
    "                                     \"callbacks\" = list(callback_early_stopping(patience = 30),callback_model_checkpoint(filepath='./models/temperature/CNN10.h5',monitor='val_loss', save_best_only=TRUE))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once trained, the CNN model is used to generate the high-resolution projections of temperature. Similarly as for precipitation, we apply for this task the `prepareNewData.keras` and `downscalePredict.keras` functions, saving the `mean` and `log_var` parameters in the `pred` object. The code below allows for obtaining deterministic projections (note that only the `mean` parameter is needed in this case) and saving them in `.nc` format with the `grid2nc` function.\n",
    "\n",
    "***Note:*** If stochastic projections are wanted, please uncomment the corresponding lines.\n",
    "he projections as `.nc`using function `grid2nc`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: Running this cell takes about 1 hour and 20 minutes\n",
    "\n",
    "lapply(1:length(dName), FUN = function(zz) {\n",
    "  lapply(c(\"h\",\"ef\",\"mf\",\"ff\"), FUN = function(sc) {\n",
    "    load(paste0(\"./Data/xn\",sc,\"_\",dName[zz],\".rda\"))\n",
    "    xyt <- prepareNewData.keras(xn,xyT)\n",
    "    pred <- downscalePredict.keras(newdata = xyt,C4R.template = y,clear.session = TRUE,loss = \"gaussianLoss\",\n",
    "                                   model = list(\"filepath\" = './models/temperature/CNN10.h5',\"custom_objects\" = c(\"custom_loss\" = gaussianLoss(last.connection = \"dense\")))) \n",
    "    \n",
    "    ## Deterministic version \n",
    "    pred <- subsetGrid(pred, var = \"mean\")\n",
    "    pred$Variable$varName <- \"tas\"\n",
    "    grid2nc(pred,NetCDFOutFile = paste0(\"./Data/temperature/CNN_\",sc,\"_\",dName[zz],\".nc4\"))               \n",
    "      \n",
    "    k_clear_session()  \n",
    "  })    \n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Dynamical climate models\n",
    "To assess the credibility of DeepESD, it is compared against two different ensembles of dynamical models (see the *Technical Validation* section), the first/second of them formed by Global/Regional Climate Models (GCMs/RCMs). Since DeepESD covers only land, we start by creating a 0.5º land-sea mask which will be later applied to eliminate sea points from both GCMs and RCMs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask <- gridArithmetics(subsetGrid(y,year = 1990),0) %>% gridArithmetics(1,operator = \"+\") %>% climatology()\n",
    "grid2nc(mask, NetCDFOutFile = \"./Data/mask.nc4\")  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Ensemble of Global Climate Models (GCMs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We perform a loop over the temporal periods of interest (1975-2005 for the historical scenario plus 2006-2040, 2041-2070 and 2071-2100 for RCP8.5) and save the GCM ensemble as netCDF files (`grid2nc` function) in a multi-member C4R object. All GCMs are interpolated to our target 0.5º resolution (E-OBS grid), using conservative remapping. To do this interpolation, we rely on the `cdo` library and use function `system` to invoke the OS command. Please note that you can install the `cdo` library with conda, by typing `conda install cdo` in a terminal. Finally, sea points are removed by applying the land-sea mask we have previously created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Please select one:   -----------------------------------------------------------------------\n",
    "# Parameter Setting for precipitation\n",
    "variable <- \"precip\"\n",
    "var <- \"pr\"\n",
    "# Parameter Setting for temperature\n",
    "variable <- \"temperature\"\n",
    "var <- \"tas\"\n",
    "## --------------------------------------------------------------------------------------------\n",
    "\n",
    "# We loop over the temporal periods of interest\n",
    "lapply(c(\"h\",\"ef\",\"mf\",\"ff\"), FUN = function(z) {\n",
    "  if (z == \"h\")  years <- 1975:2005  # historical\n",
    "  if (z == \"ef\") years <- 2006:2040  # RCP8.5\n",
    "  if (z == \"mf\") years <- 2041:2070  # RCP8.5\n",
    "  if (z == \"ff\") years <- 2071:2100  # RCP8.5\n",
    "  if (z == \"h\") {d <- dh} else {d <- df} \n",
    "    \n",
    "  # We loop over the GCM labels  \n",
    "  lapply(1:length(d), FUN = function(zzz) {    \n",
    "    # Load the data and interpolate to the target resolution with interpGrid  \n",
    "    yy <- loadGridData(dataset = d[zzz],var = var,years = years,lonLim = c(-10,30),latLim = c(34,74))\n",
    "    grid2nc(yy, NetCDFOutFile = \"./aux.nc4\") # we save the GCM cropped to the European domain, and save as .nc in an auxiliary variable \n",
    "    system(paste0(\"cdo remapcon,\", \"./Data/mask.nc4\", \" \", \"./aux.nc4\", \" \", \"./aux2.nc4\"))  # We use system function to call the cdo library and interpolate the grid using conservative remapping\n",
    "    yy <- loadGridData(\"./aux2.nc4\", var = var) # we load the interpolated GCM field\n",
    "    file.remove(c(\"./aux.nc4\",\"./aux2.nc4\"))  \n",
    "    # Apply land-sea mask \n",
    "    yy <- lapply(1:getShape(yy,\"time\"), FUN = function(z) gridArithmetics(subsetDimension(yy, dimension = \"time\", indices = z),mask)) %>% bindGrid(dimension = \"time\")            \n",
    "    # Save the GCM members    \n",
    "    grid2nc(yy,NetCDFOutFile = paste0(\"./Data/\",variable,\"/y_\",z,\"_\",dName[zzz],\".nc4\"))\n",
    "  }) \n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Ensemble of Regional Climate Models (RCMs) \n",
    "In this section we form an ensemble of EURO-CORDEX RCMs which can be easily loaded from the Santander CDS by using the appropiate labels (see the block below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Labels for the historical scenario\n",
    "dh <- c(\"CORDEX-EUR-44_CCCma-CanESM2_historical_r1i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_CNRM-CERFACS-CNRM-CM5_historical_r1i1p1_ETH-CLMcom-CCLM5-0-6_v1\",\n",
    "        \"CORDEX-EUR-44_CNRM-CERFACS-CNRM-CM5_historical_r1i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_MPI-M-MPI-ESM-LR_historical_r1i1p1_CLMcom-CCLM4-8-17_v1\",\n",
    "        \"CORDEX-EUR-44_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1\",\n",
    "        \"CORDEX-EUR-44_NCC-NorESM1-M_historical_r1i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_NOAA-GFDL-GFDL-ESM2M_historical_r1i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_ICHEC-EC-EARTH_historical_r12i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_ICHEC-EC-EARTH_historical_r12i1p1_ETH-CLMcom-CCLM5-0-6_v1\",\n",
    "        \"CORDEX-EUR-44_IPSL-IPSL-CM5A-MR_historical_r1i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_IPSL-IPSL-CM5A-MR_historical_r1i1p1_IPSL-INERIS-WRF331F_v1\")\n",
    "# Labels for the RCP8.5 scenario\n",
    "df <- c(\"CORDEX-EUR-44_CCCma-CanESM2_rcp85_r1i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_CNRM-CERFACS-CNRM-CM5_rcp85_r1i1p1_ETH-CLMcom-CCLM5-0-6_v1\",\n",
    "        \"CORDEX-EUR-44_CNRM-CERFACS-CNRM-CM5_rcp85_r1i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_MPI-M-MPI-ESM-LR_rcp85_r1i1p1_CLMcom-CCLM4-8-17_v1\",\n",
    "        \"CORDEX-EUR-44_MPI-M-MPI-ESM-LR_rcp85_r1i1p1_MPI-CSC-REMO2009_v1\",\n",
    "        \"CORDEX-EUR-44_NCC-NorESM1-M_rcp85_r1i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_NOAA-GFDL-GFDL-ESM2M_rcp85_r1i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_ICHEC-EC-EARTH_rcp85_r12i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_ICHEC-EC-EARTH_rcp85_r12i1p1_ETH-CLMcom-CCLM5-0-6_v1\",\n",
    "        \"CORDEX-EUR-44_IPSL-IPSL-CM5A-MR_rcp85_r1i1p1_SMHI-RCA4_v1\",\n",
    "        \"CORDEX-EUR-44_IPSL-IPSL-CM5A-MR_rcp85_r1i1p1_IPSL-INERIS-WRF331F_v1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We perfrom a loop over the temporal periods of interest (1975-2005 for the historical scenario plus 2006-2040, 2041-2070 and 2071-2100 for RCP8.5) and save the RCM ensemble as netCDF files (`grid2nc` function) in a multi-member C4R object. All RCMs are interpolated to our target 0.5º resolution (E-OBS grid) and sea points are removed by applying the land-sea mask we have previously created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Please select one:   -----------------------------------------------------------------------\n",
    "# Precipitation\n",
    "variable <- \"precip\"\n",
    "var <- \"pr\"\n",
    "## Temperature\n",
    "variable <- \"temperature\"\n",
    "var <- \"tas\"\n",
    "## --------------------------------------------------------------------------------------------\n",
    "\n",
    "# We loop over the temporal periods of interest\n",
    "lapply(c(\"h\",\"ef\",\"mf\",\"ff\"), FUN = function(z) {\n",
    "  if (z == \"h\")  years <- 1975:2005  # historical\n",
    "  if (z == \"ef\") years <- 2006:2040  # RCP8.5\n",
    "  if (z == \"mf\") years <- 2041:2070  # RCP8.5\n",
    "  if (z == \"ff\") years <- 2071:2099  # RCP8.5\n",
    "  if (z == \"h\") {d <- dh} else {d <- df}\n",
    "    \n",
    "  # We loop over the RCM labels  \n",
    "  yy <- lapply(1:length(d), FUN = function(zzz) { \n",
    "        # Load the data and interpolate to the target resolution with interpGrid\n",
    "        yy <- loadGridData(dataset = d[zzz],var = var,years = years,lonLim = c(-10,30),latLim = c(34,74)) %>% interpGrid(getGrid(y))\n",
    "        # Apply land-sea mask\n",
    "        yy <- lapply(1:getShape(yy,\"time\"), FUN = function(z) gridArithmetics(subsetDimension(yy, dimension = \"time\", indices = z),mask)) %>% bindGrid(dimension = \"time\")  \n",
    "        # Save the RCM members\n",
    "        grid2nc(yy,NetCDFOutFile = paste0(\"./Data/\",variable,\"/yRCM_\",z,\"_member\",zzz,\".nc4\"))\n",
    "  })\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Results\n",
    "This section provides the code needed to reproduce the figures presented in the manuscript. Note we mostly rely on the `visualizeR` package for plotting since it supports both spatial maps and temporal series.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Ensemble mean and bias with respect to E-OBS\n",
    "Figure 1 in the manuscript shows the climatology of the different ensembles built (GCM, RCM and DeepESD), along with the corresponding mean error (bias) with respect to the observed pattern in the historical period. We start thus by computing the climatology of the different contributing members forming each ensemble and saving them as netCDF files in the working directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Please select one:   -----------------------------------------------------------------------\n",
    "# Precipitation\n",
    "variable <- \"precip\"\n",
    "var <- \"pr\"\n",
    "# Temperature\n",
    "variable <- \"temperature\"\n",
    "var <- \"tas\"\n",
    "## --------------------------------------------------------------------------------------------\n",
    "\n",
    "dates <- list(start = \"2006-01-01 12:00:00 GMT\", end = \"2041-01-01 12:00:00 GMT\") # because a member of the RCM ensemble misses a value on the 01-Jan-2006, so to preserve temporal consistency in the ensemble we add this date to the ensemble mean metadata\n",
    "dName1 <- c(\"CanESM2\",\"CNRM-CM5\",\"MPI-ESM-MR\",\"MPI-ESM-LR\",\"NorESM1\",\"GFDL\",\"EC-Earth\",\"IPSL\")\n",
    "dName2 <- 1:11\n",
    "nn <- c(\"GCM\",\"RCM\",\"CNN\")\n",
    "lapply(1:length(nn), FUN = function(zz) {\n",
    "  if (nn[zz] == \"RCM\") {dName <- dName2} else {dName <- dName1}  \n",
    "  lapply(c(\"h\",\"ef\",\"mf\",\"ff\"), FUN = function(z) {  \n",
    "    pred <- lapply(1:length(dName), FUN = function(zzz) { \n",
    "      if (zz == 1) path <- paste0(\"./Data/\",variable,\"/y_\",z,\"_\",dName[zzz],\".nc4\") # GCM ensemble\n",
    "      if (zz == 2) path <- paste0(\"./Data/\",variable,\"/yRCM_\",z,\"_member\",zzz,\".nc4\") # RCM ensemble           \n",
    "      if (zz == 3) path <- paste0(\"./Data/\",variable,\"/CNN_\",z,\"_\",dName[zzz],\".nc4\") # DeepESD ensemble   \n",
    "      grid <- loadGridData(dataset = path,var = var)\n",
    "      grid <- valueIndex(grid, index.code = \"Mean\")$Index  # computing the mean of each member \n",
    "      if (z == \"ef\") grid$Dates <- dates \n",
    "      return(grid)  \n",
    "    }) %>% bindGrid(dimension = \"member\")  # bind the member means in a single C4R object along the `member` dimension\n",
    "    \n",
    "    pred$InitializationDates <- NULL\n",
    "    # Saving the ensemble mean in netCDF format  \n",
    "    grid2nc(pred,NetCDFOutFile = paste0(\"./Data/\",variable,\"/\",nn[zz],\"_\",z,\"_ensemble.nc4\"))\n",
    "    pred <- NULL  \n",
    "  })\n",
    "})  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the mean climatology for each ensemble is obtained with the `aggregateGrid` function (note the aggregation is done along the member dimension) from `transformeR`. Afterwards, we can already use `spatialPlot` to plot the corresponding spatial pattern. The resulting figure is saved in `pdf` format in the path indicated in the `pdfOutput` object.\n",
    "\n",
    "***Note:*** Depending on the target variable of interest the user should comment/uncomment the appropriate lines at the beginning of the cell, which define the plotting parameters better suited for precipitation and temperature. This applies to the rest of the notebook from now on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Please select one:   -----------------------------------------------------------------------\n",
    "# Precipitation\n",
    "cb <- c(\"#FFFFFF\",brewer.pal(n = 9, \"BuPu\"))\n",
    "cb <- cb %>% colorRampPalette()\n",
    "at <- seq(0,8,0.5)    \n",
    "units <- \"mm/day\"\n",
    "pdfOutput <- \"./figures/ensembleMean_pr.pdf\" \n",
    "var <- \"precip\"\n",
    "# Temperature\n",
    "cb <- c(\"#FFFFFF\",brewer.pal(n = 9, \"OrRd\"))\n",
    "cb <- cb %>% colorRampPalette()\n",
    "at <- seq(-5, 20,2.5) \n",
    "units <- \"ºC\"\n",
    "pdfOutput <- \"./figures/ensembleMean_tas.pdf\" \n",
    "var <- \"temperature\"\n",
    "## --------------------------------------------------------------------------------------------\n",
    "\n",
    "nn <- c(\"GCM\",\"RCM\",\"CNN\")\n",
    "figs <- lapply(1:length(nn), FUN = function(z) {\n",
    "  # We store in `grid` object the ensemble of climatologies  \n",
    "  grid <- loadGridData(paste0(\"./Data/\",var,\"/\",nn[z],\"_h_ensemble.nc4\"), var = \"Mean\")\n",
    "  # Compute the ensemble mean  \n",
    "  gridMean <-  aggregateGrid(grid,aggr.mem = list(FUN = \"mean\", na.rm = TRUE)) \n",
    "  # We depict the ensemble mean with spatialPlot function  \n",
    "  spatialPlot(gridMean,\n",
    "              backdrop.theme = \"coastline\",\n",
    "              main = paste0(\"Ensemble Mean (\",units,\") - \",nn[z]),\n",
    "              ylab = \"1975-2005\",\n",
    "              col.regions = cb,\n",
    "              at = at,\n",
    "              set.min = at[1], set.max = at[length(at)])\n",
    "}) \n",
    "pdf(pdfOutput, width = 15, height = 10)   \n",
    "grid.arrange(grobs = figs, ncol = 3)                     \n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot the bias with respecto to the observed (i.e. E-OBS) climatology. Again, we rely on `spatialPlot` to depict the spatial fields, and `aggregateGrid` and `gridArithmetics`, to compute the ensemble mean and its bias, respectively. The resulting figures are saved in `pdf` format in the path indicated by `pdfOutput`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Please select one:   -----------------------------------------------------------------------\n",
    "# Precipitation\n",
    "cb <- brewer.pal(n = 11, \"BrBG\")\n",
    "cb[6] <- \"#FFFFFF\"; cb <- cb %>% colorRampPalette()\n",
    "at <- c(seq(-2, -0.5,0.5),-0.25,0.25,seq(0.5, 2,0.5))    \n",
    "units <- \"mm/day\"\n",
    "pdfOutput <- \"./figures/bias_pr.pdf\" \n",
    "var <- \"precip\"\n",
    "# Temperature\n",
    "cb <- rev(brewer.pal(n = 11, \"RdBu\"))\n",
    "cb[6] <- \"#FFFFFF\"; cb <- cb %>% colorRampPalette()\n",
    "at <- c(seq(-2, -0.5,0.5),-0.25,0.25, seq(0.5,2,0.5)) \n",
    "units <- \"ºC\"\n",
    "pdfOutput <- \"./figures/bias_tas.pdf\"\n",
    "var <- \"temperature\"\n",
    "## --------------------------------------------------------------------------------------------\n",
    "\n",
    "nn <- c(\"GCM\",\"RCM\",\"CNN\")\n",
    "figs <- lapply(1:length(nn), FUN = function(z) {\n",
    "  # Compute the ensemble mean  \n",
    "  grid <- loadGridData(paste0(\"./Data/\",var,\"/\",nn[z],\"_h_ensemble.nc4\"), var = \"Mean\") %>% aggregateGrid(aggr.mem = list(FUN = \"mean\", na.rm = TRUE))\n",
    "  # Compute the bias with respect to the observed temporal climatology for the same period\n",
    "  grid  %<>% gridArithmetics(valueIndex(y, index.code = \"Mean\")$Index,operator = \"-\")\n",
    "  # Depict the bias of the ensemble mean\n",
    "  spatialPlot(grid,\n",
    "              backdrop.theme = \"coastline\",\n",
    "              main = paste0(\"Bias Ensemble Mean (\",units,\") - \",nn[z]),\n",
    "              ylab = \"1975-2005\",\n",
    "              col.regions = cb,\n",
    "              at = at,\n",
    "              set.min = at[1], set.max = at[length(at)])  \n",
    "}) \n",
    "pdf(pdfOutput, width = 15, height = 10)   \n",
    "grid.arrange(grobs = figs, ncol = 3)                     \n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Climate change signals\n",
    "To produce Figure 2 in the manuscript we perform a loop, for each ensemble (GCM, RCM and DeepESD), over the different RCP8.5 periods of interest and sequentially compute the difference between the future climatology and the historical one. These climate change signals are then averaged along the member dimension using the `aggregateGrid` function. The resulting figures are saved in `pdf` format in the path indicated by `pdfOutput`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Please select one:   -----------------------------------------------------------------------\n",
    "# Precipitation\n",
    "cb <- brewer.pal(n = 11, \"BrBG\")\n",
    "cb[6] <- \"#FFFFFF\"; cb <- cb %>% colorRampPalette()\n",
    "at <- c(seq(-1, -0.25,0.25),-0.125,0.125,seq(0.25, 1,0.25)) \n",
    "pdfOutput <- \"./figures/deltas_pr.pdf\" \n",
    "var <- \"precip\"\n",
    "# Temperature\n",
    "cb <- c(\"#FFFFFF\",brewer.pal(n = 11, \"OrRd\"))\n",
    "at <- c(seq(0,4,0.5),5,6)\n",
    "pdfOutput <- \"./figures/deltas_tas.pdf\"\n",
    "var <- \"temperature\"\n",
    "## --------------------------------------------------------------------------------------------\n",
    "\n",
    "nn <- c(\"GCM\",\"RCM\",\"CNN\")\n",
    "figs <- lapply(c(\"ef\",\"mf\",\"ff\"), FUN = function(z) {  \n",
    "    lapply(1:length(nn), FUN = function(zz) {\n",
    "    gridh <- loadGridData(paste0(\"./Data/\",var,\"/\",nn[zz],\"_h_ensemble.nc4\"), var = \"Mean\") \n",
    "    gridf <- loadGridData(paste0(\"./Data/\",var,\"/\",nn[zz],\"_\",z,\"_ensemble.nc4\"), var = \"Mean\") \n",
    "    grid <- gridArithmetics(gridf,gridh,operator = \"-\")\n",
    "    gridMean <-  aggregateGrid(grid,aggr.mem = list(FUN = \"mean\", na.rm = TRUE)) \n",
    "\n",
    "    if (z == \"ef\") period <- c(\"2006-2040\")\n",
    "    if (z == \"mf\") period <- c(\"2041-2070\")\n",
    "    if (z == \"ff\") period <- c(\"2071-2100\")\n",
    "    spatialPlot(gridMean,\n",
    "                backdrop.theme = \"coastline\",\n",
    "                main = paste(\"CC. signal wrt 1975-2005\"),\n",
    "                ylab = period,\n",
    "                col.regions = cb,\n",
    "                at = at,\n",
    "                set.min = at[1], set.max = at[length(at)])  \n",
    "  }) \n",
    "}) %>% unlist(recursive = FALSE)   \n",
    "pdf(pdfOutput, width = 15, height = 10)   \n",
    "grid.arrange(grobs = figs, ncol = 3)                     \n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Time-series\n",
    "The following block of code allows for plotting the time-series of the climate change signals. For precipitation (temperature), we perform a loop over the validation metrics of interest: R01, SDII (Mean). For details about these metrics please see the manuscript or type `show.indices()` in a new cell. At each iteration of the loop we define a `doCall.args` list which contains the `aggr.y` arguments needed for the `aggregateGrid` function (note the validation is done at an annual basis), which is finally passed t  `do.call`. At the end of the loop, the resulting figures are saved in `pdf` format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Please select one:   -----------------------------------------------------------------------\n",
    "# Precipitation\n",
    "indices <- c(\"R01\",\"SDII\")\n",
    "variable  <- \"precip\"\n",
    "var <- \"pr\"\n",
    "# Temperature\n",
    "indices <- c(\"Mean\")\n",
    "variable <- \"temperature\"\n",
    "var <- \"tas\"\n",
    "\n",
    "## --------------------------------------------------------------------------------------------\n",
    "figs <- lapply(indices, FUN = function(zz) {\n",
    "  doCall.args <- list() \n",
    "  doCall.args[[\"aggr.y\"]] <- list()\n",
    "  \n",
    "  # The R01 do.call arguments  \n",
    "  if (zz == \"R01\")  {\n",
    "    doCall.args[[\"aggr.y\"]][[\"FUN\"]] <- \"index.freq\"\n",
    "    doCall.args[[\"aggr.y\"]][[\"freq.type\"]] <- \"rel\"\n",
    "    doCall.args[[\"aggr.y\"]][[\"condition\"]] <- \"GE\"\n",
    "    doCall.args[[\"aggr.y\"]][[\"threshold\"]] <- 1\n",
    "    ylim <- c(0.24,0.54)  \n",
    "  }\n",
    "  # The SDII do.call arguments  \n",
    "  if (zz == \"SDII\"){\n",
    "    doCall.args[[\"aggr.y\"]][[\"FUN\"]] <- \"index.meanGE\"\n",
    "    doCall.args[[\"aggr.y\"]][[\"threshold\"]] <- 1\n",
    "    ylim <- c(2,9)  \n",
    "  } \n",
    "  # The Mean do.call arguments    \n",
    "  if (zz == \"Mean\"){\n",
    "    doCall.args[[\"aggr.y\"]][[\"FUN\"]]   <- \"mean\"\n",
    "    doCall.args[[\"aggr.y\"]][[\"na.rm\"]] <- TRUE\n",
    "    ylim <- c(0,18)  \n",
    "  }    \n",
    "    \n",
    "  # We compute the index for the GCM ensemble. To do this, we loop over the temporal periods and then bind the serie along the time dimension with bindGrid function    \n",
    "  pred1 <- lapply(c(\"h\",\"ef\",\"mf\",\"ff\"), FUN = function(z) { \n",
    "    lapply(1:length(dName), FUN = function(zzz) {            \n",
    "      doCall.args[[\"grid\"]] <- loadGridData(dataset = paste0(\"./Data/\",variable,\"/y_\",z,\"_\",dName[zzz],\".nc4\"),var = var)\n",
    "      do.call(\"aggregateGrid\",doCall.args)\n",
    "    }) %>% bindGrid(dimension = \"member\")\n",
    "  }) %>% bindGrid(dimension = \"time\") \n",
    "  \n",
    "  # We compute the index for the RCM ensemble. To do this, we loop over the temporal periods and then bind the serie along the time dimension with bindGrid function\n",
    "  pred3 <- lapply(c(\"h\",\"ef\",\"mf\",\"ff\"), FUN = function(z) { \n",
    "    lapply(1:11, FUN = function(zzz) {            \n",
    "      doCall.args[[\"grid\"]] <- loadGridData(dataset = paste0(\"./Data/\",variable,\"/yRCM_\",z,\"_member\",zzz,\".nc4\"),var = var)\n",
    "      do.call(\"aggregateGrid\",doCall.args)\n",
    "    }) %>% bindGrid(dimension = \"member\")\n",
    "  }) %>% bindGrid(dimension = \"time\")\n",
    "    \n",
    "    # We compute the index for the DeepESD ensemble. To do this, we loop over the temporal periods and then bind the serie along the time dimension with bindGrid function    \n",
    "  pred2 <- lapply(c(\"h\",\"ef\",\"mf\",\"ff\"), FUN = function(z) { \n",
    "    lapply(1:length(dName), FUN = function(zzz) {            \n",
    "      doCall.args[[\"grid\"]] <- loadGridData(dataset = paste0(\"./Data/\",variable,\"/CNN_\",z,\"_\",dName[zzz],\".nc4\"),var = var)\n",
    "      do.call(\"aggregateGrid\",doCall.args)\n",
    "    }) %>% bindGrid(dimension = \"member\")\n",
    "  }) %>% bindGrid(dimension = \"time\")\n",
    "      \n",
    "  # We compute the index for the observed temporal serie \n",
    "  doCall.args[[\"grid\"]] <- y\n",
    "  y <- do.call(\"aggregateGrid\",doCall.args)\n",
    "  \n",
    "  # We call temporalPlot to plot the times-series \n",
    "  temporalPlot(\"OBS\" = y,\"GCM\" = pred1,\"RCM\" = pred3,\"CNN\" = pred2, cols = c(\"black\",\"red\",\"blue\",\"green\"),xyplot.custom = list(ylim = ylim))       \n",
    "})\n",
    "\n",
    "# Saving the resulting figures in .pdf format\n",
    "pdf(paste0(\"./figures/serie_\",var,\".pdf\"), width = 15, height = 4)\n",
    "grid.arrange(grobs = figs, ncol = 3)  \n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Technical specifications\n",
    "Please note this notebook was run on a machine with the following technical specifications:\n",
    "\n",
    "- Operating system: Ubuntu 18.04.3 LTS (64 bits)\n",
    "- Memory: 60 GiB\n",
    "- Processor: 2x Intel(R) Xeon(R) CPU E5-2670 0 @ 2.60GHz (16 cores, 32 threads)"
   ]
  }
 ],
 "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.6.1"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}