-
Notifications
You must be signed in to change notification settings - Fork 10
/
fhempel.r
105 lines (103 loc) · 4.01 KB
/
fhempel.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
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
################################################################################
# FUNCTION TO APPLY THE CALIBRATION PROCEDURE DESCRIBED IN:
# Hempel et al. A trend-preserving bias correction - the ISI-MIP approach.
# Earth Syst Dynam, 2013;4:219-236.
#
# ARGUMENTS:
# - obs: DATE + SINGLE OBSERVED HISTORICAL SERIES
# - mod: DATE + OPTIONALLY MULTIPLE MODELLED SERIES
# - add: LOGICAL. IF TRUE, THE ADDITIVE CORRECTION IS APPLIED
# - mult: LOGICAL. IF TRUE, THE MULTIPLICATIVE CORRECTION IS APPLIED
# - output: RETURN THE SERIES ("series") OR THE PARAMETERS ("correction")
#
# Author: Antonio Gasparrini - GNU General Public License (version 3)
# Update: 24 April 2023
################################################################################
fhempel <- function(obs ,mod, add=TRUE, mult=TRUE, output="series") {
#
# CHECK THE OUTPUT
output <- match.arg(output,c("series","correction"))
#
# CHECK CONSISTENCY
if(ncol(obs)!=2) stop("obs must have two variables (date and series")
if(ncol(mod)<2) stop("mod must have at least two variables (date and series")
if(!any(class(obs[[1]])%in%"Date") | !any(class(mod[[1]])%in%"Date"))
stop("Class of first variable of 'obs' and 'mod' must be 'Date'")
#
# APPLY TO EACH MODELLED SERIES
out <- lapply(seq(ncol(mod)-1),function(j) {
#
# IDENTIFY PERIOD WITH NO MISSING
ind <- obs[[1]][obs[[1]] %in% mod[[1]]]
notna <- complete.cases(obs[obs[[1]]%in%ind,2],mod[mod[[1]]%in%ind,j+1])
indobs <- seq(nrow(obs))[obs[[1]]%in%ind][notna]
indmod <- seq(nrow(mod))[mod[[1]]%in%ind][notna]
if(length(indobs)==0) stop("no common non-missing days in 'obs' and 'mod'")
#
# EXTRACT MONTH AND YEAR
month <- as.numeric(format(obs[[1]][indobs],format="%m"))
if(length(unique(month))!=12)
stop("some months are not reprensented in the overlapping period")
year <- as.numeric(format(obs[[1]][indobs],format="%Y"))
monthyear <- factor(paste(year,month,sep="-"))
#
# COMPUTE ADDITIVE CORRECTION
mavgobs <- tapply(obs[[2]][indobs],month,mean,na.rm=T)
mavgmod <- tapply(mod[[j+1]][indmod],month,mean,na.rm=T)
C <- mavgobs - mavgmod
if(!add) C[] <- 0
#
# RESIDUALS FROM MONTHLY/YEAR AVERAGES, THEN MULTIPLICATIVE CORRECTION
myavgobs <- tapply(obs[[2]][indobs],monthyear,mean,na.rm=T)
myavgmod <- tapply(mod[[j+1]][indmod],monthyear,mean,na.rm=T)
resobs <- obs[[2]][indobs] - myavgobs[monthyear]
resmod <- mod[[j+1]][indmod] - myavgmod[monthyear]
B <- sapply(seq(12),
function(m) coef(lm(sort(resobs[month==m])~0+sort(resmod[month==m]))))
if(!mult) B[] <- 1
#
# RETURN CORRECTION IF STATED
if(output=="correction") return(matrix(c(C,B),ncol=2,
dimnames=list(unique(months(mod[,1],abbr=TRUE)),c("add","mult"))))
#
# OBTAIN DAY, MONTH AND YEAR FROM DATE (ORIGINAL SERIES)
day <- as.numeric(format(mod[[1]],format="%d"))
month <- as.numeric(format(mod[[1]],format="%m"))
year <- as.numeric(format(mod[[1]],format="%Y"))
monthyear <- factor(paste(year,month,sep="-"))
#
# DERIVE THE QUANTITIES TO REMOVE DISCONTINUITIES FROM CORRECTION
# (SEE Hempel et al, PAGE 228)
nday <- tapply(mod[[j+1]],monthyear,length)
d <- (day-1)/(nday[monthyear]-1)-0.5
dm <- 0.5*(abs(d)-d)
d0 <- 1-abs(d)
dp <- 0.5*(abs(d)+d)
#
# DERIVE THE CORRECTIONS ACCOUNTING FOR DISCONTINUITIES
Cm <- C[c(12,1:11)] ; Bm <- B[c(12,1:11)]
Cp <- C[c(2:12,1)] ; Bp <- B[c(2:12,1)]
C <- dm*Cm[month] + d0*C[month] + dp*Cp[month]
B <- dm*Bm[month] + d0*B[month] + dp*Bp[month]
#
# DERIVE THE CORRECTED SERIES
myavgout <- tapply(mod[[j+1]],monthyear,mean,na.rm=T)
resout <- mod[[j+1]] - myavgout[monthyear]
#
# RETURN CORRECTED SERIES IF STATED
return(myavgout[monthyear]+C+resout*B)
})
#
# RETURN CORRECTION IF STATED
if(output=="correction") {
names(out) <- names(mod)[-1]
return(out)
}
#
# ADD DATE AND RENAME
out <- cbind(mod[,1],as.data.frame(do.call(cbind,out)))
dimnames(out) <- dimnames(mod)
#
# RETURN
return(out)
}