---
author: "Louis-Philippe Caron and Núria Pérez-Zanón"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteEngine{knitr::knitr}
  %\VignetteIndexEntry{Most Likely Tercile}
  %\usepackage[utf8]{inputenc}
---


Computing and displaying the most likely tercile of a seasonal forecast
========================

In this example, we will use 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. 

### 1- Required packages 

To run this vignette, the following R packages should be installed and loaded:

```r
library(CSTools)
library(s2dverification)
library(multiApply)
library(zeallot)
```

### 2 - We then define a few parameters.

We first define the region for the forecast. 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).
```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
forecastsys <- 'system5c3s'
obs <- 'erainterim'
grid <- "256x128"
clim_var = 'tas'
```

### 3 - We load the data.
```r
c(exp,obs) %<-% CST_Load(var = clim_var, exp = forecastsys,
                        obs = obs, sdates = dateseq, leadtimemin = mon1, leadtimemax = monf,
                        lonmin = lon_min, lonmax = lon_max, latmin = lat_min, latmax = lat_max, 
                        storefreq = "monthly", sampleperiod = 1,  nmember = 10, 
                        output = "lonlat", method = "bilinear", 
                        grid = paste("r", grid, sep = ""))
```

### 4 - We extract the latitude and longitude of the data.
```r
Lat <- exp$lat
Lon <- exp$lon
```
#
### 5 - We extract the anomalies from the forecasts and the observations.
```r
c(ano_exp, ano_obs) %<-% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE)
```

### 6 - We compute the seasonal mean of forecasts and observations.
```r
ano_exp$data <- Mean1Dim(ano_exp$data,4)
ano_obs$data <- Mean1Dim(ano_obs$data,4)
```


### 7 - We calculate the probability of each tercile for the latest forecast (2020) using the function ProbBins in s2dv.
```r
PB <- ProbBins(ano_exp$data, fcyr = numyears, thr = c(1/3, 2/3), quantile = TRUE, 
               posdates = 3, posdim = 2, compPeriod = "Without fcyr")
```


### 8 - We plot the most likely quantile.
```
prob_map <- Mean1Dim(Mean1Dim(Mean1Dim(PB,3),2),2)
CSTools::PlotMostLikelyQuantileMap(probs = prob_map, lon = Lon, lat = Lat,coast_width=1.5, legend_scale = 0.8, 
                                   toptitle = paste0('Most likely tercile - ',variable,' - ECMWF System5 - JJA 2020'),
                                   width = 12, height = 7, fileout = paste0(clim_var, '_most_likely_tercile.png'))
```
<img src="./Figures/MostLikelyTercile_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).

### 9 - First, we evaluate and plot the RPSS.
```r
v2 <- aperm(drop(ano_obs$data),c(3,2,1))
v2 <- v2[,,1:(numyears-2)]
v1 <- aperm(drop(ano_exp$data),c(4,3,2,1))
v1 <- v1[,,1:(numyears-2),]

RPSS <- veriApply('FairRpss',fcst=v1,obs=v2,ensdim=4,tdim=3,prob =c(1/3,2/3))

PlotEquiMap(RPSS[[1]], lat=Lat,lon=Lon, brks=seq(-1,1,by=0.1),filled.continents=F, 
            fileout = paste0(figures_path, clim_var, '_RPSS.png'))
```
<img src="./Figures/MostLikelyTercile_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 will mask the areas in blue.

### 10 - From the RPSS, we create a mask. Regions with RPSS<=0 will be masked. 
```r
mask_rpss <- RPSS[[1]]
mask_rpss[RPSS[[1]] <= 0] <- 1
mask_rpss[is.na(RPSS[[1]])] <- 1
mask_rpss[RPSS[[1]] > 0] <- 0
```

### 11 - Finally, we plot the latest forecast, as in the previous step, but add the mask we just created.
```r
CSTools::PlotMostLikelyQuantileMap(probs = prob_map, lon = Lon, lat = Lat,coast_width=1.5, legend_scale = 0.8, 
                                   mask = t(mask_rpss),
                                   toptitle = paste0('Most likely tercile - ',variable,' - ECMWF System5 - JJA 2020'),
                                   width = 12, height = 7, 
                                   fileout = paste0(clim_var, '_most_likely_tercile_mask.png'))
```
<img src="./Figures/MostLikelyTercile_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 area where the system doesn't have sufficient skill over the verification period.