Nothing
#' Get data directory
#'
#' Find out where the environmental data is stored
#'
#' @usage get_datadir(datadir)
#'
#' @param datadir character. Directory as passed to \code{\link{load_layers}}.
#' This does not change the data directory used by \code{\link{load_layers}}
#' but only shows the result of passing a specific datadir. If \code{NULL} is
#' passed then the \code{sdmpredictors_datadir} option is read. To set this
#' run \code{options(sdmpredictors_datadir = "<your preferred directory>")} in
#' every session or in your .RProfile.
#'
#' @return Path to the data directory.
#' @seealso \code{\link{load_layers}}
#' @keywords internal
get_datadir <- function(datadir) {
if(is.null(datadir)) {
datadir <- getOption("sdmpredictors_datadir")
if(is.null(datadir)) {
datadir <- file.path(tempdir(), "sdmpredictors")
warning("file.path(tempdir(), \"sdmpredictors\") will be used as datadir, set options(sdmpredictors_datadir=\"<directory>\") to avoid re-downloading the data in every session or set the datadir parameter in load_layers")
} else {
datadir <- normalizePath(paste0(datadir), winslash = "/", mustWork = NA)
}
}
if(!dir.exists(datadir)) {
dir.create(datadir, recursive = TRUE)
}
datadir
}
#' Project and saves a layer from WGS84 (EPGS:4326) to Behrmann Equal Area (ESRI:54017) coordinate reference system
#'
#' @usage equalarea_project(path)
#'
#' @param path the path to the WGS84 layer
#'
#' @return the path to the Behrmann Equal Area layer
#' @seealso \code{\link{load_layers}}
#' @keywords internal
equalarea_project <- function(path){
if(grepl(".zip", path, fixed = TRUE)){
vsizip <- "/vsizip/"
out_path <- gsub(".zip", ".tif", path, fixed = TRUE)
}else{
vsizip <- ""
out_path <- gsub("_lonlat.tif", ".tif", path, fixed = TRUE)
}
stopifnot(out_path != path)
if(!file.exists(out_path)){
r <- raster::raster(paste0(vsizip, path))
message(paste0("Projecting ", path, " from native WGS84 to Behrmann Equal Areas. This might take a few minutes."))
if(grepl("/FW_", path, ignore.case = FALSE, fixed = TRUE)){res = 815}else{res = 7000} # Quickfix. Ideally: get res from .data
out <- raster::projectRaster(r, crs = sdmpredictors::equalareaproj, method = "ngb", res = res)
raster::writeRaster(out, out_path)
}
stopifnot(as.character(raster(out_path)@crs) == as.character(sdmpredictors::equalareaproj))
return(out_path)
}
#' Load layers
#'
#' Method to load rasters from disk or from the internet. By default a
#' RasterStack is returned but this is only possible When all rasters have the
#' same spatial extent and resolution.
#'
#' @usage load_layers(layercodes, equalarea=FALSE, rasterstack=TRUE,
#' datadir=NULL)
#'
#' @param layercodes character vector or dataframe. Layer_codes of the layers to
#' be loaded or dataframe with a "layer_code" column.
#' @param equalarea logical. If \code{TRUE} then layers are loaded with a
#' Behrmann cylindrical equal-area projection (\code{\link{equalareaproj}}),
#' otherwise unprojected (\code{\link{lonlatproj}}). Default is \code{FALSE}.
#' @param rasterstack logical. If \code{TRUE} (default value) then the result is
#' a \code{\link[raster]{stack}} otherwise a list of rasters is returned.
#' @param datadir character. Directory where you want to store the data. If
#' \code{NULL} is passed (default) then the \code{sdmpredictors_datadir}
#' option is read. To set this run \code{options(sdmpredictors_datadir="<your
#' preferred directory>")} in every session or add it to your .RProfile.
#'
#' @return RasterStack or list of raster
#' @examples \dontrun{
#' # warning using tempdir() implies that data will be downloaded again in the
#' # next R session
#' env <- load_layers("BO_calcite", datadir = tempdir())
#' }
#' @export
#' @seealso \code{\link{list_layers}}, \code{\link{layer_stats}},
#' \code{\link{layers_correlation}}
load_layers <- function(layercodes, equalarea = FALSE, rasterstack = TRUE, datadir = NULL) {
if(is.na(equalarea) || !is.logical(equalarea) || length(equalarea) != 1) {
stop("equalarea should be TRUE or FALSE")
}
if(is.data.frame(layercodes)) {
layercodes <- layercodes$layer_code
}
info <- get_layers_info(layercodes)
counts <- table(info$common$time)
if(length(unique(info$common$cellsize_equalarea)) > 1) {
stop("Loading layers with different cellsize is not supported")
} else if (sum(counts) != NROW(layercodes)) {
layers <- info$common$layer_code
stop(paste0("Following layer codes where not recognized: ",
paste0(layercodes[!(layercodes %in% layers)], collapse = ", ")))
}
if(max(counts) != NROW(layercodes)) {
warning("Layers from different eras (current, future, paleo) are being loaded together")
}
if(gdal_is_lower_than_3()){
warning("GDAL is lower than version 3. Consider updating GDAL to avoid errors.")
}
datadir <- get_datadir(datadir)
get_layerpath <- function(layercode) {
layer_url <- subset(info$common, info$common$layer_code == layercode)$layer_url
if(grepl(".zip", layer_url, ignore.case = FALSE, fixed = TRUE)){
ext <- "_lonlat.zip"
}else{ext <- "_lonlat.tif"}
path <- paste0(datadir, "/", layercode, ext)
if(!file.exists(path)) {
ok <- -1
# clean up of download failed
on.exit({
if(ok != 0 && file.exists(path)) {
file.remove(path)
}
})
ok <- utils::download.file(layer_url, path, method = "auto", quiet = FALSE, mode = "wb")
}
ifelse(file.exists(path), path, NA)
}
paths <- sapply(layercodes, get_layerpath)
if(equalarea){
customSupressWarning <- function(w){if(any(grepl( "point", w))){invokeRestart( "muffleWarning")}}
paths <- withCallingHandlers(sapply(paths, equalarea_project), warning = customSupressWarning)
}
if(rasterstack) {
logical_zip <- grepl(".zip", paths, fixed = TRUE)
if(all(logical_zip) | !any(logical_zip)){
if(all(logical_zip)){
vsizip <- "/vsizip/"
}else if(!any(logical_zip)){
vsizip <- ""
}
st <- raster::stack(paste0(vsizip, paths))
if("layer" %in% names(st)){
names(st) <- layercodes
}else(
names(st) <- sub("_lonlat$", "", names(st))
)
return(st)
}else{
stop("Rasterstack for zipped and non-zipped files not supported. Try `rasterstack = FALSE`")
}
} else {
return(lapply(paths, function(path) {
if(grepl(".zip", path, fixed = TRUE)){
vsizip <- "/vsizip/"
}else{
vsizip <- ""
}
r <- raster::raster(paste0(vsizip, path))
if("layer" %in% names(r)){
names(r) <- layercodes
}else(
names(r) <- sub("_lonlat$", "", names(r))
)
r}))
}
}
#' Longitude/latitude coordinate reference system (EPSG:4326), used when using
#' load_layers with equal_area = FALSE
#' @export
lonlatproj <- sp::CRS("+proj=longlat +datum=WGS84 +no_defs")
#' World Behrmann equal area coordinate reference system (ESRI:54017), used when
#' using load_layers with equal_area = TRUE
#' @export
equalareaproj <- sp::CRS("+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
#' Is gdal v3 or more?
#' @noRd
gdal_is_lower_than_3 <- function(){
vgdal <- terra::gdal()
vgdal <- strsplit(vgdal, ".", fixed = TRUE)[[1]][1]
vgdal <- as.numeric(vgdal)
is_less_than_3 <- vgdal < 3
is_less_than_3
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.