R/daymet.R

In bmcnellis/RSFIA: FIA Tree Mortality Analysis

#' @include generics.R
#' @include nc_helper.R
#'
#' @title Class, methods, and functions for pulling/analyzing daymet data
#'
#' @description
#'
#' Wrapper/processing functions for daymet data.
#'
#' @details
#'
#' Process steps are as follows:
#'
#' 1) Download daymet files using the daymetr package
#'
#' 2) Store daymet data as .nc files using the ncdf4 package
#'
#' 3) Process daymet data into dataframes useful for the statisical analysis
#'
#' @author Brandon McNellis
#'
#' @name daymet
#' @rdname daymet
NULL
#' An S4 class for daymet data
#'
#' @slot site Site metadata from daymet header
#' @slot elevation Numeric elevation vector, pulled from daymet header
#' @slot tile Tile metadata from daymet headers
#'
#' @rdname daymet
daymet <- setClass(
  'daymet',
  slots = list(
    site = 'character',
    elevation = 'numeric',
    tile = 'numeric'
  ),
  contains = 'nc_helper')
#' @export
setValidity('daymet', function(object) {
  errors <- character()

  # Returns
  if (length(errors) == 0) {
    TRUE
  } else {
    errors
  }
})
#' @rdname daymet
#' @export
setMethod('initialize',
          signature(.Object = 'daymet'),
          function(.Object, ...) {
            params <- list(...)

            # time_units
            .Object@time_units <- 'years'

            # variables
            .Object@variables <- DaymetMeta()$vars

            # daymet header info:
            .Object@elevation <- as.numeric(rep(NA, length(.Object@coords)))
            .Object@site <- as.character(rep(NA, length(.Object@coords)))
            .Object@tile <- as.numeric(rep(NA, length(.Object@coords)))

            .Object <- callNextMethod()
            mt <- validObject(.Object)
            if (isTRUE(mt)) {
              return(.Object)
            } else {
              return(mt)
            }
          }
)
#' @rdname daymet
#' @docType methods
#' @export
setMethod('SetupDataFile',
          signature(object = 'daymet'),
          function(object, overwrite = F, backup = F) {
            dwd <- object@nc_dir
            fname <- paste0(object@nc_dir, '/', object@file_name)

            if (file.exists(fname)) {
              if (overwrite) {
                cat('\nOVERWRITING = T\n')
                if (backup) {
                  cat('\nBacking up .nc file...')
                  tof <- paste0(object@nc_dir, '/backup/', object@file_name)
                  file.copy(from = fname, to = tof)
                }
                file.remove(fname)
              } else {
                stop('.nc file already exists')
              }
            }
            if (!(dir.exists(dwd))) {
              dir.create(dwd)
            }

            # Create .nc object
            in_dim <- list(
              ncdf4::ncdim_def('sample', units = '', vals = object@sample, unlim = T),
              ncdf4::ncdim_def('yday', units = '', vals = seq(365), unlim = F),
              ncdf4::ncdim_def('time', units = object@time_units, vals = object@time, unlim = T)
            )

            vl <- vector('list', length(object@variables))
            for (i in seq_along(vl)) {
              ii <- object@variables[i]
              vl[[i]] <- ncdf4::ncvar_def(ii, '', in_dim, NA)
            }

            cat('\nWriting nc file...\n')
            ncdf4::nc_create(fname, vl)
            invisible(object)
          }
)
#' @rdname daymet
#' @export
setMethod("FillArray",
          signature(object = "daymet"),
          function(object, df, nc = NULL) {
            validObject(object)

            if (is.null(nc)) {
              fname <- paste0(object@nc_dir, '/', object@file_name)
              nc <- ncdf4::nc_open(fname, write = T)
            }

            t0 <- which(object@time == unique(df$time))
            tt <- ifelse(length(t0) == 1, t0, 1)
            s0 <- which(object@sample == unique(df$sample))
            ss <- ifelse(length(s0) == 1, s0, 1)
            for (j in seq(365)) {
              st <- c(ss, j, tt)
              ct <- c(length(unique(df$sample)), 1, length(unique(df$time)))
              dfj <- df[which(df$yday == j), ]

              for (i in seq(ncol(dfj))) {
                ii <- colnames(dfj)[i]
                if (ii %in% object@variables) {
                  vali <- dfj[[ii]] # assumes sample is first column
                  #cat('\rj:', j, ' ii:', ii, ' st:', st, ' ct:', ct, ' nval:', length(vali))
                  ncdf4::ncvar_put(nc = nc, varid = ii, start = st, count = ct, vals = vali)
                } else {
                  next
                }
              } # end i
            }

            invisible()
          }
)
#' Download daymet data into an nc file
#' @rdname daymet
#' @export
DownloadDaymet <- function(object) {
  stopifnot(validObject(object))

  # get slots
  time <- object@time
  cc0 <- object@coords
  lon <- sapply(cc0, function(x) x[1])
  lat <- sapply(cc0, function(x) x[2])

  fname <- paste0(object@nc_dir, '/', object@file_name)
  if (!(file.exists(fname))) {
    SetupDataFile(object)
  }
  daymet_nc <- ncdf4::nc_open(fname, write = T, readunlim = T)

  pt <- object@n_fill
  sp <- object@sample
  for (i in pt:length(cc0)) {
    ii0 <- lon[i]
    ii1 <- lat[i]
    dmet_i <- daymetr::download_daymet('Daymet', ii1, ii0, time[1], time[length(time)], silent = T)
    dmdf <- dmet_i$data
    dmdf$sample <- rep(sp[i], nrow(dmdf))
    dmdf[, ] <- dmdf[, c(1, 10, 2, 3, 4, 5, 6, 7, 8, 9)]
    colnames(dmdf) <- c('time', 'sample', 'yday', DaymetMeta()$vars)


    FillArray(object, df = dmdf, nc = daymet_nc)

    object@n_fill <- object@n_fill + 1L
    object@elevation[i] <- dmet_i$altitude
    object@site[i] <- dmet_i$site
    object@tile[i] <- dmet_i$tile

    cat('\r', format(i / length(cc0) * 100, digits = 2, nsmall = 3), '%, i:', i)
  }

  ncdf4::nc_close(daymet_nc)
  object <- UpdateTimestamp(object)
  return(object)

}
#' @rdname daymet
#' @export
DaymetMeta <- function() {
  vars <- c('dayl', 'prcp', 'srad', 'swe', 'tmax', 'tmin', 'vp')
  units <- c('s-1', 'mm day-1', 'W m-2', 'kg m-2', 'deg C', 'deg C', 'Pa')
  lnames <- c('day length', 'precipitation', 'shortwave radiation',
              'snow water equivalent', 'maximum air temperature',
              'minimum air temperature', 'water vapor pressure')
  return(data.frame(vars, units, lnames, stringsAsFactors = F))
}