Perfect Prognosis Approach: Application to Seasonal Forecasts
downscaleR package for statistical downscaling, Santander MetGroup
Last update 18/08/2015
In this worked example, we will perform analog downscaling of the System4 seasonal forecast of summer (JJA) maximum surface temperature on six weather stations of the GSN network in Spain. NCEP reanalysis data will be used for training the model and then we will perform the test on the System4 model of 15 members.
1 Loading and preparing data
Here, test data loading is done accessing the ECOMS User Data Gateway via the loadECOMS function from package ecomsUDG.Racess (registration and log in to the UDG Thredds Administration Portal is needed). Further examples are given in the ECOMS-UDG wiki.
1.1 Observation data
In this example, VALUE ECA&D station data will be used as reference observation:
download.file("http://meteo.unican.es/work/downscaler/data/VALUE_ECA_86_v2.tar.gz",
destfile = "mydirectory/VALUE_ECA_86_v2.tar.gz")
# Extract files from the tar.gz file
untar("mydirectory/VALUE_ECA_86_v2.tar.gz", exdir = "mydirectory")
# Data inventory
value <- "mydirectory/VALUE_ECA_86_v2"
di <-(dataInventory(value))Next, we load station data within a given geographical bounding box, thus, the lonLim and latLim arguments are filled with a vector of length two, defining the corners of the bounding box. For instance:
obs <- loadStationData(dataset = value,
var = c( "tmean"),
lonLim = c(-11,5),
latLim = c(35,46),
season= 6:8,
years = 1991:2000)## [2015-10-28 10:17:53] Loading data ...
## [2015-10-28 10:17:54] Retrieving metadata ...
## [2015-10-28 10:17:54] Done.str(obs)## List of 5
## $ Variable:List of 1
## ..$ varName: chr "tmean"
## $ Data : num [1:920, 1:12] 17 16.6 18.6 16.1 14.2 17 15.5 17.6 18.6 17.3 ...
## ..- attr(*, "dimensions")= chr [1:2] "time" "station"
## $ xyCoords: num [1:12, 1:2] -6.73 -9.15 -6.83 -4.49 -4.01 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:12] "000212" "000214" "000229" "000231" ...
## .. ..$ : chr [1:2] "longitude" "latitude"
## $ Dates :List of 2
## ..$ start: chr [1:920] "1991-06-01 00:00:00" "1991-06-02 00:00:00" "1991-06-03 00:00:00" "1991-06-04 00:00:00" ...
## ..$ end : chr [1:920] "1991-06-02 00:00:00" "1991-06-03 00:00:00" "1991-06-04 00:00:00" "1991-06-05 00:00:00" ...
## $ Metadata:List of 4
## ..$ station_id: int [1:12] 212 214 229 231 232 234 236 355 800 1394 ...
## ..$ name : chr [1:12] "BRAGANCA" "LISBOA-GEOFISICA" "BADAJOZ-TALAVERALAREAL" "MALAGA" ...
## ..$ altitude : int [1:12] 690 77 185 7 1894 251 44 1567 151 370 ...
## ..$ source : chr [1:12] "ECA&D" "ECA&D" "ECA&D" "ECA&D" ...We can plot the annual time series using function getYearsAsINDEX. Here we plot the time series for each loaded station (legend = station ID):
# Plotting of annual time series
yrs <- getYearsAsINDEX(obs)
plot.data <- sapply(1:ncol(obs$Data), function(x) {tapply(obs$Data[ ,x], INDEX = yrs, FUN = mean, na.rm = TRUE)})
colnames(plot.data) <- obs$Metadata$station_id
plot(unique(yrs), plot.data[ ,1], ylim = c(floor(min(plot.data)) - 5, ceiling(max(plot.data))), ty = "b", ylab = "tasmax", xlab = "year")
for (i in 2:ncol(plot.data)) {
lines(unique(yrs), plot.data[ ,i], col = i, ty = 'b')
}
legend("bottom", colnames(plot.data), col = 1:6, pch = 21, ncol = 2)
title("Mean maximum temperature summer (JJA)")
1.2 Predictors
In this example, the NCEP reanalysis data will be used, that can be loaded accesing through the ECOMS User Data Gateway (ECOMS-UDG) or also accessing the netcdf files remotely through the corresponding ncml file (see section 2.2. Accessing and loading grid data). However, in this case, we are going to use a subdomain of the NCEP dataset encompassing the period 1961-2010 for the Iberian Peninsula, which is available in a tar.gz file for download:
download.file("http://meteo.unican.es/work/downscaler/data/Iberia_NCEP.tar.gz",
destfile = "data/datasets/Iberia_NCEP.tar.gz")
# Extract files from the tar.gz file
untar("mydirectory/Iberia_NCEP.tar.gz", exdir = "mydirectory")
# Define the path to the ncml file
ncep <- "mydirectory/Iberia_NCEP/Iberia_NCEP.ncml"The set of predictors (a.k.a. atmospheric pattern) P5 is next loaded for the training period. A single variable for a specific spatiotemporal domain is what we term a field in downscaleR terminology. In this case, for the same spatio-temporal domain we are loading two different variables, that can be merged together forming a multifield. A multifield can be loaded from a given dataset using the loadMultiField function, that takes care of regridding to ensure the spatial matching of the different input fields. (NOTE: it is also possible to load the different fields sepparately and then join them using the multifield constructor makeMultiField).
In this example, we are interested in the variables of the atmospheric pattern P5 (mean surface temperature and sea-level pressure) whose standard names are “tas” and “psl” respectively, as indicated in the vocabulary (see help("vocabulary") for details).
pred <- loadMultiField(ncep, vars = c("tas","psl"),
lonLim = c(-11,5),
latLim = c(35,46),
season = 6:8,
years = 1991:2000)## [2015-10-28 10:17:55] Loading predictor 1 (tas) out of 2
## [2015-10-28 10:17:56] Loading predictor 2 (psl) out of 2
## [2015-10-28 10:17:56] Using the original grid as no 'new.grid' has been introduced
## [2015-10-28 10:17:56] Done.# This is how the multifield looks like:
plotMeanField(pred)
Next, a PCA is performed on the calibration data retaining just the first 10 PCs (see also section 4.3. Principal Components (and EOFs)).
pred.pca <- prinComp(pred, n.eofs = 10)## [2015-10-28 10:17:56] Performing PC analysis on 2 variables plus their combination...
## [2015-10-28 10:17:56] Done1.3 Test data (Seasonal forecasting)
Next, seasonal forecasting data from the System4 model of 15 members (dataset = “System4 seasonal 15”) is loaded from the ECOMS User Data Gateway, considering the domain of analysis, the target season and the test period 2001-2010.
library(ecomsUDG.Raccess)
loginECOMS_UDG(username = "usr", password = "pswrd")
psl.s4 <- loadECOMS(dataset = "System4_seasonal_15",
var = "psl",
lonLim = c(-11,5),
latLim = c(35,46),
season = 6:8,
years = 2001:2010,
time = "DD",
aggr.d = "mean",
leadMonth = 2)
tas.s4 <- loadECOMS(dataset = "System4_seasonal_15",
var = "tas",
lonLim = c(-11,5),
latLim = c(35,46),
season = 6:8,
years = 2001:2010,
time = "DD",
aggr.d = "mean",
leadMonth = 2)# Example: visualization of a multimember
plotMeanField(psl.s4, multi.member = TRUE)
In the case of the predictors a multifield is loaded directly with downscaleR from the local dataset, in this other case, a single field (predictor) of test data is loaded each time from the ECOMS User Data Gateway, thus, we need to create a multifield joining the predictors with function makeMultiField. Prior to multifield creation, a re-gridding is performed with function interpGridData to the same grid of the predictors (function getGrid).
# Regridding to match predictors grid
psl.s4.2 <- interpData(psl.s4, new.Coordinates = getGrid(pred))## [2015-10-28 10:17:58] Performing multi-member nearest interpolation... may take a while
## [2015-10-28 10:17:58] Interpolating member 1 out of 15
## [2015-10-28 10:17:58] Interpolating member 2 out of 15
## [2015-10-28 10:17:59] Interpolating member 3 out of 15
## [2015-10-28 10:17:59] Interpolating member 4 out of 15
## [2015-10-28 10:17:59] Interpolating member 5 out of 15
## [2015-10-28 10:18:00] Interpolating member 6 out of 15
## [2015-10-28 10:18:00] Interpolating member 7 out of 15
## [2015-10-28 10:18:00] Interpolating member 8 out of 15
## [2015-10-28 10:18:00] Interpolating member 9 out of 15
## [2015-10-28 10:18:01] Interpolating member 10 out of 15
## [2015-10-28 10:18:01] Interpolating member 11 out of 15
## [2015-10-28 10:18:01] Interpolating member 12 out of 15
## [2015-10-28 10:18:02] Interpolating member 13 out of 15
## [2015-10-28 10:18:02] Interpolating member 14 out of 15
## [2015-10-28 10:18:02] Interpolating member 15 out of 15
## [2015-10-28 10:18:03] Donetas.s4.2 <- interpData(tas.s4, new.Coordinates = getGrid(pred))## [2015-10-28 10:18:03] Performing multi-member nearest interpolation... may take a while
## [2015-10-28 10:18:03] Interpolating member 1 out of 15
## [2015-10-28 10:18:03] Interpolating member 2 out of 15
## [2015-10-28 10:18:03] Interpolating member 3 out of 15
## [2015-10-28 10:18:03] Interpolating member 4 out of 15
## [2015-10-28 10:18:04] Interpolating member 5 out of 15
## [2015-10-28 10:18:04] Interpolating member 6 out of 15
## [2015-10-28 10:18:04] Interpolating member 7 out of 15
## [2015-10-28 10:18:05] Interpolating member 8 out of 15
## [2015-10-28 10:18:05] Interpolating member 9 out of 15
## [2015-10-28 10:18:05] Interpolating member 10 out of 15
## [2015-10-28 10:18:06] Interpolating member 11 out of 15
## [2015-10-28 10:18:06] Interpolating member 12 out of 15
## [2015-10-28 10:18:06] Interpolating member 13 out of 15
## [2015-10-28 10:18:06] Interpolating member 14 out of 15
## [2015-10-28 10:18:07] Interpolating member 15 out of 15
## [2015-10-28 10:18:07] Done# Multifield creation
mf.test <- makeMultiField(psl.s4.2, tas.s4.2)2 Analog downscaling
The analog method is applied with the analog function to downscale the System4 seasonal forecast. The arguments follow the order observations > predictors > test data. Further additional arguments control the different options of the analog construction (type help("analogs") for more details):
an <- analogs(obs = obs, pred = pred.pca, sim = mf.test, n.neigh = 15, sel.fun = "random")## Variables in sim reordered to match pred
## [2015-10-28 10:18:07] Calculating analogs ...
## [2015-10-28 10:18:11] Done.2.1 Validation of results
We can visualize the downscaled forecast skill using the tercileValidation function. To this aim, we load the observation data in the test period:
val <- loadStationData(dataset = value,
var = "tmax",
lonLim = c(-11,5),
latLim = c(35,46),
season = 6:8,
years = 2001:2010)## [2015-10-28 10:18:12] Loading data ...
## [2015-10-28 10:18:12] Retrieving metadata ...
## [2015-10-28 10:18:12] Done.# This is a plot that performs the mean of all stations
tercileValidation(an, obs = val)## Warning in tercileValidation(an, obs = val): The results presented are the
## mean of all input stations
# This example enables the visualization of results for a particular station (e.g. Braganca):
tercileValidation(an, obs = val, stationId = "212")
title("Braganca")
print(sessionInfo())## R version 3.2.2 (2015-08-14)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.3 LTS
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_ES.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=es_ES.UTF-8 LC_NAME=es_ES.UTF-8
## [9] LC_ADDRESS=es_ES.UTF-8 LC_TELEPHONE=es_ES.UTF-8
## [11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=es_ES.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] downscaleR_0.8-3 maps_2.3-10 downscaleR.java_1.0-0
## [4] rJava_0.9-7
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.0 knitr_1.10.5 magrittr_1.5
## [4] MASS_7.3-44 evd_2.3-0 munsell_0.4.2
## [7] colorspace_1.2-6 stringr_1.0.0 plyr_1.8.3
## [10] CircStats_0.2-4 fields_8.2-1 tools_3.2.2
## [13] grid_3.2.2 spam_1.0-1 dtw_1.17-1
## [16] htmltools_0.2.6 yaml_2.1.13 abind_1.4-3
## [19] digest_0.6.8 formatR_1.2 bitops_1.0-6
## [22] RCurl_1.95-4.7 evaluate_0.7 rmarkdown_0.6.1
## [25] proxy_0.4-15 stringi_0.4-1 scales_0.2.5
## [28] boot_1.3-17 verification_1.42