= 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:


{{{#!text/R
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")
                 time = "DD")
}}}

Messages on-screen inform about the loading process. Note that the selection of `leadMonth = 0` will give a NOTE:

{{{
[2015-05-15 14:50:46] Defining homogeneization parameters for variable "tas"
NOTE: 'leadMonth = 0' selected
[2015-05-15 14:50:46] Defining geo-location parameters
[2015-05-15 14:50:46] Defining initialization time parameters
NOTE: Daily aggregation will be computed from 6-hourly data
[2015-05-15 14:50:49] Retrieving data subset ...
[2015-05-15 14:50:56] Done
}}}

{{{#!text/R
plotMeanField(ref)
title(main = "Lead month 0 forecast of July 2001")
# This is the spatial mean of the reference field
ref.field <- apply(ref$Data, MARGIN = c(3,2), FUN = mean, na.rm = TRUE)
}}}

[[Image(image-20140903-175339.png)]]


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:


{{{#!text/R
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:

{{{#!text/R
# The library sp needs to be installed to do this example:
require(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 Kathryn!)
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)
}}}

Which returns the new spatial object class:

{{{
[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:

{{{#!text/R
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))
}}}

[[Image(image-20150130-185813.png)]]

The results show how a increasing lead month leads to a negative bias of the forecast, demonstrating that the mean state of a variable of a forecast is not stationary along the runtime dimension.

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

{{{#!text/R
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")
}}}

[[Image(image-20140903-181520.png)]]