Skip to content

Commit

Permalink
Merge pull request #73 from geco-bern/fix-cru-vpd
Browse files Browse the repository at this point in the history 
Fix CRU vpd (when vapr wasn't requested) and WFDEI WFDE5 bias-correction
  • Loading branch information
@fabern
fabern authored Jan 8, 2025
2 parents 845dc95 + 7dc7c2b commit adc1bad
Show file tree
Hide file tree
Showing 3 changed files with 225 additions and 272 deletions.
203 changes: 72 additions & 131 deletions R/ingest.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ ingest <- function(
patm_base <-patm_mean <- month <- tavg <-temp <- temp_fine <-
tmax <- tmax_fine <- tmin <- tmin_fine <- prec <- prec_fine <-
days_in_month <- rain <- snow <- srad <- srad_fine <- ppfd <-
ppfd_fine <- wind <- wind_fine <- qair <- vap <- vapr <- vapr_fine <-
ppfd_fine <- wind <- wind_fine <- qair <- vapr <- vapr_fine <-
ilon <- data <- yy <- mm <- co2_avg <- year <- . <- bias <- NULL

# Check: all sites are distinct wrt name, lon and lat
Expand Down Expand Up @@ -240,12 +240,12 @@ ingest <- function(

# this returns a flat data frame with data from all sites
ddf <- ingest_globalfields(
siteinfo,
source = source,
dir = dir,
getvars = getvars,
siteinfo = siteinfo,
source = source,
getvars = getvars,
dir = dir,
timescale = timescale,
verbose = FALSE
verbose = FALSE
)

# redo ingest_globalfields() if some sites were not extracted
Expand Down Expand Up @@ -309,24 +309,25 @@ ingest <- function(
}

# bias-correct atmospheric pressure - per default
if (!is.null(getvars)){
if ("patm" %in% getvars){

df_patm_base <- siteinfo %>%
dplyr::select(sitename, elv) %>%
mutate(patm_base = calc_patm(elv))

# scale patm with a factor so that mean(patm) corresponds to patm_base:
ddf <- ddf %>%
group_by(sitename) %>%
summarise(patm_mean = mean(patm, na.rm = TRUE)) %>%
left_join(df_patm_base, by = "sitename") %>%
mutate(scale = patm_base / patm_mean) %>%
right_join(ddf, by = "sitename") %>%
mutate(patm = patm * scale) %>%
dplyr::select(-patm_base, -elv, -patm_mean, -scale)

if ("patm" %in% getvars){
if (all(is.na(siteinfo$elv))){
stop("Aborting. Argument elv is missing.")
}

df_patm_base <- siteinfo %>%
dplyr::select(sitename, elv) %>%
mutate(patm_base = calc_patm(elv))

# scale patm with a factor so that mean(patm) corresponds to patm_base:
ddf <- ddf %>%
group_by(sitename) %>%
summarise(patm_mean = mean(patm, na.rm = TRUE)) %>%
left_join(df_patm_base, by = "sitename") %>%
mutate(scale = patm_base / patm_mean) %>%
right_join(ddf, by = "sitename") %>%
mutate(patm = patm * scale) %>%
dplyr::select(-patm_base, -elv, -patm_mean, -scale)

}

# bias-correct other variables form worldclim or (only vpd) other sources
Expand All @@ -346,13 +347,15 @@ ingest <- function(
if ("swin" %in% getvars){rlang::inform("Bias Correction: Not yet implemented for swin.")}
if ("lwin" %in% getvars){rlang::inform("Bias Correction: Not yet implemented for lwin.")}

df_fine <- ingest_globalfields(siteinfo,
source = "worldclim",
dir = settings$dir_bias,
getvars = NULL,
timescale = NULL,
verbose = FALSE,
layer = unique(getvars_wc))
df_fine <- ingest_globalfields(
siteinfo,
source = "worldclim",
dir = settings$dir_bias,
getvars = NULL,
timescale = NULL,
verbose = FALSE,
layer = unique(getvars_wc)
)

# Bias correction for temperature: subtract difference
if ("tavg" %in% getvars_wc){
Expand All @@ -366,9 +369,7 @@ ingest <- function(
mutate(month = as.integer(month)) %>%
rename(temp_fine = tavg) %>%
right_join(ddf %>%
dplyr::filter(
lubridate::year(date) %in% year_start_wc:year_end_wc
) %>%
dplyr::filter(lubridate::year(date) %in% year_start_wc:year_end_wc) %>%
mutate(month = lubridate::month(date)) %>%
group_by(sitename, month) %>%
summarise(temp = mean(temp, na.rm = TRUE)),
Expand All @@ -386,7 +387,7 @@ ingest <- function(
dplyr::select(-bias, -month)
}

# Bias correction for temperature: subtract difference
# Bias correction for minimum temperature: subtract difference
if ("tmin" %in% getvars_wc){
if (source == "cru"){ # no tmin or tmax in wwfd
df_bias <- df_fine %>%
Expand All @@ -399,28 +400,24 @@ ingest <- function(
mutate(month = as.integer(month)) %>%
rename(tmin_fine = tmin) %>%
right_join(ddf %>%
dplyr::filter(
lubridate::year(date) %in% year_start_wc:year_end_wc
) %>%
mutate(month = lubridate::month(date)) %>%
group_by(sitename, month) %>%
summarise(tmin = mean(tmin, na.rm = TRUE)),
dplyr::filter(lubridate::year(date) %in% year_start_wc:year_end_wc) %>%
mutate(month = lubridate::month(date)) %>%
group_by(sitename, month) %>%
summarise(tmin = mean(tmin, na.rm = TRUE)),
by = c("sitename", "month")) %>%
mutate(bias = tmin - tmin_fine) %>%
dplyr::select(-tmin, -tmin_fine)

# correct bias by month
ddf <- ddf %>%
mutate(month = lubridate::month(date)) %>%
left_join(
df_bias %>% dplyr::select(sitename, month, bias),
by = c("sitename", "month")) %>%
left_join(df_bias %>% dplyr::select(sitename, month, bias),
by = c("sitename", "month")) %>%
arrange(sitename, date) %>%
mutate(tmin = ifelse(is.na(bias), tmin, tmin - bias)) %>%
dplyr::select(-bias, -month)
}
}

}

# Bias correction for temperature: subtract difference
if ("tmax" %in% getvars_wc){
Expand All @@ -435,13 +432,11 @@ ingest <- function(
mutate(month = as.integer(month)) %>%
rename(tmax_fine = tmax) %>%
right_join(ddf %>%
dplyr::filter(
lubridate::year(date) %in% year_start_wc:year_end_wc
) %>%
mutate(month = lubridate::month(date)) %>%
group_by(sitename, month) %>%
summarise(tmax = mean(tmax, na.rm = TRUE)),
by = c("sitename", "month")) %>%
dplyr::filter(lubridate::year(date) %in% year_start_wc:year_end_wc) %>%
mutate(month = lubridate::month(date)) %>%
group_by(sitename, month) %>%
summarise(tmax = mean(tmax, na.rm = TRUE)),
by = c("sitename", "month")) %>%
mutate(bias = tmax - tmax_fine) %>%
dplyr::select(-tmax, -tmax_fine)

Expand All @@ -450,8 +445,7 @@ ingest <- function(
mutate(month = lubridate::month(date)) %>%
left_join(
df_bias %>% dplyr::select(sitename, month, bias),
by = c("sitename", "month")
) %>%
by = c("sitename", "month")) %>%
arrange(sitename, date) %>%
mutate(tmax = ifelse(is.na(bias), tmax, tmax - bias)) %>%
dplyr::select(-bias, -month)
Expand All @@ -462,7 +456,11 @@ ingest <- function(
if ("prec" %in% getvars_wc){
df_bias <- df_fine %>%
dplyr::select(sitename, starts_with("prec_")) %>%
tidyr::pivot_longer(cols = starts_with("prec_"), names_to = "month", values_to = "prec", names_prefix = "prec_") %>%
tidyr::pivot_longer(
cols = starts_with("prec_"),
names_to = "month",
values_to = "prec",
names_prefix = "prec_") %>%
mutate(month = as.integer(month)) %>%
rename(prec_fine = prec) %>%
mutate(prec_fine = prec_fine / lubridate::days_in_month(month)) %>% # mm/month -> mm/d
Expand All @@ -478,6 +476,7 @@ ingest <- function(

# correct bias by month
if (source == "watch_wfdei" || source == "wfde5"){
# scaling also snow and rain rates
ddf <- ddf %>%
mutate(month = lubridate::month(date)) %>%
left_join(df_bias %>% dplyr::select(sitename, month, scale), by = c("sitename", "month")) %>%
Expand All @@ -487,7 +486,6 @@ ingest <- function(
rain = ifelse(is.na(scale), rain, rain * scale),
snow = ifelse(is.na(scale), snow, snow * scale)) %>%
dplyr::select(-scale, -month)

} else {
ddf <- ddf %>%
mutate(month = lubridate::month(date)) %>%
Expand All @@ -496,17 +494,19 @@ ingest <- function(
mutate(scale = ifelse(is.infinite(scale), 0, scale)) %>%
mutate(prec = ifelse(is.na(scale), prec, prec * scale)) %>%
dplyr::select(-scale, -month)

}

}

# Bias correction for shortwave radiation: scale by ratio
if ("srad" %in% getvars_wc){
kfFEC <- 2.04
df_bias <- df_fine %>%
dplyr::select(sitename, starts_with("srad_")) %>%
tidyr::pivot_longer(cols = starts_with("srad_"), names_to = "month", values_to = "srad", names_prefix = "srad_") %>%
tidyr::pivot_longer(
cols = starts_with("srad_"),
names_to = "month",
values_to = "srad",
names_prefix = "srad_") %>%
mutate(month = as.integer(month)) %>%
rename(srad_fine = srad) %>%
mutate(ppfd_fine = 1e3 * srad_fine * kfFEC * 1.0e-6 / (60 * 60 * 24) ) %>% # kJ m-2 day-1 -> mol m−2 s−1 PAR
Expand All @@ -533,16 +533,15 @@ ingest <- function(
if ("wind" %in% getvars_wc){
df_bias <- df_fine %>%
dplyr::select(sitename, starts_with("wind_")) %>%
tidyr::pivot_longer(cols = starts_with("wind_"),
names_to = "month",
values_to = "wind",
names_prefix = "wind_") %>%
tidyr::pivot_longer(
cols = starts_with("wind_"),
names_to = "month",
values_to = "wind",
names_prefix = "wind_") %>%
mutate(month = as.integer(month)) %>%
rename(wind_fine = wind) %>%
right_join(ddf %>%
dplyr::filter(
lubridate::year(date) %in% year_start_wc:year_end_wc
) %>%
dplyr::filter(lubridate::year(date) %in% year_start_wc:year_end_wc) %>%
mutate(month = lubridate::month(date)) %>%
group_by(sitename, month) %>%
summarise(wind = mean(wind, na.rm = TRUE)),
Expand All @@ -559,31 +558,12 @@ ingest <- function(
mutate(scale = ifelse(is.infinite(scale), 0, scale)) %>%
mutate(wind = ifelse(is.na(scale), wind, wind * scale)) %>%
dplyr::select(-scale, -month)

}

# Bias correction for relative humidity (actually vapour pressure): scale
if ("vapr" %in% getvars_wc){

# calculate vapour pressure from specific humidity - needed for bias correction with worldclim data
if (source == "watch_wfdei" || source == "wfde5"){
# specific humidity (qair, g g-1) is read, convert to vapour pressure (vapr, Pa)
ddf <- ddf %>%
rowwise() %>%
dplyr::mutate(
vapr = calc_vp(qair = qair,
patm = patm)
) %>%
ungroup()

} else if (source == "cru"){
# vapour pressure is read from file, convert from hPa to Pa
ddf <- ddf %>%
dplyr::mutate(vapr = 1e2 * vap) %>%
dplyr::select(-vap)

}
if ("vapr" %in% getvars_wc){ # i.e. equivalent with: "if (vpd" %in% getvars)"

# a) correct vapor pressure (vapr, Pa) with worldclim data
df_bias <- df_fine %>%
dplyr::select(sitename, starts_with("vapr_")) %>%
tidyr::pivot_longer(
Expand Down Expand Up @@ -612,30 +592,21 @@ ingest <- function(
mutate(vapr = ifelse(is.na(scale), vapr, vapr * scale)) %>%
dplyr::select(-scale, -month)

}

# Calculate vapour pressure deficit from specific humidity
if ("vpd" %in% getvars){

# b) re-calculate vapour pressure deficit VPD from bias-corrected vapor pressure
# and bias-correcte tmep, tmin, or tmax.
# (Note: this overwrites the vpd that has been computed within `ingest_globalfields()`)
if (source == "watch_wfdei" || source == "wfde5"){
# use daily mean temperature
ddf <- ddf %>%
rowwise() %>%
dplyr::mutate(
vapr = calc_vp(
qair = qair,
patm = patm
),
vpd = calc_vpd(eact = vapr, tc = temp)) %>%
dplyr::mutate(vpd = calc_vpd(eact = vapr, tc = temp)) %>%
ungroup()

} else if (source == "cru"){
# use daily minimum and maximum temperatures
ddf <- ddf %>%
rowwise() %>%
dplyr::mutate(
vpd = calc_vpd(eact = vapr, tmin = tmin, tmax = tmax)
) %>%
dplyr::mutate(vpd = calc_vpd(eact = vapr, tmin = tmin, tmax = tmax)) %>%
ungroup()
}
}
Expand All @@ -646,36 +617,6 @@ ingest <- function(

}

} else {

# Calculate vapour pressure deficit from specific humidity
# this calculates this variable for cases where there is
# no bias correction

if ("vpd" %in% getvars){

if (source == "watch_wfdei" || source == "wfde5"){
# use daily mean temperature
ddf <- ddf %>%
rowwise() %>%
dplyr::mutate(
vapr = calc_vp(qair = qair, patm = patm),
vpd = calc_vpd(eact = vapr, tc = temp)
) %>%
ungroup()

} else if (source == "cru"){
# use daily minimum and maximum temperatures
ddf <- ddf %>%
rowwise() %>%
dplyr::mutate(
vpd = calc_vpd(eact = vapr, tmin = tmin, tmax = tmax)
) %>%
ungroup()
}

}

}

} else if (source == "gee"){
Expand Down
Loading

0 comments on commit adc1bad

Please sign in to comment.