Newer
Older
author: "Nuria Perez"
date: "`r Sys.Date()`"
output: html_document
vignette: >
%\VignetteEngine{R.rsp::md}
%\VignetteIndexEntry{Multivariate RMSE}
%\usepackage[utf8]{inputenc}
Multivariate Root Mean Square Error (RMSE)
------------------------------------------
To run this vignette, the next R packages should be installed and loaded:
```r
library(s2dverification)
library(startR)
library(RColorBrewer)
```
The function hosted in the *CSTools* gitlab project should be loaded by running:
devtools::install_git('https://earth.bsc.es/gitlab/external/cstools', branch = 'master')
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
```
### 1.- Load data
In this example, the seasonal temperature and precipitation forecasts, initialized in november, will be used to assess the glosea5 seasonal forecasting system from the Met Office, by computing the multivariate RMSE for both temperature and precipitation.
The parameters defined are the initializing month and the variables:
```{r cars}
mth = '11'
temp = 'tas'
precip = 'prlr'
```
The simulations available for this model cover the period 1992-2012. So, the starting and ending dates can be defined by running the following lines:
```r
ini <- 1992
fin <- 2012
start <- as.Date(paste(ini, mth, "01", sep = ""), "%Y%m%d")
end <- as.Date(paste(fin, mth, "01", sep = ""), "%Y%m%d")
dateseq <- format(seq(start, end, by = "year"), "%Y%m%d")
```
The grid in which all data will be interpolated should be also specified. The observational dataset used in this example is the EraInterim.
```r
grid <- "256x128"
obs <- "erainterim"
```
Using the `Load` function from **s2dverification package**, the data available in our data store can be loaded. The following lines show how this function can be used. Here, the data is loaded from a previous saved `.RData` file:
Ask deborah.verfaillie@bsc.es for the data to run the recipe.
```r
#glosea5 <- list(path = '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/
# $VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc')
# Temp <- Load(temp, exp = list(glosea5), obs = obs, sdates = dateseq, leadtimemin = 2,
# leadtimemax = 4, lonmin = -20, lonmax = 70, latmin = 25, latmax = 75,
# storefreq = "monthly", sampleperiod = 1, nmember = 9,
# output = "lonlat", method = "bilinear",
# grid = paste("r", grid, sep = ""))
# Precip <- Load(precip, exp = list(glosea5), obs = obs, sdates = dateseq,
# leadtimemin = 2, leadtimemax = 4, lonmin = -20, lonmax = 70, latmin = 25,
# latmax = 75, storefreq = "monthly", sampleperiod = 1, nmember = 9,
# output = "lonlat", method = "bilinear", grid = paste("r", grid, sep = ""))
# save(Temp, Precip, file = "tas_prlr_toydata.RData")
load(file = "tas_prlr_toydata.RData")
```
In this example the first two elements of `Temp` and `Precip` correspond to arrays containing the simulated and the observed data:
```r
> dim(Temp$mod)
dataset member sdate ftime lat lon
1 9 21 3 35 64
> dim(Temp$obs)
dataset member sdate ftime lat lon
1 1 21 3 35 64
```
Latitudes and longitudes of the common grid can be saved:
```r
Lat <- Temp$lat
Lon <- Temp$lon
```
A list containing the `Temp` and `Precip` lists should be put together to serve as input of the function to compute multivariate RMSEs.
Furthermore, some weights can be applied to the difference variables based on their relative importance (if no weights are given, a value of 1 is automatically assigned to each variable). For this example, we'll give a weight of 2 to the temperature dataset and a weight of 1 to the precipitation dataset:
```r
StartData <- list(Temp, Precip)
weight <- c(2, 1)
```
### 2.- Computing and plotting multivariate RMSEs
The multivariate RMSE gives an indication of the forecast performance (RMSE) for multiple variables simultaneously. Variables can be weighted based on their relative importance.
It is obtained by running the `CST_MultivarRMSE` function:
mvrmse <- CST_MultivarRMSE(StartData, weight)
The function `CST_MultivarRMSE` returns the multivariate RMSE value for 2 or more variables. The output is a CSTool object containing a list of elements being the first one the mvrmse:
```r
> str(mvrmse)
$ mvrmse: num [1, 1, 1, 1:35, 1:64] 0.764 0.8 0.67 0.662 0.615 ...
```
The following lines plot the multivariate RMSE
```r
PlotEquiMap(mvrmse$mvrmse, lon = Lon, lat = Lat, filled.continents = FALSE,
toptitle = "Multivariate RMSE tas, prlr 1992 - 2012", colNA = "white",
bar_limits = c(0,2.5), cols = brewer.pal(n=5,name='Reds'),
fileout = "./MultivarRMSE_gloseas5_tas_prlr_1992-2012.png")
```
