diff --git a/R/sample_data_st.R b/R/sample_data_st.R index e276a8c9fce9d12f246a7399b2c98b07131ce015..cd33ee069d270640c08a27c0a8dadb4b2fb4c0b2 100644 --- a/R/sample_data_st.R +++ b/R/sample_data_st.R @@ -23,7 +23,7 @@ #' lonmin <- -12 #' latmax <- 48 #' latmin <- 27 -#' lonlat_temp_st$exp <- CST_Start(dat = repos_exp, +#' lonlat_temp_st$exp <- CST_Start(dataset = repos_exp, #' var = 'tas', #' member = indices(1:15), #' sdate = sdates, @@ -39,15 +39,19 @@ #' return_vars = list(lat = NULL, #' lon = NULL, ftime = 'sdate'), #' retrieve = TRUE) +#' #' dates <- c(paste0(2000, c(11, 12)), paste0(2001, c('01', 11, 12)), #' paste0(2002, c('01', 11, 12)), paste0(2003, c('01', 11, 12)), #' paste0(2004, c('01', 11, 12)), paste0(2005, c('01', 11, 12)), 200601) #' dates <- sapply(dates, function(x) {paste0(x, '01')}) #' dates <- as.POSIXct(dates, format = '%Y%m%d', 'UTC') #' dim(dates) <- c(ftime = 3, sdate = 6) +#' +#' dates <- t(dates) +#' names(dim(dates)) <- c('sdate', 'ftime') #' #' path.obs <- '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$var$_$date$.nc' -#' lonlat_temp_st$obs <- CST_Start(dat = path.obs, +#' lonlat_temp_st$obs <- CST_Start(dataset = path.obs, #' var = 'tas', #' date = unique(format(dates, '%Y%m')), #' ftime = values(dates), @@ -71,6 +75,16 @@ #' lat = NULL, #' ftime = 'date'), #' retrieve = TRUE) +#' +#' library(lubridate) +#' dates_exp <- lonlat_temp_st$exp$attrs$Dates +#' lonlat_temp_st$exp$attrs$Dates <- floor_date(ymd_hms(dates_exp), unit = "months") +#' dim(lonlat_temp_st$exp$attrs$Dates) <- dim(dates_exp) +#' +#' dates_obs <- lonlat_temp_st$obs$attrs$Dates +#' lonlat_temp_st$obs$attrs$Dates <- floor_date(ymd_hms(dates_obs), unit = "months") +#' dim(lonlat_temp_st$obs$attrs$Dates) <- dim(dates_obs) +#' #'} #' #'@name lonlat_temp_st @@ -103,7 +117,7 @@ NULL #' lonmin <- 6 #' lonmax <- 9 #' -#' lonlat_prec_st <- CST_Start(dat = path, +#' lonlat_prec_st <- CST_Start(dataset = path, #' var = 'prlr', #' member = indices(1:6), #' sdate = sdates, diff --git a/data/lonlat_prec_st.rda b/data/lonlat_prec_st.rda index 405bd2297ddfac596eb65b96458b2f8f14df1975..02e8e28453a1fbbf3c3019b3f7b4fd6b43c28dd7 100644 Binary files a/data/lonlat_prec_st.rda and b/data/lonlat_prec_st.rda differ diff --git a/data/lonlat_temp_st.rda b/data/lonlat_temp_st.rda index b9f8c2bf1116adfc67fc7bef518250b7d6575b82..cad51c5546cb13c7dd3b556890388f93959e74bb 100644 Binary files a/data/lonlat_temp_st.rda and b/data/lonlat_temp_st.rda differ diff --git a/man/lonlat_prec_st.Rd b/man/lonlat_prec_st.Rd index 3f0de74af5d14b385170e6b6ce788236840d984e..817156319de7f4fb82522e78d807b4ae3397f0d4 100644 --- a/man/lonlat_prec_st.Rd +++ b/man/lonlat_prec_st.Rd @@ -28,7 +28,7 @@ extracted from forecast data available in the MEDSCOPE internal archive. lonmin <- 6 lonmax <- 9 - lonlat_prec_st <- CST_Start(dat = path, + lonlat_prec_st <- CST_Start(dataset = path, var = 'prlr', member = indices(1:6), sdate = sdates, diff --git a/man/lonlat_temp_st.Rd b/man/lonlat_temp_st.Rd index 3ac3e048eda8e1c251eaea29c7d28a4989ff316e..8dc554f2656a7d7d967d3683b66a49eb78dc2072 100644 --- a/man/lonlat_temp_st.Rd +++ b/man/lonlat_temp_st.Rd @@ -29,7 +29,7 @@ next. Note that `CST_Start` internally calls `startR::Start` and then uses lonmin <- -12 latmax <- 48 latmin <- 27 - lonlat_temp_st$exp <- CST_Start(dat = repos_exp, + lonlat_temp_st$exp <- CST_Start(dataset = repos_exp, var = 'tas', member = indices(1:15), sdate = sdates, @@ -45,6 +45,7 @@ next. Note that `CST_Start` internally calls `startR::Start` and then uses return_vars = list(lat = NULL, lon = NULL, ftime = 'sdate'), retrieve = TRUE) + dates <- c(paste0(2000, c(11, 12)), paste0(2001, c('01', 11, 12)), paste0(2002, c('01', 11, 12)), paste0(2003, c('01', 11, 12)), paste0(2004, c('01', 11, 12)), paste0(2005, c('01', 11, 12)), 200601) @@ -52,8 +53,11 @@ next. Note that `CST_Start` internally calls `startR::Start` and then uses dates <- as.POSIXct(dates, format = '%Y%m%d', 'UTC') dim(dates) <- c(ftime = 3, sdate = 6) + dates <- t(dates) + names(dim(dates)) <- c('sdate', 'ftime') + path.obs <- '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$var$_$date$.nc' - lonlat_temp_st$obs <- CST_Start(dat = path.obs, + lonlat_temp_st$obs <- CST_Start(dataset = path.obs, var = 'tas', date = unique(format(dates, '%Y%m')), ftime = values(dates), @@ -77,6 +81,16 @@ next. Note that `CST_Start` internally calls `startR::Start` and then uses lat = NULL, ftime = 'date'), retrieve = TRUE) + + library(lubridate) + dates_exp <- lonlat_temp_st$exp$attrs$Dates + lonlat_temp_st$exp$attrs$Dates <- floor_date(ymd_hms(dates_exp), unit = "months") + dim(lonlat_temp_st$exp$attrs$Dates) <- dim(dates_exp) + + dates_obs <- lonlat_temp_st$obs$attrs$Dates + lonlat_temp_st$obs$attrs$Dates <- floor_date(ymd_hms(dates_obs), unit = "months") + dim(lonlat_temp_st$obs$attrs$Dates) <- dim(dates_obs) + } } \author{ diff --git a/vignettes/Analogs_vignette.Rmd b/vignettes/Analogs_vignette.Rmd index 674dccac5e42d414b1b611b83c705fa06d4c39ab..4a8953f6819f0d76cedc3f696dd12f2d8b87f2a3 100644 --- a/vignettes/Analogs_vignette.Rmd +++ b/vignettes/Analogs_vignette.Rmd @@ -3,7 +3,7 @@ title: "Analogs based on large scale for downscaling" author: "M. Carmen Alvarez-Castro and M. del Mar Chaves-Montero (CMCC, Italy)" date: "November 2020" revisor: "Eva Rifà" -revision date: "January 2023" +revision date: "October 2023" output: rmarkdown::html_vignette vignette: > %\VignetteEngine{knitr::knitr} @@ -21,8 +21,7 @@ In this example, the seasonal temperature forecasts, initialized in october, wil ## 1. Introduction of the function -For instance if we want to perform a temperature donwscaling in Balearic Island for October we will get a daily series of temperature with 1 analog per day, -the best analog. How we define the best analog for a certain day? This function offers three options for that: +For instance if we want to perform a temperature donwscaling in Balearic Island for October we will get a daily series of temperature with 1 analog per day, the best analog. How we define the best analog for a certain day? This function offers three options for that: (1) The day with the minimum Euclidean distance in a large scale field: using i.e. pressure or geopotencial height as variables and North Atlantic region as large scale region. The Atmospheric circulation pattern in the North Atlantic (LargeScale) has an important role in the climate in Spain (LocalScale). The function will find the day with the most similar pattern in atmospheric circulation in the database (obs, slp in ERA5) to the day of interest (exp, slp in model). Once the date of the best analog is found, the function takes the associated temperature to that day (obsVar, tas in ERA5), with a subset of the region of interest (Balearic Island). @@ -39,7 +38,7 @@ Two datasets are used to illustrate how to use the function. The first one could ### Example 1: using data from CSTools -After loading **CSTools** package on the R session, the user will have access to the sample data `lonlat_temp` and `lonlat_prec`. +After loading **CSTools** package on the R session, the user will have access to the sample data created with using `CST_Start`: lonlat_temp_st` and `lonlat_prec_st`. *Note: If it is the first time using CSTools, install the package by running `install.packages("CSTools")`. @@ -50,32 +49,24 @@ library(CSTools) After exploring the data, the user can directly run the Analogs downscaling method using the 'Large_dis' metric: ``` -class(lonlat_temp$exp) -names(lonlat_temp$obs) -dim(lonlat_temp$obs$data) -dim(lonlat_temp$exp$data) -head(lonlat_temp$exp$attrs$Dates) +class(lonlat_temp_st$exp) +names(lonlat_temp_st$obs) +dim(lonlat_temp_st$obs$data) +dim(lonlat_temp_st$exp$data) +lonlat_temp_st$exp$attrs$Dates +lonlat_temp_st$obs$attrs$Dates ``` -There are 15 ensemble members available in the `exp` data set, 6 starting dates and 3 forecast times, which refer to monthly values during 3 months following starting dates on November 1st in the years 2000, 2001, 2002, 2003, 2004 and 2005. +There are 15 ensemble members available in the `exp` data set, 6 starting dates and 3 forecast times, which refer to monthly values during 3 months following starting dates on November in the years 2000, 2001, 2002, 2003, 2004 and 2005. ``` -exp1 <- lonlat_temp$exp -exp1$data <- exp1$data[, , 1, 1, , , drop = FALSE] -exp1$attrs$Dates <- exp1$attrs$Dates[1] +exp1 <- CST_Subset(x = lonlat_temp_st$exp, along = c('sdate', 'ftime'), indices = list(1, 1)) +down_1 <- CST_Analogs(expL = exp1, obsL = lonlat_temp_st$obs) -down_1 <- CST_Analogs(expL = exp1, obsL = lonlat_temp$obs) +exp2 <- CST_Subset(x = lonlat_temp_st$exp, along = c('sdate', 'ftime'), indices = list(1, 2)) +down_2 <- CST_Analogs(expL = exp2, obsL = lonlat_temp_st$obs) -exp2 <- lonlat_temp$exp -exp2$data <- exp2$data[, , 1, 2, , , drop = FALSE] -exp2$attrs$Dates <- exp2$attrs$Dates[2] - -down_2 <- CST_Analogs(expL = exp2, obsL = lonlat_temp$obs) - -exp3 <- lonlat_temp$exp -exp3$data <- exp3$data[, , 1, 3, , , drop = FALSE] -exp3$attrs$Dates <- exp3$attrs$Dates[3] - -down_3 <- CST_Analogs(expL = exp3, obsL = lonlat_temp$obs) +exp3 <- CST_Subset(x = lonlat_temp_st$exp, along = c('sdate', 'ftime'), indices = list(1, 3)) +down_3 <- CST_Analogs(expL = exp3, obsL = lonlat_temp_st$obs) ``` The visualization of the first three time steps for the ensemble mean of the forecast initialized the 1st of Noveber 2000 can be done using the package **s2dv**: @@ -97,13 +88,12 @@ PlotLayout(PlotEquiMap, c('lat', 'lon'), width = 10, height = 4) ``` - ![](./Figures/Analogs1.png) The user can also request extra Analogs and the information: ``` -down <- CST_Analogs(expL = exp1, obsL = lonlat_temp$obs, +down <- CST_Analogs(expL = exp1, obsL = lonlat_temp_st$obs, nAnalogs = 2, AnalogsInfo = TRUE) ``` @@ -121,7 +111,7 @@ The last command run concludes that the best analog of the ensemble 15 correspon ``` PlotLayout(PlotEquiMap, c('lat', 'lon'), var = list(down$data$fields[1, , , 15], - lonlat_temp$obs$data[1, 1, 5, 1, , ]), nrow = 1, ncol = 2, + lonlat_temp_st$obs$data[1, 1, 5, 1, , ]), nrow = 1, ncol = 2, lon = down$coords$lon, lat = down$coords$lat, filled.continents = FALSE, titles = c("Downscaled 2000-11-01", "Observed 2004-11-01"), units = 'T(K)', width = 7, height = 4) @@ -131,71 +121,118 @@ PlotLayout(PlotEquiMap, c('lat', 'lon'), var = list(down$data$fields[1, , , 15], As expected, they are exatly the same. -### Exemple 2: Load data using CST_Load +### Exemple 2: Load data using CST_Start -In this case, the spatial field of a single forecast day will be downscale using Analogs in this example. This will allow illustrating how to use CST_Load to retrieve observations separated from simulations. To explore other options, see other CSTools vignettes as well as `CST_Load` documentation. +In this case, the spatial field of a single forecast day will be downscale using Analogs in this example. This will allow illustrating how to use `CST_Start` to retrieve observations separated from simulations. To explore other options, see other CSTools vignettes as well as `CST_Start` documentation and [startR](https://cran.r-project.org/web/packages/startR/index.html) package. The simulations available for the desired model cover the period 1993-2016. Here, the 15th of October 2000 (for the simulation initialized in the 1st of October 2000), will be downscaled. For ERA5 from 1979 to the present days. For this example we will just use October days from 2000 to 2006, so, the starting dates can be defined by running the following lines: ``` start <- as.Date(paste(2000, 10, "01", sep = ""), "%Y%m%d") end <- as.Date(paste(2006, 10, "01", sep = ""), "%Y%m%d") -dateseq <- format(seq(start, end, by = "year"), "%Y%m%d") +dates <- as.POSIXct(seq(start, end, by = "year"), format = '%Y%m%d', 'UTC') ``` -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. The experimental datasets are interpolated to the ERA5 grid by specifying the 'grid' parameter while ERA5 doesn't need to be interpolated. While parameter leadtimemax is set to 1 for the experimental dataset, it is set to 31 for the observations, returning the daily observations for October for the years requested in 'sdate' (2000-2006). Download the data to run the recipe under the HTTPS: downloads.cmcc.bo.it/d_chaves/ANALOGS/data_for_Analogs.Rdat or ask carmen.alvarez-castro at cmcc.it or nuria.perez at bsc.es. - -``` -exp <- list(name = 'ECMWF_system4_m1', - path = file.path("/esarchive/exp/ecmwf/system4_m1/", - "$STORE_FREQ$_mean/$VAR_NAME$_*/$VAR_NAME$_$START_DATE$.nc")) -obs <- list(name = 'ERA5', - path = file.path("/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/", - "$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc")) - -expTAS <- CST_Load(var = 'tas', exp = list(exp), obs = NULL, - sdates = '20001001', latmin = 22, latmax = 70, - lonmin = -80, lonmax = 50, output = 'lonlat', - storefreq = 'daily', nmember = 15, leadtimemin = 15, - leadtimemax = 15, method = "bilinear", grid = 'r1440x721', - nprocs = 1) -obsTAS <- CST_Load(var = 'tas', exp = NULL, obs = list(obs), - sdates = dateseq, leadtimemax = 31, - latmin = 22, latmax = 70, - lonmin = -80, lonmax = 50, output = 'lonlat', - nprocs = 1, storefreq = "daily", nmember = 1) - -expPSL <- CST_Load(var = 'psl', exp = list(exp), obs = NULL, - sdates = '20001001', latmin = 22, latmax = 70, - lonmin = -80, lonmax = 50, output = 'lonlat', - storefreq = 'daily', nmember = 15, leadtimemin = 15, - leadtimemax = 15, method = "bilinear", grid = 'r1440x721', - nprocs = 1) -obsPSL <- CST_Load(var = 'psl', exp = NULL, obs = list(obs), - sdates = dateseq, leadtimemax = 31, - latmin = 22, latmax = 70, - lonmin = -80, lonmax = 50, output = 'lonlat', - nprocs = 1, storefreq = "daily", nmember = 1) - -``` - -*Note: `CST_Load` allows to load the data simultaneously for 'exp' and 'obs' already formatted to have the same dimensions as in this example. However, it is possible to request separated 'obs' and 'exp'. In this second case, the observations could be return in a continous time series instead of being split in start dates and forecast time.* - -The s2dv_cube objects `expTAS`,`obsTAS`, `expPSL` and `obsPSL` are now loaded in the R enviroment. The first two elements correspond to the experimental and observed data for temperature and the other two are the equivalent for the SLP data. - -Loading the data using `CST_Load` allows to obtain two lists, one for the experimental data and another for the observe data, with the same elements and compatible dimensions of the data element: - +Using the `CST_Start` function from **CSTools package**, the data available in our data store can be loaded. The following lines show how this function can be used. The experimental datasets are interpolated to the ERA5 grid by specifying the 'grid' parameter while ERA5 doesn't need to be interpolated. While dimension 'time' is set to 1 for the experimental dataset, it is set to 31 for the observations, returning the daily observations for October for the years requested in 'sdate' (2000-2006). Download the data to run the recipe under the HTTPS: downloads.cmcc.bo.it/d_chaves/ANALOGS/data_for_Analogs.Rdat or ask carmen.alvarez-castro at cmcc.it or nuria.perez at bsc.es. + +``` +exp_path <- paste0('/esarchive/exp/ecmwf/system4_m1/daily_mean/', + '$var$_f6h/$var$_$sdate$.nc') +obs_path <- paste0('/esarchive/recon/ecmwf/era5/daily_mean/', + '$var$_f1h-r1440x721cds/$var$_$sdate$.nc') + +date_exp <- '20001001' +lonmax <- 50 +lonmin <- -80 +latmax <- 70 +latmin <- 22 + +expTAS <- CST_Start(dataset = exp_path, + var = 'tas', + member = indices(1:15), + sdate = '20001001', + ftime = indices(15), + lat = values(list(latmin, latmax)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lonmin, lonmax)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + member = c('member', 'ensemble'), + ftime = c('ftime', 'time')), + transform = CDORemapper, + transform_params = list(grid = 'r1440x721', + method = 'bilinear'), + transform_vars = c('lat', 'lon'), + return_vars = list(lat = NULL, + lon = NULL, ftime = 'sdate'), + retrieve = TRUE) + +expPSL <- CST_Start(dataset = exp_path, + var = 'psl', + member = indices(1:15), + sdate = '20001001', + ftime = indices(15), + lat = values(list(latmin, latmax)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lonmin, lonmax)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + member = c('member', 'ensemble'), + ftime = c('ftime', 'time')), + transform = CDORemapper, + transform_params = list(grid = 'r1440x721', + method = 'bilinear'), + transform_vars = c('lat', 'lon'), + return_vars = list(lat = NULL, + lon = NULL, ftime = 'sdate'), + retrieve = TRUE) + +obsTAS <- CST_Start(dataset = obs_path, + var = 'tas', + sdate = unique(format(dates, '%Y%m')), + ftime = indices(1:31), + lat = values(list(latmin, latmax)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lonmin, lonmax)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + ftime = c('ftime', 'time')), + return_vars = list(lat = NULL, + lon = NULL, ftime = 'sdate'), + retrieve = TRUE) + +obsPSL <- CST_Start(dataset = obs_path, + var = 'psl', + sdate = unique(format(dates, '%Y%m')), + ftime = indices(1:31), + lat = values(list(latmin, latmax)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lonmin, lonmax)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + ftime = c('ftime', 'time')), + return_vars = list(lat = NULL, + lon = NULL, ftime = 'sdate'), + retrieve = TRUE) +``` + +The 's2dv_cube' objects `expTAS`,`obsTAS`, `expPSL` and `obsPSL` are now loaded in the R environment. The first two elements correspond to the experimental and observed data for temperature and the other two are the equivalent for the SLP data. + +Loading the data using `CST_Start` allows to obtain two lists, one for the experimental data and another for the observe data, with the same elements and compatible dimensions of the data element: ``` dim(expTAS$data) -# dataset member sdate ftime lat lon -# 1 15 1 1 193 521 +# dataset var member sdate ftime lat lon +# 1 1 15 1 1 193 521 dim(obsTAS$data) -# dataset member sdate ftime lat lon -# 1 1 7 31 193 521 +# dataset var sdate ftime lat lon +# 1 1 7 31 193 521 ``` - #### Two variables and criteria Large [scale] Distance: The aim is to downscale the temperature field of the simulation for the 15th of October 2000 but looking at the pressure pattern: @@ -210,11 +247,11 @@ Some warnings could appear indicating information about undefining parameters. I ``` names(down1$data) dim(down1$data$field) -# nAnalogs lat lon member time -# 3 193 521 15 1 +# nAnalogs lat lon member +# 3 193 521 15 dim(down1$data$dates) -# nAnalogs member time -# 3 15 1 +# nAnalogs member +# 3 15 down1$data$dates[1,1] # "07-10-2005" ``` @@ -222,14 +259,14 @@ Now, we can visualize the output: ``` PlotLayout(PlotEquiMap, c('lat', 'lon'), - var = list(expPSL$data[1, 1, 1, 1, , ], obsPSL$data[1, 1, 1, 15, , ], + var = list(expPSL$data[1, 1, 1, 1, 1, , ], obsPSL$data[1, 1, 1, 15, , ], obsPSL$data[1, 1, 6, 7, , ]), lon = obsPSL$coords$lon, lat = obsPSL$coords$lat, filled.continents = FALSE, titles = c('Exp PSL 15-10-2000','Obs PSL 15-10-2000', 'Obs PSL 7-10-2005'), toptitle = 'First member', ncol = 3, nrow = 1, width = 10, height = 4) PlotLayout(PlotEquiMap, c('lat', 'lon'), var = list( - expTAS$data[1, 1, 1, 1, , ], obsTAS$data[1, 1, 1, 15, , ], + expTAS$data[1, 1, 1, 1, 1, , ], obsTAS$data[1, 1, 1, 15, , ], down1$data$field[1, , , 1], obsTAS$data[1, 1, 6, 7, , ]), lon = obsTAS$coords$lon, lat = obsTAS$coords$lat, filled.continents = FALSE, titles = c('Exp TAS 15-10-2000', 'Obs TAS 15-10-2000', @@ -252,8 +289,8 @@ The aim is to downscale the temperature simulation of the 15th of October 2000, ``` region <- c(lonmin = 0, lonmax = 5, latmin = 38.5, latmax = 40.5) -expPSL$data <- expPSL$data[1, 1, 1, 1, , ] -expTAS$data <- expTAS$data[1, 1, 1, 1, , ] +expPSL$data <- expPSL$data[1, 1, 1, 1, 1, , ] +expTAS$data <- expTAS$data[1, 1, 1, 1, 1, , ] down2 <- CST_Analogs(expL = expPSL, obsL = obsPSL, AnalogsInfo = TRUE, criteria = "Local_dist", # nAnalogs = 50, obsVar = obsTAS, expVar = expTAS, @@ -262,11 +299,11 @@ down2 <- CST_Analogs(expL = expPSL, obsL = obsPSL, AnalogsInfo = TRUE, The parameter 'nAnalogs' doesn't correspond to the number of Analogs returned, but to the number of the best observations to use in the comparison between large and local scale. -In this case, when looking to a large scale pattern and also to local scale pattern the best analog for the first member is the 13th of October 2001: +In this case, when looking to a large scale pattern and also to local scale pattern the best analog for the first member is the 7th of October 2001: ``` down2$data$dates[2] -# [1] "13-10-2001" +# "13-10-2001" ``` ``` @@ -276,7 +313,7 @@ var = list(expTAS$data, obsTAS$data[1, 1, 1, 15, , ], down2$data$field[1, , ], SelBox(obsTAS$data[1, 1, 2, 13, , ], lon = as.vector(obsTAS$coords$lon), lat = as.vector(obsTAS$coords$lat), - region)$data) + region, londim = 'lon', latdim = 'lat')$data) PlotLayout(PlotEquiMap, c('lat', 'lon'), var = var, special_args = list(list(lon = expTAS$coords$lon, lat = expTAS$coords$lat), @@ -311,9 +348,11 @@ down3$data$dates[3] ``` ``` -var = list(down3$data$field[1, , ], SelBox(obsTAS$data[1, 1, 2, 10, , ], - lon = as.vector(obsTAS$coords$lon), lat = as.vector(obsTAS$coords$lat), - region)$data) +var = list(down3$data$field[1, , ], + SelBox(obsTAS$data[1, 1, 2, 10, , ], + lon = as.vector(obsTAS$coords$lon), + lat = as.vector(obsTAS$coords$lat), + region, londim = 'lon', latdim = 'lat')$data) PlotLayout(PlotEquiMap, c('lat', 'lon'), var = var, lon = down3$coords$lon, lat = down3$coords$lat, diff --git a/vignettes/Data_Considerations.Rmd b/vignettes/Data_Considerations.Rmd index 979e1c751727edfab63b9c767ca098170986f926..e6e207241a73b4a139969df8cb7db0a0ea85af26 100644 --- a/vignettes/Data_Considerations.Rmd +++ b/vignettes/Data_Considerations.Rmd @@ -1,6 +1,8 @@ --- author: "Nuria Perez" date: "`r Sys.Date()`" +revisor: "Eva Rifà" +revision date: "October 2023" output: rmarkdown::html_vignette vignette: > %\VignetteEngine{knitr::knitr} @@ -49,6 +51,7 @@ All CSTools functions have been developed following the same guidelines. The mai A reasonable important doubt that a new user may have at this point is: what 's2dv_cube' object is? 's2dv_cube' is a class of an object storing the data and metadata in several elements: + $data element is an N-dimensional array with named dimensions containing the data (e.g.: temperature values), + + $dims vector with the dimensions of $data + $coords is a named list with elements of the coordinates vectors corresponding to the dimensions of the $data, + $attrs is a named list with elements corresponding to attributes of the object. It has the following elements: + $Variable is a list with the variable name in element $varName and with the metadata of all the variables in the $metadata element, @@ -59,9 +62,9 @@ It is possible to visualize an example of the structure of 's2dv_cube' object by ``` library(CSTools) -class(lonlat_temp$exp) # check the class of the object lonlat_temp$exp -names(lonlat_temp$exp) # shows the names of the elements in the object lonlat_temp$exp -str(lonlat_temp$exp) # shows the full structure of the object lonlat_temp$exp +class(lonlat_temp_st$exp) # check the class of the object lonlat_temp$exp +names(lonlat_temp_st$exp) # shows the names of the elements in the object lonlat_temp$exp +str(lonlat_temp_st$exp) # shows the full structure of the object lonlat_temp$exp ``` ### 3. Data storage recommendations @@ -75,9 +78,11 @@ CSTools main objective is to share state-of-the-arts post-processing methods wit - CST_Load can perform spatial averages over a defined region or return the lat-lon grid and - CST_Load can read from files using multiple parallel processes among other possibilites. -If you plan to use CST_Load, we have developed guidelines to download and formatting the data. See [CDS_Seasonal_Downloader](https://earth.bsc.es/gitlab/es/cds-seasonal-downloader). +CSTools also has the function `CST_Start` from [startR](https://cran.r-project.org/web/packages/startR/index.html) that is more flexible than `CST_Load`. We recommend to use `CST_Start` since it's more efficient and flexible. -There are alternatives to CST_Load function, for instance, the user can: +If you plan to use `CST_Load` or `CST_Start`, we have developed guidelines to download and formatting the data. See [CDS_Seasonal_Downloader](https://earth.bsc.es/gitlab/es/cds-seasonal-downloader). + +There are alternatives to these functions, for instance, the user can: 1) use another tool to read the data from files (e.g.: ncdf4, easyNDCF, startR packages) and then convert it to the class 's2dv_cube' with `s2dv.cube()` function or 2) If they keep facing problems to convert the data to that class, they can just skip it and work with the functions without the prefix 'CST_'. In this case, they will be able to work with the basic class 'array'. @@ -108,6 +113,55 @@ c(exp, obs) %<-% CST_Load(var = 'sfcWind', grid = "r360x180") ``` +### 5. CST_Start example +``` +path_exp <- paste0('/esarchive/exp/meteofrance/system6c3s/monthly_mean/', + '$var$_f6h/$var$_$sdate$.nc') + +sdates <- sapply(1993:2012, function(x) paste0(x, '0501')) + +lonmax <- 60.5 +lonmin <- -19 +latmax <- 79.5 +latmin <- 0 +exp <- CST_Start(dataset = path_exp, + var = 'sfcWind', + ensemble = indices(1:9), + sdate = sdates, + time = indices(1:3), + latitude = values(list(latmin, latmax)), + latitude_reorder = Sort(decreasing = TRUE), + longitude = values(list(lonmin, lonmax)), + longitude_reorder = CircularSort(0, 360), + synonims = list(longitude = c('lon', 'longitude'), + latitude = c('lat', 'latitude')), + return_vars = list(latitude = NULL, + longitude = NULL, time = 'sdate'), + retrieve = TRUE) + +path_obs <- paste0('/esarchive/recon/ecmwf/erainterim/daily_mean/', + '$var$_f6h/$var$_$sdate$.nc') +dates <- as.POSIXct(sdates, format = '%Y%m%d', 'UTC') +obs <- CST_Start(dataset = path_obs, + var = 'sfcWind', + sdate = unique(format(dates, '%Y%m')), + time = indices(2:4), + latitude = values(list(latmin, latmax)), + latitude_reorder = Sort(decreasing = TRUE), + longitude = values(list(lonmin, lonmax)), + longitude_reorder = CircularSort(0, 360), + synonims = list(longitude = c('lon', 'longitude'), + latitude = c('lat', 'latitude')), + transform = CDORemapper, + transform_extra_cells = 2, + transform_params = list(grid = 'r360x181', + method = 'conservative'), + transform_vars = c('latitude', 'longitude'), + return_vars = list(longitude = NULL, + latitude = NULL, + time = 'sdate'), + retrieve = TRUE) +``` Extra lines to see the size of the objects and visualize the data: @@ -118,11 +172,10 @@ object_size(exp) object_size(obs) # 3.09 MB library(s2dv) -PlotEquiMap(exp$data[1,1,1,1,,], lon = exp$coords$lon, lat= exp$coords$lat, +PlotEquiMap(exp$data[1,1,1,1,1,,], lon = exp$coords$longitude, lat= exp$coords$latitude, filled.continents = FALSE, fileout = "Meteofrance_r360x180.png") ``` - ![Meteofrance](../vignettes/Figures/Meteofrance_r360x180.png) ### Managing big datasets and memory issues diff --git a/vignettes/ENSclustering_vignette.Rmd b/vignettes/ENSclustering_vignette.Rmd index 073a32100d71cf19991b1081eaa7a88ac5ce7d04..4eb96e7fad24d7c5de5991b659228db6ef08a8b9 100644 --- a/vignettes/ENSclustering_vignette.Rmd +++ b/vignettes/ENSclustering_vignette.Rmd @@ -1,6 +1,8 @@ --- author: "Ignazio Giuntoli and Federico Fabiano - CNR-ISAC" date: "`r Sys.Date()`" +revisor: "Eva Rifà" +revision date: "October 2023" output: rmarkdown::html_vignette vignette: > %\VignetteEngine{knitr::knitr} @@ -42,15 +44,15 @@ For our example we will use the sample seasonal temperature data provided within Data can be loaded as follows: ```r -datalist <- lonlat_temp$exp +datalist <- lonlat_temp_st$exp ``` The data will has the following dimension: ```r dim(datalist$data) -dataset member sdate ftime lat lon - 1 15 6 3 22 53 +dataset var member sdate ftime lat lon + 1 1 15 6 3 22 53 ``` Therefore the number of members is 15, the number of start dates is 6, while the forecast time steps are 3. The lat and lon dimensions refer to a 22x53 grid. @@ -101,11 +103,11 @@ results$freq [4,] 15.55556 ``` -Further, the cluster number to which each 'member - start-date' pair is assigned can be displayed by quering '$cluster' in 'results' as shown below (members (15) are in row and the start-dates (6) are in column (i.e. 15*6 pairs). +Further, the cluster number to which each 'member - start-date' pair is assigned can be displayed by quering '$cluster' in 'results' as shown below (members (15) are in row and the start-dates (6) are in column (i.e. 15*6 pairs)). ```r results$cluster - +, , 1, 1 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3 2 1 1 4 2 [2,] 2 4 4 2 3 1 @@ -131,16 +133,16 @@ To achieve this we can get an idea of how the clustering has performed by lookin ```r dim(results$composites) -cluster lat lon dataset - 4 22 53 1 +cluster lat lon dataset var + 4 22 53 1 1 ``` while the 'repr_field' argument of 'results' provides the spatial pattern of the member lying closes to the centroid (has the same dimensions as those of 'composites'): ```r dim(results$repr_field) -cluster lat lon dataset - 4 22 53 1 +cluster lat lon dataset var + 4 22 53 1 1 ``` Finally, the pairs 'member - start-dates' to be picked as representative for each cluster are found in the 'closest_member' argument. @@ -169,12 +171,12 @@ The following lines produce a multiplot of the representative temperature anomal These are actually the closest realizations (member - start-date pair) to the cluster centroids noted as 'repr_field' above. ```r -EnsMean <- MeanDims(datalist $data, c('member', 'sdate', 'ftime')) -EnsMean <- InsertDim(Reorder(EnsMean, c("lat", "lon", "dataset")), +EnsMean <- MeanDims(datalist$data, c('member', 'sdate', 'ftime')) +EnsMean <- InsertDim(Reorder(EnsMean, c("lat", "lon", "dataset", "var")), posdim = 1, lendim = 4, name = 'cluster') PlotLayout(PlotEquiMap, plot_dims = c("lat", "lon"), - var = results$repr_field[,,,1] - EnsMean[,,,1], + var = results$repr_field[,,,1,1] - EnsMean[,,,1,1], lon = results$lon, lat = results$lat, filled.continents = FALSE, titles = c("1","2","3","4"), brks = seq(-2, 2, 0.5), fileout = "EnsClus_4clus_both_mem_std_Fig1.png") @@ -188,16 +190,18 @@ The lines below produce a multiplot of the temperature anomaly patterns for each ```r ExpMean <- MeanDims(datalist$data, 'ftime') -EnsMean <- InsertDim(InsertDim(InsertDim(EnsMean[1,,,], +EnsMean <- InsertDim(InsertDim(InsertDim(InsertDim(EnsMean[1,,,,], posdim = 1, lendim = 6, name = 'sdate'), posdim = 1, lendim = 15, name = 'member'), - posdim = 1, lendim = 1, name = 'dataset') -ExpMeanSd <- Reorder(ExpMean - EnsMean, c('dataset', 'sdate', 'member' , 'lat', 'lon')) + posdim = 1, lendim = 1, name = 'var'), + posdim = 1, lendim = 1, name = 'dataset') + +ExpMeanSd <- Reorder(ExpMean - EnsMean, c('dataset', 'var', 'sdate', 'member' , 'lat', 'lon')) PlotLayout(PlotEquiMap, plot_dims = c("lat", "lon"), - var = ExpMeanSd[, , 1:4, , ], title_scale = 0.7, + var = ExpMeanSd[, , , 1:4, , ], title_scale = 0.7, ncol = 6, nrow = 4, row_titles = paste('member' , 1:4), col_titles = paste('sdate', 1:6), lon = results$lon, lat = results$lat, filled.continents = FALSE, - titles = as.character(t(results$cluster[1:4, 1:6,])), brks = seq(-2, 2, 0.5), + titles = as.character(t(results$cluster[1:4, 1:6,,])), brks = seq(-2, 2, 0.5), width = 24, height = 20, size_units = 'cm', fileout = "EnsClus_4clus_both_mem_std_Fig2.png") ``` diff --git a/vignettes/Figures/Analogs1.png b/vignettes/Figures/Analogs1.png index 5b4f05a94584422a6e58a638888be481a197b3b0..fa60c4e61ae52e46c6c3850bb68baa13376d99f6 100644 Binary files a/vignettes/Figures/Analogs1.png and b/vignettes/Figures/Analogs1.png differ diff --git a/vignettes/Figures/Analogs2.png b/vignettes/Figures/Analogs2.png index eb67ce3f2b5a3b405a1c36cfbf5ec466d170e73d..3b37a4c020ab6dc1e785701502e83148dcddd621 100644 Binary files a/vignettes/Figures/Analogs2.png and b/vignettes/Figures/Analogs2.png differ diff --git a/vignettes/Figures/Analogs3.png b/vignettes/Figures/Analogs3.png index 465510727fdd046a35c8318400df73f18c57d2a0..08097fdea8f0598b36d9b6db381cbda7994b8917 100644 Binary files a/vignettes/Figures/Analogs3.png and b/vignettes/Figures/Analogs3.png differ diff --git a/vignettes/Figures/Analogs5.png b/vignettes/Figures/Analogs5.png index 79907c8065e536f3a7e70536ee8b79efc9eed71d..c6343cae34495baed514525f15c068d5aa433815 100644 Binary files a/vignettes/Figures/Analogs5.png and b/vignettes/Figures/Analogs5.png differ diff --git a/vignettes/Figures/Analogs6.png b/vignettes/Figures/Analogs6.png index abbc0d7a93988ffd4f2350ba746474188896214e..e509245b8e1c3fb454dc8e541ea96b7f058cb27d 100644 Binary files a/vignettes/Figures/Analogs6.png and b/vignettes/Figures/Analogs6.png differ diff --git a/vignettes/Figures/EnsClus_4clus_both_mem_std_Fig1.png b/vignettes/Figures/EnsClus_4clus_both_mem_std_Fig1.png index 123fdb7795fc91926af753592d66c21d60f2b092..9e3a64ba6484eefcd6d82beceff60f314b37b353 100644 Binary files a/vignettes/Figures/EnsClus_4clus_both_mem_std_Fig1.png and b/vignettes/Figures/EnsClus_4clus_both_mem_std_Fig1.png differ diff --git a/vignettes/Figures/EnsClus_4clus_both_mem_std_Fig2.png b/vignettes/Figures/EnsClus_4clus_both_mem_std_Fig2.png index 2c05166cfa25e349166f301957fa23a47e96efcd..fecc8d94feb81cf0e6baa877200e9e3a98851669 100644 Binary files a/vignettes/Figures/EnsClus_4clus_both_mem_std_Fig2.png and b/vignettes/Figures/EnsClus_4clus_both_mem_std_Fig2.png differ diff --git a/vignettes/Figures/Meteofrance_r360x180.png b/vignettes/Figures/Meteofrance_r360x180.png index 1438bd493288c10cf80b1796660e97ba51e5df56..0cfc53e224e8737fc9c1afd3ceec8f29b9344bf0 100644 Binary files a/vignettes/Figures/Meteofrance_r360x180.png and b/vignettes/Figures/Meteofrance_r360x180.png differ diff --git a/vignettes/Figures/MostLikelyTercile_fig1.png b/vignettes/Figures/MostLikelyTercile_fig1.png index bd282ed25fda569eaca3c29616e0e5c9d471eda5..5f5b4322607ec7d6a098cc1f9ca40e5a5dce9f57 100644 Binary files a/vignettes/Figures/MostLikelyTercile_fig1.png and b/vignettes/Figures/MostLikelyTercile_fig1.png differ diff --git a/vignettes/Figures/MostLikelyTercile_fig2.png b/vignettes/Figures/MostLikelyTercile_fig2.png index b96854ad8a9a3ded41e0d3be2a148eda4b2c47b5..62f6a5f05a559fb83b4dce1e101789c2c35a15cb 100644 Binary files a/vignettes/Figures/MostLikelyTercile_fig2.png and b/vignettes/Figures/MostLikelyTercile_fig2.png differ diff --git a/vignettes/Figures/MostLikelyTercile_fig3.png b/vignettes/Figures/MostLikelyTercile_fig3.png index 1fa8460723a07d2abb08251758a585a2c3bcf7fc..72bb26efdc85f9fa87e98f8e0ffc244acf4f922c 100644 Binary files a/vignettes/Figures/MostLikelyTercile_fig3.png and b/vignettes/Figures/MostLikelyTercile_fig3.png differ diff --git a/vignettes/Figures/MultiModelSkill_cor_tas_1993-2012.png b/vignettes/Figures/MultiModelSkill_cor_tas_1993-2012.png index d87618198c05bcc4ec3fa3e9cf26dd9c5432bb55..962eedb6b65c280cc4dc7f49aca53ecdf254d39f 100644 Binary files a/vignettes/Figures/MultiModelSkill_cor_tas_1993-2012.png and b/vignettes/Figures/MultiModelSkill_cor_tas_1993-2012.png differ diff --git a/vignettes/Figures/MultiModelSkill_rms_tas_1993-2012.png b/vignettes/Figures/MultiModelSkill_rms_tas_1993-2012.png index 199ab861e264bcb25a48959b8dfcbefb179dc473..c0ec379c01cc9cbe767a4cf46558e526723d25ee 100644 Binary files a/vignettes/Figures/MultiModelSkill_rms_tas_1993-2012.png and b/vignettes/Figures/MultiModelSkill_rms_tas_1993-2012.png differ diff --git a/vignettes/Figures/MultiModelSkill_rmsss_tas_1993-2012.png b/vignettes/Figures/MultiModelSkill_rmsss_tas_1993-2012.png index c8ab0c65bee47265e5b1c0990ff72ce4a5f56cf2..87f8c888489ddba9289596e14a8f3897b3b3be2b 100644 Binary files a/vignettes/Figures/MultiModelSkill_rmsss_tas_1993-2012.png and b/vignettes/Figures/MultiModelSkill_rmsss_tas_1993-2012.png differ diff --git a/vignettes/Figures/MultivarRMSE_gloseas5_tas_prlr_1993-2012.png b/vignettes/Figures/MultivarRMSE_gloseas5_tas_prlr_1993-2012.png index 14d7055e00914ddf144262d07dc8af2b6d1c8c3e..31f4a2723982ebd2f71f2d3a2e2e6aa9e63143b3 100644 Binary files a/vignettes/Figures/MultivarRMSE_gloseas5_tas_prlr_1993-2012.png and b/vignettes/Figures/MultivarRMSE_gloseas5_tas_prlr_1993-2012.png differ diff --git a/vignettes/Figures/PlotForecastPDF_ex1.png b/vignettes/Figures/PlotForecastPDF_ex1.png index 516bd717e792da8c885d42257ae4d9fcb72291b1..6fd220afde64fbe9c8b7a7e5847eddfcf51f08ae 100644 Binary files a/vignettes/Figures/PlotForecastPDF_ex1.png and b/vignettes/Figures/PlotForecastPDF_ex1.png differ diff --git a/vignettes/Figures/PlotForecastPDF_ex2.png b/vignettes/Figures/PlotForecastPDF_ex2.png index c7f9d9e434180700fa858117c8cc52bf5835957c..f427fb1d090c5f1743e87659e2d1e4b6aa8e8a05 100644 Binary files a/vignettes/Figures/PlotForecastPDF_ex2.png and b/vignettes/Figures/PlotForecastPDF_ex2.png differ diff --git a/vignettes/Figures/PlotForecastPDF_ex3.png b/vignettes/Figures/PlotForecastPDF_ex3.png index d8eec5bedfc2696266972d2eb61a8eae5e514764..5fe55ad04d31e621d39c254ba7c97ea91766a9a7 100644 Binary files a/vignettes/Figures/PlotForecastPDF_ex3.png and b/vignettes/Figures/PlotForecastPDF_ex3.png differ diff --git a/vignettes/Figures/PlotForecastPDF_ex4.png b/vignettes/Figures/PlotForecastPDF_ex4.png index 7254ee2e5a21740739659f5ab5b87b4927f5383a..327375794de21bb3badd7af6233fde6bf37a6072 100644 Binary files a/vignettes/Figures/PlotForecastPDF_ex4.png and b/vignettes/Figures/PlotForecastPDF_ex4.png differ diff --git a/vignettes/Figures/RainFARM_fig1.png b/vignettes/Figures/RainFARM_fig1.png index 9c80d8fb80e3e6fc05af907ccf811193e6e7ae59..1eaede22326e9a35ce8440d86ce1222fc36ed29e 100644 Binary files a/vignettes/Figures/RainFARM_fig1.png and b/vignettes/Figures/RainFARM_fig1.png differ diff --git a/vignettes/Figures/observed_regimes.png b/vignettes/Figures/observed_regimes.png index 678ac72aad7fd3f7e55e01abcbd90372c838b3c0..e323757d43337a036cf0c1dd60cbe843cfb4159b 100644 Binary files a/vignettes/Figures/observed_regimes.png and b/vignettes/Figures/observed_regimes.png differ diff --git a/vignettes/Figures/predicted_regimes.png b/vignettes/Figures/predicted_regimes.png index 9f69484f5c97a041454330247f967b5c912c2cff..a2e1c3e062614fff59f5ce5043f4dd8e2feb445d 100644 Binary files a/vignettes/Figures/predicted_regimes.png and b/vignettes/Figures/predicted_regimes.png differ diff --git a/vignettes/MostLikelyTercile_vignette.Rmd b/vignettes/MostLikelyTercile_vignette.Rmd index aa9e998e1465afb5d9f9cdafb00c28e09d3de370..ded05182ac9bbe7bd865a0232f034d2b1aa03c0e 100644 --- a/vignettes/MostLikelyTercile_vignette.Rmd +++ b/vignettes/MostLikelyTercile_vignette.Rmd @@ -1,6 +1,8 @@ --- 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 = indices(1:10), + sdate = dateseq, + ftime = indices(2:4), + lat = values(list(lat_min, lat_max)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lon_min, lon_max)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + member = c('member', 'ensemble'), + ftime = c('ftime', 'time')), + transform = 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 = values(list(lat_min, lat_max)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lon_min, lon_max)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + ftime = c('ftime', 'time')), + transform = 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) diff --git a/vignettes/MultiModelSkill_vignette.Rmd b/vignettes/MultiModelSkill_vignette.Rmd index e927d4328ebf2d3668f09eac417d6d4bcc6afcdb..dd35e256512dfcb3dc084953d4a476d1f81f132d 100644 --- a/vignettes/MultiModelSkill_vignette.Rmd +++ b/vignettes/MultiModelSkill_vignette.Rmd @@ -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,68 @@ 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 = indices(1:4), sdate = dateseq, ftime = indices(2:4), + lat = values(list(latmin, latmax)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lonmin, lonmax)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), lat = c('lat', 'latitude'), + member = c('member', 'ensemble'), ftime = c('ftime', 'time')), + transform = 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 = values(dates_exp), + ftime_across = 'date', + ftime_var = 'ftime', + merge_across_dims = TRUE, + split_multiselected_dims = TRUE, + lat = values(list(latmin, latmax)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lonmin, lonmax)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + ftime = c('ftime', 'time')), + transform = 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 +122,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 +131,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 +154,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 +175,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 +190,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 +220,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 +244,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', diff --git a/vignettes/MultivarRMSE_vignette.Rmd b/vignettes/MultivarRMSE_vignette.Rmd index 744350a8a083a9b0fdbe1a4e75cfec78c45951d6..73f08fe33695c5cc67a9f5a4e11b0199df0a1501 100644 --- a/vignettes/MultivarRMSE_vignette.Rmd +++ b/vignettes/MultivarRMSE_vignette.Rmd @@ -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 = indices(1:9), + sdate = dateseq, + ftime = indices(2:4), + lat = values(list(latmin, latmax)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lonmin, lonmax)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + member = c('member', 'ensemble'), + ftime = c('ftime', 'time')), + transform = 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 = values(dates_exp), + ftime_across = 'date', + ftime_var = 'ftime', + merge_across_dims = TRUE, + split_multiselected_dims = TRUE, + lat = values(list(latmin, latmax)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lonmin, lonmax)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + ftime = c('ftime', 'time')), + transform = 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 = indices(1:9), + sdate = dateseq, + ftime = indices(2:4), + lat = values(list(latmin, latmax)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lonmin, lonmax)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + member = c('member', 'ensemble'), + ftime = c('ftime', 'time')), + transform = 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 = values(dates_exp), + ftime_across = 'date', + ftime_var = 'ftime', + merge_across_dims = TRUE, + split_multiselected_dims = TRUE, + lat = values(list(latmin, latmax)), + lat_reorder = Sort(decreasing = TRUE), + lon = values(list(lonmin, lonmax)), + lon_reorder = CircularSort(0, 360), + synonims = list(lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), + ftime = c('ftime', 'time')), + transform = 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 diff --git a/vignettes/PlotForecastPDF.Rmd b/vignettes/PlotForecastPDF.Rmd index bafbe7d6d001f804329b9219eb5aa823e71e1d9b..34bcb223847d7b767427df28cb01c7f3c1f77e96 100644 --- a/vignettes/PlotForecastPDF.Rmd +++ b/vignettes/PlotForecastPDF.Rmd @@ -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", diff --git a/vignettes/RainFARM_vignette.Rmd b/vignettes/RainFARM_vignette.Rmd index 28ab753ded112a1d670fb8a00abc6dc9d0f8e0b8..a51d75cb9160e14b5b1247ee1a4039486a2f1279 100644 --- a/vignettes/RainFARM_vignette.Rmd +++ b/vignettes/RainFARM_vignette.Rmd @@ -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