Nothing
#' Weather Data Object
#'
#' @description
#' Data structure containing weather data for a given site for several years.
#'
#' @details
#' All fields representing weather variables are vectors of length 365 times N,
#' where *N* is the number of years for which weather data is stored. In
#' other words, every variable has one value for each of the 365 of each of
#' the *N* years.
#'
#' # Weather inputs
#' The weather input file should be organized as space separated columns
#' with a `year` column and at least the following parameters as headers
#' (further columns are ignored):
#'
#' ```{r child = "vignettes/children/weather_parameters.Rmd"}
#' ```
#'
#' These parameters are stored in this object in the respective `PARAM_vec`
#' fields.
#'
#' # Snow model
#'
#' The precipitation and temperature inputs are used in order to estimate the
#' snow cover for each day by use of a snow model.
#' The employed model is as formulated by Kokkonen et al. 2006 and makes use
#' of parameters from Rango and Martinec, 1995.
#'
#' @field weather_file Name of provided weather data file.
#' @field years numeric Integer representation of the contained years.
#' @field vec_size Length of the PARAM_vec vectors, which is equal to *number
#' of contained years* times 365.
#' @field year_vec Vector of length *vec_size*, holding the year for the
#' respective index.
#' @field W A list generated by `get_weather_for_year()` which contains
#' weather data only for a given year. The keys in the list are:
#' - aCO2 (atmospheric CO2 concentration in ppm)
#' - year
#' - DOY
#' - Ta
#' - Ta_sm (smoothed daily average temperature)
#' - PAR
#' - PP
#' - PET
#' - liquidP
#' - melt
#' - snow
#' - ndays (number of days in this year)
#'
#' @references
#' \insertRef{rango1995RevisitingDegreeDayMethod}{growR}
#'
#' \insertRef{kokkonen2006ConstructionDegreedaySnow}{growR}
#'
#' @md
#' @export
WeatherData = R6Class(
"WeatherData",
public = list(
#-Public-attributes-----------------------------------------------------------
weather_file = NULL,
years = NULL,
day_length_vec = NULL,
liquidP_vec = NULL,
melt_vec = NULL,
snow_vec = NULL,
year_vec = NULL,
DOY_vec = NULL,
Ta_vec = NULL,
PP_vec = NULL,
PAR_vec = NULL,
PET_vec = NULL,
vec_size = NULL,
Ta_smooth = NULL,
W = NULL,
#-Public-methods--------------------------------------------------------------
#' @description Create a new WeatherData object.
#'
#' @param weather_file string Path to file containing the weather data to
#' be read.
#' @param years numeric Vector of years for which the weather is to be
#' extracted.
#'
#' @seealso [WeatherData]`$read_weather()`
initialize = function(weather_file = NULL, years = NULL) {
if (!is.null(weather_file)) {
self$read_weather(weather_file, years)
}
},
#' @description Read weather data from supplied *weather_file*.
#'
#' @param weather_file Path to or name of file containing weather data.
#' @param years Years for which the weather is to be extracted.
#' Default (NULL) is to read in all found years.
#'
read_weather = function(weather_file, years = NULL) {
self$weather_file = weather_file
self$years = years
# Load weather data
if (file.exists(weather_file)) {
logger(sprintf("Loading weather data from %s.", weather_file),
level = DEBUG)
} else {
stop(sprintf("Weather file `%s` not found.", weather_file))
}
weather = read.table(weather_file, header = TRUE)
weather = self$ensure_file_integrity(weather)
# Only consider relevant years and omit leap days
selector = weather$DOY < 366
if (!is.null(years)) {
selector = selector & weather$year %in% years
} else {
self$years = unique(weather$year)
}
# Warn if no matching data was found
if (!any(selector)) {
message = "No matching data found for specified years (%s) in `%s`."
years_string = paste(years, collapse = ", ")
stop(sprintf(message, years_string, weather_file))
}
year_vec = weather$year[selector]
DOY_vec = weather$DOY[selector]
Ta_vec = weather$Ta[selector]
PP_vec = weather$precip[selector]
PAR_vec = weather$PAR[selector]
PET_vec = weather$ET0[selector]
vec_size = length(year_vec)
#-Tweaks,-corrections-and-further-weather-variables-----------------------
# Smooth time series of air temperature, as needed to determine the
# start of the growing season.
Ta_smooth = numeric(vec_size)
weather_smooth_window = 6
for (k in 1:vec_size) {
Ta_smooth[k] = mean(Ta_vec[max(1,(k - weather_smooth_window)):k])
}
## Estimate snow cover
# melt temperature [C] (Rango and Martinec, 1995)
T_melt = -1.
# melt coefficient [mm C-1 d-1] (Rango and Martinec, 1995)
C_melt = 3.
# freeze coefficient [mm C-1 d-1] (Kokkonen et al, 2006)
C_freeze = 0.05
# liquid water retention coeff. (Kokkonen et al, 2006)
C_retention = 0.25
liquidP_vec = (1./(1. + exp(-1.5*(Ta_vec - 2.)))) * PP_vec
solidP_vec = PP_vec - liquidP_vec
melt_vec = numeric(vec_size)
freeze_vec = numeric(vec_size)
snow_vec = numeric(vec_size)
melt_vec[1] = 0.
freeze_vec[1] = 0.
snow_vec[1] = 0.
for (j in 2:vec_size) {
# :TODO: Is time_step needed here?
time_step = 1
melt_vec[j] = ifelse(snow_vec[j - 1] > 0. & Ta_vec[j] >= T_melt,
min(snow_vec[j - 1]/time_step,
C_melt * (Ta_vec[j] - T_melt)),
0.)
freeze_vec[j] = ifelse(Ta_vec[j] < T_melt,
min(liquidP_vec[j],
C_freeze * (T_melt - Ta_vec[j])),
0.)
dsnow_dt = solidP_vec[j] + freeze_vec[j] - melt_vec[j]
snow_vec[j] = max(0., snow_vec[j - 1] + dsnow_dt * time_step)
}
self[["liquidP_vec"]] = liquidP_vec
self[["melt_vec"]] = melt_vec
self[["snow_vec"]] = snow_vec
self[["year_vec"]] = year_vec
self[["DOY_vec"]] = DOY_vec
self[["Ta_vec"]] = Ta_vec
self[["PP_vec"]] = PP_vec
self[["PAR_vec"]] = PAR_vec
self[["PET_vec"]] = PET_vec
self[["vec_size"]] = vec_size
self[["Ta_smooth"]] = Ta_smooth
},
#' @description Check if supplied input file is formatted correctly.
#'
#' Check if required column names are present and fix NA entries.
#'
#' @param weather data.table of the read input file with `header = TRUE`.
#'
ensure_file_integrity = function(weather) {
required = c("year", "DOY", "Ta", "precip", "PAR", "ET0")
data_name = sprintf("the supplied weather file `%s`", self$weather_file)
ensure_table_columns(required, weather, data_name = data_name)
# Hard fix NA values
weather[is.na(weather)] = 0.0
return(weather)
},
#' @description Calculate the expected length of day based on a site's
#' geographical latitude.
#'
#' @param latitude numeric; geographical latitude in degrees.
#'
calculate_day_length = function(latitude) {
day_length_vec = numeric(self[["vec_size"]])
# Convert to radians
angular_velocity = 2. * pi / 365.
latitude = latitude * pi / 180.
# The solar declination starts from 0 degree in the March equinox
# (March 20th), increases to its maximaum value of the Earth's tilt
# angle of 23.4 degree in the Summer solstice of June 20/21. From there
# it falls again, crossing 0 during the September equinox on September
# 22/23, before approaching its minimum value of -23.4 degree on the
# Winter solstice on December 21 or 22.
# DOY for March 20th is 31 + 28 + 20 = 79
# 23.4 degree = 0.408 RAD
# IRelDst = 1. + 0.033*cos(angular_velocity*DOY_vec)
solar_declination = 0.408 *
sin(angular_velocity * (self[["DOY_vec"]] - 79))
# Use the sunset equation to calculate the hour-angle of sunset
sunset_hourangle = acos(-tan(latitude)*tan(solar_declination))
# Convert hour-angle to hours of sunshine (factor two since sunset and
# sunrise are symmetric around noon for the hour-angle). Angular
# velocity of Earth around itself is roughly 15 degrees per hour.
# day_length_vec = 24. * sunset_hourangle / pi
day_length_vec = 2 * sunset_hourangle / (15*pi/180)
self[["day_length_vec"]] = day_length_vec
return(day_length_vec)
},
#' @description Extract state variables to the weather data for given
#' year and return them as a list.
#'
#' @param year integer Year for which to extract weather data.
#' @return W List containing the keys aCO2, year, DOY, Ta, Ta_sm, PAR,
#' PP, PET, liquidP, melt, snow, ndays.
#'
get_weather_for_year = function(year) {
iW = which(self$year_vec == year)
if (!any(iW)) {
stop(sprintf("No weather data found for year %s.", year))
}
W = list()
W[["aCO2"]] = atmospheric_CO2(year)
W[["year"]] = self$year_vec[iW]
W[["DOY"]] = self$DOY_vec[iW]
W[["Ta"]] = self$Ta_vec[iW]
W[["Ta_sm"]] = self$Ta_smooth[iW]
W[["PAR"]] = self$PAR_vec[iW]
W[["PP"]] = self$PP_vec[iW]
W[["PET"]] = self$PET_vec[iW]
W[["liquidP"]] = self$liquidP_vec[iW]
W[["melt"]] = self$melt_vec[iW]
W[["snow"]] = self$snow_vec[iW]
W[["ndays"]] = length(W[["year"]])
self[["W"]] = W
return(W)
}
)
)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.