2018_drought4R_NCOMMS.Rmd

---
title: "Bias-corrected future SPEI projections with climate4R"
subtitle: "A worked example using Euro-CORDEX projections"
author: "J. Bedia"
date: "`r Sys.Date()`"
encoding: "UTF8"
urlcolor: blue
output:
  html_document:
    fig_caption: yes
    highlight: pygments
    number_sections: yes
    theme: readable
    toc: yes
    toc_float: yes
  pdf_document:
    highlight: pygments  
    toc: yes
    pandoc_args: [
      "--number-sections",
      "--number-offset=0"]
    latex_engine: xelatex
documentclass: article
abstract: This is a worked example on how to obtain bias-corrected future drought index projections using the climate4R framework for climate data access and analysis. The example is based on RCM data from EURO-CORDEX (Jacob _et al._ 2014), using the CRU-TS4 gridded observations (Harris and Jones 2017) as reference. The example describes the main steps to obtain the future SPEI projections, as used in related paper by Turco _et al._ (2018). This notebook is also available in the following link <http://www.meteo.unican.es/work/climate4r/drought4R/drought4R_notebook.html>
---

\fontfamily{cmr}
\fontsize{11}{22}
\selectfont


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      highlight = TRUE,
                      message = FALSE,
                      fig.align = "center",
                      tidy = FALSE,
                      fig.width = 7,
                      cache = TRUE)
```


# Introduction and required R packages

## The *climate4R* framework

*climate4R* (Iturbide _et al._, submitted) is a bundle of R packages for transparent climate data access, post-processing (including bias correction and downscaling) and visualization. *climate4R* builds on two main data structures (grid and station, including metadata) to deal with gridded and point data from observations, reanalysis, seasonal forecasts and climate projections. Thus, it considers ensemble members as a basic dimension of the data structures. *climate4R* exploits NetCDF-Java and allows accessing local and remote (OPeNDAP) data sources, including the User Data Gateway (UDG), a THREDDS-based service for a variety of widely used datasets (e.g. reanalysis, CMIP5, CORDEX). This provides a unique comprehensive framework for end-to-end sectoral applications, favouring the reproducibility of scientific outcomes (see e.g. Cofiño _et al._ 2018).

Compatibility with some external packages has been achieved by appropriate two-way bridging functions, enhancing *climate4R* with sector-specific functionalities: model validation, extreme climate indices, etc. In this notebook we illustrate the extension of *climate4R* to calculate drought indices. In particular, here we use [`drought4R`](https://github.com/SantanderMetGroup/drought4R), a wrapping package tailoring the SPI, SPEI and evapotranspiration routines for monthly data implemented in the R package `SPEI` (Beguería and Serrano 2017) to the needs of *climate4R* users handling large climate model datasets (multi-decadal projections, seasonal predictions...) and gridded observations or reanalysis.

## Package availability

All the required packages are publicly available through the [Santander MetGroup GitHub repository](https://github.com/SantanderMetGroup). Further details and access points can be obtained in the [climate4R site](http://www.meteo.unican.es/climate4r)

*climate4R* is formed by four main packages `loadeR` (Cofiño _et al._ 2018), `transformeR` (Iturbide _et al._, submitted), `downscaleR` (Bedia _et al._ 2018) and `visualizeR` (Frías _et al._ 2018), performing the tasks of data loading, manipulation and transformation, bias correction and downscaling and data visualization respectively. These are next loaded:

Next, installation procedure is shown. The use of the `install_github` function from package `devtools` (Wickham _et al._ 2018) is recommended. Note that the installation paths point to the stable versions used to produce this notebook. The full reproducibility of the results can’t be guaranteed with other versions.

```{r,eval=FALSE}
devtools::install_github("SantanderMetGroup/loadeR.java",
                         "SantanderMetGroup/loadeR@v1.4.5",
                         "SantanderMetGroup/transformeR@v1.4.3",
                         "SantanderMetGroup/downscaleR@v3.0.1",
                         "SantanderMetGroup/visualizeR@v1.2.1")
```


In addition, two "add-on" packages seamlessly integrated in *climate4R* are used for SPEI calculation (`drought4R`) and unit conversion (`convertR`). The latter relies on the Unidata's [UDUNITS](https://www.unidata.ucar.edu/software/udunits) software libraries for handling physical quantity units. Specific installation instructions are available in their respective GitHub sites in case any problems arise:

```{r}
devtools::install_github("SantanderMetGroup/drought4R@v0.1.0",
                         "SantanderMetGroup/convertR@v0.1.2")
```

Finally, the pipe operator (`%>%`) in package `magrittr` (Bache and Wickham 2014) will be used throughout the examples to concatenate command calls in a convenient way (the package is available in CRAN and can be installed using the R base function `install.packages`):

```{r}
library(magrittr)
```


# Data access

## Baseline observations (CRU 4.0, period 1971-2000) 

The CRU TS gridded observations (Harris and Jones 2017) can be downloaded from the [data access site](http://catalogue.ceda.ac.uk/uuid/58a8802721c94c66ae45c3baa4d814d0). For brevity, the process to prepare these raw data to be accessible using `loadeR` is skipped (the interested reader has worked examples in other sources, see e.g. Iturbide _et al._ (in press), or the examples in the [loadeR's wiki](https://github.com/SantanderMetGroup/loadeR/wiki). 

Thus, the data have been already prepared in a convenient format to be used in this example, and are available on-line. The following function can be used for a straightforward download into the current R session:

```{r}
my_readRDS <- function(file.url) {
    tmpfile <- tempfile()
    utils::download.file(url = file.url, destfile = tmpfile)
    readRDS(tmpfile)
}
```



```{r,echo=FALSE}
# load("~/workspace/COLLABORATIONS/2017_CORDEX_SPEI/data/CRU_tp_monthly_1971_2011.Rdata", verbose = TRUE)
# pr.cru <- subsetGrid(tp.cru, years = 1971:2001)
# # ref.grid <- getGrid(pr.cru)
# load("~/workspace/COLLABORATIONS/2017_CORDEX_SPEI/data/CRU_tasmax_monthly_1971_2011.Rdata", verbose = TRUE)
#  tasmax.cru <- subsetGrid(tasmax.cru, years = 1971:2001)
#  load("~/workspace/COLLABORATIONS/2017_CORDEX_SPEI/data/CRU_tasmin_monthly_1971_2011.Rdata", verbose = TRUE)
#  tasmin.cru <- subsetGrid(tasmin.cru, years = 1971:2001)
# saveRDS(pr.cru, file = "data/CRU_precip.rds", compress = "xz")
# saveRDS(tasmin.cru, file = "data/CRU_tasmin.rds", compress = "xz")
# saveRDS(tasmax.cru, file = "data/CRU_tasmax.rds", compress = "xz")
pr.cru <- readRDS("data/CRU_precip.rds")
tasmin.cru <- readRDS("data/CRU_tasmin.rds")
tasmax.cru <- readRDS("data/CRU_tasmax.rds")
```

We next load the three observational datasets used:

```{r,eval=FALSE}
tasmin.cru <- my_readRDS("http://www.meteo.unican.es/work/UDG/drought4R/CRU_tasmin.rds")
tasmax.cru <- my_readRDS("http://www.meteo.unican.es/work/UDG/drought4R/CRU_tasmax.rds")
pr.cru <- my_readRDS("http://www.meteo.unican.es/work/UDG/drought4R/CRU_tpr.rds")
```




The reference grid of the CRU data is retained in order to interpolate the RCM data afterwards:

```{r,message=TRUE}
library(transformeR)
ref.grid <- getGrid(pr.cru)
```

Next, the climatology of one of the reference observations just loaded is shown:

```{r, message=TRUE}
library(visualizeR)
climatology(tasmin.cru) %>% spatialPlot(backdrop.theme = "countries",
                                        rev.colors = TRUE,
                                        main = "CRU-TS4.0 tasmin annual climatology (1971-2000)")
```


## Accessing Euro-CORDEX data from the User Data Gateway (and bias correction)

In this section we illustrate the use of *climate4R* as the data access layer of the UDG. The added value of using the `loadeR` package tools for data access directly pointing to the UDG server are here exemplified.


### Historical experiment data (1971-2000)


The UDG service requires (free) registration to accept the data policies of the different data providers (http://www.meteo.unican.es/udg-wiki). Prior to data access, authentication with valid UDG credentials is required for the current R session in order to access the UDG. Once a valid user name and password have been issued, the authentication can be done in one step within the R session using the `loginUDG` function from `loadeR`:


```{r,message=TRUE}
library(loadeR)
```


```{r,eval=FALSE}
loginUDG(username = "jdoe", password = "****")
```

To get a quick overview of the available datasets, the function `UDG.datasets` prints a summary. For all datasets included in the UDG and listed by the function, the name of the target dataset can be used as a valid entry for the argument dataset in `loadGridData`, instead of the full URL. Next, data from the CORDEX historical scenario available at UDG are loaded. In this example, for simplicity and faster calculations, the 0.44 regular degree grid will be used (note that the 0.11º simulations are also available at UDG).

```{r}
models <- UDG.datasets()
```

There are many available datasets. Here, for illustration we will use one of the GCM/RCM couplings of Euro-CORDEX used by Turco _et al._ 2018, in particular, the EC-Earth/RACMO22. However, for brevity and minimal computational cost, the 0.44 degree horizontal resolution will be used, instead of the 11 degree used in the paper. Note that the latter would be accessed just by replacing the "EUR44"" prefix by "EUR11". Using pattern matching we identify the dataset:

```{r}
grep("EUR44.*EC-EARTH.*RACMO22", models$name, value = TRUE)
```

The following call to the function `loadGridData` retrieves the historical simulation considering all the months of the year (`season = 1:12`, is omitted because it is the default), minimum temperature (`var = "tasmin"`), considering a Euro-Mediterranean spatial domain (`lonLim`and `latLim`) for the period 1971-2000 (argument `years`). Furthermore, monthly mean data is requested, via the corresponding temporal aggregation arguments. More details in `help("loadGridData", package = "loadeR")`


```{r}
lonLim <- c(-10, 30)
latLim <- c(34, 48)
```

```{r,eval=FALSE}
tasmin.hist <- loadGridData(dataset = "CORDEX-EUR44_EC-EARTH_r1i1p1_historical_RACMO22E_v1",
                            var = "tasmin",
                            lonLim = lonLim,
                            latLim = latLim,
                            years = 1971:2000,
                            time = "DD",
                            aggr.d = "mean",
                            aggr.m = "mean")
```

Next, we ensure that the data are in degress Celsius (model data are originally stored in Kelvin in most cases). To this aim, the function `udConvertGrid` from package `convertR` is used:

```{r,eval=FALSE}
library(convertR)
tasmin.hist <- udConvertGrid(tasmin.hist, new.units = "degC")
```

Finally, the data are regridded to the reference regular grid from the CRU (the original one is a rotated grid). This is done using the `interpGrid` function from package `transformeR`:

```{r,eval=FALSE}
tasmin.hist <- interpGrid(tasmin.hist, new.coordinates = ref.grid)
```

For brevity, in the following the above steps (data loading + unit conversion + regridding) will be concatenated using the `%>%` operator in the same call. 


```{r,echo=FALSE}
tasmin.hist <- readRDS(file = "data/tasmin_EC-EARTH_historical_r1i1p1_RACMO22E.rds")
```

The function `biasCorrection` of the package `downscaleR` allows applying a number of standard bias correction techniques within the *climate4R* framework. In particular, when dealing with monthly data, the common bias correction technique is the (additive and/or multiplicative) local scaling method. It is next applied to the historical data, using as reference the CRU observations:

```{r,message=TRUE}
library(downscaleR)
```


```{r}
tasmin.hist.corr <- biasCorrection(tasmin.hist,
                                   y = tasmin.cru,
                                   method = "scaling",
                                   scaling.type = "additive") %>% redim(drop = TRUE)
```

The same steps are next undertaken with maximum temperature (`var = "tasmax"`):


```{r,eval=FALSE}
tasmax.hist <- loadGridData(dataset = "CORDEX-EUR44_EC-EARTH_r1i1p1_historical_RACMO22E_v1",
                            var = "tasmax",
                            dictionary = dic, 
                            lonLim = lonLim,
                            latLim = latLim,
                            years = 1971:2000,
                            time = "DD",
                            aggr.d = "mean",
                            aggr.m = "mean") %>% udConvertGrid(new.units = "degC") %>% interpGrid(new.coordinates = ref.grid)
```

```{r,echo=FALSE}
tasmax.hist <- readRDS(file = "data/tasmax_EC-EARTH_historical_r1i1p1_RACMO22E.rds")
```



```{r}
tasmax.hist.corr <- biasCorrection(tasmax.hist,
                                   y = tasmax.cru,
                                   method = "scaling",
                                   scaling.type = "additive") %>% redim(drop = TRUE)
```

and precipitation (`var = "pr"`). Note that the monthly aggregation function is now set to `aggr.m = "sum"` (i.e., total accumulated monthly precipitation), as opposite to `"mean"` monthly temperature:

```{r,eval=FALSE}
pr.hist <- loadGridData(dataset = "CORDEX-EUR44_EC-EARTH_r1i1p1_historical_RACMO22E_v1",
                        var = "pr",
                        lonLim = lonLim,
                        latLim = latLim,
                        years = 1971:2000,
                        time = "DD",
                        aggr.d = "sum",
                        aggr.m = "sum") %>% udConvertGrid(new.units = "degC") %>% interpGrid(new.coordinates = ref.grid)
```

```{r,echo=FALSE}
pr.hist <- readRDS(file = "data/pr_EC-EARTH_historical_r1i1p1_RACMO22E.rds")
```


```{r}
pr.hist.corr <- biasCorrection(pr.hist,
                               y = pr.cru,
                               method = "scaling",
                               scaling.type = "multiplicative") %>% redim(drop = TRUE)
```


### RCP 8.5 experiment data

In a similar vein, the future (`years = 2010:2100`) data from the RCP 8.5 are loaded and transformed. Note that the value of the argument `dataset` is changed to point to the RCP 8.5 experiment dataset (`dataset = "CORDEX-EUR44_EC-EARTH_r1i1p1_rcp85_RACMO22E_v1"`)

```{r,eval=FALSE}
rcp85 <- "CORDEX-EUR44_EC-EARTH_r1i1p1_rcp85_RACMO22E_v1" 
tasmin.85 <- loadGridData(dataset = rcp85,
                          var = "tasmin",
                          lonLim = lonLim,
                          latLim = latLim,
                          years = 2010:2100,
                          time = "DD",
                          aggr.d = "mean",
                          aggr.m = "mean") %>% udConvertGrid(new.units = "degC") %>% interpGrid(new.coordinates = ref.grid)
```

```{r,echo=FALSE}
tasmin.85 <- readRDS(file = "data/tasmin_EC-EARTH_rcp85_r1i1p1_RACMO22E.rds")
```

Note that for the future projection correction, the calibration is undertaken considering the parameters estimated from the training period (1971-2000), using the historical simulation (`x = tasmin.hist`), and then applied to the new future slice of RCP 8.5 (passed via the `newdata` argument):

```{r}
tasmin.85.corr <- biasCorrection(y = tasmin.cru,
                                 x = tasmin.hist,
                                 newdata = tasmin.85,
                                 method = "scaling",
                                 scaling.type = "additive") %>% redim(drop = TRUE)
```

The same is done with maximum temperature:

```{r,eval=FALSE}
tasmax.85 <- loadGridData(dataset = rcp85,
                          var = "tasmax", 
                          lonLim = lonLim,
                          latLim = latLim,
                          years = 2010:2100,
                          time = "DD",
                          aggr.d = "mean",
                          aggr.m = "mean") %>% udConvertGrid(new.units = "degC") %>% interpGrid(new.coordinates = ref.grid)
```

```{r,echo=FALSE}
tasmax.85 <- readRDS(file = "data/tasmax_EC-EARTH_rcp85_r1i1p1_RACMO22E.rds")
```

```{r}
tasmax.85.corr <- biasCorrection(y = tasmax.cru,
                                 x = tasmax.hist,
                                 newdata = tasmax.85,
                                 method = "scaling",
                                 scaling.type = "additive") %>% redim(drop = TRUE)
```

and precipitation: 

```{r,eval=FALSE}
pr.85 <- loadGridData(dataset = rcp85,
                      var = "pr",
                      dictionary = dic, 
                      lonLim = lonLim,
                      latLim = latLim,
                      years = 2010:2100,
                      time = "DD",
                      aggr.d = "sum",
                      aggr.m = "sum") %>% interpGrid(new.coordinates = ref.grid)
```

```{r,echo=FALSE}
pr.85 <- readRDS(file = "data/pr_EC-EARTH_rcp85_r1i1p1_RACMO22E.rds")
```


```{r}
pr.85.corr <- biasCorrection(y = pr.cru, x = pr.hist, newdata = pr.85, method = "scaling",
                             scaling.type = "multiplicative") %>% redim(drop = TRUE)
```

# Calculatig future SPEI projections

Once both the historical and the RCP 8.5 projections have been bias-corrected using the CRU observations as reference, these can be joined along time into a single object, so a continuous time series for SPEI can be computed (however, note that there is a gap between the end of the reference period in 2001 and the start of transient period in 2010. This will be transparently handled by the involver *climate4R* functions involved). This is achieved by the function `bindGrid` from package `transformeR`:

```{r}
tx <- bindGrid(tasmax.hist.corr, tasmax.85.corr, dimension = "time") %>% redim(drop = TRUE)
tn <- bindGrid(tasmin.hist.corr, tasmin.85.corr, dimension = "time") %>% redim(drop = TRUE)
pr <- bindGrid(pr.hist.corr, pr.85.corr, dimension = "time") %>% redim(drop = TRUE)
```

Potential Evapotranspiration need to be calculated prior to computing SPEI. The function `petGrid` of the `drought4R` package is a wrapper of the `hargreaves` function from package `SPEI`, implementing the Hargreaves method for estimation of Potential Evapo-transpiration (see `help(hargreaves, package = "SPEI")` for further details):

```{r}
library(drought4R)
pet.har.85 <- petGrid(tasmin = tn, tasmax = tx, pr = pr, method = "hargreaves")
```

Finally, `speiGrid` is used, a wrapper to the `spei` function from package `SPEI`. Note that the index is computed considering the historical reference period (1971-2000), although the SPEI time series is calculated afterwards over the transient period.

This is done by passing the function the argument `ref.start` and `ref.end`, providing the temporal extent for calibrating the index:

```{r}
spei.rcp85 <- speiGrid(et0.grid = pet.har.85,
                       pr.grid = pr,
                       scale = 12,
                       ref.start = c(1971, 1),
                       ref.end = c(2000, 12),
                       na.rm = TRUE)
```

In the next plot, a time series for the nearest grid point to Madrid (Spain) shows the RCP 8.5 SPEI projections obtained. Note the use of `subsetGrid` to remove the historical period before plotting:

```{r}
m <- subsetGrid(spei.rcp85, lonLim = -3.43, latLim = 40.23, years = 2010:2100)
plot(ts(m$Data, start = c(2010,1), frequency = 12),
     main = "SPEI-12 (Hargreaves) RCP 8.5 Projection\nCalibration period: 1971-2000",
     ylab = "SPEI-12")
grid()
abline(h = 0)
mtext("Madrid, 3.43W / 40.23N")
```

# References

* Bache, S.M. and Wickham, H. (2014). magrittr: A Forward-Pipe Operator for R. R package version 1.5. https://CRAN.R-project.org/package=magrittr

* Bedia, J., Gutiérrez, J.M., Herrera, S., Iturbide, M., Manzanas, R. and Medina, J.B. (2018). downscaleR: An R package for bias correction and statistical downscaling. R package version 3.0.2. https://github.com/SantanderMetGroup/downscaleR/wiki

* Beguería, S. and Vicente-Serrano, S.M., 2017. SPEI: Calculation of the Standardised Precipitation-Evapotranspiration Index. R package version 1.7. https://CRAN.R-project.org/package=SPEI

* Cofiño A, Bedia J, Iturbide M, Vega M, Herrera S, Fernandez J, Frias M, Manzanas R and Gutierrez J, 2018. The
ECOMS User Data Gateway: Towards seasonal forecast data provision and research reproducibility in the era of
Climate Services. Climate Services. http://doi.org/10.1016/j.cliser.2017.07.001

* Frias, M.D., Iturbide, M., Manzanas, R., Bedia, J., Fernández, J., Herrera, S., Cofiño, A.S., Gutiérrez, J.M., 2018. An R package to visualize and communicate uncertainty in seasonal climate prediction. Environmental Modelling & Software 99, 101–110. https://doi.org/10.1016/j.envsoft.2017.09.008

* Harris, I.C., Jones, P.D. (2017): CRU TS4.01: Climatic Research Unit (CRU) Time-Series (TS) version 4.01 of high-resolution gridded data of month-by-month variation in climate (Jan. 1901- Dec. 2016). University of East Anglia Climatic Research Unit; Centre for Environmental Data Analysis, 04 December 2017. doi:10.5285/58a8802721c94c66ae45c3baa4d814d0. http://dx.doi.org/10.5285/58a8802721c94c66ae45c3baa4d814d0

* Iturbide, M., Bedia, J., Herrera, S., Bano-Medina, J., Fernández, J., Frı́as, M., Manzanas, R., San-Martı́n, D.,
Cimadevilla, E., Cofiño, A., Gutiérrez, J., _submitted_. climate4R: An R-based framework for Climate Data Access,
Post-processing and Bias Correction. Submitted to Environmental Modelling and Software.

* Jacob, D. _et al._, 2014. EURO-CORDEX: new high-resolution climate change projections for European impact research. doi:10.1007/s10113-013-0499-2.

* Turco, M., Cánovas, J.J.R., Bedia, J., Jerez, S., Montávez, J.P., Llasat, M.C. and Provenzale, A., 2018. Exacerbated fires in Mediterranean Europe due to anthropogenic warming projected with non-stationary climate-fire models. Nature Communications, in press


# Session information

```{r}
sessionInfo() %>% print()
```