Newer
Older
Weather regime analysis
========================
Weather regimes are a set of patterns able to characterize the different circulation structures occurring each month/season. This vignette aims to illustrate a step-by-step example of how to perform a weather regime assessment using s2dverification functionalities.
### 1- Load s2dverification
The package version including the weather regimes functions can be installed and loaded with the following commands:
devtools::install_git('https://earth.bsc.es/gitlab/es/s2dverification',
library(s2dverification)
- The documentation of the package can be found here:
https://cran.r-project.org/web/packages/s2dverification/s2dverification.pdf
### 2- Load files with the `Load ()` function included in the package.
The full description of `Load()` and examples of its use can be found [*here*] (https://earth.bsc.es/gitlab/es/s2dverification/blob/develop-vignettes/vignettes/data_retrieval.md).
The data employed in this example are described below.
- Sea level pressure (psl): this has been selected as the circulation variable, however other variables such as geopotential at 500 hPa can be also used.
- Region: Euro-Atlantic domain [85.5ºW-45ºE; 27-81ºN].
- Datasets: seasonal predictions from ECMWF System 4 ([**Molteni et al. 2011**] (https://www.ecmwf.int/sites/default/files/elibrary/2011/11209-new-ecmwf-seasonal-forecast-system-system-4.pdf)) and ERA-Interim reanalysis ([**Dee et al. 2011**] (http://onlinelibrary.wiley.com/doi/10.1002/qj.828/pdf)) as a reference dataset.
- Period: 1991-2010. Only 20 years have been selected for illustrative purposes, but the full hindcast period could be used for the analysis.
sdates <- paste0(1991:2010, '1201')
data <- Load('psl', exp = 'system4_m1',
obs = 'erainterim',
sdates = sdates,
storefreq='daily',
leadtimemin = 1, leadtimemax = 31,
latmin = 27, latmax = 81,
lonmin = 274.5, lonmax = 45,output = 'lonlat')
```
The object returned by `Load()` contains, among others, two arrays with the
requested data, one for the experimental data and the other for observational
```r
> dim(data$mod)
dataset member sdate ftime lat lon
> dim(data$obs)
dataset member sdate ftime lat lon
### 3- Daily anomalies based on a smoothed climatology
The weather regimes classification is based on daily anomalies, which have been computed by following these steps:
- Daily climatology is computed with the `Clim()` function.
```r
clim<-Clim(var_exp = data$mod,var_obs = data$obs,memb=F)
- The LOESS filter has been applied to the climatology to remove the short-term variability and retain the annual cycle. The parameter `loess_span` should be adapted to the period used to compute the climatology (e.g. season, month, week,...). In this example we are using monthly data, so we have selected `loess_span=1`.
Loess<-function(clim,loess_span){
data<-data.frame(ensmean=clim,day=1:length(clim))
loess_filt<-loess(ensmean~day,data,span=loess_span)
output<-predict(loess_filt)
return(output)
}
clim_smoothed_obs<-aperm(apply(clim$clim_obs,c(1:length(dim(clim$clim_obs)))[-which(names(dim(clim$clim_obs))=='ftime')],Loess,loess_span=1),c(2,1,3,4))
clim_smoothed_exp<-aperm(apply(clim$clim_exp,c(1:length(dim(clim$clim_exp)))[-which(names(dim(clim$clim_exp))=='ftime')],Loess,loess_span=1),c(2,1,3,4))
- Anomalies are computed as the deviations of the sea level pressure from the smoothed climatologies with the `Ano()` function.
anom_obs<-Ano(data$obs,clim_smoothed_obs)
anom_exp<-Ano(data$mod,clim_smoothed_exp)
```
### 4- Weather regimes in observations
`WeatherRegimes()` function is used to define the clusters based on the sea level pressure anomalies from ERA-Interim. This function is based on the [*kmeans function*] (http://stat.ethz.ch/R-manual/R-devel/library/stats/html/kmeans.html)
from the stats R package. In this example we have made different assumptions: four clusters (`ncenters=4`) will be produced and the Empirical orthogonal functions are not used to filter the data (`EOFS=FALSE`) just to take into account the extreme values. More details about the methodology can be found in Cortesi et al. 2018 (submitted).
WR_obs<-WeatherRegime(data=anom_obs,EOFS=FALSE,lat=data$lat,ncenters = 4)
`WeatherRegime()` provides a list with 6 elements (`$composite`, `$pvalue`, `$cluster`,`$center`, `$cluster_length`, `$cluster_values`, `$persistence` and `$frequency`) which are the needed parameters for the weather regimes assessment. Further details about the outputs provided by the `WeatherRegime()` function can be found in the package documentation ([*WeatherRegime documentation*] (https://earth.bsc.es/gitlab/es/s2dverification/blob/develop-Regimes/man/WeatherRegime.Rd)).
### 5- Visualisation of the observed weather regimes
To plot the composite maps of each regime and the mean frequencies of each cluster, we have employed the `PlotLayout()` and `PlotEquiMap()` functions available in s2dverifcation. The object `WR_obs$composite` is divided by 100 to change from Pa to hPa. As the `WR_obs$frequency` provides the monthly frequencies the climatological frequencies are obtained as the average across the 20 years of the monthly frequencies. Note that these frequencies could change as a consequence of the randomness inherent to the iterative processes involved in the k-means.
clim_frequencies<-paste0('freq = ',round(Mean1Dim(WR_obs$frequency,1),1),'%')
PlotLayout(PlotEquiMap,c(1,2),lon=data$lon,lat=data$lat,var=WR_obs$composite/100,
titles=paste0(paste0('Cluster ',1:4),' (',clim_frequencies,' )'),filled.continents=F,
axelab=F,draw_separators = T,subsampleg = 1,brks=seq(-16,16,by=2),
Verónica Torralba-Fernández
committed
bar_extra_labels = c(2,0,0,0),fileout='observed_regimes.png')
Verónica Torralba-Fernández
committed
<img src="observed_regimes.png" />
### 6- Weather regimes in the predictions
Predicted anomalies for each day, month, member and lead time are matched with the observed clusters (obtained in step 4). The assignment of the anomalies to a pre-defined set of clusters guarantees that the predicted weather regimes have very similar spatial structures to the observed regimes, which is an essential requirement for the verification of weather regimes. This is an example of how to produce a set of weather regimes based on the predictions that can be verified with the observational dataset, but this approach can be also used in an operational context for which the probability of occurence of each cluster could be estimated.
The observational pre-defined clusters are slightly adjusted to be used as a reference.
reference<-drop(WR_obs$composite)
names(dim(reference))<-c('lat','lon','nclust')
```
The matching is based on the minimization of Eucledian distance `method='distance'`, but it can also be also done in terms of spatial correlation `method='ACC'`. However the computational efficiency is superior for the distance method.
```r
WR_exp<-RegimesAssign(var_ano=anom_exp,ref_maps=reference,lats=data$lat,method='distance')
`RegimesAssign()` provides a list with 4 elements (`$composite`, `$pvalue`, `$cluster` and `$frequency`) with similar information to that provided by the `WeatherRegime()` function. Further details about the output can be found in the package documentation ([*RegimesAssign documentation*] (https://earth.bsc.es/gitlab/es/s2dverification/blob/develop-Regimes/man/RegimesAssign.Rd)).
### 7- Visualisation of the predicted weather regimes
Verónica Torralba-Fernández
committed
The outputs of `RegimesAssign()` have been represented to be compared with those obtained for the observational classification (step 5).
PlotLayout(PlotEquiMap,c(1,2),lon=data$lon,lat=data$lat,var=WR_exp$composite/100,
titles=paste0(paste0('Cluster ',1:4),' (',paste0('freq = ',round(WR_exp$frequency,1),'%'),' )'),
filled.continents=F,
axelab=F,draw_separators = T,subsampleg = 1,brks=seq(-16,16,by=2),
bar_extra_labels = c(2,0,0,0),fileout='predicted_regimes.png')
<img src="predicted_regimes.png" />