wiki:udg/ecoms/RPackage/examples/drift

Analysing model drift in South-western Iberia

In this practice, we will analyse the model drift by using the forecast of daily mean (time = "DD", aggr.d = "mean") surface temperature (var = "tas") for July (season = 7) 2006 (years = 2006) considering 6 different forecast (lead) months, from January to June. The lead month 0 (i.e., the initialization of July itself) will be the reference from which we will compute the anomalies (leadMonth = 0). In this example, we will consider the first member (members = 1) of the CFSv2 hindcast as reference:

ref <- loadECOMS(dataset = "CFSv2_seasonal", 
                 var = "tas",
                 members = 1,
                 lonLim = c(-10,-1),
                 latLim = c(36,40),
                 season = 7,
                 years = 2006,
                 leadMonth = 0,
                 time = "DD",
                 aggr.d = "mean")

## [2016-05-13 14:28:46] Defining homogeneization parameters for variable "tas"
## [2016-05-13 14:28:46] Opening dataset...
## [2016-05-13 14:28:56] The dataset was successfuly opened
## [2016-05-13 14:28:56] Defining geo-location parameters
## NOTE: 'leadMonth = 0' selected
## [2016-05-13 14:28:57] Defining initialization time parameters
## NOTE: Daily aggregation will be computed from 6-hourly data
## [2016-05-13 14:28:58] Retrieving data subset ...
## [2016-05-13 14:29:04] Done

Messages on-screen inform about the loading process. The selection of leadMonth = 0 will give a NOTE.

library(downscaleR)
plotMeanGrid(ref)
title(main = "Lead month 0 forecast of July 2001")
# This is the climatology of the reference field
ref.field <- apply(ref$Data, MARGIN = c(3,2), FUN = mean, na.rm = TRUE)

Next, we will load the forecast of the target variable recursively for lead month values from 1 to 6 (i.e.: the initializations from January to June). The different objects are arranged in a list:

cfs.list <- lapply(1:6, function(lead.month) {
      loadECOMS(dataset = "CFSv2_seasonal",
                var = "tas",
                members = 1,
                lonLim = c(-10,-1),
                latLim = c(36,40),
                season = 7,
                years = 2006,
                leadMonth = lead.month,
                time = "DD",
                aggr.d = "mean")
      }
)

In order to visualize the departures of each lead month from the reference in the same range of values, we will use the spplot method of package sp for plotting spatial objects. To this aim, we will first compute the multi-member spatial mean for each lead month forecast, and then we will arrange the data in a matrix of 6 columns (one for each month), and x * y rows, as follows:

# The library sp needs to be installed to do this example:
if (!require(sp)) install.packages(sp)
# Matrix of anomalies between lead month and reference 
aux.mat <- sapply(1:length(cfs.list), function(i) {apply(cfs.list[[i]]$Data, MARGIN = c(3,2), FUN = mean, na.rm = TRUE) - ref.field})
# 2D coordinates
xy <- expand.grid(ref$xyCoords$x, ref$xyCoords$y)
# This step ensures regularity of the CFS grid, which is not perfectly regular:
xy.coords <- coordinates(points2grid(points = SpatialPoints(xy), tolerance = .003))
# Trick to adequately adjust the coordinate ordering (thanks to Kathryn Nicklin for spotting a previous error here)
xy.coords[ ,2] <- rev(xy.coords[ ,2])
# Now we create a data.frame with the coordinates X-Y in the first two columns and the mean anomalies in the next 6 columns:
df <- cbind.data.frame(xy.coords, aux.mat)
names(df) <- c("x","y",paste("LeadMonth_",1:6, sep = ""))
str(df)
## 'data.frame':        55 obs. of  8 variables:
##  $ x          : num  -10.31 -9.38 -8.44 -7.5 -6.56 ...
##  $ y          : num  36.4 36.4 36.4 36.4 36.4 ...
##  $ LeadMonth_1: num  0.0596 0.1601 0.5955 0.9068 1.2601 ...
##  $ LeadMonth_2: num  -0.1215 0.0509 0.3622 0.5444 0.6977 ...
##  $ LeadMonth_3: num  -0.48 -0.359 -0.402 -0.693 -0.967 ...
##  $ LeadMonth_4: num  0.22303 0.27295 0.23869 0.04764 -0.00855 ...
##  $ LeadMonth_5: num  -1.212 -0.986 -0.682 -0.587 -0.584 ...
##  $ LeadMonth_6: num  -1.314 -0.732 -0.32 -0.5 -0.98 ...
coordinates(df) <- c(1,2)
gridded(df) <- TRUE
class(df)
## [1] "SpatialPixelsDataFrame"
## attr(,"package")
## [1] "sp"

In the next lines we use apply the spplot method of package sp, generating a lattice-type map. In first place, we will also load a SpatialLines dataset remotely stored at Santander Met Group server, in order to represent the coastline in the lattice map generated as a geographical reference:

load(url("http://meteo.unican.es/work/downscaler/aux/wlines.rda"), verbose = TRUE)
l1 <- list("sp.lines", wlines)
spplot(df,
       as.table = TRUE,
       col.regions = colorRampPalette(c("blue","white","red")),
       at = seq(-5.25,5.25,.25),
       scales = list(draw = TRUE),
       sp.layout = list(l1))

The results show how a varying bias of the forecast depending on the lead time, i.e., the model drift.

Finally, we display the spatial mean of the anomalies w.r.t. the reference for each lead month considered using a barplot:

barplot(colMeans(df@data), names.arg = abbreviate(names(df)), xlab = "lead month", ylab = "anomaly (ºC)")
title(main = "Mean forecast bias w.r.t. the lead-month 0 initialization")
mtext("Member 1")

Last modified 5 years ago Last modified on May 13, 2016 2:58:49 PM

Attachments (3)

Download all attachments as: .zip