R/weather.R

In growR: Implementation of the Vegetation Model ModVege

#' 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)
    }
  )
)