Bias_Correction_practice
Learning goals:
The goal of this practice is to apply and analyze the bias correction techniques shown in the previous demonstration to calibrate the direct output of the seasonal forecasting. To this aim, the user should follow the next steps:
- Loading data of various types from local and remote locations using the downscaleR and ecomsUDG.Raccess tools to this aim.
- Homogenize the data obtained in the previous step.
- Apply the bias correction techniques and analyze their effect on the direct output of the model in terms of the differences between both, raw and calibrated, forecast and/or the spread of the predictions.
- Visualization of the bias corrected forecast skill and basic assessment of the quality of the predictions
Aim of the practice:
The objective of this practice is to work with and to analyze the different bias correction techniques included in the loadECOMS and downscaleR packages. To this aim, we will use the functions, observations and forecast used in previous demonstration and practices.
Note that in this case we will use the datasets included in the ECOMS-UDG data server and then we must to load/install the library ecomsUDG.Raccess and to be logged using the username and password of the Thredds Administration Panel (TAP):
devtools::install_github(c("SantanderMetGroup/downscaleR.java@stable","SantanderMetGroup/downscaleR@stable","SantanderMetGroup/ecomsUDG.Raccess@stable"))
library(downscaleR)
library(ecomsUDG.Raccess)
loginECOMS_UDG("username", "password")Data needed and methods:
- Target variable: Mean surface daily temperature (tmean)
- Observations: WATCH-Forcing-Data-ERA-Interim (WFDEI)
- Global Circulation Model: The NCEP Climate Forecast System Version 2 (CFSv2).
- Target season: Boreal winter (DJF)
- Lead Month: 3 months, which corresponds with the initialization of September.
- Calibration (training) period: 1991-2000
- Test period: 2001-2010
- Bias Correction methods: All the methods applicable for mean temperature (unbiasing, qq-map, ISI-MIP, etc.)
Practice Steps:
1 Loading the observational dataset and use some graphics to make a first preliminary analisys of the observations. In this case we will use as observations the WFDEI, which is a predefined dataset of the ECOMS User Data Gateway (UDG) Data Server. Load the grid points within the domain of analysis for the calibration period and the target season using the appropriate arguments of the loadECOMS function.
# For some graphics it could be useful the library fields
library(fields) # e.g: using the fields library
# Loading the observations
obs <- loadECOMS(dataset = "WFDEI", var = "tas", lonLim = c(-15,35), latLim = c(32, 75), season = c(12,1,2), years = 1991:2000)
# Maps: Mean and standard deviation of the observations
par(mfrow = c(1,2))
plotMeanField(obs, multi.member = FALSE)
obsStd <- t(apply(obs$Data, FUN = sd, MARGIN = c(2,3), na.rm = TRUE))
image.plot(obs$xyCoords$x,obs$xyCoords$y,obsStd, asp = 1, main = "Standard Deviation")
world(add = TRUE)
2 Load the seasonal forecast from the ECOMS User Data Gateway, considering again the domain of analysis, the target season and the calibration period. Explore the returned objects and estimate the bias of the forecast in the control period.
# Loading the seasonal forecast
prd <- loadECOMS(dataset = "CFSv2_seasonal_16", var = "tas", members = 1, lonLim = c(-15,35), latLim = c(32, 75), season = c(12,1,2), years = 1991:2000, leadMonth = 3, time = "DD")
# REMEMBER: We should interpolate the observations to the grid of model: we use the method "nearest" and the getGrid function to ensure spatial consistency:
obs <- interpGridData(obs, new.grid = getGrid(prd), method = "nearest")
par(mfrow = c(1,3))
plotMeanField(obs, multi.member = FALSE)
plotMeanField(prd, multi.member = FALSE)
obsMean <- t(apply(obs$Data, FUN = sd, MARGIN = c(2,3), na.rm = TRUE))
prdMean <- t(apply(prd$Data, FUN = sd, MARGIN = c(2,3), na.rm = TRUE))
image.plot(obs$xyCoords$x,obs$xyCoords$y,prdMean-obsMean, asp = 1, main = "Bias")
world(add = TRUE)
3 Load the seasonal forecast from the ECOMS User Data Gateway, considering again the domain of analysis, the target season and an unobserved/test period. Is there any difference between the control and the test periods?
# Loading an unobserved periods of the simulation:
sim <- loadECOMS(dataset = "CFSv2_seasonal_16", var = "tas", members = 1, lonLim = c(-15,35), latLim = c(32, 75), season = c(12,1,2), years = 2001:2010, leadMonth = 3, time = "DD")
simMean <- t(apply(sim$Data, FUN = mean, MARGIN = c(2,3), na.rm = TRUE))
prdMean <- t(apply(prd$Data, FUN = mean, MARGIN = c(2,3), na.rm = TRUE))
simStd <- t(apply(sim$Data, FUN = sd, MARGIN = c(2,3), na.rm = TRUE))
prdStd <- t(apply(prd$Data, FUN = sd, MARGIN = c(2,3), na.rm = TRUE))
par(mfrow = c(2,3))
image.plot(prd$xyCoords$x,prd$xyCoords$y,prdMean, asp = 1, main = "Climatology")
world(add = TRUE)
image.plot(prd$xyCoords$x,prd$xyCoords$y,simMean, asp = 1, main = "Climatology")
world(add = TRUE)
image.plot(prd$xyCoords$x,prd$xyCoords$y,simMean-prdMean, asp = 1, main = "Delta")
world(add = TRUE)
image.plot(prd$xyCoords$x,prd$xyCoords$y,prdStd, asp = 1, main = "Climatology")
world(add = TRUE)
image.plot(prd$xyCoords$x,prd$xyCoords$y,simStd, asp = 1, main = "Climatology")
world(add = TRUE)
image.plot(prd$xyCoords$x,prd$xyCoords$y,simStd-prdStd, asp = 1, main = "Delta")
world(add = TRUE)
4 Apply the bias correction techniques and compare the different results.
# In the case of the temperature these are the methods available:
sim1 <- biasCorrection (obs, prd, sim, method = "unbiasing")# bias
sim2 <- biasCorrection (obs, prd, sim, method = "qqmap")# qq-mapping
sim3 <- isimip(obs, prd, sim)# isimip
par(mfrow = c(2,2))
plotMeanField(sim, multi.member = FALSE)
simMean <- t(apply(sim1$Data, FUN = mean, MARGIN = c(2,3), na.rm = TRUE))
image.plot(prd$xyCoords$x,prd$xyCoords$y,simMean, asp = 1, main = "Unbiasing")
world(add = TRUE)
simMean <- t(apply(sim2$Data, FUN = mean, MARGIN = c(2,3), na.rm = TRUE))
image.plot(prd$xyCoords$x,prd$xyCoords$y,simMean, asp = 1, main = "QQ-Mapping")
world(add = TRUE)
simMean <- t(apply(sim3$Data, FUN = mean, MARGIN = c(2,3), na.rm = TRUE))
image.plot(prd$xyCoords$x,prd$xyCoords$y,simMean, asp = 1, main = "ISI-MIP")
world(add = TRUE)
We compare the results in terms of the mean and variability differences
simMean <- t(apply(sim$Data, FUN = mean, MARGIN = c(2,3), na.rm = TRUE))
simStd <- t(apply(sim$Data, FUN = sd, MARGIN = c(2,3), na.rm = TRUE))
par(mfrow = c(2,4))
image.plot(prd$xyCoords$x,prd$xyCoords$y,simMean, asp = 1, main = "Climatology")
world(add = TRUE)
sim1Mean <- t(apply(sim1$Data, FUN = mean, MARGIN = c(2,3), na.rm = TRUE))
image.plot(prd$xyCoords$x,prd$xyCoords$y,sim1Mean-simMean, asp = 1, main = "Unbiasing")
world(add = TRUE)
sim1Mean <- t(apply(sim2$Data, FUN = mean, MARGIN = c(2,3), na.rm = TRUE))
image.plot(prd$xyCoords$x,prd$xyCoords$y,sim1Mean-simMean, asp = 1, main = "QQ-Mapping")
world(add = TRUE)
sim1Mean <- t(apply(sim3$Data, FUN = mean, MARGIN = c(2,3), na.rm = TRUE))
image.plot(prd$xyCoords$x,prd$xyCoords$y,sim1Mean-simMean, asp = 1, main = "ISI-MIP")
world(add = TRUE)
image.plot(prd$xyCoords$x,prd$xyCoords$y,simStd, asp = 1, main = "Standard Deviation")
world(add = TRUE)
sim1Std <- t(apply(sim1$Data, FUN = sd, MARGIN = c(2,3), na.rm = TRUE))
image.plot(prd$xyCoords$x,prd$xyCoords$y,sim1Std/simStd, asp = 1, main = "Unbiasing")
world(add = TRUE)
sim1Std <- t(apply(sim2$Data, FUN = sd, MARGIN = c(2,3), na.rm = TRUE))
image.plot(prd$xyCoords$x,prd$xyCoords$y,sim1Std/simStd, asp = 1, main = "QQ-Mapping")
world(add = TRUE)
sim1Std <- t(apply(sim3$Data, FUN = sd, MARGIN = c(2,3), na.rm = TRUE))
image.plot(prd$xyCoords$x,prd$xyCoords$y,sim1Std/simStd, asp = 1, main = "ISI-MIP")
world(add = TRUE)
In order to show the spread of the spatial mean of the ensemble of bias correction methodologies we could draw the plume of the corrected seasonal forecast with the function polygon, including the multi-method mean (black line) and the direct output (red line).
spatialMeanSim <- apply(sim$Data, MARGIN = 1, FUN = mean, na.rm = TRUE)
spatialMean <- matrix(data = NA, ncol = 3, nrow = length(spatialMeanSim))
spatialMean[,1] <- apply(sim1$Data, MARGIN = 1, FUN = mean, na.rm = TRUE)
spatialMean[,2] <- apply(sim2$Data, MARGIN = 1, FUN = mean, na.rm = TRUE)
spatialMean[,3] <- apply(sim3$Data, MARGIN = 1, FUN = mean, na.rm = TRUE)
ens.mean <- rowMeans(spatialMean)
dates <- as.POSIXlt(sim$Dates$start, tz="GMT")
par(mfrow = c(1,1))
plot(dates, ens.mean, ylim = range(spatialMean), ty = 'l', ylab = "Spatial Mean", xlab = "time")
polygon(x = c(dates, rev(dates)), y = c(apply(spatialMean, MARGIN = 1, FUN = min, na.rm = TRUE), rev(apply(spatialMean, MARGIN = 1, FUN = max, na.rm = TRUE))), border = "transparent", col = rgb(0,0,1,.4))
lines(dates, spatialMeanSim, col = "red")
5 Extend the previous analysis to the multi-member case.
# We can consider a multi-member forecast: loading the test/control period
prd <- loadECOMS(dataset = "CFSv2_seasonal_16", var = "tas", members = 1:9, lonLim = c(-15,35), latLim = c(32, 75), season = c(12,1,2), years = 1991:2000, leadMonth = 3, time = "DD")
# Function plotMeanField let us to draw all the members in a same figure:
plotMeanField(prd, multi.member = TRUE)
The multi-member approach let us to analyze the skill of the raw forecast and compare it with the skill obtained when the bias correction is applied. To this aim we can use the tercileValidation function:
# In the case of the temperature these are the methods available:
prd1 <- biasCorrection (obs, prd, prd, method = "unbiasing")# bias
prd2 <- biasCorrection (obs, prd, prd, method = "qqmap")# qq-mapping
prd3 <- isimip(obs, prd, prd)# isimip
# This is a plot that performs the mean of all stations
tercileValidation(prd, obs = obs)# Raw skill
tercileValidation(prd1, obs = obs)
We apply the correction in the test period and plot the different members:
# Loading an unobserved periods of the simulation:
sim <- loadECOMS(dataset = "CFSv2_seasonal_16", var = "tas", members = 1:9, lonLim = c(-15,35), latLim = c(32, 75), season = c(12,1,2), years = 2001:2010, leadMonth = 3, time = "DD")
plotMeanField(sim, multi.member = TRUE)
# In the case of the temperature these are the methods available:
sim1 <- biasCorrection (obs, prd, sim, method = "unbiasing")# bias
plotMeanField(sim1, multi.member = TRUE)
Key points: Identify the effect of the different bias correction techniques and point out their main advantages and shortcomings.
