2018_climate4R_example2.Rmd

---
title: "Full code for Example 2 of the paper `climate4R: An R-based Open Framework for Reproducible Climate Data Access and Post-processing'"
author: "M. Iturbide, J. Bedia, S. Herrera, J. Baño-Medina, J. Fernández, M. D. Frías, R. Manzanas, D. San Martín, E. Cimadevilla, A.S. Cofiño, J. M. Gutiérrez"
date: "`r Sys.Date()`"
csl: elsarticle.csl
header-includes:
  - \usepackage[font={small}]{caption}
output: 
    rmarkdown::pdf_document:
        fig_caption: yes
        toc: yes
        pandoc_args: [
      "--number-sections",
      "--number-offset=0"
    ] 
vignette: >
  %\VignetteIndexEntry{mopa within the climate4R ecosystem}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
urlcolor: blue
---

```{r set, results='hide', message=FALSE, echo=FALSE}
 knitr::opts_chunk$set(fig.width = 6, fig.height = 4, cache = TRUE, cache.path = "./cache/ex2/", fig.path = "./cache/ex2/figs") 
```

```{r, echo=FALSE, message=FALSE, warning=FALSE, cache=FALSE}
options(java.parameters = "-Xmx8000m")
```


# Introduction

This worked example contains the full code that reproduces the 2nd example of the paper "climate4R: An R-based Framework for Climate Data Access, Post-processing and Bias Correction" (Sec. 6 of the manuscript). The main section is divided in additional subsections to help with the understanding of the different code chunks. All operations hereinafter are performed with the core packages of climate4R, excepting (1) package installation and (2) the creation of color palettes, for which packages `devtools` and `RColorBrewer` are used respectively:

(1) Package installation:
```{r, eval=FALSE}
library(devtools)
install_github(c("SantanderMetGroup/loadeR",
                 "SantanderMetGroup/loadeR.java",
                 "SantanderMetGroup/transformeR",
                 "SantanderMetGroup/visualizeR",
                 "SantanderMetGroup/downscaleR",
                 "SantanderMetGroup/climate4R.climdex")
```

```{r, message=FALSE, warning=FALSE}
library(loadeR)
library(transformeR)
library(visualizeR)
library(downscaleR)
library(climate4R.climdex)
```


(2) Brewer palettes:
```{r palettes, message=FALSE, warning=FALSE}
library(RColorBrewer)
colstx <- rev(brewer.pal(n = 9, "Spectral"))
colsindex <- rev(brewer.pal(n = 9, "RdYlBu"))
colsdelta <- brewer.pal(n = 9, "Reds")
colsbias <- brewer.pal(n = 9, "PiYG")
colssd <- brewer.pal(n = 9, "Blues")
```

***
NOTE: see also [2018_climate4R_example1.pdf](https://github.com/SantanderMetGroup/notebooks/blob/master/2018_climate4R_example1.pdf) for a better understanding of this document.

***

# Example 2: CORDEX Ensembles via the User Data Gateway

## Loading via the User Data Gateway

We define the domain of the Iberian Peninsula with the following bounding coordinates:

```{r boundary, message=FALSE, warning=FALSE}
lon <- c(-10, 5)
lat <- c(36, 44)
```

The UDG service requires (free) [registration](http://www.meteo.unican.es/udg-wiki) to accept the data policies of the different data providers. Once a valid user name and password have been issued, the authentication must be done within the R session before data loading with function `loginUDG`:

```{r loginudg, eval=FALSE}
loginUDG(username = "userUDG", password = "pswrdUDG")
```

If the data is to be loaded from the UDG, we can use function `UDG.datasets` to print the inventory of the available public and harmonized **UDG** datasets, where the name, type and url are specified (see Table 1 in the manuscript). Additionally, we can use the "name" (`UDG.datasets()["name"]`) of the desired dataset instead of passing the complete url to `loadGridData`. If this is the case, we do not need to create a dictionary, since the data is harmonized by default. Nevertheless, we still can use the complete url (`UDG.datasets()["url"]`) to access the data in its original form. 

For example, if we are interested in loading observations from the **E-OBS** dataset at 0.25 degrees resolution we can filter the names returned by `UDG.datasets` passing an appropriate pattern to the optional argument `pattern`:

```{r eobs, message=FALSE, warning=FALSE}
models <- UDG.datasets(pattern = "E-OBS.*0.25")
eobs <- models$name
```


Object `eobs` contains the name of the dataset and is passed to `loadGridData` for loading maximum temperature using the standard name "tasmax" (and without using the dictionary argument, see [2018_climate4R_example1.pdf](https://github.com/SantanderMetGroup/notebooks/blob/master/2018_climate4R_example1.pdf)):

```{r login, echo = FALSE, message=FALSE, warning=FALSE, cache=FALSE}
source("/media/maialen/work/WORK/creds")
```


```{r loadeobs, message=FALSE, warning=FALSE}
TX <- loadGridData(eobs,
                         var = "tasmax",
                         season = 1:12,
                         lonLim = lon,
                         latLim = lat,
                         years = 1971:2000)
```

To load historical CORDEX data for the Iberian Peninsula we also filtered the names in `UDG.datasets()` with an appropriate pattern:

```{r cordexens, message=FALSE, warning=FALSE}
models <- UDG.datasets(pattern = "CORDEX-EUR44.*historical")
ensemble.h <- models$name[1:6]
```

Unlike the first example of the paper (see [2018_climate4R_example1.pdf](https://github.com/SantanderMetGroup/notebooks/blob/master/2018_climate4R_example1.pdf)), here we considered 6 regional climate models (RCMs) from EURO-CORDEX, thus, object `ensemble.h` contains 6 names. Everything can be loaded in a single step by combining function `loadGridData` with `lapply`:

```{r loadhist, message=FALSE, warning=FALSE}
TXh.list <- lapply(ensemble.h, function(x) 
                              loadGridData(dataset = x,
                              var = "tasmax",
                              season = 1:12,
                              lonLim = lon,
                              latLim = lat,
                              years = 1971:2000))
```

We repeat the operation for the RCP8.5 scenario:

```{r loadrcp, message=FALSE, warning=FALSE}
ensemble.f <- UDG.datasets(pattern = "CORDEX-EUR44.*rcp85")$name[1:6]
TXf.list <- lapply(ensemble.f, function(x) 
                          loadGridData(dataset = x,
                          var = "tasmax",
                          season = 1:12,
                          lonLim = lon,
                          latLim = lat,
                          years = 2071:2100))
```

As a result, we obtain the following harmonized grids:

* `TX`: a single grid for E-OBS and reference period 1971-2000
* `TXh.list`: a list of 6 grids for the reference period 1971-2000 (historical scenario)
* `TXf.list`: a list of 6 grids for the future period 2071-2100 (RCP8.5 scenario).


## Working with multi-model ensembles

climate4R functionalities allow working with model ensembles through the `member` dimension. Thus, we can aggregate the list of `grid`s created before (`TXh.list` and `TXf.list`) along the `member` dimension to obtain a single `grid`. However, we need to check temporal and spatial consistency among the different models. In this case, a temporal inconsistency exists in our ensemble, since two of the models contain less calendar days. We can check this easily with function `getShape`:

```{r checkintersecttime, message=FALSE, warning=FALSE}
lapply(TXh.list, function(x) getShape(x))
```

Function `intersectGrid` performs time subsetting of a collection of grids, to the dates they have in common:

```{r intersecttime, message=FALSE, warning=FALSE}
# Temporal intersection
TXh.list <- intersectGrid(TXh.list, type = "temporal",
                                          which.return = 1:6) 
TXf.list <- intersectGrid(TXf.list, type = "temporal",
                                          which.return = 1:6)
```

If there is not spatial consistency, we can use `interpGrid` to interpolate all `grid`s to the same spatial structure. Despite this not being the case, here we perform interpolation of CORDEX data to the E-OBS spatial grid in order to apply a land-sea mask (by means of function `gridArithmetics`) and produce comparable map figures.

```{r interp, message=FALSE, warning=FALSE}
# Interpolation 
TXh.list <- lapply(TXh.list, function(x) interpGrid(x, getGrid(TX)))
TXf.list <- lapply(TXf.list, function(x) interpGrid(x, getGrid(TX)))
```

```{r applymask, message=FALSE, warning=FALSE}
# Create mask
m <- TX$Data[1,,]*0
mask.hist <- array(dim = c(getShape(TXh.list[[1]])["time"], dim(m)))
for (i in 1:dim(mask.hist)[1]) mask.hist[i,,] <- m 
mask.rcp <- array(dim = c(getShape(TXf.list[[1]])["time"], dim(m)))
for (i in 1:dim(mask.rcp)[1]) mask.rcp[i,,] <- m 

# Apply mask
TXh.list <- lapply(TXh.list, function(x) 
  gridArithmetics(x, mask.hist, operator = "+"))
TXf.list <- lapply(TXf.list, function(x) 
  gridArithmetics(x, mask.rcp, operator = "+"))
```

Finally, we crate the multi-member grids for the historical and RCP8.5 scenarios of CORDEX (objects `TXh.ens` and `TXf.ens`):

```{r multimember, message=FALSE, warning=FALSE}
# Create a multimember grid
TXh.ens <- bindGrid(TXh.list, dimension = "member")
TXf.ens <- bindGrid(TXf.list, dimension = "member")
```

Note that function `spatialPlot` recognizes a multi-member grid and displays a map for each member. For instance, next we plot the ensemble of the RCP8.5 scenario to generate Figure \ref{fig:fig8a} (Fig. 8(top) in the manuscript)):

```{r fig8a, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig8a}Maximum temperature (ºC) in Iberia for an ensemble of 6 CORDEX RCMs under the RCP8.5 scenario and for future period 2071-2100. Fig. 8(top) in the manuscript."}
spatialPlot(climatology(TXf.ens), at = seq(5, 33, 1), backdrop.theme = "countries",
            col.regions = colorRampPalette(colstx), layout = c(3, 2), as.table = TRUE)
```


In order to calculate the bias of each model with respect to the observations and generate Figure \ref{fig:fig8b} (Fig. 8(below) in the manuscript) we apply function `aggregateGrid`, `gridArithmetics`, `bindGrid` and `spatialPlot` as follows:

```{r calcbias, message=FALSE, warning=FALSE}
# Create a multimember grid
TXh.list.ann <- lapply(TXh.list, function(x) 
  aggregateGrid(x, aggr.y = list(FUN = "mean", na.rm = TRUE)))
TX.ann <- aggregateGrid(TX, aggr.y = list(FUN = "mean", na.rm = TRUE))
TXh.list.bias <- lapply(TXh.list.ann, function(x) 
  gridArithmetics(x, TX.ann, operator = "-"))
```

```{r null, echo=FALSE, message=FALSE, warning=FALSE}
tx_iberia.ann <- NULL
txHist_iberia.ann <- NULL
```


```{r multimemberbias, message=FALSE, warning=FALSE}
TXh.bias.ens <- bindGrid(TXh.list.bias, dimension = "member")
```

```{r fig8b, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig8b}Bias of the maximum temperature (ºC) in Iberia for an ensemble of 6 CORDEX RCMs w.r.t. E-OBS in the historical period 1971-2000. Fig. 8(below) in the manuscript."}
spatialPlot(climatology(TXh.bias.ens), at = seq(-10, 10, 1), backdrop.theme = "countries",
            col.regions = colorRampPalette(colsbias), layout = c(3, 2), as.table = TRUE)
```


We can use function `aggregateGrid` to for example calculate the multi-member mean and deviation of the ensemble:

```{r aggrmeansd, message=FALSE, warning=FALSE}
TXf.ens.mean <- aggregateGrid(climatology(TXf.ens), aggr.mem = list(FUN = mean, na.rm = TRUE))
TXf.ens.sd <- aggregateGrid(climatology(TXf.ens), aggr.mem = list(FUN = sd, na.rm = TRUE))
```

Next we generate the corresponding Figures \ref{fig:fig8c} and \ref{fig:fig8d} (Not shown in the manuscript).

```{r fig8c, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig8c} Ensemble mean of maximum temperature (ºC) in Iberia for  6 CORDEX RCMs under the RCP8.5 scenario and for future period 2071-2100. Not shown in the manuscript."}
spatialPlot(TXf.ens.mean, at = seq(5, 33, 1),
            col.regions = colorRampPalette(colstx), backdrop.theme = "countries")
```
```{r fig8d, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig8d} Ensemble standard deviation (sd) of maximum temperature (ºC) in Iberia for 6 CORDEX RCMs under the RCP8.5 scenario and for future period 2071-2100. Not shown in the manuscript."}
spatialPlot(TXf.ens.sd, at = seq(0, 6, .5),
            col.regions = colorRampPalette(colssd), backdrop.theme = "countries")
```


## ETCCDI index calculation (SU) from raw data

Next the raw SU index (summer days) is calculated for the whole ensemble and future period 2071-2100 (object `SUf.ens`), in a single line with function `climdexGrid`:

```{r climdexraw, message=FALSE, warning=FALSE}
SUf.ens <- climdexGrid(tx = TXf.ens, index.code = "SU")
```

We calculate the ensemble mean and deviation using function `aggregateGrid`:

```{r aggrmeansdSU, message=FALSE, warning=FALSE}
SUf.ens.mean <- aggregateGrid(climatology(SUf.ens), aggr.mem = list(FUN = mean, na.rm = TRUE))
SUf.ens.sd <- aggregateGrid(climatology(SUf.ens), aggr.mem = list(FUN = sd, na.rm = TRUE))
```

And finally, we plot the results to generate Figures \ref{fig:fig6top1} and \ref{fig:fig6top2} (Figs. 6(top, left) and 6(top, right) in the manuscript:

```{r fig6top1, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig6top1} Summer days in Iberia for the future period 2071-2100 computed from the original RCM daily maximum temperature data. The figure shows the ensemble mean. Fig. 6(top, left) in the manuscript."}
spatialPlot(SUf.ens.mean, backdrop.theme = "countries", 
            at = seq(0, 260, 10), col.regions = colorRampPalette(colsindex))
```
```{r fig6top2, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig6top2} Standard deviation of summer days in Iberia for the future period 2071-2100 computed from the original RCM daily maximum temperature data. The figure shows the spread of the ensemble. Fig. 6(top, right) in the manuscript."}
spatialPlot(SUf.ens.sd, set.max = 50, backdrop.theme = "countries", 
            at = seq(0, 50, 5), col.regions = colorRampPalette(colssd))
```


## Bias correction of the maximum temperature

Here we apply Empirical Quantile Mapping (EQM) for correcting the bias of each ensemble member (by default argument `join.members = FALSE`) by considering a moving correction window of 30 days to correct each 7-day time interval:

```{r biascorrectTX, message=FALSE, warning=FALSE}
TXf.ens.bc <- biasCorrection(TX, 
                              TXh.ens, 
                              TXf.ens, 
                              window = c(30, 7), 
                              extrapolation = "constant",
                              method = "eqm")
```

As done before, we calculate the ensemble mean and deviation of the maximum temperature with function `aggregateGrid`:

```{r aggrmeansedBC, message=FALSE, warning=FALSE}
TXf.ens.bc.mean <- aggregateGrid(TXf.ens.bc, aggr.mem = list(FUN = mean, na.rm = TRUE))
TXf.ens.bc.sd <- aggregateGrid(TXf.ens.bc, aggr.mem = list(FUN = sd, na.rm = TRUE))
```

We plot the results to generate Figures \ref{fig:fig8e} and \ref{fig:fig8f} (Not shown in the manuscript):

```{r fig8e, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig8e} Ensemble mean of bias corrected maximum temperature (ºC) in Iberia for  6 CORDEX RCMs under the RCP8.5 scenario and for future period 2071-2100. Not shown in the manuscript."}
spatialPlot(climatology(TXf.ens.bc.mean), backdrop.theme = "countries", 
            at = seq(5, 33, 1), col.regions = colorRampPalette(colstx))
```
```{r fig8f, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig8f} Ensemble standard deviation (sd) of bias corrected maximum temperature (ºC) in Iberia for 6 CORDEX RCMs under the RCP8.5 scenario and for future period 2071-2100. Not shown in the manuscript."}
spatialPlot(climatology(TXf.ens.bc.sd), backdrop.theme = "countries", 
            at = seq(0, 6, .5), col.regions = colorRampPalette(colssd))
```

## ETCCDI index calculation (SU) from bias corrected data

Next we apply function `climdexGrid` over the bias corrected maximum temperature (object `TXf.ens.bc`) to obtain the corrected SU index for the future period.

```{r climdexBC, message=FALSE, warning=FALSE}
SUf.ens.bc <- climdexGrid(tx = TXf.ens.bc, index.code = "SU")
```

Again, the ensemble mean and deviation is calculated:

```{r aggrmeansdSUbc, message=FALSE, warning=FALSE}
SUf.ens.bc.mean <- aggregateGrid(SUf.ens.bc, aggr.mem = list(FUN = mean, na.rm = TRUE))
SUf.ens.bc.sd <- aggregateGrid(SUf.ens.bc, aggr.mem = list(FUN = sd, na.rm = TRUE))
```

By plotting these results Figures \ref{fig:fig6below1} and \ref{fig:fig6below2} are generated (Figs. 6(below, left) and 6(below, right) in the manuscript):

```{r fig6below1, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig6below1} Summer days in Iberia for the future period 2071-2100 computed from the bias corrected RCM daily maximum temperature data. The figure shows the ensemble mean. Fig. 6(below, left) in the manuscript."}
spatialPlot(climatology(SUf.ens.bc.mean), backdrop.theme = "countries", 
            at = seq(0, 230, 10), col.regions = colorRampPalette(colsindex))
```
```{r fig6below2, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig6below2} Standard deviation of summer days in Iberia for the future period 2071-2100 computed from bias corrected RCM daily maximum temperature data. The figure shows the spread of the ensemble. Fig. 6(below, right) in the manuscript."}
spatialPlot(climatology(SUf.ens.bc.sd), backdrop.theme = "countries", 
            at = seq(0, 50, 5), col.regions = colorRampPalette(colssd))
```



In order to generate Figure \ref{fig:fig7} (Fig. 7 in the manuscript), we next calculate the SU index from observational data (`SU`) and the historical model ensemble (`SUh.ens`):

```{r climdexSUall, message=FALSE, warning=FALSE}
#A single location
SU <- climdexGrid(tx = TX, index.code = "SU")
SUh.ens <- climdexGrid(tx = TXh.ens, index.code = "SU")
```

For visualization purposes, the 1rst memeber of CORDEX data is extracted and everything (E-OBS, the ensemble of CORDEX, and the 1st member of CORDEX) is listed in object `SU.Z` for plotting:

```{r preparetemporalplot, message=FALSE, warning=FALSE}
SU.Z.cdx <- list(SUh.ens, SUf.ens, SUf.ens.bc)
SU.Z.1m <- lapply(SU.Z.cdx, function(x)
                        subsetGrid(x, members = 1))
SU.Z <- c(list(SU), SU.Z.cdx, SU.Z.1m)
```

Finally `temporalPlot` is used to generate Figure  \ref{fig:fig7} (Fig. 7 in the manuscript). The nearest point to Zaragoza (Spain) is selected via arguments `latLim` and `lonLim`:

```{r fig7, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig7} Annual summer days time series for a single gridbox (the one closest to Zaragoza, Spain) computed form (red) the original RCM daily maximum temperature data, and (blue) daily maximum temperature bias corrected data using E-OBS (black). When it comes to CORDEX data, continuous lines correspond to the ensemble mean and the shadowed area to the range (uncertainty). Dashed lines correspond to the 1st member of the ensemble. Fig. 7 in the manuscript."}
cols  <- c("black", rep(c("red", "red", "blue"), 2))
temporalPlot(SU.Z,
             latLim = 41.64, 
             lonLim = -0.89,
             cols = cols,
             lty = rep(c(1,3), each = 4),
             lwd = 0.8,
             xyplot.custom = 
               list(ylim = c(70, 220), ylab = "",
                                  key = list(space = "top", 
                                             lines = list(lty = c(rep(1,4), 3), 
                                                          col = c(cols[1:4], cols[1]),
                                                          lwd = 0.8), 
                                             text = list(c("E-OBS", "CDX_hist", 
                                                    "CDX_rcp85", "CDX_rcp85_corrected",
                                                           "1st_member"),
                                                         cex = .7), 
                                             columns = 3, rows = 3)))
```

```{r endpart1, warning=FALSE, message=FALSE, echo=FALSE}
rm(list=c("SU.Z", "SU","SUh.ens","SU.Z.eobs","SU.Z.cdx","SU.Z.1m","SUf.ens.bc.mean","TXh.ens","TXf.ens","TXf.ens.bc"))
```

## Calculation of the ETCCDI index CDD (Consecutive Dry Days)

This section replicates the procedures shown top but for calculating a precipitation based index, i.e. CDD (consecutive dry days). Function `climdexShow` displays an overview of the available ETCDDI indices:

```{r showindexes, eval= FALSE, message=FALSE, warning=FALSE}
climdexShow()[,1:6]
```


### Data loading and multi-member ensemble building

In order to calculate the CDD index, we need to load daily precipitation data. Therefore, we repeat the loading operation (function `loadGridData`) shown in the previous section (this time `var = "pr"`) for observational data (E-OBS) and historical and RCP8.5 projection data (6 models from CORDEX).

```{r authenticationhide2, echo = FALSE, message=FALSE, warning=FALSE, cache=FALSE}
source("/media/maialen/work/WORK/creds")
```



```{r eobsurlupdate, message=FALSE, warning=FALSE, echo=FALSE, eval=FALSE}
eobs.pr<-"http://opendap.knmi.nl/knmi/thredds/dodsC/e-obs_0.25regular/rr_0.25deg_reg_v16.0.nc"
```


```{r dicupdate, message=FALSE, warning=FALSE, eval=FALSE, echo=FALSE}
write(c("identifier,short_name,time_step,lower_time_bound,upper_time_bound,cell_method,
        offset,scale,deaccum,derived,interface",
             "tasmax,tx,24h,0,24,max,0,1,0,0,",
        "pr,rr,24h,0,24,max,0,1,0,0,"), "dicEOBS.dic")
```

```{r loadpreobs, message=FALSE, warning=FALSE}
pr <- loadGridData(eobs,
                         var = "pr",
                         season = 1:12,
                         lonLim = lon,
                         latLim = lat,
                         years = 1971:2000)
```

```{r dicupdateCDX, message=FALSE, warning=FALSE, eval=FALSE, echo=FALSE}
write(c("identifier,short_name,time_step,lower_time_bound,upper_time_bound,cell_method,
        offset,scale,deaccum,derived,interface",
             "tasmax,tasmax,24h,0,24,max,-273.15,1,0,0,",
             "pr,pr,24h,0,24,max,0,86400,0,0,"), "dicCDX.dic")
```

```{r loadprHistRcp, message=FALSE, warning=FALSE}
prh <- lapply(ensemble.h, function(x) 
                              loadGridData(dataset = x,
                              var = "pr",
                              season = 1:12,
                              lonLim = lon,
                              latLim = lat,
                              years = 1971:2000))
prf <- lapply(ensemble.f, function(x) 
                          loadGridData(dataset = x,
                          var = "pr",
                          season = 1:12,
                          lonLim = lon,
                          latLim = lat,
                          years = 2071:2100))
prh.t <- intersectGrid(prh, type = "temporal", which.return = 1:6)
prf.t <- intersectGrid(prf, type = "temporal", which.return = 1:6)
prh.ens <- bindGrid(prh.t, dimension = "member")
prf.ens <- bindGrid(prf.t, dimension = "member")
```

Next the CDD index is calculated (function `climdexGrid`), and regridded (function `interpGrid`) to the spatial structure given by the observational data (`getGrid(pr)`). Then, the ensemble mean (object `CDDf.ens.mean`) and deviation (object `CDDf.ens.sd`) are calculated with `aggregateGrid`.

```{r climdexinterpaggrmask, warning=FALSE, message=FALSE}
CDDf.ens <- climdexGrid(pr = prf.ens, index.code = "CDD")
CDDf.ens.interp <- interpGrid(CDDf.ens, getGrid(pr))
CDDf.ens.mean <- aggregateGrid(CDDf.ens.interp, 
                               aggr.mem = list(FUN = "mean", na.rm = TRUE))
CDDf.ens.sd <- aggregateGrid(CDDf.ens.interp, 
                             aggr.mem = list(FUN = "sd", na.rm = TRUE))
```

We also create and apply the land-sea mask:

```{r cddapplymask, warning=FALSE, message=FALSE}
m <- pr$Data[1,,]*0
mask.cdd <- array(dim = c(getShape(CDDf.ens.mean)["time"], dim(m)))
for (i in 1:dim(mask.cdd)[1]) mask.cdd[i,,] <- m 
CDDf.ens.mean <- gridArithmetics(CDDf.ens.mean, mask.cdd, operator = "+")
CDDf.ens.sd <- gridArithmetics(CDDf.ens.sd, mask.cdd, operator = "+")
```

Finally, we plot the result with function `spatialPlot` to generate Figures \ref{fig:fig9} and \ref{fig:fig10} (not shown in the manuscript).

```{r fig9, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig9} Consecutive dry days (CDD) in Iberia for the future period 2071-2100 computed from raw RCM daily precipitation data. The figure shows the the ensemble mean. Not shown in the manuscript."}
spatialPlot(climatology(CDDf.ens.mean), backdrop.theme = "countries", at = seq(0, 60, 5), 
            col.regions = colorRampPalette(colsindex))
```
```{r fig10, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig10} Standard deviation of consecutive dry days (CDD) in Iberia for the future period 2071-2100 computed from raw RCM daily precipitation data. The figure shows the spread of the ensemble. Not shown in the manuscript."}
spatialPlot(climatology(CDDf.ens.sd), backdrop.theme = "countries", at = seq(0, 60, 5), 
            col.regions = colorRampPalette(colssd))
```

### Bias correction of the precipitation

Here, we use again the EQM method for bias correcting precipitation data to subsequently calculate the corrected CDD index (objects `CDDf.ens.bc`, `CDDf.ens.bc.mean`, `CDDf.ens.bc.sd`).

```{r biascorrclimdexmemaggr, warning=FALSE, message=FALSE}
prf.ens.bc <- biasCorrection(y = pr, 
                                    x = prh.ens, 
                                    newdata = prf.ens,
                                    precipitation = TRUE,
                                    window = c(30, 7), 
                                    extrapolation = "constant",
                                    method = "eqm",
                                    wet.threshold = 0.1)
CDDf.ens.bc <- climdexGrid(pr = prf.ens.bc, index.code = "CDD")
CDDf.ens.bc.mean <- aggregateGrid(CDDf.ens.bc, 
                                  aggr.mem = list(FUN = "mean", na.rm = TRUE))
CDDf.ens.bc.sd <- aggregateGrid(CDDf.ens.bc, 
                                aggr.mem = list(FUN = "sd", na.rm = TRUE))
```

This results are plotted to generate Figures \ref{fig:fig11} and \ref{fig:fig12} (not shown in the manuscript):

```{r fig11, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig11} Consecutive dry days (CDD) in Iberia for the future period 2071-2100 computed from bias corrected RCM daily precipitation data. The figure shows the the ensemble mean. Not shown in the manuscript."}
spatialPlot(climatology(CDDf.ens.bc.mean), backdrop.theme = "countries", 
            at = seq(0, 60, 5), col.regions = colorRampPalette(colsindex))
```

```{r fig12, message=FALSE, warning=FALSE, fig.cap="\\label{fig:fig12} Standard deviation of consecutive dry days (CDD) in Iberia for the future period 2071-2100 computed from bias corrected RCM daily precipitation data. The figure shows the spread of the ensemble. Not shown in the manuscript."}
spatialPlot(climatology(CDDf.ens.bc.sd), backdrop.theme = "countries", 
            at = seq(0, 60, 5), col.regions = colorRampPalette(colssd))
```

In this case, unlike the SU index, the uncertainty is not reduced before and after bias correcting the original variable (compare Figs. \ref{fig:fig10} and \ref{fig:fig12}). Nevertheless the pattern of the index is significantly modified (compare Figs. \ref{fig:fig9} and \ref{fig:fig11}).

# Other available material

* [2018_climate4R_example1.pdf](https://github.com/SantanderMetGroup/notebooks/blob/devel/2018_climate4R_example1.pdf) contains the full code for **Example 1** of the paper `climate4R: An Ecosystem of R packages for Climate Data Access, Post-processing and Bias Correction'.

* Find more worked examples on the utilization of climate4R packages in their respective GitHub **wiki**-s at [https://github.com/SantanderMetGroup](https://github.com/SantanderMetGroup):

    + [loadeR: https://github.com/SantanderMetGroup/loadeR/wiki](https://github.com/SantanderMetGroup/loadeR/wiki)
    + [transformeR: https://github.com/SantanderMetGroup/transformeR/wiki](https://github.com/SantanderMetGroup/transformeR/wiki)
    + [downscaleR: https://github.com/SantanderMetGroup/downscaleR/wiki](https://github.com/SantanderMetGroup/downscaleR/wiki)
    + [visualizeR: https://github.com/SantanderMetGroup/visualizeR/wiki](https://github.com/SantanderMetGroup/visualizeR/wiki)