Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
Multi-model Skill Assessment
-----------------------------------------
**reference**: Mishra, N., Prodhomme, C., & Guemas, V. (2018). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29–31. <https://link.springer.com/article/10.1007%2Fs00382-018-4404-z>
A version of s2dverification under development should be loaded by running:
```r
devtools::install_git('https://earth.bsc.es/gitlab/es/s2dverification',
branch = 'mostlikely')
library(s2dverification)
```
The functions belonging to the future *CSTools* package would load to your R session, by running the following line after cloning the GitLab repository:
```r
source("./R/AnoMultiMetric.R")
```
*Notice that it's suppose that your working directory is the GitLab project *medscope* you've already cloned.
*Also notice that to clone the repository, you should run in your terminal:
`git clone https://earth.bsc.es/gitlab/external/cstools.git`
to create the folder containing all functions
### 1.- Load data
In this case, the seasonal temperature forecasted, initialized in November, will be used to assess the EUROSIP multi-model seasonal forecasting system consists of a number of independent coupled seasonal forecasting systems integrated into a common framework. From September 2012, the systems include those from ECMWF, the Met Office, Météo-France and NCEP.
The parameters defined are the initializating month and the variable:
```{r cars}
mth = '11'
clim_var = 'tas'
```
The simulations available for these models are covering 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, shows how this function can be used. However, the data is loaded from a previous saved `.RData` file:
Ask nuria.perez at bsc.es to achieve 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')
#StartData <- Load(clim_var, exp = list(glosea5, list(name = 'ecmwf/system4_m1'),
# list(name = 'meteofrance/system5_m1')),
# obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4,
# lonmin = -15, lonmax = 45, latmin = 25, latmax = 50,
# storefreq = "monthly", sampleperiod = 1, nmember = 9,
# output = "lonlat", method = "bilinear",
# grid = paste("r", grid, sep = ""))
#save(StartData, StartData, file = "tas_toydata.RData")
#load(file = "./data/tas_toydata.RData")
```
In this example the first two elements of `StartData` corresponds to arrays containing the simulated and the observed data:
```r
> dim(StartData$mod)
dataset member sdate ftime lat lon
3 9 21 3 18 43
> dim(StartData$obs)
dataset member sdate ftime lat lon
1 1 21 3 18 43
```
Latitudes and longitudes of the common grid can be saved:
```r
Lat <- StartData$lat
Lon <- StartData$lon
```
### 2.- Computing and plotting Anomaly Correlation Coefficient
The Anomaly Correlation Coefficient (ACC) is the most widely used skill metric for Seasonal Climate Forecast quality (Mishra et al., 2018). The ACC is obtained by running the `AnoMultiMetric` function defining the parameter 'metric' as correlation. Furthermore, the function includes the option of computing the Multi-Model Mean ensemble (MMM).
```r
AnomDJF <- AnoMultiMetric(data = StartData, metric = 'correlation',
multimodel = TRUE)
```
The first element of the list returned by the function `AnoMultiMetric` contains the anomaly data for the simulations and the observations, where the second is the correlations (including the correlation for the MMM in the latest position) and the third the names of the models being compared.
```r
> str(AnomDJF)
List of 3
$ ano :List of 2
..$ ano_exp: num [1:3, 1:9, 1:21, 1:3, 1:18, 1:43] 0.637 -0.295 -0.623 1.419 0.464 ...
..$ ano_obs: num [1, 1, 1:21, 1:3, 1:18, 1:43] -0.488 0.709 2.253 -1.842 -2.125 ...
$ metric: num [1:4, 1, 1:4, 1:18, 1:43] 0.4766 0.4937 0.6845 0.5677 0.0565 ...
$ names : chr [1:4] "exp1" "ecmwf/system4_m1" "meteofrance/system5_m1" "MMM"
```
In the second element of the list `AnomDJF`, `metric`, the third dimension contains the lower limit of the 95% confidence interval, the correlation, the upper limit of the 95% confidence interval and the 95% significance level given by a one-sided T-test.
```r
corre <- AnomDJF$metric[ , , 2, , ]
names(dim(corre)) <- c("maps", "lat", "lon")
```
To obtain a spatial plot with a scale from -1 to 1 value of correlation for the model with the highest correlation for each grid point, the following lines should be run:
```r
PlotCombinedMap(corre, lon = Lon, lat = Lat, map_select_fun = max,
display_range = c(-1, 1), map_dim = 'maps',
legend_scale = 0.5, brks = 11,
cols = list(c('white', 'darkblue'),
c('white', 'darkred'),
c('white', 'darkorange'),
c('white', 'black')),
bar_titles = c('glosea5', AnomDJF$names[-1]),
fileout = "./vignettes/Figures/MultiModelSkill_cor_as_1992-2012.png",
width = 14, height = 8)
```
The next figure is the map of the maximum positive Anomaly Correlation Coefficient (ACC) among the three individual models from EUROSIP and the multimodel ensemble. ACC for each model is calculated between their respective predicted ensemble mean anomalies and the anomalies of the observed temperature obtained from ERAINT for winter (DJF) seasons over the period 1992-2012. Blue, red, yellow and black colors indicate that the maximum correlation is obtained for GloSea5, ECMWF, MF and the Multi-Model Mean respectively (similar to figure 3 in Mishra et al., 2018).
![Max Skills Correlation](../vignettes/Figures/MultiModelSkill_cor_tas_1992-2012.png)
### 3.- Computing and plotting Root Mean Square error (RMS)
The same function can be used to compute the RMS error by defining the parameter `metric` as 'rms'.
```r
AnomDJF <- AnoMultiMetric(data = StartData, metric = 'rms',
multimodel = TRUE)
RMS <- AnomDJF$metric[ , , 2, , ]
```
The following lines are necessary to obtain the plot which visualizes the best model given this metric for each grid point.
```r
names(dim(RMS)) <- c("maps", "lat", "lon")
PlotCombinedMap(RMS, lon = Lon, lat = Lat, map_select_fun = min,
display_range = c(0, ceiling(max(abs(RMS)))), map_dim = 'maps',
legend_scale = 0.5, brks = 11,
cols = list(c('darkblue', 'white'),
c('darkred', 'white'),
c('darkorange', 'white'),
c('black', 'white')),
bar_titles = c('glosea5', AnomDJF$names[-1]),
fileout = "./vignettes/Figures/MultiModelSkill_rms_tas_1992-2012.png",
width = 14, height = 8)
```
![Max Skills RMS](../vignettes/Figures/MultiModelSkill_rms_tas_1992-2012.png)
### 4.- Computing and plotting Root Mean Square error Skill Scores (RMSSS)
By running the following lines a plot for the best model given the RMSSS is obtained.
When parameter `metric` is defined as `rmsss`, the RMSSS are stored in the first position on the third dimension of the `metric` component in the AnoMultiMetric output.
Notice that the perfect RMSSS is 1 and the parameter `map_select_fun` from `PlotCombinedMap` function (see *s2dverification R package*) has been defined in order to select the best model.
```r
AnomDJF <- AnoMultiMetric(data = StartData, metric = 'rmsss',
multimodel = TRUE)
RMSSS <- AnomDJF$metric[ , , 1, , ]
names(dim(RMSSS)) <- c("maps", "lat", "lon")
PlotCombinedMap(RMSSS, lon = Lon, lat = Lat,
map_select_fun = function(x) {x[which.min(abs(x - 1))]},
display_range = c(0,
ceiling(max(abs(RMSSS)))), map_dim = 'maps',
legend_scale = 0.5, brks = 11,
cols = list(c('white', 'darkblue'),
c('white', 'darkred'),
c('white', 'darkorange'),
c('white', 'black')),
bar_titles = c('glosea5', AnomDJF$names[-1]),
fileout = "./vignettes/Figures/MultiModelSkill_rmsss_tas_1992-2012.png",
width = 14, height = 8)
```
![Max Skills RMSSS](../vignettes/Figures/MultiModelSkill_rmsss_tas_1992-2012.png)