R/extractEDEN.R

In troyhill/EvergladesEBM: Everglades Ecosystem-Based Management

#' extractEDEN
#'
#' @description This function extracts EDEN data for spatial locations of interest. This function can take a very long time depending on the number of regions-of-interest and their size, and if EDEN data is downloaded within the function call. NOTE: This is based on availability of data on the home page of EDEN. If data aren't there, they aren't accessible for this code.
#' 
#' @param targetLocations      locations to extract data from; should be class SpatVector or sf
#' @param targetLocationNames  option to specify the name of target locations (e.g., pts$gage)
#' @param EDEN_data            The name of the list of EDEN data generated from, e.g., lapply(dateRange, fireHydro::getEDEN, returnType = 'raster') or a rasterStack as generated by fireHydro::getAnnualEden(2020). Specifically, object is a list of lists with (1) a date (YYYYMMDD) and (2) EDEN data (must be a rasterStack; in `getEDEN()` use argument `returnType = 'raster'`). 
#'
#' @return a dataframe formatted identically to the traceDataLong output from extractSimData()
#'
#' @importFrom terra extract
#' @importFrom terra crs
#' @importFrom terra rast
#' @importFrom terra vect
#' @importFrom terra project
#' @importFrom terra plot
#' 
#' @examples  
#' 
#' \dontrun{
#' ### requires internet connection to download EDEN data
#' 
#' ### set regions of interest
#' locs   <- IRMap[[2]]
#' loc    <- locs[locs$INDICATOR %in% 118:119, ]
#' 
#' dateVec <- format(seq.Date(from = as.Date("20201101", format = "%Y%m%d"), 
#'             to = as.Date("202012017", format = "%Y%m%d"), by = "day"), format = "%Y%m%d")
#'             
#' ### pull and process EDEN data so it's retained in the global environment
#' EDEN_data  <- lapply(dateVec, fireHydro::getEDEN, returnType = 'raster')
#' input_eden <- list(date = do.call(rbind, lapply(EDEN_data, function(x) {x[[1]]})),
#'                 data = do.call(stack, lapply(EDEN_data, function(x) {x[[2]]})))
#' 
#' EDEN_by_IR <- extractEDEN(targetLocations     = loc,
#'                            targetLocationNames = loc$NAME,
#'                            EDEN_data           = input_eden)
#'               
#' }
#' 
#' @export
#' 



extractEDEN  <- function(targetLocations,
                         targetLocationNames,
                         EDEN_data # the name of the list of EDEN data generated from lapply(dateRange, fireHydro::getEDEN).
) {
  # if (length(needEDEN) == 1 && typeof(needEDEN) == logical && needEDEN) {
  #   test <- lapply(dateRange, fireHydro::getEDEN)
  # } else if(is.list(needEDEN)) {
  #   test <- needEDEN
  # } else {
  #   stop("needEDEN argument not recognized. It should be either 'TRUE' (signaling that EDEN data should be pulled) or a list of EDEN data.")
  # }
  test <- EDEN_data
  ### convert targetLocations to sf if necessary
  if(any(class(targetLocations) %in% c("SpatialPolygonsDataFrame", "SpatialPointsDataFrame", 'sf'))) {
    targetLocations <- terra::vect(targetLocations)
  }
  ### convert EDEN data to terra, if necessary
  if(any(class(test$data) %in% c("RasterBrick", "RasterStack", "Raster", "raster"))) {
    ### if test is a raster:
    test$data <- terra::rast(test$data)
  } else if (!any(class(test$data) %in% c("SpatRaster"))) {
    stop("EDEN data must be in terra (SpatRaster) or raster format. Pro-tip: use `returnType = 'raster'` argument in fireHydro::getEDEN()")
  }
    ### check/transform CRS
    if (!identical(terra::crs(targetLocations, proj = TRUE), terra::crs(test$data, proj = TRUE))) {
      targetLocations <- terra::project(targetLocations, terra::crs(test$data, proj = TRUE))
      terra::plot(test$data[[1]], axes = FALSE, legend = FALSE, main = 'reprojected targetLocations')
      terra::plot(targetLocations, add = TRUE)
      cat("Coordinate reference systems did not match. Conversion has been performed internally but may be incorrect (and result in NAs). \n")
    }
    
    ### pull dates
    dateRange <- format(seq.Date(from = as.Date(min(test$date, na.rm = TRUE), format = "%Y%m%d"),
                               to = as.Date(max(test$date, na.rm = TRUE), format = "%Y%m%d"), by = "day"), format = "%Y%m%d")
    test$data <- subset(test$data, match(dateRange, format(x = test$date, format = "%Y%m%d")))
    test$date <- test$date[match(dateRange, format(x = test$date, format = "%Y%m%d"))]
    # calculate mean for each polygon
    # warning appears with NAs and stops function: "In .getRat(x, ratvalues, ratnames, rattypes) : NAs introduced by coercion"
    withCallingHandlers({
      r.vals <- terra::extract(x = test$data, y = targetLocations, fun = mean, na.rm = TRUE)
    }, warning=function(w) {
      if (startsWith(conditionMessage(w), "NAs introduced by coercion"))
        invokeRestart("muffleWarning")
    })
    
    ###
    r.vals2         <- data.frame(t(r.vals))
    r.vals2         <- r.vals2[-1, ]
    if (length(targetLocations) == 1) {
      r.vals2 <- data.frame(
        mean    = r.vals2,
        cluster = rep("observed", times = length(r.vals2)))
      names(r.vals2)[1] <- targetLocationNames
    } else {
      names(r.vals2)  <- targetLocationNames
      r.vals2$cluster <- "observed"
    }
    r.vals2$date    <- dateRange[1:nrow(r.vals2)]
    
    ### convert long to wide
    dat <- stats::reshape(r.vals2, direction = "long", 
                          varying = list(names(r.vals2)[1:length(targetLocationNames)]), 
                          v.names = "ave", 
                          idvar = c("cluster", "date"), 
                          timevar = "name", 
                          times = targetLocationNames)
    dat$stdev   <- 0

    
    ### sf usage removed 20201217
  #   if (any(class(test$data) %in% "sf")) {
  #   for (i in 1:length(targetLocations)) {
  #     IR.tmp <- sapply(test, function(x) {
  #       x.sf  <- x$data
  #       a.sub <- sf::st_intersection(x.sf, sf::st_set_crs(sf::st_as_sf(targetLocations[i, ]), sf::st_crs(x.sf)))
  #       mean(a.sub$WaterDepth, na.rm = TRUE)
  #     })
  #     dat.tmp <- data.frame(cluster  = "observed",
  #                           name     = targetLocationNames[i],
  #                           date     = as.Date(as.character(sapply(test, "[[", 1)), format = "%Y%m%d"),
  #                           ave      = IR.tmp, 
  #                           stdev    = 0)
  #     if (i == 1) {
  #       dat <- dat.tmp
  #     } else {
  #       dat <- rbind(dat, dat.tmp)
  #     }
  #   }
  # }
  
  dat$date <- as.Date(dat$date, format = "%Y%m%d")
  
  invisible(dat)
}