diff --git a/vignettes/Figures/MostLikelyTercile_fig3.png b/vignettes/Figures/MostLikelyTercile_fig3.png new file mode 100644 index 0000000000000000000000000000000000000000..edf922e1a3071e653c35f862e7fbfa0173ce8dc8 Binary files /dev/null and b/vignettes/Figures/MostLikelyTercile_fig3.png differ diff --git a/vignettes/Figures/VizForecastPDF_ex1.png b/vignettes/Figures/VizForecastPDF_ex1.png new file mode 100644 index 0000000000000000000000000000000000000000..0fddba29a26775647f1ee454f8f2c920a07c3dab Binary files /dev/null and b/vignettes/Figures/VizForecastPDF_ex1.png differ diff --git a/vignettes/Figures/VizForecastPDF_ex2.png b/vignettes/Figures/VizForecastPDF_ex2.png new file mode 100644 index 0000000000000000000000000000000000000000..8288d3a1c9bb64669b136fcfb07820e9f3dc58ce Binary files /dev/null and b/vignettes/Figures/VizForecastPDF_ex2.png differ diff --git a/vignettes/Figures/VizForecastPDF_ex3.png b/vignettes/Figures/VizForecastPDF_ex3.png new file mode 100644 index 0000000000000000000000000000000000000000..3525ad0ccef39302421dbcd5e8d47b3e0a4c549e Binary files /dev/null and b/vignettes/Figures/VizForecastPDF_ex3.png differ diff --git a/vignettes/Figures/VizForecastPDF_ex4.png b/vignettes/Figures/VizForecastPDF_ex4.png new file mode 100644 index 0000000000000000000000000000000000000000..95d6d3bcc051abdc1c055a6e221953f81f648e71 Binary files /dev/null and b/vignettes/Figures/VizForecastPDF_ex4.png differ diff --git a/vignettes/Figures/most_likely_tercile_fig1.png b/vignettes/Figures/most_likely_tercile_fig1.png new file mode 100644 index 0000000000000000000000000000000000000000..0150a7f2cf1fa75388598d3e83031450252309c0 Binary files /dev/null and b/vignettes/Figures/most_likely_tercile_fig1.png differ diff --git a/vignettes/Figures/most_likely_tercile_fig2.png b/vignettes/Figures/most_likely_tercile_fig2.png new file mode 100644 index 0000000000000000000000000000000000000000..5ba969e4dba881f34a5364707327524a84d56989 Binary files /dev/null and b/vignettes/Figures/most_likely_tercile_fig2.png differ diff --git a/vignettes/Figures/multi_model_skill_cor_tas_1993-2012.png b/vignettes/Figures/multi_model_skill_cor_tas_1993-2012.png new file mode 100644 index 0000000000000000000000000000000000000000..cad7506c9a9a442180cc891fd8dcae6001db32a7 Binary files /dev/null and b/vignettes/Figures/multi_model_skill_cor_tas_1993-2012.png differ diff --git a/vignettes/Figures/multi_model_skill_rms_tas_1993-2012.png b/vignettes/Figures/multi_model_skill_rms_tas_1993-2012.png new file mode 100644 index 0000000000000000000000000000000000000000..764a6fc20dc86668addec44db2e556c80c538b94 Binary files /dev/null and b/vignettes/Figures/multi_model_skill_rms_tas_1993-2012.png differ diff --git a/vignettes/Figures/multi_model_skill_rmsss_tas_1993-2012.png b/vignettes/Figures/multi_model_skill_rmsss_tas_1993-2012.png new file mode 100644 index 0000000000000000000000000000000000000000..5c4b4a1b6d0574b3e0489105435f6bcaa9f1d911 Binary files /dev/null and b/vignettes/Figures/multi_model_skill_rmsss_tas_1993-2012.png differ diff --git a/vignettes/forecast_pdf_viz.Rmd b/vignettes/forecast_pdf_viz.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..d27d039d6bf9468000c2cb42bbd2b49d09cb13c2 --- /dev/null +++ b/vignettes/forecast_pdf_viz.Rmd @@ -0,0 +1,89 @@ +--- +author: "Francesc Roura-Adserias and Llorenç Lledó" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteEngine{knitr::knitr} + %\VignetteIndexEntry{Viz Forecast PDFs} + %\usepackage[utf8]{inputenc} +--- + +Visualize Forecast PDFs (Probability Distibution Functions) +------------------------------------------ + +Library *esviz*, should be installed and loaded. + +``` +library(esviz) +``` + +### 1. A simple example + +The first step is to put your forecasts in an appropriate format. For this vignette we generate some random values from two normal distributions. The VizForecastPDF by default will plot the ensemble members, the estimated density distributions and the tercile probabilities. + +```{r,fig.show = 'hide',warning=F} +fcst <- data.frame(fcst1 = rnorm(mean = 25, sd = 3, n = 30), + fcst2 = rnorm(mean = 23, sd = 4.5, n = 30)) +VizForecastPDF(fcst, tercile.limits = c(20, 26)) +``` + +![Example 1](./Figures/VizForecastPDF_ex1.png) + +Input data can also be provided as an two-dimensional array, as far as one of the dimensions is named 'member' or adjusted in 'memb_dim' parameter: + +```{r,fig.show = 'hide',warning=F} +fcst <- array(rnorm(mean = 25, sd = 2, n = 90), dim = c(member = 30, 3)) +VizForecastPDF(fcst, tercile.limits = c(23, 27)) +``` + +### 2. Customizing the appearance of your plots +Some parameters allow to customize your plot by changing the title, the forecast labels, the variable name and units, or the colors. + +```{r,fig.show = 'hide',warning=F} +fcst <- data.frame(fcst1 = rnorm(mean = 25, sd = 3, n = 30), + fcst2 = rnorm(mean = 23, sd = 4.5, n = 30)) +VizForecastPDF(fcst, tercile.limits = c(20, 26), var.name = "Temperature (ºC)", + title = "Forecasts valid for 2019-01-01 at Sunny Hills", + fcst.names = c("model a", "model b"), + color.set = "s2s4e") +``` +![Example 2](./Figures/VizForecastPDF_ex2.png) + +### 3. Adding extremes and observed values +Optionally, we can include the probability of extreme values or the actually observed values. The tercile limits, extreme limits and observation values can be specified for each panel separately. + +```{r,fig.show = 'hide',warning=F} +fcst <- data.frame(fcst1 = rnorm(mean = 25, sd = 3, n = 30), + fcst2 = rnorm(mean = 28, sd = 4.5, n = 30), fcst3 = rnorm(mean = 17, sd = 3, n = 30)) +VizForecastPDF(fcst, tercile.limits = rbind(c(20, 26), c(22, 28), c(15, 22)), + var.name = "Temperature (ºC)", title = "Forecasts at Sunny Hills", + fcst.names = c("January", "February", "March"), obs = c(21, 24, 17), + extreme.limits = rbind(c(18, 28), c(20, 30), c(12, 24)), + color.set = "s2s4e") +``` + +![Example 3](./Figures/VizForecastPDF_ex3.png) + +### 4. Saving your plot to a file +VizForecastPDF uses ggplot2, so you can save the output of the function to a variable and operate with it as a ggplot2 object. For instance, you can save it to a file: + +``` +library(ggplot2) +fcst <- array(rnorm(mean = 25, sd = 2, n = 90), dim = c(member = 30, 3)) +plot <- VizForecastPDF(fcst, tercile.limits = c(23, 27)) +ggsave("outfile.pdf", plot, width = 7, height = 5) +``` + +### 5. A reproducible example using lonlat_temp_st +This final example uses the sample lonlat data from CSTools. It is suitable for checking reproducibility of results. + +```{r,fig.show = 'hide',warning=F} +fcst <- data.frame(fcst1 = lonlat_temp_st$exp$data[1,1,,1,1,1,1] - 273.15, + fcst2 = lonlat_temp_st$exp$data[1,1,,1,2,1,1] - 273.15) +VizForecastPDF(fcst, tercile.limits = c(5, 7), extreme.limits = c(4, 8), + var.name = "Temperature (ºC)", + title = "Forecasts initialized on Nov 2000 at sample Mediterranean region", + fcst.names = c("November", "December")) +``` + +![Example 4](./Figures/VizForecastPDF_ex4.png) diff --git a/vignettes/most_likely_tercile.Rmd b/vignettes/most_likely_tercile.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..6d59117f2ce4a1575fc8be48162922b54b7a41bf --- /dev/null +++ b/vignettes/most_likely_tercile.Rmd @@ -0,0 +1,255 @@ +--- +author: "Louis-Philippe Caron and Núria Pérez-Zanón" +date: "`r Sys.Date()`" +revisor: "Eva Rifà" +revision date: "October 2023" +output: rmarkdown::html_vignette +vignette: > + %\VignetteEngine{knitr::knitr} + %\VignetteIndexEntry{Most Likely Terciles} + %\usepackage[utf8]{inputenc} +--- + + +Computing and displaying the most likely tercile of a seasonal forecast +======================== + +In this example, we will use **esviz** and **CSTools** to visualize a probabilistic forecast (most likely tercile) of summer near-surface temperature produced by ECMWF System5. The (re-)forecasts used are initilized on May 1st for the period 1981-2020. The target for the forecast is June-August (JJA) 2020. The forecast data are taken from the Copernicus Climate Data Store. + + +### 1. Preliminary setup + +Library *esviz*, should be installed and loaded. + +``` +library(esviz) +``` + +To run this vignette, the following R packages should be installed and loaded: + + +``` +library(CSTools) +library(s2dv) +library(zeallot) +library(easyVerification) +library(ClimProjDiags) +``` + +### 2. Loading the data + + +We first define a few parameters. We start by defining the region we are interested in. In this example, we focus on the Mediterranean region. + + +```r +lat_min = 25 +lat_max = 52 +lon_min = -10 +lon_max = 40 +``` + +If the boundaries are not specified, the domain will be the entire globe. + +We also define the start dates for the hindcasts/forecast (in this case, May 1st 1981-2020) and create a sequence of dates that will be required by the load function. + + +```r +ini <- 1981 +fin <- 2020 +numyears <- fin - ini +1 +mth = '05' +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") +``` + +We then define the target months for the forecast (i.e. JJA). The months are given relative to the start date (May in this case) considering that monthly simulations are being analyzed. + + +```r +mon1 <- 2 +monf <- 4 +``` + +Finally, we define the forecast system, an observational reference, the variable of interest and the common grid onto which to interpolate. + + +```r +clim_var = 'tas' +``` + +Finally, the data are loaded using `CST_Start`: + +```r +repos_exp <- paste0('/esarchive/exp/ecmwf/system5c3s/monthly_mean/', + '$var$_f6h/$var$_$sdate$.nc') +exp <- CST_Start(dataset = repos_exp, + var = clim_var, + member = startR::indices(1:10), + sdate = dateseq, + ftime = startR::indices(2:4), + lat = startR::values(list(lat_min, lat_max)), + lat_reorder = startR::Sort(decreasing = TRUE), + lon = startR::values(list(lon_min, lon_max)), + lon_reorder = startR::CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + member = c('member', 'ensemble'), + ftime = c('ftime', 'time')), + transform = startR::CDORemapper, + transform_extra_cells = 2, + transform_params = list(grid = 'r256x128', + method = 'bilinear'), + transform_vars = c('lat', 'lon'), + return_vars = list(lat = NULL, + lon = NULL, ftime = 'sdate'), + retrieve = TRUE) + +# Give the correct time values (identical as the netCDF files) +dates_obs <- c(paste0(ini:fin, '-06-30 18:00:00'), + paste0(ini:fin, '-07-31 18:00:00'), + paste0(ini:fin, '-08-31 18:00:00')) +dates_obs <- as.POSIXct(dates_obs, tz = "UTC") +dim(dates_obs) <- c(sdate = 40, ftime = 3) + +date_arr <- array(format(dates_obs, '%Y%m'), dim = c(sdate = 40, ftime = 3)) + +repos_obs <- paste0('/esarchive/recon/ecmwf/erainterim/monthly_mean/', + '$var$/$var$_$date$.nc') + +obs <- CST_Start(dataset = repos_obs, + var = clim_var, + date = date_arr, + split_multiselected_dims = TRUE, + lat = startR::values(list(lat_min, lat_max)), + lat_reorder = startR::Sort(decreasing = TRUE), + lon = startR::values(list(lon_min, lon_max)), + lon_reorder = startR::CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + ftime = c('ftime', 'time')), + transform = startR::CDORemapper, + transform_extra_cells = 2, + transform_params = list(grid = 'r256x128', + method = 'bilinear'), + transform_vars = c('lat', 'lon'), + return_vars = list(lon = NULL, + lat = NULL, + ftime = 'date'), + retrieve = TRUE) +``` + +Loading the data using CST_Start returns two objects, one for the experimental data and another one for the observe data, with the same elements and compatible dimensions of the data element: + + +```r +> dim(exp$data) +dataset var member sdate ftime lat lon + 1 1 10 40 3 19 36 +> dim(obs$data) +dataset var sdate ftime lat lon + 1 1 40 3 19 36 +``` + + +The latitude and longitude are saved for later use: + + +```r +Lat <- exp$coords$lat +Lon <- exp$coords$lon +``` + +### 3. Computing probabilities + +First, anomalies of forecast and observations are computed using cross-validation on individual members: + +``` +c(Ano_Exp, Ano_Obs) %<-% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE) +``` + +The seasonal mean of both forecasts and observations are computed by averaging over the ftime dimension. + +```r +Ano_Exp$data <- MeanDims(Ano_Exp$data, 'ftime') +Ano_Obs$data <- MeanDims(Ano_Obs$data, 'ftime') +``` + +Finally, the probabilities of each tercile are computed by evaluating which tercile is forecasted by each ensemble member for the latest forecast (2020) using the function `ProbBins` in **s2dv** and then averaging the results along the member dimension to obtain the probability of each tercile. + +```r +PB <- ProbBins(Ano_Exp$data, fcyr = numyears, thr = c(1/3, 2/3), compPeriod = "Without fcyr") +prob_map <- MeanDims(PB, c('sdate', 'member', 'dataset', 'var')) +``` + +### 4. Visualization with VizMostLikelyQuantileMap + + +We then visualize the most likely quantile using the **esviz** function `VizMostLikelyQuantileMap`. + + +``` +VizMostLikelyQuantileMap(probs = prob_map, lon = Lon, lat = Lat, + coast_width = 1.5, legend_scale = 0.5, + toptitle = paste0('Most likely tercile - ', clim_var, + ' - ECMWF System5 - JJA 2020'), + width = 10, height = 8) +``` + +![](./Figures/most_likely_tercile_fig1.png) + +The forecast calls for above average temperature over most of the Mediterranean basin and near average temperature for some smaller regions as well. But can this forecast be trusted? + +For this, it is useful evaluate the skill of the system at forecasting near surface temperature over the period for which hindcasts are available. We can then use this information to mask the regions for which the system doesn't have skill. + +In order to do this, we will first calculate the ranked probability skill score (RPSS) and then exclude/mask from the forecasts the regions for which the RPSS is smaller or equal to 0 (no improvement with respect to climatology). + + +### 5. Computing Skill Score + + +First, we evaluate and plot the RPSS. Therefore, we use `RPSS` metric included in CST_MultiMetric from function from **CSTools** package which requires to remove missing values from latest start dates: + + +```r +Ano_Exp <- CST_Subset(Ano_Exp, along = 'sdate', indices = 1:38) +Ano_Obs <- CST_Subset(Ano_Obs, along = 'sdate', indices = 1:38) +Ano_Obs <- CST_InsertDim(Ano_Obs, posdim = 3, lendim = 1, name = "member") + +RPSS <- CST_MultiMetric(Ano_Exp, Ano_Obs, metric = 'rpss', multimodel = FALSE) + +VizEquiMap(RPSS$data[[1]], lat = Lat, lon = Lon, brks = seq(-1, 1, by = 0.1), + filled.continents = FALSE) +``` + +![](./Figures/most_likely_tercile_fig2.png) + + +Areas displayed in red (RPSS > 0) are areas for which the forecast system shows skill above climatology whereas areas in blue (such as a large part of the Iberian Peninsula) are areas for which the model does not. We thus want to mask the areas currently displayed in blue. + +### 6. Simultaneous visualization of probabilities and skill scores + +From the RPSS, we create a mask: regions with RPSS <= 0 will be masked. + + +```r +mask_rpss <- ifelse((RPSS$data$skillscore <= 0) | is.na(RPSS$data$skillscore), 1, 0) +``` + +Finally, we plot the latest forecast, as in the previous step, but add the mask we just created. + + +```r +VizMostLikelyQuantileMap(probs = prob_map, lon = Lon, lat = Lat, coast_width = 1.5, + legend_scale = 0.5, mask = mask_rpss[ , , ,], + toptitle = paste('Most likely tercile -', clim_var, + '- ECMWF System5 - JJA 2020'), + width = 10, height = 8) +``` + +![](./Figures/most_likely_tercile_fig3.png) + +We obtain the same figure as before, but this time, we only display the areas for which the model has skill at forecasting the right tercile. The gray regions represents areas where the system doesn't have sufficient skill over the verification period. + + diff --git a/vignettes/multimodel_skill_assessment.Rmd b/vignettes/multimodel_skill_assessment.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..e0b16e070937dcee78b2f2be11223eb8d6b04cb6 --- /dev/null +++ b/vignettes/multimodel_skill_assessment.Rmd @@ -0,0 +1,265 @@ +--- +author: "Nuria Perez" +date: "`r Sys.Date()`" +revisor: "Eva Rifà" +revision date: "October 2023" +output: rmarkdown::html_vignette +vignette: > + %\VignetteEngine{knitr::knitr} + %\VignetteIndexEntry{Multi-model Skill Assessment} + %\usepackage[utf8]{inputenc} +--- + +Multi-model Skill Assessment and Visualization +----------------------------------------- + +In this example we will do a multi-model skill assesment and visualization using **esviz** and **CSTools**. + + +**Reference**: Mishra, N., Prodhomme, C., & Guemas, V. (2018). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31. + +Library *esviz*, should be installed and loaded. + +``` +library(esviz) +``` + +The following R packages should be loaded by running: + +``` +library(s2dv) +library(zeallot) +library(CSTools) +``` + + +### 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, Meteo-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 1993-2012. So, the starting and ending dates can be defined by running the following lines: + +``` +ini <- 1993 +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") +``` + +The grid in which all data will be interpolated needs to be specified within the `CST_Start` call (256x128 grid). The observational dataset used in this example is the EraInterim. + +Using the `CST_Start` function, 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 +lonmin = -20 +lonmax = 70 +latmin = 25 +latmax = 75 +repos1 <- "/esarchive/exp/glosea5/glosea5c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc" +repos2 <- "/esarchive/exp/ecmwf/system4_m1/monthly_mean/$var$_f6h/$var$_$sdate$01.nc" +repos3 <- "/esarchive/exp/meteofrance/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$01.nc" + +exp <- CST_Start(dataset = list(list(name = 'glosea5c3s', path = repos1), + list(name = 'ecmwf/system4_m1', path = repos2), + list(name = 'meteofrance/system5_m1', path = repos3)), + var = clim_var, + member = startR::indices(1:4), + sdate = dateseq, + ftime = startR::indices(2:4), + lat = startR::values(list(latmin, latmax)), + lat_reorder = startR::Sort(decreasing = TRUE), + lon = startR::values(list(lonmin, lonmax)), + lon_reorder = startR::CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), lat = c('lat', 'latitude'), + member = c('member', 'ensemble'), ftime = c('ftime', 'time')), + transform = startR::CDORemapper, + transform_extra_cells = 2, + transform_params = list(grid = 'r256x128', method = 'bilinear'), + transform_vars = c('lat', 'lon'), + return_vars = list(lat = 'dataset', lon = 'dataset', ftime = 'sdate'), + retrieve = TRUE) + +dates_exp <- exp$attrs$Dates +repos_obs <- "/esarchive/recon/ecmwf/erainterim/monthly_mean/$var$/$var$_$date$.nc" +obs <- CST_Start(dataset = list(list(name = 'erainterim', path = repos_obs)), + var = clim_var, + date = unique(format(dates_exp, '%Y%m')), + ftime = startR::values(dates_exp), + ftime_across = 'date', + ftime_var = 'ftime', + merge_across_dims = TRUE, + split_multiselected_dims = TRUE, + lat = startR::values(list(latmin, latmax)), + lat_reorder = startR::Sort(decreasing = TRUE), + lon = startR::values(list(lonmin, lonmax)), + lon_reorder = startR::CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + ftime = c('ftime', 'time')), + transform = startR::CDORemapper, + transform_extra_cells = 2, + transform_params = list(grid = 'r256x128', + method = 'bilinear'), + transform_vars = c('lat', 'lon'), + return_vars = list(lon = NULL, + lat = NULL, + ftime = 'date'), + retrieve = TRUE) + +# save(exp, obs, file = "../tas_toydata.RData") + +# Or use the following line to load the file provided in .RData format: +# load(file = "./tas_toydata.RData") +``` + +There should be two new elements loaded in the R working environment: `exp` and `obs`, containing the experimental and the observed data for temperature. It is possible to check that they are of class `s2dv_cube` by running: + +``` +class(exp) +class(obs) +``` + +The corresponding data is saved in the element `data` of each object, while other relevant information is saved in different elements, such as `lat` and `lon`: + +```r +> dim(exp$data) +dataset var member sdate ftime lat lon + 3 1 9 20 3 35 64 +> dim(obs$data) +dataset var sdate ftime lat lon + 1 1 20 3 35 64 +Lat <- exp$coords$lat +Lon <- exp$coords$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). + +First step is to compute the anomalies over the loaded data applying cross validation technique on individual members by running: + +```r +c(ano_exp, ano_obs) %<-% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE) +``` + +The dimensions are preserved: + +``` +> str(ano_exp$data) + num [1:20, 1:3, 1:9, 1, 1:3, 1:35, 1:64] -1.399 -0.046 -0.133 0.361 -5.696 ... +> str(ano_obs$data) + num [1:20, 1, 1, 1, 1:3, 1:35, 1:64] 1.556 1.397 -0.346 -5.99 -0.273 ... +``` + +The ACC is obtained by running the `CST_MultiMetric` function defining the parameter 'metric' as correlation. The function also includes the option of computing the Multi-Model Mean ensemble (MMM). + + +```r +ano_obs <- CST_InsertDim(ano_obs, posdim = 3, lendim = 1, name = "member") +AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'correlation', + multimodel = TRUE) +``` + +The output of the function `CST_MultiMetric` is a object of class `s2dv_cube`, it contains the result of the metric, in this case correlation, in the `data` element (including the correlation for the MMM in the latest position). +While other relevant data is being stored in the corresponding element of the object: + + +```r +> str(AnomDJF$data) +List of 4 + $ corr : num [1:4, 1, 1, 1:35, 1:64] 0.3061 0.4401 0.0821 0.2086 0.1948 ... + $ p.val : num [1:4, 1, 1, 1:35, 1:64] 0.0947 0.0261 0.3653 0.1887 0.2052 ... + $ conf.lower: num [1:4, 1, 1, 1:35, 1:64] -0.15782 -0.00297 -0.37399 -0.25768 -0.27106 ... + $ conf.upper: num [1:4, 1, 1, 1:35, 1:64] 0.659 0.739 0.506 0.596 0.587 ... +> names(AnomDJF) +[1] "data" "dims" "coords" "attrs" +> AnomDJF$attrs$Datasets +[1] "glosea5" "ecmwf/system4_m1" "meteofrance/system5_m1" "erainterim" +``` + +In the element $data of the `AnomDJF` object is a list of object for the metric and its statistics: correlation, p-value, the lower limit of the 95% confidence interval and the upper limit of the 95% confidence interval and the 95% significance level given by a one-sided T-test. + +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 +VizCombinedMap(AnomDJF$data$corr[,1,1,,], lon = Lon, lat = Lat, map_select_fun = max, + display_range = c(0, 1), map_dim = 'nexp', + legend_scale = 0.5, brks = 11, + cols = list(c('white', 'black'), + c('white', 'darkblue'), + c('white', 'darkred'), + c('white', 'darkorange')), + bar_titles = c("MMM", AnomDJF$attrs$Datasets), + 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 1993-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](./Figures/multi_model_skill_cor_tas_1993-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 <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rms', + multimodel = TRUE) +``` + +The following lines are necessary to obtain the plot which visualizes the best model given this metric for each grid point. + +```r +VizCombinedMap(AnomDJF$data$rms[,1,1,,], lon = Lon, lat = Lat, map_select_fun = min, + display_range = c(0, ceiling(max(abs(AnomDJF$data$rms)))), map_dim = 'nexp', + legend_scale = 0.5, brks = 11, + cols = list(c('black', 'white'), + c('darkblue', 'white'), + c('darkred', 'white'), + c('darkorange', 'white')), + bar_titles = c("MMM", AnomDJF$attrs$Datasets), + width = 14, height = 8) +``` + +![Max Skills RMS](./Figures/multi_model_skill_rms_tas_1993-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 `VizCombinedMap` function (see *esviz R package documentation*) has been defined in order to select the best model. + + +```r +AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rmsss', + multimodel = TRUE) + +VizCombinedMap(AnomDJF$data$rmsss[,1,1,,], lon = Lon, lat = Lat, + map_select_fun = function(x) {x[which.min(abs(x - 1))]}, + display_range = c(0, + ceiling(max(abs(AnomDJF$data$rmsss)))), map_dim = 'nexp', + legend_scale = 0.5, brks = 11, + cols = list(c('white', 'black'), + c('white', 'darkblue'), + c('white', 'darkred'), + c('white', 'darkorange')), + bar_titles = c("MMM", AnomDJF$attrs$Datasets), + width = 14, height = 8) +``` + +![Max Skills RMSSS](./Figures/multi_model_skill_rmsss_tas_1993-2012.png)