2020_Bano_GMD_FULL.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Configuration and Intercomparison of Deep Learning Neural Models for Statistical Downscaling *(code for full reproducibility)*\n",
    "### *Geoscientific Model Development preprint*\n",
    "### J. Baño-Medina, R. Manzanas and J. M. Gutiérrez\n",
    "https://doi.org/10.5194/gmd-2019-278\n",
    "\n",
    "GitHub repository at https://github.com/SantanderMetGroup/DeepDownscaling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "This notebook reproduces the results presented in the paper *Configuration and Intercomparison of Deep Learning Neural Models for Statistical Downscaling* by *J. Baño-Medina, R. Manzanas and J. M. Gutiérrez*, in *Geoscientific Model Development* (https://doi.org/10.5194/gmd-2019-278). In particular, the code developed herein concerns the downscaling of temperature and precipitation and therefore, their particularities are treated throughout the notebook in two different sections. Note that the programming lenguage is `R` and the technical specifications of the machine can be found at the end of the notebook. **The notebook takes around 5 days to fully reproduce the results. For this reason we refer the reader to a reduced version of the notebook in which we only consider a particular model configuration ([2020_Bano_GMD.ipynb](https://github.com/SantanderMetGroup/DeepDownscaling))**.\n",
    "\n",
    "**Note:** This notebook is written in the free programming language `R` and builds on the `R` framework [`climate4R`](https://github.com/SantanderMetGroup/climate4R) (C4R hereafter, conda and docker installations available), 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, via the [`downscaleR`](https://github.com/SantanderMetGroup/downscaleR) package; [Bedia et al. 2020](https://www.geosci-model-dev-discuss.net/gmd-2019-224/)) and visualization. The interested reader is referred to [Iturbide et al. 2019](https://www.sciencedirect.com/science/article/pii/S1364815218303049?via%3Dihub). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## 1. Loading libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "\n",
    "All the working steps rely on the [climate4R](https://github.com/SantanderMetGroup/climate4R) (conda and docker installations available), which is a set of libraries especifically developed to handle climate data (`loadeR`,`transformeR`,`downscaleR`, `visualizeR` and `climate4R.value`). In this study, [C4R](https://github.com/SantanderMetGroup/climate4R) is used for loading and post-processing, downscaling, validation and visualization. Different sectorial [notebooks](https://github.com/SantanderMetGroup/notebooks) are available illustrating the use of C4R functions. \n",
    "\n",
    "Deep learning models are indcluded as an extension of the downscaleR package: [`downscaleR.keras`](https://github.com/SantanderMetGroup/downscaleR.keras) which integrates *keras* in the C4R framework.\n",
    "\n",
    "Here we use a **specific version of the conda C4R installer (v1.3.0)** containing the package versions needed to fully reproduce the results of this notebook and including also the [`downscaleR.keras`](https://github.com/SantanderMetGroup/downscaleR.keras) package which is not included in the standard C4R conda distribution (visit this [site](https://docs.conda.io/projects/conda/en/latest/user-guide/install/index.html) for more information on how to obtain conda). Thus, below we illustrate how to create a new environment with conda to install the C4R libraries (please visit the [`climate4R`](https://github.com/SantanderMetGroup/climate4R) github repository for more details regarding the installation procedures). Note that the packages can be also installed individually, but this is more time consuming than using the conda installation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## --- type this in terminal to proceed with the installation of C4R with conda\n",
    "\n",
    "# conda create --name myEnvironment  # create a new environment\n",
    "# conda activate myEnvironment  # activate the environment\n",
    "# conda install -c conda-forge -c defaults -c santandermetgroup climate4r=1.3.0  # install climate4R v1.3.0 in the new environment\n",
    "# R # type 'R' to initialize R software"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once in R, we load the libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# C4R framework and auxiliary packages\n",
    "\n",
    "library(loadeR) # version 1.6.1      \n",
    "library(transformeR) # version 1.7.4\n",
    "library(downscaleR) # version 3.1.3\n",
    "library(visualizeR) # version 1.5.1\n",
    "library(climate4R.value) # version 0.0.2 (also relies on VALUE version 2.2.1)\n",
    "library(magrittr)\n",
    "library(gridExtra)\n",
    "library(RColorBrewer)\n",
    "library(sp)            \n",
    "library(downscaleR.keras) # version 0.0.2 (relies on keras version 2.2.2 and tensorflow version 2.0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "In order to avoid possible errors while running the notebook, you have to set the path to your desired working directory and create two files named \"Data\" and \"models\", that will contain the downscaled predictions and the trained deep models, respectively. Moreover, as we perform 2 distinct studies, one for precipitation and other for temperature, you should create 2 new directories named \"precip\" and \"temperature\" within the previous created directories (i.e., \"Data\" and \"models\"). An example of the latter would be \"personalpath/Data/temperature\". The predictions and models infered will be automatically saved in these folders and therefore not creating them will end into saving errors across the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "path = \"\"\n",
    "setwd(path)\n",
    "dir.create(\"Data\")\n",
    "dir.create(\"Data/precip/\")\n",
    "dir.create(\"models\")\n",
    "dir.create(\"models/temperature/\")\n",
    "dir.create(\"models/precip/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## 2. Loading data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "As explained in the paper, we have considered 20 large-scale variables from the ERA-Interim reanalysis as predictors and surface temperature from E-OBS as predictand. All these data can be loaded remotely from the [Santander Climate Data Service](http://meteo.unican.es/cds) (register [here](https://meteo.unican.es/trac/wiki/udg) freely to get a user), which provides access to various kinds of climate datasets (global and regional climate models, reanalysis, observations...). We will use the C4R packages [`loadeR`](https://github.com/SantanderMetGroup/loadeR) and [`transformeR`](https://github.com/SantanderMetGroup/transformeR) to load and postprocess the required information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "loginUDG(username = \"\", password = \"\") # login into the Santander CDS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We find the label associated to ERA-Interim via the `UDG.datasets()` function of loadeR: “ECMWF_ERA-Interim-ESD”. Then we load the predictors by calling `loadGridData` from loadeR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "variables <- c(\"z@500\",\"z@700\",\"z@850\",\"z@1000\",\n",
    "               \"hus@500\",\"hus@700\",\"hus@850\",\"hus@1000\",\n",
    "               \"ta@500\",\"ta@700\",\"ta@850\",\"ta@1000\",\n",
    "               \"ua@500\",\"ua@700\",\"ua@850\",\"ua@1000\",\n",
    "               \"va@500\",\"va@700\",\"va@850\",\"va@1000\")\n",
    "x <- lapply(variables, function(x) {\n",
    "  loadGridData(dataset = \"ECMWF_ERA-Interim-ESD\",\n",
    "               var = x,\n",
    "               lonLim = c(-10,32), # 22 points en total\n",
    "               latLim = c(36,72),  # 19 points en total\n",
    "               years = 1979:2008)\n",
    "}) %>% makeMultiGrid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We split into train (i.e., 1979-2002) and test (i.e., 2003-2008) datasets by using the `subsetGrid` function from transformeR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "xT <- subsetGrid(x,years = 1979:2002)\n",
    "xt <- subsetGrid(x,years = 2003:2008)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## 3. Temperature\n",
    "In this section we present the code needed to downscale temperature.  Once the predictors were loaded above we proceed to download the predictand dataset: E-OBS version 14 at a resolution of 0.5º. The E-OBS dataset is also accesible through the UDG portal. Thus, we load the temperature by calling again `loadGridData` and split it into train and test periods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "y <- loadGridData(dataset = \"E-OBS_v14_0.50regular\",\n",
    "                  var = \"tas\",\n",
    "                  lonLim = c(-10,32),\n",
    "                  latLim = c(36,72), \n",
    "                  years = 1979:2008)\n",
    "\n",
    "yT <- subsetGrid(y,years = 1979:2002)\n",
    "yt <- subsetGrid(y,years = 2003:2008)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We can take a look at the grid's resolutions of ERA-Interim (2º) and E-OBS (0.5º) in order to better visualize the gap we try to bridge with the downscaling (Figure 1 of the manuscript)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "cb <- colorRampPalette(brewer.pal(9, \"OrRd\"))(80)\n",
    "coords_x <- expand.grid(xt$xyCoords$x,xt$xyCoords$y) ; names(coords_x) <- c(\"x\",\"y\")\n",
    "coords_y <- expand.grid(yt$xyCoords$x,yt$xyCoords$y) ; names(coords_y) <- c(\"x\",\"y\")\n",
    "colsindex <- rev(brewer.pal(n = 9, \"RdBu\"))\n",
    "cb2 <- colorRampPalette(colsindex)\n",
    "\n",
    "pplot <- list()\n",
    "pplot[[1]] <- spatialPlot(climatology(subsetGrid(xt,var = \"ta@1000\")), backdrop.theme = \"coastline\",\n",
    "                          main = \"Temperature 1000 hPa (ERA-Interim)\",\n",
    "                          col.regions = cb2,\n",
    "                          at = seq(-3, 15, 1),\n",
    "                          set.min = -3, set.max = 15, colorkey = TRUE,\n",
    "                          sp.layout = list(list(SpatialPoints(coords_x), \n",
    "                                                first = FALSE, col = \"black\", \n",
    "                                                pch = 20, cex = 1)))\n",
    "pplot[[2]] <- spatialPlot(climatology(yt), backdrop.theme = \"coastline\", \n",
    "                          main = \"Temperature (E-OBS)\",\n",
    "                          col.regions = cb,\n",
    "                          at = seq(-3, 15, 1),\n",
    "                          set.min = -3, set.max = 15, colorkey = TRUE,\n",
    "                          sp.layout = list(list(SpatialPoints(coords_y), \n",
    "                                                first = FALSE, col = \"black\", \n",
    "                                                pch = 20, cex = 1)))\n",
    "\n",
    "lay = rbind(c(1,2))\n",
    "grid.arrange(grobs = pplot, layout_matrix = lay)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We can visualize some statistics of the train and test distributions, such as the climatology, or the percentiles 02th and 98th in order to gain knowledge about the observed data (Figure 2 of the manuscript). To compute the statistics we use `valueIndex` and `valueMeasure` functions from `climate4R.value` , a wrapper to the [`VALUE`](https://github.com/SantanderMetGroup/VALUE) package, which was developed in the [COST action VALUE](http://www.value-cost.eu/) for validation purposes. If interested in knowing more about the plotting options used here, the reader is referred to the documentation of the `spatialPlot` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "cb <- colorRampPalette(brewer.pal(9, \"OrRd\"))(80)\n",
    "colsindex <- rev(brewer.pal(n = 9, \"RdBu\"))\n",
    "cb2 <- colorRampPalette(colsindex)\n",
    "\n",
    "pplot <- at <- list()\n",
    "n1 <- 0; n2 <- 3\n",
    "indexNames <- c(\"Climatology\", \"P02\", \"P98\")\n",
    "for (indexName in indexNames) {\n",
    "  if (indexName == \"Climatology\") {\n",
    "    indexTrain <- valueIndex(yT,index.code = \"Mean\")$Index %>% redim()  \n",
    "    indexTest <- valueIndex(yt,index.code = \"Mean\")$Index %>% redim()\n",
    "    at[[1]] <- seq(-3, 15, 1); at[[2]] <- seq(-2, 2, 0.1)\n",
    "  }\n",
    "  if (indexName == \"P02\") {\n",
    "    indexTrain <- valueIndex(yT,index.code = \"P02\")$Index %>% redim()  \n",
    "    indexTest <- valueIndex(yt,index.code = \"P02\")$Index %>% redim()\n",
    "    at[[1]] <- seq(-20, 10, 1); at[[2]] <- seq(-2, 2, 0.1)\n",
    "  }\n",
    "  if (indexName == \"P98\") {\n",
    "    indexTrain <- valueIndex(yT,index.code = \"P98\")$Index %>% redim()  \n",
    "    indexTest <- valueIndex(yt,index.code = \"P98\")$Index %>% redim()\n",
    "    at[[1]] <- seq(10, 30, 1); at[[2]] <- seq(-2, 2, 0.1)\n",
    "  }\n",
    "  \n",
    "  for (i in 1:2) {\n",
    "    if (i == 1) {\n",
    "      dataset <- \"(train)\"; index <- indexTrain; n1 <- n1 + 1; n <- n1\n",
    "      value <- index$Data; colorbar <- cb\n",
    "    }\n",
    "    if (i == 2) {\n",
    "      indexTest <- gridArithmetics(indexTest,indexTrain,operator = \"-\")\n",
    "      dataset <- \"(test bias)\"; index <- indexTest; n2 <- n2 + 1; n <- n2\n",
    "      value <- abs(index$Data); colorbar <- cb2\n",
    "    }\n",
    "    pplot[[n]] <- spatialPlot(climatology(index), backdrop.theme = \"coastline\", \n",
    "                              main = paste(indexName,paste0(dataset,\":\"),\n",
    "                                           round(mean(value, na.rm = TRUE),digits = 2)),\n",
    "                              col.regions = colorbar,\n",
    "                              at = at[[i]],\n",
    "                              set.min = at[[i]][1], set.max = at[[i]][length(at[[i]])], \n",
    "                              colorkey = TRUE)\n",
    "  }\n",
    "}\n",
    "\n",
    "lay = rbind(c(1,2,3),\n",
    "            c(4,5,6))\n",
    "grid.arrange(grobs = pplot, layout_matrix = lay)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Once the data is loaded we standardize the predictors by calling `scaleGrid` function of transformeR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "xt <- scaleGrid(xt,xT, type = \"standardize\", spatial.frame = \"gridbox\") %>% redim(drop = TRUE)\n",
    "xT <- scaleGrid(xT, type = \"standardize\", spatial.frame = \"gridbox\") %>% redim(drop = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### 3.1 Downscaling with generalized linear models (temperature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "To downscale via generalized linear models (GLM) we rely on the `downscaleR` package of `C4R`. In particular, we use the `downscaleChunk` function of `downscaleR`. In the case of temperature, the generalized linear model has a gaussian family with link identity which is, in fact, an ordinary least squares regression. Therefore, we input to the function the predictor (x), the number of local predictors to be used (neighbours), the predictand (y) and the test set to apply the infered relationship as newdata. Finally we save the predictions which will be loaded during validation. Note that `downscaleChunk` temporarily creates .rda files in your working directory, containing the predictions per chunk. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "glmName <- c(\"glm1\",\"glm4\")\n",
    "neighs <- c(1,4)\n",
    "lapply(1:length(glmName), FUN = function(z) {\n",
    "  pred <- downscaleChunk(x = xT, y = yT, newdata = list(xt),\n",
    "                         method = \"GLM\", family = \"gaussian\", simulate = FALSE,\n",
    "                         prepareData.args = list(local.predictors = list(n=neighs[z], \n",
    "                                                                         vars = getVarNames(xT))))[[2]] %>% \n",
    "    redim(drop = TRUE)\n",
    "  save(pred,file = paste0(\"./Data/temperature/predictions_\",glmName[z],\".rda\"))\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### 3.2 Downscaling with deep neural networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "In the following code we define a function containing the deep learning topologies intercompared in the study (see Table 2 and Figure 3 of the manuscript for more information about these architectures)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "deepName <- c(\"CNN-LM\",\"CNN1\",\"CNN10\",\"CNN-PR\",\"CNNdense\")\n",
    "architectures <- function(architecture,input_shape,output_shape) {\n",
    "  if (architecture == \"CNN-LM\") {\n",
    "    inputs <- layer_input(shape = input_shape)\n",
    "    x = inputs\n",
    "    l1 = layer_conv_2d(x ,filters = 50, kernel_size = c(3,3), activation = 'linear', padding = \"same\")\n",
    "    l2 = layer_conv_2d(l1,filters = 25, kernel_size = c(3,3), activation = 'linear', padding = \"same\")\n",
    "    l3 = layer_conv_2d(l2,filters = 1, kernel_size = c(3,3), activation = 'linear', padding = \"same\")\n",
    "    l4 = layer_flatten(l3)\n",
    "    outputs = layer_dense(l4,units = output_shape)\n",
    "    model <- keras_model(inputs = inputs, outputs = outputs)\n",
    "  }\n",
    "  \n",
    "  if (architecture == \"CNN1\") {\n",
    "    inputs <- layer_input(shape = input_shape)\n",
    "    x = inputs\n",
    "    l1 = layer_conv_2d(x ,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",
    "    outputs = layer_dense(l4,units = output_shape)\n",
    "    model <- keras_model(inputs = inputs, outputs = outputs)\n",
    "  }\n",
    "  \n",
    "  if (architecture == \"CNN10\") {\n",
    "    inputs <- layer_input(shape = input_shape)\n",
    "    x = inputs\n",
    "    l1 = layer_conv_2d(x ,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",
    "    outputs = layer_dense(l4,units = output_shape)\n",
    "    model <- keras_model(inputs = inputs, outputs = outputs)\n",
    "  }\n",
    "  \n",
    "  if (architecture == \"CNN-PR\") {\n",
    "    inputs <- layer_input(shape = input_shape)\n",
    "    x = inputs\n",
    "    l1 = layer_conv_2d(x ,filters = 10, 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 = 50, kernel_size = c(3,3), activation = 'relu', padding = \"valid\")\n",
    "    l4 = layer_flatten(l3)\n",
    "    outputs = layer_dense(l4,units = output_shape)\n",
    "    model <- keras_model(inputs = inputs, outputs = outputs)\n",
    "  }\n",
    "  \n",
    "  if (architecture == \"CNNdense\") {\n",
    "    inputs <- layer_input(shape = input_shape)\n",
    "    x = inputs\n",
    "    l1 = layer_conv_2d(x ,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",
    "    l5 = layer_dense(l4,units = 50, activation = \"relu\")\n",
    "    l6 = layer_dense(l5,units = 50, activation = \"relu\")\n",
    "    outputs = layer_dense(l6,units = output_shape)\n",
    "    model <- keras_model(inputs = inputs, outputs = outputs)\n",
    "  }\n",
    "  \n",
    "  return(model)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We prepare the predictor and predictand datasets for integration with keras with the functions `prepareData.keras` and `prepareNewData.keras`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "xy.T <- prepareData.keras(xT,yT,\n",
    "                          first.connection = \"conv\",\n",
    "                          last.connection = \"dense\",\n",
    "                          channels = \"last\")\n",
    "xy.t <- prepareNewData.keras(xt,xy.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We loop over the topologies to train the deep models and predict over the test set. Unlike GLMs where there was a model per gridpoint, deep models perform multi-task, downscaling to all sites at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "lapply(1:length(deepName), FUN = function(z){\n",
    "  model <- architectures(architecture = deepName[z],\n",
    "                         input_shape = dim(xy.T$x.global)[-1],\n",
    "                         output_shape = dim(xy.T$y$Data)[2])\n",
    "  downscaleTrain.keras(obj = xy.T,\n",
    "                       model = model,\n",
    "                       clear.session = TRUE,\n",
    "                       compile.args = list(\"loss\" = \"mse\",\n",
    "                                           \"optimizer\" = optimizer_adam(lr = 0.0001)),\n",
    "                       fit.args = list(\"batch_size\" = 100,\n",
    "                                       \"epochs\" = 1000,\n",
    "                                       \"validation_split\" = 0.1,\n",
    "                                       \"verbose\" = 1,\n",
    "                                       \"callbacks\" = list(callback_early_stopping(patience = 30),\n",
    "                                                 callback_model_checkpoint(\n",
    "                                                 filepath=paste0('./models/temperature/',deepName[z],'.h5'),\n",
    "                                                                  monitor='val_loss', save_best_only=TRUE))))\n",
    "  pred <- downscalePredict.keras(newdata = xy.t,\n",
    "                                 model = list(\"filepath\" = \n",
    "                                                paste0(\"./models/temperature/\",deepName[z],\".h5\")),\n",
    "                                 C4R.template = yT,\n",
    "                                 clear.session = TRUE)\n",
    "  save(pred,file = paste0(\"./Data/temperature/predictions_\",deepName[z],\".rda\"))\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### 3.3 Validation of results\n",
    "In this code, we calculate the validation indices by using the library `climate4R.value` of `C4R`. In particular the indices used are: the Root Mean Squared Error (RMSE), the deseasonal perason correlation, the biases of the climatology and of the percentile 2th and 98th and the ratio of the standard deviations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "models <- c(\"glm1\",\"glm4\",\n",
    "            \"CNN-LM\",\"CNN1\",\"CNN10\",\n",
    "            \"CNN-PR\",\"CNNdense\")\n",
    "measures <- c(\"ts.RMSE\",\"ts.rp\",\"ratio\",rep(\"bias\",6))\n",
    "index <- c(rep(NA,2),\"sd\",\"Mean\",\"P02\",\"P98\",\"AC1\",\n",
    "           \"WarmAnnualMaxSpell\",\"ColdAnnualMaxSpell\")\n",
    "yt2 <- scaleGrid(yt,time.frame = \"daily\",window.width = 31) %>% redim(drop=TRUE)\n",
    "validation.list <- lapply(1:length(measures), FUN = function(z) {\n",
    "  lapply(1:length(models), FUN = function(zz){\n",
    "    args <- list()\n",
    "    load(paste0(\"./Data/temperature/predictions_\",models[zz],\".rda\"))\n",
    "    if (any(measures[z] == c(\"ts.rp\",\"ratio\"))) {\n",
    "      pred2 <- scaleGrid(pred,time.frame = \"daily\",window.width = 31) %>% redim(drop=TRUE)\n",
    "      args[[\"y\"]] <- yt2; args[[\"x\"]] <- pred2\n",
    "    } else {\n",
    "      args[[\"y\"]] <- yt; args[[\"x\"]] <- pred\n",
    "    }\n",
    "    args[[\"measure.code\"]] <- measures[z]\n",
    "    if (!is.na(index[z])) args[[\"index.code\"]] <- index[z]\n",
    "    do.call(\"valueMeasure\",args)$Measure\n",
    "  }) %>% makeMultiGrid()\n",
    "})\n",
    "save(validation.list, file = \"./Data/temperature/validation.rda\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Once the validation indices are calculated, we represent the results in boxplots (Figure 4 of the manuscript)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "ylabs <- c(\"RMSE (ºC)\",\"Cor. deseasonal\",\n",
    "           \"Standard dev. ratio\",\"bias (ºC)\",\n",
    "           \"bias P02 (ºC)\",\"bias P98 (ºC)\",\n",
    "           \"bias AC1\", \"bias WAMS (days)\",\n",
    "           \"bias CAMS (days)\")\n",
    "par(mfrow = c(3,3))\n",
    "lapply(1:length(validation.list), FUN = function(z) {\n",
    "  if (z == 1) {ylim <- c(0.75,1.75)}\n",
    "  if (z == 2) {ylim <- c(0.9,1)}\n",
    "  if (z == 3) {ylim <- c(0.9,1.05)}\n",
    "  if (z == 4) {ylim <- c(-0.2,0.5)}\n",
    "  if (z == 5) {ylim <- c(-0.5,1.5)}\n",
    "  if (z == 6) {ylim <- c(-1,1)}\n",
    "  if (z == 7) {ylim <- c(-0.1,0.1)}\n",
    "  if (any(z == c(8,9))) {ylim <- c(-3,3)}\n",
    "  index <- (validation.list[[z]] %>% redim(drop = TRUE))$Data\n",
    "  dim(index) <- c(nrow(index),prod(dim(index)[2:3]))\n",
    "  indLand <- (!apply(index,MARGIN = 2,anyNA)) %>% which()\n",
    "  index <- index[,indLand] %>% t()\n",
    "  mglm4 <- median(index[,2],na.rm = TRUE)\n",
    "  perc <- apply(index,MARGIN = 2,FUN = function(z) quantile(z,probs = c(0.1,0.9)))\n",
    "  boxplot(index, outline = FALSE, ylim = ylim, range = 0.0001, ylab = ylabs[z], asp = 1)\n",
    "  lines(c(0,8),c(mglm4,mglm4), col = \"red\")\n",
    "  for (i in 1:ncol(index)) lines(c(i,i),perc[,i], lty = 2)\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "In order to obtain a spatial representation of the validation indices computed above we use the function `spatialPlot` of visualizeR. In particular, we plot the deseasonal correlation and the biases of the percentile 2th, 98th and the mean for the glm1, glm4 and CNN10 models (Figure 5 of the manuscript)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "ylabs <- c(\"glm1\",\"glm4\",NA,NA,\"CNN10\")\n",
    "mains <- c(NA,\"Cor. deseasonal\",NA,\"bias (ºC)\",\"bias P02 (ºC)\",\"bias P98 (ºC)\")\n",
    "cb <- colorRampPalette(brewer.pal(9, \"OrRd\"))(80)\n",
    "colsindex <- rev(brewer.pal(n = 9, \"RdBu\"))\n",
    "cb2 <- colorRampPalette(colsindex)\n",
    "validation.plots <- lapply(c(2,4,5,6),FUN = function(z) {\n",
    "  lapply(c(1,2,5),FUN = function(zz) {\n",
    "    if (z == 2) {\n",
    "      at <- seq(0.85, 1, 0.005); colorbar <- cb\n",
    "    } else {\n",
    "      at <- seq(-2, 2, 0.1); colorbar <- cb2\n",
    "    }\n",
    "    index <- subsetDimension(validation.list[[z]],dimension = \"var\",indices = zz) %>% \n",
    "      redim(drop = TRUE)\n",
    "    spatialPlot(index, backdrop.theme = \"coastline\",\n",
    "                ylab = ylabs[zz],\n",
    "                main = paste(mains[z],\n",
    "                             round(mean(abs(index$Data), na.rm = TRUE), digits = 2)),\n",
    "                col.regions = colorbar,\n",
    "                at = at,\n",
    "                set.min = at[1], set.max = at[length(at)], colorkey = TRUE)\n",
    "  })\n",
    "})\n",
    "lay = cbind(1:3,4:6,7:9,10:12)\n",
    "grid.arrange(grobs = unlist(validation.plots,recursive = FALSE), layout_matrix = lay)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To provide more insight into the extrapolation abilities of the statistical downscaling methods tested herein, we compute the bias of the frequency of days above the percentile 95th (over the train period) in the observed and predicted (GLM4 and CNN10) time series with respect to the same frequency in the observed train values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grids <- c(\"yt\",\"predictions_glm1\",\"predictions_glm4\",\"predictions_CNN10\")\n",
    "gridsNames <- c(\"Test\",\"GLM1\",\"GLM4\",\"CNN10\")\n",
    "p.plots <- lapply(c(0.95,0.99), FUN = function(zzz){\n",
    "  pX <- apply(yT$Data,MARGIN = c(2,3),FUN = function(z) quantile(z,probs = zzz,na.rm=TRUE))\n",
    "  if (zzz == 0.95) at <- seq(-0.1, 0.1, length.out = 30)\n",
    "  if (zzz == 0.99) at <- seq(-0.05, 0.05, length.out = 30)\n",
    "  lapply(1:length(grids),FUN = function(zz) {\n",
    "    grid <- get(load(paste0(\"./Data/temperature/\",grids[zz],\".rda\")))\n",
    "    freq.t <- grid\n",
    "    dimNames <- attr(grid$Data,\"dimensions\")\n",
    "    days.t <- lapply(1:getShape(grid,\"time\"),FUN = function(z){\n",
    "      grid$Data[z,,] > pX \n",
    "    }) %>% abind::abind(along = 0)\n",
    "    freq.t$Data <- apply(days.t,MARGIN = c(2,3),FUN = function(z) mean(z,na.rm = TRUE))\n",
    "    attr(freq.t$Data,\"dimensions\") <- dimNames[-1]\n",
    "    freq.t <- gridArithmetics(freq.t,1-zzz,operator = \"-\")\n",
    "    spatialPlot(freq.t,\n",
    "                backdrop.theme = \"coastline\",\n",
    "                col.regions = rev(brewer.pal(n = 9, \"RdBu\")) %>% colorRampPalette(),\n",
    "                main = gridsNames[zz],\n",
    "                at = at,\n",
    "                set.min = at[[1]], set.max = at[[length(at)]], \n",
    "                colorkey = TRUE)\n",
    "  })\n",
    "})\n",
    "# P95\n",
    "grid.arrange(grobs = p.plots[[1]], ncol = 4)\n",
    "# P99\n",
    "grid.arrange(grobs = p.plots[[2]], ncol = 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## 4. Precipitation\n",
    "In this section we present the code needed to downscale precipitation. Though the steps taken are very similar to those of temperature there are some particularities that are good to mention. We start by loading the precipitation using `loadGridData` and to convert to binary the precipitation (values greater than 1mm/day of rain are set to 1 and the rest to 0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "y <- loadGridData(dataset = \"E-OBS_v14_0.50regular\",\n",
    "                  var = \"pr\",lonLim = c(-10,32),\n",
    "                  latLim = c(36,72), \n",
    "                  years = 1979:2008)\n",
    "\n",
    "yT <- subsetGrid(y,years = 1979:2002)\n",
    "yT_bin <- binaryGrid(yT,threshold = 1,condition = \"GT\")\n",
    "yt <- subsetGrid(y,years = 2003:2008)\n",
    "yt_bin <- binaryGrid(yt,threshold = 1,condition = \"GT\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We can take a look at the grid resolutions of ERA-Interim (2º) and E-OBS (0.5º) by plotting the specific humidity at 1000hPa and the observed precipitation (Figure 1 of the manuscript). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "colsindex <- brewer.pal(n = 9, \"BrBG\")\n",
    "cb <- colorRampPalette(colsindex)\n",
    "coords_x <- expand.grid(xt$xyCoords$x,xt$xyCoords$y) ; names(coords_x) <- c(\"x\",\"y\")\n",
    "coords_y <- expand.grid(yt$xyCoords$x,yt$xyCoords$y) ; names(coords_y) <- c(\"x\",\"y\")\n",
    "\n",
    "pplot <- list()\n",
    "pplot[[1]] <- spatialPlot(climatology(subsetGrid(xt,var = \"hus@1000\")), backdrop.theme = \"coastline\",\n",
    "                          main = \"Q1000 (ERA-Interim)\",\n",
    "                          col.regions = cb,\n",
    "                          at = seq(0,0.01, 0.0001),\n",
    "                          set.min = 0, set.max = 0.01, colorkey = TRUE,\n",
    "                          sp.layout = list(list(SpatialPoints(coords_x), \n",
    "                                                first = FALSE, col = \"black\", \n",
    "                                                pch = 20, cex = 1)))\n",
    "pplot[[2]] <- spatialPlot(climatology(yt), backdrop.theme = \"coastline\", \n",
    "                          main = \"Precipitation (E-OBS)\",\n",
    "                          col.regions = cb,\n",
    "                          at = seq(0, 4, 0.1),\n",
    "                          set.min = 0, set.max = 4, colorkey = TRUE,\n",
    "                          sp.layout = list(list(SpatialPoints(coords_y), \n",
    "                                                first = FALSE, col = \"black\", \n",
    "                                                pch = 20, cex = 1)))\n",
    "\n",
    "lay = rbind(c(1,2))\n",
    "grid.arrange(grobs = pplot, layout_matrix = lay)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We can visualize some statistics of the train and test distributions, such as the climatology, or the percentiles 2th and 98th in order to gain knowledge about the observed data (Figure 2 of the manuscript). To compute the statistics we use the library `climate4R.value`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "pplot <- at <- list()\n",
    "n1 <- 0; n2 <- 3\n",
    "indexNames <- c(\"Climatology\", \"Frequency of rain\", \"P98\")\n",
    "for (indexName in indexNames) {\n",
    "  if (indexName == \"Climatology\") {\n",
    "    indexTrain <- valueIndex(yT,index.code = \"Mean\")$Index %>% redim()  \n",
    "    indexTest <- valueIndex(yt,index.code = \"Mean\")$Index %>% redim()\n",
    "    at[[1]] <- seq(0, 4, 0.1); at[[2]] <- seq(-1, 1, 0.1)\n",
    "  }\n",
    "  if (indexName == \"Frequency of rain\") {\n",
    "    indexTrain <- valueIndex(yT_bin,index.code = \"Mean\")$Index %>% redim()  \n",
    "    indexTest <- valueIndex(yt_bin,index.code = \"Mean\")$Index %>% redim()\n",
    "    at[[1]] <- seq(0, 0.5, 0.01); at[[2]] <- seq(-0.1, 0.1, 0.01)\n",
    "  }\n",
    "  if (indexName == \"P98\") {\n",
    "    indexTrain <- valueIndex(yT,index.code = \"P98\")$Index %>% redim()  \n",
    "    indexTest <- valueIndex(yt,index.code = \"P98\")$Index %>% redim()\n",
    "    at[[1]] <- seq(10, 20, 0.25); at[[2]] <- seq(-5, 5, 0.2)\n",
    "  }\n",
    "  \n",
    "  for (i in 1:2) {\n",
    "    if (i == 1) {\n",
    "      dataset <- \"(train)\"; index <- indexTrain; n1 <- n1 + 1; n <- n1\n",
    "      }\n",
    "    if (i == 2) {\n",
    "      indexTest <- gridArithmetics(indexTest,indexTrain,operator = \"-\")\n",
    "      dataset <- \"(test bias)\"; index <- indexTest; n2 <- n2 + 1; n <- n2\n",
    "      }\n",
    "    pplot[[n]] <- spatialPlot(climatology(index), backdrop.theme = \"coastline\", \n",
    "                              main = paste(indexName,paste0(dataset,\":\"),\n",
    "                                           round(mean(abs(index$Data), na.rm = TRUE),digits = 2)),\n",
    "                              col.regions = cb,\n",
    "                              at = at[[i]],\n",
    "                              set.min = at[[i]][1], set.max = at[[i]][length(at[[i]])], \n",
    "                              colorkey = TRUE)\n",
    "  }\n",
    "}\n",
    "\n",
    "lay = rbind(c(1,2,3),\n",
    "            c(4,5,6))\n",
    "grid.arrange(grobs = pplot, layout_matrix = lay)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### 4.1 Downscaling with generalized linear models (precipitation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "In the case of precipitation, there are 2 generalized linear models (GLMs): one to predict the occurrence of precipitation with binomial family and link logit and another to predict tha rainfall amount based on a gamma family and link logarithmic. We train the model and predict on the test set using the package `downscaleR`. In particular we train two configurations using as predictor data the closest neighbour (glm1) or the 4 closest neigbours (glm4). We obtain stochastic predictions by setting the parameter `simulate` to `TRUE`. We refer the reader to the [downscaleR wiki](https://github.com/SantanderMetGroup/downscaleR/wiki) for a detailed explanation about the working of the downscaleR functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "simulateName <- c(\"deterministic\",\"stochastic\")\n",
    "glmName <- c(\"glm1\",\"glm4\")\n",
    "neighs <- c(1,4)\n",
    "y.ocu <- binaryGrid(yT,condition = \"GT\",threshold = 1)\n",
    "y.rest <- gridArithmetics(yT,1,operator = \"-\")\n",
    "simulateGLM <- c(FALSE,TRUE)\n",
    "lapply(1:length(glmName), FUN = function(z){\n",
    "  lapply(1:length(simulateGLM),FUN = function(zz) {\n",
    "    pred <- downscaleChunk(x = xT, y = y.ocu, newdata = list(xt),\n",
    "                           method = \"GLM\", \n",
    "                           family = binomial(link = \"logit\"), \n",
    "                           simulate = simulateGLM[zz],\n",
    "                           prepareData.args = list(local.predictors = list(n=neighs[z],\n",
    "                                                                           vars = getVarNames(xT)))\n",
    "    )\n",
    "    pred_ocu_train <- pred[[1]] %>% redim(drop = TRUE)\n",
    "    pred_ocu <- pred[[2]] %>% redim(drop = TRUE)\n",
    "    rm(pred)\n",
    "    pred_amo <- downscaleChunk(x = xT, y = y.rest, newdata = list(xt),\n",
    "                               method = \"GLM\", \n",
    "                               family = Gamma(link = \"log\"), \n",
    "                               simulate = simulateGLM[zz],\n",
    "                               condition = \"GT\", threshold = 0,\n",
    "                               prepareData.args = list(local.predictors = list(n=neighs[z], \n",
    "                                                                               vars = getVarNames(xT))))[[2]] %>% \n",
    "      gridArithmetics(1,operator = \"+\") %>% redim(drop = TRUE)\n",
    "    pred_bin <- binaryGrid(pred_ocu,ref.obs = yT_bin,ref.pred = pred_ocu_train); rm(pred_ocu_train)\n",
    "    save(pred_bin,pred_ocu,pred_amo,\n",
    "         file = paste0(\"./Data/precip/predictions_\",simulateName[zz],\"_\",glmName[z],\".rda\"))\n",
    "  })\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### 4.2 Downscaling with deep neural networks (precipitation)\n",
    "As done with temperature we define a function containing the topologies intercompared in the study. The main difference with respect to the downscaling of temperature is that for precipitation we infer the probability of rain and the shape and scale parameters of a *Gamma* distribution, and therefore there are three output parameters per predictand's gridpoint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "deepName <- c(\"CNN-LM\",\"CNN1\",\"CNN10\",\"CNN-PR\",\"CNNdense\")\n",
    "architectures <- function(architecture,input_shape,output_shape) {\n",
    "  if (architecture == \"CNN-LM\") {\n",
    "    inputs <- layer_input(shape = input_shape)\n",
    "    x = inputs\n",
    "    l1 = layer_conv_2d(x ,filters = 50, kernel_size = c(3,3), activation = 'linear', padding = \"same\")\n",
    "    l2 = layer_conv_2d(l1,filters = 25, kernel_size = c(3,3), activation = 'linear', padding = \"same\")\n",
    "    l3 = layer_conv_2d(l2,filters = 1, kernel_size = c(3,3), activation = 'linear', padding = \"same\")\n",
    "    l4 = layer_flatten(l3)\n",
    "    parameter1 = layer_dense(l4,units = output_shape, activation = \"sigmoid\")\n",
    "    parameter2 = layer_dense(l4,units = output_shape)\n",
    "    parameter3 = layer_dense(l4,units = output_shape)\n",
    "    outputs = layer_concatenate(list(parameter1,parameter2,parameter3))\n",
    "    model <- keras_model(inputs = inputs, outputs = outputs)\n",
    "  }\n",
    "  \n",
    "  if (architecture == \"CNN1\") {\n",
    "    inputs <- layer_input(shape = input_shape)\n",
    "    x = inputs\n",
    "    l1 = layer_conv_2d(x ,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",
    "    parameter1 = layer_dense(l4,units = output_shape, activation = \"sigmoid\")\n",
    "    parameter2 = layer_dense(l4,units = output_shape)\n",
    "    parameter3 = layer_dense(l4,units = output_shape)\n",
    "    outputs = layer_concatenate(list(parameter1,parameter2,parameter3))\n",
    "    model <- keras_model(inputs = inputs, outputs = outputs)\n",
    "  }\n",
    "  \n",
    "  if (architecture == \"CNN10\") {\n",
    "    inputs <- layer_input(shape = input_shape)\n",
    "    x = inputs\n",
    "    l1 = layer_conv_2d(x ,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",
    "    parameter1 = layer_dense(l4,units = output_shape, activation = \"sigmoid\")\n",
    "    parameter2 = layer_dense(l4,units = output_shape)\n",
    "    parameter3 = layer_dense(l4,units = output_shape)\n",
    "    outputs = layer_concatenate(list(parameter1,parameter2,parameter3))\n",
    "    model <- keras_model(inputs = inputs, outputs = outputs)\n",
    "  }\n",
    "  \n",
    "  if (architecture == \"CNN-PR\") {\n",
    "    inputs <- layer_input(shape = input_shape)\n",
    "    x = inputs\n",
    "    l1 = layer_conv_2d(x ,filters = 10, 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 = 50, kernel_size = c(3,3), activation = 'relu', padding = \"valid\")\n",
    "    l4 = layer_flatten(l3)\n",
    "    parameter1 = layer_dense(l4,units = output_shape, activation = \"sigmoid\")\n",
    "    parameter2 = layer_dense(l4,units = output_shape)\n",
    "    parameter3 = layer_dense(l4,units = output_shape)\n",
    "    outputs = layer_concatenate(list(parameter1,parameter2,parameter3))\n",
    "    model <- keras_model(inputs = inputs, outputs = outputs)\n",
    "  }\n",
    "  \n",
    "  if (architecture == \"CNNdense\") {\n",
    "    inputs <- layer_input(shape = input_shape)\n",
    "    x = inputs\n",
    "    l1 = layer_conv_2d(x ,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",
    "    l5 = layer_dense(l4,units = 50, activation = \"relu\")\n",
    "    l6 = layer_dense(l5,units = 50, activation = \"relu\")\n",
    "    parameter1 = layer_dense(l6,units = output_shape, activation = \"sigmoid\")\n",
    "    parameter2 = layer_dense(l6,units = output_shape)\n",
    "    parameter3 = layer_dense(l6,units = output_shape)\n",
    "    outputs = layer_concatenate(list(parameter1,parameter2,parameter3))\n",
    "    model <- keras_model(inputs = inputs, outputs = outputs)\n",
    "  }\n",
    "  \n",
    "  return(model)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We prepare the predictor and predictand datasets. We substract 1 to the precipitation (*yT*) to center the conditional Gamma distribution in 0 and convert to 0 the negative values with binaryGrid in order to avoid issues with logarithms when minimizing the loss function. This will be later added to the prediction output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "xy.T <- prepareData.keras(xT,binaryGrid(gridArithmetics(yT,1,operator = \"-\"),\n",
    "                                        condition = \"GE\",\n",
    "                                        threshold = 0,\n",
    "                                        partial = TRUE),\n",
    "                          first.connection = \"conv\",\n",
    "                          last.connection = \"dense\",\n",
    "                          channels = \"last\")\n",
    "xy.tT <- prepareNewData.keras(xT,xy.T)\n",
    "xy.t <- prepareNewData.keras(xt,xy.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We loop over the topologies to train the deep models and predict over the test set. Unlike GLMs where there was two models, one for the occurrence of rain and other for the amount of rain, with deep learning we minimize the negative log-likelihood of a *Bernouilli Gamma* distribution and therefore both the occurrence and quantity of rain for a given day are derived from these infered parameters. The custom loss function `bernouilliGamma.loss_function` is part of `downscaleR.keras`. In addition we use the function `bernouilliGamma.statistics` (also from `downscaleR.keras`) to compute the deterministic (i.e., expectance of the conditional distribution) or the stochastic (i.e., sample from the conditional distribution) prediction. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "simulateName <- c(\"deterministic\",\"stochastic\")\n",
    "simulateDeep <- c(FALSE,TRUE)\n",
    "lapply(1:length(deepName), FUN = function(z){\n",
    "  model <- architectures(architecture = deepName[z],\n",
    "                         input_shape = dim(xy.T$x.global)[-1],\n",
    "                         output_shape = dim(xy.T$y$Data)[2])\n",
    "  downscaleTrain.keras(obj = xy.T,\n",
    "             model = model,\n",
    "             clear.session = TRUE,\n",
    "             compile.args = list(\"loss\" = bernouilliGamma.loss_function(last.connection = \"dense\"),\n",
    "                                           \"optimizer\" = optimizer_adam(lr = 0.0001)),\n",
    "                       fit.args = list(\"batch_size\" = 100,\n",
    "                            \"epochs\" = 1000,\n",
    "                            \"validation_split\" = 0.1,\n",
    "                            \"verbose\" = 1,\n",
    "                            \"callbacks\" = list(callback_early_stopping(patience = 30),\n",
    "                                callback_model_checkpoint(filepath=paste0('./models/precip/',deepName[z],'.h5'),\n",
    "                                monitor='val_loss', save_best_only=TRUE))))\n",
    "  lapply(1:length(simulateDeep),FUN = function(zz) {\n",
    "    pred_ocu_train <- downscalePredict.keras(newdata = xy.tT,\n",
    "                                             model = list(\"filepath\" = \n",
    "                                                      paste0(\"./models/precip/\",deepName[z],\".h5\"), \n",
    "                                                     \"custom_objects\" = \n",
    "                                                     c(\"custom_loss\" = \n",
    "                                                      bernouilliGamma.loss_function(\n",
    "                                                        last.connection = \"dense\"))),\n",
    "                                             C4R.template = yT,\n",
    "                                             clear.session = TRUE) %>% \n",
    "      subsetGrid(var = \"pr1\")\n",
    "    pred <- downscalePredict.keras(newdata = xy.t,\n",
    "                                   model = list(\"filepath\" = \n",
    "                                              paste0(\"./models/precip/\",deepName[z],\".h5\"), \n",
    "                                              \"custom_objects\" = \n",
    "                                              c(\"custom_loss\" = \n",
    "                                              bernouilliGamma.loss_function(last.connection = \"dense\"))),\n",
    "                                   C4R.template = yT,\n",
    "                                   clear.session = TRUE) \n",
    "    pred <- bernouilliGamma.statistics(p = subsetGrid(pred,var = \"pr1\"),\n",
    "                                       alpha = subsetGrid(pred,var = \"pr2\"),\n",
    "                                       beta = subsetGrid(pred,var = \"pr3\"),\n",
    "                                       simulate = simulateDeep[zz],\n",
    "                                       bias = 1)\n",
    "    pred_ocu <- subsetGrid(pred,var = \"probOfRain\") %>% redim(drop = TRUE)\n",
    "    pred_amo <- subsetGrid(pred,var = \"amountOfRain\") %>% redim(drop = TRUE)\n",
    "    pred_bin <- binaryGrid(pred_ocu,ref.obs = yT_bin,ref.pred = pred_ocu_train); rm(pred_ocu_train)\n",
    "    save(pred_bin,pred_ocu,pred_amo,file = \n",
    "            paste0(\"./Data/precip/predictions_\",simulateName[zz],\"_\",deepName[z],\".rda\"))\n",
    "  })\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### 4.3 Validation of results (precipitation)\n",
    "In this code, we calculate the validation indices by using the library `climate4R.value` of `C4R`. In particular the indices used are: the Root Mean Squared Error (RMSE), the deseasonal perason correlation, the biases of the climatology and of the percentile 2th and 98th and the ratio of the standard deviations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "simulateName <- c(rep(\"deterministic\",5),\"stochastic\",rep(\"deterministic\",3))\n",
    "models <- c(\"glm1\",\"glm4\",\n",
    "               \"CNN-LM\",\"CNN1\",\"CNN10\",\n",
    "               \"CNN-PR\",\"CNNdense\")\n",
    "measures <- c(\"ts.rocss\",\"ts.RMSE\",\"ts.rs\",rep(\"biasRel\",3),rep(\"bias\",3))\n",
    "index <- c(rep(NA,3),\"Mean\",rep(\"P98\",2),\"AnnualCycleRelAmp\",\n",
    "           \"WetAnnualMaxSpell\",\"DryAnnualMaxSpell\")\n",
    "validation.list <- lapply(1:length(measures), FUN = function(z) {\n",
    "  lapply(1:length(models), FUN = function(zz){\n",
    "    args <- list()\n",
    "    load(paste0(\"./Data/precip/predictions_\",simulateName[z],\"_\",models[zz],\".rda\"))\n",
    "    if (simulateName[z] == \"deterministic\") {\n",
    "      pred <- gridArithmetics(pred_bin,pred_amo,operator = \"*\")\n",
    "      if (measures[z] == \"ts.rocss\") {\n",
    "        args[[\"y\"]] <- yt_bin; args[[\"x\"]] <- pred_ocu\n",
    "      } else if (measures[z] == \"ts.RMSE\") {\n",
    "        args[[\"y\"]] <- yt; args[[\"x\"]] <- pred_amo\n",
    "        args[[\"condition\"]] = \"GT\"; args[[\"threshold\"]] = 1; args[[\"which.wetdays\"]] = \"Observation\"  \n",
    "      } else {\n",
    "        args[[\"y\"]] <- yt; args[[\"x\"]] <- pred\n",
    "      }\n",
    "    } else {\n",
    "      pred <- gridArithmetics(pred_ocu,pred_amo,operator = \"*\")\n",
    "      args[[\"y\"]] <- yt; args[[\"x\"]] <- pred\n",
    "    }\n",
    "    args[[\"measure.code\"]] <- measures[z]\n",
    "    if (!is.na(index[z])) args[[\"index.code\"]] <- index[z]\n",
    "    do.call(\"valueMeasure\",args)$Measure\n",
    "  }) %>% makeMultiGrid()\n",
    "})\n",
    "save(validation.list, file = \"./Data/precip/validation.rda\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Once the validation indices are calculated, we represent the results in boxplots (Figure 6 of the manuscript)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "par(mfrow = c(3,3)) \n",
    "ylabs <- c(\"ROCSS\",\"RMSE (wet days, mm)\",\n",
    "           \"Spearman Corr.\",\"biasRel(%)\",\n",
    "           \"biasRel P98 (DET, %)\",\"biasRel P98 (STO, %)\",\n",
    "           \"Annual Cycle Rel. Amplitude\", \"biasRel WetAMS (days)\",\n",
    "           \"biasRel DryAMS (days)\")\n",
    "lapply(1:length(validation.list), FUN = function(z) {\n",
    "  if (z == 1) {ylim <- c(0.65,0.9)}\n",
    "  if (z == 2) {ylim <- c(3,6.5)}\n",
    "  if (z == 3) {ylim <- c(0.5,0.8)}\n",
    "  if (z == 4) {ylim <- c(-0.2,0.2)}\n",
    "  if (z == 5) {ylim <- c(-0.4,0.0)}\n",
    "  if (z == 6) {ylim <- c(-0.2,0.2)}\n",
    "  if (z == 7) {ylim <- c(-1,1)}\n",
    "  if (any(z == c(8,9))) {ylim <- c(-1,1)}\n",
    "  index <- (validation.list[[z]] %>% redim(drop = TRUE))$Data\n",
    "  dim(index) <- c(nrow(index),prod(dim(index)[2:3]))\n",
    "  indLand <- (!apply(index,MARGIN = 2,anyNA)) %>% which()\n",
    "  index <- index[,indLand] %>% t()\n",
    "  mglm4 <- median(index[,2],na.rm = TRUE)\n",
    "  perc <- apply(index,MARGIN = 2,FUN = function(z) quantile(z,probs = c(0.1,0.9)))\n",
    "  boxplot(index, outline = FALSE, ylim = ylim, range = 0.0001, ylab = ylabs[z], asp = 1)\n",
    "  lines(c(0,8),c(mglm4,mglm4), col = \"red\")\n",
    "  for (i in 1:ncol(index)) lines(c(i,i),perc[,i], lty = 2)\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "In order to obtain a spatial representation of the validation indices computed above we use the function `spatialPlot` of `visualizeR`. In particular, we plot the ROCSS, the spearman correlation and the relative biases of the mean and the P98 for the glm1, glm4 and CNN1 models (Figure 7 of the manuscript)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "ylabs <- c(\"glm1\",\"glm4\",NA,\"CNN1\")\n",
    "mains <- c(\"ROCSS\",NA,\"Spearman Corr.\",\"biasRel\",\"biasRel P98\")\n",
    "cb <- colorRampPalette(brewer.pal(9, \"OrRd\"))(80)\n",
    "colsindex <- rev(brewer.pal(n = 9, \"RdBu\"))\n",
    "cb2 <- colorRampPalette(colsindex)\n",
    "validation.plots <- lapply(c(1,3,4,5),FUN = function(z) {\n",
    "  lapply(c(1,2,4),FUN = function(zz) {\n",
    "    if (z == 1) {\n",
    "      at <- seq(0.5, 1, 0.01); colorbar <- cb\n",
    "    } else if (z == 3) {\n",
    "      at <- seq(0.5, 1, 0.02); colorbar <- cb\n",
    "    } else {\n",
    "      at <- seq(-0.5, 0.5, 0.01); colorbar <- cb2\n",
    "    }\n",
    "    index <- subsetDimension(validation.list[[z]],dimension = \"var\",indices = zz) %>% redim(drop = TRUE)\n",
    "    spatialPlot(index, backdrop.theme = \"coastline\",\n",
    "                ylab = ylabs[zz],\n",
    "                main = paste(mains[z],\n",
    "                             round(mean(abs(index$Data), na.rm = TRUE), digits = 2)),\n",
    "                col.regions = colorbar,\n",
    "                at = at,\n",
    "                set.min = at[1], set.max = at[length(at)], colorkey = TRUE)\n",
    "  })\n",
    "})\n",
    "lay = cbind(1:3,4:6,7:9,10:12)\n",
    "grid.arrange(grobs = unlist(validation.plots,recursive = FALSE), layout_matrix = lay)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## 5. Technical aspects\n",
    "To perform all the stages involved in this study we relied on the local machine described below.\n",
    "  \n",
    "1. Local Machine (HP-ProDesk-600-G2-MT)\n",
    "     + Operating system: ubuntu 4.15.0-72-generic\n",
    "     + Memory: 15.6 GiB\n",
    "     + Processor: Intel® Core™ i7-6700 CPU @ 3.40GHz × 8 \n",
    "     + SO: 64 bits\n",
    "     + Disc: 235.1 GiB"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}