Created
December 5, 2022 09:13
-
-
Save khufkens/1119a1db7ef8066977ec6a81909ab407 to your computer and use it in GitHub Desktop.
Format rsofun drivers
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#' Format and run p-model drivers on Euler | |
#' | |
#' Takes site information for a single | |
#' or multiple sites and grabs all data | |
#' required for a p-model run in rsofun. | |
#' | |
#' Parameter settings are provided as | |
#' arguments, but could be altered after | |
#' the fact if desired. | |
#' | |
#' @param siteinfo data frame using minimum information required | |
#' being five columns: sitename, lon, lat, year_start, year_end | |
#' @param params_siml simulation parameters (preset) | |
#' @param params_modl model parameters (preset) | |
#' @param df_soiltexture soil data specifics (preset) | |
#' @param bias_correction bolean TRUE/FALSE, when using high resolution | |
#' locations set the bias correction to TRUE, for global runs leave FALSE | |
#' (default) | |
#' @param verbose provide verbose output (default = FALSE) | |
#' @param run_model shoud model evaluation be run, which returns evaluated | |
#' data instead of driver files? | |
#' | |
#' @return returns an rsofun compatible driver file for the provided | |
#' sites | |
format_drivers <- function( | |
siteinfo, | |
params_siml = list( | |
spinup = TRUE, # to bring soil moisture to steady state | |
spinupyears = 10, # 10 is enough for soil moisture. | |
recycle = 1, # number of years recycled during spinup | |
soilmstress = FALSE, # soil moisture stress function is included | |
tempstress = FALSE, # temperature stress function is included | |
calc_aet_fapar_vpd = FALSE, # set to FALSE - should be dropped again | |
in_ppfd = TRUE, # if available from forcing files, set to TRUE | |
in_netrad = FALSE, # if available from forcing files, set to TRUE | |
outdt = 1, | |
ltre = FALSE, | |
ltne = FALSE, | |
ltrd = FALSE, | |
ltnd = FALSE, | |
lgr3 = TRUE, | |
lgn3 = FALSE, | |
lgr4 = FALSE | |
), | |
params_modl = list( | |
kphio = 0.09423773, | |
soilm_par_a = 0.33349283, | |
soilm_par_b = 1.45602286 | |
), | |
df_soiltexture = bind_rows( | |
top = tibble( | |
layer = "top", | |
fsand = 0.4, | |
fclay = 0.3, | |
forg = 0.1, | |
fgravel = 0.1), | |
bottom = tibble( | |
layer = "bottom", | |
fsand = 0.4, | |
fclay = 0.3, | |
forg = 0.1, | |
fgravel = 0.1) | |
), | |
bias_correction = FALSE, | |
verbose = TRUE, | |
run_model = FALSE | |
){ | |
#---- startup checks ---- | |
# bail if not on euler | |
if(!grepl('eu-', Sys.info()['nodename'])){ | |
stop("You are not on Euler, source data unavailable - abort abort abort!") | |
} | |
# bail if not on euler | |
if(!dir.exists("~/data")){ | |
stop("Data path is not linked, create a soft link in your home directory | |
to setup a link to the (CES) Euler data storage: | |
ln -s /cluster/work/climate/bestocke/data data | |
") | |
} | |
# check format of the siteinfo | |
names(siteinfo) %in% c("sitename","lon","lat","start_year","end_year","elv") | |
if(!dir.exists("~/data")){ | |
stop("Data path is not linked, create a soft link in your home directory | |
to setup a link to the (CES) Euler data storage: | |
ln -s /cluster/work/climate/bestocke/data data | |
") | |
} | |
# some feedback on the processing | |
if(verbose){ | |
message("Running on Euler, data linkages in place. Proceeding ....") | |
} | |
#---- complement siteinfo with WHC based on S_CWDX80 ---- | |
# some feedback on the processing | |
if(verbose){ | |
message("Processing WHC data ....") | |
} | |
# bail if not on euler | |
if(grepl('eu-', Sys.info()['nodename'])){ | |
if(verbose){ | |
message("Processing cwdx80 data on server ....") | |
} | |
if(bias_correction){ | |
filn <- "~/data/mct_data/cwdx80.nc" | |
siteinfo <- siteinfo %>% | |
left_join(rbeni::extract_nc( | |
dplyr::select(siteinfo, | |
sitename, | |
lon, | |
lat), | |
filn) %>% | |
unnest(data) %>% | |
rename(whc = V1), | |
by = c("sitename", "lon", "lat") | |
) | |
} else { | |
if(verbose){ | |
message("Processing cwdx80 data from repo ....") | |
} | |
# use local data for global runs | |
filn <- "data/mct/whc_0.5deg.tif" | |
if(!file.exists(filn)){ | |
stop("Please generate 0.5 degree WHC values using: | |
data-raw/mct/00_downsample_hwc.R in the project | |
directory") | |
} | |
siteinfo <- siteinfo %>% | |
left_join(rbeni::extract_nc( | |
dplyr::select(siteinfo, | |
sitename, | |
lon, | |
lat), | |
filn) %>% | |
unnest(data) %>% | |
rename(whc = V1), | |
by = c("sitename", "lon", "lat") | |
) | |
} | |
## fill gaps in whc | |
whc_median <- median(siteinfo$whc, na.rm = TRUE) | |
siteinfo <- siteinfo %>% | |
mutate( | |
whc = ifelse( | |
is.na(whc), | |
whc_median, | |
whc | |
)) | |
} | |
#---- grabbing WATCH data ---- | |
# some feedback on the processing | |
if(verbose){ | |
message("Processing WATCH data ....") | |
} | |
if(bias_correction){ | |
ddf_watch <- ingest( | |
siteinfo = siteinfo, | |
source = "watch_wfdei", | |
getvars = c("temp", "prec", "ppfd", "vpd", "patm"), | |
dir = "~/data/watch_wfdei/", | |
settings = list( | |
correct_bias = "worldclim", | |
dir_bias = "~/data/worldclim") | |
) | |
} else { | |
if(verbose){ | |
message("Processing WATCH data coarse levels ....") | |
} | |
ddf_watch <- ingest( | |
siteinfo = siteinfo, | |
source = "watch_wfdei", | |
getvars = c("temp", "prec", "ppfd", "vpd", "patm"), | |
dir = "~/data/watch_wfdei/", | |
settings = list(correct_bias = NULL) | |
) | |
if(verbose){ | |
print(head(ddf_watch$data[[1]])) | |
} | |
} | |
# memory intensive, purge memory | |
gc() | |
#---- Processing CRU data ---- | |
if(verbose){ | |
message("Processing CRU data ....") | |
} | |
# adjust this with your local path | |
ddf_cru <- ingest( | |
siteinfo = siteinfo, | |
source = "cru", | |
getvars = "ccov", | |
dir = "~/data/cru/ts_4.01/", | |
settings = list(correct_bias = NULL) | |
) | |
# memory intensive, purge memory | |
gc() | |
#---- Merging climate data ---- | |
if(verbose){ | |
message("Merging climate data ....") | |
} | |
# merge all climate drivers into | |
# one format | |
ddf_meteo <- ddf_watch %>% | |
tidyr::unnest(data) %>% | |
left_join( | |
ddf_cru %>% | |
tidyr::unnest(data), | |
by = c("sitename", "date") | |
) %>% | |
group_by(sitename) %>% | |
tidyr::nest() | |
#---- Append CO2 data ---- | |
if(verbose){ | |
message("Append CO2 data ....") | |
} | |
# grab the CO2 data matching date ranges | |
df_co2 <- ingest( | |
siteinfo, | |
source = "co2_cmip", | |
verbose = FALSE, | |
dir = "~/data/co2/" | |
) | |
#---- Append FAPAR data ---- | |
if(verbose){ | |
message("Append FAPAR data ....") | |
} | |
# grab the FAPAR data | |
ddf_fapar_unity <- ingest( | |
siteinfo = siteinfo, | |
source = "fapar_unity" | |
) | |
#---- Format p-model driver data ---- | |
if(verbose){ | |
message("Combining all driver data ....") | |
} | |
output <- collect_drivers_sofun( | |
site_info = siteinfo, | |
params_siml = params_siml, | |
meteo = ddf_meteo, | |
fapar = ddf_fapar_unity, | |
co2 = df_co2, | |
params_soil = df_soiltexture | |
) | |
if (run_model){ | |
message("Running model on formatted data ....") | |
output <- process_pmodel_data( | |
file = output | |
) | |
} | |
# return data, either a driver | |
# or processed output | |
return(output) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment