2020_Lamb_ClimDyn.Rmd

---
title: 'Circulation Weather Typing with `climate4R`: GCM evaluation with Reanalysis using Lamb Weather Types'
author: J. A. Fernández, A. Casanueva, J. Bedia \& J. Fernández
date: '`r Sys.Date()`'
output:
  html_document:
    fig_caption: yes
    highlight: pygments
    number_sections: yes
    theme: readable
    toc: yes
    toc_float: yes
  pdf_document:
    fig_caption: yes
    highlight: pygments
    latex_engine: pdflatex
    pandoc_args:
    - --number-sections
    - --number-offset=0
    toc: yes
encoding: UTF8
documentclass: article
subtitle: Paper notebook - submitted to Climate Dynamics
abstract: This is an example notebook illustrating the calculations undertaken in the paper. It is not intended to provide full reproducibility of the results, but a sample on how to achieve this using the [climate4R framework](https://github.com/SantanderMetGroup/climate4R), used in the paper. To this aim, we provide an example using public datasets only, namely the NCEP-NCAR Reanalisys1 (NNRP) and the CMIP5 simulations of the EC-EARTH model, considering the historical and RCP8.5 experiments and the ensemble member _r12i1pi_.
urlcolor: blue
---

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



```{r include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      highlight = TRUE,
                      message = FALSE,
                      fig.align = "center",
                      tidy = FALSE,
                      eval = TRUE,
                      fig.width = 7,
                      cache = TRUE,
                      cache.path = "./cache_html/",
                      fig.path = "./cache_html/figs")
# Thanks to this:
# https://www.r-bloggers.com/wrapper-of-knitrinclude_graphics-to-handle-urls-pdf-outputs/
# https://github.com/liao961120/linguisticsdown/blob/master/R/include_graphics2.R
  
include_graphics2 <- function(path, alt_path = NULL, handler = function(path) knitr::asis_output(paste('View', tools::file_ext(path), 'at', path)), ...) {
  if (knitr::is_latex_output()) {
    return(include_graphics_latex(path, alt_path, handler, ...))
  } else {
    return(knitr::include_graphics(path, ...))
  }
}
include_graphics_latex <- function(path, alt_path = NULL, handler = function(path) knitr::asis_output(paste('View', tools::file_ext(path), 'at', path)), ...) {
  # URL
  if (grepl('^https?://', path)) {
     ifelse(use_alt_path(path, alt_path),
            path <- alt_path,
            return(handler(path)))
    ## Download Figure
    dir_path <- paste0('downloadFigs4latex_',
                       tools::file_path_sans_ext(knitr::current_input()))
    if (!dir.exists(dir_path)) dir.create(dir_path)
    file_path <- paste0(dir_path, '/',
                        knitr::opts_current$get()$label, '.',
                        tools::file_ext(path))
    download.file(path, destfile = file_path)
    path <- file_path
  }
  # Local files
  else {
     ifelse(use_alt_path(path, alt_path),
            path <- alt_path,
            return(handler(path)))
  }
  # Insert Figure
  return(knitr::include_graphics(path, ...))
}
use_alt_path <- function(path, alt_path) {
  # Invalid img ext & no alt provided: Don't include in File
  if (inval_latex_img(path) && is.null(alt_path)) return(FALSE)
  # Invalid img ext with alt provided: insert alt-figure
  if (inval_latex_img(path) && !is.null(alt_path)) {
      stopifnot(!inval_latex_img(alt_path))
      return(TRUE)
    }
}
inval_latex_img <- function(path) {
  invalid_ext <- c('svg', 'SVG', 'GIF', 'gif')
  return(tools::file_ext(path) %in% invalid_ext)
}
```

```{r, eval=TRUE, echo = FALSE, cache = FALSE}
# library("rticles")
# library("rmarkdown")
# rmarkdown::draft(file = "2019_downscaleR_GMD.Rmd",
#                  template = "copernicus_article",
#                  package = "rticles", edit = FALSE)
# rmarkdown::render(input = "2019_downscaleR_GMD/2019_downscaleR_GMD.Rmd")
```

# Used packages

To ensure the reproducibility of the paper results as accurately as possible, it is recommended to install the package versions used to compile this notebook. The appropriate package versions are indicated here through their version tags using the `devtools` package function `install_github` (Wickham _et al._ 2020), or alternatively, their commit has:

```{r, eval=FALSE}
devtools::install_github(c("SantanderMetGroup/loadeR.java@v1.1.1",
                           "SantanderMetGroup/climate4R.UDG@0.1.1",
                           "SantanderMetGroup/loadeR@1.6.1",
                           "SantanderMetGroup/transformeR@7005f67",
                           "SantanderMetGroup/visualizeR@v1.6.0"))
```

Alternatively, and updated image of the packages can be installed using the [conda recipe for climate4R](https://github.com/SantanderMetGroup/climate4R/tree/master/conda).  


## Cloud computing with the climate4R Hub

Furthermore, there is a [docker](https://github.com/SantanderMetGroup/climate4R/tree/master/docker) `climate4R` installation available. The docker file also includes the [jupyter](https://jupyter.readthedocs.io/en/latest) framework enabling a direct usage of `climate4R` via the **climate4R Hub**, a cloud-based computing facility to run `climate4R` notebooks on the cloud using the [IFCA/CSIC Cloud Services](https://ifca.unican.es/en-us/research/advanced-computing-and-e-science)).

The `climate4R` packages used in this experiment are next loaded:

```{r,eval=TRUE,message=FALSE}
require(loadeR)
require(transformeR)
require(visualizeR) 
```

Additional packages will be used for convenience. For instance, the package `magrittr` (Bache and Wickham 2014) allows to conveniently concatenate functions via the pipe operator `%>%`. In addition, the `philentropy` package is used for calculating KL Divergences (Drost, 2018). `sp` (Pebesma and Bivand 2005, Bivand _et al._ 2013) and `lattice` (Sarkar 2008) are used for plotting.  

```{r,eval=TRUE,message=FALSE}
require(magrittr)
require(philentropy)
require(sp)
require(lattice)
```


# Analysis of Lamb Weather Types from NCEP reanalysis {#exp1}

We next load the required datasets. First of all, the NCEP-NCAR reanalysis1 (NNRP) will be loaded, using to this aim the Santander MetGroup Climate Data Service (User Data Gateway, UDG), and the package loadeR for remote access. 

Note that prior to remotely accessing the UDG, login is required. To obtain credentials, pease visit the Thredds Administration Panel ([TAP](http://www.meteo.unican.es/udg-tap/home)). These components are further described in Cofiño _et al._ 2018 and Iturbide _et al._ 2019.


```{r, echo = FALSE}
source("~/workspace/jb")
```

```{r}
loginUDG(username, password)
```

The `lonLim` and `latLim` vectors are used in the following to consider the Euro-CORDEX domain as bounding box for data load. We consider the period 1981-2010, following the WMO guidelines on the calculation of climate normals (WMO, 2017). 

```{r, eval = TRUE}
wmo.years <- 1981:2010
lonLim = c(-45, 66)
latLim = c(22, 73)
```

## Loading reanalysis data {#reanalysis}

```{r,eval=FALSE,message=FALSE}
var <- "slp"
dataset <- "http://meteo.unican.es/tds5/dodsC/ncepReanalysis1/ncepReanalysis1_4xDaily.ncml"
ncep <- loadGridData(dataset = dataset,
                     var = var,
                     lonLim = lonLim,
                     latLim = latLim,
                     season = c(12, 1:11),
                     years = wmo.years,
                     time = "DD",
                     aggr.d = "mean")
```

```{r,echo=FALSE}
# save(ncep, file = "ncep.Rdata")
load("ncep.Rdata")
```

The function `clusterGrid` of package `transformeR` is the workhorse for the application of [clustering methods](https://github.com/SantanderMetGroup/transformeR/wiki/Clustering) to climate datasets. Here, we indicate the Lamb Weater Typing through argument `type = "lamb"`. The default options are fine to compute the LWTs as presented in this paper.


```{r,eval=TRUE,message=FALSE}
clusters.ncep <- clusterGrid(ncep, type = "lamb")

## Figure 1: SpatialPlot of annual climatologies from 8-LWTs-subset of NCEP:

wt.names <-  c("A", "ANE", "AE", "ASE", "AS", "ASW", "AW", "ANW", "AN",
               "NE",  "E", "SE",  "S",  "SW",  "W",  "NW",  "N", 
               "C", "CNE", "CE", "CSE", "CS", "CSW", "CW", "CNW", "CN")

#grid of points from Lamb "cross":
centerlon = -5
centerlat = 55
lon.array <- rep(centerlon, times = 16) + c(-5, 5, -15, -5, 5, 15, -15, -5, 5, 15,
                                            -15, -5, 5, 15, -5, 5)
lat.array <- rep(centerlat, times = 16) + c(10, 10, 5, 5, 5, 5, 0, 0, 0, 0,
                                            -5, -5, -5, -5, -10, -10)
coords <- cbind(lon.array, lat.array) %>% sp::SpatialPoints()
l.points <- list("sp.points", coords, col = 1)

subsetLWT <- c(1, 18, 15, 14, 16, 13, 7, 17)
names.subset <- c("A", "C", "W", "SW", "NW", "S", "AW", "N")
freqLWT <- getWT(clusters.ncep) %>% names() %>% table() %>% prop.table()

#To set freqLWTs[i] = 0 when LWTs 'i' does not occur: 
freqLWT <- freqLWT[match(wt.names, names(freqLWT))]
freqLWT[which(is.na(freqLWT))] <- 0
names(freqLWT) <- wt.names
freqLWT <- round(freqLWT * 100, 2)

names.attr <- paste0(names.subset, ": ", freqLWT[subsetLWT], "%")

LWTs.list <- lapply(1:8, function(x) {
  suppressMessages(climatology(subsetGrid(clusters.ncep, cluster = subsetLWT[x])))
})

LWTs.mg <- makeMultiGrid(LWTs.list, skip.temporal.check = TRUE)
dev.new()
breaks <- seq(99300, 103300, 200)
colorkey.labels <- c(99300, 99800, 100300, 100800, 101300, 101800, 102300, 102800, 133000)
```

The function `spatialPlot` from package `visualizeR` (Frías _et al._ 2018) is a wrapper for the `spplot` method in package `sp`, thus accepting the many possible arguments of the lattice framework for fine tuning of the plot. Here, we reproduce the paper Fig. 1 as accurately as possible:

```{r,fig.width=12,fig.height=14,fig.cap='Composite maps of Lamb Weather Types (LWTs) derived from MSLP (Pa) from the NNRP reanalysis for the period 1981-2010. A subset of the 8 (out of 26) most frequent LWTs annually is displayed. Sub-panels are labelled with their LWT abbreviation (frequency in % in parenthesis) and sorted in decreasing frequency order from top to bottom and from left to right. Colorbar is centered on average sea-level atmospheric pressure (reds are highs andblues are lows). Lamb’s cross coordinates are also indicated over the British Isles domain.'}
visualizeR::spatialPlot(LWTs.mg,
                        sp.layout = list(l.points),
                        backdrop.theme = "coastline",
                        rev.colors = TRUE,
                        main = "Lamb WTs from ERA-Interim (1981-2010)",
                        useRaster = TRUE,
                        set.min = min(breaks),
                        set.max = max(breaks),
                        at = breaks, 
                        colorkey = list(space = 'bottom',
                                        labels = list(at = seq(99300, 103300, 500), 
                                                      labels = colorkey.labels)),
                        layout = c(2,4),
                        as.table = TRUE,
                        names.attr = names.attr,
                        contour = TRUE,
                        lty = 3)
```

Next, the LWT frequencies as captured by the NNRN reanalysis are also displayed as a barplot, similar to paper Fig. 2:

```{r, fig.width=12, fig.cap='Comparison of the seasonal relative frequencies of Lamb Weather Types (LWTs) obtained from the NNRP reanalysis. The LWTs are  sorted  in decreasing order of their annual frequencies.'}
# Definition of seasons:

seasons <- list(
  DJF = c(12,1,2),
  MAM = c(3,4,5),
  JJA = c(6,7,8),
  SON = c(9,10,11)
)

ncep.freqs.LWTs <- lapply(1:length(seasons), function(x) {
  grid <- subsetGrid(clusters.ncep, season = seasons[[x]])
  freqLWT <- getWT(grid) %>% names() %>% table() %>% prop.table()
  # Missing LWTs are set to zero frequency: 
  freqLWT <- freqLWT[match(wt.names, names(freqLWT))]
  freqLWT[which(is.na(freqLWT))] <- 0
  names(freqLWT) <- wt.names
  freqLWT <- round(freqLWT * 100, 2)
})

seasonal.freqs <- lapply(1:length(seasons), function(x) {
  sort.int(ncep.freqs.LWTs[[x]], decreasing = TRUE)
})
  
seasonal.freqs.mat <- matrix(as.numeric(unlist(seasonal.freqs)),
                             ncol = length(wt.names),
                             byrow = TRUE,
                             dimnames = list(c("DJF", "MAM", "JJA", "SON"), wt.names))

layout(matrix(c(rep(1, 6), 2), ncol = 1)) 
par(mai = rep(0.6, 4))
bar.colors <- c("#612C69", "#459ED5", "#FCCF61", rgb(0.3, 0.9, 0.4, 0.6))
bp <- barplot(seasonal.freqs.mat, 
              beside = TRUE, col = bar.colors, border = NA, las = 1,
              ylab = "freq. [%]", ylim = c(0, 21),
              main = "NCEP-NCAR Reanalysis (NNRP) LWT frequencies")

par(mai = c(0, 0, 0, 0))
plot.new()
legend(legend = c("DJF", "MAM","JJA", "SON") , 
       fill = c("#612C69", "#459ED5", "#FCCF61", rgb(0.3, 0.9, 0.4, 0.6)), 
       "center", horiz = TRUE, border = "transparent", bty = "n")
```

# Evaluation of EC-EARTH vs. NCEP

## Lamb WTs of EC-Earth: 

### Loading GCM data

```{r,echo=FALSE}
load("ec-earth.Rdata")
```

```{r,eval=FALSE,message=FALSE}
historical = "http://meteo.unican.es/tds5/dodsC/cmip5/EC-EARTH/EC-EARTH/historical/day/ec-earth_ec-earth_historical_r12i1p1.ncml"
rcp8.5 = "http://meteo.unican.es/tds5/dodsC/cmip5/EC-EARTH/EC-EARTH/rcp85/day/ec-earth_ec-earth_rcp85_r12i1p1.ncml"
var <- "psl"

grid1 <- loadGridData(dataset = historical,
                      var = var,
                      lonLim = lonLim,
                      latLim = latLim,
                      years = 1980:2005,
                      time = "DD", 
                      aggr.d = "mean")

# merge of 5 years from rcp8.5 as specified in section 2.1 of the paper:

grid2 <- loadGridData(dataset = rcp8.5,
                      var = var,
                      lonLim = lonLim,
                      latLim = latLim,
                      years = 2006:2010,
                      time = "DD", 
                      aggr.d = "mean")
```

Both datasets (historical and RCP8.5) are next joined with `bindGrid` along their time dimension:

```{r}
ec.earth <- bindGrid(grid1, grid2, dimension = "time")
ec.earth <- subsetGrid(ec.earth, season = c(12,1:11))
```

```{r,echo=FALSE,eval=FALSE}
# save(grid1, grid2, file = "ec-earth.Rdata")
```


The LWTs from the EC-EARTH model are next computed:

```{r,eval=TRUE,message=FALSE}
wts.ec.earth <- clusterGrid(ec.earth, type = "lamb")
```


## Relative Biases between EC-Earth and NCEP seasonal LWTs

```{r,eval=TRUE,message=FALSE,fig.cap='Seasonal relative biases of the eight main LWT frequencies as simulated by EC-EARTH, compared against the NNRP reanalysis.'}
ec.earth.freqs.LWTs <- lapply(1:length(seasons), function(x) {
  grid <- subsetGrid(wts.ec.earth, season = seasons[[x]])
  freqLWT <- getWT(grid) %>% names() %>% table() %>% prop.table()
  # Non existing LWTs are assigned a zero probability  
  freqLWT <- freqLWT[match(wt.names, names(freqLWT))]
  freqLWT[which(is.na(freqLWT))] <- 0
  names(freqLWT) <- wt.names
  freqLWT <- round(freqLWT*100, 2)
})
names(ec.earth.freqs.LWTs) <- c("DJF", "MAM", "JJA", "SON")

rel.bias <- lapply(1:length(seasons), function(x) {
    diff.freqs <- ec.earth.freqs.LWTs[[x]][subsetLWT] - ncep.freqs.LWTs[[x]][subsetLWT]
    diff.freqs/ncep.freqs.LWTs[[x]][subsetLWT]
})
names(rel.bias) <- c("DJF", "MAM", "JJA", "SON")

rel.bias.mat <- matrix(as.numeric(unlist(rel.bias)),
                       ncol = length(subsetLWT), byrow = TRUE,
                       dimnames = list(c("DJF", "MAM", "JJA", "SON"),
                                       names.subset))

RColorBrewer::brewer.pal(n = 9, "RdBu") %>% rev()
pcolors <- RColorBrewer::brewer.pal(n = 9, "RdBu") %>% colorRampPalette()
levelplot(rel.bias.mat,
          ylab = "LWTs", xlab = "Season",
          main = "DJF - LWT Relative Bias",
          col.regions = rev(pcolors(201)),
          at = seq(-0.5, 0.5, 0.05))
```

## Transition probabilities

### TPM of the NNRP reanalysis

The following helper functions from this public repository are used to calculate the tansition probability matrices (TPM) and the associated score (TPMS). Note that the plot depends on the `image.plot` function from package `fields` (Nychka _et al._ 2017), that is also loaded.

```{r, message=FALSE}
require(fields)
source("https://raw.githubusercontent.com/juanferngran/TFM/master/R/tprobPlot.R")
source("https://raw.githubusercontent.com/juanferngran/TFM/master/R/transitionProb.R")
source("https://raw.githubusercontent.com/juanferngran/TFM/master/R/transitionProb.pvalue.R")
source("https://raw.githubusercontent.com/juanferngran/TFM/master/R/transitionProbMatrixScore.R")
```

```{r,fig.asp=1,fig.width=9,fig.cap='Transition probability matrix of the Lamb Weather Type classification, as produced by the NNRP reanalysis for the period 1981-2010.'}
tprobPlot(tprob.matrix = transitionProb(clusters.ncep), title = "Reference TPM (NNRP reanalysis)")
```


### TPM of the EC-EARTH GCM 

```{r,fig.asp=1,fig.width=9,fig.cap='Transition Probability Matrix of CMIP5 EC-EARTH model, considering the Lamb Weather Type classification for the period 1981-2010, using as reference the NNRP reanalysis.'}
tprobPlot(tprob.matrix = transitionProb(wts.ec.earth),
          pval.matrix = transitionProb.test(obs.grid = clusters.ncep,
                                            gcm.grid = wts.ec.earth),
          tprob.ref = transitionProb(clusters.ncep),
          title = "CMIP5 EC-EARTH Transition Probability Matrix")
```

### TPM Score calculation

The TPMS provides an overall score of the similarity between the NNRP and EC-EARTH transition probability measures:

```{r}
TPMS(obs.wt.grid = clusters.ncep, gcm.wt.grid = wts.ec.earth, include.nonexisting = TRUE)
```


## Kullback-Leibler Divergence between EC-Earth and NCEP annual/seasonal LWTs

Next, the KL divergence is computed, used in this study as a measure of departure between the reanalysis and the GCM representation of the Lamb Weather typing classification. To this aim, we use the from package `phylentropy` (Drost, 2018).

```{r,eval=TRUE,message=FALSE}
seasons <- list(
  DJF = c(12,1,2),
  MAM = c(3,4,5),
  JJA = c(6,7,8),
  SON = c(9,10,11),
  year = c(12,1:11)
)

ncep.abs.freqs <- lapply(1:length(seasons), function(x) {
  grid <- subsetGrid(clusters.ncep, season = seasons[[x]])
  freqLWT <- getWT(grid) %>% names() %>% table()
  #To set freqLWTs[i] = 0 when LWTs 'i' does not occur: 
  freqLWT <- freqLWT[match(wt.names, names(freqLWT))]
  freqLWT[which(is.na(freqLWT))] <- 0
  names(freqLWT) <- wt.names
  return(freqLWT)
})
names(ncep.abs.freqs) <- c("DJF", "MAM", "JJA", "SON", "Annual")

ec.earth.abs.freqs <- lapply(1:length(seasons), function(x) {
  grid <- subsetGrid(wts.ec.earth, season = seasons[[x]])
  freqLWT <- attr(x = grid, which = "wt.index") %>% table()
  # Non existing LWTs are assigned a zero probability 
  freqLWT <- freqLWT[match(1:26, names(freqLWT))]
  freqLWT[which(is.na(freqLWT))] <- 0
  names(freqLWT) <- wt.names
  return(freqLWT)
})
names(ec.earth.abs.freqs) <- c("DJF", "MAM", "JJA", "SON", "Annual")

divergence <- sapply(1:length(seasons), function(x) {
  vector1 <- ec.earth.abs.freqs[[x]]
  vector2 <- ncep.abs.freqs[[x]]
  x <- rbind(vector1, vector2)
  philentropy::KL(x, unit = "log",  est.prob = "empirical")
})
names(divergence) <- c("DJF", "MAM", "JJA", "SON", "Annual")
print(divergence)
```

# References
 * Cofiño, A.S., Bedia, J., Iturbide, M., Vega, M., Herrera, S., Fernández, J., Frías, M.D., Manzanas, R., Gutiérrez, J.M., 2018. The ECOMS User Data Gateway: Towards seasonal forecast data provision and research reproducibility in the era of Climate Services. Climate Services 9, 33–43. https://doi.org/10.1016/j.cliser.2017.07.001
 * Drost HG., 2018. Philentropy: Information Theory and Distance Quantification with R. Journal of Open Source Software.
  doi:10.21105/joss.00765
 * Frías, 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
 * Iturbide, M., Bedia, J., Herrera, S., Baño-Medina, J., Fernández, J., Frías, M.D., Manzanas, R., San-Martín, D., Cimadevilla, E., Cofiño, A.S., Gutiérrez, J.M., 2019. The R-based climate4R open framework for reproducible climate data access and post-processing. Environmental Modelling & Software 111, 42–54. https://doi.org/10.1016/j.envsoft.2018.09.009
 * Milton Bache, S. and Wickham, H., 2014. magrittr: A Forward-Pipe Operator for R. R package version 1.5.
  https://CRAN.R-project.org/package=magrittr
 * Nychka, D., Furrer, R., Paige, J. and Sain, S., 2017. “fields: Tools for spatial data.” doi: 10.5065/D6W957CT
(URL: https://doi.org/10.5065/D6W957CT), R package version 10.3, <URL: https://github.com/NCAR/Fields>.
 * Pebesma, E.J., R.S. Bivand, 2005. Classes and methods for spatial data in R. R News 5 (2),
  https://cran.r-project.org/doc/Rnews/.
 * Roger S. Bivand, Edzer Pebesma, Virgilio Gomez-Rubio, 2013. Applied spatial data analysis with R, Second edition.
  Springer, NY. https://asdar-book.org/
 * Sarkar, Deepayan, 2008. Lattice: Multivariate Data Visualization with R. Springer, New York. ISBN 978-0-387-75968-5
 * Wickham, H., Hester, J. and Chang, W., 2020. devtools: Tools to Make Developing R Packages Easier. R package
  version 2.3.0. https://CRAN.R-project.org/package=devtools

# Session info

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