R/calc_bioclim_raster.R

In sitscholl/rebecka_package: Utility functions used in the REBECKA project

#' Bioclimatic index calculation
#'
#' This function allows you to calculate bioclimatic indices based on the overlay of raster datasets.
#' @param tx Raster brick with daily rasters of maximum temperature. The names of the rasters have to be the same as the names in 'tn' and must correspond to the day which the values refer to.
#' @param tn Raser brick with daily rasters of minimum temperature. The names of the rasters have to be the same as the names in 'tx' and must correspond to the day which the values refer to.
#' @param index.calc Name of the index that you want to calculate. There are some predefined indices within this function (Currently: Winkler, Huglin, GSTavg). If your index is not within this list you have to supply the necessary parameters manually and can leave out this parameter.
#' @param frml A formula object that describes how to calculate the daily index values. The left side of the formula will be used to name the resulting RasterLayer and the right side to calculate the daily rasters.
#' @param agg.func A function that will be used to aggregate the daily rasters calculated by 'frml' for every year
#' @param period A sequence of numbers that indicates the day of years that should be used to aggregate the index
#' @param set.zero Are values below 0 possible or should they be converted to 0?
#' @keywords bioclimatic index
#' @export
#' @import data.table
#' @examples
#' calc_bioclim_raster(tx = tmax_brick,
#'   tn = tmin_brick,
#'   frml = winkler ~ ((T_max + T_min) / 2) - 10,
#'   agg.func = sum,
#'   period = 91:304,
#'   set.zero = T)

calc_bioclim_raster_old <- function(tx, tn, index.calc,
                                frml, agg.func, period, set.zero) {


  calc_daily_raster <- function(i) {

    T_max <- tx[[i]]
    T_min <- tn[[i]]

    index.daily <- eval(parse(text = frml.string))

    if (set.zero) {
      index.daily[index.daily < 0] <- 0
    }

    index.daily
  }
  get_index_params <- function(x) {

    index.param.list <- list('huglin' = list('frml' = huglin ~ (((T_mean - 10) + (T_max - 10)) / 2) * 1.05,
                                             'agg.func' = sum,
                                             'period' = 91:273,
                                             'set.zero' = T),
                             'winkler' = list('frml' = winkler ~ T_mean - 10,
                                              'agg.func' = sum,
                                              'period' = 91:304,
                                              'set.zero' = T),
                             'GSTavg' = list('frml' = GSTavg ~ T_mean,
                                             'agg.func' = mean,
                                             'period' = 91:304,
                                             'set.zero' = F),
                             'rad_sum' = list('frml' = rad_sum ~ SolRad,
                                              'agg.func' = sum,
                                              'period' = 91:304,
                                              'set.zero' = F))

    if (!x %in% names(index.param.list)) {
      stop('Index name has not been found in parameter list')
    }

    return(index.param.list[[x]])
  }


  #Check if both raster bricks are of the same length
  if (dim(tn)[3] != dim(tx)[3]) {
    stop('The two raster bricks are not of the same length!')
  }

  #Check if both raster bricks contain the same days
  if (!all(names(tn) == names(tx))) {
    stop('The two raster bricks have not the same raster names')
  }

  raster.names <- c(names(tn), names(tx))
  #Regex expression that finds every character between two numbers as well as
  #at the beginning and end of numbers
  #which is then used as separating string for Date generation
  sep.string <- unique(unlist(str_extract_all(raster.names, '(?<=\\d)([^0-9])(?=\\d)')))
  start.string <- unique(unlist(str_extract_all(raster.names, '^[^0-9]*')))
  end.string <- unique(unlist(str_extract_all(raster.names, '[^0-9]+$')))

  #Make sure all the names have the same strings
  if (length(sep.string) > 1 | length(start.string) > 1 | length(end.string) > 1) {
    stop('No unique string found for raster names!')
  }

  #Transform raster names to dates and extract years and day of years
  date.format <- str_c(start.string,
                       str_c('%Y', '%m', '%d', sep = sep.string),
                       end.string)
  names.date <- as.Date(names(tn), format = date.format)
  years <- as.numeric(format(names.date, '%Y'))
  days <- as.numeric(format(names.date, '%j'))

  #If index name is supplied get parameters from default list
  if (!missing(index.calc)) {
    params <- get.index.params(index.calc)
    frml <- params[['frml']]
    agg.func <- params[['agg.func']]
    period <- params[['period']]
    set.zero <- params[['set.zero']]
  }

  index.name <- all.vars(frml[[2]])
  frml.string <- deparse(frml[[3]])

  #Iterate over the years and calculate the index for every year
  #The results are collected in a stack of raster files
  result.stack <- stack()
  for (year in unique(years)) {

    #Extract the index of dates in names.date that are within period
    raster.stack.index <- which(years == year & days %in% period)

    #Calculate daily rasters of the index
    index.list <- lapply(raster.stack.index, calc_daily_raster)
    index.brick <- brick(index.list)

    #Summarize the daily rasters
    index.final <- calc(index.brick,
                        fun = agg.func)

    index.final <- projectRaster(index.final,
                                 crs="+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
    names(index.final) <- str_c(index.name, '_', year)

    result.stack <- stack(result.stack, index.final)
  }
  result.stack
}



#' Radiation sum calculation
#'
#' This function allows you to calculate the sum of solar radiation based on the overlay of raster datasets.
#' @param rad Raster brick with daily rasters of radiation energy input. The names of the rasters must correspond to the day which the values refer to.
#' @param period A sequence of numbers that indicates the day of years that should be used to aggregate the index
#' @keywords radiation, raster
#' @export
#' @import data.table
#' @examples
#' calc_rad_sum_raster(rad = solrad_brick,
#'   period = 91:304)

calc_rad_sum_raster_old <- function(rad, period = 91:304) {


  calc_daily_raster <- function(i) {

    SolRad <- rad[[i]]

    index.daily <- eval(parse(text = frml.string))

    index.daily
  }

  frml = rad_sum ~ SolRad
  agg.func = sum
  period = period

  raster.names <- names(rad)
  #Regex expression that finds every character between two numbers as well as
  #at the beginning and end of numbers
  #which is then used as separating string for Date generation
  sep.string <- unique(unlist(stringr::str_extract_all(raster.names, '(?<=\\d)([^0-9])(?=\\d)')))
  start.string <- unique(unlist(stringr::str_extract_all(raster.names, '^[^0-9]*')))
  end.string <- unique(unlist(stringr::str_extract_all(raster.names, '[^0-9]+$')))

  #Make sure all the names have the same strings
  if (length(sep.string) > 1 | length(start.string) > 1 | length(end.string) > 1) {
    stop('No unique string found for raster names!')
  }

  #Transform raster names to dates and extract years and day of years
  date.format <- stringr::str_c(start.string,
                                stringr::str_c('%Y', '%m', '%d', sep = sep.string),
                                end.string)
  names.date <- as.Date(raster.names, format = date.format)
  years <- as.numeric(format(names.date, '%Y'))
  days <- as.numeric(format(names.date, '%j'))

  index.name <- all.vars(frml[[2]])
  frml.string <- deparse(frml[[3]])

  #Iterate over the years and calculate the index for every year
  #The results are collected in a stack of raster files
  result.stack <- raster::stack()
  years.unique <- unique(years)
  total <- length(years.unique)
  for (i in seq_along(years.unique)) {

    year <- years.unique[i]

    if (i == 1) {
      pb <- txtProgressBar(min = 0, max = total, style = 3)
    }

    #Extract the index of dates in names.date that are within period
    raster.stack.index <- which(years == year & days %in% period)

    #Calculate daily rasters of the index
    index.list <- lapply(raster.stack.index, calc_daily_raster)
    index.brick <- raster::brick(index.list)

    #Summarize the daily rasters
    index.final <- raster::calc(index.brick,
                                fun = agg.func)

    index.final <- raster::projectRaster(index.final,
                                         crs="+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
    names(index.final) <- stringr::str_c(index.name, '_', year)

    result.stack <- raster::stack(result.stack, index.final)

    setTxtProgressBar(pb, i)
  }
  result.stack
}