From 4f5397bd4e906698e38c47054f52dd941038accb Mon Sep 17 00:00:00 2001 From: Beni Stocker Date: Thu, 6 Jun 2024 13:48:07 +0200 Subject: [PATCH] successful test at single site: correct forcing, p-model runs --- DESCRIPTION | 1 + R/calc_patm.R | 47 ++++++++++ R/calc_vp.R | 42 ++++----- R/calc_vpd.R | 58 ++++++------ R/grsofun_run.R | 133 +++++++++++++++++++++++++--- vignettes/global_rsofun_example.Rmd | 1 + 6 files changed, 221 insertions(+), 61 deletions(-) create mode 100644 R/calc_patm.R diff --git a/DESCRIPTION b/DESCRIPTION index 1781e5f..90a4a4d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,6 +13,7 @@ Description: A wrapper for global rsofun runs License: AGPL-3 Encoding: UTF-8 Imports: + zoo, dplyr, tidyr, magrittr, diff --git a/R/calc_patm.R b/R/calc_patm.R new file mode 100644 index 0000000..9123e8d --- /dev/null +++ b/R/calc_patm.R @@ -0,0 +1,47 @@ +#' Calculates atmospheric pressure +#' +#' Calculates atmospheric pressure as a function of elevation, by default assuming +#' standard atmosphere (101325 Pa at sea level) +#' +#' @param elv Elevation above sea-level (m.a.s.l.) +#' @param patm0 (Optional) Atmospheric pressure at sea level (Pa), defaults to 101325 Pa. +#' +#' @details The elevation-dependence of atmospheric pressure is computed by +#' assuming a linear decrease in temperature with elevation and a mean +#' adiabatic lapse rate (Berberan-Santos et al., 1997): +#' \deqn{ +#' p(z) = p0 ( 1 - Lz / TK0) ^ ( g M / (RL) ) +#' } +#' where \eqn{z} is the elevation above mean sea level (m, argument \code{elv}), +#' \eqn{g} is the gravity constant (9.80665 m s-2), \eqn{p0} is the atmospheric +#' pressure at 0 m a.s.l. (argument \code{patm0}, defaults to 101325 Pa), +#' \eqn{L} is the mean adiabatic lapse rate (0.0065 K m-2), +#' \eqn{M} is the molecular weight for dry air (0.028963 kg mol-1), +#' \eqn{R} is the universal gas constant (8.3145 J mol-1 K-1), and \eqn{TK0} +#' is the standard temperature (298.15 K, corresponds to 25 deg C). +#' +#' @return A numeric value for \eqn{p} +#' +#' @examples print("Standard atmospheric pressure, in Pa, corrected for 1000 m.a.s.l.:") +#' print(calc_patm(1000)) +#' +#' @references Allen, R. G., Pereira, L. S., Raes, D., Smith, M.: +#' FAO Irrigation and Drainage Paper No. 56, Food and +#' Agriculture Organization of the United Nations, 1998 +#' +#' @export +#' +calc_patm <- function( elv, patm0 = 101325 ){ + + # Define constants: + kTo <- 298.15 # base temperature, K (Prentice, unpublished) + kL <- 0.0065 # adiabiatic temperature lapse rate, K/m (Allen, 1973) + kG <- 9.80665 # gravitational acceleration, m/s^2 (Allen, 1973) + kR <- 8.3145 # universal gas constant, J/mol/K (Allen, 1973) + kMa <- 0.028963 # molecular weight of dry air, kg/mol (Tsilingiris, 2008) + + # Convert elevation to pressure, Pa: + patm <- patm0*(1.0 - kL*elv/kTo)^(kG*kMa/(kR*kL)) + + return(patm) +} diff --git a/R/calc_vp.R b/R/calc_vp.R index 2394c0c..c30ba78 100644 --- a/R/calc_vp.R +++ b/R/calc_vp.R @@ -4,49 +4,49 @@ #' #' @param qair Air specific humidity (g g-1) #' @param patm Atmospehric pressure (Pa) -#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing +#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing #' for calculating it based on standard sea level pressure) #' #' @return vapor pressure (Pa) #' @export -#' +#' calc_vp <- function( qair, patm = NA, elv = NA ) { - - # Ref: Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure - # Calculation Methods, in Evaporation and Evapotranspiration: + + # Ref: Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure + # Calculation Methods, in Evaporation and Evapotranspiration: # Measurements and Estimations, Springer, London. # vpd = 0.611*exp[ (17.27 tc)/(tc + 237.3) ] - ea # where: # tc = average daily air temperature, deg C # eact = actual vapor pressure, Pa - - # calculate atmopheric pressure (Pa) assuming standard + + # calculate atmopheric pressure (Pa) assuming standard # conditions at sea level (elv=0) # if (is.na(elv) && is.na(patm)){ - # + # # warning( # " - # calc_vp(): Either patm or elv must be provided + # calc_vp(): Either patm or elv must be provided # if eact is not given. # ") # vp <- NA - # + # # } else { patm <- ifelse(is.na(patm), calc_patm(elv), patm) - + # Calculate VP. vp <- calc_vp_inst(qair = qair, patm = patm) # } return( vp ) - + } #' Instantaneous vapour pressure @@ -63,26 +63,26 @@ calc_vp_inst <- function( qair, patm ){ - - # Ref: Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure - # Calculation Methods, in Evaporation and Evapotranspiration: + + # Ref: Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure + # Calculation Methods, in Evaporation and Evapotranspiration: # Measurements and Estimations, Springer, London. # vpd = 0.611*exp[ (17.27 tc)/(tc + 237.3) ] - ea # where: # tc = average daily air temperature, deg C # ea = actual vapor pressure, Pa - + kR = 8.3143 # universal gas constant, J/mol/K (Allen, 1973) kMv = 18.02 # molecular weight of water vapor, g/mol (Tsilingiris, 2008) kMa = 28.963 # molecular weight of dry air, g/mol (Tsilingiris, 2008) - + # calculate the mass mixing ratio of water vapor to dry air (dimensionless) wair <- qair / (1 - qair) - - # calculate water vapor pressure + + # calculate water vapor pressure rv <- kR / kMv rd <- kR / kMa eact = patm * wair * rv / (rd + wair * rv) - + return( eact ) -} \ No newline at end of file +} diff --git a/R/calc_vpd.R b/R/calc_vpd.R index 71b6909..285cfd2 100644 --- a/R/calc_vpd.R +++ b/R/calc_vpd.R @@ -6,15 +6,15 @@ #' @param qair Air specific humidity (g g-1) #' @param eact Water vapour pressure (Pa) #' @param tc temperature, deg C -#' @param tmin (optional) min daily air temp, deg C +#' @param tmin (optional) min daily air temp, deg C #' @param tmax (optional) max daily air temp, deg C #' @param patm Atmospehric pressure (Pa) -#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing +#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing #' for calculating it based on standard sea level pressure) #' #' @return vapor pressure deficit (Pa) #' @export -#' +#' calc_vpd <- function( qair = NA, eact = NA, @@ -24,10 +24,10 @@ calc_vpd <- function( patm = NA, elv = NA ) { - + ##----------------------------------------------------------------------- - ## Ref: Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure - ## Calculation Methods, in Evaporation and Evapotranspiration: + ## Ref: Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure + ## Calculation Methods, in Evaporation and Evapotranspiration: ## Measurements and Estimations, Springer, London. ## vpd = 0.611*exp[ (17.27 tc)/(tc + 237.3) ] - ea ## where: @@ -37,18 +37,18 @@ calc_vpd <- function( ## calculate atmopheric pressure (Pa) assuming standard conditions at sea level (elv=0) # if (is.na(elv) && is.na(patm) && is.na(eact)){ - # + # # warning("calc_vpd(): Either patm or elv must be provided if eact is not given.") # vpd <- NA - # + # # } else { - + # if (is.na(eact)){ patm <- ifelse(is.na(patm), - ingestr::calc_patm(elv), + calc_patm(elv), patm) # } - + ## Calculate VPD as mean of VPD based on Tmin and VPD based on Tmax if they are availble. ## Otherwise, use just tc for calculating VPD. vpd <- ifelse( @@ -62,15 +62,15 @@ calc_vpd <- function( } #' Calculate instantenous VPD from ambient conditions -#' +#' #' Follows Abtew and Meleese (2013), Ch. 5 Vapor Pressure Calculation Methods, #' in Evaporation and Evapotranspiration -#' +#' #' @param qair Air specific humidity (g g-1) #' @param eact Water vapour pressure (Pa) #' @param tc temperature, deg C #' @param patm Atmospehric pressure (Pa) -#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing +#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing #' for calculating it based on standard sea level pressure) #' #' @return instantenous VPD from ambient conditions @@ -78,15 +78,15 @@ calc_vpd <- function( calc_vpd_inst <- function( qair=NA, - eact=NA, + eact=NA, tc=NA, patm=NA, elv=NA ) { - + ##----------------------------------------------------------------------- - ## Ref: Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure - ## Calculation Methods, in Evaporation and Evapotranspiration: + ## Ref: Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure + ## Calculation Methods, in Evaporation and Evapotranspiration: ## Measurements and Estimations, Springer, London. ## vpd = 0.611*exp[ (17.27 tc)/(tc + 237.3) ] - ea ## where: @@ -101,13 +101,13 @@ calc_vpd_inst <- function( calc_eact(qair, patm), eact ) - + ## calculate saturation water vapour pressure in Pa esat <- 611.0 * exp( (17.27 * tc)/(tc + 237.3) ) - + ## calculate VPD in units of Pa - vpd <- ( esat - eact ) - + vpd <- ( esat - eact ) + ## this empirical equation may lead to negative values for VPD ## (happens very rarely). assume positive... vpd <- ifelse( @@ -115,24 +115,24 @@ calc_vpd_inst <- function( 0, vpd ) - + return( vpd ) } calc_eact <- function(qair, patm){ - + # kTo <- 288.15 # base temperature, K (Prentice, unpublished) kR <- 8.3143 # universal gas constant, J/mol/K (Allen, 1973) kMv <- 18.02 # molecular weight of water vapor, g/mol (Tsilingiris, 2008) kMa <- 28.963 # molecular weight of dry air, g/mol (Tsilingiris, 2008) - + ## calculate the mass mixing ratio of water vapor to dry air (dimensionless) wair <- qair / (1 - qair) - - ## calculate water vapor pressure + + ## calculate water vapor pressure rv <- kR / kMv rd <- kR / kMa eact <- patm * wair * rv / (rd + wair * rv) - + return(eact) -} \ No newline at end of file +} diff --git a/R/grsofun_run.R b/R/grsofun_run.R index c984d3b..5e41b77 100644 --- a/R/grsofun_run.R +++ b/R/grsofun_run.R @@ -65,7 +65,8 @@ grsofun_run_byilon <- function(ilon, par, settings){ # That's where it all happens... message(paste("Running rsofun for ilon =", ilon)) - # xxx test: for DE-Tha (DE-Tha lon = 13.6, lat = 51.0, elv = 380 m), use ilon = 388 (lon = 13.75, lat = 50.75, sitename = grid_ilon_388_ilat_174) + # xxx test: ilon = 388 + # for DE-Tha (DE-Tha lon = 13.6, lat = 51.0, elv = 380 m), use ilon = 388 (lon = 13.75, lat = 50.75, sitename = grid_ilon_388_ilat_174) # get land mask - variable name hard coded df <- readr::read_rds(paste0(settings$dir_landmask_tidy, "WFDEI-elevation_ilon_", ilon, ".rds")) |> @@ -103,6 +104,7 @@ grsofun_run_byilon <- function(ilon, par, settings){ ) |> # convert units and rename + dplyr::rowwise() |> dplyr::mutate( Tair = Tair - 273.15, # K -> deg C ppfd = SWdown * kfFEC * 1.0e-6, # W m-2 -> mol m-2 s-1 @@ -118,7 +120,6 @@ grsofun_run_byilon <- function(ilon, par, settings){ netrad = NA, ccov = 0.5, co2 = 400, - fapar = 1, ) |> dplyr::select( lon, @@ -132,7 +133,6 @@ grsofun_run_byilon <- function(ilon, par, settings){ ccov, snow = Snowf, co2, - fapar, patm = PSurf, tmin = Tair, tmax = Tair @@ -141,19 +141,114 @@ grsofun_run_byilon <- function(ilon, par, settings){ dplyr::group_by(lon, lat) |> tidyr::nest() + # # xxx test + # df_climate |> + # filter(lat == 50.75) |> + # unnest(data) |> + # ggplot(aes(date, vpd)) + + # geom_line() + # + # df_climate |> + # filter(lat == 50.75) |> + # unnest(data) |> + # visdat::vis_miss() + } if (settings$source_fapar == "modis"){ - df <- read_rds(paste0(settings$dir_fapar_tidy, "MODIS-C006_MOD15A2_LAI_FPAR_zmaw_ilon_", ilon, ".rds")) |> - tidyr::unnest(data) - # match years read from climate - dates <- df_forcing$data[[1]]$date - year_end <- max(lubridate::year(dates[lubridate::month(dates) == 12 & lubridate::mday(dates) == 31])) - year_start <- min(lubridate::year(dates[lubridate::month(dates) == 1 & lubridate::mday(dates) == 1])) + # read monthly fAPAR data + df_fapar_mon <- read_rds(paste0(settings$dir_fapar_tidy, "MODIS-C006_MOD15A2_LAI_FPAR_zmaw_ilon_", ilon, ".rds")) |> + mutate(data = purrr::map(data, ~rename(., date = time))) |> + + # xxx try + filter(lat == 50.75) + + # df_fapar_mon$data[[1]] |> + # ggplot(aes(date, fpar)) + + # geom_line() + + # linearly interpolate monthly fAPAR values to daily + dates <- df_fapar_mon$data[[1]]$date + year_start <- min(lubridate::year(dates[lubridate::month(dates) == 1])) + year_end <- max(lubridate::year(dates[lubridate::month(dates) == 12])) + + ddf <- tibble( + date = seq( + from = lubridate::ymd(paste0(year_start, "-01-01")), + to = lubridate::ymd(paste0(year_end, "-12-31")), + by = "days" + )) + + interpolate2daily_fpar <- function(df, ddf){ + # linearly interpolate (leaves training NAs) + ddf <- ddf |> + dplyr::left_join( + df, + by = "date" + ) |> + dplyr::mutate( + fpar_daily = zoo::na.approx(fpar, na.rm = FALSE) + ) + + # fill remaining with mean seasonal cycle + meandf <- ddf |> + dplyr::mutate(doy = lubridate::yday(date)) |> + dplyr::group_by(doy) |> + dplyr::summarise(fpar_meandoy = mean(fpar_daily, na.rm = TRUE)) + + ddf <- ddf |> + dplyr::mutate(doy = lubridate::yday(date)) |> + dplyr::left_join( + meandf, + by = "doy" + ) |> + dplyr::mutate(fpar_daily = ifelse(is.na(fpar_daily), fpar_meandoy, fpar_daily)) |> + dplyr::select(-fpar_meandoy, -doy) + + return(ddf) + } + + df_fapar <- df_fapar_mon |> + mutate(data = purrr::map(data, ~interpolate2daily_fpar(., ddf))) |> + mutate(data = purrr::map(data, ~dplyr::select(., -fpar))) |> + mutate(data = purrr::map(data, ~dplyr::rename(., fapar = fpar_daily))) + + # df_fapar$data[[1]] |> + # ggplot(aes(date, fpar)) + + # geom_line() + + # # reduce to match years available from climate forcing + # dates <- df_climate$data[[1]]$date + # year_start <- min(lubridate::year(dates[lubridate::month(dates) == 1 & lubridate::mday(dates) == 1])) + # year_end <- max(lubridate::year(dates[lubridate::month(dates) == 12 & lubridate::mday(dates) == 31])) + } + # combine to the most limiting time span among the climate and fapar time series + df_forcing <- df_climate |> + tidyr::unnest(data) |> + dplyr::inner_join( + df_fapar |> + tidyr::unnest(data), + by = c("lon", "lat", "date") + ) |> + dplyr::group_by(lon, lat) |> + tidyr::nest() + + # # xxx test + # df_forcing |> + # filter(lat == 50.75) |> + # unnest(data) |> + # ggplot(aes(date, vpd)) + + # geom_line() + # + # df_forcing |> + # filter(lat == 50.75) |> + # unnest(data) |> + # visdat::vis_miss() + # populate driver object if (settings$model == "pmodel"){ @@ -161,7 +256,7 @@ grsofun_run_byilon <- function(ilon, par, settings){ # merge with forcing time series left_join( - df_climate |> + df_forcing |> dplyr::rename(forcing = data), join_by(lon, lat) ) |> @@ -221,6 +316,18 @@ grsofun_run_byilon <- function(ilon, par, settings){ ) } + # xxx test (near DE-Tha) + df |> + filter(sitename == "grid_ilon_388_ilat_174") |> + unnest(forcing) |> + ggplot(aes(date, vpd)) + + geom_line() + + df |> + filter(sitename == "grid_ilon_388_ilat_174") |> + unnest(forcing) |> + visdat::vis_miss() + # run model out <- rsofun::runread_pmodel_f( drivers = df |> @@ -228,7 +335,11 @@ grsofun_run_byilon <- function(ilon, par, settings){ par = par ) - # xxx test (near DE-Tha) + # xxx test + out$data[[1]] |> + ggplot(aes(date, fapar)) + + geom_line() + out$data[[1]] |> ggplot(aes(date, gpp)) + geom_line() diff --git a/vignettes/global_rsofun_example.Rmd b/vignettes/global_rsofun_example.Rmd index c9ebf04..bd19072 100644 --- a/vignettes/global_rsofun_example.Rmd +++ b/vignettes/global_rsofun_example.Rmd @@ -13,6 +13,7 @@ vignette: > knitr::opts_chunk$set(echo = TRUE) source_files <- list.files(here::here("R/"), "*.R$") # locate all .R files purrr::map(paste0(here::here("R/"), source_files), ~source(.)) +library(tidyverse) ``` ## Specify settings