Newer
Older
Example
=======
Next you can see an example of usage of s2dverification spanning its
four modules ([Data retrieval](data_retrieval.md),
[Basic statistics](basic_statistics.md), [Verification](verification.md)
and [Visualisation](visualisation.md)).
The goal of the example is to load and analyse the skill of two sample
forecasts of the near-surface air temperature over the Atlanticregion and
assess their skill against a reference reanalysis.
Loading the data
----------------
First the package is loaded and attached.
Then a list is built with information on the location of each dataset to load.
```r
library(s2dverification)
expA <- list(name = 'experimentA',
path = file.path('/path/to/experiments/$EXP_NAME$/monthly_mean',
'$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc'))
expB <- list(name = 'experimentB',
path = file.path('/path/to/experiments/$EXP_NAME$/monthly_mean',
'$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc'))
obsX <- list(name = 'observationX',
path = file.path('/path/to/observations/$OBS_NAME$/monthly_mean',
'$VAR_NAME$/$VAR_NAME$_$YEAR$$MONTH$.nc'))
```
Then the data is loaded with `Load()` providing the previously built lists
-to specify the desired datasets- and other parameters to select the Earth
surface region of interest, starting dates and forecast time steps to load
data from.
In this example, the requested format is 2-dimensional: all the loaded data
will be remapped onto the specified common 'grid' via CDO libraries.
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
```r
sdates <- paste0(1991:1993, '1101')
data <- Load('tas', list(expA, expB), list(obsX),
sdates = sdates,
leadtimemin = 2, leadtimemax = 7,
latmin = -10, latmax = 60,
lonmin = 100, lonmax = 250,
output = 'lonlat', grid = 't106grid',
method = 'distance-weighted')
## * The load call you issued is:
## * Load(var = "tas", exp = list(structure(list(name = "experimentA", path =
## * "/path/to/experiments/$EXP_NAME$/monthly_mean/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc"),
## * .Names = c("name", "path")), structure(list(name =
## * "experimentB", path =
## * "/path/to/experiments/$EXP_NAME$/monthly_mean/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc"),
## * .Names = c("name", "path"))), obs = list(structure(list(name =
## * "observationX", path =
## * "/path/to/observations/$OBS_NAME$/monthly_mean/$VAR_NAME$/$VAR_NAME$_$YEAR$$MONTH$.nc"),
## * .Names = c("name", "path"))), sdates = c("19911101",
## * "19921101", "19931101"), grid = "t106grid", output = "lonlat",
## * storefreq = "monthly", ...)
## * See the full call in '$load_parameters' after Load() finishes.
## * Fetching first experimental files to work out 'var_exp' size...
## * Exploring dimensions... /path/to/experiments/experimentA/monthly_mean/tas/tas_19911101.nc
## * Success. Detected dimensions of experimental data: 2, 5, 3, 6, 63, 134
## * Fetching first observational files to work out 'var_obs' size...
## * Exploring dimensions... /path/to/observations/observationX/monthly_mean/tas/tas_199112.nc
## * Success. Detected dimensions of observational data: 1, 1, 3, 6, 63, 134
## * Will now proceed to read and process 24 data files:
## * /path/to/experiments/experimentA/monthly_mean/tas/tas_19911101.nc
## * /path/to/experiments/experimentA/monthly_mean/tas/tas_19921101.nc
## * /path/to/experiments/experimentA/monthly_mean/tas/tas_19931101.nc
## * /path/to/experiments/experimentB/monthly_mean/tas/tas_19911101.nc
## * /path/to/experiments/experimentB/monthly_mean/tas/tas_19921101.nc
## * /path/to/experiments/experimentB/monthly_mean/tas/tas_19931101.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199112.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199201.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199202.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199203.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199204.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199205.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199212.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199301.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199302.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199303.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199304.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199305.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199312.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199401.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199402.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199403.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199404.nc
## * /path/to/observations/observationX/monthly_mean/tas/tas_199405.nc
## * Total size of requested data: 13372128 bytes.
## * - Experimental data: ( 2 x 5 x 3 x 6 x 63 x 134 ) x 8 bytes = 12156480 bytes.
## * - Observational data: ( 1 x 1 x 3 x 6 x 63 x 134 ) x 8 bytes = 1215648 bytes.
## * If size of requested data is close to or above the free shared RAM memory, R will crash.
## starting worker pid=14895 on localhost:11924 at 15:11:37.477
## starting worker pid=14916 on localhost:11924 at 15:11:37.716
## starting worker pid=14932 on localhost:11924 at 15:11:37.891
## starting worker pid=14948 on localhost:11924 at 15:11:38.110
## starting worker pid=14964 on localhost:11924 at 15:11:38.360
## starting worker pid=14980 on localhost:11924 at 15:11:38.540
## starting worker pid=14996 on localhost:11924 at 15:11:38.718
## starting worker pid=15012 on localhost:11924 at 15:11:38.908
## * Loading... This may take several minutes...
## * Progress: 0% + 33.33% + 33.33% + 33.33%
```
The output consists of two arrays of data (experimental and observational
data) with labelled dimensions, a list of loaded files, a list of not found
files and a call stamp to exactly reproduce as needed, among others.
See [**Data
retrieval**](data_retrieval.md)
for a full explanation of the capabilities and outputs of `Load()`.
Keep in mind that the two arrays obtained from `Load()` will have the
following dimensions in the following order:
c(n. of datasets, n. of members, n. of starting dates, n. of lead-times, n. of latitudes, n. of longitudes)
Checking and summarizing data
-----------------------------
Next, the data of the first member, starting date and leadtime of each dataset
can be visualised on an equidistant projection with `PlotEquiMap()` to check
that the loaded data is as expected.
```r
PlotEquiMap(data$mod[1, 1, 1, 1, , ], data$lon, data$lat)
PlotEquiMap(data$mod[2, 1, 1, 1, , ], data$lon, data$lat)
PlotEquiMap(data$obs[1, 1, 1, 1, , ], data$lon, data$lat)
```
<img src="equi_map_raw_all.png" />
See the full code used to obtain this figure in
[Snippet 1](snippets.md#snippet1).
The values for some starting dates at a single grid point of an experimental
dataset can be plotted together with the observations with the function
`PlotAno()`:
```r
mod <- data$mod[, , , , 30, 60, drop = FALSE]
dim(mod) <- dim(mod)[1:4]
obs <- data$obs[, , , , 30, 60, drop = FALSE]
dim(obs) <- dim(obs)[1:4]
PlotAno(mod, obs, gsub('1101', '1201', sdates),
toptitle = paste0("Experiment ", c("A", "B"),
": raw 'tas' at 166.5ºE, 27.47ºN."),
ytitle = "K",
fileout = c('raw_expA_obsX.eps', 'raw_expB_obsX.eps'))
<img src="raw_expA_obsX.png" width="600" />
<img src="raw_expB_obsX.png" width="600" />
Each coloured region represents data corresponding to a single starting date.
The bold line represents the mean value and the thin lines represent the values
of each ensemble member. The upper and lower limits of each regions are
defined, respectively, by the maximum and minimum values among the members.
The black line represents the reference reanalysis values.
We can see in the plot how the 'tas' reaches a minimum in February/March.
We can work out the exact coordinates of the selected grid cell as follows:
```r
data$lon[60]
## 166.5
data$lat[30]
At this point, additional numerical checks could be carried out to ensure the
data retrieval stage has been sucessful.
The next common step is to perform basic statistics on the raw data such as
computing the climatologies with bias correction to account for the drift
dependency on the initial conditions in the verification stage.
This is straightforward with the functions `Clim()` and `Ano()`. Before,
however, it is usual to average the data over the region of interest instead
of focusing on a single grid cell. This can be done by directly requesting the
data to `Load()` with the 'areave' output type, which will apply area weights
and calculate total average:
```r
data <- Load('tas', list(expA, expB), list(obsX),
sdates = sdates,
leadtimemin = 2, leadtimemax = 7,
latmin = -10, latmax = 60,
lonmin = 100, lonmax = 250,
output = 'areave', grid = 't106grid',
method = 'distance-weighted')
```
It is also common at this point to use `Smoothing()` to calculate running
means along the lead-times dimension or `Trend()` to compute and remove the
trend of the raw data. ¿¿??
Then we calculate per-pair climatologies with bias correction with `Clim()`:
```r
clim <- Clim(mod, obs)
PlotClim(clim$clim_exp, clim$clim_obs, monini = 12,
toptitle = "Experiment A, B, Observation X: 'tas' over Atlantic region.",
ytitle = "K", fileout = 'clim_expA_expB_obsX.eps')
Each line in this plot represents the climatology of each member of the
experimental datasets. A single climatology could be obtained providing the
parameter memb = FALSE to `Clim()`.
We can obtain the anomalies with `Ano()`, which simply substracts the
climatologies to the raw data:
```r
ano_mod <- Ano(mod, clim$clim_exp)
ano_obs <- Ano(obs, clim$clim_obs)
PlotAno(ano_mod, ano_obs, gsub('1101', '1201', sdates),
toptitle = paste0(c("Experiment A", "Experiment B"),
": bias corrected 'tas' anomalies over Atlantic region."),
ytitle = 'K',
fileout = c('ano_expA_obsX.eps', 'ano_expB_obsX.eps'))
<img src="ano_expA_obsX.png" width="600" />
<img src="ano_expB_obsX.png" width="600" />
To complete the bias correction we would need to add the observed climatologies
to these anomalies. The package, however, uses as working units the anomalies
and observed climatologies separately.
Other basic statistics such as trends, detrended anomalies, seasonal means,
running means, EOF/PCA and empirical model generation are implemented in the
package. Check [**Basic statistics**](basic_statistics.md) for a thorough
explanation.