Skip to content
GitLab
Hide whitespace changes
Inline Side-by-side
vignettes/Figures/RainFARM_fig1.png
View replaced file @ 8173c9b3 View file @ 0450e139
vignettes/Figures/RainFARM_fig1.png

31 KB | W: | H:

vignettes/Figures/RainFARM_fig1.png

31.2 KB | W: | H:

vignettes/Figures/RainFARM_fig1.png
vignettes/Figures/RainFARM_fig1.png
vignettes/Figures/RainFARM_fig1.png
vignettes/Figures/RainFARM_fig1.png
  • 2-up
  • Swipe
  • Onion skin
vignettes/Figures/observed_regimes.png
View replaced file @ 8173c9b3 View file @ 0450e139
vignettes/Figures/observed_regimes.png

34.9 KB | W: | H:

vignettes/Figures/observed_regimes.png

39.7 KB | W: | H:

vignettes/Figures/observed_regimes.png
vignettes/Figures/observed_regimes.png
vignettes/Figures/observed_regimes.png
vignettes/Figures/observed_regimes.png
  • 2-up
  • Swipe
  • Onion skin
vignettes/Figures/predicted_regimes.png
View replaced file @ 8173c9b3 View file @ 0450e139
vignettes/Figures/predicted_regimes.png

34.1 KB | W: | H:

vignettes/Figures/predicted_regimes.png

39.4 KB | W: | H:

vignettes/Figures/predicted_regimes.png
vignettes/Figures/predicted_regimes.png
vignettes/Figures/predicted_regimes.png
vignettes/Figures/predicted_regimes.png
  • 2-up
  • Swipe
  • Onion skin
vignettes/MostLikelyTercile_vignette.Rmd
View file @ 0450e139
---
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}
......@@ -66,40 +68,84 @@ 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'
```
Finally, the data are loaded using `CST_Load`:
Finally, the data are loaded using `CST_Start`:
```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 = ""))
```
Loading the data using CST_Load 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:
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 member sdate ftime lat lon
1 10 40 3 19 36
dataset var member sdate ftime lat lon
1 1 10 40 3 19 36
> dim(obs$data)
dataset member sdate ftime lat lon
1 1 40 3 19 36
dataset var sdate ftime lat lon
1 1 40 3 19 36
```
......@@ -130,7 +176,7 @@ Finally, the probabilities of each tercile are computed by evaluating which terc
```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'))
prob_map <- MeanDims(PB, c('sdate', 'member', 'dataset', 'var'))
```
### 4. Visualization with PlotMostLikelyQuantileMap
......@@ -163,8 +209,10 @@ First, we evaluate and plot the RPSS. Therefore, we use `RPSS` metric included i
```r
Ano_Exp$data <- Subset(Ano_Exp$data, along = 'sdate', indices = 1:38)
Ano_Obs$data <- Subset(Ano_Obs$data, along = 'sdate', indices = 1:38)
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)
PlotEquiMap(RPSS$data[[1]], lat = Lat, lon = Lon, brks = seq(-1, 1, by = 0.1),
......@@ -190,7 +238,7 @@ Finally, we plot the latest forecast, as in the previous step, but add the mask
```r
PlotMostLikelyQuantileMap(probs = prob_map, lon = Lon, lat = Lat, coast_width = 1.5,
legend_scale = 0.5, mask = mask_rpss[ , , 1],
legend_scale = 0.5, mask = mask_rpss[ , , ,],
toptitle = paste('Most likely tercile -', clim_var,
'- ECMWF System5 - JJA 2020'),
width = 10, height = 8)
......
vignettes/MultiModelSkill_vignette.Rmd
View file @ 0450e139
......@@ -2,7 +2,7 @@
author: "Nuria Perez"
date: "`r Sys.Date()`"
revisor: "Eva Rifà"
revision date: "March 2023"
revision date: "October 2023"
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::knitr}
......@@ -20,6 +20,7 @@ The R package s2dv should be loaded by running:
```r
library(s2dv)
library(zeallot)
```
Library *CSTools*, should be installed from CRAN and loaded:
......@@ -51,35 +52,70 @@ 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%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.
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 `CST_Load` 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:
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
require(zeallot)
glosea5 <- '/esarchive/exp/glosea5/glosea5c3s/$STORE_FREQ$_mean/$VAR_NAME$_f6h/$VAR_NAME$_$YEAR$$MONTH$.nc'
c(exp, obs) %<-%
CST_Load(var = clim_var, exp = list(list(name = 'glosea5', path = glosea5),
list(name = 'ecmwf/system4_m1'),
list(name = 'meteofrance/system5_m1')),
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 = ""))
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:
......@@ -88,7 +124,6 @@ c(exp, obs) %<-%
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 `sd2v_cube` by running:
```
class(exp)
class(obs)
......@@ -98,10 +133,10 @@ The corresponding data is saved in the element `data` of each object, while othe
```r
> dim(exp$data)
dataset member sdate ftime lat lon
3 9 20 3 35 64
dataset var member sdate ftime lat lon
3 1 9 20 3 35 64
> dim(obs$data)
dataset member sdate ftime lat lon
dataset var sdate ftime lat lon
1 1 20 3 35 64
Lat <- exp$coords$lat
Lon <- exp$coords$lon
......@@ -121,15 +156,16 @@ The dimensions are preserved:
```
> str(ano_exp$data)
num [1:20, 1:3, 1:9, 1:3, 1:35, 1:64] -1.3958 -0.0484 -0.1326 0.3621 -5.6905 ...
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:3, 1:35, 1:64] 1.551 1.393 -0.344 -5.986 -0.27 ...
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)
```
......@@ -141,10 +177,10 @@ While other relevant data is being stored in the corresponding element of the ob
```r
> str(AnomDJF$data)
List of 4
$ corr : num [1:4, 1, 1:35, 1:64] 0.584 0.649 0.131 0.565 0.484 ...
$ p.val : num [1:4, 1, 1:35, 1:64] 0.0034 0.000989 0.291589 0.00475 0.015262 ...
$ conf.lower: num [1:4, 1, 1:35, 1:64] 0.192 0.289 -0.331 0.163 0.053 ...
$ conf.upper: num [1:4, 1, 1:35, 1:64] 0.816 0.848 0.542 0.806 0.763 ...
$ 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
......@@ -156,7 +192,7 @@ In the element $data of the `AnomDJF` object is a list of object for the metric
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(AnomDJF$data$corr[,1,,], lon = Lon, lat = Lat, map_select_fun = max,
PlotCombinedMap(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'),
......@@ -186,7 +222,7 @@ AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rms',
The following lines are necessary to obtain the plot which visualizes the best model given this metric for each grid point.
```r
PlotCombinedMap(AnomDJF$data$rms[,1,,], lon = Lon, lat = Lat, map_select_fun = min,
PlotCombinedMap(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'),
......@@ -210,9 +246,9 @@ Notice that the perfect RMSSS is 1 and the parameter `map_select_fun` from `Plo
```r
AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rmsss',
multimodel = TRUE)
multimodel = TRUE)
PlotCombinedMap(AnomDJF$data$rmsss[,1,,], lon = Lon, lat = Lat,
PlotCombinedMap(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',
......
vignettes/MultivarRMSE_vignette.Rmd
View file @ 0450e139
......@@ -2,7 +2,7 @@
author: "Deborah Verfaillie"
date: "`r Sys.Date()`"
revisor: "Eva Rifà"
revision date: "March 2023"
revision date: "October 2023"
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::knitr}
......@@ -20,6 +20,7 @@ To run this vignette, the next R packages should be installed and loaded:
```r
library(s2dv)
library(RColorBrewer)
library(zeallot)
```
Library *CSTools*, should be installed from CRAN and loaded:
......@@ -61,30 +62,114 @@ grid <- "256x128"
obs <- "erainterim"
```
Using the `CST_Load` function from **CSTool 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:
Using the `CST_Start` function from **CSTool 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 nuria.perez at bsc.es for the data to run the recipe.
```r
require(zeallot)
glosea5 <- '/esarchive/exp/glosea5/glosea5c3s/$STORE_FREQ$_mean/$VAR_NAME$_f6h/$VAR_NAME$_$YEAR$$MONTH$.nc'
c(exp_T, obs_T) %<-%
CST_Load(var = temp, exp = list(list(name = 'glosea5', path = glosea5)),
obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4,
latmin = 25, latmax = 75, lonmin = -20, lonmax = 70, output = 'lonlat',
nprocs = 1, storefreq = "monthly", sampleperiod = 1, nmember = 9,
method = "bilinear", grid = paste("r", grid, sep = ""))
glosea5 <- '/esarchive/exp/glosea5/glosea5c3s/$STORE_FREQ$_mean/$VAR_NAME$_f24h/$VAR_NAME$_$YEAR$$MONTH$.nc'
c(exp_P, obs_P) %<-%
CST_Load(var = precip, exp = list(list(name = 'glosea5', path = glosea5)),
obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4,
latmin = 25, latmax = 75, lonmin = -20, lonmax = 70, output = 'lonlat',
nprocs = 1, storefreq = "monthly", sampleperiod = 1, nmember = 9,
method = "bilinear", grid = paste("r", grid, sep = ""))
latmin = 25
latmax = 75
lonmin = -20
lonmax = 70
dateseq <- format(seq(start, end, by = "year"), "%Y%m")
repos1 <- "/esarchive/exp/glosea5/glosea5c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc"
exp_T <- CST_Start(dataset = list(list(name = 'glosea5c3s', path = repos1)),
var = temp,
member = startR::indices(1:9),
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 = NULL,
lon = NULL, ftime = 'sdate'),
retrieve = TRUE)
dates_exp <- exp_T$attrs$Dates
repos2 <- "/esarchive/recon/ecmwf/erainterim/monthly_mean/$var$/$var$_$date$.nc"
obs_T <- CST_Start(dataset = list(list(name = 'erainterim', path = repos2)),
var = temp,
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)
repos3 <- "/esarchive/exp/glosea5/glosea5c3s/monthly_mean/$var$_f24h/$var$_$sdate$.nc"
exp_P <- CST_Start(dataset = list(list(name = 'glosea5c3s', path = repos3)),
var = precip,
member = startR::indices(1:9),
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 = NULL,
lon = NULL, ftime = 'sdate'),
retrieve = TRUE)
dates_exp <- exp_P$attrs$Dates
obs_P <- CST_Start(dataset = list(list(name = 'erainterim', path = repos2)),
var = precip,
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_T, obs_T, exp_P, obs_P, file = "./tas_prlr_toydata.RData")
# Or use the following line to load the file provided in .RData format:
# load(file = "./tas_prlr_toydata.RData")
......@@ -105,11 +190,11 @@ Loading the data using `CST_Load` allows to obtain two lists, one for the experi
```
> dim(exp_T$data)
dataset member sdate ftime lat lon
1 9 20 3 35 64
dataset var member sdate ftime lat lon
1 1 9 20 3 35 64
> dim(obs_T$data)
dataset member sdate ftime lat lon
1 1 20 3 35 64
dataset var sdate ftime lat lon
1 1 20 3 35 64
```
Latitudes and longitudes of the common grid can be saved:
......@@ -131,9 +216,9 @@ The original dimensions are preserved and the anomalies are stored in the `data`
```
> str(ano_exp_T$data)
num [1:20, 1, 1:9, 1:3, 1:35, 1:64] -1.3958 -0.0484 -0.1326 0.3621 -5.6905 ...
num [1:20, 1, 1:9, 1, 1:3, 1:35, 1:64] -1.399 -0.046 -0.133 0.361 -5.696 ...
> str(ano_obs_T$data)
num [1:20, 1, 1, 1:3, 1:35, 1:64] 1.551 1.393 -0.344 -5.986 -0.27 ...
num [1:20, 1, 1, 1, 1:3, 1:35, 1:64] 1.556 1.397 -0.346 -5.99 -0.273 ...
```
Two lists containing the experiment ,`ano_exp`, and the observation, `ano_obs`, lists should be put together to serve as input of the function to compute multivariate RMSEs.
......@@ -154,6 +239,8 @@ It is obtained by running the `CST_MultivarRMSE` function:
```r
ano_obs[[1]] <- CST_InsertDim(ano_obs[[1]], posdim = 3, lendim = 1, name = "member")
ano_obs[[2]] <- CST_InsertDim(ano_obs[[2]], posdim = 3, lendim = 1, name = "member")
mvrmse <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight)
```
......@@ -163,51 +250,185 @@ The function `CST_MultivarRMSE` returns the multivariate RMSE value for 2 or mor
```r
> class(mvrmse)
> str(mvrmse$data)
num [1, 1, 1:35, 1:64] 1002261 1034354 1041180 1034907 1238147 ...
num [1, 1, 1, 1:35, 1:64] 806916 832753 838254 833206 996828 ...
> str(mvrmse$attrs$Variable)
List of 2
$ varName : chr [1:2] "tas" "prlr"
$ metadata:List of 6
..$ tas :List of 7
.. ..$ use_dictionary : logi FALSE
.. ..$ units : chr "K"
.. ..$ longname : chr "2 metre temperature"
.. ..$ description : chr "none"
.. ..$ daily_agg_cellfun : chr "none"
.. ..$ monthly_agg_cellfun: chr "none"
.. ..$ verification_time : chr "none"
..$ lon : num [1:64(1d)] 0 1.41 2.81 4.22 5.62 ...
.. ..- attr(*, "cdo_grid_name")= chr "r256x128"
.. ..- attr(*, "data_across_gw")= logi TRUE
.. ..- attr(*, "array_across_gw")= logi FALSE
.. ..- attr(*, "first_lon")= num 340
.. ..- attr(*, "last_lon")= num 68.9
.. ..- attr(*, "projection")= chr "none"
..$ lat : num [1:35(1d)] 73.8 72.4 71 69.6 68.2 ...
.. ..- attr(*, "cdo_grid_name")= chr "r256x128"
.. ..- attr(*, "first_lat")= num [1(1d)] 26
.. ..- attr(*, "last_lat")= num [1(1d)] 73.8
.. ..- attr(*, "projection")= chr "none"
..$ prlr:List of 7
.. ..$ use_dictionary : logi FALSE
.. ..$ units : chr "m s-1"
.. ..$ longname : chr "Total precipitation"
.. ..$ description : chr "none"
.. ..$ daily_agg_cellfun : chr "none"
.. ..$ monthly_agg_cellfun: chr "none"
.. ..$ verification_time : chr "none"
..$ lon : num [1:64(1d)] 0 1.41 2.81 4.22 5.62 ...
.. ..- attr(*, "cdo_grid_name")= chr "r256x128"
.. ..- attr(*, "data_across_gw")= logi TRUE
.. ..- attr(*, "array_across_gw")= logi FALSE
.. ..- attr(*, "first_lon")= num 340
.. ..- attr(*, "last_lon")= num 68.9
.. ..- attr(*, "projection")= chr "none"
..$ lat : num [1:35(1d)] 73.8 72.4 71 69.6 68.2 ...
.. ..- attr(*, "cdo_grid_name")= chr "r256x128"
.. ..- attr(*, "first_lat")= num [1(1d)] 26
.. ..- attr(*, "last_lat")= num [1(1d)] 73.8
.. ..- attr(*, "projection")= chr "none"
$ metadata:List of 8
..$ lat : num [1:35(1d)] 73.8 72.4 71 69.6 68.2 ...
..$ lon : num [1:64(1d)] 0 1.41 2.81 4.22 5.62 ...
..$ ftime: POSIXct[1:60], format: "1993-12-16 00:00:00" "1994-12-16 00:00:00" ...
..$ tas :List of 11
.. ..$ prec : chr "float"
.. ..$ units : chr "K"
.. ..$ dim :List of 4
.. .. ..$ :List of 10
.. .. .. ..$ name : chr "lon"
.. .. .. ..$ len : int 64
.. .. .. ..$ unlim : logi FALSE
.. .. .. ..$ group_index : int 1
.. .. .. ..$ group_id : int 65536
.. .. .. ..$ id : int 1
.. .. .. ..$ dimvarid :List of 5
.. .. .. .. ..$ id : int 1
.. .. .. .. ..$ group_index: int 1
.. .. .. .. ..$ group_id : int 65536
.. .. .. .. ..$ list_index : num -1
.. .. .. .. ..$ isdimvar : logi TRUE
.. .. .. .. ..- attr(*, "class")= chr "ncid4"
.. .. .. ..$ units : chr "degrees_east"
.. .. .. ..$ vals : num [1:64(1d)] -19.7 -18.3 -16.9 -15.5 -14.1 ...
.. .. .. ..$ create_dimvar: logi TRUE
.. .. .. ..- attr(*, "class")= chr "ncdim4"
.. .. ..$ :List of 10
.. .. .. ..$ name : chr "lat"
.. .. .. ..$ len : int 35
.. .. .. ..$ unlim : logi FALSE
.. .. .. ..$ group_index : int 1
.. .. .. ..$ group_id : int 65536
.. .. .. ..$ id : int 0
.. .. .. ..$ dimvarid :List of 5
.. .. .. .. ..$ id : int 0
.. .. .. .. ..$ group_index: int 1
.. .. .. .. ..$ group_id : int 65536
.. .. .. .. ..$ list_index : num -1
.. .. .. .. ..$ isdimvar : logi TRUE
.. .. .. .. ..- attr(*, "class")= chr "ncid4"
.. .. .. ..$ units : chr "degrees_north"
.. .. .. ..$ vals : num [1:35(1d)] 26 27.4 28.8 30.2 31.6 ...
.. .. .. ..$ create_dimvar: logi TRUE
.. .. .. ..- attr(*, "class")= chr "ncdim4"
.. .. ..$ :List of 10
.. .. .. ..$ name : chr "ensemble"
.. .. .. ..$ len : int 14
.. .. .. ..$ unlim : logi FALSE
.. .. .. ..$ group_index : int 1
.. .. .. ..$ group_id : int 65536
.. .. .. ..$ id : int 3
.. .. .. ..$ dimvarid :List of 5
.. .. .. .. ..$ id : int -1
.. .. .. .. ..$ group_index: int 1
.. .. .. .. ..$ group_id : int 65536
.. .. .. .. ..$ list_index : num -1
.. .. .. .. ..$ isdimvar : logi TRUE
.. .. .. .. ..- attr(*, "class")= chr "ncid4"
.. .. .. ..$ vals : int [1:14] 1 2 3 4 5 6 7 8 9 10 ...
.. .. .. ..$ units : chr ""
.. .. .. ..$ create_dimvar: logi FALSE
.. .. .. ..- attr(*, "class")= chr "ncdim4"
.. .. ..$ :List of 11
.. .. .. ..$ name : chr "time"
.. .. .. ..$ len : int 7
.. .. .. ..$ unlim : logi TRUE
.. .. .. ..$ group_index : int 1
.. .. .. ..$ group_id : int 65536
.. .. .. ..$ id : int 2
.. .. .. ..$ dimvarid :List of 5
.. .. .. .. ..$ id : int 2
.. .. .. .. ..$ group_index: int 1
.. .. .. .. ..$ group_id : int 65536
.. .. .. .. ..$ list_index : num -1
.. .. .. .. ..$ isdimvar : logi TRUE
.. .. .. .. ..- attr(*, "class")= chr "ncid4"
.. .. .. ..$ units : chr "months since 1993-11-15 12:00:00"
.. .. .. ..$ calendar : chr "proleptic_gregorian"
.. .. .. ..$ vals : num [1:7(1d)] 0 1 2 3 4 5 6
.. .. .. ..$ create_dimvar: logi TRUE
.. .. .. ..- attr(*, "class")= chr "ncdim4"
.. ..$ unlim : logi TRUE
.. ..$ make_missing_value: logi TRUE
.. ..$ missval : num -9e+33
.. ..$ hasAddOffset : logi FALSE
.. ..$ hasScaleFact : logi FALSE
.. ..$ table : int 128
.. ..$ _FillValue : num -9e+33
.. ..$ missing_value : num -9e+33
..$ lat : num [1:35(1d)] 73.8 72.4 71 69.6 68.2 ...
..$ lon : num [1:64(1d)] 0 1.41 2.81 4.22 5.62 ...
..$ ftime: POSIXct[1:60], format: "1993-12-16 00:00:00" "1994-12-16 00:00:00" ...
..$ prlr :List of 9
.. ..$ prec : chr "float"
.. ..$ units : chr "m s-1"
.. ..$ dim :List of 4
.. .. ..$ :List of 10
.. .. .. ..$ name : chr "lon"
.. .. .. ..$ len : int 64
.. .. .. ..$ unlim : logi FALSE
.. .. .. ..$ group_index : int 1
.. .. .. ..$ group_id : int 65536
.. .. .. ..$ id : int 1
.. .. .. ..$ dimvarid :List of 5
.. .. .. .. ..$ id : int 1
.. .. .. .. ..$ group_index: int 1
.. .. .. .. ..$ group_id : int 65536
.. .. .. .. ..$ list_index : num -1
.. .. .. .. ..$ isdimvar : logi TRUE
.. .. .. .. ..- attr(*, "class")= chr "ncid4"
.. .. .. ..$ units : chr "degrees_east"
.. .. .. ..$ vals : num [1:64(1d)] -19.7 -18.3 -16.9 -15.5 -14.1 ...
.. .. .. ..$ create_dimvar: logi TRUE
.. .. .. ..- attr(*, "class")= chr "ncdim4"
.. .. ..$ :List of 10
.. .. .. ..$ name : chr "lat"
.. .. .. ..$ len : int 35
.. .. .. ..$ unlim : logi FALSE
.. .. .. ..$ group_index : int 1
.. .. .. ..$ group_id : int 65536
.. .. .. ..$ id : int 0
.. .. .. ..$ dimvarid :List of 5
.. .. .. .. ..$ id : int 0
.. .. .. .. ..$ group_index: int 1
.. .. .. .. ..$ group_id : int 65536
.. .. .. .. ..$ list_index : num -1
.. .. .. .. ..$ isdimvar : logi TRUE
.. .. .. .. ..- attr(*, "class")= chr "ncid4"
.. .. .. ..$ units : chr "degrees_north"
.. .. .. ..$ vals : num [1:35(1d)] 26 27.4 28.8 30.2 31.6 ...
.. .. .. ..$ create_dimvar: logi TRUE
.. .. .. ..- attr(*, "class")= chr "ncdim4"
.. .. ..$ :List of 10
.. .. .. ..$ name : chr "ensemble"
.. .. .. ..$ len : int 14
.. .. .. ..$ unlim : logi FALSE
.. .. .. ..$ group_index : int 1
.. .. .. ..$ group_id : int 65536
.. .. .. ..$ id : int 3
.. .. .. ..$ dimvarid :List of 5
.. .. .. .. ..$ id : int -1
.. .. .. .. ..$ group_index: int 1
.. .. .. .. ..$ group_id : int 65536
.. .. .. .. ..$ list_index : num -1
.. .. .. .. ..$ isdimvar : logi TRUE
.. .. .. .. ..- attr(*, "class")= chr "ncid4"
.. .. .. ..$ vals : int [1:14] 1 2 3 4 5 6 7 8 9 10 ...
.. .. .. ..$ units : chr ""
.. .. .. ..$ create_dimvar: logi FALSE
.. .. .. ..- attr(*, "class")= chr "ncdim4"
.. .. ..$ :List of 11
.. .. .. ..$ name : chr "time"
.. .. .. ..$ len : int 7
.. .. .. ..$ unlim : logi TRUE
.. .. .. ..$ group_index : int 1
.. .. .. ..$ group_id : int 65536
.. .. .. ..$ id : int 2
.. .. .. ..$ dimvarid :List of 5
.. .. .. .. ..$ id : int 2
.. .. .. .. ..$ group_index: int 1
.. .. .. .. ..$ group_id : int 65536
.. .. .. .. ..$ list_index : num -1
.. .. .. .. ..$ isdimvar : logi TRUE
.. .. .. .. ..- attr(*, "class")= chr "ncid4"
.. .. .. ..$ units : chr "months since 1993-11-15 12:00:00"
.. .. .. ..$ calendar : chr "proleptic_gregorian"
.. .. .. ..$ vals : num [1:7(1d)] 0 1 2 3 4 5 6
.. .. .. ..$ create_dimvar: logi TRUE
.. .. .. ..- attr(*, "class")= chr "ncdim4"
.. ..$ unlim : logi TRUE
.. ..$ make_missing_value: logi FALSE
.. ..$ missval : num 1e+30
.. ..$ hasAddOffset : logi FALSE
.. ..$ hasScaleFact : logi FALSE
.. ..$ table : int 128
```
The following lines plot the multivariate RMSE
......
vignettes/PlotForecastPDF.Rmd
View file @ 0450e139
......@@ -18,7 +18,7 @@ Library *CSTools*, should be installed from CRAN and loaded:
library(CSTools)
```
### 1.- A simple example
### 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 PlotForecastPDF by default will plot the ensemble members, the estimated density distributions and the tercile probabilities.
......@@ -37,7 +37,7 @@ fcst <- array(rnorm(mean = 25, sd = 2, n = 90), dim = c(member = 30, 3))
PlotForecastPDF(fcst, tercile.limits = c(23, 27))
```
### 2.- Customizing the appearance of your plots
### 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}
......@@ -50,7 +50,7 @@ PlotForecastPDF(fcst, tercile.limits = c(20, 26), var.name = "Temperature (ºC)"
```
![Example 2](./Figures/PlotForecastPDF_ex2.png)
### 3.- Adding extremes and observed values
### 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}
......@@ -65,7 +65,7 @@ PlotForecastPDF(fcst, tercile.limits = rbind(c(20, 26), c(22, 28), c(15, 22)),
![Example 3](./Figures/PlotForecastPDF_ex3.png)
### 4.- Saving your plot to a file
### 4. Saving your plot to a file
PlotForecastPDF 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:
```
......@@ -75,12 +75,12 @@ plot <- PlotForecastPDF(fcst, tercile.limits = c(23, 27))
ggsave("outfile.pdf", plot, width = 7, height = 5)
```
### 5.- A reproducible example using lonlat_temp
### 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$exp$data[1,,1,1,1,1] - 273.15,
fcst2 = lonlat_temp$exp$data[1,,1,2,1,1] - 273.15)
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)
PlotForecastPDF(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",
......
vignettes/RainFARM_vignette.Rmd
View file @ 0450e139
......@@ -2,6 +2,8 @@
title: "Rainfall Filtered Autoregressive Model (RainFARM) precipitation downscaling"
author: "Jost von Hardenberg (ISAC-CNR)"
date: "26/03/2019"
revisor: "Eva Rifà"
revision date: "October 2023"
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::knitr}
......@@ -38,23 +40,22 @@ library(CSTools)
We use test data provided by CSTools to load a seasonal precipitation forecast:
```{r}
library(CSTools)
exp <- lonlat_prec
exp <- lonlat_prec_st
```
This gives us a CSTools object `exp`, containing an element `exp$data` with dimensions:
```{r}
dim(exp$data)
# dataset member sdate ftime lat lon
# 1 6 3 31 4 4
# dataset var member sdate ftime lat lon
# 1 1 6 3 31 4 4
```
There are 6 ensemble members available in the data set, 3 starting dates and 31 forecast times, which refer to daily values in the month of March following starting dates on November 1st in the years 2010, 2011, 2012. Please notice that RainFARM (in this version) only accepts square domains, possibly with an even number of pixels on each side, so we always need to select an appropriate cutout. Also, there are time and memory limitations when a large ensemble of downscaled realizations is generated with RainFARM, so that selecting a smaller target area is advised. On the other hand, if spectral slopes are to be determined from the large scales we will still need enough resolution to allow this estimation. In this example we have preselected a 4x4 pixel cutout at resolution 1 degree in a smaller area lon=[6,9], lat=[44,47] covering Northern Italy.
```{r}
ilon <- which(exp$coords$lon %in% 5:12)
ilat <- which(exp$coords$lat %in% 40:47 )
exp$data <- exp$data[ , , , , ilon, ilat, drop = FALSE]
names(dim(exp$data)) <- names(dim(lonlat_prec$data))
ilat <- which(exp$coords$lat %in% 40:47)
exp$data <- exp$data[, , , , , ilat, ilon, drop = FALSE]
names(dim(exp$data)) <- names(dim(lonlat_prec_st$data))
exp$coords$lon <- exp$coords$lon[ilon]
exp$coords$lat <- exp$coords$lat[ilat]
```
......@@ -65,13 +66,14 @@ Our goal is to downscale with RainFARM these data from the resolution of 1 degre
RainFARM can compute automatically its only free parameter, i.e. the spatial spectral slope, from the large-scale field (here only with size 4x4 pixel, but in general we reccomend selecting at least 8x8 pixels).
In this example we would like to compute this slope as an average over the _member_ and _ftime_ dimensions, while we will use different slopes for the remaining _dataset_ and _sdate_ dimensions (a different choice may be more appropriate in a real application). To obtain this we specify the parameter `time_dim = c("member", "ftime")`. The slope is computed starting from the wavenumber corresponding to the box, `kmin = 1`. We create 3 stochastic realizations for each dataset, member, starting date and forecast time with `nens = 5`. The command to donwscale and the resulting fields are:
```{r}
```
exp_down <- CST_RainFARM(exp, nf = 20, kmin = 1, nens = 3,
time_dim = c("member", "ftime"))
dim(exp_down$data)
# dataset member realization sdate ftime lat lon
# 1 6 3 3 31 80 80
# dataset var member realization sdate ftime lat lon
# 1 1 6 3 3 31 80 80
str(exp_down$coords$lon)
# num [1:80] 5.53 5.58 5.62 5.67 5.72 ...
str(exp_down$coords$lat)
......@@ -95,13 +97,13 @@ png("Figures/RainFARM_fig1.png", width = 640, height = 365)
par(mfrow = c(1,2))
-->
```{r}
a <- exp$data[1, 1, 1, 17, , ] * 86400 * 1000
a <- exp$data[1, 1, 1, 1, 17, , ] * 86400 * 1000
a[a > 60] <- 60
image(exp$coords$lon, rev(exp$coords$lat), t(apply(a, 2, rev)), xlab = "lon", ylab = "lat",
col = rev(terrain.colors(20)), zlim = c(0,60))
map("world", add = TRUE)
title(main = "pr 17/03/2010 original")
a <- exp_down$data[1, 1, 1, 1, 17, , ] * 86400 * 1000
a <- exp_down$data[1, 1, 1, 1, 1, 17, , ] * 86400 * 1000
a[a > 60] <- 60
image(exp_down$coords$lon, rev(exp_down$coords$lat), t(apply(a, 2, rev)), xlab = "lon", ylab = "lat",
col = rev(terrain.colors(20)), zlim = c(0, 60))
......@@ -137,8 +139,8 @@ From a single realization and time it is not possible to see that a more realist
<!--
```{r}
exp_down1 <- exp_down$data[1, , , 1, , , ]
exp_down_weights1 <- exp_down_weights$data[1, , , 1, , , ]
exp_down1 <- exp_down$data[, , , , , , , 1]
exp_down_weights1 <- exp_down_weights$data[, , , , , , , 1]
dim(exp_down1) <- c(member = 6 * 3 * 31, lat = 80, lon = 80)
dim(exp_down_weights1) <- c(member = 6 * 3 * 31, lat = 80, lon = 80)
ad <- apply(exp_down1, c(2, 3), mean)
......@@ -146,7 +148,7 @@ adw <- apply(exp_down_weights1, c(2, 3), mean);
png("Figures/RainFARM_fig2.png", width = 640, height = 243)
par(mfrow = c(1,3))
a <- exp_down_weights$data[1, 1, 1, 1, 17, , ] * 86400 * 1000
a <- exp_down_weights$data[1, 1, 1, 1, 17, , ,1] * 86400 * 1000
a[a > 60] <- 60
image(exp_down$coords$lon, rev(exp_down$coords$lat), t(apply(a, 2, rev)), xlab = "lon",
ylab = "lat", col = rev(terrain.colors(20)), zlim = c(0, 60))
......@@ -176,11 +178,23 @@ The only free parameter of the RainFARM method is represented by the spatial spe
```{r}
slopes <- CST_RFSlope(exp, time_dim = c("member", "ftime"))
dim(slopes)
# dataset sdate
# 1 3
slopes
[,1] [,2] [,3]
[1,] 1.09958 1.768862 1.190185
# dataset var sdate
# 1 1 3
# slopes
# , , 1
# [,1]
# [1,] 1.09957
# , , 2
# [,1]
# [1,] 1.768861
# , , 3
# [,1]
# [1,] 1.190176
```
which return an array of spectral slopes, one for each "dataset" and starting date "sdate".
......
vignettes/WeatherRegimes_vignette.Rmd
View file @ 0450e139
---
author: "Verónica Torralba, Nicola Cortesi 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}
......@@ -35,37 +37,78 @@ The data employed in this example are described below.
```r
repos_exp <- paste0('/esarchive/exp/ecmwf/system4_m1/daily_mean/$var$_f6h/',
'$var$_$sdate$.nc')
sdates <- paste0(1991:2010, '1201')
c(exp, obs) %<-% CST_Load(var = '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')
latmin = 27
latmax = 81
lonmin = 274.5
lonmax = 45
exp <- CST_Start(dataset = repos_exp,
var = 'psl',
member = "all",
sdate = sdates,
ftime = startR::indices(1:31),
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')),
return_vars = list(lat = NULL,
lon = NULL, ftime = 'sdate'),
retrieve = TRUE)
dates0 <- exp$attrs$Dates
repos_obs <-"/esarchive/recon/ecmwf/erainterim/daily_mean/$var$_f6h/$var$_$date$.nc"
obs <- CST_Start(dataset = repos_obs,
var = 'psl',
date = unique(format(dates0, '%Y%m')),
ftime = startR::values(dates0),
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')),
return_vars = list(lon = NULL,
lat = NULL,
ftime = 'date'),
retrieve = TRUE)
```
Notice that you need the files to be stored locally in your computer or server with correct configuration file. If you are interested into run this vignette, contact nuria.perez at bsc.es to get a data sample.
The objects returned by `CST_Load()` are s2v_cube class. They contains among others, the array with the requested data.
The objects returned by `CST_Start()` are `s2v_cube` class. They contains among others, the array with the requested data.
```r
> dim(exp$data)
dataset member sdate ftime lat lon
1 15 20 31 77 186
dataset var member sdate ftime lat lon
1 1 15 20 31 77 186
> dim(obs$data)
dataset member sdate ftime lat lon
1 1 20 31 77 186
dataset var sdate ftime lat lon
1 1 20 31 77 186
```
### 3- Daily anomalies based on a smoothed climatology)
### 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:
```r
c(ano_exp, ano_obs) %<-% CST_Anomaly(exp = exp, obs = obs, filter_span = 1)
c(ano_exp, ano_obs) %<-% CST_Anomaly(exp = exp, obs = obs, filter_span = 1,
dat_dim = c('dataset', 'member'))
```
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 length of the period used to compute the climatology (e.g. season, month, week,...). In this example we are using daily data, so we have selected `loess_span = 1`.
The LOESS filter has been applied to the climatology to remove the short-term variability and retain the annual cycle. The parameter `filter_span` should be adapted to the length of the period used to compute the climatology (e.g. season, month, week,...). In this example we are using daily data, so we have selected `filter_span = 1`.
### 4- Weather regimes in observations
......@@ -75,23 +118,23 @@ from the stats R package. In this example we have made different assumptions: fo
```r
WR_obs <- CST_WeatherRegimes(data = ano_obs, EOFs = FALSE, ncenters = 4)
```
```
`CST_WeatherRegime()` provides a s2dv_cube object with several elements. `$data` the 4 weather regimes composites are stored while `$statistics` contains extra information (`$pvalue`, `$cluster`, `$persistence` and `$frequency`) which are the needed parameters for the weather regimes assessment. Further details about the outputs provided by the `CST_WeatherRegime()` function can be found in the package documentation or typing `?CST_WeatherRegimes` in the R session.
`CST_WeatherRegime()` provides a `s2dv_cube` object with several elements. `$data` the 4 weather regimes composites are stored while `$statistics` contains extra information (`$pvalue`, `$cluster`, `$persistence` and `$frequency`) which are the needed parameters for the weather regimes assessment. Further details about the outputs provided by the `CST_WeatherRegime()` function can be found in the package documentation or typing `?CST_WeatherRegimes` in the R session.
### 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 s2dv. The object `WR_obs$data` is divided by 100 to change from Pa to hPa. As the `WR_obs$statistics$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 slightly change as a consequence of the randomness inherent to the iterative processes involved in the k-means.
```r
clim_frequencies <- paste0('freq = ', round(Mean1Dim(WR_obs$statistics$frequency, 1), 1), '%')
```r
clim_frequencies <- paste0('freq = ', round(MeanDims(WR_obs$statistics$frequency, 1), 1), '%')
PlotLayout(PlotEquiMap, c(1, 2), lon = obs$coords$lon, lat = obs$coords$lat,
var = WR_obs$data / 100,
titles = paste0(paste0('Cluster ', 1:4), ' (', clim_frequencies,' )'),
filled.continents = FALSE,
axelab = FALSE, draw_separators = TRUE, subsampleg = 1,
brks = seq(-16, 16, by = 2),
bar_extra_labels = c(2, 0, 0, 0), fileout = './Figures/observed_regimes.png')
bar_extra_labels = c(2, 0, 0, 0))
```
<img src="./Figures/observed_regimes.png" />
......@@ -109,10 +152,8 @@ PlotTriangles4Categories(freq_obs, toptitle = 'Persistence',
xtitle = 'Start Dates', ytitle = '', xlab = FALSE,
ylabels = substr(sdates, 1, 4), cex_leg = 0.6,
lab_legend = c('AR', 'NAO-', 'BL', 'NAO+'), figure.width = .7)
```
<img src="./Figures/Obs_Persistence.png" width = "400px"/>
### 7- Weather regimes in the predictions
......@@ -121,30 +162,27 @@ Predicted anomalies for each day, month, member and lead time are matched with t
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 <- CST_RegimesAssign(data = ano_exp, ref_maps = WR_obs, method = 'distance', composite = TRUE, memb = TRUE)
WR_exp <- CST_RegimesAssign(data = ano_exp, ref_maps = WR_obs,
method = 'distance', composite = TRUE, memb = TRUE)
```
`CST_RegimesAssign()` provides a object of 's2dv_cube' class which lists several elements. The composites are stored in the $data element which in the case of `memb = TRUE` corresponds to the mean composites for all member. `$pvalue`, `$cluster` and `$frequency` are stored in element `$statistics. This elements contain similar information to that provided by the `WeatherRegime()` function. Further details about the output can be found in the package documentation.
`CST_RegimesAssign()` provides a object of 's2dv_cube' class which lists several elements. The composites are stored in the $data element which in the case of `memb = TRUE` corresponds to the mean composites for all member. `$pvalue`, `$cluster` and `$frequency` are stored in element `$statistics`. This elements contain similar information to that provided by the `WeatherRegime()` function. Further details about the output can be found in the package documentation.
### 8- Visualisation of the predicted weather regimes
The outputs of `RegimesAssign()` have been represented to be compared with those obtained for the observational classification (step 5).
```r
PlotLayout(PlotEquiMap, c(1, 2),lon = exp$coords$lon, lat = exp$coords$lat,
PlotLayout(PlotEquiMap, c(1, 2), lon = exp$coords$lon, lat = exp$coords$lat,
var = WR_exp$data/100,
titles = paste0(paste0('Cluster ',1:4), ' (',paste0('freq = ',
round(WR_exp$statistics$frequency,1),'%'),' )'),
filled.continents = FALSE,
axelab = FALSE, draw_separators = TRUE,
subsampleg = 1, brks = seq(-16, 16, by = 2),
bar_extra_labels = c(2, 0, 0, 0), fileout = './Figures/predicted_regimes.png')
bar_extra_labels = c(2, 0, 0, 0))
```
<img src="./Figures/predicted_regimes.png" />
Observed and predicted weather regimes are very similar although their frequencies are slightly different. Cluster 1 is the Atlantic Ridge and cluster 3 the Blocking pattern, while cluster 4 and 2 are the positive and negative phases of the NAO. This patterns can change depending on the period analyzed.
Observed and predicted weather regimes are very similar although their frequencies are slightly different. Cluster 4 is the Atlantic Ridge and cluster 2 the Blocking pattern, while cluster 3 and 1 are the positive and negative phases of the NAO. This patterns can change depending on the period analyzed.