https://github.com/cran/FedData
Tip revision: 88660a732f4b9b57d13b7505620f5fe11cc60566 authored by R. Kyle Bocinsky on 03 October 2016, 22:45:45 UTC
version 2.3.0
version 2.3.0
Tip revision: 88660a7
DAYMET_FUNCTIONS.R
#' Download and crop the 1-km DAYMET daily weather dataset.
#'
#' \code{get_daymet} returns a \code{RasterBrick} of weather data cropped to a given
#' template study area.
#'
#' @param template A Raster* or Spatial* object to serve
#' as a template for cropping.
#' @param label A character string naming the study area.
#' @param elements A character vector of elemets to extract.\cr
#' The available elements are:\cr
#' dayl = Duration of the daylight period in seconds per day. This calculation is based on the period of the day during which the sun is above a hypothetical flat horizon.\cr
#' prcp = Daily total precipitation in millimeters per day, sum of all forms converted to water-equivalent. Precipitation occurrence on any given day may be ascertained.\cr
#' srad = Incident shortwave radiation flux density in watts per square meter, taken as an average over the daylight period of the day. NOTE: Daily total radiation (MJ/m2/day) can be calculated as follows: ((srad (W/m2) * dayl (s/day)) / l,000,000)\cr
#' swe = Snow water equivalent in kilograms per square meter. The amount of water contained within the snowpack.\cr
#' tmax = Daily maximum 2-meter air temperature in degrees Celsius.\cr
#' tmin = Daily minimum 2-meter air temperature in degrees Celsius.\cr
#' vp = Water vapor pressure in pascals. Daily average partial pressure of water vapor.\cr
#' @param years A numeric vector of years to extract.
#' @param raw.dir A character string indicating where raw downloaded files should be put.
#' The directory will be created if missing. Defaults to "./RAW/DAYMET/".
#' @param extraction.dir A character string indicating where the extracted and cropped DEM should be put.
#' The directory will be created if missing. Defaults to "./EXTRACTIONS/DAYMET/".
#' @param force.redo If an extraction for this template and label already exists, should a new one be created?
#' @return A named list of \code{RasterBrick}s of weather data cropped to the extent of the template.
#' @importFrom lubridate year
#' @importFrom magrittr %>% %$%
#' @importFrom readr read_rds write_rds
#' @importFrom raster brick projection mosaic writeRaster
#' @importFrom sp spTransform %over%
#' @importFrom foreach foreach %do% %dopar%
#' @importFrom utils data
#' @export
#' @examples
#' \dontrun{
#' # Extract data for the Village Ecodynamics Project "VEPIIN" study area:
#' # http://village.anth.wsu.edu
#' vepPolygon <- polygon_from_extent(raster::extent(672800,740000,4102000,4170000),
#' proj4string="+proj=utm +datum=NAD83 +zone=12")
#'
#' # Get the DAYMET (North America only)
#' # Returns a list of raster bricks
#' DAYMET <- get_daymet(template=vepPolygon,
#' label="VEPIIN",
#' elements = c("prcp","tmin","tmax"),
#' years = 1980:1985)
#'
#' # Plot with raster::plot
#' plot(DAYMET$tmin$X1985.10.23)
#' }
get_daymet <- function(template, label, elements=NULL, years = NULL, raw.dir="./RAW/DAYMET/", extraction.dir=paste0("./EXTRACTIONS/",label,"/DAYMET/"), force.redo=F){
dir.create(raw.dir, showWarnings = FALSE, recursive = TRUE)
dir.create(extraction.dir, showWarnings = FALSE, recursive = TRUE)
all.elements <- c("dayl", "prcp", "srad", "swe", "tmax", "tmin", "vp")
if(is.null(elements))
elements <- all.elements
elements <- tolower(elements)
missing.elements <- setdiff(elements,all.elements)
if(length(missing.elements)>0) warning("Elements not available: ",paste(missing.elements,collapse = ", "))
elements <- setdiff(elements,missing.elements)
if(length(elements) == 0)
stop("No elements available")
all.years <- 1980:(lubridate::year(Sys.time())-1)
if(is.null(years))
years <- all.years
missing.years <- setdiff(years,all.years)
if(length(missing.years)>0) warning("Years not available: ",paste(missing.years,collapse = ", "))
years <- setdiff(years,missing.years)
if(length(years) == 0)
stop("No years available")
out.files <- paste0(extraction.dir,"/",label,"_DAYMET_",elements,"_",min(years),"-",max(years),".tif")
if(!force.redo & all(file.exists(out.files)) & file.exists(paste0(extraction.dir,"/",label,"_DAYMET_layer_names.Rds"))){
extracted.DAYMET <- out.files %>% lapply(FUN = function(x){
out <- raster::brick(x)
names(out) <- readr::read_rds(paste0(extraction.dir,"/",label,"_DAYMET_layer_names.Rds"))
return(out)
})
names(extracted.DAYMET) <- elements
return(extracted.DAYMET)
}
data(tiles, envir = environment())
template.latlon <- template %>%
sp::spTransform(raster::projection(tiles))
tile.ids <- (template.latlon %over% tiles)$TileID
message("Area of interest includes ",length(tile.ids)," DAYMET tile(s).")
# Download and crop tiles
tiles <- lapply(tile.ids,function(tile){
return(FedData::get_daymet_tile(template=template, elements = elements, years = years, tileID = tile, raw.dir=raw.dir))
})
names(tiles) <- tile.ids
# Mosaic all tiles
if(length(tiles)>1){
message('Mosaicking DAYMET tiles.')
foreach::foreach(element=elements) %do% {
utils::flush.console()
these.tiles <- lapply(tiles,"[[",element)
these.tiles$fun <- mean
these.tiles <- do.call(raster::mosaic, these.tiles)
gc()
}
}else{
tiles <- tiles[[1]]
}
tiles %>%
mapply(x = ., y = names(tiles), FUN = function(x,y){
raster::writeRaster(x, paste0(extraction.dir,"/",label,"_DAYMET_",y,"_",min(years),"-",max(years),".tif"), datatype="FLT4S", options=c("COMPRESS=DEFLATE", "ZLEVEL=9", "INTERLEAVE=BAND"), overwrite=T, setStatistics=FALSE)
})
tiles[[1]] %>%
names() %>%
readr::write_rds(paste0(extraction.dir,"/",label,"_DAYMET_layer_names.Rds"))
return(tiles)
}
#' Download a netcdf tile from the 1-km DAYMET daily weather dataset.
#'
#' Tiles are specified by a tile ID, year, and element;
#' Tiles are downloaded in the NetCDF format. \code{download_daymet_tile} returns the path to the downloaded NetCDF tile files.
#'
#' @param tileID A numeric indicating the DAYMET tile ID number.
#' @param elements A character vector of elemets to extract.\cr
#' The available elements are:\cr
#' dayl = Duration of the daylight period in seconds per day. This calculation is based on the period of the day during which the sun is above a hypothetical flat horizon.\cr
#' prcp = Daily total precipitation in millimeters per day, sum of all forms converted to water-equivalent. Precipitation occurrence on any given day may be ascertained.\cr
#' srad = Incident shortwave radiation flux density in watts per square meter, taken as an average over the daylight period of the day. NOTE: Daily total radiation (MJ/m2/day) can be calculated as follows: ((srad (W/m2) * dayl (s/day)) / l,000,000)\cr
#' swe = Snow water equivalent in kilograms per square meter. The amount of water contained within the snowpack.\cr
#' tmax = Daily maximum 2-meter air temperature in degrees Celsius.\cr
#' tmin = Daily minimum 2-meter air temperature in degrees Celsius.\cr
#' vp = Water vapor pressure in pascals. Daily average partial pressure of water vapor.\cr
#' @param years A numeric vector of years to extract.
#' @param raw.dir A character string indicating where raw downloaded files should be put.
#' The directory will be created if missing. Defaults to "./RAW/DAYMET/".
#' @return A named list of character vectors, each representing the full local paths of the tile downloads.
#' @export
#' @keywords internal
#' @importFrom doParallel registerDoParallel stopImplicitCluster
#' @importFrom foreach foreach %dopar% %do% %:%
download_daymet_tile <- function(tileID, elements, years, raw.dir){
doParallel::registerDoParallel()
out <- foreach::foreach(element = elements) %:%
foreach::foreach(year = sort(years), .combine='c') %dopar% {
destdir <- paste0(raw.dir,'/',tileID,'/',year)
dir.create(destdir, recursive = TRUE, showWarnings = FALSE)
url <- paste0('http://thredds.daac.ornl.gov/thredds/fileServer/ornldaac/1328/tiles/',year,'/',tileID,'_',year,'/',element,'.nc')
download_data(url=url, destdir=destdir, timestamping = FALSE, nc = TRUE)
paste0(destdir,'/',basename(url)) %>% normalizePath(mustWork = T)
}
doParallel::stopImplicitCluster()
names(out) <- elements
return(out)
}
#' Download and crop a netcdf tile from the 1-km DAYMET daily weather dataset.
#'
#' \code{get_daymet_tile} returns a list of \code{RasterBrick}s---one for each element---cropped within the specified \code{template}.
#' If template is not provided, returns the entire DAYMET tile.
#'
#' @param template A Raster* or Spatial* object to serve
#' as a template for cropping. If missing, entire tile is returned.
#' @param tileID A numeric indicating the DAYMET tile ID number.
#' @param elements A character vector of elemets to extract.\cr
#' The available elements are:\cr
#' dayl = Duration of the daylight period in seconds per day. This calculation is based on the period of the day during which the sun is above a hypothetical flat horizon.\cr
#' prcp = Daily total precipitation in millimeters per day, sum of all forms converted to water-equivalent. Precipitation occurrence on any given day may be ascertained.\cr
#' srad = Incident shortwave radiation flux density in watts per square meter, taken as an average over the daylight period of the day. NOTE: Daily total radiation (MJ/m2/day) can be calculated as follows: ((srad (W/m2) * dayl (s/day)) / l,000,000)\cr
#' swe = Snow water equivalent in kilograms per square meter. The amount of water contained within the snowpack.\cr
#' tmax = Daily maximum 2-meter air temperature in degrees Celsius.\cr
#' tmin = Daily minimum 2-meter air temperature in degrees Celsius.\cr
#' vp = Water vapor pressure in pascals. Daily average partial pressure of water vapor.\cr
#' @param years A numeric vector of years to extract.
#' @param raw.dir A character string indicating where raw downloaded files should be put.
#' The directory will be created if missing. Defaults to "./RAW/NED/".
#' @return A a list of \code{RasterBrick}s---one for each element---cropped within the specified \code{template}.
#' @export
#' @keywords internal
#' @importFrom doParallel registerDoParallel stopImplicitCluster
#' @importFrom foreach foreach %do% %dopar%
#' @importFrom raster projection crop stack
#' @importFrom sp spTransform CRS
#' @importFrom lubridate year
#' @importFrom magrittr %<>% %>%
get_daymet_tile <- function(template, tileID, elements = NULL, years = NULL, raw.dir){
tmpdir <- tempfile()
if (!dir.create(tmpdir))
stop("failed to create my temporary directory")
message("(Down)Loading DAYMET tile ",tileID)
all.elements <- c("dayl", "prcp", "srad", "swe", "tmax", "tmin", "vp")
elements <- tolower(elements)
if(is.null(elements))
elements <- all.elements
missing.elements <- setdiff(elements,all.elements)
if(length(missing.elements)>0) warning("Elements not available: ",paste(missing.elements,collapse = ", "))
elements <- setdiff(elements,missing.elements)
if(length(elements) == 0)
stop("No elements available")
all.years <- 1980:(lubridate::year(Sys.time())-1)
if(is.null(years))
years <- all.years
missing.years <- setdiff(years,all.years)
if(length(missing.years)>0) warning("Years not available: ",paste(missing.years,collapse = ", "))
years <- setdiff(years,missing.years)
if(length(years) == 0)
stop("No years available")
files <- download_daymet_tile(tileID = tileID, elements = elements, years = years, raw.dir=raw.dir)
doParallel::registerDoParallel()
tiles <- foreach::foreach(element=files) %dopar% {
tile <- foreach::foreach(file=element) %do%
raster::brick(file)
tile %<>% raster::stack(quick = TRUE)
if(!is.null(template)){
tile <- tryCatch(tile %>%
raster::crop(template %>%
sp::spTransform(tile %>%
raster::projection() %>%
sp::CRS()), snap="out"),
error=function(e){
tile %>%
raster::crop(template %>%
sp::spTransform(tile %>%
raster::projection() %>%
sp::CRS()))
})
}
}
doParallel::stopImplicitCluster()
names(tiles) <- elements
unlink(tmpdir, recursive = TRUE)
return(tiles)
}
#' The DAYMET tiles SpatialPolygonsDataFrame.
#'
#' A dataset containing the DAYMET tiles.
#'
#' @format A SpatialPolygonsDataFrame with 1060 features and 5 variables:
#' \describe{
#' \item{TileID}{the numeric identifier of the tile}
#' \item{XMin}{the minimum longitude of the tile}
#' \item{XMax}{the maximum longitude of the tile}
#' \item{YMin}{the minimum latitude of the tile}
#' \item{YMax}{the maximum latitude of the tile}
#' }
#' @source \url{https://github.com/khufkens/daymetr/blob/master/data/DAYMET_grid.RData}
"tiles"