diff --git a/.gitignore b/.gitignore index 9222578..aad7bcb 100644 --- a/.gitignore +++ b/.gitignore @@ -15,3 +15,5 @@ packrat/ .Trashes ehthumbs.db Thumbs.db +packrat/lib*/ +.lintr diff --git a/DESCRIPTION b/DESCRIPTION index d0c3b9f..672ef31 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: KrigR Type: Package Title: Downloading, Aggregating, and Kriging of ECMWF CDS-Data -Version: 0.3.0 +Version: 0.9.0 Authors@R: as.person(c( "Erik Kusch [aut, cre]", "Richard Davy [aut]" @@ -9,13 +9,13 @@ Authors@R: as.person(c( Description: An R Package for downloading, preprocessing, and statistical downscaling of data provided by the European Centre for Medium‐Range Weather Forecasts (ECMWF). KrigR contains functions for: - Downloading ECMWF data directly from within R - - Downloading USGS GMTED 2010 elevation data - Working towards also implementing support for HWSD data + - Downloading USGS GMTED 2010 and HWSD elevation data - Preparing covariate data for Kriging - Kriging spatial input to desired output using user-specified covariates License: MIT + file LICENSE Encoding: UTF-8 LazyData: true -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Imports: ecmwfr, httr, @@ -29,5 +29,9 @@ Imports: tools, progress, doSNOW, - pbapply + pbapply, + tidyr, + ggplot2, + viridis, + cowplot Depends: R (>= 4.0.0) diff --git a/LICENSE b/LICENSE deleted file mode 100644 index ddf5d79..0000000 --- a/LICENSE +++ /dev/null @@ -1,2 +0,0 @@ -YEAR: 2020 -COPYRIGHT HOLDER: Erik Kusch diff --git a/NAMESPACE b/NAMESPACE index ca020df..694aa8a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,32 +15,48 @@ export(Meta.List) export(Meta.QuickFacts) export(Meta.Read) export(Meta.Variables) -export(SummarizeRaster) -export(buffer_Points) +export(Plot.BioClim) +export(Plot.Covariates) +export(Plot.Kriged) +export(Plot.SpatRast) export(download_DEM) export(download_ERA) export(krigR) -export(mask_Shape) importFrom(automap,autoKrige) +importFrom(cowplot,plot_grid) importFrom(doSNOW,registerDoSNOW) +importFrom(ecmwfr,wf_delete) importFrom(ecmwfr,wf_get_key) importFrom(ecmwfr,wf_request) importFrom(ecmwfr,wf_set_key) importFrom(ecmwfr,wf_transfer) importFrom(foreach,"%dopar%") importFrom(foreach,foreach) +importFrom(ggplot2,aes) +importFrom(ggplot2,facet_wrap) +importFrom(ggplot2,geom_raster) +importFrom(ggplot2,geom_sf) +importFrom(ggplot2,ggplot) +importFrom(ggplot2,labs) +importFrom(ggplot2,scale_fill_gradientn) +importFrom(ggplot2,theme) +importFrom(ggplot2,theme_bw) +importFrom(ggplot2,unit) importFrom(httr,DELETE) importFrom(httr,GET) importFrom(httr,add_headers) importFrom(httr,authenticate) importFrom(httr,progress) importFrom(httr,write_disk) +importFrom(lubridate,date) importFrom(lubridate,days_in_month) +importFrom(lubridate,year) importFrom(methods,getClass) importFrom(ncdf4,nc_close) importFrom(ncdf4,nc_open) importFrom(ncdf4,ncatt_get) importFrom(ncdf4,ncatt_put) +importFrom(parallel,detectCores) importFrom(parallel,makeCluster) importFrom(parallel,stopCluster) importFrom(pbapply,pblapply) @@ -59,11 +75,13 @@ importFrom(terra,crop) importFrom(terra,crs) importFrom(terra,ext) importFrom(terra,mask) +importFrom(terra,mean) importFrom(terra,metags) importFrom(terra,nlyr) importFrom(terra,rast) importFrom(terra,res) importFrom(terra,resample) +importFrom(terra,sources) importFrom(terra,subset) importFrom(terra,tapp) importFrom(terra,time) @@ -72,4 +90,9 @@ importFrom(terra,values) importFrom(terra,varnames) importFrom(terra,writeCDF) importFrom(terra,writeRaster) +importFrom(tidyr,gather) importFrom(tools,file_path_sans_ext) +importFrom(viridis,cividis) +importFrom(viridis,inferno) +importFrom(viridis,mako) +importFrom(viridis,viridis) diff --git a/R/BioClim.R b/R/BioClim.R index 0e4212d..202091a 100644 --- a/R/BioClim.R +++ b/R/BioClim.R @@ -2,334 +2,440 @@ #' #' This function queries download of required essential climate variables from the [Climate Data Store](https://cds.climate.copernicus.eu/#!/home) hosted by the [Copernicus Climate Change Service (C3S)](https://cds.climate.copernicus.eu/about-c3s) for retrieval of climate data and subsequent calculation of bioclimatic variables for user-defined regions and time-frames. #' -#' @param Water_Var ERA5(Land)-contained climate variable targeting water availability information. Recommended values: "volumetric_soil_water_layer_1", "total_precipitation". -#' @param DataSet Which ERA5 data set to download data from. 'era5' or 'era5-land'. +#' @param Temperature_Var CDS-contained climate variable targeting temperature information. Recommended values: "2m_temperature". +#' @param Temperature_DataSet Character. Which dataset to query data from. See currently supported datasets by calling \code{\link{Meta.List}}. For now, this function is conceptualised to support "reanalysis-era5-land". +#' @param Temperature_Type Either NA or Character. Which kind of sub-type to query per data set. See \code{\link{Meta.QucikFacts}} for options per dataset. +#' @param Water_Var CDS-contained climate variable targeting water availability information. Recommended values: "volumetric_soil_water_layer_X" (where X is an integer of either 1, 2, 3, 4), "total_precipitation". +#' @param Water_DataSet Character. Which dataset to query water availability data from. See currently supported datasets by calling \code{\link{Meta.List}}. For now, this function is conceptualised to support "reanalysis-era5-land". +#' @param Water_Type Either NA or Character. Which kind of sub-type to query per water availability data set. See \code{\link{Meta.QucikFacts}} for options per dataset. #' @param Y_start Year ('YYYY') at which to start time series of downloaded data. #' @param Y_end Year ('YYYY') at which to stop time series of downloaded data. -#' @param T_res Temporal resolution from which to obtain minimum and maximum values of temperature. 'hour' or 'day' -#' @param Extent Optional, download data according to rectangular bounding box. specify as extent() object or as a raster, a SpatialPolygonsDataFrame object, or a data.frame opbject. If Extent is a SpatialPolygonsDataFrame, this will be treated as a shapefile and the output will be cropped and masked to this shapefile. If Extent is a data.frame of geo-referenced point records, it needs to contain Lat and Lon columns as well as a non-repeating ID-column. -#' @param Buffer Optional. Identifies how big a rectangular buffer to draw around points if Extent is a data frame of points. Buffer is expressed as centessimal degrees. -#' @param ID Optional. Identifies which column in Extent to use for creation of individual buffers if Extent is a data.frame. -#' @param Dir Directory specifying where to download data to. -#' @param FileName A file name for the netcdf produced. +#' @param TZone Character. Time zone in which to represent and evaluate time dimension of data. See the output of OlsonNames() for a full overview of supported specifications. Default is UTC. +#' @param Extent Optional, download data according to desired spatial specification. If missing/unspecified, total area of queried data set is used. Can be specified either as a raster object, an sf object, a terra object, or a data.frame. If Extent is a raster or terra object, data will be queried according to rectangular extent thereof. If Extent is an sf (MULTI-)POLYGON object, this will be treated as a shapefile and the output will be cropped and masked to this shapefile. If Extent is a data.frame of geo-referenced point records, it needs to contain Lat and Lon columns around which a buffered shapefile will be created using the Buffer argument. +#' @param Buffer Optional, Numeric. Identifies how big a circular buffer to draw around points if Extent is a data.frame of points. Buffer is expressed as centessimal degrees. +#' @param Dir Character/Directory Pointer. Directory specifying where to download data to. +#' @param FileName Character. A file name for the produced file. +#' @param FileExtension Character. A file extension for the produced file. Supported values are ".nc" (default) and ".tif" (better support for metadata). +#' @param Compression Integer between 1 to 9. Applied to final .nc file that the function writes to hard drive. Same as compression argument in terra::writeCDF(). Ignored if FileExtension = ".tif". #' @param API_Key Character; ECMWF cds API key. #' @param API_User Character; ECMWF cds user number. -#' @param verbose Optional, logical. Whether to report progress of the function in the console or not. -#' @param Keep_Raw Logical. Whether to keep monthly netcdf files of raw data aggregated to temporal resolution of `T_res`. Default FALSE. +#' @param TChunkSize Numeric. Number of layers to bundle in each individual download. Default is 6000 to adhere to most restrictive CDS limits: https://cds.climate.copernicus.eu/live/limits. +#' @param Cores Numeric. How many cores to use when carrying out temporal aggregation. Default is 1. +#' @param verbose Logical. Whether to print/message function progress in console or not. +#' @param Keep_Raw Logical. Whether to retain raw downloaded data or not. Default is FALSE. #' @param Keep_Monthly Logical. Whether to keep monthly netcdf files of raw data aggregated to temporal resolution of months. Default FALSE. -#' @param Cores Numeric. How many cores to use.^This can speed up downloads of long time-series. If you want output to your console during the process, use Cores = 1. Parallel processing is carried out when Cores is bigger than 1. Default is 1. -#' @param TryDown Optional, numeric. How often to attempt the download of each individual file that the function queries from the server. This is to circumvent having to restart the entire function when encountering connectivity issues. -#' @param TimeOut Numeric. The timeout for each download in seconds. Default 36000 seconds (10 hours). -#' @param SingularDL Logical. Whether to force download of data in one call to CDS or automatically break download requests into individual monthly downloads. Default is FALSE. -#' @return A raster object containing the downloaded ERA5(-Land) data, and a NETCDF (.nc) file in the specified directory. +#' @param TryDown Optional, numeric. Legacy, ignored when querying data from new CDS (https://cds-beta.climate.copernicus.eu/; this happens when the package version of ecmwfr is >= 2.0.0). How often to attempt the download of each individual file that the function queries from the CDS. This is to circumvent having to restart the entire function when encountering connectivity issues. +#' @param TimeOut Numeric. Legacy, ignored when querying data from new CDS (https://cds-beta.climate.copernicus.eu/; this happens when the package version of ecmwfr is >= 2.0.0). The timeout for each download in seconds. Default 36000 seconds (10 hours). +#' +#' @return A SpatRaster object containing the queried bioclimatic data, and a NETCDF (.nc) file in the specified directory. +#' +#' The SpatRaster contains metadata/attributes as a named vector that can be retrieved with terra::metags(...): +#' \itemize{ +#' \item{Citation}{ - A string which to use for in-line citation of the data product obtained with BioClim}. +#' \item{KrigRCall.X}{ - Arguments passed to the BioClim function that produced the file (API credentials are omitted from these metadata)}. +#' } +#' +#' @importFrom tools file_path_sans_ext +#' @importFrom lubridate year +#' @importFrom lubridate date +#' @importFrom stringr str_pad +#' @importFrom terra rast +#' @importFrom terra time +#' @importFrom terra nlyr +#' @importFrom terra metags +#' @importFrom terra writeRaster +#' @importFrom terra tapp +#' @importFrom terra app +#' @importFrom terra mean +#' @importFrom terra values +#' @importFrom terra subset +#' #' @examples #' \dontrun{ -#' BioClim (Water_Var = "volumetric_soil_water_layer_1", # could also be total_precipitation -#' Y_start = 1981, -#' Y_end = 2015, -#' T_res = "day", -#' Extent = extent(11.8,15.1,50.1,51.7), -#' DataSet = "era5-land", -#' Dir = getwd(), -#' verbose = TRUE, -#' FileName = "NULL", -#' Keep_Raw = FALSE, -#' Keep_Monthly = FALSE, -#' API_User = API_User, -#' API_Key = API_Key) +#' CN_ext <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +#' CN_BC <- BioClim( +#' Temperature_Var = "2m_temperature", +#' Temperature_DataSet = "reanalysis-era5-land", +#' Temperature_Type = NA, +#' Water_Var = "volumetric_soil_water_layer_1", +#' Water_DataSet = "reanalysis-era5-land-monthly-means", +#' Water_Type = "monthly_averaged_reanalysis", +#' Y_start = 1970, +#' Y_end = 1979, +#' TZone = "CET", +#' Extent = CN_ext, +#' Dir = getwd(), +#' FileName = "CN_BC", +#' FileExtension = ".nc", +#' Compression = 9, +#' API_User = API_User, +#' API_Key = API_Key, +#' TChunkSize = 6000, TryDown = 10, TimeOut = 36000, +#' Cores = parallel::detectCores(), +#' verbose = TRUE, +#' Keep_Raw = FALSE, Keep_Monthly = FALSE +#' ) #' } #' #' @export -BioClim <- function(Water_Var = "volumetric_soil_water_layer_1", # could also be total_precipitation - Y_start = 1981, - Y_end = 2015, - T_res = "day", - Extent = extent(-180,180,-90,90), - DataSet = "era5-land", - Dir = getwd(), - verbose = TRUE, - FileName = "NULL", - Keep_Raw = FALSE, - Keep_Monthly = FALSE, - Buffer = .5, - ID = "ID", - API_User = API_User, - API_Key = API_Key, - Cores = 1, - TryDown = 10, - TimeOut = 36000, - SingularDL = FALSE){ - - stop("Function currently deprecated as KrigR undergoes major re-development. Please use the stable release to gain access to this functionality.") - - Vars <- c("2m_temperature", Water_Var) - - if(Y_end == year(Sys.Date())){ - stop("Please note that calculation of bioclimatic variables requires consideration of data sets spanning full years. Therefore, the current year cannot be included in the calculation of bioclimatic variables.") +BioClim <- function( + # nolint: cyclocomp_linter. + Temperature_Var = "2m_temperature", + Temperature_DataSet = "reanalysis-era5-land", + Temperature_Type = NA, + Water_Var = "volumetric_soil_water_layer_1", # could also be total_precipitation + Water_DataSet = "reanalysis-era5-land-monthly-means", + Water_Type = "monthly_averaged_reanalysis", + Y_start, Y_end, TZone = "UTC", # time-window, default set to range of dataset-type + Extent, # spatial limitation, default set to range of dataset-type + Buffer = 0.5, # point buffering if desired + Dir = getwd(), FileName, FileExtension = ".nc", Compression = 9, # file storing + API_User, API_Key, # API credentials + TChunkSize = 6000, TryDown = 10, TimeOut = 36000, # Calls to CDS + Cores = 1, # parallelisation + verbose = TRUE, # verbosity + Keep_Raw = FALSE, Keep_Monthly = FALSE # continued file storage + ) { + ## Catching Most Frequent Issues =============== + on.exit(closeAllConnections()) + #--- File Name and Extension + ### check if file name has been specified + if (!exists("FileName")) { + stop("Please provide a value for the FileName argument.") } - - ####### SHAPE HANDLING ####### - if(class(Extent) == "data.frame"){ # if we have been given point data - Extent <- KrigR:::buffer_Points(Points = Extent, Buffer = Buffer, ID = ID) + FileName <- paste0(file_path_sans_ext(FileName), FileExtension) + if (!(FileExtension %in% c(".nc", ".tif"))) { + stop("Please specify a FileExtension of either '.tif' or '.nc'") } - if(class(Extent) == "Raster" | class(Extent) == "SpatialPolygonsDataFrame" | class(Extent) == "SpatialPolygons"){ # sanity check: ensure it is a raster of Spatialpolygonsdataframe object if not an extent object - if(class(Extent) == "SpatialPolygonsDataFrame" | class(Extent) == "SpatialPolygons"){ # shape check - Shape <- Extent # save the shapefile for later masking - } # end of shape check - } # end of sanity check - - - ####### DATA RETRIEVAL ####### - ### DATE HANDLER ---- - Down_start <- lubridate::date(paste0(Y_start, "-01-01")) - Down_end <- lubridate::date(paste0(Y_end, "-12-31")) - T_seq <- seq(Down_start, Down_end, by = "month") - Y_seq <- year(T_seq) - M_seq <- str_pad(month(T_seq), 2, "left", 0) - - ## LOOP FOR EACH MONTH (immediate reducing of raster layers for storage purposes) or in one-go (reducing of download calls) - looptext <- " - # check name (file(s) that will be written) - if(SingularDL){ - CheckName <- file.path(Dir, paste0(Var_down, '-', Fun_vec, 'MonthlyBC.nc')) - TempName <- file.path(Dir, paste0(Var_down, '_Temporary.nc')) - }else{ - CheckName <- file.path(Dir, paste0(Var_down, '-', Fun_vec, '-', Y_seq[Down_Iter], '_', M_seq[Down_Iter], 'MonthlyBC.nc')) - TempName <- file.path(Dir, paste0(Var_down, '_Temporary_', Y_seq[Down_Iter], '_', M_seq[Down_Iter], '.nc')) + #--- Time Zone + if (!(TZone %in% OlsonNames())) { + stop("The TZone argument you have specified is not supported. Please refer to OlsonNames() for an overview of all supported specifications.") } - # DATA CHECK (skip this iteration if data is already downloaded) - if(all(file.exists(CheckName))){ - if(verbose){ - message(ifelse(SingularDL, paste(Var_down, 'already processed'), paste0(Var_down, ' already processed for ', M_seq[Down_Iter], '/', Y_seq[Down_Iter]))) + #--- File Check + Meta_vec <- c(Temperature_Var, Temperature_DataSet, Temperature_Type, Water_Var, Water_DataSet, Water_Type, Y_start, Y_end, TZone, as.character(quote(Extent)), Buffer, Dir, FileName) + names(Meta_vec) <- c("Temperature_Var", "Temperature_DataSet", "Temperature_Type", "Water_Var", "Water_DataSet", "Water_Type", "Y_start", "Y_end", "TZone", "Extent", "Buffer", "Dir", "FileName") + Meta_vec <- c( + "Citation" = paste0("Bioclimatic variables obtained with KrigR (DOI:10.1088/1748-9326/ac48b3) on ", Sys.time()), + "KrigRCall" = Meta_vec + ) + FCheck <- Check.File(FName = FileName, Dir = Dir, loadFun = terra::rast, load = TRUE, verbose = TRUE) + if (!is.null(FCheck)) { + names(FCheck) <- paste0("BIO", 1:19) + if (FileExtension == ".nc") { + FCheck <- Meta.NC(NC = FCheck, FName = file.path(Dir, FileName), Attrs = Meta_vec, Read = TRUE) } - next() + terra::time(FCheck) <- as.POSIXct(terra::time(FCheck), tz = TZone) # assign the correct time zone, when loading from disk, time zone is set to UTC + names(FCheck) <- paste0("BIO", 1:19) + return(FCheck) + } + + #--- API Credentials + ### checking if API User and Key have been supplied + if (exists("API_User") + exists("API_Key") != 2) { + stop("Please provide a value for the API_User and API_Key arguments.") } - ## DATE HANDLER - if(SingularDL){ - month_start <- paste0(Y_start, '-01-01') - month_end <- paste0(Y_end, '-12-31') - }else{ - month_start <- lubridate::date(paste(Y_seq[Down_Iter], M_seq[Down_Iter], '01', sep = '-')) # set start date - month_end <- month_start+(lubridate::days_in_month(month_start)-1) # find end date depending on days in current month + ### making API_User into a character string + API_User <- as.character(API_User) + + #--- Variable and product list + Vars_ls <- list( + list( + Var = Temperature_Var, + DataSet = Temperature_DataSet, + Type = Temperature_Type + ), + list( + Var = Water_Var, + DataSet = Water_DataSet, + Type = Water_Type + ) + ) + + #--- DataSets and Types + MetaCheck <- lapply(lapply(Vars_ls, "[[", "DataSet"), Meta.QuickFacts) + + ## Projection mismatch + ProjCheck <- sapply(1:length(MetaCheck), FUN = function(x) { + MetaCheck[[x]]$Projection + }) + if (length(unique(ProjCheck)) > 1) { + stop("The CDS data products you have specified are stored in different projections. Please select data products via the *_DataSet arguments that use the same spatial projection.") } - ## ensuring that water var is pulled at monthly resolution - if(Var_Iter == 2 ){ - T_resDL <- 'month' - }else{ - T_resDL <- T_res + ## extent checking + ExtCheck <- lapply(1:length(MetaCheck), FUN = function(x) { + ext(MetaCheck[[x]]$CDSArguments$area) + }) + if (!all.equal(ExtCheck[[1]], ExtCheck[[2]])) { + stop("Please specify CDS products via the *_DataSet arguments which have matching spatial coverage.") } - ## DOWNLOAD - if(file.exists(TempName)){ - if(verbose){ - message(ifelse(SingularDL, paste(Var_down, 'already downloaded'), paste0(Var_down, ' already downloaded for ', M_seq[Down_Iter], '/', Y_seq[Down_Iter]))) - } - Temp_Ras <- stack(TempName) - }else{ - Temp_Ras <- download_ERA( - Variable = Var_down, - DataSet = DataSet, - Type = 'reanalysis', - DateStart = month_start, - DateStop = month_end, - TResolution = T_resDL, - TStep = 1, - FUN = AggrFUN, - Extent = Extent, - Dir = Dir, - FileName = strsplit(TempName, '/')[[1]][length(strsplit(TempName, '/')[[1]])], - API_User = API_User, - API_Key = API_Key, - verbose = verbose, - PrecipFix = PrecipFix, - TryDown = TryDown, - TimeOut = TimeOut, - SingularDL = SingularDL - ) - } - ## PROCESSING - if(Var_Iter == 1){ - Counter <- 1 - for(Iter_fun in Fun_vec){ - Save_Ras <- stackApply(Temp_Ras, - indices = month(seq(as.POSIXct(paste(month_start, '00:00:00')), as.POSIXct(paste(month_end, '23:00:00')), by=T_res)), - fun = Iter_fun) - - if(Iter_fun == 'sum' & exists('Shape')){ - range <- KrigR:::mask_Shape(base.map = Save_Ras[[1]], Shape = Shape) - Save_Ras <- mask(Save_Ras, range) - } - terra::writeCDF(x = as(brick(Save_Ras), 'SpatRaster'), filename = CheckName[Counter], overwrite = TRUE) - # writeRaster(x = Save_Ras, - # filename = CheckName[Counter], - # format = 'CDF', overwrite = TRUE) - Counter <- Counter + 1 - } - }else{ - file.copy(from = TempName, to = CheckName) + ## resolution mismatch + ResCheck <- diff( + sapply(1:length(MetaCheck), FUN = function(x) { + MetaCheck[[x]]$SpatialResolution + }) + ) + if (ResCheck != 0) { + stop("Please specify CDS products via the *_DataSet arguments which have matching spatial resolution.") + } + + #--- Year Specificiation + if (Y_end == year(Sys.Date())) { + stop("Please note that calculation of bioclimatic variables requires consideration of data sets spanning full years. Therefore, the current year cannot be included in the calculation of bioclimatic variables.") } - ## DELETING RAW - if(!isTRUE(Keep_Raw)){ - unlink(TempName) - } - " - ### DATA DOWNLOAD ---- - for(Var_Iter in 1:length(Vars)){ ## LOOP FOR EACH VARIABLE + ## Data Retrieval =============== + message("###### Retrieving raw data ######") + #--- Date handling + Down_start <- paste0(Y_start, "-01-01 00:00") + Down_end <- paste0(Y_end, "-12-31 23:00") + M_seq <- seq(lubridate::date(Down_start), lubridate::date(Down_end), by = "month") + 15 # this is to ensure proper downloading of monthly recorded data despite TZone adjustments + #--- Data download + RawNames <- paste(file_path_sans_ext(FileName), sapply(Vars_ls, "[[", "Var"), "RAW", sep = "_") + Raw_data <- lapply(1:length(Vars_ls), FUN = function(Var_Iter) { ## VARIABLE IDENTIFICATION - Var_down <- Vars[Var_Iter] - if(Var_down == "total_precipitation"){PrecipFix <- TRUE}else{PrecipFix <- FALSE} + Var_down <- Vars_ls[[Var_Iter]]$Var + CumulVar <- ifelse(startsWith(prefix = "total_", Var_down), TRUE, FALSE) ## PROCESSING FUNCTIONS - Fun_vec <- 'mean' # for all water variables that are not total precip - if(Var_down == 'total_precipitation'){Fun_vec <- 'sum'} - if(Var_down == '2m_temperature'){Fun_vec <- c('min', 'mean', 'max')} - - if(Var_down == 'total_precipitation'){ + if (startsWith(prefix = "total_", Var_down)) { AggrFUN <- sum - }else{ + } else { # for all variables that are not total precip AggrFUN <- mean } - if(SingularDL){ - n_down <- 1 - Cores = 1 - }else{ - n_down <- length(T_seq) + + ## CDownloadS CALL + Raw_rast <- CDownloadS( + Variable = Var_down, + CumulVar = CumulVar, + DataSet = Vars_ls[[Var_Iter]]$DataSet, + Type = Vars_ls[[Var_Iter]]$Type, + DateStart = Down_start, + DateStop = Down_end, + TZone = TZone, + TResolution = ifelse(MetaCheck[[Var_Iter]]$TResolution == "month", "month", "day"), + FUN = AggrFUN, + Extent = Extent, + Buffer = Buffer, + Dir = Dir, + FileName = RawNames[Var_Iter], + FileExtension = FileExtension, + Compression = Compression, + API_User = API_User, + API_Key = API_Key, + TryDown = TryDown, + TimeOut = TimeOut, + TChunkSize = TChunkSize, + Cores = Cores, + verbose = verbose, + Keep_Raw = FALSE + ) + Raw_rast + }) + + ## Data Processing =============== + message("###### Prcoessing raw data ######") + #--- Month summaries + if (startsWith(prefix = "total_", Water_Var)) { + AggrFUN <- sum + AggrFunName <- "sum" + } else { # for all variables that are not total precip + AggrFUN <- mean + AggrFunName <- "mean" + } + MonthNames <- paste0(paste(file_path_sans_ext(FileName), c(rep("Temperature", 3), "Water"), "Monthly", c("mean", "min", "max", AggrFunName), sep = "_"), ".nc") + Monthly <- list( + list( + Source = 1, + AggrFUN = mean, + Name = MonthNames[1] + ), + list( + Source = 1, + AggrFUN = min, + Name = MonthNames[2] + ), + list( + Source = 1, + AggrFUN = max, + Name = MonthNames[3] + ), + list( + Source = 2, + AggrFUN = AggrFUN, + Name = MonthNames[4] + ) + ) + names(Monthly) <- c("T_mean", "T_min", "T_max", "W_mean") + Monthly <- lapply(Monthly, FUN = function(Monthly_Iter) { + if (verbose) { + print(paste("Producing", Monthly_Iter$Name)) } - n_downrep <- n_down - if(Var_down == 'total_precipitation' & !SingularDL){ - n_downrep <- length(T_seq)*2} - if(verbose){ - message(paste("The KrigR::BioClim() function is going to stage", n_downrep, "download(s) for", Var_down, "data now.")) + MetaMonthly_vec <- metags(Raw_data[[Monthly_Iter$Source]]) + MetaMonthly_vec["KrigRCall.TResolution"] <- "month" + MetaMonthly_vec["KrigRCall.FUN"] <- "Monthly_Iter$AggrFUN" + FCheck <- Check.File(FName = Monthly_Iter$Name, Dir = Dir, loadFun = terra::rast, load = TRUE, verbose = TRUE) + if (!is.null(FCheck)) { + if (FileExtension == ".nc") { + FCheck <- Meta.NC(NC = FCheck, FName = file.path(Dir, Monthly_Iter$Name), Attrs = MetaMonthly_vec, Read = TRUE) + } + terra::time(FCheck) <- as.POSIXct(terra::time(FCheck), tz = TZone) # assign the correct time zone, when loading from disk, time zone is set to UTC + return(FCheck) } + Monthly_rast <- Temporal.Aggr( + CDS_rast = Raw_data[[Monthly_Iter$Source]], + BaseResolution = "day", BaseStep = 1, + TResolution = "month", TStep = 1, FUN = Monthly_Iter$AggrFUN, + Cores = Cores, QueryTargetSteps = NULL, TZone = TZone, verbose = FALSE + ) + terra::metags(Monthly_rast) <- MetaMonthly_vec + if (FileExtension == ".tif") { + terra::writeRaster(Monthly_rast, filename = file.path(Dir, Monthly_Iter$Name)) + Monthly_rast <- terra::rast(filename = file.path(Dir, Monthly_Iter$Name)) + } + if (FileExtension == ".nc") { + Monthly_rast <- Meta.NC( + NC = Monthly_rast, FName = file.path(Dir, Monthly_Iter$Name), + Attrs = MetaMonthly_vec, Write = TRUE, + Compression = Compression + ) + terra::time(Monthly_rast) <- as.POSIXct(terra::time(Monthly_rast), tz = TZone) # assign the correct time zone, when loading from disk, time zone is set to UTC + } + Monthly_rast + }) + names(Monthly) <- c("T_mean", "T_min", "T_max", "W_mean") - if(Cores > 1){ # Cores check: if parallel processing has been specified - ForeachObjects <- c("Var_down", "Var_Iter", "Dir", "Y_seq", "M_seq", "DataSet", "PrecipFix", "API_User", "API_Key", "T_res", "Extent", "Keep_Raw", "Fun_vec", "TryDown", "TimeOut", "SingularDL", "Var_down", "Fun_vec", "AggrFUN", "verbose", "Down_start") - pb <- txtProgressBar(max = n_down, style = 3) - progress <- function(n){setTxtProgressBar(pb, n)} - opts <- list(progress = progress) - cl <- makeCluster(Cores) # Assuming Cores node cluster - registerDoSNOW(cl) # registering cores - foreach(Down_Iter = 1:n_down, - .packages = c("KrigR"), # import packages necessary to each itteration - .export = ForeachObjects, - .options.snow = opts) %:% when(!file.exists(file.path(Dir, paste0(Var_down, '-', Fun_vec[length(Fun_vec)], '-', Y_seq[Down_Iter], '_', M_seq[Down_Iter], 'MonthlyBC.nc')))) %dopar% { - eval(parse(text=looptext)) - } # end of parallel kriging loop - close(pb) - stopCluster(cl) # close down cluster - }else{ # if non-parallel processing has been specified - for(Down_Iter in 1:n_down){eval(parse(text=looptext))} - } # end of non-parallel loop - } # end of variable loop - - ### DATA LOADING ---- - setwd(Dir) - Tair_min <- raster::stack(list.files(path = Dir, pattern = paste0(Vars[1], "-min"))) - Tair_mean <- raster::stack(list.files(path = Dir, pattern = paste0(Vars[1], "-mean"))) - Tair_max <- raster::stack(list.files(path = Dir, pattern = paste0(Vars[1], "-max"))) - Water <- raster::stack(list.files(path = Dir, pattern = paste0(Vars[2], "-", Fun_vec))) - - ### QUARTER COMPUTATION ---- - Tair_mean_quarter <- stackApply(Tair_mean, indices = rep(1:4, each = 3, length.out = nlayers(Tair_mean)), mean) - if(Water_Var == "total_precipitation"){ - Water_quarter <- stackApply(Water, indices = rep(1:4, each = 3, length.out = nlayers(Water)), sum) - }else{ - Water_quarter <- stackApply(Water, indices = rep(1:4, each = 3, length.out = nlayers(Water)), mean) - } - if(Water_Var == "total_precipitation" & exists("Shape")){ - range <- KrigR:::mask_Shape(base.map = Water_quarter[[1]], Shape = Shape) - Water_quarter <- mask(Water_quarter, range) - } + #--- Quarter summaries + T_mean_quarter <- tapp(x = Monthly$T_mean, index = rep(1:4, each = 3, length.out = nlyr(Monthly$T_mean)), fun = mean, cores = Cores) + W_mean_quarter <- tapp(x = Monthly$W_mean, index = rep(1:4, each = 3, length.out = nlyr(Monthly$W_mean)), fun = AggrFUN, cores = Cores) - ####### BIOCLIMATIC VARIABLES ####### + ## Bioclimatic Variables =============== + message("###### Calculating aggregate metrics ######") ### BIO1 = Annual Mean Temperature ---- - BIO1 <- raster::mean(Tair_mean) + BIO1 <- mean(Monthly$T_mean) ### BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp)) ---- - BIO2 <- raster::mean(Tair_max-Tair_min, na.rm = TRUE) + BIO2 <- mean(Monthly$T_max - Monthly$T_min, na.rm = TRUE) ### BIO4 = Temperature Seasonality (standard deviation * 100) ---- - BIO4 <- calc(Tair_mean, sd)*100 + BIO4 <- app(Monthly$T_mean, sd) * 100 + + ### Monthly Temperature ---- + T_mean <- values(Monthly$T_mean) + T_max <- values(Monthly$T_max) + T_min <- values(Monthly$T_min) + BIO5 <- BIO6 <- T_mean_quarter[[1]] + bio5vals <- bio6vals <- values(BIO5)[, 1] + + for (i in seq_along(bio5vals)) { + if (is.na(bio5vals[i])) { + next() + } + bio5vals[i] <- T_max[i, which.max(T_mean[i, ])] + bio6vals[i] <- T_min[i, which.min(T_mean[i, ])] + } - ### BIO5 = Max Temperature of Warmest Month ---- - BIO5 <- max(Tair_max) + #### BIO5 = Max Temperature of Warmest Month ---- + terra::values(BIO5) <- bio5vals - ### BIO6 = Min Temperature of Coldest Month ---- - BIO6 <- min(Tair_min) + #### BIO6 = Min Temperature of Coldest Month ---- + terra::values(BIO6) <- bio6vals ### BIO7 = Temperature Annual Range (BIO5-BIO6) ---- - BIO7 <- BIO5-BIO6 + BIO7 <- BIO5 - BIO6 ### BIO3 = Isothermality (BIO2/BIO7) (*100) ---- - BIO3 <- BIO2/BIO7*100 - - ### BIO8 = Mean Temperature of Wettest Quarter ---- - BIO8 <- raster::stackSelect(Tair_mean_quarter, raster::which.max(Water_quarter)) - - ### BIO9 = Mean Temperature of Driest Quarter ---- - BIO9 <- raster::stackSelect(Tair_mean_quarter, raster::which.min(Water_quarter)) + BIO3 <- BIO2 / BIO7 * 100 ### BIO10 = Mean Temperature of Warmest Quarter ---- - BIO10 <- raster::stackSelect(Tair_mean_quarter, raster::which.max(Tair_mean_quarter)) + BIO10 <- max(T_mean_quarter) ### BIO11 = Mean Temperature of Coldest Quarter ---- - BIO11 <- raster::stackSelect(Tair_mean_quarter, raster::which.min(Tair_mean_quarter)) + BIO11 <- min(T_mean_quarter) ### BIO12 = Annual Precipitation ---- - if(Water_Var == "total_precipitation"){ - BIO12 <- sum(Water)/(nlayers(Water)/12) - }else{ - BIO12 <- raster::mean(Water) - } - if(Water_Var == "total_precipitation" & exists("Shape")){ - range <- KrigR:::mask_Shape(base.map = BIO12, Shape = Shape) - BIO12 <- mask(BIO12, range) + if (startsWith(prefix = "total_", Water_Var)) { + BIO12 <- sum(Monthly$W_mean) / (nlyr(Monthly$W_mean) / 12) + } else { + BIO12 <- mean(Monthly$W_mean) } ### BIO13 = Precipitation of Wettest Month ---- - BIO13 <- max(Water) - if(Water_Var == "total_precipitation" & exists("Shape")){ - range <- KrigR:::mask_Shape(base.map = BIO13, Shape = Shape) - BIO13 <- mask(BIO13, range) - } + BIO13 <- max(Monthly$W_mean) ### BIO14 = Precipitation of Driest Month ---- - BIO14 <- min(Water) - if(Water_Var == "total_precipitation" & exists("Shape")){ - range <- KrigR:::mask_Shape(base.map = BIO14, Shape = Shape) - BIO14 <- mask(BIO14, range) - } + BIO14 <- min(Monthly$W_mean) ### BIO15 = Precipitation Seasonality (Coefficient of Variation) ---- - BIO15 <- calc(Water, sd)/BIO12 * 100 + BIO15 <- app(Monthly$W_mean, sd) / BIO12 * 100 ### BIO16 = Precipitation of Wettest Quarter ---- - BIO16 <- max(Water_quarter) + BIO16 <- max(W_mean_quarter) ### BIO17 = Precipitation of Driest Quarter ---- - BIO17 <- min(Water_quarter) + BIO17 <- min(W_mean_quarter) - ### BIO18 = Precipitation of Warmest Quarter ---- - BIO18 <- raster::stackSelect(Water_quarter, raster::which.max(Tair_mean_quarter)) + ### Quarterly ---- + T_q <- values(T_mean_quarter) + W_q <- values(W_mean_quarter) + BIO8 <- BIO9 <- BIO18 <- BIO19 <- T_mean_quarter[[1]] + bio8vals <- bio9vals <- bio18vals <- bio19vals <- values(BIO8)[, 1] - ### BIO19 = Precipitation of Coldest Quarter ---- - BIO19 <- raster::stackSelect(Water_quarter, raster::which.min(Tair_mean_quarter)) + for (i in seq_along(bio8vals)) { + if (is.na(bio8vals[i])) { + next() + } + bio8vals[i] <- T_q[i, which.max(W_q[i, ])] + bio9vals[i] <- T_q[i, which.min(W_q[i, ])] + bio18vals[i] <- W_q[i, which.max(T_q[i, ])] + bio19vals[i] <- W_q[i, which.min(T_q[i, ])] + } - ####### EXPORT ####### - BIO_Ras <-raster::stack(BIO1, BIO2, BIO3, BIO4, BIO5, BIO6, BIO7, BIO8, BIO9, - BIO10, BIO11, BIO12, BIO13, BIO14, BIO15, BIO16, BIO17, BIO18, BIO19) - names(BIO_Ras) <- paste0("BIO", 1:19) + #### BIO8 = Mean Temperature of Wettest Quarter ---- + terra::values(BIO8) <- bio8vals - terra::writeCDF(x = as(brick(BIO_Ras), 'SpatRaster'), filename = file.path(Dir, FileName), overwrite = TRUE) + #### BIO9 = Mean Temperature of Driest Quarter ---- + terra::values(BIO9) <- bio9vals - # writeRaster(BIO_Ras, file.path(Dir, FileName), format = "CDF", overwrite = TRUE) - if(!isTRUE(Keep_Monthly)){ - RM_fs <- list.files(Dir, pattern = "MonthlyBC.nc") - unlink(RM_fs) + ### BIO18 = Precipitation of Warmest Quarter ---- + terra::values(BIO18) <- bio18vals + + ### BIO19 = Precipitation of Coldest Quarter ---- + terra::values(BIO19) <- bio19vals + + ## Data Export =============== + message("###### Data Export ######") + BIO_rast <- c( + BIO1, BIO2, BIO3, BIO4, BIO5, BIO6, BIO7, BIO8, BIO9, + BIO10, BIO11, BIO12, BIO13, BIO14, BIO15, BIO16, BIO17, BIO18, BIO19 + ) + names(BIO_rast) <- paste0("BIO", 1:19) + terra::metags(BIO_rast) <- Meta_vec + + ### write file + if (FileExtension == ".tif") { + terra::writeRaster(BIO_rast, filename = file.path(Dir, FileName)) + BIO_rast <- terra::rast(filename = file.path(Dir, FileName)) } - return(BIO_Ras) + if (FileExtension == ".nc") { + BIO_rast <- Meta.NC( + NC = BIO_rast, FName = file.path(Dir, FileName), + Attrs = Meta_vec, Write = TRUE, + Compression = Compression + ) + } + names(BIO_rast) <- paste0("BIO", 1:19) + + ### cleaning up + if (!Keep_Monthly) { + unlink(file.path(Dir, MonthNames)) + } + if (!Keep_Raw) { + unlink(file.path(Dir, paste0(RawNames, ".nc"))) + } + + return(BIO_rast) } diff --git a/R/CDSAPI.R b/R/CDSAPI.R index c0f03d1..47aeaa1 100644 --- a/R/CDSAPI.R +++ b/R/CDSAPI.R @@ -13,14 +13,37 @@ #' #' @seealso \code{\link{Make.Request}}, \code{\link{Execute.Requests}}. #' -Register.Credentials <- function(API_User, API_Key){ - API_Service = "cds" - KeyRegisterCheck <- tryCatch(ecmwfr::wf_get_key(user = API_User, service = API_Service), - error = function(e){e}) - if(any(class(KeyRegisterCheck) == "simpleError")){ - ecmwfr::wf_set_key(user = API_User, - key = as.character(API_Key), - service = API_Service) +Register.Credentials <- function(API_User, API_Key) { + if (packageVersion("ecmwfr") < "2.0.0") { + warning("You are using an ecmwfr (a KrigR dependency) of version < 2.0.0. This causes queries to be directed to the old CDS service (https://cds.climate.copernicus.eu) instead of the new CDS (https://cds-beta.climate.copernicus.eu/). You may want to update ecmwfr to the latest version to ensure reliable downloads of CDS data going forward.") + API_Service <- "cds" + KeyRegisterCheck <- tryCatch(ecmwfr::wf_get_key(user = API_User, service = API_Service), + error = function(e) { + e + } + ) + if (any(class(KeyRegisterCheck) == "simpleError")) { + ecmwfr::wf_set_key( + user = API_User, + key = as.character(API_Key), + service = API_Service + ) + } + } else { + if (!grepl("@", API_User)) { + stop("With the adoption of the new CDS (https://cds-beta.climate.copernicus.eu/), API_User must be you E-mail registered with the new CDS.") + } + KeyRegisterCheck <- tryCatch(ecmwfr::wf_get_key(user = API_User), + error = function(e) { + e + } + ) + if (any(class(KeyRegisterCheck) == "simpleError")) { + ecmwfr::wf_set_key( + user = API_User, + key = as.character(API_Key) + ) + } } } ### FORMING CDS Requests ======================================================= @@ -39,6 +62,7 @@ Register.Credentials <- function(API_User, API_Key){ #' @param verbose Logical. Whether to print/message function progress in console or not. #' @param API_User Character. CDS API User #' @param API_Key Character. CDS API Key +#' @param TimeOut Numeric. Legacy, ignored when querying data from new CDS (https://cds-beta.climate.copernicus.eu/; this happens when the package version of ecmwfr is >= 2.0.0). The timeout for each download in seconds. Default 36000 seconds (10 hours). #' #' @importFrom ecmwfr wf_request #' @@ -48,74 +72,76 @@ Register.Credentials <- function(API_User, API_Key){ #' Make.Request <- function(QueryTimeWindows, QueryDataSet, QueryType, QueryVariable, QueryTimes, QueryExtent, QueryFormat, Dir = getwd(), verbose = TRUE, - API_User, API_Key){ + API_User, API_Key, TimeOut = 36000) { #' Make list of CDS Requests - Requests_ls <- lapply(1:length(QueryTimeWindows), FUN = function(requestID){ + Requests_ls <- lapply(1:length(QueryTimeWindows), FUN = function(requestID) { FName <- paste("TEMP", QueryVariable, stringr::str_pad(requestID, 5, "left", "0"), sep = "_") - if(grepl("month", QueryType)){ # monthly data needs to be specified with year, month fields - list('dataset_short_name' = QueryDataSet, - 'product_type' = QueryType, - 'variable' = QueryVariable, - 'year' = unique(as.numeric(format(as.POSIXct(QueryTimeWindows[[requestID]]), "%Y"))), - 'month' = unique(as.numeric(format(QueryTimeWindows[[requestID]], "%m"))), - 'time' = QueryTimes, - 'area' = QueryExtent, - 'format' = QueryFormat, - 'target' = FName - # , - # 'grid' = QueryGrid + if (grepl("month", QueryType)) { # monthly data needs to be specified with year, month fields + list( + "dataset_short_name" = QueryDataSet, + "product_type" = QueryType, + "variable" = QueryVariable, + "year" = unique(as.numeric(format(as.POSIXct(QueryTimeWindows[[requestID]]), "%Y"))), + "month" = unique(as.numeric(format(QueryTimeWindows[[requestID]], "%m"))), + "time" = QueryTimes, + "area" = QueryExtent, + "format" = QueryFormat, + "target" = FName ) - }else{ - list('dataset_short_name' = QueryDataSet, - 'product_type' = QueryType, - 'variable' = QueryVariable, - 'date' = paste0( - head(QueryTimeWindows[[requestID]], n = 1), - "/", - tail(QueryTimeWindows[[requestID]], n = 1)), - 'time' = QueryTimes, - 'area' = QueryExtent, - 'format' = QueryFormat, - 'target' = FName - # , - # 'grid' = QueryGrid + } else { + list( + "dataset_short_name" = QueryDataSet, + "product_type" = QueryType, + "variable" = QueryVariable, + "date" = paste0( + head(QueryTimeWindows[[requestID]], n = 1), + "/", + tail(QueryTimeWindows[[requestID]], n = 1) + ), + "time" = QueryTimes, + "area" = QueryExtent, + "format" = QueryFormat, + "target" = FName ) } - }) ## making list names useful for request execution updates to console Iterators <- paste0("[", 1:length(Requests_ls), "/", length(Requests_ls), "] ") FNames <- unlist(lapply(Requests_ls, "[[", "target")) Dates <- unlist(lapply(lapply(Requests_ls, "[[", "date"), gsub, pattern = "/", replacement = " - ")) - if(length(Dates) == 0){ # this happens for monthly data queries - Dates <- unlist(lapply(lapply(Requests_ls, "[[", "year"), FUN = function(x){ + if (length(Dates) == 0) { # this happens for monthly data queries + Dates <- unlist(lapply(lapply(Requests_ls, "[[", "year"), FUN = function(x) { paste0(head(x, 1), " - ", tail(x, 1)) })) } names(Requests_ls) <- paste0(Iterators, FNames, " (UTC: ", Dates, ")") ## check if files are already present - FCheck <- sapply(FNames, Check.File, Dir = Dir, loadFun = "terra::rast", load = FALSE, - verbose = FALSE) - if(length(names(unlist(FCheck))) > 0){ + FCheck <- sapply(FNames, Check.File, + Dir = Dir, loadFun = "terra::rast", load = FALSE, + verbose = FALSE + ) + if (length(names(unlist(FCheck))) > 0) { Requests_ls[match(names(unlist(FCheck)), FNames)] <- NA } Requests_ls - if(verbose){print("## Staging CDS Requests")} - for(requestID in 1:length(Requests_ls)){ ## looping over CDS requests - if(verbose){print(names(Requests_ls)[requestID])} - if(class(Requests_ls[[requestID]]) == "logical"){ - # if(verbose){print("File with this name is already present.")} + if (verbose) { + print("## Staging CDS Requests") + } + for (requestID in 1:length(Requests_ls)) { ## looping over CDS requests + if (verbose) { + print(names(Requests_ls)[requestID]) + } + if (class(Requests_ls[[requestID]]) == "logical") { next() } - API_request <- suppressMessages({ - ecmwfr::wf_request(user = API_User, - request = Requests_ls[[requestID]], - transfer = FALSE, - path = Dir, - verbose = verbose, - time_out = TimeOut) - }) + API_request <- ecmwfr::wf_request( + user = API_User, + request = Requests_ls[[requestID]], + transfer = FALSE, + path = Dir, + verbose = FALSE + ) Requests_ls[[requestID]]$API_request <- API_request } Requests_ls @@ -137,69 +163,154 @@ Make.Request <- function(QueryTimeWindows, QueryDataSet, QueryType, QueryVariabl #' @importFrom httr DELETE #' @importFrom httr authenticate #' @importFrom httr add_headers +#' @importFrom ecmwfr wf_delete #' #' @return No R object. Resulting files of CDS query/queries in signated directory. #' #' @seealso \code{\link{Register.Credentials}}, \code{\link{Make.Request}}. #' -Execute.Requests <- function(Requests_ls, Dir, API_User, API_Key, TryDown, verbose = TRUE){ - if(verbose){print("## Listening for CDS Requests")} - for(requestID in 1:length(Requests_ls)){ ## looping over CDS requests - if(verbose){print(names(Requests_ls)[requestID])} - if(class(Requests_ls[[requestID]]) == "logical"){ - if(verbose){print("File with this name is already present.")} +Execute.Requests <- function(Requests_ls, Dir, API_User, API_Key, TryDown, verbose = TRUE) { # nolint: cyclocomp_linter. + if (verbose) { + print("## Listening for CDS Requests") + } + + for (requestID in 1:length(Requests_ls)) { ## looping over CDS requests + if (verbose) { + print(names(Requests_ls)[requestID]) + } + + if (class(Requests_ls[[requestID]]) == "logical") { + if (verbose) { + print("File with this name is already present.") + } next() } API_request <- Requests_ls[[requestID]]$API_request - FileDown <- list(state = "queued") - Down_try <- 0 - while(FileDown$state != "completed" & Down_try <= TryDown){ - ## console output that shows the status of the request on CDS - if(verbose){ - if(FileDown$state == "queued"){ - for(rep_iter in 1:10){ - cat(rep(" ", 100)) - flush.console() - cat('\r', "Waiting for CDS to start processing the query", rep(".", rep_iter)) - flush.console() - Sys.sleep(0.5) + + ## old CDS + if (packageVersion("ecmwfr") < "2.0.0") { + FileDown <- list(state = "queued") + Down_try <- 0 + while (FileDown$state != "completed" && Down_try <= TryDown) { + ## console output that shows the status of the request on CDS + if (verbose) { + if (FileDown$state == "queued") { + for (rep_iter in 1:10) { + cat(rep(" ", 100)) + flush.console() + cat("\r", "Waiting for CDS to start processing the query", rep(".", rep_iter)) + flush.console() + Sys.sleep(0.25) + } + } + if (FileDown$state == "running") { + for (rep_iter in 1:10) { + cat(rep(" ", 100)) + flush.console() + cat("\r", "CDS is processing the query", rep(".", rep_iter)) + flush.console() + Sys.sleep(0.25) + } + } + } + ## download file for current request when ready + FileDown <- tryCatch( + ecmwfr::wf_transfer( + url = API_request$get_url(), + user = API_User, + service = "cds", + verbose = TRUE, + path = Dir, + filename = API_request$get_request()$target + ), + error = function(e) { + e } + ) + if (Down_try == TryDown) { + stop("Download of CDS query result continues to fail after ", Down_try, " trys. The most recent error message is: \n", FileDown, "Assess issues at https://cds.climate.copernicus.eu/cdsapp#!/yourrequests.") + } + if (any(class(FileDown) == "simpleError")) { + FileDown <- list(state = "queued") + Down_try <- Down_try + 1 + } + } + if (FileDown$state == "completed") { + delete <- httr::DELETE( + API_request$get_url(), + httr::authenticate(API_User, API_Key), + httr::add_headers( + "Accept" = "application/json", + "Content-Type" = "application/json" + ) + ) + rm(delete) + } + } else { ## new CDS!!! + FileDown <- list(state = "accepted") + API_request <- API_request$update_status(verbose = FALSE) + + while (FileDown$state != "successful") { + ## console output that shows the status of the request on CDS + if (verbose) { + if (FileDown$state == "accepted") { + for (rep_iter in 1:10) { + cat(rep(" ", 100)) + flush.console() + cat("\r", "Waiting for CDS to start processing the query", rep(".", rep_iter)) + flush.console() + Sys.sleep(0.25) + } + } + if (FileDown$state == "running") { + for (rep_iter in 1:10) { + cat(rep(" ", 100)) + flush.console() + cat("\r", "CDS is processing the query", rep(".", rep_iter)) + flush.console() + Sys.sleep(0.25) + } + } + } + + API_request <- API_request$update_status(verbose = FALSE) + FileDown$state <- API_request$get_status() + + if (API_request$is_failed()) { + stop("Query failed on CDS. Assess issues at https://cds.climate.copernicus.eu/cdsapp#!/yourrequests.") } - if(FileDown$state == "running"){ - for(rep_iter in 1:10){ + + if (FileDown$state == "successful") { + for (rep_iter in 1:10) { cat(rep(" ", 100)) flush.console() - cat('\r', "CDS is processing the query", rep(".", rep_iter)) + cat("\r", "CDS finished processing the query. Download starting soon.", rep(".", rep_iter)) flush.console() - Sys.sleep(0.5) + Sys.sleep(0.25) } + Download_CDS <- capture.output( + ecmwfr::wf_transfer( + url = API_request$get_url(), + user = API_User, + verbose = TRUE, + path = Dir, + filename = API_request$get_request()$target + ), + type = "message" + ) + rm(Download_CDS) + + ## purge request and check succes of doing so + checkdeletion <- capture.output( + wf_delete( + url = API_request$get_url(), + user = API_User + ), + type = "message" + ) + rm(checkdeletion) } } - ## download file for current request when ready - FileDown <- tryCatch(ecmwfr::wf_transfer(url = API_request$get_url(), - user = API_User, - service = "cds", - verbose = TRUE, - path = Dir, - filename = API_request$get_request()$target), - error = function(e){e} - ) - if(Down_try == TryDown){ - stop("Download of CDS query result continues to fail after ", Down_try, " trys. The most recent error message is: \n", FileDown, "Assess issues at https://cds.climate.copernicus.eu/cdsapp#!/yourrequests.") - } - if(any(class(FileDown) == "simpleError")){ - FileDown <- list(state = "queued") - Down_try <- Down_try+1 - } - } - if(FileDown$state == "completed"){ - delete <- httr::DELETE( - API_request$get_url(), - httr::authenticate(API_User, API_Key), - httr::add_headers( - "Accept" = "application/json", - "Content-Type" = "application/json") - ) - } - } -} + } # new CDS end + } # request loop +} # function diff --git a/R/CDownloadS.R b/R/CDownloadS.R index 1ed6150..b008ea2 100644 --- a/R/CDownloadS.R +++ b/R/CDownloadS.R @@ -17,15 +17,15 @@ #' @param Dir Character/Directory Pointer. Directory specifying where to download data to. #' @param FileName Character. A file name for the produced file. #' @param FileExtension Character. A file extension for the produced file. Supported values are ".nc" (default) and ".tif" (better support for metadata). +#' @param Compression Integer between 1 to 9. Applied to final .nc file that the function writes to hard drive. Same as compression argument in terra::writeCDF(). Ignored if FileExtension = ".tif". #' @param API_Key Character; ECMWF cds API key. #' @param API_User Character; ECMWF cds user number. -#' @param TryDown Optional, numeric. How often to attempt the download of each individual file that the function queries from the CDS. This is to circumvent having to restart the entire function when encountering connectivity issues. -#' @param TimeOut Numeric. The timeout for each download in seconds. Default 36000 seconds (10 hours). #' @param TChunkSize Numeric. Number of layers to bundle in each individual download. Default is 6000 to adhere to most restrictive CDS limits: https://cds.climate.copernicus.eu/live/limits. #' @param Cores Numeric. How many cores to use when carrying out temporal aggregation. Default is 1. #' @param verbose Logical. Whether to print/message function progress in console or not. #' @param Keep_Raw Logical. Whether to retain raw downloaded data or not. Default is FALSE. -#' @param Save_Final Logical. Whether to write the final SpatRaster to the hard drive. Default is TRUE. +#' @param TryDown Optional, numeric. Legacy, ignored when querying data from new CDS (https://cds-beta.climate.copernicus.eu/; this happens when the package version of ecmwfr is >= 2.0.0). How often to attempt the download of each individual file that the function queries from the CDS. This is to circumvent having to restart the entire function when encountering connectivity issues. +#' @param TimeOut Numeric. Legacy, ignored when querying data from new CDS (https://cds-beta.climate.copernicus.eu/; this happens when the package version of ecmwfr is >= 2.0.0). The timeout for each download in seconds. Default 36000 seconds (10 hours). #' #' @importFrom tools file_path_sans_ext #' @importFrom terra rast @@ -47,52 +47,52 @@ #' #' \strong{ATTENTION:} If data is loaded again from disk at a later point with a different function, take note that the time zone will have to be set anew and existing time parameters in the .nc contents will need to be converted to the desired time zone. Likewise, citation and KrigR-call metadata will not be loaded properly from a .nc when loading data through a different function. CDownloads() handles these .nc specific issues when loading .nc files created previously with CDownloadS from disk. #' -#' @seealso \code{\link{Meta.List}}, \code{\link{Meta.Variables}}, \code{\link{Meta.QuickFacts}}. +#' @seealso \code{\link{Meta.List}}, \code{\link{Meta.Variables}}, \code{\link{Meta.QuickFacts}}, \code{\link{Plot.SpatRast}}. #' #' @examples #' \dontrun{ #' ## Raw data for one month of full globe #' RawGlobe_rast <- CDownloadS( -#' Variable = "2m_temperature", -#' DataSet = "reanalysis-era5-land-monthly-means", -#' Type = "monthly_averaged_reanalysis", -#' # time-window, default set to range of dataset-type -#' DateStart = "1995-01-01 00:00", -#' DateStop = "1995-01-01 23:00", -#' TZone = "CET", -#' # temporal aggregation -#' TResolution = "month", -#' TStep = 1, -#' # file storing -#' FileName = "RawGlobe", -#' # API credentials -#' API_User = API_User, -#' API_Key = API_Key +#' Variable = "2m_temperature", +#' DataSet = "reanalysis-era5-land-monthly-means", +#' Type = "monthly_averaged_reanalysis", +#' # time-window, default set to range of dataset-type +#' DateStart = "1995-01-01 00:00", +#' DateStop = "1995-01-01 23:00", +#' TZone = "CET", +#' # temporal aggregation +#' TResolution = "month", +#' TStep = 1, +#' # file storing +#' FileName = "RawGlobe", +#' # API credentials +#' API_User = API_User, +#' API_Key = API_Key #' ) -#' terra::plot(RawGlobe_rast) +#' Plot.SpatRast(RawGlobe_rast) #' #' ## Monthly air temperature aggregated to bi-annual maximum by SpatRaster -#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package="KrigR")) +#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) #' BiAnnAirTemp_rast <- CDownloadS( -#' Variable = "2m_temperature", -#' DataSet = "reanalysis-era5-land-monthly-means", -#' Type = "monthly_averaged_reanalysis", -#' # time-window, default set to range of dataset-type -#' DateStart = "1995-01-01 00:00", -#' DateStop = "1996-12-31 23:00", -#' TZone = "EET", -#' # temporal aggregation -#' TResolution = "year", -#' TStep = 2, -#' # spatial -#' Extent = CDS_rast, -#' # file storing -#' FileName = "BiAnnAirTemp", -#' # API credentials -#' API_User = API_User, -#' API_Key = API_Key +#' Variable = "2m_temperature", +#' DataSet = "reanalysis-era5-land-monthly-means", +#' Type = "monthly_averaged_reanalysis", +#' # time-window, default set to range of dataset-type +#' DateStart = "1995-01-01 00:00", +#' DateStop = "1996-12-31 23:00", +#' TZone = "EET", +#' # temporal aggregation +#' TResolution = "year", +#' TStep = 2, +#' # spatial +#' Extent = CDS_rast, +#' # file storing +#' FileName = "BiAnnAirTemp", +#' # API credentials +#' API_User = API_User, +#' API_Key = API_Key #' ) -#' terra::plot(BiAnnAirTemp_rast) +#' Plot.SpatRast(BiAnnAirTemp_rast) #' #' ## Hourly back-calculated precipitation aggregated to daily averages by shapefiles #' data("Jotunheimen_poly") @@ -114,122 +114,142 @@ #' API_User = API_User, #' API_Key = API_Key #' ) -#' terra::plot(DailyBackCPrecip_rast) +#' Plot.SpatRast(DailyBackCPrecip_rast, SF = Jotunheimen_poly, Legend = "Precipitation [m]", COL = rev(viridis::cividis(100))) #' #' ## 6-hourly ensemble member spread sum for air temperature by buffered points #' data("Mountains_df") #' EnsembleSpreadSum6hour_rast <- CDownloadS( -#' Variable = "2m_temperature", -#' DataSet = "reanalysis-era5-single-levels", -#' Type = "ensemble_spread", -#' # time-window, default set to range of dataset-type -#' DateStart = "1995-01-01 00:00:00", -#' DateStop = "1995-01-01 21:00:00", -#' TZone = "UTC", -#' # temporal aggregation -#' TResolution = "hour", -#' TStep = 6, -#' FUN = sum, -#' # spatial -#' Extent = Mountains_df, -#' Buffer = 0.2, -#' # file storing -#' FileName = "EnsembleSpreadSum6hour", -#' FileExtension = ".tif", -#' # API credentials -#' API_User = API_User, -#' API_Key = API_Key, -#' Keep_Raw = TRUE +#' Variable = "2m_temperature", +#' DataSet = "reanalysis-era5-single-levels", +#' Type = "ensemble_spread", +#' # time-window, default set to range of dataset-type +#' DateStart = "1995-01-01 00:00:00", +#' DateStop = "1995-01-01 21:00:00", +#' TZone = "UTC", +#' # temporal aggregation +#' TResolution = "hour", +#' TStep = 6, +#' FUN = sum, +#' # spatial +#' Extent = Mountains_df, +#' Buffer = 0.2, +#' # file storing +#' FileName = "EnsembleSpreadSum6hour", +#' FileExtension = ".tif", +#' # API credentials +#' API_User = API_User, +#' API_Key = API_Key, +#' Keep_Raw = TRUE #' ) -#' terra::plot(EnsembleSpreadSum6hour_rast) +#' Plot.SpatRast(EnsembleSpreadSum6hour_rast, Legend = "Air Temperature Uncertainty [K]") #' } #' @export -CDownloadS <- function(Variable = NULL, # which variable +CDownloadS <- function(Variable = NULL, # which variable # nolint: cyclocomp_linter. CumulVar = FALSE, # cumulative variable? DataSet = "reanalysis-era5-land", # data set Type = NA, # type of data set DateStart, DateStop, TZone = "UTC", # time-window, default set to range of dataset-type - TResolution = "month", TStep = 1, FUN = 'mean', # temporal aggregation + TResolution = "month", TStep = 1, FUN = "mean", # temporal aggregation Extent, # spatial limitation, default set to range of dataset-type Buffer = 0.5, # point buffering if desired - Dir = getwd(), FileName, FileExtension = ".nc", # file storing + Dir = getwd(), FileName, FileExtension = ".nc", Compression = 9, # file storing API_User, API_Key, # API credentials TryDown = 10, TimeOut = 36000, # Calls to CDS TChunkSize = 6000, Cores = 1, # parallelisation verbose = TRUE, # verbosity - Keep_Raw = FALSE, - Save_Final = TRUE - ){ + Keep_Raw = FALSE) { ## Catching Most Frequent Issues =============== + on.exit(closeAllConnections()) #--- API Credentials ### checking if API User and Key have been supplied - if(exists("API_User") + exists("API_Key") != 2){ + if (exists("API_User") + exists("API_Key") != 2) { stop("Please provide a value for the API_User and API_Key arguments.") } ### making API_User into a character string API_User <- as.character(API_User) #--- Data set & Type ### checking if a supported data set has been queried - if(!(DataSet %in% Meta.List())){ - stop("Please specify a supported dataset as the DataSet argument. Your options are:", - "\n", paste(Meta.List(), collapse = (" \n"))) - } - if(!(Type %in% Meta.QuickFacts(dataset = DataSet)$Type)){ - stop("Please specify a Type argument that is supported by your chosen data set. Your options are:", - "\n", paste(Meta.QuickFacts(dataset = DataSet)$Type, collapse = (" \n")), - "\n !! If you are seeing an NA on the above line, note that this is not an error. Please specify NA as the Type.") + if (!(DataSet %in% Meta.List())) { + stop( + "Please specify a supported dataset as the DataSet argument. Your options are:", + "\n", paste(Meta.List(), collapse = (" \n")) + ) + } + if (!(Type %in% Meta.QuickFacts(dataset = DataSet)$Type)) { + stop( + "Please specify a Type argument that is supported by your chosen data set. Your options are:", + "\n", paste(Meta.QuickFacts(dataset = DataSet)$Type, collapse = (" \n")), + "\n !! If you are seeing an NA on the above line, note that this is not an error. Please specify NA as the Type." + ) } #--- File Name and Extension ### check if file name has been specified - if(!exists("FileName")){stop("Please provide a value for the FileName argument.")} + if (!exists("FileName")) { + stop("Please provide a value for the FileName argument.") + } FileName <- paste0(file_path_sans_ext(FileName), FileExtension) - if(!(FileExtension %in% c(".nc", ".tif"))){stop("Please specify a FileExtension of either '.tif' or '.nc'")} + if (!(FileExtension %in% c(".nc", ".tif"))) { + stop("Please specify a FileExtension of either '.tif' or '.nc'") + } #--- Time Zone - if(!(TZone %in% OlsonNames())){stop("The TZone argument you have specified is not supported. Please refer to OlsonNames() for an overview of all supported specifications.")} + if (!(TZone %in% OlsonNames())) { + stop("The TZone argument you have specified is not supported. Please refer to OlsonNames() for an overview of all supported specifications.") + } ## The Request ================================= - if(verbose){message("###### CDS Request & Data Download")} + if (verbose) { + message("###### CDS Request & Data Download") + } ### Building ===== - if(verbose){print("Building request")} + if (verbose) { + print("Building request") + } #--- Name resolving - VarPos <- which(rowSums(Meta.Variables(dataset = DataSet)[,1:2] == Variable) != 0) - QueryVariable <- Meta.Variables(dataset = DataSet)[VarPos,"CDSname"] + VarPos <- which(rowSums(Meta.Variables(dataset = DataSet)[, 1:2] == Variable) != 0) + QueryVariable <- Meta.Variables(dataset = DataSet)[VarPos, "CDSname"] #--- Extent resolving; formatting as SpatExtent object - if(missing(Extent)){Extent <- ext(Meta.QuickFacts(dataset = DataSet)$CDSArguments$area)} ## assign maximum extent for dataset if not specified - if(class(Extent)[1] == "data.frame"){ - Extent <- Buffer.pts(USER_pts = Make.SpatialPoints(USER_df = Extent), - USER_buffer = Buffer) + if (missing(Extent)) { + Extent <- ext(Meta.QuickFacts(dataset = DataSet)$CDSArguments$area) + } ## assign maximum extent for dataset if not specified + if (class(Extent)[1] == "data.frame") { + Extent <- Buffer.pts( + USER_pts = Make.SpatialPoints(USER_df = Extent), + USER_buffer = Buffer + ) } QuerySpace <- Ext.Check(Extent) - QueryExtent <- QuerySpace$Ext[c(4,1,3,2)] #N,W,S,E + QueryExtent <- QuerySpace$Ext[c(4, 1, 3, 2)] # N,W,S,E Extent <- QuerySpace$SpatialObj # terra/sf version of input extent to be used for easy cropping and masking #--- Base Dataset Information BaseResolution <- Meta.QuickFacts(dataset = DataSet)$TResolution BaseStep <- Meta.QuickFacts(dataset = DataSet)$TStep[ - na.omit(match(Type, Meta.QuickFacts(dataset = DataSet)$Type))] + na.omit(match(Type, Meta.QuickFacts(dataset = DataSet)$Type)) + ] BaseStart <- Meta.QuickFacts(dataset = DataSet)$TStart - if(BaseResolution == "hour" & CumulVar){ - DateStop <- as.character(as.POSIXct(DateStop, tz = TZone)+1*60*60*24) # add one day to hourly pulls when cumulVar is turned on as an extra layer of data is needed for proper backcalculation + if (BaseResolution == "hour" && CumulVar) { + DateStopIn <- as.POSIXct(DateStop, tz = TZone) + DateStop <- as.character(as.POSIXct(DateStop, tz = TZone) + 1 * 60 * 60 * 24) # add one day to hourly pulls when cumulVar is turned on as an extra layer of data is needed for proper backcalculation } #--- Time windows Dates_df <- Make.UTC(DatesVec = c( as.POSIXct(DateStart, tz = TZone), - as.POSIXct(DateStop, tz = TZone)) - ) - QueryTimeWindows <- Make.RequestWindows(Dates_df = Dates_df, - BaseTResolution = BaseResolution, - BaseTStep = 24/BaseStep, - BaseTStart = BaseStart, - TChunkSize = TChunkSize, - DataSet = DataSet + as.POSIXct(DateStop, tz = TZone) + )) + QueryTimeWindows <- Make.RequestWindows( + Dates_df = Dates_df, + BaseTResolution = BaseResolution, + BaseTStep = 24 / BaseStep, + BaseTStart = BaseStart, + TChunkSize = TChunkSize, + DataSet = DataSet ) QueryTimes <- QueryTimeWindows$QueryTimes QueryTimeWindows <- QueryTimeWindows$QueryTimeWindows @@ -246,23 +266,25 @@ CDownloadS <- function(Variable = NULL, # which variable ) ### Checking ===== - if(verbose){print("Checking request validity")} + if (verbose) { + print("Checking request validity") + } #--- Metadata check - can the queried dataset-type deliver the queried data? - MetaCheck_ls <- Meta.Check(DataSet = DataSet, - Type = Type, - VariableCheck = QueryVariable, - CumulativeCheck = CumulVar, - ExtentCheck = QueryExtent, - DateCheck = Dates_df, - AggrCheck = list(TStep, TResolution), - QueryTimes = QueryTimes) + MetaCheck_ls <- Meta.Check( + DataSet = DataSet, + Type = Type, + VariableCheck = QueryVariable, + CumulativeCheck = CumulVar, + ExtentCheck = QueryExtent, + DateCheck = Dates_df, + AggrCheck = list(TStep, TResolution), + QueryTimes = QueryTimes + ) ## File/Call metadata - KrigRCall <- match.call() - KrigRCall <- KrigRCall[!(names(KrigRCall) %in% c("API_Key", "API_User"))] - Meta_vec <- as.character(KrigRCall) - names(Meta_vec) <- names(KrigRCall) + Meta_vec <- c(Variable, CumulVar, DataSet, Type, DateStart, DateStop, TZone, TResolution, TStep, as.character(substitute(FUN)), as.character(quote(Extent)), Buffer, Dir, FileName) + names(Meta_vec) <- c("Variable", "CumulVar", "DataSet", "Type", "DateStart", "DateStop", "TZone", "TResolution", "TStep", "FUN", "Extent", "Buffer", "Dir", "FileName") Meta_vec <- c( "Citation" = paste0(MetaCheck_ls$QueryDataSet, " data (DOI:", Meta.DOI("reanalysis-era5-land-monthly-means"), ") obtained with KrigR (DOI:10.1088/1748-9326/ac48b3) on ", Sys.time()), "KrigRCall" = Meta_vec @@ -270,8 +292,8 @@ CDownloadS <- function(Variable = NULL, # which variable #--- File check, if already a file with this name present then load from disk FCheck <- Check.File(FName = FileName, Dir = Dir, loadFun = terra::rast, load = TRUE, verbose = TRUE) - if(!is.null(FCheck)){ - if(FileExtension == ".nc"){ + if (!is.null(FCheck)) { + if (FileExtension == ".nc") { FCheck <- Meta.NC(NC = FCheck, FName = file.path(Dir, FileName), Attrs = Meta_vec, Read = TRUE) } terra::time(FCheck) <- as.POSIXct(terra::time(FCheck), tz = TZone) # assign the correct time zone, when loading from disk, time zone is set to UTC @@ -279,53 +301,73 @@ CDownloadS <- function(Variable = NULL, # which variable } ### Executing ===== - if(verbose){print("Executing request - Notice that time windows may vary slightly at this step due to timezone conversions. This will be resolved automatically.")} + if (verbose) { + print("Executing request - Notice that time windows may vary slightly at this step due to timezone conversions. This will be resolved automatically.") + } #--- API credentials Register.Credentials(API_User, API_Key) #--- Make list of CDS Requests Requests_ls <- Make.Request(QueryTimeWindows, - MetaCheck_ls$QueryDataSet, MetaCheck_ls$QueryType, MetaCheck_ls$QueryVariable, - QueryTimes, QueryExtent, MetaCheck_ls$QueryFormat, - Dir, verbose = TRUE, API_User, API_Key) + MetaCheck_ls$QueryDataSet, MetaCheck_ls$QueryType, MetaCheck_ls$QueryVariable, + QueryTimes, QueryExtent, MetaCheck_ls$QueryFormat, + Dir, + verbose = TRUE, API_User, API_Key, + TimeOut = TimeOut + ) ## work an on.exit in here to allow restarting downloads themselves without new queries #--- Execution of requests Execute.Requests(Requests_ls, Dir, API_User, API_Key, TryDown, verbose = TRUE) ## The Data ================================= - if(verbose){message("###### Data Checking & Limitting & Aggregating")} + if (verbose) { + message("###### Data Checking & Limitting & Aggregating") + } TempFs <- file.path(Dir, unlist(lapply(strsplit(names(Requests_ls), " "), "[[", 2))) ## Checking ===== - if(verbose){print("Checking for known data issues")} + if (verbose) { + print("Checking for known data issues") + } #--- layers - NLyrCheck <- unlist(lapply(TempFs, FUN = function(LayerCheckIter){ + NLyrCheck <- unlist(lapply(TempFs, FUN = function(LayerCheckIter) { nlyr(rast(LayerCheckIter)) })) NLyrIssue <- which(NLyrCheck != unlist(lapply(QueryTimeWindows, length))) - if(length(NLyrIssue) > 0){ + if (length(NLyrIssue) > 0) { stop("Download of ", paste(basename(TempFs[NLyrIssue]), collapse = ", "), " produced file(s) of incorrect amount of layers. You may want to delete these files and try again. If the error persists. Please consult your queue on CDS: https://cds.climate.copernicus.eu/cdsapp#!/yourrequests. Alternatively, you may want to consult the corresponding download query/queries used behind the scenes:", paste(capture.output(str(Requests_ls[NLyrIssue])), collapse = "\n")) } #--- Loading data CDS_rast <- rast(TempFs) terra::time(CDS_rast) <- as.POSIXct(terra::time(CDS_rast), tz = TZone) # assign time in queried timezone - ## subset to desired time - CDS_rast <- CDS_rast[[which(!(terra::time(CDS_rast) < Dates_df$IN[1] | terra::time(CDS_rast) > Dates_df$IN[2]))]] ## Spatial ===== - if(verbose){print("Spatial Limiting")} + if (verbose) { + print("Spatial Limiting") + } CDS_rast <- Handle.Spatial(CDS_rast, Extent) ## Temporal ===== - if(verbose){print("Temporal Aggregation")} #--- Cumulative Fix - CDS_rast <- Temporal.Cumul(CDS_rast, CumulVar, BaseResolution, BaseStep, TZone) + CDS_rast <- Temporal.Cumul(CDS_rast, CumulVar, BaseResolution, BaseStep, TZone, verbose) + + #--- Subset to desired time, happens here to allow for correct disaggregation of cumulative variables in previous step + terra::time(CDS_rast) <- as.POSIXct(terra::time(CDS_rast), tz = TZone) # assign time in queried timezone + if (exists("DateStopIn")) { + Dates_df$IN[2] <- DateStopIn + } + CDS_rast <- CDS_rast[[which((terra::time(CDS_rast) < Dates_df$IN[1]) + (terra::time(CDS_rast) > Dates_df$IN[2]) == 0)]] + terra::time(CDS_rast) <- as.POSIXct(terra::time(CDS_rast), tz = TZone) # assign time in queried timezone #--- Temporal aggregation - CDS_rast <- Temporal.Aggr(CDS_rast, BaseResolution, BaseStep, - TResolution, TStep, FUN, Cores, QueryTargetSteps, TZone) + CDS_rast <- Temporal.Aggr( + CDS_rast, BaseResolution, BaseStep, + TResolution, TStep, FUN, Cores, QueryTargetSteps, TZone, verbose + ) ## Exports ================================= - if(verbose){message("###### Data Export & Return")} + if (verbose) { + message("###### Data Export & Return") + } ### Assign additional information terra::varnames(CDS_rast) <- MetaCheck_ls$QueryVariable @@ -333,21 +375,25 @@ CDownloadS <- function(Variable = NULL, # which variable terra::metags(CDS_rast) <- Meta_vec ### write file - if(Save_Final){ - if(FileExtension == ".tif"){ - terra::writeRaster(CDS_rast, filename = file.path(Dir, FileName)) - } - if(FileExtension == ".nc"){ - CDS_rast <- Meta.NC(NC = CDS_rast, FName = file.path(Dir, FileName), - Attrs = terra::metags(CDS_rast), Write = TRUE) - } + if (FileExtension == ".tif") { + terra::writeRaster(CDS_rast, filename = file.path(Dir, FileName)) + CDS_rast <- terra::rast(filename = file.path(Dir, FileName)) + } + if (FileExtension == ".nc") { + CDS_rast <- Meta.NC( + NC = CDS_rast, FName = file.path(Dir, FileName), + Attrs = Meta_vec, Write = TRUE, + Compression = Compression + ) + terra::time(CDS_rast) <- as.POSIXct(terra::time(CDS_rast), tz = TZone) # assign the correct time zone, when loading from disk, time zone is set to UTC } ### unlink temporary files - if(!Keep_Raw){ + if (!Keep_Raw) { unlink(TempFs) } ### return object + closeAllConnections() return(CDS_rast) } diff --git a/R/Checks.R b/R/Checks.R index 0213ba7..50556d5 100644 --- a/R/Checks.R +++ b/R/Checks.R @@ -13,23 +13,29 @@ #' #' @examples #' KrigR::Check.File( -#' FName = basename(system.file("extdata", "CentralNorway.nc", package="KrigR")), -#' Dir = dirname(system.file("extdata", "CentralNorway.nc", package="KrigR")), -#' loadFun = terra::rast -#' ) +#' FName = basename(system.file("extdata", "CentralNorway.nc", package = "KrigR")), +#' Dir = dirname(system.file("extdata", "CentralNorway.nc", package = "KrigR")), +#' loadFun = terra::rast +#' ) #' @export -Check.File <- function(FName, Dir = getwd(), loadFun, load = TRUE, verbose = TRUE){ +Check.File <- function(FName, Dir = getwd(), loadFun, load = TRUE, verbose = TRUE) { FNAME <- file.path(Dir, FName) file <- NULL - if(file.exists(FNAME)){ - if(verbose){print(paste0("A file with the name ", FName, " already exists in ", Dir, - "."))} - if(load){ - if(verbose){print("Loading this file for you from the disk.")} + if (file.exists(FNAME)) { + if (verbose) { + print(paste0( + "A file with the name ", FName, " already exists in ", Dir, + "." + )) + } + if (load) { + if (verbose) { + print("Loading this file for you from the disk.") + } file <- sapply(FNAME, loadFun)[[1]] - }else{ + } else { file <- "Present. Not Loaded." - } + } } return(file) } @@ -54,65 +60,65 @@ Check.File <- function(FName, Dir = getwd(), loadFun, load = TRUE, verbose = TRU #' #' @seealso \code{\link{Kriging}}, \code{\link{KrigingCovariateSetup}}. #' -Check.Krig <- function(Data, CovariatesCoarse, CovariatesFine, KrigingEquation){ +Check.Krig <- function(Data, CovariatesCoarse, CovariatesFine, KrigingEquation) { ## Resolutions =============== - if(terra::res(CovariatesFine)[1] < terra::res(Data)[1]/10){ + if (terra::res(CovariatesFine)[1] < terra::res(Data)[1] / 10) { warning("It is not recommended to use kriging for statistical downscaling of more than one order of magnitude. You are currently attempting this. Kriging will proceed.") } - if(all.equal(terra::res(CovariatesCoarse)[1], terra::res(Data)[1]) != TRUE){ - stop(paste0("The resolution of your data (", terra::res(Data)[1], ") does not match the resolution of your covariate data (", terra::res(CovariatesCoarse)[1], ") used for training the kriging model. Kriging can't be performed!" )) + if (all.equal(terra::res(CovariatesCoarse)[1], terra::res(Data)[1]) != TRUE) { + stop(paste0("The resolution of your data (", terra::res(Data)[1], ") does not match the resolution of your covariate data (", terra::res(CovariatesCoarse)[1], ") used for training the kriging model. Kriging can't be performed!")) } ## Extent =============== - if(terra::ext(Data) == terra::ext(-180, 180, -90, 90)){ + if (terra::ext(Data) == terra::ext(-180, 180, -90, 90)) { stop("You are attempting to use kriging at a global extent. For reasons of computational expense and identity of relationships between covariates and variables not being homogenous across the globe, this is not recommended. Instead, try kriging of latitude bands if global kriging is really your goal.") } - if(!all.equal(terra::ext(CovariatesCoarse), terra::ext(Data))){ + if (!all.equal(terra::ext(CovariatesCoarse), terra::ext(Data))) { stop("The extents of your data and training covariates don't match. Kriging can't be performed!") } ## Data Availability =============== DataSkips <- NULL # data layers without enough data to be skipped in kriging Data_vals <- base::colSums(matrix(!is.na(terra::values(Data)), ncol = terra::nlyr(Data))) # a value of 0 indicates a layer only made of NAs - if(length(which(Data_vals < 2)) > 0){ - if(length(which(Data_vals < 2)) != terra::nlyr(Data)){ - stop(paste0("Layer(s) ", paste(which(Data_vals == 0), collapse=", "), " of your data do(es) not contain enough data. Kriging cannot be performed. Usually, increasing the extent of kriging can fix this issue.")) + if (length(which(Data_vals < 2)) > 0) { + if (length(which(Data_vals < 2)) != terra::nlyr(Data)) { + stop(paste0("Layer(s) ", paste(which(Data_vals == 0), collapse = ", "), " of your data do(es) not contain enough data. Kriging cannot be performed. Usually, increasing the extent of kriging can fix this issue.")) DataSkips <- which(Data_vals < 2) - }else{ + } else { stop("Your Data does not contain enough values. Kriging cannot be performed. Usually, increasing the extent of kriging can fix this issue.") } } CovCo_vals <- base::colSums(matrix(!is.na(terra::values(CovariatesCoarse)), ncol = terra::nlyr(CovariatesCoarse))) # a value of 0 indicates a layer only made of NAs - if(length(which(CovCo_vals < 2)) > 0){ - if(length(which(CovCo_vals < 2)) != terra::nlyr(CovariatesCoarse)){ - warning(paste0("Layer(s) ", paste(which(CovCo_vals < 2), collapse=", "), " of your covariates at training resolution do(es) not contain enough data. This/these layer(s) is/are dropped. The Kriging equation might get altered.")) + if (length(which(CovCo_vals < 2)) > 0) { + if (length(which(CovCo_vals < 2)) != terra::nlyr(CovariatesCoarse)) { + warning(paste0("Layer(s) ", paste(which(CovCo_vals < 2), collapse = ", "), " of your covariates at training resolution do(es) not contain enough data. This/these layer(s) is/are dropped. The Kriging equation might get altered.")) CovariatesCoarse <- CovariatesCoarse[[-which(CovCo_vals < 2)]] - }else{ + } else { stop("Your covariate data at training resolution does not contain enough values. Kriging can't be performed!") } } CovFin_vals <- base::colSums(matrix(!is.na(terra::values(CovariatesFine)), ncol = terra::nlyr(CovariatesFine))) # a value of 0 indicates a layer only made of NAs - if(length(which(CovFin_vals < 2)) > 0){ - if(length(which(CovFin_vals < 2)) != terra::nlyr(CovariatesFine)){ - warning(paste0("Layer(s) ", paste(which(CovFin_vals == 0), collapse=", "), " of your covariates at target resolution do(es) not contain enough data. This/these layer(s) is/are dropped.")) + if (length(which(CovFin_vals < 2)) > 0) { + if (length(which(CovFin_vals < 2)) != terra::nlyr(CovariatesFine)) { + warning(paste0("Layer(s) ", paste(which(CovFin_vals == 0), collapse = ", "), " of your covariates at target resolution do(es) not contain enough data. This/these layer(s) is/are dropped.")) CovariatesFine <- CovariatesFine[[-which(CovFin_vals < 2)]] - }else{ + } else { stop("Your covariate data at target resolution does not contain enough values. Kriging can't be performed!") } } ## Kriging Equation =============== Terms <- unlist(strsplit(labels(terms(KrigingEquation)), split = ":")) # identify parameters called to in formula Terms_Required <- unique(Terms) # isolate double-references (e.g. due to ":" indexing for interactions) - Terms_Present <- Reduce(intersect, list(Terms_Required, terra::varnames(CovariatesCoarse), terra::varnames(CovariatesFine))) # identify the terms that are available and required - if(sum(Terms_Required %in% Terms_Present) != length(Terms_Required)){ - if(length(Terms_Present) == 0){ # if none of the specified terms were found - KrigingEquation <- paste0("Data ~ ", paste(terra::varnames(CovariatesCoarse), collapse = "+")) + Terms_Present <- Reduce(intersect, list(Terms_Required, names(CovariatesCoarse), names(CovariatesFine))) # identify the terms that are available and required + if (sum(Terms_Required %in% Terms_Present) != length(Terms_Required)) { + if (length(Terms_Present) == 0) { # if none of the specified terms were found + KrigingEquation <- paste0("Data ~ ", paste(names(CovariatesCoarse), collapse = "+")) warn <- paste("None of the terms specified in your KrigingEquation are present in the covariate data sets. The KrigingEquation has been altered to include all available terms in a linear model:", "\n\n", KrigingEquation) - }else{ # at least some of the specified terms were found + } else { # at least some of the specified terms were found KrigingEquation <- paste0("Data ~ ", paste(Terms_Present, collapse = "+")) warn <- paste("Not all of the terms specified in your KrigingEquation are present in the covariate data sets. The KrigingEquation has been altered to include all available and specified terms in a linear model:", "\n\n", KrigingEquation) } - Cotinue <- menu(c("Yes", "No"), title=paste0(warn, ". \n\nDo you wish to continue using the new formula?")) - if(Cotinue == 2){ # break operation if user doesn't want this + Cotinue <- menu(c("Yes", "No"), title = paste0(warn, ". \n\nDo you wish to continue using the new formula?")) + if (Cotinue == 2) { # break operation if user doesn't want this stop("Kriging terminated by user due to formula issues.") } } diff --git a/R/CovariateSetup.R b/R/CovariateSetup.R index e84c4ef..b25af21 100644 --- a/R/CovariateSetup.R +++ b/R/CovariateSetup.R @@ -1,16 +1,18 @@ #' Preparing Covariate Data for Use in Kriging #' -#' This function is used to setup products of covariate data ready for use in Kriging. This functiuonality can either be applied to user-supplied covariate data or ready-made data products such as the Harmised World Soil Data Base and the median statistic of the Global Multi-resolution Terrain Elevation Data (GMTED2010; available at \url{https://topotools.cr.usgs.gov/gmted_viewer/gmted2010_global_grids.php}). In case of the latter, the data is downloaded at 30 arc-sec latitude/longitude grid cells and subsequently resampled to match training and target resolutions specified by the user. +#' This function is used to setup products of covariate data ready for use in Kriging. This functionality can either be applied to user-supplied covariate data or ready-made data products such as the global dataset of soil hydraulic and thermal parameters for earth system modeling (available at \url{http://globalchange.bnu.edu.cn/research/soil4.jsp}) and the median statistic of the Global Multi-resolution Terrain Elevation Data (GMTED2010; available at \url{https://topotools.cr.usgs.gov/gmted_viewer/gmted2010_global_grids.php}). In case of the latter, the data is downloaded at 30 arc-sec latitude/longitude grid cells and subsequently resampled to match training and target resolutions specified by the user. #' #' @param Training A SpatRaster file containing the data which is to be downscaled. Covariate data will be resampled to match this. #' @param Target Either numeric or a SpatRaster. If numeric, a single number representing the target resolution for the kriging step (i.e. wich resolution to downscale to). If a SpatRaster, data that the covariates and kriged products should align with. In case of a numeric input, covariate data is aggregated as closely as possible to desired resolution. If a SpatRaster, covariate data is resampled to match desired output directly. -#' @param Covariates Either character or a SpatRaster. If character, obtain frequently used and provably useful covariate data (i.e., GMTED2010 and HWSD) and prepare for use in Kriging. Supported character values are "GMTED2010" and "HWSD". Note that currently, HWSD data download is not functional. If a SpatRaster, a user-supplied set of covariate data to be prepared for use in Kriging. +#' @param FilePrefix Character. A file name prefix for the produced files. +#' @param Covariates Either character or a SpatRaster. If character, obtain frequently used and provably useful covariate data (i.e., GMTED2010 and soil data) and prepare for use in Kriging. Supported character values are "GMTED2010", "tksat", "tkdry", "csol", "k_s", "lambda", "psi", and "theta_s". If a SpatRaster, a user-supplied set of covariate data to be prepared for use in Kriging. #' @param Source Character. Only comes into effect when Covariates argument is specified as a character. Whether to attempt download of covariate data from the official sources (Source = "Origin") or a static copy of the data set on a private drive (Source = "Drive"). Default is "Origin". #' @param Extent Optional, prepare covariate data according to desired spatial specification. If missing/unspecified, maximal area of supplied data and covariat sets is used. Can be specified either as a raster object, an sf object, a terra object, or a data.frame. If Extent is a raster or terra object, covariates will be prepared according to rectangular extent thereof. If Extent is an sf (MULTI-)POLYGON object, this will be treated as a shapefile and the output will be cropped and masked to this shapefile. If Extent is a data.frame of geo-referenced point records, it needs to contain Lat and Lon columns around which a buffered shapefile will be created using the Buffer argument. #' @param Buffer Optional, Numeric. Identifies how big a circular buffer to draw around points if Extent is a data.frame of points. Buffer is expressed as centessimal degrees. #' @param Dir Character/Directory Pointer. Directory specifying where to download data to. #' @param Keep_Global Logical. Only comes into effect when Covariates argument is specified as a character. Whether to retain raw downloaded covariate data or not. Default is FALSE. #' @param FileExtension Character. A file extension for the produced files. Supported values are ".nc" (default) and ".tif" (better support for metadata). +#' @param Compression Integer between 1 to 9. Applied to final .nc file that the function writes to hard drive. Same as compression argument in terra::writeCDF(). Ignored if FileExtension = ".tif". #' #' @importFrom httr GET #' @importFrom httr write_disk @@ -25,6 +27,8 @@ #' @importFrom terra writeRaster #' @importFrom terra writeCDF #' @importFrom terra varnames +#' @importFrom terra sources +#' @importFrom tools file_path_sans_ext #' #' @return A list containing two SpatRaster objects (Training and Target) ready to be used as covariates for kriging, and two files called Covariates_Target and Covariates_Train in the specified directory. #' @@ -33,98 +37,127 @@ #' \item{Citation}{ - A string which to use for in-line citation of the data product.} #' } #' -#' @seealso \code{\link{Kriging}}. +#' @seealso \code{\link{Kriging}}, \code{\link{Plot.Covariates}}. #' #' @examples #' \dontrun{ #' ## Rectangular Covariate data according to input data -#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package="KrigR")) -#' Covariates_ls <- CovariateSetup(Training = CDS_rast, -#' Target = 0.01, -#' Covariates = "GMTED2010", -#' Keep_Global = TRUE, -#' FileExtension = ".nc") -#' terra::plot(Covariates_ls[[1]]) -#' terra::plot(Covariates_ls[[2]]) +#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +#' Covariates_ls <- CovariateSetup( +#' Training = CDS_rast, +#' Target = 0.01, +#' Covariates = "GMTED2010", +#' Keep_Global = TRUE, +#' FileExtension = ".nc" +#' ) +#' Plot.Covariates(Covariates_ls) #' #' ## Shapefile-limited covariate data #' data("Jotunheimen_poly") -#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package="KrigR")) -#' Covariates_ls <- CovariateSetup(Training = CDS_rast, -#' Target = 0.01, -#' Covariates = "GMTED2010", -#' Extent = Jotunheimen_poly, -#' Keep_Global = TRUE, -#' FileExtension = ".nc") -#' terra::plot(Covariates_ls[[1]]) -#' terra::plot(Covariates_ls[[2]]) +#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +#' Covariates_ls <- CovariateSetup( +#' Training = CDS_rast, +#' Target = 0.01, +#' Covariates = "GMTED2010", +#' Extent = Jotunheimen_poly, +#' Keep_Global = TRUE, +#' FileExtension = ".nc" +#' ) +#' Plot.Covariates(Covariates_ls, SF = Jotunheimen_poly) #' #' ## buffered-point-limited covariate data #' data("Mountains_df") -#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package="KrigR")) -#' Covariates_ls <- CovariateSetup(Training = CDS_rast, -#' Target = 0.01, -#' Covariates = "GMTED2010", -#' Extent = Mountains_df, -#' Buffer = 0.2, -#' Keep_Global = TRUE, -#' FileExtension = ".nc") -#' terra::plot(Covariates_ls[[1]]) -#' terra::plot(Covariates_ls[[2]]) +#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +#' Covariates_ls <- CovariateSetup( +#' Training = CDS_rast, +#' Target = 0.01, +#' Covariates = c("tksat", "tkdry", "csol", "k_s", "lambda", "psi", "theta_s"), +#' Source = "Drive", +#' Extent = Mountains_df, +#' Buffer = 0.2, +#' Keep_Global = TRUE, +#' FileExtension = ".nc" +#' ) +#' Plot.Covariates(Covariates_ls) #' } #' @export -CovariateSetup <- function(Training, - Target, - Covariates = "GMTED2010", - Source = "Origin", - Extent, - Buffer = 0.5, - Dir = getwd(), - Keep_Global = FALSE, - FileExtension = ".nc" - ){ +CovariateSetup <- function(Training, # nolint: cyclocomp_linter. + Target, + FilePrefix = "", + Covariates = "GMTED2010", + Source = "Origin", + Extent, + Buffer = 0.5, + Dir = getwd(), + Keep_Global = FALSE, + FileExtension = ".nc", + Compression = 9) { ## Catching Most Frequent Issues =============== - if(class(Covariates) == "character"){ - if(sum(!(Covariates %in% c("GMTED2010", "HWSD"))) > 0){ - stop("Please specify a valid covariate data set. You may supply either the character string 'GMTED2010' or 'HWSD', or a SpatRaster object.") + if (class(Covariates) == "character") { + if (sum(!(Covariates %in% c("GMTED2010", "tksat", "tkdry", "csol", "k_s", "lambda", "psi", "theta_s"))) > 0) { + stop("Please specify a valid covariate data set. You may supply either a character string (allowed values are: 'GMTED2010', 'tksat', 'tkdry', 'csol', 'k_s', 'lambda', 'psi' and 'theta_s') or a SpatRaster object.") } - if("HWSD" %in% Covariates){ - stop("HWSD download currently not supported. We are working on it.") - } - if(length(Source) > 1){ + if (length(Source) > 1) { stop("Please specify only one source type for the covariate data that will be downloaded. You may specify either 'Origin' or 'Drive'.") } - if(!(Source %in% c("Origin", "Drive"))){ + if (!(Source %in% c("Origin", "Drive"))) { stop("Please specify a valid source type for the covariate data that will be downloaded. You may specify either 'Origin' or 'Drive'.") } } ## Links for Data download =============== Links_df <- data.frame( - Cov = rep(c("GMTED2010", "HWSD"), each = 2), - Source = rep(c("Origin", "Drive")), - UnzippedFile = rep(c("mn30_grd/w001001.adf", "NULL"), each = 2), - DOI = rep(c("10.3133/ofr20111073", "NULL"), each = 2), + Cov = rep(c( + "GMTED2010", + "tksat", "tkdry", "csol", "k_s", "lambda", "psi", "theta_s" + ), each = 2), + Source = rep(c("Origin", "Drive"), 8), + UnzippedFile = rep(c( + "mn30_grd/w001001.adf", + paste0(c("tksatu", "tkdry", "csol", "k_s", "lambda", "psi_s", "theta_s"), "_l1.nc") + ), each = 2), + DOI = rep(c("10.3133/ofr20111073", "10.1175/JHM-D-12-0149.1"), c(2, 14)), Link = c( "https://edcintl.cr.usgs.gov/downloads/sciweb1/shared/topo/downloads/GMTED/Grid_ZipFiles/mn30_grd.zip", # Link to GMTED2010 - "https://www.dropbox.com/s/whkje7jc401xuwx/GMTED2010.zip?raw=1", # Link to DropBox with GMTED2010 - NULL, # Link to Harmonized World Soil Database v2.0 - NULL # Link to DropBox with GMTED2010 + "https://www.dropbox.com/scl/fi/io099a5lrxawdc6v9wy0b/GMTED2010.zip?rlkey=izqwy2uhebvh4ypgaiwkiuk6l&st=br70b3jt&dl=1", # Link to DropBox with GMTED2010 + ## tksat + "http://globalchange.bnu.edu.cn/download/data/worldptf/tksat.zip", # Origin + "https://www.dropbox.com/scl/fi/k8ly97961yukgsw3dyzxy/tksat.zip?rlkey=v03z7vjx122nhj156xd3bnfyj&st=piieeigq&dl=1", # DropBox + ## tkdry + "http://globalchange.bnu.edu.cn/download/data/worldptf/tkdry.zip", # Origin + "https://www.dropbox.com/scl/fi/866sny1x7393a36naqfez/tkdry.zip?rlkey=pgqulsdoz2p3qznyvo5eqt7yw&st=s01p43us&dl=1", # DropBox + ## csol + "http://globalchange.bnu.edu.cn/download/data/worldptf/csol.zip", # Origin + "https://www.dropbox.com/scl/fi/a8cbzy8zn38zw6ou8l580/csol.zip?rlkey=3nv321henvysppav8l3ln5ja6&st=40sjwly9&dl=1", # DropBox + ## k_s + "http://globalchange.bnu.edu.cn/download/data/worldptf/k_s.zip", # Origin + "https://www.dropbox.com/scl/fi/87kak98esfe4b3srd8vwk/k_s.zip?rlkey=2ktmnl1rhakbx9xvroxb2emcz&st=goevwluk&dl=1", # DropBox + ## lambda + "http://globalchange.bnu.edu.cn/download/data/worldptf/lambda.zip", # Origin + "https://www.dropbox.com/scl/fi/uqp3g3pcak96qxkvnc1uw/lambda.zip?rlkey=tchnkjxpgf4wz0xpfmna7h8we&st=njbd097a&dl=1", # DropBox + ## psi + "http://globalchange.bnu.edu.cn/download/data/worldptf/psi.zip", # Origin + "https://www.dropbox.com/scl/fi/mrh4jql8c5nytli5n32um/psi.zip?rlkey=km89j9t3rum83o56d81qf8nyi&st=sqd7ua6v&dl=1", # DropBox + ## theta_s + "http://globalchange.bnu.edu.cn/download/data/worldptf/theta_s.zip", # Origin + "https://www.dropbox.com/scl/fi/9a2z1k0i91awqczgpi378/theta_s.zip?rlkey=zanege1p7tevp1b1k01vb277b&st=48neb5te&dl=1" # DropBox ) ) ## Data Download (skipped if own SpatRaster supplied) =============== - if(class(Covariates) == "character"){ + if (class(Covariates) == "character") { + message("###### Downloading global covariate data") CovariatesIn <- Covariates ### Directory for raw files Dir.Covs <- file.path(Dir, "CovariateSetup") - if(!dir.exists(Dir.Covs)){dir.create(Dir.Covs)} + if (!dir.exists(Dir.Covs)) { + dir.create(Dir.Covs) + } ### Data downloads - Data_ls <- lapply(Covariates, FUN = function(Cov_iter){ - + Data_ls <- lapply(Covariates, FUN = function(Cov_iter) { #### Figure out which link to use - Match_vec <- sapply(1:nrow(Links_df), FUN = function(x){ + Match_vec <- sapply(1:nrow(Links_df), FUN = function(x) { sum(Links_df$Cov[x] == Cov_iter, Links_df$Source[x] == Source) }) @@ -141,44 +174,73 @@ CovariateSetup <- function(Training, #### Check if data is already present Data <- Check.File(FName = FName, Dir = Dir, loadFun = terra::rast, load = TRUE, verbose = FALSE) - if(!is.null(Data) & FileExtension == ".nc"){ + if (!is.null(Data) & FileExtension == ".nc") { Data <- Meta.NC(NC = Data, FName = file.path(Dir, FName), Attrs = Meta_vec, Read = TRUE) } - if(is.null(Data)){ + if (is.null(Data)) { #### Downloading data - message("Downloading ", Name, " covariate data.") # inform user of download in console + print(paste(Name, "covariate data.")) # inform user of download in console httr::GET(Link, - httr::write_disk(file.path(Dir.Covs, paste0(Name, ".zip"))), - httr::progress(), overwrite = TRUE) + httr::write_disk(file.path(Dir.Covs, paste0(Name, ".zip"))), + httr::progress(), + overwrite = TRUE + ) #### Unzipping data - unzip(file.path(Dir.Covs, paste0(Name, ".zip")), # which file to unzip - exdir = Dir.Covs) # where to unzip to + if (Name == "GMTED2010") { + unzip(file.path(Dir.Covs, paste0(Name, ".zip")), # which file to unzip + exdir = Dir.Covs + ) # where to unzip to + } else { + unzip(file.path(Dir.Covs, paste0(Name, ".zip")), # which file to unzip + files = UnzippedFile, exdir = Dir.Covs + ) # where to unzip to + } + #### Loading data Data <- terra::rast(file.path(Dir.Covs, UnzippedFile)) - if(class(terra::values(Data)[,1]) == "integer"){ - print("Reformatting integer data into continuous numeric data.") - Data <- Data+0 # +0 to avoid integer reading in faulty way, https://gis.stackexchange.com/questions/398061/reading-rasters-in-r-using-terra-package + if (class(terra::values(Data)[, 1]) == "integer") { + print("Reformatting integer data into continuous numeric data. Necessary for GMTED2010 data.") + Data <- Data + 0 # +0 to avoid integer reading in faulty way, https://gis.stackexchange.com/questions/398061/reading-rasters-in-r-using-terra-package } terra::metags(Data) <- Meta_vec terra::varnames(Data) <- Name #### Saving data as single file - if(FileExtension == ".tif"){ + if (FileExtension == ".tif") { terra::writeRaster(Data, filename = file.path(Dir, FName)) } - if(FileExtension == ".nc"){ - Data <- Meta.NC(NC = Data, FName = file.path(Dir, FName), - Attrs = terra::metags(Data), Write = TRUE) + if (FileExtension == ".nc") { + Data <- Meta.NC( + NC = Data, FName = file.path(Dir, FName), + Attrs = terra::metags(Data), Write = TRUE, + Compression = Compression + ) } - }else{ - message("Raw ", Name, " covariate data already downloaded.") + Data <- Check.File(FName = FName, Dir = Dir, loadFun = terra::rast, load = TRUE, verbose = FALSE) + if (FileExtension == ".nc") { + Data <- Meta.NC(NC = Data, FName = file.path(Dir, FName), Attrs = Meta_vec, Read = TRUE) + } + } else { + print(paste(Name, "covariate data already downloaded.")) } Data }) + + if ("GMTED2010" %in% CovariatesIn && any(c("tksat", "tkdry", "csol", "k_s", "lambda", "psi", "theta_s") %in% CovariatesIn)) { + message("###### Aligning data from different sources with one another") + MinExt <- apply(abs(do.call(rbind, lapply(Data_ls, FUN = function(x) { + as.vector(ext(x)) + }))), 2, min) + Data_ls <- lapply(Data_ls, FUN = function(x) { + print(basename(tools::file_path_sans_ext(terra::sources(x)))) + crop(x, ext(-MinExt[1], MinExt[2], -MinExt[3], MinExt[4])) # I can simply to -/+ here because this data is global + }) + } + Covariates <- do.call(c, Data_ls) unlink(Dir.Covs, recursive = TRUE) } @@ -186,67 +248,75 @@ CovariateSetup <- function(Training, ## Spatial Limitting =============== ### Extent Handling - if(missing("Extent")){ ## assign maximum extent of supplied data and covariates (only when no extent is specified) - Extent <-terra::ext( + if (missing("Extent")) { ## assign maximum extent of supplied data and covariates (only when no extent is specified) + Extent <- terra::ext( ifelse(terra::ext(Training)[1] > terra::ext(Covariates)[1], terra::ext(Training)[1], terra::ext(Covariates)[1]), ifelse(terra::ext(Training)[2] < terra::ext(Covariates)[2], terra::ext(Training)[2], terra::ext(Covariates)[2]), ifelse(terra::ext(Training)[3] > terra::ext(Covariates)[3], terra::ext(Training)[3], terra::ext(Covariates)[3]), ifelse(terra::ext(Training)[4] < terra::ext(Covariates)[4], terra::ext(Training)[4], terra::ext(Covariates)[4]) ) - } - if(class(Extent)[1] == "data.frame"){ - Extent <- Buffer.pts(USER_pts = Make.SpatialPoints(USER_df = Extent), - USER_buffer = Buffer) + } + if (class(Extent)[1] == "data.frame") { + Extent <- Buffer.pts( + USER_pts = Make.SpatialPoints(USER_df = Extent), + USER_buffer = Buffer + ) } QuerySpace <- Ext.Check(Extent) Extent <- QuerySpace$SpatialObj # terra/sf version of input extent to be used for easy cropping and masking ## Spatial Aggregation/Resampling =============== + ### Cropping and Masking + Covariates <- Handle.Spatial(Covariates, ext(Extent)) + ### Sanity Check - if(class(Target) == "numeric"){ + if (class(Target) == "numeric") { Target_res <- Target[1] - }else{ + } else { Target_res <- terra::res(Target) } - if(Target_res < terra::res(Covariates)[1]){ - stop(paste0("You have specified resolution(s) to be finer than ", res(GMTED2010_ras), " (native GMTED2010 reslution). Please download higher-resolution DEM data instead.")) + if (Target_res < terra::res(Covariates)[1]) { + stop(paste0("You have specified resolution(s) to be finer than ", res(Covariates), " (native covariate reslution). Please provide higher-resolution data instead.")) } ### Resampling + message("###### Resampling Data") Cov_train <- terra::resample(Covariates, Training) - if(class(Extent)[1] == "SpatRaster"){ + if (class(Extent)[1] == "SpatRaster") { Cov_target <- terra::resample(Covariates, Extent) - }else{ - Cov_target <- suppressWarnings(terra::aggregate(Covariates, fact = Target_res[1]/terra::res(Covariates)[1])) + } else { + Cov_target <- suppressWarnings(terra::aggregate(Covariates, fact = Target_res[1] / terra::res(Covariates)[1])) } + ### Cropping and Masking - Training <- Handle.Spatial(BASE = Training, Shape = Extent) Cov_train <- Handle.Spatial(Cov_train, Extent) Cov_target <- Handle.Spatial(Cov_target, Extent) ## Data Saving & Export =============== - TrainName <- file.path(Dir, paste0("Covariates_Train", FileExtension)) - TargetName <- file.path(Dir, paste0("Covariates_Target", FileExtension)) - if(FileExtension == ".tif"){ + TrainName <- file.path(Dir, paste0(FilePrefix, "Covariates_Train", FileExtension)) + TargetName <- file.path(Dir, paste0(FilePrefix, "Covariates_Target", FileExtension)) + if (FileExtension == ".tif") { terra::writeRaster(x = Cov_train, filename = TrainName, overwrite = TRUE) terra::writeRaster(x = Cov_target, filename = TargetName, overwrite = TRUE) } - if(FileExtension == ".nc"){ - terra::writeCDF(x = Cov_train, filename = TrainName, overwrite = TRUE) - terra::writeCDF(x = Cov_target, filename = TargetName, overwrite = TRUE) + if (FileExtension == ".nc") { + terra::writeCDF(x = Cov_train, filename = TrainName, overwrite = TRUE, Compression = Compression) + terra::writeCDF(x = Cov_target, filename = TargetName, overwrite = TRUE, Compression = Compression) } ## Cleaning up files =============== - if(!Keep_Global){ # cleanup check - unlink(list.files(Dir, pattern = paste0(CovariatesIn, FileExtension), full.names = TRUE)) + if (!Keep_Global) { # cleanup check + unlink(file.path(Dir, paste0(CovariatesIn, FileExtension))) } ## Return data =============== TrainRet <- terra::rast(TrainName) TargetRet <- terra::rast(TargetName) - terra::varnames(TrainRet) <- terra::varnames(TargetRet) <- VarNames + names(TrainRet) <- names(TargetRet) <- VarNames return( - list(Training = TrainRet, - Target = TargetRet) - ) + list( + Training = TrainRet, + Target = TargetRet + ) + ) } diff --git a/R/download.R b/R/DEPRECATED.R similarity index 60% rename from R/download.R rename to R/DEPRECATED.R index 2004da9..96ee66d 100644 --- a/R/download.R +++ b/R/DEPRECATED.R @@ -57,7 +57,7 @@ download_ERA <- function(Variable = NULL, PrecipFix = FALSE, Type = "reanalysis" Cores = 1, TimeOut = 36000, SingularDL = FALSE, ...) { - stop("Function currently deprecated as KrigR undergoes major re-development. Please use the stable release to gain access to this functionality.") + stop("This function is deprecated and has been superceeded by the CDownloadS() function which offers improvements over this implementation. To get started using the new function, you may want to consult the relevant documentation by calling ?CDownloadS(). Please note that this deprecated function will be removed alltogether when KrigR reaches version 1.0.0.") if(verbose){message("download_ERA() is starting. Depending on your specifications, this can take a significant time.")} @@ -558,7 +558,7 @@ download_DEM <- function(Train_ras = NULL, Shape = NULL, Buffer = 0.5, ID = "ID", Dir = getwd(), Keep_Temporary = FALSE, Source = "USGS"){ - stop("Function currently deprecated as KrigR undergoes major re-development. Please use the stable release to gain access to this functionality.") + stop("This function is deprecated and has been superceeded by the CovariateSetup() function which offers improvements over this implementation. To get started using the new function, you may want to consult the relevant documentation by calling ?CovariateSetup(). Please note that this deprecated function will be removed alltogether when KrigR reaches version 1.0.0.") ### PREPARATION ----- Extent <- extent(Train_ras) # extract extent for later cropping @@ -635,3 +635,329 @@ download_DEM <- function(Train_ras = NULL, return(list(GMTED2010Train_ras, GMTED2010Target_ras)) } + +#' (multi-core) Kriging +#' +#' This function statistically downscales input data using covariate data and the kriging methodology. The function can be run in two ways: +#' \enumerate{ +#' \item \strong{By Itself}: Use the arguments Data, Covariates_coarse, Covariates_fine when you already have raster files for your data which is to be downscaled as well as covariate raster data. +#' \item \strong{From Scratch}: Use the arguments Variable, Type, DataSet, DateStart, DateStop, TResolution, TStep, Extent, Dir, FileName, API_Key, API_User, and arget_res. By doing so, krigR will call the functions download_ERA() and download_DEM() for one coherent kriging workflow. Note that this process does not work when targetting UERRA data. +#' } +#' Use optional arguments such as Dir, FileName, Keep_Temporary, SingularTry, KrigingEquation and Cores for ease of use, substitution of non-GMTED2010 covariates, and parallel processing. +#' +#' @param Data Raster file which is to be downscaled. +#' @param Covariates_coarse Raster file containing covariates at training resolution. +#' @param Covariates_fine Raster file containing covariates at target resolution. +#' @param KrigingEquation Formula or character string specifying which covariates to use and how. Layer names in Covariates_coarse and Covariates_fine need to match Parameters in this formula. Needs to start with "X ~ ". X can read anything you like. +#' @param Dir Optional. Directory specifying where to place final kriged product. Default is current working directory. +#' @param FileName Optional. A file name for the netcdf produced. Default is a combination parameters in the function call. +#' @param Keep_Temporary Logical, whether to delete individual kriging products of layers in Data after processing. Default is TRUE. +#' @param Cores Numeric. How many cores to use. If you want output to your console during the process, use Cores == 1. Parallel processing is carried out when Cores is bigger than 1. Default is detecting all cores of your machine. +#' @param SingularTry Numeric. How often to try kriging of each layer of the input. This usually gets around issues of singular covariance matrices in the kriging process, but takes some time. Default is 10 +#' @param Variable Optional, calls download_ERA(). ERA5(Land)-contained climate variable. +#' @param PrecipFix Optional. Era5(-land) total precipitation is recorded in cumulative steps per hour from the 00:00 time mark per day. Setting PrecipFix to TRUE converts these into records which represent the total precipitation per hour. Monthly records in Era5(-land) express the average daily total precipitation. Setting this argument to TRUE multiplies monthly records by the number of days per the respective month(s) to get to total precipitation records instead of average. Default is FALSE. +#' @param Type Optional. Whether to download reanalysis ('reanalysis') or ensemble ('ensemble_members', 'ensemble_mean', or 'ensemble_spread') data. Passed on to download_ERA. +#' @param DataSet Optional. Which ERA5 data set to download data from. 'era5' or 'era5-land'. Passed on to download_ERA. +#' @param DateStart Optional. Date ('YYYY-MM-DD') at which to start time series of downloaded data. Passed on to download_ERA. +#' @param DateStop Optional. Date ('YYYY-MM-DD') at which to stop time series of downloaded data. Passed on to download_ERA. +#' @param TResolution Optional. Temporal resolution of final product. hour', 'day', 'month'. Passed on to download_ERA. +#' @param TStep Optional. Which time steps (numeric) to consider for temporal resolution. Passed on to download_ERA. +#' @param FUN Optional. A raster calculation argument as passed to `raster::stackApply()`. This controls what kind of data to obtain for temporal aggregates of reanalysis data. Specify 'mean' (default) for mean values, 'min' for minimum values, and 'max' for maximum values, among others. +#' @param Extent Optional, download data according to rectangular bounding box. specify as extent() object or as a raster, a SpatialPolygonsDataFrame object, or a data.frame object. If Extent is a SpatialPolygonsDataFrame, this will be treated as a shapefile and the output will be cropped and masked to this shapefile. If Extent is a data.frame of geo-referenced point records, it needs to contain Lat and Lon columns as well as a non-repeating ID-column. Passed on to download_ERA and download_DEM. +#' @param Buffer Optional. Identifies how big a rectangular buffer to draw around points if Extent is a data frame of points. Buffer is expressed as centessimal degrees. Passed on to download_ERA and download_DEM. +#' @param ID Optional. Identifies which column in Extent to use for creation of individual buffers if Extent is a data.frame. Passed on to download_ERA and download_DEM. +#' @param Target_res Optional. The target resolution for the kriging step (i.e. which resolution to downscale to). An object as specified/produced by raster::res(). Passed on to download_DEM. +#' @param Source Optional, character. Whether to attempt download from the official USGS data viewer (Source = "USGS") or a static copy of the data set on a private drive (Source = "Drive"). Default is "USGS". Use this if the USGS viewer is unavailable. Passed on to download_DEM. +#' @param API_Key Optional. ECMWF cds API key. Passed on to download_ERA. +#' @param API_User Optional. ECMWF cds user number. Passed on to download_ERA. +#' @param nmax Optional. Controls local kriging. Number of nearest observations to be used kriging of each observation. Default is to use all available (Inf). You can specify as a number (numeric). +#' @param TryDown Optional, numeric. How often to attempt the download of each individual file (if querying data download) that the function queries from the server. This is to circumvent having to restart the entire function when encountering connectivity issues. +#' @param verbose Optional, logical. Whether to report progress of data download (if queried) in the console or not. +#' @param TimeOut Numeric. The timeout for each download in seconds. Default 36000 seconds (10 hours). +#' @param SingularDL Logical. Whether to force download of data in one call to CDS or automatically break download requests into individual monthly downloads. Default is FALSE. +#' @return A list object containing the downscaled data as well as the standard error for downscaling as well as the call to the krigR function, and two NETCDF (.nc) file in the specified directory which are the two data contents of the aforementioned list. A temporary directory is populated with individual NETCDF (.nc) files throughout the runtime of krigR which is deleted upon completion if Keep_Temporary = TRUE and all layers in the Data raster object were kriged successfully. +#' @examples +#' \dontrun{ +#' ## THREE-STEP PROCESS (By Itself) +#' # Downloading ERA5-Land air temperature reanalysis data in 12-hour intervals for 02/01/1995 - 04/01/1995 (DD/MM/YYYY). API User and Key in this example are non-functional. Substitute with your user number and key to run this example. +#' Extent <- extent(c(11.8, 15.1, 50.1, 51.7)) # roughly the extent of Saxony +#' API_User <- "..." +#' API_Key <- "..." +#' State_Raw <- download_ERA( +#' Variable = "2m_temperature", +#' DataSet = "era5-land", +#' DateStart = "1995-01-02", +#' DateStop = "1995-01-04", +#' TResolution = "hour", +#' TStep = 12, +#' Extent = Extent, +#' API_User = API_User, +#' API_Key = API_Key +#' ) +#' State_Raw # a raster brick with 6 layers at resolution of ~0.1° +#' # Downloading GMTED2010-data at resolution and extent obtained by a call to download_ERA and a target resolution of .02. +#' Covs_ls <- download_DEM( +#' Train_ras = State_Raw, +#' Target_res = .02, +#' Keep_Temporary = TRUE +#' ) +#' Covs_ls # a list with two elements: (1) GMTED 2010 data at training resolution, and (2) GMTED 2010 data aggregated as close as possible to a resolution of 0.02 +#' # Kriging the data sets prepared with the previous functions. +#' State_Krig <- krigR( +#' Data = State_Raw, # data we want to krig as a raster object +#' Covariates_coarse = Covs_ls[[1]], # training covariate as a raster object +#' Covariates_fine = Covs_ls[[2]], # target covariate as a raster object +#' ) +#' +#' ## PIPELINE (From Scratch) +#' #' # Downloading ERA5-Land air temperature reanalysis data in 12-hour intervals for 02/01/1995 - 04/01/1995 (DD/MM/YYYY), downloading and preparing GMTED 2010 covariate data, and kriging. API User and Key in this example are non-functional. Substitute with your user number and key to run this example. This example produces the same output as the example above. +#' Extent <- extent(c(11.8, 15.1, 50.1, 51.7)) # roughly the extent of Saxony +#' API_User <- "..." +#' API_Key <- "..." +#' Pipe_Krig <- krigR( +#' Variable = "2m_temperature", +#' Type = "reanalysis", +#' DataSet = "era5-land", +#' DateStart = "1995-01-02", +#' DateStop = "1995-01-04", +#' TResolution = "hour", # +#' TStep = 12, +#' Extent = Extent, +#' API_User = API_User, +#' API_Key = API_Key, +#' Target_res = .02, +#' ) +#' } +#' +#' @export +krigR <- function(Data = NULL, Covariates_coarse = NULL, Covariates_fine = NULL, KrigingEquation = "ERA ~ DEM", Cores = detectCores(), Dir = getwd(), FileName, Keep_Temporary = TRUE, SingularTry = 10, Variable, PrecipFix = FALSE, Type = "reanalysis", DataSet = "era5-land", DateStart, DateStop, TResolution = "month", TStep = 1, FUN = "mean", Extent, Buffer = 0.5, ID = "ID", API_Key, API_User, Target_res, Source = "USGS", nmax = Inf, TryDown = 10, verbose = TRUE, TimeOut = 36000, SingularDL = FALSE, ...) { + stop("This function is deprecated and has been superceeded by the Kriging() function which offers improvements over this implementation. To get started using the new function, you may want to consult the relevant documentation by calling ?Kriging(). Please note that this deprecated function will be removed alltogether when KrigR reaches version 1.0.0.") + + ## CALL LIST (for storing how the function as called in the output) ---- + if (is.null(Data)) { + Data_Retrieval <- list( + Variable = Variable, + Type = Type, + PrecipFix = PrecipFix, + DataSet = DataSet, + DateStart = DateStart, + DateStop = DateStop, + TResolution = TResolution, + TStep = TStep, + Extent = Extent + ) + } else { + Data_Retrieval <- "None needed. Data was not queried via krigR function, but supplied by user." + } + ## CLIMATE DATA (call to download_ERA function if no Data set is specified) ---- + if (is.null(Data)) { # data check: if no data has been specified + Data <- download_ERA(Variable = Variable, PrecipFix = PrecipFix, Type = Type, DataSet = DataSet, DateStart = DateStart, DateStop = DateStop, TResolution = TResolution, TStep = TStep, FUN = FUN, Extent = Extent, API_User = API_User, API_Key = API_Key, Dir = Dir, TryDown = TryDown, verbose = verbose, ID = ID, Cores = Cores, TimeOut = TimeOut, SingularDL = SingularDL) + } # end of data check + + ## COVARIATE DATA (call to download_DEM function when no covariates are specified) ---- + if (is.null(Covariates_coarse) & is.null(Covariates_fine)) { # covariate check: if no covariates have been specified + if (class(Extent) == "SpatialPolygonsDataFrame" | class(Extent) == "data.frame") { # Extent check: if Extent input is a shapefile + Shape <- Extent # save shapefile for use as Shape in masking covariate data + } else { # if Extent is not a shape, then extent specification is already baked into Data + Shape <- NULL # set Shape to NULL so it is ignored in download_DEM function when masking is applied + } # end of Extent check + Covs_ls <- download_DEM(Train_ras = Data, Target_res = Target_res, Shape = Shape, Buffer = Buffer, ID = ID, Keep_Temporary = Keep_Temporary, Dir = Dir, Source = Source) + Covariates_coarse <- Covs_ls[[1]] # extract coarse covariates from download_DEM output + Covariates_fine <- Covs_ls[[2]] # extract fine covariates from download_DEM output + } # end of covariate check + + ## KRIGING FORMULA (assure that KrigingEquation is a formula object) ---- + KrigingEquation <- as.formula(KrigingEquation) + + ## CALL LIST (for storing how the function as called in the output) ---- + Call_ls <- list( + Data = SummarizeRaster(Data), + Covariates_coarse = SummarizeRaster(Covariates_coarse), + Covariates_fine = SummarizeRaster(Covariates_fine), + KrigingEquation = KrigingEquation, + Cores = Cores, + FileName = FileName, + Keep_Temporary = Keep_Temporary, + nmax = nmax, + Data_Retrieval = Data_Retrieval + ) + + ## SANITY CHECKS (step into check_Krig function to catch most common error messages) ---- + Check_Product <- check_Krig(Data = Data, CovariatesCoarse = Covariates_coarse, CovariatesFine = Covariates_fine, KrigingEquation = KrigingEquation) + KrigingEquation <- Check_Product[[1]] # extract KrigingEquation (this may have changed in check_Krig) + DataSkips <- Check_Product[[2]] # extract which layers to skip due to missing data (this is unlikely to ever come into action) + Terms <- unique(unlist(strsplit(labels(terms(KrigingEquation)), split = ":"))) # identify which layers of data are needed + + ## DATA REFORMATTING (Kriging requires spatially referenced data frames, reformatting from rasters happens here) --- + Origin <- raster::as.data.frame(Covariates_coarse, xy = TRUE) # extract covariate layers + Origin <- Origin[, c(1:2, which(colnames(Origin) %in% Terms))] # retain only columns containing terms + + Target <- raster::as.data.frame(Covariates_fine, xy = TRUE) # extract covariate layers + Target <- Target[, c(1:2, which(colnames(Target) %in% Terms))] # retain only columns containing terms + Target <- na.omit(Target) + suppressWarnings(gridded(Target) <- ~ x + y) # establish a gridded data product ready for use in kriging + Target@grid@cellsize[1] <- Target@grid@cellsize[2] # ensure that grid cells are square + + ## SET-UP TEMPORARY DIRECTORY (this is where kriged products of each layer will be saved) ---- + Dir.Temp <- file.path(Dir, paste("Kriging", FileName, sep = "_")) + if (!dir.exists(Dir.Temp)) { + dir.create(Dir.Temp) + } + + ## KRIGING SPECIFICATION (this will be parsed and evaluated in parallel and non-parallel evaluations further down) ---- + looptext <- " + OriginK <- cbind(Origin, raster::extract(x = Data[[Iter_Krige]], y = Origin[,1:2], df=TRUE)[, 2]) # combine data of current data layer with training covariate data + OriginK <- na.omit(OriginK) # get rid of NA cells + colnames(OriginK)[length(Terms)+3] <- c(terms(KrigingEquation)[[2]]) # assign column names + suppressWarnings(gridded(OriginK) <- ~x+y) # generate gridded product + OriginK@grid@cellsize[1] <- OriginK@grid@cellsize[2] # ensure that grid cells are square + + Iter_Try = 0 # number of tries set to 0 + kriging_result <- NULL + while(class(kriging_result)[1] != 'autoKrige' & Iter_Try < SingularTry){ # try kriging SingularTry times, this is because of a random process of variogram identification within the automap package that can fail on smaller datasets randomly when it isn't supposed to + try(invisible(capture.output(kriging_result <- autoKrige(formula = KrigingEquation, input_data = OriginK, new_data = Target, nmax = nmax))), silent = TRUE) + Iter_Try <- Iter_Try +1 + } + if(class(kriging_result)[1] != 'autoKrige'){ # give error if kriging fails + message(paste0('Kriging failed for layer ', Iter_Krige, '. Error message produced by autoKrige function: ', geterrmessage())) + } + + ## retransform to raster + try( # try fastest way - this fails with certain edge artefacts in meractor projection and is fixed by using rasterize + Krig_ras <- raster(x = kriging_result$krige_output, layer = 1), # extract raster from kriging product + silent = TRUE + ) + try( + Var_ras <- raster(x = kriging_result$krige_output, layer = 3), # extract raster from kriging product + silent = TRUE + ) + if(!exists('Krig_ras') & !exists('Var_ras')){ + Krig_ras <- rasterize(x = kriging_result$krige_output, y = Covariates_fine[[1]])[[2]] # extract raster from kriging product + Var_ras <- rasterize(x = kriging_result$krige_output, y = Covariates_fine)[[4]] # extract raster from kriging product + } + crs(Krig_ras) <- crs(Data) # setting the crs according to the data + crs(Var_ras) <- crs(Data) # setting the crs according to the data + + if(Cores == 1){ + Ras_Krig[[Iter_Krige]] <- Krig_ras + Ras_Var[[Iter_Krige]] <- Var_ras + } # stack kriged raster into raster list if non-parallel computing + + terra::writeCDF(x = as(brick(Krig_ras), 'SpatRaster'), filename = file.path(Dir.Temp, paste0(str_pad(Iter_Krige,4,'left','0'), '_data.nc')), overwrite = TRUE) + # writeRaster(x = Krig_ras, filename = file.path(Dir.Temp, paste0(str_pad(Iter_Krige,4,'left','0'), '_data.nc')), overwrite = TRUE, format='CDF') # save kriged raster to temporary directory + terra::writeCDF(x = as(brick(Var_ras), 'SpatRaster'), filename = file.path(Dir.Temp, paste0(str_pad(Iter_Krige,4,'left','0'), '_SE.nc')), overwrite = TRUE) + # writeRaster(x = Var_ras, filename = file.path(Dir.Temp, paste0(str_pad(Iter_Krige,4,'left','0'), '_SE.nc')), overwrite = TRUE, format='CDF') # save kriged raster to temporary directory + + if(Cores == 1){ # core check: if processing non-parallel + if(Count_Krige == 1){ # count check: if this was the first actual computation + T_End <- Sys.time() # record time at which kriging was done for current layer + Duration <- as.numeric(T_End)-as.numeric(T_Begin) # calculate how long it took to krig on layer + message(paste('Kriging of remaining ', nlayers(Data)-Iter_Krige, ' data layers should finish around: ', as.POSIXlt(T_Begin + Duration*nlayers(Data), tz = Sys.timezone(location=TRUE)), sep='')) # console output with estimate of when the kriging should be done + ProgBar <- txtProgressBar(min = 0, max = nlayers(Data), style = 3) # create progress bar when non-parallel processing + Count_Krige <- Count_Krige + 1 # raise count by one so the stimator isn't called again + } # end of count check + setTxtProgressBar(ProgBar, Iter_Krige) # update progress bar with number of current layer + } # end of core check + " + + ## KRIGING PREPARATION (establishing objects which the kriging refers to) ---- + Ras_Krig <- as.list(rep(NA, nlayers(Data))) # establish an empty list which will be filled with kriged layers + Ras_Var <- as.list(rep(NA, nlayers(Data))) # establish an empty list which will be filled with kriged layers + + if (verbose) { + message("Commencing Kriging") + } + ## DATA SKIPS (if certain layers in the data are empty and need to be skipped, this is handled here) --- + if (!is.null(DataSkips)) { # Skip check: if layers need to be skipped + for (Iter_Skip in DataSkips) { # Skip loop: loop over all layers that need to be skipped + Ras_Krig[[Iter_Skip]] <- Data[[Iter_Skip]] # add raw data (which should be empty) to list + terra::writeCDF(x = as(brick(Ras_Krig[[Iter_Skip]]), "SpatRaster"), filename = file.path(Dir.Temp, str_pad(Iter_Skip, 4, "left", "0")), overwrite = TRUE) + # writeRaster(x = Ras_Krig[[Iter_Skip]], filename = file.path(Dir.Temp, str_pad(Iter_Skip,4,'left','0')), overwrite = TRUE, format = 'CDF') # save raw layer to temporary directory, needed for loading back in when parallel processing + } # end of Skip loop + Layers_vec <- 1:nlayers(Data) # identify vector of all layers in data + Compute_Layers <- Layers_vec[which(!Layers_vec %in% DataSkips)] # identify which layers can actually be computed on + } else { # if we don't need to skip any layers + Compute_Layers <- 1:nlayers(Data) # set computing layers to all layers in data + } # end of Skip check + + + ## ACTUAL KRIGING (carry out kriging according to user specifications either in parallel or on a single core) ---- + if (Cores > 1) { # Cores check: if parallel processing has been specified + ### PARALLEL KRIGING --- + ForeachObjects <- c("Dir.Temp", "Cores", "Data", "KrigingEquation", "Origin", "Target", "Covariates_coarse", "Covariates_fine", "Terms", "SingularTry", "nmax") # objects which are needed for each kriging run and are thus handed to each cluster unit + pb <- txtProgressBar(max = length(Compute_Layers), style = 3) + progress <- function(n) { + setTxtProgressBar(pb, n) + } + opts <- list(progress = progress) + cl <- makeCluster(Cores) # Assuming Cores node cluster + registerDoSNOW(cl) # registering cores + foreach( + Iter_Krige = Compute_Layers, # kriging loop over all layers in Data, with condition (%:% when(...)) to only run if current layer is not present in Dir.Temp yet + .packages = c("raster", "stringr", "automap", "ncdf4", "rgdal", "terra"), # import packages necessary to each itteration + .export = ForeachObjects, + .options.snow = opts + ) %:% when(!paste0(str_pad(Iter_Krige, 4, "left", "0"), "_data.nc") %in% list.files(Dir.Temp)) %dopar% { # parallel kriging loop + Ras_Krig <- eval(parse(text = looptext)) # evaluate the kriging specification per cluster unit per layer + } # end of parallel kriging loop + close(pb) + stopCluster(cl) # close down cluster + Files_krig <- list.files(Dir.Temp)[grep(pattern = "_data.nc", x = list.files(Dir.Temp))] + Files_var <- list.files(Dir.Temp)[grep(pattern = "_SE.nc", x = list.files(Dir.Temp))] + for (Iter_Load in 1:length(Files_krig)) { # load loop: load data from temporary files in Dir.Temp + Ras_Krig[[Iter_Load]] <- raster(file.path(Dir.Temp, Files_krig[Iter_Load])) # load current temporary file and write contents to list of rasters + Ras_Var[[Iter_Load]] <- raster(file.path(Dir.Temp, Files_var[Iter_Load])) # load current temporary file and write contents to list of rasters + } # end of load loop + } else { # if non-parallel processing has been specified + ### NON-PARALLEL KRIGING --- + Count_Krige <- 1 # Establish count variable which is targeted in kriging specification text for producing an estimator + for (Iter_Krige in Compute_Layers) { # non-parallel kriging loop over all layers in Data + if (paste0(str_pad(Iter_Krige, 4, "left", "0"), "_data.nc") %in% list.files(Dir.Temp)) { # file check: if this file has already been produced + Ras_Krig[[Iter_Krige]] <- raster(file.path(Dir.Temp, paste0(str_pad(Iter_Krige, 4, "left", "0"), "_data.nc"))) # load already produced kriged file and save it to list of rasters + Ras_Var[[Iter_Krige]] <- raster(file.path(Dir.Temp, paste0(str_pad(Iter_Krige, 4, "left", "0"), "_SE.nc"))) + if (!exists("ProgBar")) { + ProgBar <- txtProgressBar(min = 0, max = nlayers(Data), style = 3) + } # create progress bar when non-parallel processing} + setTxtProgressBar(ProgBar, Iter_Krige) # update progress bar + next() # jump to next layer + } # end of file check + T_Begin <- Sys.time() # record system time when layer kriging starts + eval(parse(text = looptext)) # evaluate the kriging specification per layer + } # end of non-parallel kriging loop + } # end of Cores check + + ## SAVING FINAL PRODUCT ---- + if (is.null(DataSkips)) { # Skip check: if no layers needed to be skipped + # convert list of kriged layers in actual rasterbrick of kriged layers + names(Ras_Krig) <- names(Data) + if (class(Ras_Krig) != "RasterBrick") { + Ras_Krig <- brick(Ras_Krig) + } + Krig_terra <- as(Ras_Krig, "SpatRaster") + names(Krig_terra) <- names(Data) + terra::writeCDF(x = Krig_terra, filename = file.path(Dir, paste0(FileName, ".nc")), overwrite = TRUE) + # writeRaster(x = Ras_Krig, filename = file.path(Dir, FileName), overwrite = TRUE, format="CDF") # save final product as raster + # convert list of kriged layers in actual rasterbrick of kriged layers + names(Ras_Var) <- names(Data) + if (class(Ras_Var) != "RasterBrick") { + Ras_Var <- brick(Ras_Var) + } + Var_terra <- as(Ras_Var, "SpatRaster") + names(Var_terra) <- names(Data) + + terra::writeCDF(x = Var_terra, filename = file.path(Dir, paste0("SE_", paste0(FileName, ".nc"))), overwrite = TRUE) + # writeRaster(x = Ras_Var, filename = file.path(Dir, paste0("SE_",FileName)), overwrite = TRUE, format="CDF") # save final product as raster + } else { # if some layers needed to be skipped + warning(paste0("Some of the layers in your raster could not be kriged. You will find all the individual layers (kriged and not kriged) in ", Dir, ".")) + Keep_Temporary <- TRUE # keep temporary files so kriged products are not deleted + } # end of Skip check + + ### REMOVE FILES FROM HARD DRIVE --- + if (Keep_Temporary == FALSE) { # cleanup check + unlink(Dir.Temp, recursive = TRUE) + } # end of cleanup check + + Krig_ls <- list(Ras_Krig, Ras_Var, Call_ls) + names(Krig_ls) <- c("Kriging_Output", "Kriging_SE", "Call") + return(Krig_ls) # return raster or list of layers +} diff --git a/R/Kriging.R b/R/Kriging.R index fd5ad48..2cb9bfc 100644 --- a/R/Kriging.R +++ b/R/Kriging.R @@ -12,6 +12,7 @@ #' @param Dir Character/Directory Pointer. Directory specifying where to place final kriged product. Default is current working directory. #' @param FileName Character. A file name for the produced files. #' @param FileExtension Character. A file extension for the produced file. Supported values are ".nc" (default) and ".tif" (better support for metadata). +#' @param Compression Integer between 1 to 9. Applied to final .nc file that the function writes to hard drive. Same as compression argument in terra::writeCDF(). Ignored if FileExtension = ".tif". #' @param Keep_Temporary Logical, whether to delete individual kriging products of layers in Data after processing. Default is TRUE. These temporary files are stored in a newly created directory in Dir which is pre-pended with "TEMP-" and is deleted if Keep_Temporary = FALSE upon completion. #' @param verbose Optional, logical. Whether to report progress of data download (if queried) in the console or not. #' @@ -33,6 +34,7 @@ #' @importFrom automap autoKrige #' @importFrom stringr str_pad #' @importFrom doSNOW registerDoSNOW +#' @importFrom parallel detectCores #' @importFrom parallel makeCluster #' @importFrom parallel stopCluster #' @importFrom foreach %dopar% @@ -48,16 +50,17 @@ #' #' \strong{ATTENTION:} If data is loaded again from disk at a later point with a different function, take note that citation and KrigR-call metadata will not be loaded properly from a .nc when loading data through a different function. Kriging() handles these .nc specific issues when loading .nc files created previously with Kriging() from disk. #' -#' @seealso \code{\link{CovariateSetup}}. +#' @seealso \code{\link{CovariateSetup}}, \code{\link{Plot.Kriged}}. #' #' @examples #' \dontrun{ #' ## Kriging using pre-fab data with a rectangular extent and a fives layers of data with parallel processing #' ### Loading data -#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package="KrigR")) -#' Cov_train <- terra::rast(system.file("extdata", "Covariates_Train.nc", package="KrigR")) -#' Cov_target <- terra::rast(system.file("extdata", "Covariates_Target.nc", package="KrigR")) -#' terra::varnames(Cov_train) <- terra::varnames(Cov_target) <- "GMTED2010" # we must ensure that the varnames in the Covariate file match +#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +#' Cov_train <- terra::rast(system.file("extdata", "Covariates_Train.nc", package = "KrigR")) +#' Cov_target <- terra::rast(system.file("extdata", "Covariates_Target.nc", package = "KrigR")) +#' names(Cov_train) <- names(Cov_target) <- "GMTED2010" +#' #' ### kriging itself #' ExtentKrig <- Kriging( #' Data = CDS_rast, @@ -71,7 +74,7 @@ #' nmax = 40, #' verbose = TRUE #' ) -#' +#' Plot.Kriged(Krigs = ExtentKrig) #' ## Kriging using full KrigR pipeline with shapefile data #' ### Shapefile loading #' data("Jotunheimen_poly") @@ -94,19 +97,21 @@ #' ) #' #' ### Covariate preparations -#' Covariates_ls <- CovariateSetup(Training = Qsoil_rast, -#' #' Target = 0.03, -#' Covariates = "GMTED2010", # this shiuld really be HWSD -#' Extent = Jotunheimen_poly, -#' Keep_Global = TRUE) -#' terra::varnames(Covariates_ls[[1]]) <- terra::varnames(Covariates_ls[[2]]) <- "GMTED2010" # we must ensure that the varnames in the Covariate file match +#' Covariates_ls <- CovariateSetup( +#' Training = Qsoil_rast, +#' Target = 0.03, +#' Covariates = c("tksat", "tkdry", "csol", "k_s", "lambda", "psi", "theta_s"), +#' Source = "Drive", +#' Extent = Jotunheimen_poly, +#' Keep_Global = TRUE +#' ) #' #' ### kriging itself #' ShapeKrig <- Kriging( #' Data = Qsoil_rast, #' Covariates_training = Covariates_ls[[1]], #' Covariates_target = Covariates_ls[[2]], -#' Equation = "GMTED2010", +#' Equation = "tksat + tkdry + csol + k_s + lambda + psi + theta_s", #' Cores = 1, #' FileName = "KrigTest2", #' FileExtension = ".nc", @@ -114,33 +119,40 @@ #' nmax = 40, #' verbose = TRUE #' ) +#' Plot.Kriged(Krigs = ShapeKrig, SF = Jotunheimen_poly) #' } #' @export -Kriging <- function( +Kriging <- function( # nolint: cyclocomp_linter. Data, Covariates_training, Covariates_target, Equation = NULL, - Cores = detectCores(), + Cores = parallel::detectCores(), nmax = Inf, Dir = getwd(), FileName, - FileExtension, + FileExtension = ".nc", + Compression = 9, Keep_Temporary = FALSE, - verbose = TRUE - ){ + verbose = TRUE) { ## Run Preparations =============== + if (verbose) { + message("###### Checking your Kriging Specification") + } ### Number of layers in data & Progress bar KrigIterations <- nlyr(Data) # used for krig looping pb <- progress_bar$new( format = "Kriging (:current/:total) | [:bar] Elapsed: :elapsed | Remaining: :eta", - total = KrigIterations, # 100 - width = getOption("width")) - progress_layer <- 1:KrigIterations # token reported in progress bar + total = KrigIterations, # 100 + width = getOption("width") + ) + progress_layer <- 1:KrigIterations # token reported in progress bar ### CRS for assignment in loop CRS_dat <- crs(Data) ### if no equation is specified, assign additive combination of variables in training covariates - if(is.null(Equation)){Equation <- paste(terra::varnames(Covariates_training), collapse = " + ")} + if (is.null(Equation)) { + Equation <- paste(terra::names(Covariates_training), collapse = " + ") + } ### assure that KrigingEquation is a formula object KrigingEquation <- as.formula(paste("Data ~", Equation)) ### Metadata @@ -152,16 +164,18 @@ Kriging <- function( "KrigRCall" = Meta_vec ) ### Temporary Directory - Dir.Temp <- file.path(Dir, paste("TEMP-Kriging", FileName, sep="_")) - if(!dir.exists(Dir.Temp)){dir.create(Dir.Temp)} + Dir.Temp <- file.path(Dir, paste("TEMP-Kriging", FileName, sep = "_")) + if (!dir.exists(Dir.Temp)) { + dir.create(Dir.Temp) + } ### Remove extension from file name FileName <- file_path_sans_ext(FileName) ## Check if already executed once =============== FCheck1 <- Check.File(FName = paste0(FileName, "_Kriged", FileExtension), Dir = Dir, loadFun = terra::rast, load = TRUE, verbose = TRUE) FCheck2 <- Check.File(FName = paste0(FileName, "_StDev", FileExtension), Dir = Dir, loadFun = terra::rast, load = TRUE, verbose = TRUE) - if(!is.null(FCheck1)){ - if(FileExtension == ".nc"){ + if (!is.null(FCheck1)) { + if (FileExtension == ".nc") { FCheck1 <- Meta.NC(NC = FCheck1, FName = file.path(Dir, paste0(FileName, "_Kriged.nc")), Attrs = Meta_vec, Read = TRUE) FCheck2 <- Meta.NC(NC = FCheck2, FName = file.path(Dir, paste0(FileName, "_StDev.nc")), Attrs = Meta_vec, Read = TRUE) } @@ -171,25 +185,28 @@ Kriging <- function( terra::metags(FCheck1) <- terra::metags(FCheck2) <- Meta_vec Krig_ls <- list(FCheck1, FCheck2) names(Krig_ls) <- c("Prediction", "StDev") + unlink(Dir.Temp, recursive = TRUE) return(Krig_ls) } ## Catching Most Frequent Issues =============== Check_Product <- Check.Krig(Data = Data, CovariatesCoarse = Covariates_training, CovariatesFine = Covariates_target, KrigingEquation = KrigingEquation) KrigingEquation <- Check_Product[[1]] # extract KrigingEquation (this may have changed in check_Krig) - # DataSkips <- Check_Product[[2]] # extract which layers to skip due to missing data (this is unlikely to ever come into action) Terms <- unique(unlist(strsplit(labels(terms(KrigingEquation)), split = ":"))) # identify which layers of data are needed ## Data Reformatting =============== + if (verbose) { + message("###### Preparing And Reformatting your Data") + } # (Kriging requires spatially referenced data frames, reformatting from rasters happens here) ### Make Training sf object Origin <- as.data.frame(Covariates_training, xy = TRUE, na.rm = FALSE) - colnames(Origin)[-1:-2] <- terra::varnames(Covariates_training) + colnames(Origin)[-1:-2] <- names(Covariates_training) Origin <- Origin[, c(1:2, which(colnames(Origin) %in% Terms))] # retain only columns containing terms Origin <- st_as_sf(Origin, coords = c("x", "y")) ### Make Target sf object Target <- as.data.frame(Covariates_target, xy = TRUE) - colnames(Target)[-1:-2] <- terra::varnames(Covariates_target) + colnames(Target)[-1:-2] <- names(Covariates_target) Target <- Target[, c(1:2, which(colnames(Target) %in% Terms))] # retain only columns containing terms Target <- st_as_sf(Target, coords = c("x", "y")) ### Make data into data frame for handling in parallel (SpatRasters cannot be used in foreach) @@ -249,41 +266,45 @@ Kriging <- function( ## Kriging Execution =============== ## carry out kriging according to user specifications either in parallel or on a single core - if(verbose){message("Commencing Kriging")} + if (verbose) { + message("###### Commencing Kriging") + } ### multi-core kriging ---- - if(Cores > 1){ + if (Cores > 1) { ## registering cluster and progress bar for foreach cl <- makeCluster(Cores) registerDoSNOW(cl) - progress <- function(n){ + progress <- function(n) { pb$tick(tokens = list(layer = progress_layer[n])) } on.exit(stopCluster(cl)) ## executing foreach kriging ForeachObjects <- c("Dir.Temp", "Cores", "CRS_dat", "Data_df", "KrigingEquation", "Origin", "Target", "nmax", "FileExtension") # objects which are needed in Kriging - foreach(Iter_Krige = 1:KrigIterations, # kriging loop over all layers in Data, with condition (%:% when(...)) to only run if current layer is not present in Dir.Temp yet - .packages = c("terra", "sf", "stringr", "automap", "ncdf4"), # import packages necessary to each iteration - .export = ForeachObjects, - .options.snow = list(progress = progress)) %dopar% { # parallel kriging loop # %:% when(!paste0(str_pad(Iter_Krige,7,"left","0"), '_data.nc') %in% list.files(Dir.Temp)) - eval(parse(text=looptext)) # evaluate the kriging specification per cluster unit per layer - Sys.sleep(0.5) - NULL - } # end of parallel kriging loop + foreach( + Iter_Krige = 1:KrigIterations, # kriging loop over all layers in Data, with condition (%:% when(...)) to only run if current layer is not present in Dir.Temp yet + .packages = c("terra", "sf", "stringr", "automap", "ncdf4"), # import packages necessary to each iteration + .export = ForeachObjects, + .options.snow = list(progress = progress) + ) %dopar% { # parallel kriging loop # %:% when(!paste0(str_pad(Iter_Krige,7,"left","0"), '_data.nc') %in% list.files(Dir.Temp)) + eval(parse(text = looptext)) # evaluate the kriging specification per cluster unit per layer + Sys.sleep(0.5) + NULL + } # end of parallel kriging loop } ### single-core kriging ---- - if(Cores == 1){ - for(Iter_Krige in 1:KrigIterations){ - # FileExis <- paste0(str_pad(Iter_Krige,7,'left','0'), '_data', FileExtension) %in% list.files(Dir.Temp) - # if(!FileExis){ - eval(parse(text=looptext)) # evaluate the kriging specification per layer - # } + if (Cores == 1) { + for (Iter_Krige in 1:KrigIterations) { + eval(parse(text = looptext)) # evaluate the kriging specification per layer Sys.sleep(0.5) pb$tick(tokens = list(layer = progress_layer[Iter_Krige])) } } ## Data Loading and Saving =============== + if (verbose) { + message("###### Loading Temporary Files") + } ### loading kriged data back in Krig_rast <- rast(list.files(Dir.Temp, full.names = TRUE, pattern = "_data")) SE_rast <- rast(list.files(Dir.Temp, full.names = TRUE, pattern = "_StDev")) @@ -293,21 +314,37 @@ Kriging <- function( terra::units(Krig_rast) <- terra::units(SE_rast) <- terra::units(Data) terra::metags(Krig_rast) <- terra::metags(SE_rast) <- Meta_vec ### Data Saving - if(FileExtension == ".tif"){ + if (verbose) { + message("###### Saving Kriged Data") + } + if (FileExtension == ".tif") { writeRaster(Krig_rast, filename = file.path(Dir, paste0(FileName, "_Kriged", FileExtension))) + Krig_rast <- rast(filename = file.path(Dir, paste0(FileName, "_Kriged", FileExtension))) writeRaster(SE_rast, filename = file.path(Dir, paste0(FileName, "_StDev", FileExtension))) + SE_rast <- rast(filename = file.path(Dir, paste0(FileName, "_StDev", FileExtension))) } - if(FileExtension == ".nc"){ - Krig_rast <- Meta.NC(NC = Krig_rast, FName = file.path(Dir, paste0(FileName, "_Kriged", FileExtension)), - Attrs = terra::metags(Krig_rast), Write = TRUE) - SE_rast <- Meta.NC(NC = SE_rast, FName = file.path(Dir, paste0(FileName, "_STDev", FileExtension)), - Attrs = terra::metags(SE_rast), Write = TRUE) + if (FileExtension == ".nc") { + Krig_rast <- Meta.NC( + NC = Krig_rast, FName = file.path(Dir, paste0(FileName, "_Kriged", FileExtension)), + Attrs = terra::metags(Krig_rast), Write = TRUE, Compression = Compression + ) + SE_rast <- Meta.NC( + NC = SE_rast, FName = file.path(Dir, paste0(FileName, "_STDev", FileExtension)), + Attrs = terra::metags(SE_rast), Write = TRUE, Compression = Compression + ) } + Krig_rast <- rast(file.path(Dir, paste0(FileName, "_Kriged", FileExtension))) + SE_rast <- rast(file.path(Dir, paste0(FileName, "_STDev", FileExtension))) + terra::time(Krig_rast) <- terra::time(SE_rast) <- terra::time(Data) + terra::varnames(Krig_rast) <- terra::varnames(SE_rast) <- terra::varnames(Data) + terra::units(Krig_rast) <- terra::units(SE_rast) <- terra::units(Data) + terra::metags(Krig_rast) <- terra::metags(SE_rast) <- Meta_vec + ## Removing Temporary Files =============== - if(Keep_Temporary == FALSE){ # cleanup check + if (Keep_Temporary == FALSE) { # cleanup check unlink(Dir.Temp, recursive = TRUE) - } # end of cleanup check + } # end of cleanup check ## Data Return =============== Krig_ls <- list(Krig_rast, SE_rast) diff --git a/R/Meta.R b/R/Meta.R index 0857499..15dbd5f 100644 --- a/R/Meta.R +++ b/R/Meta.R @@ -9,7 +9,7 @@ #' #' @seealso \code{\link{Meta.List}}, \code{\link{Meta.Read}}, \code{\link{Meta.Variables}}, \code{\link{Meta.DOI}}, \code{\link{Meta.QuickFacts}}. #' -Meta.Register <- function(Dir = file.path(getwd(), "metadata")){ +Meta.Register <- function(Dir = file.path(getwd(), "metadata")) { sink(file = file.path(Dir, "metadata.txt")) cat(list.files(Dir, ".rds"), sep = "\n") sink() @@ -32,9 +32,8 @@ Meta.Register <- function(Dir = file.path(getwd(), "metadata")){ #' Meta.List() #' #' @export -Meta.List <- function(URL = "https://raw.githubusercontent.com/ErikKusch/KrigR/Development/metadata" - ){ - file_path_sans_ext(read.table(file.path(URL, "metadata.txt"))[,1]) +Meta.List <- function(URL = "https://raw.githubusercontent.com/ErikKusch/KrigR/Development/metadata") { + file_path_sans_ext(read.table(file.path(URL, "metadata.txt"))[, 1]) } ### READ METADATA FACTS ======================================================== @@ -54,12 +53,12 @@ Meta.List <- function(URL = "https://raw.githubusercontent.com/ErikKusch/KrigR/D #' #' @export Meta.Read <- function(URL = "https://raw.githubusercontent.com/ErikKusch/KrigR/Development/metadata", ## change this to github repo for these data once ready - dataset = "reanalysis-era5-land"){ + dataset = "reanalysis-era5-land") { load(url( paste0( - "https://github.com/ErikKusch/KrigR/blob/Development/metadata/", - dataset, - ".RData?raw=true" + "https://github.com/ErikKusch/KrigR/blob/Development/metadata/", + dataset, + ".RData?raw=true" ) )) get(ls()[ls() == gsub(dataset, pattern = "-", replacement = "_")]) @@ -87,7 +86,7 @@ Meta.Read <- function(URL = "https://raw.githubusercontent.com/ErikKusch/KrigR/D #' Meta.Variables() #' #' @export -Meta.Variables <- function(dataset = "reanalysis-era5-land"){ +Meta.Variables <- function(dataset = "reanalysis-era5-land") { Meta.Read(dataset = dataset)$Variables } @@ -106,7 +105,7 @@ Meta.Variables <- function(dataset = "reanalysis-era5-land"){ #' Meta.DOI() #' #' @export -Meta.DOI <- function(dataset = "reanalysis-era5-land"){ +Meta.DOI <- function(dataset = "reanalysis-era5-land") { Meta.Read(dataset = dataset)$Citation } @@ -138,11 +137,13 @@ Meta.DOI <- function(dataset = "reanalysis-era5-land"){ #' Meta.QuickFacts() #' #' @export -Meta.QuickFacts <- function(dataset = "reanalysis-era5-land"){ - Meta.Read(dataset = dataset)[c("DataSet", "Type", "URL", "Description", - "TResolution", "TStep", "TStart", "TEnd", - "Projection", "SpatialResolution", - "CDSArguments")] +Meta.QuickFacts <- function(dataset = "reanalysis-era5-land") { + Meta.Read(dataset = dataset)[c( + "DataSet", "Type", "URL", "Description", + "TResolution", "TStep", "TStart", "TEnd", + "Projection", "SpatialResolution", + "CDSArguments" + )] } ### CDS QUERY VALIDATION AGAINST DATA SET METADATA ============================= @@ -174,102 +175,111 @@ Meta.QuickFacts <- function(dataset = "reanalysis-era5-land"){ #' @seealso \code{\link{Meta.List}}, \code{\link{Meta.Read}}, \code{\link{Meta.Variables}}, \code{\link{Meta.DOI}}, \code{\link{Meta.QuickFacts}}. #' #' @examples -#' Meta.Check(DataSet = "reanalysis-era5-land", Type = NA, VariableCheck = "2m_temperature", CumulativeCheck = FALSE, ExtentCheck = c(53.06, 9.87, 49.89, 15.03), DateCheck = data.frame(IN = c(as.POSIXct("1995-01-01 CET"), as.POSIXct("2005-01-01 23:00:00 CET")), UTC = c(as.POSIXct("1994-12-31 23:00:00 UTC"), as.POSIXct("2005-01-01 22:00:00 UTC"))), AggrCheck = list(1, "hour"), QueryTimes = c('00:00', '03:00', '06:00', '09:00', '12:00', '15:00', '18:00', '21:00')) +#' Meta.Check(DataSet = "reanalysis-era5-land", Type = NA, VariableCheck = "2m_temperature", CumulativeCheck = FALSE, ExtentCheck = c(53.06, 9.87, 49.89, 15.03), DateCheck = data.frame(IN = c(as.POSIXct("1995-01-01 CET"), as.POSIXct("2005-01-01 23:00:00 CET")), UTC = c(as.POSIXct("1994-12-31 23:00:00 UTC"), as.POSIXct("2005-01-01 22:00:00 UTC"))), AggrCheck = list(1, "hour"), QueryTimes = c("00:00", "03:00", "06:00", "09:00", "12:00", "15:00", "18:00", "21:00")) #' -Meta.Check <- function(DataSet = "reanalysis-era5-land", Type = NA, VariableCheck, CumulativeCheck, ExtentCheck, DateCheck, AggrCheck, QueryTimes){ +Meta.Check <- function(DataSet = "reanalysis-era5-land", Type = NA, VariableCheck, CumulativeCheck, ExtentCheck, DateCheck, AggrCheck, QueryTimes) { # nolint: cyclocomp_linter. #--- Variable ### if a variable not in the data set has been specified - if(length(VariableCheck) == 0){stop("Please specify a variable provided by the data set. Your can be retrieved with the function call: ", "\n", "Meta.Variables(dataset = '", DataSet, "')")} + if (length(VariableCheck) == 0) { + stop("Please specify a variable provided by the data set. Your can be retrieved with the function call: ", "\n", "Meta.Variables(dataset = '", DataSet, "')") + } #--- Cumulative ### if the cumulative back-calculation is attempting to be applied to a non-cumulative variable CumVar <- Meta.Variables(dataset = DataSet)$Cumulative[which(Meta.Variables(dataset = DataSet)$CDSname == VariableCheck)] - if(CumulativeCheck & !CumVar){ + if (CumulativeCheck && !CumVar) { stop("You have specified to back-calculation of cumulative data for a non-cumulatively recorded variable. This would produce nonsense data. Please specify CumulVar = FALSE instead. For an overview of which variables are recorded cumulatively for the data set you are querying, please consider the function call:", "\n", "Meta.Variables(dataset = '", DataSet, "')") } + if (!CumulativeCheck && CumVar && Meta.QuickFacts(dataset = DataSet)$TResolution == "hour" && attr(DateCheck$IN[1], "tzone") != "UTC") { + warning("You have selected to download data recorded cumulatively in hourly intervals for a different timezone than UTC while not applying back-calculation/disaggregation of these cumulative records. Be aware that the resulting data will likely not be useful as the cumulative window resets at each new day in UTC.") + } #--- Extent ### if an extent outside the data product has been specified - DataExt <- ext(Meta.QuickFacts(dataset = DataSet)$CDSArguments$area)[c(4,1,3,2)] #N,W,S,E - if( + DataExt <- ext(Meta.QuickFacts(dataset = DataSet)$CDSArguments$area)[c(4, 1, 3, 2)] # N,W,S,E + if ( ( # ymax (ExtentCheck[1] > DataExt[1]) + - # xmin - (ExtentCheck[2] < DataExt[2]) + - # ymin - (ExtentCheck[3] < DataExt[3]) + - # ymax - (ExtentCheck[4] > DataExt[4]) + # xmin + (ExtentCheck[2] < DataExt[2]) + + # ymin + (ExtentCheck[3] < DataExt[3]) + + # ymax + (ExtentCheck[4] > DataExt[4]) ) != 0 - ){ - stop("Please specify an area using the Extent argument that is contained within the data set. The data set covers the area defined by the following extent:", - "\n", ext(Meta.QuickFacts(dataset = DataSet)$CDSArguments$area), " in ", Meta.QuickFacts(dataset = DataSet)$Projection) + ) { + stop( + "Please specify an area using the Extent argument that is contained within the data set. The data set covers the area defined by the following extent:", + "\n", ext(Meta.QuickFacts(dataset = DataSet)$CDSArguments$area), " in ", Meta.QuickFacts(dataset = DataSet)$Projection + ) } #--- Time Window ### check if time window is exceeded CheckStart <- DateCheck$UTC[1] < Meta.QuickFacts(dataset = DataSet)$TStart - if(class(Meta.QuickFacts(dataset = DataSet)$TEnd)[1] == "POSIXct"){ + if (class(Meta.QuickFacts(dataset = DataSet)$TEnd)[1] == "POSIXct") { CheckEnd <- (DateCheck$UTC[2] > Meta.QuickFacts(dataset = DataSet)$TEnd) - }else{ + } else { CheckEnd <- FALSE - warning("Cannot validate user-specified end date (DateStop) because specified data set is being updated regularly (", - strsplit(Meta.QuickFacts(dataset = DataSet)$TEnd, split = "; ")[[1]][2], "). User-specification may lead to an error.") + warning( + "Cannot validate user-specified end date (DateStop) because specified data set is being updated regularly (", + strsplit(Meta.QuickFacts(dataset = DataSet)$TEnd, split = "; ")[[1]][2], "). User-specification may lead to an error." + ) } - if(CheckStart + CheckEnd != 0){ - stop("The time window you have specified is not supported by the data set. The data set makes data available from ", - Meta.QuickFacts(dataset = DataSet)$TStart, " until ", Meta.QuickFacts(dataset = DataSet)$TEnd) + if (CheckStart + CheckEnd != 0) { + stop( + "The time window you have specified is not supported by the data set. The data set makes data available from ", + Meta.QuickFacts(dataset = DataSet)$TStart, " until ", Meta.QuickFacts(dataset = DataSet)$TEnd + ) } #--- Aggregation Match ### check if desired aggregation is supported SuppRes <- c("hour", "day", "month", "year") BaseStep <- BaseStep <- Meta.QuickFacts(dataset = DataSet)$TStep[ - na.omit(match(Type, Meta.QuickFacts(dataset = DataSet)$Type))] - if(Meta.QuickFacts(dataset = DataSet)$TResolution != AggrCheck[[2]] | - BaseStep != AggrCheck[[1]]){ # if this is TRUE, we need to check if aggregation works - + na.omit(match(Type, Meta.QuickFacts(dataset = DataSet)$Type)) + ] + if (Meta.QuickFacts(dataset = DataSet)$TResolution != AggrCheck[[2]] || BaseStep != AggrCheck[[1]]) { # if this is TRUE, we need to check if aggregation works ## specification of a temporal resolution finer than the data? - if(which(SuppRes == AggrCheck[[2]]) < which(SuppRes == Meta.QuickFacts(dataset = DataSet)$TResolution)){ - stop("You have specified a temporal aggregation at a scale finer than what the data set reports natively (", - Meta.QuickFacts(dataset = DataSet)$TResolution, "). Please specify the same or a coarser temporal resolution for the TResolution argument. Supported options are '", paste(SuppRes, collapse = "', '"), "'.") + if (which(SuppRes == AggrCheck[[2]]) < which(SuppRes == Meta.QuickFacts(dataset = DataSet)$TResolution)) { + stop( + "You have specified a temporal aggregation at a scale finer than what the data set reports natively (", + Meta.QuickFacts(dataset = DataSet)$TResolution, "). Please specify the same or a coarser temporal resolution for the TResolution argument. Supported options are '", paste(SuppRes, collapse = "', '"), "'." + ) } ## specification of tsteps that cannot be achieved with the data? - if(Meta.QuickFacts(dataset = DataSet)$TResolution == AggrCheck[[2]] & - ((AggrCheck[[1]] /BaseStep) %%1!=0)){ + if (Meta.QuickFacts(dataset = DataSet)$TResolution == AggrCheck[[2]] && ((AggrCheck[[1]] / BaseStep) %% 1 != 0)) { stop("You have specified a temporal aggregation that cannot be achieved with the data. When specifying the same temporal resolution as the data (you have specified TResolution = ", AggrCheck[[2]], "), the TStep must be a multiple of the base temporal resolution of the data (", BaseStep, " for DataSet = ", DataSet, " and Type = ", Type, ").") } ## specification of daily, monthly or annual aggregates but not setting tstart or tend to beginning or end of day/month/year? - if(AggrCheck[[2]] == "day" & - (as.numeric(substr(QueryTimes[1], 0, 2)) != 0 | - as.numeric(substr(QueryTimes[length(QueryTimes)], 0, 2)) != 23)){ + if (AggrCheck[[2]] == "day" && (as.numeric(substr(QueryTimes[1], 0, 2)) != 0 || as.numeric(substr(QueryTimes[length(QueryTimes)], 0, 2)) != 23)) { stop("You have specified (multi-)daily temporal aggregation but are querying a time window which does not start at 00:00 and/or does not terminate at 23:00. Please ensure that you set the argument DateStart and DateStop accordingly.") } ## these may fail when querying monthly raw data MustStartMonth <- as.POSIXct(paste( paste(format(DateCheck$IN[1], "%Y"), format(DateCheck$IN[1], "%m"), "01", sep = "-"), - "00:00:00"), tz = format(DateCheck$IN[2], "%Z")) + "00:00:00" + ), tz = format(DateCheck$IN[2], "%Z")) MustEndMonth <- as.POSIXct(paste( paste(format(DateCheck$IN[2], "%Y"), format(DateCheck$IN[2], "%m"), - days_in_month(DateCheck$IN[2]), sep = "-"), - "24:00:00"), tz = format(DateCheck$IN[2], "%Z")) - if(AggrCheck[[2]] == "month" & - (DateCheck$IN[1] != MustStartMonth | - DateCheck$IN[2] != MustEndMonth)){ + days_in_month(DateCheck$IN[2]), + sep = "-" + ), + "24:00:00" + ), tz = format(DateCheck$IN[2], "%Z")) + if (AggrCheck[[2]] == "month" && (DateCheck$IN[1] != MustStartMonth || DateCheck$IN[2] != MustEndMonth)) { stop("You have specified (multi-)monthly temporal aggregation but are querying a time window which does not start at the first day of a month at 00:00 and/or does not terminate on the last day of a month at 24:00. Please ensure that you set the argument DateStart and DateStop accordingly.") } MustStartYear <- as.POSIXct(paste( paste(format(DateCheck$IN[1], "%Y"), "01-01", sep = "-"), - "00:00:00"), tz = format(DateCheck$IN[2], "%Z") - ) + "00:00:00" + ), tz = format(DateCheck$IN[2], "%Z")) MustEndYear <- as.POSIXct(paste( paste(format(DateCheck$IN[2], "%Y"), "12-31", sep = "-"), - "23:00:00"), tz = format(DateCheck$IN[2], "%Z") - ) - if(AggrCheck[[2]] == "year" & - (DateCheck$IN[1] != MustStartYear | - DateCheck$IN[2] != MustEndYear)){ + "23:00:00" + ), tz = format(DateCheck$IN[2], "%Z")) + if (AggrCheck[[2]] == "year" && (DateCheck$IN[1] != MustStartYear || DateCheck$IN[2] != MustEndYear)) { stop("You have specified (multi-)yearly temporal aggregation but are querying a time window which does not start at the first of day of a year at 00:00 and/or does not terminate on the last day of a year at 23:00. Please ensure that you set the argument DateStart and DateStop accordingly.") } } @@ -283,7 +293,8 @@ Meta.Check <- function(DataSet = "reanalysis-era5-land", Type = NA, VariableChec QueryType = Type, QueryVariable = VariableCheck, QueryFormat = QueryFormat, - QueryUnit = Meta.Variables(dataset = DataSet)$Unit[which(Meta.Variables(dataset = DataSet)$CDSname == VariableCheck)]) + QueryUnit = Meta.Variables(dataset = DataSet)$Unit[which(Meta.Variables(dataset = DataSet)$CDSname == VariableCheck)] + ) } ### NETCDF METADATA WRITING AND READING ======================================== @@ -296,32 +307,33 @@ Meta.Check <- function(DataSet = "reanalysis-era5-land", Type = NA, VariableChec #' @param Attrs Named vector of metadata attributes #' @param Write Logical. Whether to write metadata #' @param Read Logical Whether to read metadata +#' @param Compression Integer between 1 to 9. Applied to final .nc file that the function writes to hard drive. Same as compression argument in terra::writeCDF(). Ignored if FileExtension = ".tif". #' #' @importFrom terra writeCDF +#' @importFrom terra rast #' @importFrom ncdf4 nc_open #' @importFrom ncdf4 ncatt_put #' @importFrom ncdf4 nc_close #' @importFrom ncdf4 ncatt_get #' @importFrom terra metags -#' @importFrom lubridate days_in_month #' #' @return A SpatRaster with metadata #' -Meta.NC <- function(NC, FName, Attrs, Write = FALSE, Read = FALSE){ +Meta.NC <- function(NC, FName, Attrs, Write = FALSE, Read = TRUE, Compression = 9) { ## Writing metadata - if(Write){ - writeCDF(x = NC, filename = FName) + if (Write) { + NC <- writeCDF(x = NC, filename = FName, compression = Compression) nc <- nc_open(FName, write = TRUE) - for(name in names(Attrs)) { + for (name in names(Attrs)) { ncatt_put(nc, 0, name, Attrs[[name]]) } nc_close(nc) } ## Reading metadata - if(Read){ + if (Read) { nc <- nc_open(FName) # Retrieve custom metadata - Meta <- lapply(names(Attrs), FUN = function(name){ + Meta <- lapply(names(Attrs), FUN = function(name) { ncatt_get(nc, 0, name)$value }) # Close the NetCDF file diff --git a/R/Plotting.R b/R/Plotting.R new file mode 100644 index 0000000..c3b9b40 --- /dev/null +++ b/R/Plotting.R @@ -0,0 +1,260 @@ +### Raw Data ========================================================== +#' Visualise raster data and overlay sf polygons if desired. +#' +#' Use the ggplot2 plotting engine to easily create visualisations of raster data - like the ones obtained using CDownloadS(...) - and overlay sf polygon data if desired. +#' +#' @param SpatRast SpatRast object to visualise. +#' @param SF Optional. SF object which to overlay. +#' @param Dates Optional. Character vector of labels to apply to each layer of the SpatRast. By default, the content of the terra::time() field of the supplied SpatRast object. +#' @param Legend Colour label legend. +#' @param COL Colour palette. +#' +#' @importFrom viridis inferno +#' @importFrom terra time +#' @importFrom tidyr gather +#' @importFrom ggplot2 ggplot +#' @importFrom ggplot2 aes +#' @importFrom ggplot2 geom_raster +#' @importFrom ggplot2 theme_bw +#' @importFrom ggplot2 facet_wrap +#' @importFrom ggplot2 labs +#' @importFrom ggplot2 scale_fill_gradientn +#' @importFrom ggplot2 theme +#' @importFrom ggplot2 unit +#' @importFrom ggplot2 geom_sf +#' +#' @return A ggplot2 object visualising a raster. +#' +#' @seealso \code{\link{CDownloadS}}. +#' +#' @examples +#' SpatRast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR"))[[1:2]] +#' data("Jotunheimen_poly") +#' SF <- Jotunheimen_poly +#' Plot.SpatRast(SpatRast = SpatRast, SF = SF) +#' +#' @export +Plot.SpatRast <- function(SpatRast, SF, Dates, Legend = "Air Temperature [K]", COL = viridis::inferno(100)) { + if (missing(Dates)) { + Dates <- as.character(terra::time(SpatRast)) + } + Raw_df <- as.data.frame(SpatRast, xy = TRUE) # turn raster into dataframe + colnames(Raw_df)[c(-1, -2)] <- Dates # set colnames + Raw_df <- tidyr::gather(data = Raw_df, key = Values, value = "value", colnames(Raw_df)[c(-1, -2)]) # make ggplot-ready + Raw_plot <- ggplot() + # create plot + geom_raster(data = Raw_df, aes(x = x, y = y, fill = value)) + # plot the covariate data + theme_bw() + + facet_wrap(~Values) + + labs(x = "Longitude", y = "Latitude") + # make plot more readable + scale_fill_gradientn(name = Legend, colours = COL, na.value = "transparent") + # add colour and legend + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + # reduce margins (for fusing of plots) + theme(legend.key.size = unit(1.5, "cm")) + if (!missing(SF)) { # if a shape has been designated + Raw_plot <- Raw_plot + geom_sf(data = SF, colour = "black", fill = "NA") # add shape + } + return(Raw_plot) +} # export the plot + +### Covariate Data ========================================================== +#' Visualise covariate raster data and overlay sf polygons if desired. +#' +#' Use the ggplot2 plotting engine to easily create visualisations of covariate raster data - like the ones obtained using CovariateSetup(...) - and overlay sf polygon data if desired. +#' +#' @param Covariates List of length 2. Containing training resolution covariate SpatRast in slot 1 and target resolution covariate SpatRast in slot 2. +#' @param SF Optional. SF object which to overlay. +#' @param COL Colour palette. +#' +#' @importFrom viridis cividis +#' @importFrom terra nlyr +#' @importFrom tidyr gather +#' @importFrom ggplot2 ggplot +#' @importFrom ggplot2 aes +#' @importFrom ggplot2 geom_raster +#' @importFrom ggplot2 theme_bw +#' @importFrom ggplot2 facet_wrap +#' @importFrom ggplot2 labs +#' @importFrom ggplot2 scale_fill_gradientn +#' @importFrom ggplot2 theme +#' @importFrom ggplot2 unit +#' @importFrom ggplot2 geom_sf +#' @importFrom cowplot plot_grid +#' +#' @return A ggplot2 object visualising a raster. +#' +#' @seealso \code{\link{CovariateSetup}}. +#' +#' @examples +#' Cov_train <- terra::rast(system.file("extdata", "Covariates_Train.nc", package = "KrigR")) +#' Cov_target <- terra::rast(system.file("extdata", "Covariates_Target.nc", package = "KrigR")) +#' names(Cov_train) <- names(Cov_target) <- "GMTED2010 [m]" +#' Covariates <- list(Training = Cov_train, Target = Cov_target) +#' data("Jotunheimen_poly") +#' SF <- Jotunheimen_poly +#' Plot.Covariates(Covariates = Covariates, SF = SF) +#' +#' @export +Plot.Covariates <- function(Covariates, SF, COL = viridis::cividis(100)) { + Plots_ls <- as.list(rep(NA, nlyr(Covariates[[1]]))) # create as many plots as there are covariates variables + for (Variable in 1:nlyr(Covariates[[1]])) { # loop over all covariate variables + Covariates_Iter <- list(Covariates[[1]][[Variable]], Covariates[[2]][[Variable]]) # extract the data for this variable + Covariates_Iter <- lapply(Covariates_Iter, FUN = function(x) { + Cov_df <- as.data.frame(x, xy = TRUE) # turn raster into dataframe + gather(data = Cov_df, key = Values, value = "value", colnames(Cov_df)[c(-1, -2)]) # make ggplot-ready + }) + Covariates_Iter[[1]][, 3] <- "Native" + Covariates_Iter[[2]][, 3] <- "Target" + Cov_df <- do.call(rbind, Covariates_Iter) + Plots_ls[[Variable]] <- ggplot() + # create plot + geom_raster(data = Cov_df, aes(x = x, y = y, fill = value)) + # plot the covariate data + theme_bw() + + facet_wrap(~Values) + + labs(x = "Longitude", y = "Latitude") + # make plot more readable + scale_fill_gradientn(name = names(Covariates[[1]][[Variable]]), colours = COL, na.value = "transparent") + # add colour and legend + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + # reduce margins (for fusing of plots) + theme(legend.key.size = unit(1.5, "cm")) + if (!missing(SF)) { # if a shape has been designated + Plots_ls[[Variable]] <- Plots_ls[[Variable]] + geom_sf(data = SF, colour = "black", fill = "NA") # add shape + } + # } # end of resolution loop + } # end of variable loop + if (nlyr(Covariates[[1]]) > 1) { + ggPlot <- plot_grid(plotlist = Plots_ls, ncol = 1, labels = "AUTO") # fuse the plots into one big plot + return(ggPlot) + } else { + return(Plots_ls[[1]]) + } +} # export the plot + +### Kriged Data ========================================================== +#' Visualise kriged raster data and associated uncertainties and overlay sf polygons if desired. +#' +#' Use the ggplot2 plotting engine to easily create visualisations of kriged raster data and associated uncertainties raster data - like the ones obtained using Kriging(...) - and overlay sf polygon data if desired. +#' +#' @param Krigs List of length 2. Containing kriging prediction SpatRast in slot 1 and kriging standard error SpatRast in slot 2. +#' @param SF Optional. SF object which to overlay. +#' @param Dates Optional. Character vector of labels to apply to each layer of the SpatRast. By default, the content of the terra::time() field of the supplied SpatRast objects in list. +#' @param Legend Colour label legend. +#' +#' @importFrom viridis inferno +#' @importFrom viridis viridis +#' @importFrom terra time +#' @importFrom tidyr gather +#' @importFrom ggplot2 ggplot +#' @importFrom ggplot2 aes +#' @importFrom ggplot2 geom_raster +#' @importFrom ggplot2 theme_bw +#' @importFrom ggplot2 facet_wrap +#' @importFrom ggplot2 labs +#' @importFrom ggplot2 scale_fill_gradientn +#' @importFrom ggplot2 theme +#' @importFrom ggplot2 unit +#' @importFrom ggplot2 geom_sf +#' @importFrom cowplot plot_grid +#' +#' @return A ggplot2 object visualising a raster. +#' +#' @seealso \code{\link{Kriging}}. +#' +#' @examples +#' CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +#' Cov_train <- terra::rast(system.file("extdata", "Covariates_Train.nc", package = "KrigR")) +#' Cov_target <- terra::rast(system.file("extdata", "Covariates_Target.nc", package = "KrigR")) +#' names(Cov_train) <- names(Cov_target) <- "GMTED2010 [m]" +#' +#' ### kriging itself +#' ExtentKrig <- Kriging( +#' Data = CDS_rast[[1:2]], +#' Covariates_training = Cov_train, +#' Covariates_target = Cov_target, +#' Equation = "GMTED2010", +#' Cores = 2, +#' FileName = "KrigTest1", +#' FileExtension = ".nc", +#' Keep_Temporary = TRUE, +#' nmax = 40, +#' verbose = TRUE +#' ) +#' +#' Plot.Kriged(Krigs = ExtentKrig) +#' +#' @export +Plot.Kriged <- function(Krigs, SF, Dates, Legend = "Air Temperature [K]") { + if (missing(Dates)) { + Dates <- as.character(terra::time(Krigs[[1]])) + } + Type_vec <- c("Prediction", "Standard Error") # these are the output types of krigR + Colours_ls <- list(inferno(100), rev(viridis(100))) # we want separate colours for the types + Plots_ls <- as.list(NA, NA) # this list will be filled with the output plots + for (Plot in 1:2) { # loop over both output types + Krig_df <- as.data.frame(Krigs[[Plot]], xy = TRUE) # turn raster into dataframe + colnames(Krig_df)[c(-1, -2)] <- paste(Type_vec[Plot], Dates) # set colnames + Krig_df <- gather(data = Krig_df, key = Values, value = "value", colnames(Krig_df)[c(-1, -2)]) # make ggplot-ready + Plots_ls[[Plot]] <- ggplot() + # create plot + geom_raster(data = Krig_df, aes(x = x, y = y, fill = value)) + # plot the kriged data + facet_wrap(~Values) + # split raster layers up + theme_bw() + + labs(x = "Longitude", y = "Latitude") + # make plot more readable + scale_fill_gradientn(name = Legend, colours = Colours_ls[[Plot]], na.value = "transparent") + # add colour and legend + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + # reduce margins (for fusing of plots) + theme(legend.key.size = unit(1, "cm")) + if (!missing(SF)) { # if a shape has been designated + Plots_ls[[Plot]] <- Plots_ls[[Plot]] + geom_sf(data = SF, colour = "black", fill = "NA") # add shape + } + } # end of type-loop + ggPlot <- plot_grid(plotlist = Plots_ls, ncol = 1, labels = "AUTO") # fuse the plots into one big plot + return(ggPlot) +} # export the plot + +### Bioclimatic Data ========================================================== +#' Visualise bioclimatic raster data and overlay sf polygons if desired. +#' +#' Use the ggplot2 plotting engine to easily create visualisations of biolcimatic raster data - like the ones obtained using BioClim(...) - and overlay sf polygon data if desired. +#' +#' @param BioClims SpatRast object to visualise. +#' @param Which Numeric. Which bioclimatic variable(s) to visualise. +#' @param SF Optional. SF object which to overlay. +#' @param Water_Var Optional, character. Name of water availability variable in the bioclimatic variables. +#' @param ncol Number of columns for panel arrangement of plots +#' +#' @importFrom terra nlyr +#' @importFrom viridis inferno +#' @importFrom viridis mako +#' @importFrom cowplot plot_grid +#' +#' @return A ggplot2 object visualising a raster. +#' +#' @seealso \code{\link{BioClim}}. +#' +#' @examples +#' BC_rast <- terra::rast(system.file("extdata", "CN_BC.nc", package = "KrigR")) +#' Plot.BioClim(BioClims = BC_rast, Water_Var = "Soil Moisture (0-7cm)") +#' +#' @export +Plot.BioClim <- function(BioClims, Which = 1:19, SF, Water_Var = "Water Availability", ncol = 3) { + if (missing(SF)) { + SF <- NULL + } + if (nlyr(BioClims) != 19) { + stop("The raster data you supplied does not contain 19 layers. Please supply raster data containing 19 layers corresponding to the 19 bioclimatic variables.") + } + BC_names <- c("Annual Mean Temperature", "Mean Diurnal Range", "Isothermality", "Temperature Seasonality", "Max Temperature of Warmest Month", "Min Temperature of Coldest Month", "Temperature Annual Range (BIO5-BIO6)", "Mean Temperature of Wettest Quarter", "Mean Temperature of Driest Quarter", "Mean Temperature of Warmest Quarter", "Mean Temperature of Coldest Quarter", paste("Annual", Water_Var), paste(Water_Var, "of Wettest Month"), paste(Water_Var, "of Driest Month"), paste(Water_Var, "Seasonality"), paste(Water_Var, "of Wettest Quarter"), paste(Water_Var, "of Driest Quarter"), paste(Water_Var, "of Warmest Quarter"), paste(Water_Var, "of Coldest Quarter")) + BC_names <- paste0("BIO", 1:19, " - ", BC_names) + names(BioClims) <- BC_names + + ToPlot <- BioClims[[Which]] + Colours <- list( + Temperature = inferno(1e3), + Water = mako(1e3) + ) + + Plots_ls <- lapply(1:nlyr(ToPlot), FUN = function(x) { + COL <- ifelse(grepl(Water_Var, names(ToPlot[[x]]), fixed = TRUE), "Water", "Temperature") + if (!is.null(SF)) { + Plot.SpatRast(SpatRast = ToPlot[[x]], SF = SF, Dates = names(ToPlot[[x]]), Legend = "", COL = Colours[[COL]]) + } else { + Plot.SpatRast(SpatRast = ToPlot[[x]], Dates = names(ToPlot[[x]]), Legend = " ", COL = Colours[[COL]]) + } + }) + + cowplot::plot_grid(plotlist = Plots_ls, ncol = ncol) +} # export the plot diff --git a/R/Spatial.R b/R/Spatial.R index 4bdd917..4b9a216 100644 --- a/R/Spatial.R +++ b/R/Spatial.R @@ -14,9 +14,11 @@ #' Make.SpatialPoints(Mountains_df) #' #' @export -Make.SpatialPoints <- function(USER_df){ +Make.SpatialPoints <- function(USER_df) { USER_df <- data.frame(USER_df) ## attempt to catch tibbles or data.tables - if(sum(c("Lat", "Lon") %in% colnames(USER_df)) != 2){stop("Please provide your geo-locations with a Lat and a Lon column (named exactly like such).")} + if (sum(c("Lat", "Lon") %in% colnames(USER_df)) != 2) { + stop("Please provide your geo-locations with a Lat and a Lon column (named exactly like such).") + } st_as_sf(USER_df, coords = c("Lon", "Lat"), remove = FALSE) } ### EXTENT CHECKING ============================================================ @@ -35,7 +37,7 @@ Make.SpatialPoints <- function(USER_df){ #' @return A list containg (1) a terra/sf object and (2) the corresponding SpatExtent object. #' #' @examples -#' ## raster +#' ## raster #' Check.Ext(raster::extent(c(9.87, 15.03, 49.89, 53.06))) #' ## terra #' Check.Ext(terra::ext(c(9.87, 15.03, 49.89, 53.06))) @@ -43,13 +45,13 @@ Make.SpatialPoints <- function(USER_df){ #' set.seed(42) #' nb_pt <- 10 #' dd <- data.frame(x = runif(nb_pt, 9.87, 15.03), y = runif(nb_pt, 49.89, 53.06), val = rnorm(nb_pt)) -#' sf <- sf::st_as_sf(dd, coords = c("x","y")) +#' sf <- sf::st_as_sf(dd, coords = c("x", "y")) #' Check.Ext(sf) #' ## sp #' Check.Ext(as(sf, "Spatial")) #' #' @export -Ext.Check <- function(USER_ext){ +Ext.Check <- function(USER_ext) { ## find package where USER_ext class originates class_name <- class(USER_ext) class_def <- getClass(class_name) @@ -57,31 +59,33 @@ Ext.Check <- function(USER_ext){ ## sanity check if USER_ext is supported SupportedPackages <- c("raster", "terra", "sf", "sp") - if(!(package_name %in% SupportedPackages)){ + if (!(package_name %in% SupportedPackages)) { stop("Please specify the Extent argument as an object defined either with classes found in the raster, terra, or sf packages") } ## Transform into SpatExtent class - if(package_name == "raster"){ + if (package_name == "raster") { OUT_spatialobj <- rast(USER_ext) OUT_ext <- ext(USER_ext) } - if(package_name == "terra" | package_name == "sf"){ + if (package_name == "terra" || package_name == "sf") { OUT_spatialobj <- USER_ext - if(class_name[1] == "sfc_MULTIPOLYGON"){ + if (class_name[1] == "sfc_MULTIPOLYGON") { OUT_ext <- ext(st_bbox(USER_ext)) - }else{ + } else { OUT_ext <- ext(USER_ext) } } - if(package_name == "sp"){ + if (package_name == "sp") { OUT_spatialobj <- st_as_sf(USER_ext) OUT_ext <- ext(OUT_spatialobj) } ## Round digits and return - OUT_list = list(SpatialObj = OUT_spatialobj, - Ext = round(OUT_ext, 3)) + OUT_list <- list( + SpatialObj = OUT_spatialobj, + Ext = round(OUT_ext, 3) + ) return(OUT_list) } @@ -105,7 +109,7 @@ Ext.Check <- function(USER_ext){ #' Buffer.pts(User_pts, USER_buffer = 0.5) #' #' @export -Buffer.pts <- function(USER_pts, USER_buffer = .5){ +Buffer.pts <- function(USER_pts, USER_buffer = .5) { st_as_sf(st_union(st_buffer(USER_pts, USER_buffer, endCapStyle = "SQUARE"))) } @@ -130,17 +134,17 @@ Buffer.pts <- function(USER_pts, USER_buffer = .5){ #' Mask.Shape(Jotunheimen_ras, Jotunheimen_poly) #' #' @export -Handle.Spatial <- function(BASE, Shape){ - +Handle.Spatial <- function(BASE, Shape) { ## splitting by rasterlayers if necessary to avoid error reported in https://github.com/rspatial/terra/issues/1556 - if(terra::nlyr(BASE) > 65535){ - Indices <- ceiling((1:terra::nlyr(BASE))/2e4) + if (terra::nlyr(BASE) > 65535) { + Indices <- ceiling((1:terra::nlyr(BASE)) / 2e4) r_ls <- terra::split(x = BASE, f = Indices) - ret_ls <- pblapply(r_ls, FUN = function(BASE_iter){ + ret_ls <- pblapply(r_ls, FUN = function(BASE_iter) { ret_rast <- crop(BASE_iter, ext(Shape)) - if(class(Shape)[1] == "sf"){ + if (class(Shape)[1] == "sf") { ret_rast <- mask(ret_rast, Shape, touches = TRUE) } + ret_rast }) ret_rast <- do.call(c, ret_ls) return(ret_rast) @@ -148,7 +152,7 @@ Handle.Spatial <- function(BASE, Shape){ ## regular cropping and masking for SPatRasters not exceeding layer limit ret_rast <- crop(BASE, ext(Shape)) - if(class(Shape)[1] == "sf"){ + if (class(Shape)[1] == "sf") { ret_rast <- mask(ret_rast, Shape, touches = TRUE) } return(ret_rast) diff --git a/R/Temporal.R b/R/Temporal.R index fc07162..142104d 100644 --- a/R/Temporal.R +++ b/R/Temporal.R @@ -14,9 +14,11 @@ #' Dates_df #' #' @export -Make.UTC <- function(DatesVec = NULL){ +Make.UTC <- function(DatesVec = NULL) { data.frame(IN = DatesVec, - UTC = do.call(c ,lapply(DatesVec, FUN = function(x){as.POSIXct(x, tz = "UTC")})) + UTC = do.call(c, lapply(DatesVec, FUN = function(x) { + as.POSIXct(x, tz = "UTC") + })) ) } ### QUERY SEPARATING INTO TIME WINDOWS ========================================= @@ -50,63 +52,75 @@ Make.UTC <- function(DatesVec = NULL){ #' BaseTStart = as.POSIXct("1950-01-01 00:01", tz = "UTC") #' TChunkSize = 12000) #' -Make.RequestWindows <- function(Dates_df, BaseTResolution, BaseTStep, BaseTStart, TChunkSize, DataSet){ +Make.RequestWindows <- function(Dates_df, BaseTResolution, BaseTStep, BaseTStart, TChunkSize, DataSet) { ## reformat input DateStart <- Dates_df$UTC[1] DateStop <- Dates_df$UTC[2] - if(BaseTResolution == "month"){ + if (BaseTResolution == "month") { DateStart <- as.POSIXct(paste0(format(DateStart, "%Y-"), "01-01 00:00"), tz = "UTC") # ensure that first of first month in first queried year is used for sequence creation to avoid month skips DateStop <- as.POSIXct(paste0(format(DateStop, "%Y-"), "12-31 23:00"), tz = "UTC") # ensure that last of last month in last queried year is used for sequence creation to avoid month skips } ## checking chunksize specification - if((TChunkSize/BaseTStep)%%1!=0){stop("Please specify a TChunkSize (currently = ", TChunkSize, - ") that is a multiple of the base temporal resolution of the data you queried from CDS (curently = ", BaseTStep, ").")} + if ((TChunkSize / BaseTStep) %% 1 != 0) { + stop( + "Please specify a TChunkSize (currently = ", TChunkSize, + ") that is a multiple of the base temporal resolution of the data you queried from CDS (curently = ", BaseTStep, ")." + ) + } ## checking alignment of queried data with raw data - if(BaseTResolution == "hour" & BaseTStep != 24){ + if (BaseTResolution == "hour" && BaseTStep != 24) { # when we are pulling from non-1-hourly records, check whether specified start-date aligns with date layers in raw data - StartCheck <- difftime(DateStart, Meta.QuickFacts(dataset = DataSet)$TStart, units = "hour")/BaseTStep - EndCheck <- difftime(DateStop, Meta.QuickFacts(dataset = DataSet)$TStart, units = "hour")/BaseTStep - AlignCheck <- (as.numeric(StartCheck)%%1==0 | as.numeric(EndCheck)%%1==0) + StartCheck <- difftime(DateStart, Meta.QuickFacts(dataset = DataSet)$TStart, units = "hour") / BaseTStep + EndCheck <- difftime(DateStop, Meta.QuickFacts(dataset = DataSet)$TStart, units = "hour") / BaseTStep + AlignCheck <- (as.numeric(StartCheck) %% 1 == 0 || as.numeric(EndCheck) %% 1 == 0) } ## making query time call - if(!(BaseTResolution %in% c("hour", "month"))){stop("Non-hour or -month base resolutions not supported yet")} - if(BaseTResolution == "hour"){ - if(BaseTStep == 24){ - QueryTimes <- str_pad(str_c(0:23,"00",sep=":"), 5,"left","0") ## this is used for telling CDS which layers we want per day - }else{ + if (!(BaseTResolution %in% c("hour", "month"))) { + stop("Non-hour or -month base resolutions not supported yet") + } + if (BaseTResolution == "hour") { + if (BaseTStep == 24) { + QueryTimes <- str_pad(str_c(0:23, "00", sep = ":"), 5, "left", "0") ## this is used for telling CDS which layers we want per day + } else { QueryTimes <- str_pad(str_c( - seq(from = as.numeric(format(Meta.QuickFacts(dataset = DataSet)$TStart, "%H")), - to = 23, - by = 24/BaseTStep) - ,"00",sep=":"), 5,"left","0") ## this is used for telling CDS which layers we want per day, relevant for ensemble_mean and ensemble_spread for example, which are recorded at 3-hour intervals starting at 00:00 per day + seq( + from = as.numeric(format(Meta.QuickFacts(dataset = DataSet)$TStart, "%H")), + to = 23, + by = 24 / BaseTStep + ), + "00", + sep = ":" + ), 5, "left", "0") ## this is used for telling CDS which layers we want per day, relevant for ensemble_mean and ensemble_spread for example, which are recorded at 3-hour intervals starting at 00:00 per day } } - if(BaseTResolution == "month"){ + if (BaseTResolution == "month") { QueryTimes <- "00:00" ## this is used for telling CDS which layers we want per day } ## check alignment with non-1-BaseTStep data products - if(exists("AlignCheck")){ - if(!AlignCheck){ - stop("You have specified download of a data set whose raw layers are provided at a temporal resolution = ", BaseTResolution, " at intervalstime steps = ", BaseTStep, ".", - "\n Either one or both of the the time-window defining dates (DateStart and DateStop arguments) you have specified, once converted to UTC (", DateStart, " and ", DateStop, ") do not align with the structure of the raw data which requires querying of data to start and terminate at any of the following UTC hours of the day: ", paste(QueryTimes, collapse = "; "), ". Please adjust your date specification accordingly.") + if (exists("AlignCheck")) { + if (!AlignCheck) { + stop( + "You have specified download of a data set whose raw layers are provided at a temporal resolution = ", BaseTResolution, " at intervalstime steps = ", BaseTStep, ".", + "\n Either one or both of the the time-window defining dates (DateStart and DateStop arguments) you have specified, once converted to UTC (", DateStart, " and ", DateStop, ") do not align with the structure of the raw data which requires querying of data to start and terminate at any of the following UTC hours of the day: ", paste(QueryTimes, collapse = "; "), ". Please adjust your date specification accordingly." + ) } } ## making request ranges - if(BaseTResolution == "month"){ + if (BaseTResolution == "month") { BaseTStep <- 1 # do not repeat each month, hence set this to 1 } T_RequestRange <- seq(from = DateStart, to = DateStop, by = BaseTResolution) - # T_RequestRange <- as.POSIXct(paste(rep(unique(T_RequestDates), each = length(QueryTimes)), QueryTimes), tz = "UTC") T_RequestDates <- as.Date(rep(unique(format(T_RequestRange, "%Y-%m-%d")), each = BaseTStep)) - list(QueryTimeWindows = split(T_RequestDates, ceiling(seq_along(T_RequestDates)/TChunkSize)), - QueryTimes = QueryTimes - ) + list( + QueryTimeWindows = split(T_RequestDates, ceiling(seq_along(T_RequestDates) / TChunkSize)), + QueryTimes = QueryTimes + ) } ### BACK-CALCULATION OF CUMULATIVE VARIABLES =================================== @@ -119,50 +133,61 @@ Make.RequestWindows <- function(Dates_df, BaseTResolution, BaseTStep, BaseTStart #' @param BaseResolution Character. Base temporal resolution of data set #' @param BaseStep Numeric. Base time step of data set #' @param TZone Character. Time zone for queried data. +#' @param verbose Logical. Whether to print/message function progress in console or not. #' #' @importFrom terra rast #' @importFrom terra nlyr #' @importFrom terra subset #' @importFrom terra time #' @importFrom lubridate days_in_month +#' @importFrom pbapply pblapply #' #' @return A SpatRaster #' -Temporal.Cumul <- function(CDS_rast, CumulVar, BaseResolution, BaseStep, TZone){ +Temporal.Cumul <- function(CDS_rast, CumulVar, BaseResolution, BaseStep, TZone, verbose = TRUE) { # nolint: cyclocomp_linter. Era5_ras <- CDS_rast - if(CumulVar & BaseResolution == "hour"){ - if(BaseStep != 1){stop("Back-calculation of hourly cumulative variables only supported for 1-hour interval data. The data you have specified reports hourly data in intervals of ", BaseStep, ".")} + if (verbose && CumulVar) { + print("Disaggregation of cumulative records") + } + if (CumulVar && BaseResolution == "hour") { + if (BaseStep != 1) { + stop("Back-calculation of hourly cumulative variables only supported for 1-hour interval data. The data you have specified reports hourly data in intervals of ", BaseStep, ".") + } ## removing non-needed layers - RemovalLyr <- c(1, (nlyr(Era5_ras)-22):nlyr(Era5_ras)) # need to remove first layer and last 23 for backcalculation - Era5_ras <- subset(Era5_ras, RemovalLyr, negate=TRUE) + RemovalLyr <- c(1, (nlyr(Era5_ras) - 22):nlyr(Era5_ras)) # need to remove first layer and last 23 for backcalculation + Era5_ras <- subset(Era5_ras, RemovalLyr, negate = TRUE) ## back-calculation - counter <- 1 - Era5_ls <- as.list(rep(NA, nlyr(Era5_ras))) - names(Era5_ls) <- terra::time(Era5_ras) - for(i in 1:nlyr(Era5_ras)){ - if(counter > 24){counter <- 1} - if(counter == 1){ - Era5_ls[[i]] <- Era5_ras[[i]] - StartI <- i - } - if(counter == 24){ - Era5_ls[[i]] <- Era5_ras[[i]]-sum(rast(Era5_ls[StartI:(StartI+counter-2)])) + #' break apart sequence by UTC days and apply back-calculation per day in pblapply loop, for loop for each hour in each day + DataDays <- ceiling(1:nlyr(Era5_ras) / 24) + DissagDays <- unique(DataDays) + Era5_ls <- pblapply(DissagDays, FUN = function(DissagDay_Iter) { + counter <- 1 + Interior_ras <- Era5_ras[[which(DataDays == DissagDay_Iter)]] + Interior_ls <- as.list(rep(NA, nlyr(Interior_ras))) + names(Interior_ls) <- terra::time(Interior_ras) + for (i in 1:nlyr(Interior_ras)) { + if (counter == 1) { + Interior_ls[[i]] <- Interior_ras[[i]] + } + if (counter == 24) { + Interior_ls[[i]] <- Interior_ras[[i]] - sum(rast(Interior_ls[1:(1 + counter - 2)])) + } + if (counter != 24 & counter != 1) { + Interior_ls[[i]] <- Interior_ras[[i + 1]] - Interior_ras[[i]] + } + counter <- counter + 1 } - if(counter != 24 & counter != 1){ - Era5_ls[[i]] <- Era5_ras[[i+1]] - Era5_ras[[i]] - } - counter <- counter + 1 - } + rast(Interior_ls) + }) ## finishing off object Ret_ras <- rast(Era5_ls) - terra::time(Ret_ras) <- as.POSIXct(terra::time(Era5_ras), tz = TZone) - 60*60 # back-dating to be in-line with regular specifications Era5_ras <- Ret_ras warning("You toggled on the CumulVar option in the function call. Hourly records have been converted from cumulative aggregates to individual hourly records.") } ## multiply by number of days per month - if(CumulVar & BaseResolution == "month"){ + if (CumulVar && BaseResolution == "month") { Days_in_Month_vec <- days_in_month(terra::time(CDS_rast)) - if(grepl("ensemble_members", Type)){ + if (grepl("ensemble_members", Type)) { Days_in_Month_vec <- rep(Days_in_Month_vec, each = 10) } Era5_ras <- Era5_ras * Days_in_Month_vec @@ -185,6 +210,7 @@ Temporal.Cumul <- function(CDS_rast, CumulVar, BaseResolution, BaseStep, TZone){ #' @param Cores Numeric. Number of cores for parallel processing #' @param QueryTargetSteps Character. Target resolution steps #' @param TZone Character. Time zone for queried data. +#' @param verbose Logical. Whether to print/message function progress in console or not. #' #' @importFrom terra time #' @importFrom terra tapp @@ -193,63 +219,64 @@ Temporal.Cumul <- function(CDS_rast, CumulVar, BaseResolution, BaseStep, TZone){ #' @return A SpatRaster #' Temporal.Aggr <- function(CDS_rast, BaseResolution, BaseStep, - TResolution, TStep, FUN, Cores, QueryTargetSteps, TZone){ - if(BaseResolution == TResolution & BaseStep == TStep){ + TResolution, TStep, FUN, Cores, QueryTargetSteps, TZone, verbose = TRUE) { + if (verbose) { + print("Temporal Aggregation") + } + if (BaseResolution == TResolution && BaseStep == TStep) { Final_rast <- CDS_rast # no temporal aggregation needed - }else{ - - - TimeDiff <- sapply(terra::time(CDS_rast), FUN = function(xDate){ - length(seq(from = terra::time(CDS_rast)[1], - to = xDate, - by = TResolution))-1 + } else { + TimeDiff <- sapply(terra::time(CDS_rast), FUN = function(xDate) { + length(seq( + from = terra::time(CDS_rast)[1], + to = xDate, + by = TResolution + )) - 1 }) - AggrIndex <- floor(TimeDiff/TStep)+1 + AggrIndex <- floor(TimeDiff / TStep) + 1 Form <- substr(TResolution, 1, 1) Form <- ifelse(Form %in% c("h", "y"), toupper(Form), Form) LayerFormat <- format(terra::time(CDS_rast), paste0("%", Form)) - # LayerMatches <- match(LayerFormat, QueryTargetSteps) - # - # if(grepl("Factor =", QueryTargetSteps)){ # happens for hourly ensemble data - # Factor <- as.numeric(sub("Ensembling at base resolution, Factor = ", "", QueryTargetSteps)) - # AggrIndex <- rep(seq(from = 1, to = length(LayerFormat)/(Factor)), each = Factor) - # }else{ - # AggrIndex <- ceiling(LayerMatches/TStep) - # } - if(length(unique(AggrIndex)) == 1){ ## this is to avoid a warning message thrown by terra - Final_rast <- app(x = CDS_rast, - cores = Cores, - fun = FUN) - }else{ - Final_rast <- tapp(x = CDS_rast, - index = AggrIndex, - cores = Cores, - fun = FUN) + if (length(unique(AggrIndex)) == 1) { ## this is to avoid a warning message thrown by terra + Final_rast <- app( + x = CDS_rast, + cores = Cores, + fun = FUN + ) + } else { + Final_rast <- tapp( + x = CDS_rast, + index = AggrIndex, + cores = Cores, + fun = FUN + ) } - - - if(TResolution == "year"){ + if (TResolution == "year") { terra::time(Final_rast) <- as.POSIXct( paste0(LayerFormat[!duplicated(AggrIndex)], "-01-01"), - tz = TZone) + tz = TZone + ) } - if(TResolution == "month"){ + if (TResolution == "month") { terra::time(Final_rast) <- as.POSIXct( paste0(format(terra::time(CDS_rast)[!duplicated(AggrIndex)], "%Y-%m"), "-01"), - tz = TZone) + tz = TZone + ) } - if(TResolution == "day"){ + if (TResolution == "day") { terra::time(Final_rast) <- as.POSIXct( format(terra::time(CDS_rast)[!duplicated(AggrIndex)], "%Y-%m-%d"), - tz = TZone) + tz = TZone + ) } - if(TResolution == "hour"){ + if (TResolution == "hour") { terra::time(Final_rast) <- as.POSIXct( terra::time(CDS_rast)[!duplicated(AggrIndex)], - tz = TZone) + tz = TZone + ) } } return(Final_rast) @@ -271,22 +298,20 @@ Temporal.Aggr <- function(CDS_rast, BaseResolution, BaseStep, #' @return Character - target resolution formatted steps in data. #' TemporalAggregation.Check <- function( - QuerySeries, - DateStart, - DateStop, - TResolution, - BaseTResolution, - TStep, - BaseTStep -){ - + QuerySeries, + DateStart, + DateStop, + TResolution, + BaseTResolution, + TStep, + BaseTStep +) { ## check clean division - if(BaseTResolution == TResolution){ ## this comes into play for hourly aggregates of ensemble data - if((TStep/BaseTStep) %%1!=0){ - stop("Your specified time range does not allow for a clean integration of your selected time steps. You specified a time series of raw data with a length of ", length(QueryTargetFormat), " (", BaseTResolution, " intervals of length ", BaseTStep, "). Applying your desired temporal aggregation of ", TResolution, " intervals of length ", TStep, " works out to ", round(TStep/BaseTStep, 3), " intervals. Please fix this so the specified time range can be cleanly divided into aggregation intervals.") + if (BaseTResolution == TResolution) { ## this comes into play for hourly aggregates of ensemble data + if ((TStep / BaseTStep) %% 1 != 0) { + stop("Your specified time range does not allow for a clean integration of your selected time steps. You specified a time series of raw data with a length of ", length(QueryTargetFormat), " (", BaseTResolution, " intervals of length ", BaseTStep, "). Applying your desired temporal aggregation of ", TResolution, " intervals of length ", TStep, " works out to ", round(TStep / BaseTStep, 3), " intervals. Please fix this so the specified time range can be cleanly divided into aggregation intervals.") } - QueryTargetSteps <- paste("Ensembling at base resolution, Factor =", TStep/BaseTStep) - # print("ensemble detected") + QueryTargetSteps <- paste("Ensembling at base resolution, Factor =", TStep / BaseTStep) return(QueryTargetSteps) } @@ -300,7 +325,7 @@ TemporalAggregation.Check <- function( QueryTargetFormat <- format(as.POSIXct(QuerySeries, tz = "UTC"), paste0("%", Form)) QueryTargetSteps <- unique(QueryTargetFormat) - if((length(QueryTargetSteps) / TStep) %%1!=0){ + if ((length(QueryTargetSteps) / TStep) %% 1 != 0) { stop("Your specified time range does not allow for a clean integration of your selected time steps. You specified a time series of raw data with a length of ", length(QueryTargetFormat), " (", BaseTResolution, " intervals of length ", BaseTStep, "). Applying your desired temporal aggregation of ", TResolution, " intervals of length ", TStep, " works out to ", round(length(QueryTargetSteps) / TStep, 3), " intervals. Please fix this so the specified time range can be cleanly divided into aggregation intervals.") } return(QueryTargetSteps) diff --git a/R/data.R b/R/data.R index 159df78..38daa8e 100644 --- a/R/data.R +++ b/R/data.R @@ -13,11 +13,3 @@ #' @format a sf POLYGON #' "Jotunheimen_poly" - -#' #' Daily Air-Temperature Raster across Southern and Central Norway -#' #' -#' #' A SpatRaster object containing five layers of daily mean air-temperature data for the time-period 1995-01-01 to 1995-01-5 sourced from ERA5-Land. -#' # -#' #' @format a SpatRaster -#' #' @references Muñoz Sabater, J. (2019): ERA5-Land hourly data from 1950 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS). DOI: 10.24381/cds.e2161bac (Accessed on 2024-06-24) -#' "CDS_rast" diff --git a/R/krigr.R b/R/krigr.R deleted file mode 100644 index 45f5e74..0000000 --- a/R/krigr.R +++ /dev/null @@ -1,308 +0,0 @@ -#' (multi-core) Kriging -#' -#' This function statistically downscales input data using covariate data and the kriging methodology. The function can be run in two ways: -#' \enumerate{ -#' \item \strong{By Itself}: Use the arguments Data, Covariates_coarse, Covariates_fine when you already have raster files for your data which is to be downscaled as well as covariate raster data. -#' \item \strong{From Scratch}: Use the arguments Variable, Type, DataSet, DateStart, DateStop, TResolution, TStep, Extent, Dir, FileName, API_Key, API_User, and arget_res. By doing so, krigR will call the functions download_ERA() and download_DEM() for one coherent kriging workflow. Note that this process does not work when targetting UERRA data. -#' } -#' Use optional arguments such as Dir, FileName, Keep_Temporary, SingularTry, KrigingEquation and Cores for ease of use, substitution of non-GMTED2010 covariates, and parallel processing. -#' -#' @param Data Raster file which is to be downscaled. -#' @param Covariates_coarse Raster file containing covariates at training resolution. -#' @param Covariates_fine Raster file containing covariates at target resolution. -#' @param KrigingEquation Formula or character string specifying which covariates to use and how. Layer names in Covariates_coarse and Covariates_fine need to match Parameters in this formula. Needs to start with "X ~ ". X can read anything you like. -#' @param Dir Optional. Directory specifying where to place final kriged product. Default is current working directory. -#' @param FileName Optional. A file name for the netcdf produced. Default is a combination parameters in the function call. -#' @param Keep_Temporary Logical, whether to delete individual kriging products of layers in Data after processing. Default is TRUE. -#' @param Cores Numeric. How many cores to use. If you want output to your console during the process, use Cores == 1. Parallel processing is carried out when Cores is bigger than 1. Default is detecting all cores of your machine. -#' @param SingularTry Numeric. How often to try kriging of each layer of the input. This usually gets around issues of singular covariance matrices in the kriging process, but takes some time. Default is 10 -#' @param Variable Optional, calls download_ERA(). ERA5(Land)-contained climate variable. -#' @param PrecipFix Optional. Era5(-land) total precipitation is recorded in cumulative steps per hour from the 00:00 time mark per day. Setting PrecipFix to TRUE converts these into records which represent the total precipitation per hour. Monthly records in Era5(-land) express the average daily total precipitation. Setting this argument to TRUE multiplies monthly records by the number of days per the respective month(s) to get to total precipitation records instead of average. Default is FALSE. -#' @param Type Optional. Whether to download reanalysis ('reanalysis') or ensemble ('ensemble_members', 'ensemble_mean', or 'ensemble_spread') data. Passed on to download_ERA. -#' @param DataSet Optional. Which ERA5 data set to download data from. 'era5' or 'era5-land'. Passed on to download_ERA. -#' @param DateStart Optional. Date ('YYYY-MM-DD') at which to start time series of downloaded data. Passed on to download_ERA. -#' @param DateStop Optional. Date ('YYYY-MM-DD') at which to stop time series of downloaded data. Passed on to download_ERA. -#' @param TResolution Optional. Temporal resolution of final product. hour', 'day', 'month'. Passed on to download_ERA. -#' @param TStep Optional. Which time steps (numeric) to consider for temporal resolution. Passed on to download_ERA. -#' @param FUN Optional. A raster calculation argument as passed to `raster::stackApply()`. This controls what kind of data to obtain for temporal aggregates of reanalysis data. Specify 'mean' (default) for mean values, 'min' for minimum values, and 'max' for maximum values, among others. -#' @param Extent Optional, download data according to rectangular bounding box. specify as extent() object or as a raster, a SpatialPolygonsDataFrame object, or a data.frame object. If Extent is a SpatialPolygonsDataFrame, this will be treated as a shapefile and the output will be cropped and masked to this shapefile. If Extent is a data.frame of geo-referenced point records, it needs to contain Lat and Lon columns as well as a non-repeating ID-column. Passed on to download_ERA and download_DEM. -#' @param Buffer Optional. Identifies how big a rectangular buffer to draw around points if Extent is a data frame of points. Buffer is expressed as centessimal degrees. Passed on to download_ERA and download_DEM. -#' @param ID Optional. Identifies which column in Extent to use for creation of individual buffers if Extent is a data.frame. Passed on to download_ERA and download_DEM. -#' @param Target_res Optional. The target resolution for the kriging step (i.e. which resolution to downscale to). An object as specified/produced by raster::res(). Passed on to download_DEM. -#' @param Source Optional, character. Whether to attempt download from the official USGS data viewer (Source = "USGS") or a static copy of the data set on a private drive (Source = "Drive"). Default is "USGS". Use this if the USGS viewer is unavailable. Passed on to download_DEM. -#' @param API_Key Optional. ECMWF cds API key. Passed on to download_ERA. -#' @param API_User Optional. ECMWF cds user number. Passed on to download_ERA. -#' @param nmax Optional. Controls local kriging. Number of nearest observations to be used kriging of each observation. Default is to use all available (Inf). You can specify as a number (numeric). -#' @param TryDown Optional, numeric. How often to attempt the download of each individual file (if querying data download) that the function queries from the server. This is to circumvent having to restart the entire function when encountering connectivity issues. -#' @param verbose Optional, logical. Whether to report progress of data download (if queried) in the console or not. -#' @param TimeOut Numeric. The timeout for each download in seconds. Default 36000 seconds (10 hours). -#' @param SingularDL Logical. Whether to force download of data in one call to CDS or automatically break download requests into individual monthly downloads. Default is FALSE. -#' @return A list object containing the downscaled data as well as the standard error for downscaling as well as the call to the krigR function, and two NETCDF (.nc) file in the specified directory which are the two data contents of the aforementioned list. A temporary directory is populated with individual NETCDF (.nc) files throughout the runtime of krigR which is deleted upon completion if Keep_Temporary = TRUE and all layers in the Data raster object were kriged successfully. -#' @examples -#' \dontrun{ -#' ## THREE-STEP PROCESS (By Itself) -#' # Downloading ERA5-Land air temperature reanalysis data in 12-hour intervals for 02/01/1995 - 04/01/1995 (DD/MM/YYYY). API User and Key in this example are non-functional. Substitute with your user number and key to run this example. -#' Extent <- extent(c(11.8,15.1,50.1,51.7)) # roughly the extent of Saxony -#' API_User <- "..." -#' API_Key <- "..." -#' State_Raw <- download_ERA( -#' Variable = "2m_temperature", -#' DataSet = "era5-land", -#' DateStart = "1995-01-02", -#' DateStop = "1995-01-04", -#' TResolution = "hour", -#' TStep = 12, -#' Extent = Extent, -#' API_User = API_User, -#' API_Key = API_Key -#' ) -#' State_Raw # a raster brick with 6 layers at resolution of ~0.1° -#' # Downloading GMTED2010-data at resolution and extent obtained by a call to download_ERA and a target resolution of .02. -#' Covs_ls <- download_DEM( -#' Train_ras = State_Raw, -#' Target_res = .02, -#' Keep_Temporary = TRUE -#' ) -#' Covs_ls # a list with two elements: (1) GMTED 2010 data at training resolution, and (2) GMTED 2010 data aggregated as close as possible to a resolution of 0.02 -#' # Kriging the data sets prepared with the previous functions. -#' State_Krig <- krigR( -#' Data = State_Raw, # data we want to krig as a raster object -#' Covariates_coarse = Covs_ls[[1]], # training covariate as a raster object -#' Covariates_fine = Covs_ls[[2]], # target covariate as a raster object -#' ) -#' -#' ## PIPELINE (From Scratch) -#' #' # Downloading ERA5-Land air temperature reanalysis data in 12-hour intervals for 02/01/1995 - 04/01/1995 (DD/MM/YYYY), downloading and preparing GMTED 2010 covariate data, and kriging. API User and Key in this example are non-functional. Substitute with your user number and key to run this example. This example produces the same output as the example above. -#' Extent <- extent(c(11.8,15.1,50.1,51.7)) # roughly the extent of Saxony -#' API_User <- "..." -#' API_Key <- "..." -#' Pipe_Krig <- krigR( -#' Variable = "2m_temperature", -#' Type = "reanalysis", -#' DataSet = "era5-land", -#' DateStart = "1995-01-02", -#' DateStop = "1995-01-04", -#' TResolution = "hour",# -#' TStep = 12, -#' Extent = Extent, -#' API_User = API_User, -#' API_Key = API_Key, -#' Target_res = .02, -#' ) -#' } -#' -#' @export -krigR <- function(Data = NULL, Covariates_coarse = NULL, Covariates_fine = NULL, KrigingEquation = "ERA ~ DEM", Cores = detectCores(), Dir = getwd(), FileName, Keep_Temporary = TRUE, SingularTry = 10, Variable, PrecipFix = FALSE, Type = "reanalysis", DataSet = "era5-land", DateStart, DateStop, TResolution = "month", TStep = 1, FUN = 'mean', Extent, Buffer = 0.5, ID = "ID", API_Key, API_User, Target_res, Source = "USGS", nmax = Inf, TryDown = 10, verbose = TRUE, TimeOut = 36000, SingularDL = FALSE, ...){ - - stop("Function currently deprecated as KrigR undergoes major re-development. Please use the stable release to gain access to this functionality.") - - ## CALL LIST (for storing how the function as called in the output) ---- - if(is.null(Data)){ - Data_Retrieval <- list(Variable = Variable, - Type = Type, - PrecipFix = PrecipFix, - DataSet = DataSet, - DateStart = DateStart, - DateStop = DateStop, - TResolution = TResolution, - TStep = TStep, - Extent = Extent) - }else{ - Data_Retrieval <- "None needed. Data was not queried via krigR function, but supplied by user." - } - ## CLIMATE DATA (call to download_ERA function if no Data set is specified) ---- - if(is.null(Data)){ # data check: if no data has been specified - Data <- download_ERA(Variable = Variable, PrecipFix = PrecipFix, Type = Type, DataSet = DataSet, DateStart = DateStart, DateStop = DateStop, TResolution = TResolution, TStep = TStep, FUN = FUN, Extent = Extent, API_User = API_User, API_Key = API_Key, Dir = Dir, TryDown = TryDown, verbose = verbose, ID = ID, Cores = Cores, TimeOut = TimeOut, SingularDL = SingularDL) - } # end of data check - - ## COVARIATE DATA (call to download_DEM function when no covariates are specified) ---- - if(is.null(Covariates_coarse) & is.null(Covariates_fine)){ # covariate check: if no covariates have been specified - if(class(Extent) == "SpatialPolygonsDataFrame" | class(Extent) == "data.frame"){ # Extent check: if Extent input is a shapefile - Shape <- Extent # save shapefile for use as Shape in masking covariate data - }else{ # if Extent is not a shape, then extent specification is already baked into Data - Shape <- NULL # set Shape to NULL so it is ignored in download_DEM function when masking is applied - } # end of Extent check - Covs_ls <- download_DEM(Train_ras = Data, Target_res = Target_res, Shape = Shape, Buffer = Buffer, ID = ID, Keep_Temporary = Keep_Temporary, Dir = Dir, Source = Source) - Covariates_coarse <- Covs_ls[[1]] # extract coarse covariates from download_DEM output - Covariates_fine <- Covs_ls[[2]] # extract fine covariates from download_DEM output - } # end of covariate check - - ## KRIGING FORMULA (assure that KrigingEquation is a formula object) ---- - KrigingEquation <- as.formula(KrigingEquation) - - ## CALL LIST (for storing how the function as called in the output) ---- - Call_ls <- list(Data = SummarizeRaster(Data), - Covariates_coarse = SummarizeRaster(Covariates_coarse), - Covariates_fine = SummarizeRaster(Covariates_fine), - KrigingEquation = KrigingEquation, - Cores = Cores, - FileName = FileName, - Keep_Temporary = Keep_Temporary, - nmax = nmax, - Data_Retrieval = Data_Retrieval) - - ## SANITY CHECKS (step into check_Krig function to catch most common error messages) ---- - Check_Product <- check_Krig(Data = Data, CovariatesCoarse = Covariates_coarse, CovariatesFine = Covariates_fine, KrigingEquation = KrigingEquation) - KrigingEquation <- Check_Product[[1]] # extract KrigingEquation (this may have changed in check_Krig) - DataSkips <- Check_Product[[2]] # extract which layers to skip due to missing data (this is unlikely to ever come into action) - Terms <- unique(unlist(strsplit(labels(terms(KrigingEquation)), split = ":"))) # identify which layers of data are needed - - ## DATA REFORMATTING (Kriging requires spatially referenced data frames, reformatting from rasters happens here) --- - Origin <- raster::as.data.frame(Covariates_coarse, xy = TRUE) # extract covariate layers - Origin <- Origin[, c(1:2, which(colnames(Origin) %in% Terms))] # retain only columns containing terms - - Target <- raster::as.data.frame(Covariates_fine, xy = TRUE) # extract covariate layers - Target <- Target[, c(1:2, which(colnames(Target) %in% Terms))] # retain only columns containing terms - Target <- na.omit(Target) - suppressWarnings(gridded(Target) <- ~x+y) # establish a gridded data product ready for use in kriging - Target@grid@cellsize[1] <- Target@grid@cellsize[2] # ensure that grid cells are square - - ## SET-UP TEMPORARY DIRECTORY (this is where kriged products of each layer will be saved) ---- - Dir.Temp <- file.path(Dir, paste("Kriging", FileName, sep="_")) - if(!dir.exists(Dir.Temp)){dir.create(Dir.Temp)} - - ## KRIGING SPECIFICATION (this will be parsed and evaluated in parallel and non-parallel evaluations further down) ---- - looptext <- " - OriginK <- cbind(Origin, raster::extract(x = Data[[Iter_Krige]], y = Origin[,1:2], df=TRUE)[, 2]) # combine data of current data layer with training covariate data - OriginK <- na.omit(OriginK) # get rid of NA cells - colnames(OriginK)[length(Terms)+3] <- c(terms(KrigingEquation)[[2]]) # assign column names - suppressWarnings(gridded(OriginK) <- ~x+y) # generate gridded product - OriginK@grid@cellsize[1] <- OriginK@grid@cellsize[2] # ensure that grid cells are square - - Iter_Try = 0 # number of tries set to 0 - kriging_result <- NULL - while(class(kriging_result)[1] != 'autoKrige' & Iter_Try < SingularTry){ # try kriging SingularTry times, this is because of a random process of variogram identification within the automap package that can fail on smaller datasets randomly when it isn't supposed to - try(invisible(capture.output(kriging_result <- autoKrige(formula = KrigingEquation, input_data = OriginK, new_data = Target, nmax = nmax))), silent = TRUE) - Iter_Try <- Iter_Try +1 - } - if(class(kriging_result)[1] != 'autoKrige'){ # give error if kriging fails - message(paste0('Kriging failed for layer ', Iter_Krige, '. Error message produced by autoKrige function: ', geterrmessage())) - } - - ## retransform to raster - try( # try fastest way - this fails with certain edge artefacts in meractor projection and is fixed by using rasterize - Krig_ras <- raster(x = kriging_result$krige_output, layer = 1), # extract raster from kriging product - silent = TRUE - ) - try( - Var_ras <- raster(x = kriging_result$krige_output, layer = 3), # extract raster from kriging product - silent = TRUE - ) - if(!exists('Krig_ras') & !exists('Var_ras')){ - Krig_ras <- rasterize(x = kriging_result$krige_output, y = Covariates_fine[[1]])[[2]] # extract raster from kriging product - Var_ras <- rasterize(x = kriging_result$krige_output, y = Covariates_fine)[[4]] # extract raster from kriging product - } - crs(Krig_ras) <- crs(Data) # setting the crs according to the data - crs(Var_ras) <- crs(Data) # setting the crs according to the data - - if(Cores == 1){ - Ras_Krig[[Iter_Krige]] <- Krig_ras - Ras_Var[[Iter_Krige]] <- Var_ras - } # stack kriged raster into raster list if non-parallel computing - - terra::writeCDF(x = as(brick(Krig_ras), 'SpatRaster'), filename = file.path(Dir.Temp, paste0(str_pad(Iter_Krige,4,'left','0'), '_data.nc')), overwrite = TRUE) - # writeRaster(x = Krig_ras, filename = file.path(Dir.Temp, paste0(str_pad(Iter_Krige,4,'left','0'), '_data.nc')), overwrite = TRUE, format='CDF') # save kriged raster to temporary directory - terra::writeCDF(x = as(brick(Var_ras), 'SpatRaster'), filename = file.path(Dir.Temp, paste0(str_pad(Iter_Krige,4,'left','0'), '_SE.nc')), overwrite = TRUE) - # writeRaster(x = Var_ras, filename = file.path(Dir.Temp, paste0(str_pad(Iter_Krige,4,'left','0'), '_SE.nc')), overwrite = TRUE, format='CDF') # save kriged raster to temporary directory - - if(Cores == 1){ # core check: if processing non-parallel - if(Count_Krige == 1){ # count check: if this was the first actual computation - T_End <- Sys.time() # record time at which kriging was done for current layer - Duration <- as.numeric(T_End)-as.numeric(T_Begin) # calculate how long it took to krig on layer - message(paste('Kriging of remaining ', nlayers(Data)-Iter_Krige, ' data layers should finish around: ', as.POSIXlt(T_Begin + Duration*nlayers(Data), tz = Sys.timezone(location=TRUE)), sep='')) # console output with estimate of when the kriging should be done - ProgBar <- txtProgressBar(min = 0, max = nlayers(Data), style = 3) # create progress bar when non-parallel processing - Count_Krige <- Count_Krige + 1 # raise count by one so the stimator isn't called again - } # end of count check - setTxtProgressBar(ProgBar, Iter_Krige) # update progress bar with number of current layer - } # end of core check - " - - ## KRIGING PREPARATION (establishing objects which the kriging refers to) ---- - Ras_Krig <- as.list(rep(NA, nlayers(Data))) # establish an empty list which will be filled with kriged layers - Ras_Var <- as.list(rep(NA, nlayers(Data))) # establish an empty list which will be filled with kriged layers - - if(verbose){message("Commencing Kriging")} - ## DATA SKIPS (if certain layers in the data are empty and need to be skipped, this is handled here) --- - if(!is.null(DataSkips)){ # Skip check: if layers need to be skipped - for(Iter_Skip in DataSkips){ # Skip loop: loop over all layers that need to be skipped - Ras_Krig[[Iter_Skip]] <- Data[[Iter_Skip]] # add raw data (which should be empty) to list - terra::writeCDF(x = as(brick(Ras_Krig[[Iter_Skip]]), 'SpatRaster'), filename = file.path(Dir.Temp, str_pad(Iter_Skip,4,'left','0')), overwrite = TRUE) - # writeRaster(x = Ras_Krig[[Iter_Skip]], filename = file.path(Dir.Temp, str_pad(Iter_Skip,4,'left','0')), overwrite = TRUE, format = 'CDF') # save raw layer to temporary directory, needed for loading back in when parallel processing - } # end of Skip loop - Layers_vec <- 1:nlayers(Data) # identify vector of all layers in data - Compute_Layers <- Layers_vec[which(!Layers_vec %in% DataSkips)] # identify which layers can actually be computed on - }else{ # if we don't need to skip any layers - Compute_Layers <- 1:nlayers(Data) # set computing layers to all layers in data - } # end of Skip check - - - ## ACTUAL KRIGING (carry out kriging according to user specifications either in parallel or on a single core) ---- - if(Cores > 1){ # Cores check: if parallel processing has been specified - ### PARALLEL KRIGING --- - ForeachObjects <- c("Dir.Temp", "Cores", "Data", "KrigingEquation", "Origin", "Target", "Covariates_coarse", "Covariates_fine", "Terms", "SingularTry", "nmax") # objects which are needed for each kriging run and are thus handed to each cluster unit - pb <- txtProgressBar(max = length(Compute_Layers), style = 3) - progress <- function(n){setTxtProgressBar(pb, n)} - opts <- list(progress = progress) - cl <- makeCluster(Cores) # Assuming Cores node cluster - registerDoSNOW(cl) # registering cores - foreach(Iter_Krige = Compute_Layers, # kriging loop over all layers in Data, with condition (%:% when(...)) to only run if current layer is not present in Dir.Temp yet - .packages = c("raster", "stringr", "automap", "ncdf4", "rgdal", "terra"), # import packages necessary to each itteration - .export = ForeachObjects, - .options.snow = opts) %:% when(!paste0(str_pad(Iter_Krige,4,"left","0"), '_data.nc') %in% list.files(Dir.Temp)) %dopar% { # parallel kriging loop - Ras_Krig <- eval(parse(text=looptext)) # evaluate the kriging specification per cluster unit per layer - } # end of parallel kriging loop - close(pb) - stopCluster(cl) # close down cluster - Files_krig <- list.files(Dir.Temp)[grep(pattern = "_data.nc", x = list.files(Dir.Temp))] - Files_var <- list.files(Dir.Temp)[grep(pattern = "_SE.nc", x = list.files(Dir.Temp))] - for(Iter_Load in 1:length(Files_krig)){ # load loop: load data from temporary files in Dir.Temp - Ras_Krig[[Iter_Load]] <- raster(file.path(Dir.Temp, Files_krig[Iter_Load])) # load current temporary file and write contents to list of rasters - Ras_Var[[Iter_Load]] <- raster(file.path(Dir.Temp, Files_var[Iter_Load])) # load current temporary file and write contents to list of rasters - } # end of load loop - }else{ # if non-parallel processing has been specified - ### NON-PARALLEL KRIGING --- - Count_Krige <- 1 # Establish count variable which is targeted in kriging specification text for producing an estimator - for(Iter_Krige in Compute_Layers){ # non-parallel kriging loop over all layers in Data - if(paste0(str_pad(Iter_Krige,4,'left','0'), '_data.nc') %in% list.files(Dir.Temp)){ # file check: if this file has already been produced - Ras_Krig[[Iter_Krige]] <- raster(file.path(Dir.Temp, paste0(str_pad(Iter_Krige,4,'left','0'), '_data.nc'))) # load already produced kriged file and save it to list of rasters - Ras_Var[[Iter_Krige]] <- raster(file.path(Dir.Temp, paste0(str_pad(Iter_Krige,4,'left','0'), '_SE.nc'))) - if(!exists("ProgBar")){ProgBar <- txtProgressBar(min = 0, max = nlayers(Data), style = 3)} # create progress bar when non-parallel processing} - setTxtProgressBar(ProgBar, Iter_Krige) # update progress bar - next() # jump to next layer - } # end of file check - T_Begin <- Sys.time() # record system time when layer kriging starts - eval(parse(text=looptext)) # evaluate the kriging specification per layer - } # end of non-parallel kriging loop - } # end of Cores check - - ## SAVING FINAL PRODUCT ---- - if(is.null(DataSkips)){ # Skip check: if no layers needed to be skipped - # convert list of kriged layers in actual rasterbrick of kriged layers - names(Ras_Krig) <- names(Data) - if(class(Ras_Krig) != "RasterBrick"){Ras_Krig <- brick(Ras_Krig)} - Krig_terra <- as(Ras_Krig, "SpatRaster") - names(Krig_terra) <- names(Data) - terra::writeCDF(x = Krig_terra, filename = file.path(Dir, paste0(FileName, ".nc")), overwrite = TRUE) - # writeRaster(x = Ras_Krig, filename = file.path(Dir, FileName), overwrite = TRUE, format="CDF") # save final product as raster - # convert list of kriged layers in actual rasterbrick of kriged layers - names(Ras_Var) <- names(Data) - if(class(Ras_Var) != "RasterBrick"){Ras_Var <- brick(Ras_Var)} - Var_terra <- as(Ras_Var, "SpatRaster") - names(Var_terra) <- names(Data) - - terra::writeCDF(x = Var_terra, filename = file.path(Dir, paste0("SE_", paste0(FileName, ".nc"))), overwrite = TRUE) - # writeRaster(x = Ras_Var, filename = file.path(Dir, paste0("SE_",FileName)), overwrite = TRUE, format="CDF") # save final product as raster - }else{ # if some layers needed to be skipped - warning(paste0("Some of the layers in your raster could not be kriged. You will find all the individual layers (kriged and not kriged) in ", Dir, ".")) - Keep_Temporary <- TRUE # keep temporary files so kriged products are not deleted - } # end of Skip check - - ### REMOVE FILES FROM HARD DRIVE --- - if(Keep_Temporary == FALSE){ # cleanup check - unlink(Dir.Temp, recursive = TRUE) - } # end of cleanup check - - Krig_ls <- list(Ras_Krig, Ras_Var, Call_ls) - names(Krig_ls) <- c("Kriging_Output", "Kriging_SE", "Call") - return(Krig_ls) # return raster or list of layers -} diff --git a/R/misc.R b/R/misc.R deleted file mode 100644 index 9d31452..0000000 --- a/R/misc.R +++ /dev/null @@ -1,92 +0,0 @@ -#' Summary of Raster file characteristics -#' -#' This function is called upon in the krigR function and summarizes Raster characteristics without carrying along the raster file itself. This is used to create lists tracking calls to the function krigR without bloating them too much. -#' -#' @param Object_ras A raster object. -#' -#' @return A list containing information about the input raster. -#' -#'@export -SummarizeRaster <- function(Object_ras = NULL){ - Summary_ls <- list(Class = class(Object_ras), - Dimensions = list(nrow = nrow(Object_ras), - ncol = ncol(Object_ras), - ncell = ncell(Object_ras)), - Extent = Object_ras@extent, - CRS = crs(Object_ras), - layers = names(Object_ras)) - return(Summary_ls) -} - -#' Square Buffers Around Point Data -#' -#' Allow for drawing of buffer zones around point-location data for downloading and kriging of spatial data around point-locations. Overlapping individual buffers are merged. -#' -#' @param Points A data.frame containing geo-referenced points with Lat and Lon columns -#' @param Buffer Identifies how big a rectangular buffer to draw around points. Expressed as centessimal degrees. -#' @param ID Identifies which column in to use for creation of individual buffers. -#' -#' @return A shape made up of individual square buffers around point-location input. -#' -#' @export -buffer_Points <- function(Points = NULL, Buffer = .5, ID = "ID"){ - # set the radius for the plots - radius <- Buffer # radius in meters - # define the plot edges based upon the plot radius. - yPlus <- Points$Lat+radius - xPlus <- Points$Lon+radius - yMinus <- Points$Lat-radius - xMinus <- Points$Lon-radius - # calculate polygon coordinates for each plot centroid. - square=cbind(xMinus,yPlus, # NW corner - xPlus, yPlus, # NE corner - xPlus,yMinus, # SE corner - xMinus,yMinus, # SW corner - xMinus,yPlus) # NW corner again - close ploygon - # Extract the plot ID information - if(any(is.na(Points[,ID]))){stop("Your ID column in your point-location data frame contains NA values. We cannot process these.")} - ID = Points[,ID] - # create spatial polygons from coordinates - polys <- SpatialPolygons(mapply(function(poly, id) - { - xy <- matrix(poly, ncol=2, byrow=TRUE) - Polygons(list(Polygon(xy)), ID=id) - }, - split(square, row(square)), ID), - proj4string = CRS(as.character("+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")) - ) -} - - - -#' Range Masking with Edge Support -#' -#' Creating a raster mask identifying all cells in the original raster (`base.map`) which are at least partially covered by the supplied shapefile (`Shape`). -#' -#' @param base.map A raster within which coverage should be identified -#' @param Shape A polygon(-collection) whose coverage of the raster object is to be found. -#' -#' @return A raster layer. -#' -#' @export -mask_Shape <- function(base.map = NULL, Shape = NULL){ - base.map[] <- NA - stars.base.map <- stars::st_as_stars(base.map) - # Subset shape file - select.ranges <- sf::st_as_sf(Shape) - # Cast polygon as lines instead - select.ranges.lines <- sf::st_cast(select.ranges, "MULTILINESTRING") - select.ranges.lines$STARS <- 1 - # Get centroids (FAST!) - range <- fasterize(select.ranges, base.map, fun = "first", background = 0) - # Get edges (slower than fasterize but faster than rasterize) - range.edges <- stars::st_rasterize(select.ranges.lines, stars.base.map, options = "ALL_TOUCHED=TRUE") - if(class(as.vector(range.edges[[1]])) == "stars"| class(as.vector(range.edges[[1]])) == "numeric"){ - range.edges <- as.vector(range.edges[[1]]) - range.edges <- ifelse(is.na(range.edges), 0, 1) - # Merge - range[] <- ifelse(range[] + range.edges, 1, 0) - } - range[range==0] <- NA # set all cells which the shape doesn't touch to NA - return(range) -} diff --git a/R/onLoad.R b/R/onLoad.R new file mode 100644 index 0000000..f8f5028 --- /dev/null +++ b/R/onLoad.R @@ -0,0 +1,14 @@ +# .onLoad <- function(libname, pkgname) +# { +# library.dynam("mclust", pkgname, libname) +# } +.onAttach <- function(lib, pkg) { + # startup message + msg <- paste0("This is KrigR (version ", packageVersion("KrigR"), "). Note that as of version 0.5.0, considerable changes from the development path have been merged with the main branch leading to the deprecation of some older and outdated functionality. If you have used a version of KrigR prior to 0.5.0 before, we strongly recommend you re-familiarise yourself with the complete suite of KrigR. This message will keep showing until KrigR version 1.0.0 is achieved.") + packageStartupMessage(msg) + # invisible() +} + + +usethis::use_citation() +# Please cite it as Kusch, E., & Davy, R. (2022). KrigR-a tool for downloading and statistically downscaling climate reanalysis data. Environmental Research Letters, 17(2). https://doi.org/10.1088/1748-9326/ac48b3" diff --git a/README.md b/README.md index 77990a1..0b8ba27 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,15 @@ # KrigR An R Package for downloading, preprocessing, and statistical downscaling of the European Centre for Medium-range Weather Forecasts ReAnalysis 5 (ERA5) family provided by the [European Centre for Medium‐Range Weather Forecasts (ECMWF)](https://www.ecmwf.int/).The package interfaces with the [Climate Data Store](https://cds.climate.copernicus.eu/#!/home) hosted by the [Copernicus Climate Change Service (C3S)](https://cds.climate.copernicus.eu/about-c3s) for retrieval of climate data. - + + + KrigR contains functions for: - - Downloading Era5(-Land) data directly from within `R` via a wrapper function for the [`ecmwfr` package](https://github.com/bluegreen-labs/ecmwfr) - - Downloading [USGS GMTED 2010](https://www.usgs.gov/core-science-systems/eros/coastal-changes-and-impacts/gmted2010?qt-science_support_page_related_con=0#qt-science_support_page_related_con) elevation data - - Kriging spatial input to desired output using user-specified covariates via a wrapper for the [`automap` package](https://github.com/cran/automap) - - Downloading and Kriging Era5(Land) data using USGS GMTED 2010 elevation as covariate data in one function call +- Downloading Era5(-Land) data directly from within `R` via a wrapper function for the [`ecmwfr` package](https://github.com/bluegreen-labs/ecmwfr) +- Downloading interpolation covaraite data: + - [USGS GMTED 2010](https://www.usgs.gov/core-science-systems/eros/coastal-changes-and-impacts/gmted2010?qt-science_support_page_related_con=0#qt-science_support_page_related_con) elevation data + - [soil hydraulic and thermal parameters](http://globalchange.bnu.edu.cn/research/soil4.jsp) data +- Preparing covariate data for use in statistical interpolation +- Kriging spatial input to desired output using user-specified covariates via a wrapper for the [`automap` package](https://github.com/cran/automap) **NOTE:** All kriging functionality can be parallelised for faster computing using the `Cores` argument. @@ -13,24 +17,16 @@ KrigR contains functions for: KrigR has been [published here](https://iopscience.iop.org/article/10.1088/1748-9326/ac48b3). # Installation -A KrigR dependency, rgdal, is no longer available at CRAN. You need to download and install it separately as follows: -```{r} -remotes::install_github("https://github.com/cran/rgdal") -``` - -I will work to remove this dependency as soon as possible and hope to achieve this (alongside additional KrigR development) in 2024. - KrigR is not yet on CRAN, so it needs to be installed as such: ```{r} -Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS="true") devtools::install_github("https://github.com/ErikKusch/KrigR") library(KrigR) ``` -Users require personal API-access tokens which can be obtained [here](https://cds.climate.copernicus.eu/api-how-to). +Users require personal API-access tokens which can be obtained [here](https://accounts.ecmwf.int/auth/realms/ecmwf/login-actions/registration?client_id=cds&tab_id=VkbipqjwuIQ). # Workshop Material and Further Information -We have put together a comprehensive [workshop](https://www.erikkusch.com/courses/krigr/) that walks you through the functionality of the KrigR package. This workshop was presented during the BIOCHANGE Methods Workshop Series at Aarhus University. A recording of this presentation can be found on [YouTube](https://www.youtube.com/watch?v=wwb107L4wVw&ab_channel=ErikKusch). For any additional information on the project please refer to the project [website](https://www.erikkusch.com/project/krigr/). Upates will be available on said website as well as the KrigR [twitter profile](https://twitter.com/ERAKrigR) +I have put together a comprehensive [workshop](https://www.erikkusch.com/courses/krigr/) that walks you through the functionality of the KrigR package. For any additional information on the project please refer to the project [website](https://www.erikkusch.com/project/krigr/). Upates will be available on said website as well as the KrigR [twitter profile](https://twitter.com/ERAKrigR) # Abstract Advances in climate science have rendered obsolete the gridded observation data widely used in downstream applications. Novel climate reanalysis products outperform legacy data products in accuracy, temporal resolution, and provision of uncertainty metrics. Consequently, there is an urgent need to develop a workflow through which to integrate these improved data into biological analyses. The ERA5 product family (ERA5 and ERA5-Land) are the latest and most advanced global reanalysis products created by the European Center for Medium-range Weather Forecasting (ECMWF). These data products offer up to 83 essential climate variables (ECVs) at hourly intervals for the time-period of 1981 to today with preliminary back-extensions being available for 1950-1981. Spatial resolutions range from 30x30km (ERA5) to 11x11km (ERA5-Land) and can be statistically downscaled to study-requirements at finer spatial resolutions. Kriging is one such method to interpolate data to finer resolutions and has the advantages that one can leverage additional covariate information and obtain the uncertainty associated with the downscaling. diff --git a/data-raw/DATASET.R b/data-raw/DATASET.R index dc8d7be..43e31c7 100644 --- a/data-raw/DATASET.R +++ b/data-raw/DATASET.R @@ -24,10 +24,36 @@ CDS_rast <- CDownloadS( Dir = file.path(getwd(), "inst", "extdata"), FileName = "CentralNorway", API_User = API_User, - API_Key = API_Key) + API_Key = API_Key +) usethis::use_data(CDS_rast) ## Jotunheimen_poly Jotunheimen_poly <- sf::st_read("data-raw/Shape/Shape-polygon.shp") -sf::st_transform(Jotunheimen_poly, crs =sf::st_crs(CDS_rast)) +sf::st_transform(Jotunheimen_poly, crs = sf::st_crs(CDS_rast)) usethis::use_data(Jotunheimen_poly) + +## Bioclimatic variables +CN_BC <- BioClim( + Temperature_Var = "2m_temperature", + Temperature_DataSet = "reanalysis-era5-land", + Temperature_Type = NA, + Water_Var = "volumetric_soil_water_layer_1", + Water_DataSet = "reanalysis-era5-land-monthly-means", + Water_Type = "monthly_averaged_reanalysis", + Y_start = 1970, + Y_end = 1979, + TZone = "CET", + Extent = CN_ext, + Dir = getwd(), + FileName = "CN_BC", + FileExtension = ".nc", + Compression = 9, + API_User = API_User, + API_Key = API_Key, + TChunkSize = 6000, TryDown = 10, TimeOut = 36000, + Cores = parallel::detectCores(), + verbose = TRUE, + Keep_Raw = FALSE, Keep_Monthly = FALSE +) +usethis::use_data(CN_BC) diff --git a/inst/CITATION b/inst/CITATION new file mode 100644 index 0000000..83850d3 --- /dev/null +++ b/inst/CITATION @@ -0,0 +1,11 @@ +bibentry( + bibtype = "Article", + title = "KrigR-a tool for downloading and statistically downscaling climate reanalysis data", + author = "Erik Kusch and Richard Davy", + journal = "Environmental Research Letters", + year = "2022", + volume = "17", + number = "2", + pages = "", + doi = "10.1088/1748-9326/ac48b3" +) \ No newline at end of file diff --git a/inst/extdata/CN_BC.nc b/inst/extdata/CN_BC.nc new file mode 100644 index 0000000..688078f Binary files /dev/null and b/inst/extdata/CN_BC.nc differ diff --git a/inst/extdata/KrigRLogo.png b/inst/extdata/KrigRLogo.png new file mode 100644 index 0000000..9f2483f Binary files /dev/null and b/inst/extdata/KrigRLogo.png differ diff --git a/man/BioClim.Rd b/man/BioClim.Rd index 34d809a..e80552e 100644 --- a/man/BioClim.Rd +++ b/man/BioClim.Rd @@ -5,87 +5,118 @@ \title{Computation of Bioclimatic Variables} \usage{ BioClim( + Temperature_Var = "2m_temperature", + Temperature_DataSet = "reanalysis-era5-land", + Temperature_Type = NA, Water_Var = "volumetric_soil_water_layer_1", - Y_start = 1981, - Y_end = 2015, - T_res = "day", - Extent = extent(-180, 180, -90, 90), - DataSet = "era5-land", - Dir = getwd(), - verbose = TRUE, - FileName = "NULL", - Keep_Raw = FALSE, - Keep_Monthly = FALSE, + Water_DataSet = "reanalysis-era5-land-monthly-means", + Water_Type = "monthly_averaged_reanalysis", + Y_start, + Y_end, + TZone = "UTC", + Extent, Buffer = 0.5, - ID = "ID", - API_User = API_User, - API_Key = API_Key, - Cores = 1, + Dir = getwd(), + FileName, + FileExtension = ".nc", + Compression = 9, + API_User, + API_Key, + TChunkSize = 6000, TryDown = 10, TimeOut = 36000, - SingularDL = FALSE + Cores = 1, + verbose = TRUE, + Keep_Raw = FALSE, + Keep_Monthly = FALSE ) } \arguments{ -\item{Water_Var}{ERA5(Land)-contained climate variable targeting water availability information. Recommended values: "volumetric_soil_water_layer_1", "total_precipitation".} +\item{Temperature_Var}{CDS-contained climate variable targeting temperature information. Recommended values: "2m_temperature".} -\item{Y_start}{Year ('YYYY') at which to start time series of downloaded data.} +\item{Temperature_DataSet}{Character. Which dataset to query data from. See currently supported datasets by calling \code{\link{Meta.List}}. For now, this function is conceptualised to support "reanalysis-era5-land".} -\item{Y_end}{Year ('YYYY') at which to stop time series of downloaded data.} +\item{Temperature_Type}{Either NA or Character. Which kind of sub-type to query per data set. See \code{\link{Meta.QucikFacts}} for options per dataset.} -\item{T_res}{Temporal resolution from which to obtain minimum and maximum values of temperature. 'hour' or 'day'} +\item{Water_Var}{CDS-contained climate variable targeting water availability information. Recommended values: "volumetric_soil_water_layer_X" (where X is an integer of either 1, 2, 3, 4), "total_precipitation".} -\item{Extent}{Optional, download data according to rectangular bounding box. specify as extent() object or as a raster, a SpatialPolygonsDataFrame object, or a data.frame opbject. If Extent is a SpatialPolygonsDataFrame, this will be treated as a shapefile and the output will be cropped and masked to this shapefile. If Extent is a data.frame of geo-referenced point records, it needs to contain Lat and Lon columns as well as a non-repeating ID-column.} +\item{Water_DataSet}{Character. Which dataset to query water availability data from. See currently supported datasets by calling \code{\link{Meta.List}}. For now, this function is conceptualised to support "reanalysis-era5-land".} -\item{DataSet}{Which ERA5 data set to download data from. 'era5' or 'era5-land'.} +\item{Water_Type}{Either NA or Character. Which kind of sub-type to query per water availability data set. See \code{\link{Meta.QucikFacts}} for options per dataset.} -\item{Dir}{Directory specifying where to download data to.} +\item{Y_start}{Year ('YYYY') at which to start time series of downloaded data.} -\item{verbose}{Optional, logical. Whether to report progress of the function in the console or not.} +\item{Y_end}{Year ('YYYY') at which to stop time series of downloaded data.} -\item{FileName}{A file name for the netcdf produced.} +\item{TZone}{Character. Time zone in which to represent and evaluate time dimension of data. See the output of OlsonNames() for a full overview of supported specifications. Default is UTC.} -\item{Keep_Raw}{Logical. Whether to keep monthly netcdf files of raw data aggregated to temporal resolution of `T_res`. Default FALSE.} +\item{Extent}{Optional, download data according to desired spatial specification. If missing/unspecified, total area of queried data set is used. Can be specified either as a raster object, an sf object, a terra object, or a data.frame. If Extent is a raster or terra object, data will be queried according to rectangular extent thereof. If Extent is an sf (MULTI-)POLYGON object, this will be treated as a shapefile and the output will be cropped and masked to this shapefile. If Extent is a data.frame of geo-referenced point records, it needs to contain Lat and Lon columns around which a buffered shapefile will be created using the Buffer argument.} -\item{Keep_Monthly}{Logical. Whether to keep monthly netcdf files of raw data aggregated to temporal resolution of months. Default FALSE.} +\item{Buffer}{Optional, Numeric. Identifies how big a circular buffer to draw around points if Extent is a data.frame of points. Buffer is expressed as centessimal degrees.} -\item{Buffer}{Optional. Identifies how big a rectangular buffer to draw around points if Extent is a data frame of points. Buffer is expressed as centessimal degrees.} +\item{Dir}{Character/Directory Pointer. Directory specifying where to download data to.} -\item{ID}{Optional. Identifies which column in Extent to use for creation of individual buffers if Extent is a data.frame.} +\item{FileName}{Character. A file name for the produced file.} + +\item{FileExtension}{Character. A file extension for the produced file. Supported values are ".nc" (default) and ".tif" (better support for metadata).} + +\item{Compression}{Integer between 1 to 9. Applied to final .nc file that the function writes to hard drive. Same as compression argument in terra::writeCDF(). Ignored if FileExtension = ".tif".} \item{API_User}{Character; ECMWF cds user number.} \item{API_Key}{Character; ECMWF cds API key.} -\item{Cores}{Numeric. How many cores to use.^This can speed up downloads of long time-series. If you want output to your console during the process, use Cores = 1. Parallel processing is carried out when Cores is bigger than 1. Default is 1.} +\item{TChunkSize}{Numeric. Number of layers to bundle in each individual download. Default is 6000 to adhere to most restrictive CDS limits: https://cds.climate.copernicus.eu/live/limits.} + +\item{TryDown}{Optional, numeric. Legacy, ignored when querying data from new CDS (https://cds-beta.climate.copernicus.eu/; this happens when the package version of ecmwfr is >= 2.0.0). How often to attempt the download of each individual file that the function queries from the CDS. This is to circumvent having to restart the entire function when encountering connectivity issues.} + +\item{TimeOut}{Numeric. Legacy, ignored when querying data from new CDS (https://cds-beta.climate.copernicus.eu/; this happens when the package version of ecmwfr is >= 2.0.0). The timeout for each download in seconds. Default 36000 seconds (10 hours).} -\item{TryDown}{Optional, numeric. How often to attempt the download of each individual file that the function queries from the server. This is to circumvent having to restart the entire function when encountering connectivity issues.} +\item{Cores}{Numeric. How many cores to use when carrying out temporal aggregation. Default is 1.} -\item{TimeOut}{Numeric. The timeout for each download in seconds. Default 36000 seconds (10 hours).} +\item{verbose}{Logical. Whether to print/message function progress in console or not.} -\item{SingularDL}{Logical. Whether to force download of data in one call to CDS or automatically break download requests into individual monthly downloads. Default is FALSE.} +\item{Keep_Raw}{Logical. Whether to retain raw downloaded data or not. Default is FALSE.} + +\item{Keep_Monthly}{Logical. Whether to keep monthly netcdf files of raw data aggregated to temporal resolution of months. Default FALSE.} } \value{ -A raster object containing the downloaded ERA5(-Land) data, and a NETCDF (.nc) file in the specified directory. +A SpatRaster object containing the queried bioclimatic data, and a NETCDF (.nc) file in the specified directory. + +The SpatRaster contains metadata/attributes as a named vector that can be retrieved with terra::metags(...): +\itemize{ +\item{Citation}{ - A string which to use for in-line citation of the data product obtained with BioClim}. +\item{KrigRCall.X}{ - Arguments passed to the BioClim function that produced the file (API credentials are omitted from these metadata)}. +} } \description{ This function queries download of required essential climate variables from the [Climate Data Store](https://cds.climate.copernicus.eu/#!/home) hosted by the [Copernicus Climate Change Service (C3S)](https://cds.climate.copernicus.eu/about-c3s) for retrieval of climate data and subsequent calculation of bioclimatic variables for user-defined regions and time-frames. } \examples{ \dontrun{ -BioClim (Water_Var = "volumetric_soil_water_layer_1", # could also be total_precipitation -Y_start = 1981, -Y_end = 2015, -T_res = "day", -Extent = extent(11.8,15.1,50.1,51.7), -DataSet = "era5-land", -Dir = getwd(), -verbose = TRUE, -FileName = "NULL", -Keep_Raw = FALSE, -Keep_Monthly = FALSE, -API_User = API_User, -API_Key = API_Key) +CN_ext <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +CN_BC <- BioClim( + Temperature_Var = "2m_temperature", + Temperature_DataSet = "reanalysis-era5-land", + Temperature_Type = NA, + Water_Var = "volumetric_soil_water_layer_1", + Water_DataSet = "reanalysis-era5-land-monthly-means", + Water_Type = "monthly_averaged_reanalysis", + Y_start = 1970, + Y_end = 1979, + TZone = "CET", + Extent = CN_ext, + Dir = getwd(), + FileName = "CN_BC", + FileExtension = ".nc", + Compression = 9, + API_User = API_User, + API_Key = API_Key, + TChunkSize = 6000, TryDown = 10, TimeOut = 36000, + Cores = parallel::detectCores(), + verbose = TRUE, + Keep_Raw = FALSE, Keep_Monthly = FALSE +) } } diff --git a/man/CDownloadS.Rd b/man/CDownloadS.Rd index 4398173..5fead06 100644 --- a/man/CDownloadS.Rd +++ b/man/CDownloadS.Rd @@ -20,6 +20,7 @@ CDownloadS( Dir = getwd(), FileName, FileExtension = ".nc", + Compression = 9, API_User, API_Key, TryDown = 10, @@ -27,8 +28,7 @@ CDownloadS( TChunkSize = 6000, Cores = 1, verbose = TRUE, - Keep_Raw = FALSE, - Save_Final = TRUE + Keep_Raw = FALSE ) } \arguments{ @@ -62,13 +62,15 @@ CDownloadS( \item{FileExtension}{Character. A file extension for the produced file. Supported values are ".nc" (default) and ".tif" (better support for metadata).} +\item{Compression}{Integer between 1 to 9. Applied to final .nc file that the function writes to hard drive. Same as compression argument in terra::writeCDF(). Ignored if FileExtension = ".tif".} + \item{API_User}{Character; ECMWF cds user number.} \item{API_Key}{Character; ECMWF cds API key.} -\item{TryDown}{Optional, numeric. How often to attempt the download of each individual file that the function queries from the CDS. This is to circumvent having to restart the entire function when encountering connectivity issues.} +\item{TryDown}{Optional, numeric. Legacy, ignored when querying data from new CDS (https://cds-beta.climate.copernicus.eu/; this happens when the package version of ecmwfr is >= 2.0.0). How often to attempt the download of each individual file that the function queries from the CDS. This is to circumvent having to restart the entire function when encountering connectivity issues.} -\item{TimeOut}{Numeric. The timeout for each download in seconds. Default 36000 seconds (10 hours).} +\item{TimeOut}{Numeric. Legacy, ignored when querying data from new CDS (https://cds-beta.climate.copernicus.eu/; this happens when the package version of ecmwfr is >= 2.0.0). The timeout for each download in seconds. Default 36000 seconds (10 hours).} \item{TChunkSize}{Numeric. Number of layers to bundle in each individual download. Default is 6000 to adhere to most restrictive CDS limits: https://cds.climate.copernicus.eu/live/limits.} @@ -77,8 +79,6 @@ CDownloadS( \item{verbose}{Logical. Whether to print/message function progress in console or not.} \item{Keep_Raw}{Logical. Whether to retain raw downloaded data or not. Default is FALSE.} - -\item{Save_Final}{Logical. Whether to write the final SpatRaster to the hard drive. Default is TRUE.} } \value{ A SpatRaster object containing the downloaded, cropped/masked, and subsequently temporally aggregated data, and a file (either .nc or .tif) in the specified directory. @@ -98,46 +98,46 @@ This function is used to obtain data from the \href{https://cds.climate.copernic \dontrun{ ## Raw data for one month of full globe RawGlobe_rast <- CDownloadS( - Variable = "2m_temperature", - DataSet = "reanalysis-era5-land-monthly-means", - Type = "monthly_averaged_reanalysis", - # time-window, default set to range of dataset-type - DateStart = "1995-01-01 00:00", - DateStop = "1995-01-01 23:00", - TZone = "CET", - # temporal aggregation - TResolution = "month", - TStep = 1, - # file storing - FileName = "RawGlobe", - # API credentials - API_User = API_User, - API_Key = API_Key + Variable = "2m_temperature", + DataSet = "reanalysis-era5-land-monthly-means", + Type = "monthly_averaged_reanalysis", + # time-window, default set to range of dataset-type + DateStart = "1995-01-01 00:00", + DateStop = "1995-01-01 23:00", + TZone = "CET", + # temporal aggregation + TResolution = "month", + TStep = 1, + # file storing + FileName = "RawGlobe", + # API credentials + API_User = API_User, + API_Key = API_Key ) -terra::plot(RawGlobe_rast) +Plot.SpatRast(RawGlobe_rast) ## Monthly air temperature aggregated to bi-annual maximum by SpatRaster -CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package="KrigR")) +CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) BiAnnAirTemp_rast <- CDownloadS( - Variable = "2m_temperature", - DataSet = "reanalysis-era5-land-monthly-means", - Type = "monthly_averaged_reanalysis", - # time-window, default set to range of dataset-type - DateStart = "1995-01-01 00:00", - DateStop = "1996-12-31 23:00", - TZone = "EET", - # temporal aggregation - TResolution = "year", - TStep = 2, - # spatial - Extent = CDS_rast, - # file storing - FileName = "BiAnnAirTemp", - # API credentials - API_User = API_User, - API_Key = API_Key + Variable = "2m_temperature", + DataSet = "reanalysis-era5-land-monthly-means", + Type = "monthly_averaged_reanalysis", + # time-window, default set to range of dataset-type + DateStart = "1995-01-01 00:00", + DateStop = "1996-12-31 23:00", + TZone = "EET", + # temporal aggregation + TResolution = "year", + TStep = 2, + # spatial + Extent = CDS_rast, + # file storing + FileName = "BiAnnAirTemp", + # API credentials + API_User = API_User, + API_Key = API_Key ) -terra::plot(BiAnnAirTemp_rast) +Plot.SpatRast(BiAnnAirTemp_rast) ## Hourly back-calculated precipitation aggregated to daily averages by shapefiles data("Jotunheimen_poly") @@ -159,36 +159,36 @@ DailyBackCPrecip_rast <- CDownloadS( API_User = API_User, API_Key = API_Key ) -terra::plot(DailyBackCPrecip_rast) +Plot.SpatRast(DailyBackCPrecip_rast, SF = Jotunheimen_poly, Legend = "Precipitation [m]", COL = rev(viridis::cividis(100))) ## 6-hourly ensemble member spread sum for air temperature by buffered points data("Mountains_df") EnsembleSpreadSum6hour_rast <- CDownloadS( - Variable = "2m_temperature", - DataSet = "reanalysis-era5-single-levels", - Type = "ensemble_spread", - # time-window, default set to range of dataset-type - DateStart = "1995-01-01 00:00:00", - DateStop = "1995-01-01 21:00:00", - TZone = "UTC", - # temporal aggregation - TResolution = "hour", - TStep = 6, - FUN = sum, - # spatial - Extent = Mountains_df, - Buffer = 0.2, - # file storing - FileName = "EnsembleSpreadSum6hour", - FileExtension = ".tif", - # API credentials - API_User = API_User, - API_Key = API_Key, - Keep_Raw = TRUE + Variable = "2m_temperature", + DataSet = "reanalysis-era5-single-levels", + Type = "ensemble_spread", + # time-window, default set to range of dataset-type + DateStart = "1995-01-01 00:00:00", + DateStop = "1995-01-01 21:00:00", + TZone = "UTC", + # temporal aggregation + TResolution = "hour", + TStep = 6, + FUN = sum, + # spatial + Extent = Mountains_df, + Buffer = 0.2, + # file storing + FileName = "EnsembleSpreadSum6hour", + FileExtension = ".tif", + # API credentials + API_User = API_User, + API_Key = API_Key, + Keep_Raw = TRUE ) -terra::plot(EnsembleSpreadSum6hour_rast) +Plot.SpatRast(EnsembleSpreadSum6hour_rast, Legend = "Air Temperature Uncertainty [K]") } } \seealso{ -\code{\link{Meta.List}}, \code{\link{Meta.Variables}}, \code{\link{Meta.QuickFacts}}. +\code{\link{Meta.List}}, \code{\link{Meta.Variables}}, \code{\link{Meta.QuickFacts}}, \code{\link{Plot.SpatRast}}. } diff --git a/man/Check.File.Rd b/man/Check.File.Rd index 131d277..f48fdc5 100644 --- a/man/Check.File.Rd +++ b/man/Check.File.Rd @@ -25,8 +25,8 @@ If a file already exists in a given place, load that file } \examples{ KrigR::Check.File( - FName = basename(system.file("extdata", "CentralNorway.nc", package="KrigR")), - Dir = dirname(system.file("extdata", "CentralNorway.nc", package="KrigR")), - loadFun = terra::rast - ) + FName = basename(system.file("extdata", "CentralNorway.nc", package = "KrigR")), + Dir = dirname(system.file("extdata", "CentralNorway.nc", package = "KrigR")), + loadFun = terra::rast +) } diff --git a/man/CovariateSetup.Rd b/man/CovariateSetup.Rd index d35710b..4f831c0 100644 --- a/man/CovariateSetup.Rd +++ b/man/CovariateSetup.Rd @@ -7,13 +7,15 @@ CovariateSetup( Training, Target, + FilePrefix = "", Covariates = "GMTED2010", Source = "Origin", Extent, Buffer = 0.5, Dir = getwd(), Keep_Global = FALSE, - FileExtension = ".nc" + FileExtension = ".nc", + Compression = 9 ) } \arguments{ @@ -21,7 +23,9 @@ CovariateSetup( \item{Target}{Either numeric or a SpatRaster. If numeric, a single number representing the target resolution for the kriging step (i.e. wich resolution to downscale to). If a SpatRaster, data that the covariates and kriged products should align with. In case of a numeric input, covariate data is aggregated as closely as possible to desired resolution. If a SpatRaster, covariate data is resampled to match desired output directly.} -\item{Covariates}{Either character or a SpatRaster. If character, obtain frequently used and provably useful covariate data (i.e., GMTED2010 and HWSD) and prepare for use in Kriging. Supported character values are "GMTED2010" and "HWSD". Note that currently, HWSD data download is not functional. If a SpatRaster, a user-supplied set of covariate data to be prepared for use in Kriging.} +\item{FilePrefix}{Character. A file name prefix for the produced files.} + +\item{Covariates}{Either character or a SpatRaster. If character, obtain frequently used and provably useful covariate data (i.e., GMTED2010 and soil data) and prepare for use in Kriging. Supported character values are "GMTED2010", "tksat", "tkdry", "csol", "k_s", "lambda", "psi", and "theta_s". If a SpatRaster, a user-supplied set of covariate data to be prepared for use in Kriging.} \item{Source}{Character. Only comes into effect when Covariates argument is specified as a character. Whether to attempt download of covariate data from the official sources (Source = "Origin") or a static copy of the data set on a private drive (Source = "Drive"). Default is "Origin".} @@ -34,6 +38,8 @@ CovariateSetup( \item{Keep_Global}{Logical. Only comes into effect when Covariates argument is specified as a character. Whether to retain raw downloaded covariate data or not. Default is FALSE.} \item{FileExtension}{Character. A file extension for the produced files. Supported values are ".nc" (default) and ".tif" (better support for metadata).} + +\item{Compression}{Integer between 1 to 9. Applied to final .nc file that the function writes to hard drive. Same as compression argument in terra::writeCDF(). Ignored if FileExtension = ".tif".} } \value{ A list containing two SpatRaster objects (Training and Target) ready to be used as covariates for kriging, and two files called Covariates_Target and Covariates_Train in the specified directory. @@ -44,46 +50,50 @@ The SpatRasters produced and stored when specifying the Covariates argument as a } } \description{ -This function is used to setup products of covariate data ready for use in Kriging. This functiuonality can either be applied to user-supplied covariate data or ready-made data products such as the Harmised World Soil Data Base and the median statistic of the Global Multi-resolution Terrain Elevation Data (GMTED2010; available at \url{https://topotools.cr.usgs.gov/gmted_viewer/gmted2010_global_grids.php}). In case of the latter, the data is downloaded at 30 arc-sec latitude/longitude grid cells and subsequently resampled to match training and target resolutions specified by the user. +This function is used to setup products of covariate data ready for use in Kriging. This functionality can either be applied to user-supplied covariate data or ready-made data products such as the global dataset of soil hydraulic and thermal parameters for earth system modeling (available at \url{http://globalchange.bnu.edu.cn/research/soil4.jsp}) and the median statistic of the Global Multi-resolution Terrain Elevation Data (GMTED2010; available at \url{https://topotools.cr.usgs.gov/gmted_viewer/gmted2010_global_grids.php}). In case of the latter, the data is downloaded at 30 arc-sec latitude/longitude grid cells and subsequently resampled to match training and target resolutions specified by the user. } \examples{ \dontrun{ ## Rectangular Covariate data according to input data -CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package="KrigR")) -Covariates_ls <- CovariateSetup(Training = CDS_rast, - Target = 0.01, - Covariates = "GMTED2010", - Keep_Global = TRUE, - FileExtension = ".nc") -terra::plot(Covariates_ls[[1]]) -terra::plot(Covariates_ls[[2]]) +CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +Covariates_ls <- CovariateSetup( + Training = CDS_rast, + Target = 0.01, + Covariates = "GMTED2010", + Keep_Global = TRUE, + FileExtension = ".nc" +) +Plot.Covariates(Covariates_ls) ## Shapefile-limited covariate data data("Jotunheimen_poly") -CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package="KrigR")) -Covariates_ls <- CovariateSetup(Training = CDS_rast, - Target = 0.01, - Covariates = "GMTED2010", - Extent = Jotunheimen_poly, - Keep_Global = TRUE, - FileExtension = ".nc") -terra::plot(Covariates_ls[[1]]) -terra::plot(Covariates_ls[[2]]) +CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +Covariates_ls <- CovariateSetup( + Training = CDS_rast, + Target = 0.01, + Covariates = "GMTED2010", + Extent = Jotunheimen_poly, + Keep_Global = TRUE, + FileExtension = ".nc" +) +Plot.Covariates(Covariates_ls, SF = Jotunheimen_poly) ## buffered-point-limited covariate data data("Mountains_df") -CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package="KrigR")) -Covariates_ls <- CovariateSetup(Training = CDS_rast, - Target = 0.01, - Covariates = "GMTED2010", - Extent = Mountains_df, - Buffer = 0.2, - Keep_Global = TRUE, - FileExtension = ".nc") -terra::plot(Covariates_ls[[1]]) -terra::plot(Covariates_ls[[2]]) +CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +Covariates_ls <- CovariateSetup( + Training = CDS_rast, + Target = 0.01, + Covariates = c("tksat", "tkdry", "csol", "k_s", "lambda", "psi", "theta_s"), + Source = "Drive", + Extent = Mountains_df, + Buffer = 0.2, + Keep_Global = TRUE, + FileExtension = ".nc" +) +Plot.Covariates(Covariates_ls) } } \seealso{ -\code{\link{Kriging}}. +\code{\link{Kriging}}, \code{\link{Plot.Covariates}}. } diff --git a/man/Ext.Check.Rd b/man/Ext.Check.Rd index a4b9e9a..69cc85a 100644 --- a/man/Ext.Check.Rd +++ b/man/Ext.Check.Rd @@ -16,7 +16,7 @@ A list containg (1) a terra/sf object and (2) the corresponding SpatExtent objec Try to convert user input into (1) a terra or sf object and also read out the corresponding (2) SpatExtent object. Supports inputs of classes belonging to the packages raster, terra, sf, and sp } \examples{ - ## raster +## raster Check.Ext(raster::extent(c(9.87, 15.03, 49.89, 53.06))) ## terra Check.Ext(terra::ext(c(9.87, 15.03, 49.89, 53.06))) @@ -24,7 +24,7 @@ Check.Ext(terra::ext(c(9.87, 15.03, 49.89, 53.06))) set.seed(42) nb_pt <- 10 dd <- data.frame(x = runif(nb_pt, 9.87, 15.03), y = runif(nb_pt, 49.89, 53.06), val = rnorm(nb_pt)) -sf <- sf::st_as_sf(dd, coords = c("x","y")) +sf <- sf::st_as_sf(dd, coords = c("x", "y")) Check.Ext(sf) ## sp Check.Ext(as(sf, "Spatial")) diff --git a/man/Kriging.Rd b/man/Kriging.Rd index ae2076d..d65bd7c 100644 --- a/man/Kriging.Rd +++ b/man/Kriging.Rd @@ -9,11 +9,12 @@ Kriging( Covariates_training, Covariates_target, Equation = NULL, - Cores = detectCores(), + Cores = parallel::detectCores(), nmax = Inf, Dir = getwd(), FileName, - FileExtension, + FileExtension = ".nc", + Compression = 9, Keep_Temporary = FALSE, verbose = TRUE ) @@ -37,6 +38,8 @@ Kriging( \item{FileExtension}{Character. A file extension for the produced file. Supported values are ".nc" (default) and ".tif" (better support for metadata).} +\item{Compression}{Integer between 1 to 9. Applied to final .nc file that the function writes to hard drive. Same as compression argument in terra::writeCDF(). Ignored if FileExtension = ".tif".} + \item{Keep_Temporary}{Logical, whether to delete individual kriging products of layers in Data after processing. Default is TRUE. These temporary files are stored in a newly created directory in Dir which is pre-pended with "TEMP-" and is deleted if Keep_Temporary = FALSE upon completion.} \item{verbose}{Optional, logical. Whether to report progress of data download (if queried) in the console or not.} @@ -60,10 +63,11 @@ Use optional arguments such as Dir, Keep_Temporary, KrigingEquation, nmax and Co \dontrun{ ## Kriging using pre-fab data with a rectangular extent and a fives layers of data with parallel processing ### Loading data -CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package="KrigR")) -Cov_train <- terra::rast(system.file("extdata", "Covariates_Train.nc", package="KrigR")) -Cov_target <- terra::rast(system.file("extdata", "Covariates_Target.nc", package="KrigR")) -terra::varnames(Cov_train) <- terra::varnames(Cov_target) <- "GMTED2010" # we must ensure that the varnames in the Covariate file match +CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +Cov_train <- terra::rast(system.file("extdata", "Covariates_Train.nc", package = "KrigR")) +Cov_target <- terra::rast(system.file("extdata", "Covariates_Target.nc", package = "KrigR")) +names(Cov_train) <- names(Cov_target) <- "GMTED2010" + ### kriging itself ExtentKrig <- Kriging( Data = CDS_rast, @@ -77,7 +81,7 @@ ExtentKrig <- Kriging( nmax = 40, verbose = TRUE ) - +Plot.Kriged(Krigs = ExtentKrig) ## Kriging using full KrigR pipeline with shapefile data ### Shapefile loading data("Jotunheimen_poly") @@ -100,19 +104,21 @@ Qsoil_rast <- CDownloadS( ) ### Covariate preparations -Covariates_ls <- CovariateSetup(Training = Qsoil_rast, - #' Target = 0.03, - Covariates = "GMTED2010", # this shiuld really be HWSD - Extent = Jotunheimen_poly, - Keep_Global = TRUE) -terra::varnames(Covariates_ls[[1]]) <- terra::varnames(Covariates_ls[[2]]) <- "GMTED2010" # we must ensure that the varnames in the Covariate file match +Covariates_ls <- CovariateSetup( + Training = Qsoil_rast, + Target = 0.03, + Covariates = c("tksat", "tkdry", "csol", "k_s", "lambda", "psi", "theta_s"), + Source = "Drive", + Extent = Jotunheimen_poly, + Keep_Global = TRUE +) ### kriging itself ShapeKrig <- Kriging( Data = Qsoil_rast, Covariates_training = Covariates_ls[[1]], Covariates_target = Covariates_ls[[2]], - Equation = "GMTED2010", + Equation = "tksat + tkdry + csol + k_s + lambda + psi + theta_s", Cores = 1, FileName = "KrigTest2", FileExtension = ".nc", @@ -120,8 +126,9 @@ ShapeKrig <- Kriging( nmax = 40, verbose = TRUE ) +Plot.Kriged(Krigs = ShapeKrig, SF = Jotunheimen_poly) } } \seealso{ -\code{\link{CovariateSetup}}. +\code{\link{CovariateSetup}}, \code{\link{Plot.Kriged}}. } diff --git a/man/Make.Request.Rd b/man/Make.Request.Rd index 630fcd5..59e55d9 100644 --- a/man/Make.Request.Rd +++ b/man/Make.Request.Rd @@ -15,7 +15,8 @@ Make.Request( Dir = getwd(), verbose = TRUE, API_User, - API_Key + API_Key, + TimeOut = 36000 ) } \arguments{ @@ -40,6 +41,8 @@ Make.Request( \item{API_User}{Character. CDS API User} \item{API_Key}{Character. CDS API Key} + +\item{TimeOut}{Numeric. Legacy, ignored when querying data from new CDS (https://cds-beta.climate.copernicus.eu/; this happens when the package version of ecmwfr is >= 2.0.0). The timeout for each download in seconds. Default 36000 seconds (10 hours).} } \value{ List. Each element holding either (1) a list object representing a CDS request or (2) the value NA indicating that a file of this name is already present. diff --git a/man/Meta.Check.Rd b/man/Meta.Check.Rd index f015229..8f6c8bb 100644 --- a/man/Meta.Check.Rd +++ b/man/Meta.Check.Rd @@ -46,7 +46,7 @@ List. Contains: Read and return short overview of data set characteristics, supported types, extent, time frames and required arguments. } \examples{ -Meta.Check(DataSet = "reanalysis-era5-land", Type = NA, VariableCheck = "2m_temperature", CumulativeCheck = FALSE, ExtentCheck = c(53.06, 9.87, 49.89, 15.03), DateCheck = data.frame(IN = c(as.POSIXct("1995-01-01 CET"), as.POSIXct("2005-01-01 23:00:00 CET")), UTC = c(as.POSIXct("1994-12-31 23:00:00 UTC"), as.POSIXct("2005-01-01 22:00:00 UTC"))), AggrCheck = list(1, "hour"), QueryTimes = c('00:00', '03:00', '06:00', '09:00', '12:00', '15:00', '18:00', '21:00')) +Meta.Check(DataSet = "reanalysis-era5-land", Type = NA, VariableCheck = "2m_temperature", CumulativeCheck = FALSE, ExtentCheck = c(53.06, 9.87, 49.89, 15.03), DateCheck = data.frame(IN = c(as.POSIXct("1995-01-01 CET"), as.POSIXct("2005-01-01 23:00:00 CET")), UTC = c(as.POSIXct("1994-12-31 23:00:00 UTC"), as.POSIXct("2005-01-01 22:00:00 UTC"))), AggrCheck = list(1, "hour"), QueryTimes = c("00:00", "03:00", "06:00", "09:00", "12:00", "15:00", "18:00", "21:00")) } \seealso{ diff --git a/man/Meta.NC.Rd b/man/Meta.NC.Rd index b07ea3a..d46b90d 100644 --- a/man/Meta.NC.Rd +++ b/man/Meta.NC.Rd @@ -4,7 +4,7 @@ \alias{Meta.NC} \title{Read or write metadata into NetCDF files} \usage{ -Meta.NC(NC, FName, Attrs, Write = FALSE, Read = FALSE) +Meta.NC(NC, FName, Attrs, Write = FALSE, Read = TRUE, Compression = 9) } \arguments{ \item{NC}{SpatRaster} @@ -16,6 +16,8 @@ Meta.NC(NC, FName, Attrs, Write = FALSE, Read = FALSE) \item{Write}{Logical. Whether to write metadata} \item{Read}{Logical Whether to read metadata} + +\item{Compression}{Integer between 1 to 9. Applied to final .nc file that the function writes to hard drive. Same as compression argument in terra::writeCDF(). Ignored if FileExtension = ".tif".} } \value{ A SpatRaster with metadata diff --git a/man/Plot.BioClim.Rd b/man/Plot.BioClim.Rd new file mode 100644 index 0000000..7bf1709 --- /dev/null +++ b/man/Plot.BioClim.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Plotting.R +\name{Plot.BioClim} +\alias{Plot.BioClim} +\title{Visualise bioclimatic raster data and overlay sf polygons if desired.} +\usage{ +Plot.BioClim( + BioClims, + Which = 1:19, + SF, + Water_Var = "Water Availability", + ncol = 3 +) +} +\arguments{ +\item{BioClims}{SpatRast object to visualise.} + +\item{Which}{Numeric. Which bioclimatic variable(s) to visualise.} + +\item{SF}{Optional. SF object which to overlay.} + +\item{Water_Var}{Optional, character. Name of water availability variable in the bioclimatic variables.} + +\item{ncol}{Number of columns for panel arrangement of plots} +} +\value{ +A ggplot2 object visualising a raster. +} +\description{ +Use the ggplot2 plotting engine to easily create visualisations of biolcimatic raster data - like the ones obtained using BioClim(...) - and overlay sf polygon data if desired. +} +\examples{ +BC_rast <- terra::rast(system.file("extdata", "CN_BC.nc", package = "KrigR")) +Plot.BioClim(BioClims = BC_rast, Water_Var = "Soil Moisture (0-7cm)") + +} +\seealso{ +\code{\link{BioClim}}. +} diff --git a/man/Plot.Covariates.Rd b/man/Plot.Covariates.Rd new file mode 100644 index 0000000..e9dc8ee --- /dev/null +++ b/man/Plot.Covariates.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Plotting.R +\name{Plot.Covariates} +\alias{Plot.Covariates} +\title{Visualise covariate raster data and overlay sf polygons if desired.} +\usage{ +Plot.Covariates(Covariates, SF, COL = viridis::cividis(100)) +} +\arguments{ +\item{Covariates}{List of length 2. Containing training resolution covariate SpatRast in slot 1 and target resolution covariate SpatRast in slot 2.} + +\item{SF}{Optional. SF object which to overlay.} + +\item{COL}{Colour palette.} +} +\value{ +A ggplot2 object visualising a raster. +} +\description{ +Use the ggplot2 plotting engine to easily create visualisations of covariate raster data - like the ones obtained using CovariateSetup(...) - and overlay sf polygon data if desired. +} +\examples{ +Cov_train <- terra::rast(system.file("extdata", "Covariates_Train.nc", package = "KrigR")) +Cov_target <- terra::rast(system.file("extdata", "Covariates_Target.nc", package = "KrigR")) +names(Cov_train) <- names(Cov_target) <- "GMTED2010 [m]" +Covariates <- list(Training = Cov_train, Target = Cov_target) +data("Jotunheimen_poly") +SF <- Jotunheimen_poly +Plot.Covariates(Covariates = Covariates, SF = SF) + +} +\seealso{ +\code{\link{CovariateSetup}}. +} diff --git a/man/Plot.Kriged.Rd b/man/Plot.Kriged.Rd new file mode 100644 index 0000000..b36d35c --- /dev/null +++ b/man/Plot.Kriged.Rd @@ -0,0 +1,49 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Plotting.R +\name{Plot.Kriged} +\alias{Plot.Kriged} +\title{Visualise kriged raster data and associated uncertainties and overlay sf polygons if desired.} +\usage{ +Plot.Kriged(Krigs, SF, Dates, Legend = "Air Temperature [K]") +} +\arguments{ +\item{Krigs}{List of length 2. Containing kriging prediction SpatRast in slot 1 and kriging standard error SpatRast in slot 2.} + +\item{SF}{Optional. SF object which to overlay.} + +\item{Dates}{Optional. Character vector of labels to apply to each layer of the SpatRast. By default, the content of the terra::time() field of the supplied SpatRast objects in list.} + +\item{Legend}{Colour label legend.} +} +\value{ +A ggplot2 object visualising a raster. +} +\description{ +Use the ggplot2 plotting engine to easily create visualisations of kriged raster data and associated uncertainties raster data - like the ones obtained using Kriging(...) - and overlay sf polygon data if desired. +} +\examples{ +CDS_rast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR")) +Cov_train <- terra::rast(system.file("extdata", "Covariates_Train.nc", package = "KrigR")) +Cov_target <- terra::rast(system.file("extdata", "Covariates_Target.nc", package = "KrigR")) +names(Cov_train) <- names(Cov_target) <- "GMTED2010 [m]" + +### kriging itself +ExtentKrig <- Kriging( + Data = CDS_rast[[1:2]], + Covariates_training = Cov_train, + Covariates_target = Cov_target, + Equation = "GMTED2010", + Cores = 2, + FileName = "KrigTest1", + FileExtension = ".nc", + Keep_Temporary = TRUE, + nmax = 40, + verbose = TRUE +) + +Plot.Kriged(Krigs = ExtentKrig) + +} +\seealso{ +\code{\link{Kriging}}. +} diff --git a/man/Plot.SpatRast.Rd b/man/Plot.SpatRast.Rd new file mode 100644 index 0000000..0c57e62 --- /dev/null +++ b/man/Plot.SpatRast.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Plotting.R +\name{Plot.SpatRast} +\alias{Plot.SpatRast} +\title{Visualise raster data and overlay sf polygons if desired.} +\usage{ +Plot.SpatRast( + SpatRast, + SF, + Dates, + Legend = "Air Temperature [K]", + COL = viridis::inferno(100) +) +} +\arguments{ +\item{SpatRast}{SpatRast object to visualise.} + +\item{SF}{Optional. SF object which to overlay.} + +\item{Dates}{Optional. Character vector of labels to apply to each layer of the SpatRast. By default, the content of the terra::time() field of the supplied SpatRast object.} + +\item{Legend}{Colour label legend.} + +\item{COL}{Colour palette.} +} +\value{ +A ggplot2 object visualising a raster. +} +\description{ +Use the ggplot2 plotting engine to easily create visualisations of raster data - like the ones obtained using CDownloadS(...) - and overlay sf polygon data if desired. +} +\examples{ +SpatRast <- terra::rast(system.file("extdata", "CentralNorway.nc", package = "KrigR"))[[1:2]] +data("Jotunheimen_poly") +SF <- Jotunheimen_poly +Plot.SpatRast(SpatRast = SpatRast, SF = SF) + +} +\seealso{ +\code{\link{CDownloadS}}. +} diff --git a/man/SummarizeRaster.Rd b/man/SummarizeRaster.Rd deleted file mode 100644 index f9072bb..0000000 --- a/man/SummarizeRaster.Rd +++ /dev/null @@ -1,17 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/misc.R -\name{SummarizeRaster} -\alias{SummarizeRaster} -\title{Summary of Raster file characteristics} -\usage{ -SummarizeRaster(Object_ras = NULL) -} -\arguments{ -\item{Object_ras}{A raster object.} -} -\value{ -A list containing information about the input raster. -} -\description{ -This function is called upon in the krigR function and summarizes Raster characteristics without carrying along the raster file itself. This is used to create lists tracking calls to the function krigR without bloating them too much. -} diff --git a/man/Temporal.Aggr.Rd b/man/Temporal.Aggr.Rd index 8a0d4c9..01e23a3 100644 --- a/man/Temporal.Aggr.Rd +++ b/man/Temporal.Aggr.Rd @@ -13,7 +13,8 @@ Temporal.Aggr( FUN, Cores, QueryTargetSteps, - TZone + TZone, + verbose = TRUE ) } \arguments{ @@ -34,6 +35,8 @@ Temporal.Aggr( \item{QueryTargetSteps}{Character. Target resolution steps} \item{TZone}{Character. Time zone for queried data.} + +\item{verbose}{Logical. Whether to print/message function progress in console or not.} } \value{ A SpatRaster diff --git a/man/Temporal.Cumul.Rd b/man/Temporal.Cumul.Rd index d1de4c7..78a6822 100644 --- a/man/Temporal.Cumul.Rd +++ b/man/Temporal.Cumul.Rd @@ -4,7 +4,14 @@ \alias{Temporal.Cumul} \title{Make cumulatively stored records into sequential ones} \usage{ -Temporal.Cumul(CDS_rast, CumulVar, BaseResolution, BaseStep, TZone) +Temporal.Cumul( + CDS_rast, + CumulVar, + BaseResolution, + BaseStep, + TZone, + verbose = TRUE +) } \arguments{ \item{CDS_rast}{SpatRaster} @@ -16,9 +23,13 @@ Temporal.Cumul(CDS_rast, CumulVar, BaseResolution, BaseStep, TZone) \item{BaseStep}{Numeric. Base time step of data set} \item{TZone}{Character. Time zone for queried data.} + +\item{verbose}{Logical. Whether to print/message function progress in console or not.} } \value{ A SpatRaster + +break apart sequence by UTC days and apply back-calculation per day in pblapply loop, for loop for each hour in each day } \description{ Takes a SpatRaster of cumulatively stored records and returns a SpatRaster of sequential counterparts diff --git a/man/buffer_Points.Rd b/man/buffer_Points.Rd deleted file mode 100644 index 2684b70..0000000 --- a/man/buffer_Points.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/misc.R -\name{buffer_Points} -\alias{buffer_Points} -\title{Square Buffers Around Point Data} -\usage{ -buffer_Points(Points = NULL, Buffer = 0.5, ID = "ID") -} -\arguments{ -\item{Points}{A data.frame containing geo-referenced points with Lat and Lon columns} - -\item{Buffer}{Identifies how big a rectangular buffer to draw around points. Expressed as centessimal degrees.} - -\item{ID}{Identifies which column in to use for creation of individual buffers.} -} -\value{ -A shape made up of individual square buffers around point-location input. -} -\description{ -Allow for drawing of buffer zones around point-location data for downloading and kriging of spatial data around point-locations. Overlapping individual buffers are merged. -} diff --git a/man/download_DEM.Rd b/man/download_DEM.Rd index 88b2e9b..ab622ec 100644 --- a/man/download_DEM.Rd +++ b/man/download_DEM.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/download.R +% Please edit documentation in R/DEPRECATED.R \name{download_DEM} \alias{download_DEM} \title{Downloading DEM data from USGS servers} diff --git a/man/download_ERA.Rd b/man/download_ERA.Rd index 5b1ab03..477a1e3 100644 --- a/man/download_ERA.Rd +++ b/man/download_ERA.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/download.R +% Please edit documentation in R/DEPRECATED.R \name{download_ERA} \alias{download_ERA} \title{Downloading ERA5(Land)-data from ECMWF servers} diff --git a/man/krigR.Rd b/man/krigR.Rd index b3ba54a..8d6d9bf 100644 --- a/man/krigR.Rd +++ b/man/krigR.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/krigr.R +% Please edit documentation in R/DEPRECATED.R \name{krigR} \alias{krigR} \title{(multi-core) Kriging} @@ -114,52 +114,52 @@ Use optional arguments such as Dir, FileName, Keep_Temporary, SingularTry, Krigi \dontrun{ ## THREE-STEP PROCESS (By Itself) # Downloading ERA5-Land air temperature reanalysis data in 12-hour intervals for 02/01/1995 - 04/01/1995 (DD/MM/YYYY). API User and Key in this example are non-functional. Substitute with your user number and key to run this example. -Extent <- extent(c(11.8,15.1,50.1,51.7)) # roughly the extent of Saxony +Extent <- extent(c(11.8, 15.1, 50.1, 51.7)) # roughly the extent of Saxony API_User <- "..." API_Key <- "..." State_Raw <- download_ERA( -Variable = "2m_temperature", -DataSet = "era5-land", -DateStart = "1995-01-02", -DateStop = "1995-01-04", -TResolution = "hour", -TStep = 12, -Extent = Extent, -API_User = API_User, -API_Key = API_Key + Variable = "2m_temperature", + DataSet = "era5-land", + DateStart = "1995-01-02", + DateStop = "1995-01-04", + TResolution = "hour", + TStep = 12, + Extent = Extent, + API_User = API_User, + API_Key = API_Key ) State_Raw # a raster brick with 6 layers at resolution of ~0.1° # Downloading GMTED2010-data at resolution and extent obtained by a call to download_ERA and a target resolution of .02. Covs_ls <- download_DEM( -Train_ras = State_Raw, -Target_res = .02, -Keep_Temporary = TRUE + Train_ras = State_Raw, + Target_res = .02, + Keep_Temporary = TRUE ) Covs_ls # a list with two elements: (1) GMTED 2010 data at training resolution, and (2) GMTED 2010 data aggregated as close as possible to a resolution of 0.02 # Kriging the data sets prepared with the previous functions. State_Krig <- krigR( Data = State_Raw, # data we want to krig as a raster object -Covariates_coarse = Covs_ls[[1]], # training covariate as a raster object -Covariates_fine = Covs_ls[[2]], # target covariate as a raster object + Covariates_coarse = Covs_ls[[1]], # training covariate as a raster object + Covariates_fine = Covs_ls[[2]], # target covariate as a raster object ) ## PIPELINE (From Scratch) #' # Downloading ERA5-Land air temperature reanalysis data in 12-hour intervals for 02/01/1995 - 04/01/1995 (DD/MM/YYYY), downloading and preparing GMTED 2010 covariate data, and kriging. API User and Key in this example are non-functional. Substitute with your user number and key to run this example. This example produces the same output as the example above. -Extent <- extent(c(11.8,15.1,50.1,51.7)) # roughly the extent of Saxony +Extent <- extent(c(11.8, 15.1, 50.1, 51.7)) # roughly the extent of Saxony API_User <- "..." API_Key <- "..." Pipe_Krig <- krigR( -Variable = "2m_temperature", -Type = "reanalysis", -DataSet = "era5-land", -DateStart = "1995-01-02", -DateStop = "1995-01-04", -TResolution = "hour",# -TStep = 12, -Extent = Extent, -API_User = API_User, -API_Key = API_Key, -Target_res = .02, + Variable = "2m_temperature", + Type = "reanalysis", + DataSet = "era5-land", + DateStart = "1995-01-02", + DateStop = "1995-01-04", + TResolution = "hour", # + TStep = 12, + Extent = Extent, + API_User = API_User, + API_Key = API_Key, + Target_res = .02, ) } diff --git a/man/mask_Shape.Rd b/man/mask_Shape.Rd deleted file mode 100644 index ffd18b4..0000000 --- a/man/mask_Shape.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/misc.R -\name{mask_Shape} -\alias{mask_Shape} -\title{Range Masking with Edge Support} -\usage{ -mask_Shape(base.map = NULL, Shape = NULL) -} -\arguments{ -\item{base.map}{A raster within which coverage should be identified} - -\item{Shape}{A polygon(-collection) whose coverage of the raster object is to be found.} -} -\value{ -A raster layer. -} -\description{ -Creating a raster mask identifying all cells in the original raster (`base.map`) which are at least partially covered by the supplied shapefile (`Shape`). -}