From d23aa1065a5a2f94eddf02beaa66e91256ef593b Mon Sep 17 00:00:00 2001 From: Alex Chubaty Date: Thu, 23 May 2024 23:12:34 -0600 Subject: [PATCH 1/2] shorten long lines identified in R CMD check also use `styler::style_active_file()` to automate formatting etc. --- R/downscale_core.R | 126 ++++++++++++++++++++--------------- R/plot_timeSeries_input.R | 80 +++++++++++----------- man/climr-package.Rd | 2 +- man/data-option-lists.Rd | 6 +- man/downscale_core.Rd | 36 +++++++--- man/plot_timeSeries_input.Rd | 21 +++--- 6 files changed, 157 insertions(+), 114 deletions(-) diff --git a/R/downscale_core.R b/R/downscale_core.R index d1f5384..2dc691c 100644 --- a/R/downscale_core.R +++ b/R/downscale_core.R @@ -3,10 +3,10 @@ #' @description #' `downscale_core()` is the engine for [`downscale()`]. #' It takes user-supplied high- and low-resolution rasters as input and downscales to user-specified point locations. -#' While less user-friendly than [`downscale()`], `downscale_core()` is more flexible in that users can supply their -#' own raster inputs. For example, a user could supply their own high-resolution climate map, instead of what is -#' available in climr, as the input to `refmap`. Another example is in downscaling a uniform warming level, as shown -#' in the example for this function. +#' While less user-friendly than [`downscale()`], `downscale_core()` is more flexible in that users can supply their +#' own raster inputs. For example, a user could supply their own high-resolution climate map, instead of what is +#' available in climr, as the input to `refmap`. Another example is in downscaling a uniform warming level, as shown +#' in the example for this function. #' #' @details #' We recommend [`downscale()`] for most purposes. @@ -31,7 +31,7 @@ #' @param out_spatial logical. Should a SpatVector be returned instead of a #' `data.frame`. #' @param plot character. If `out_spatial` is TRUE, the name of a variable to plot. -#' If the variable exists in `reference`, then its reference values will also be plotted. +#' If the variable exists in `reference`, then its reference values will also be plotted. #' Otherwise, reference January total precipitation (PPT01) values will be plotted. #' Defaults to no plotting (NULL). #' @@ -49,39 +49,55 @@ #' #' @export #' @examples -#' ## +#' ## #' library(terra) -#' xyz <- data.frame(lon = runif(10, -130, -106), lat = runif(10, 37, 50), elev = runif(10), id = 1:10) +#' xyz <- data.frame( +#' lon = runif(10, -130, -106), lat = runif(10, 37, 50), +#' elev = runif(10), id = 1:10 +#' ) #' #' ## get bounding box based on input points #' thebb <- get_bb(xyz) -#' +#' #' ## get database connection #' dbCon <- data_connect() -#' -#' # obtain the climatena 1961-1990 normals for the study area. +#' +#' # obtain the climatena 1961-1990 normals for the study area. #' refmap <- input_refmap(dbCon, thebb, reference = "refmap_climatena") -#' -#' # obtain the low-resolution climate data for a single gcm, 20-year period, and ssp scenario. -#' gcm_raw <- input_gcms(dbCon, thebb, list_gcms()[3], list_ssps()[1], period = list_gcm_periods()[2]) -#' +#' +#' # obtain the low-resolution climate data for a single gcm, 20-year period, and ssp scenario. +#' gcm_raw <- input_gcms(dbCon, thebb, list_gcms()[3], list_ssps()[1], +#' period = list_gcm_periods()[2] +#' ) +#' #' # downscale the GCM data -#' gcm_downscaled <- downscale_core(xyz = xyz, refmap = refmap, gcms = gcm_raw, vars = c("MAT", "PAS")) -#' -#' # create an input of uniform warming of 2 degrees Celsius and no precipitation change, for use as a null comparison to the GCM warming +#' gcm_downscaled <- downscale_core( +#' xyz = xyz, refmap = refmap, gcms = gcm_raw, +#' vars = c("MAT", "PAS") +#' ) +#' +#' # create an input of uniform warming of 2 degrees Celsius and no precipitation change, +#' # for use as a null comparison to the GCM warming #' null <- gcm_raw #' use the gcm input object as a template #' names(null) <- "null_2C" -#' names(null[[1]]) <- sapply(strsplit(names(null[[1]]), "_"), function(x) paste("null2C", x[2], x[3], "NA", "NA", "NA", "NA", sep="_")) -#' for(var in names(null[[1]])){ values(null[[1]][[var]]) <- if(length(grep("PPT", var)==1)) 1 else 2 } #' repopulate with the null values -#' +#' names(null[[1]]) <- sapply(strsplit(names(null[[1]]), "_"), function(x) { +#' paste("null2C", x[2], x[3], "NA", "NA", "NA", "NA", sep = "_") +#' }) +#' for (var in names(null[[1]])) { +#' values(null[[1]][[var]]) <- if (length(grep("PPT", var) == 1)) 1 else 2 +#' } # repopulate with the null values +#' #' # downscale the null values for variables of interest -#' null_downscaled <- downscale_core(xyz = xyz, refmap = refmap, gcms = null, vars = c("MAT", "PAS")) +#' null_downscaled <- downscale_core( +#' xyz = xyz, refmap = refmap, gcms = null, +#' vars = c("MAT", "PAS") +#' ) #' pool::poolClose(dbCon) -#' +#' downscale_core <- function(xyz, refmap, gcms = NULL, obs = NULL, gcm_ssp_ts = NULL, - gcm_hist_ts = NULL, obs_ts = NULL, return_refperiod = TRUE, - vars = sort(sprintf(c("PPT_%02d", "Tmax_%02d", "Tmin_%02d"), sort(rep(1:12, 3)))), - ppt_lr = FALSE, nthread = 1L, out_spatial = FALSE, plot = NULL) { + gcm_hist_ts = NULL, obs_ts = NULL, return_refperiod = TRUE, + vars = sort(sprintf(c("PPT_%02d", "Tmax_%02d", "Tmin_%02d"), sort(rep(1:12, 3)))), + ppt_lr = FALSE, nthread = 1L, out_spatial = FALSE, plot = NULL) { ## checks .checkDwnsclArgs( xyz, refmap, gcms, obs, gcm_ssp_ts, gcm_hist_ts, @@ -363,7 +379,7 @@ downscale_ <- function(xyz, refmap, gcms, gcm_ssp_ts, gcm_hist_ts, labels <- nm normal_ <- res # Reshape (melt / dcast) to obtain final form - #ref_dt <- tstrsplit(nm, "_") + # ref_dt <- tstrsplit(nm, "_") ref_dt <- data.table(VAR = nm) # setDT(ref_dt) # setnames(ref_dt, c("VAR")) @@ -469,38 +485,40 @@ process_one_climaterast <- function(climaterast, res, xyz, timeseries = FALSE, # Cropping will reduce the size of data to load in memory climaterast <- crop(climaterast, ex, snap = "out") - gc(reset = TRUE) ## free unused memory - + gc(reset = TRUE) ## free unused memory + climaterast <- try(extract(x = climaterast, y = xyz[, .(lon, lat)], method = "bilinear")) - + ## we may have run out of memory if there are MANY rasters - ## attempt to get only unique raster cell values + ## attempt to get only unique raster cell values ## (i.e. xyz may be at higher res than the climaterast leading to extracting the same values many times) if (is(climaterast, "try-error")) { if (grepl("bad_alloc", climaterast)) { message("System is out of memory to extract climate values for the supplied coordinates") - stop("Insufficient memory to downscale climate data for these many points/climate layers.\n", - " Try reducing number of points/layers.") - } - } - + stop( + "Insufficient memory to downscale climate data for these many points/climate layers.\n", + " Try reducing number of points/layers." + ) + } + } + # else { Ceres not sure what this is for but it's always causing fails # stop("Climate value extraction failed.", # "\n Please contact developers with a reproducible example and the error:\n", - # climaterast) + # climaterast) # } - + # Create match set to match with res names - - - labels <- vapply( - strsplit(nm, "_"), - \(x) { - paste0(x[2:3], collapse = "_") - }, - character(1) - ) - + + + labels <- vapply( + strsplit(nm, "_"), + \(x) { + paste0(x[2:3], collapse = "_") + }, + character(1) + ) + if (type %in% c("obs")) { ## Create match set to match with res names labels <- nm @@ -525,12 +543,12 @@ process_one_climaterast <- function(climaterast, res, xyz, timeseries = FALSE, } setDT(ref_dt) - if (type %in% c("obs","obs_ts")) { + if (type %in% c("obs", "obs_ts")) { if (timeseries) { setnames(ref_dt, c("DATASET", "VAR", "MONTH", "PERIOD")) set(ref_dt, j = "variable", value = nm) } else { - setnames(ref_dt, c("VAR","MONTH")) + setnames(ref_dt, c("VAR", "MONTH")) set(ref_dt, j = "variable", value = nm) set(ref_dt, j = "PERIOD", value = "2001_2020") } @@ -686,7 +704,7 @@ unpackRasters <- function(ras) { if (!inherits(xyz$id, colTypes)) { stop("'xyz$id' must be an column of type ", paste(colTypes, collapse = ", ")) } - + return(xyz) } @@ -701,18 +719,18 @@ unpackRasters <- function(ras) { obs_ts = NULL, return_refperiod = FALSE, out_spatial = FALSE, plot = NULL, vars = list_vars()) { vars <- match.arg(vars, list_vars(), several.ok = TRUE) - + if (!return_refperiod %in% c(TRUE, FALSE)) { stop("'return_refperiod' must be TRUE or FALSE") } if (!out_spatial %in% c(TRUE, FALSE)) { stop("'out_spatial' must be TRUE or FALSE") } - + plot <- if (!is.null(plot)) { - match.arg(plot,list_vars()) + match.arg(plot, list_vars()) } - + if (!isTRUE(attr(refmap, "builder") == "climr")) { stop( "Please use `input_refmap` function to create `refmap`.", diff --git a/R/plot_timeSeries_input.R b/R/plot_timeSeries_input.R index 1647de3..bf6b0d6 100644 --- a/R/plot_timeSeries_input.R +++ b/R/plot_timeSeries_input.R @@ -1,62 +1,66 @@ #' Input data for the time series climate change plot #' #' @description -#' Input data for the [`plot_timeSeries()`] function. Since these inputs are time-consuming to generate, +#' Input data for the [`plot_timeSeries()`] function. Since these inputs are time-consuming to generate, #' the purpose of conducting the generation of the input table in a separate function is to allow users -#' to make multiple calls to [`plot_timeSeries()`] (e.g., for comparing different climate variables) +#' to make multiple calls to [`plot_timeSeries()`] (e.g., for comparing different climate variables) #' without needing to generate the inputs each time. -#' +#' #' @details #' This function generates standardized inputs for one or multiple locations at any spatial scale. -#' If multiple locations are specified, the output is the average of the climate variables for all locations. -#' -#' Downloads of GCM time series take a long time. The `plot_timeSeries_input()` function can take >1hr -#' to run for the first time it is called for a location. We are looking into ways to speed this up, but until then -#' we recommend users dedicate some time to run this function in background. Once the time series are cached, they -#' don't need to be downloaded again. +#' If multiple locations are specified, the output is the average of the climate variables for all locations. +#' +#' Downloads of GCM time series take a long time. The `plot_timeSeries_input()` function can take >1hr +#' to run for the first time it is called for a location. We are looking into ways to speed this up, but until then +#' we recommend users dedicate some time to run this function in background. Once the time series are cached, they +#' don't need to be downloaded again. #' #' @template xyz #' @inheritParams downscale #' @template vars #' #' @return `data.table` of average downscaled climate variables for all locations. -#' +#' #' @examples -#' if(FALSE){ -#' # data frame of arbitrary points -#' my_points <- data.frame(lon = c(-127.7300,-127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) -#' -#' # generate the input data -#' my_data <- plot_timeSeries_input(my_points) -#' -#' # use the input to create a plot -#' plot_timeSeries(my_data, variable1 = "Tmin_sm") +#' if (FALSE) { +#' # data frame of arbitrary points +#' my_points <- data.frame( +#' lon = c(-127.7300, -127.7500), +#' lat = c(55.34114, 55.25), +#' elev = c(711, 500), +#' id = 1:2 +#' ) +#' +#' # generate the input data +#' my_data <- plot_timeSeries_input(my_points) +#' +#' # use the input to create a plot +#' plot_timeSeries(my_data, variable1 = "Tmin_sm") #' } -#' #' +#' #' #' @export - plot_timeSeries_input <- function( - xyz, + xyz, gcms = list_gcms(), ssps = list_ssps(), max_run = 10, - obs_ts_dataset = c("cru.gpcc", "climatena"), + obs_ts_dataset = c("cru.gpcc", "climatena"), obs_years = 1901:2022, - gcm_hist_years = 1850:2014, - gcm_ssp_years = 2015:2100, - vars = list_vars() -) { - data <- downscale(xyz = xyz, - gcms = gcms, - ssps = ssps, - max_run = max_run, - obs_ts_dataset = obs_ts_dataset, - obs_years = obs_years, - gcm_hist_years = gcm_hist_years, - gcm_ssp_years = gcm_ssp_years, - vars = vars + gcm_hist_years = 1850:2014, + gcm_ssp_years = 2015:2100, + vars = list_vars()) { + data <- downscale( + xyz = xyz, + gcms = gcms, + ssps = ssps, + max_run = max_run, + obs_ts_dataset = obs_ts_dataset, + obs_years = obs_years, + gcm_hist_years = gcm_hist_years, + gcm_ssp_years = gcm_ssp_years, + vars = vars ) - data.agg <- data[, lapply(.SD, mean), by = .(GCM, SSP, RUN, PERIOD, DATASET), .SDcols = -c("id", "GCM", "SSP", "RUN", "PERIOD", "DATASET")] + data.agg <- data[, lapply(.SD, mean), by = .(GCM, SSP, RUN, PERIOD, DATASET), + .SDcols = -c("id", "GCM", "SSP", "RUN", "PERIOD", "DATASET")] return(data.agg) } - diff --git a/man/climr-package.Rd b/man/climr-package.Rd index 18d89d3..17668bd 100644 --- a/man/climr-package.Rd +++ b/man/climr-package.Rd @@ -15,7 +15,7 @@ Authors: \itemize{ \item Colin Mahony \email{Colin.Mahony@gov.bc.ca} (\href{https://orcid.org/0000-0002-6111-5675}{ORCID}) \item Bruno Tremblay \email{bruno@boostao.ca} (\href{https://orcid.org/0000-0002-2945-356X}{ORCID}) - \item Ceres Barros \email{ceres.barros@gov.bc.ca} (\href{https://orcid.org/0000-0003-4036-977X}{ORCID}) + \item Ceres Barros \email{ceres.barros@nrcan-rncan.gc.ca} (\href{https://orcid.org/0000-0003-4036-977X}{ORCID}) } Other contributors: diff --git a/man/data-option-lists.Rd b/man/data-option-lists.Rd index 59eb891..a4d8169 100644 --- a/man/data-option-lists.Rd +++ b/man/data-option-lists.Rd @@ -54,9 +54,9 @@ a character vector. \code{list_run} lists available runs for a given GCM. -\code{list_refmaps} lists available normals. +\code{list_refmaps} lists available reference maps of gridded climate normals -\code{list_obs_periods} lists available obs periods +\code{list_obs_periods} lists available observational periods \code{list_vars} lists climate variables @@ -67,7 +67,7 @@ a character vector. \code{list_gcm_hist_years} lists available years for obs projections' time series } \details{ -Currently available normals (\code{list_refmaps()}) are: +Currently available reference maps of gridded climate normals (\code{list_refmaps()}) are: \itemize{ \item "refmap_climatena" for Climate NA derived normals \item "refmap_prism" for British Columbia PRISM climatologies derived normals diff --git a/man/downscale_core.Rd b/man/downscale_core.Rd index 236c396..b934c2a 100644 --- a/man/downscale_core.Rd +++ b/man/downscale_core.Rd @@ -78,9 +78,12 @@ in the example for this function. We recommend \code{\link[=downscale]{downscale()}} for most purposes. } \examples{ -## +## library(terra) -xyz <- data.frame(lon = runif(10, -130, -106), lat = runif(10, 37, 50), elev = runif(10), id = 1:10) +xyz <- data.frame( + lon = runif(10, -130, -106), lat = runif(10, 37, 50), + elev = runif(10), id = 1:10 +) ## get bounding box based on input points thebb <- get_bb(xyz) @@ -88,23 +91,36 @@ thebb <- get_bb(xyz) ## get database connection dbCon <- data_connect() -# obtain the climatena 1961-1990 normals for the study area. +# obtain the climatena 1961-1990 normals for the study area. refmap <- input_refmap(dbCon, thebb, reference = "refmap_climatena") -# obtain the low-resolution climate data for a single gcm, 20-year period, and ssp scenario. -gcm_raw <- input_gcms(dbCon, thebb, list_gcms()[3], list_ssps()[1], period = list_gcm_periods()[2]) +# obtain the low-resolution climate data for a single gcm, 20-year period, and ssp scenario. +gcm_raw <- input_gcms(dbCon, thebb, list_gcms()[3], list_ssps()[1], + period = list_gcm_periods()[2] +) # downscale the GCM data -gcm_downscaled <- downscale_core(xyz = xyz, refmap = refmap, gcms = gcm_raw, vars = c("MAT", "PAS")) +gcm_downscaled <- downscale_core( + xyz = xyz, refmap = refmap, gcms = gcm_raw, + vars = c("MAT", "PAS") +) -# create an input of uniform warming of 2 degrees Celsius and no precipitation change, for use as a null comparison to the GCM warming +# create an input of uniform warming of 2 degrees Celsius and no precipitation change, +# for use as a null comparison to the GCM warming null <- gcm_raw #' use the gcm input object as a template names(null) <- "null_2C" -names(null[[1]]) <- sapply(strsplit(names(null[[1]]), "_"), function(x) paste("null2C", x[2], x[3], "NA", "NA", "NA", "NA", sep="_")) -for(var in names(null[[1]])){ values(null[[1]][[var]]) <- if(length(grep("PPT", var)==1)) 1 else 2 } #' repopulate with the null values +names(null[[1]]) <- sapply(strsplit(names(null[[1]]), "_"), function(x) { + paste("null2C", x[2], x[3], "NA", "NA", "NA", "NA", sep = "_") +}) +for (var in names(null[[1]])) { + values(null[[1]][[var]]) <- if (length(grep("PPT", var) == 1)) 1 else 2 +} # repopulate with the null values # downscale the null values for variables of interest -null_downscaled <- downscale_core(xyz = xyz, refmap = refmap, gcms = null, vars = c("MAT", "PAS")) +null_downscaled <- downscale_core( + xyz = xyz, refmap = refmap, gcms = null, + vars = c("MAT", "PAS") +) pool::poolClose(dbCon) } diff --git a/man/plot_timeSeries_input.Rd b/man/plot_timeSeries_input.Rd index 4a44763..93d02af 100644 --- a/man/plot_timeSeries_input.Rd +++ b/man/plot_timeSeries_input.Rd @@ -70,15 +70,20 @@ we recommend users dedicate some time to run this function in background. Once t don't need to be downloaded again. } \examples{ -if(FALSE){ -# data frame of arbitrary points -my_points <- data.frame(lon = c(-127.7300,-127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) +if (FALSE) { + # data frame of arbitrary points + my_points <- data.frame( + lon = c(-127.7300, -127.7500), + lat = c(55.34114, 55.25), + elev = c(711, 500), + id = 1:2 + ) -# generate the input data -my_data <- plot_timeSeries_input(my_points) + # generate the input data + my_data <- plot_timeSeries_input(my_points) -# use the input to create a plot -plot_timeSeries(my_data, variable1 = "Tmin_sm") + # use the input to create a plot + plot_timeSeries(my_data, variable1 = "Tmin_sm") } -#' +#' } From 336a513db69933670d5843c822190b28282841db Mon Sep 17 00:00:00 2001 From: Alex Chubaty Date: Thu, 23 May 2024 23:23:55 -0600 Subject: [PATCH 2/2] deal with long lines in plot_timeSeries example + cleanup --- R/plot_timeSeries.R | 651 +++++++++++++++++++++-------------------- man/plot_timeSeries.Rd | 75 ++--- 2 files changed, 372 insertions(+), 354 deletions(-) diff --git a/R/plot_timeSeries.R b/R/plot_timeSeries.R index 386fe1d..747e599 100644 --- a/R/plot_timeSeries.R +++ b/R/plot_timeSeries.R @@ -12,91 +12,98 @@ #' All global climate model anomalies are bias-corrected to the 1961-1990 reference period normals. #' #' @details -#' The input table `X` provides climate data for a single location or the average of multiple locations. +#' The input table `X` provides climate data for a single location or the average of multiple locations. #' The purpose of conducting the generation of the input table in a separate function is to allow users -#' to make multiple calls to [`plot_timeSeries()`] without needing to generate the inputs each time. -#' -#' Some combinations of `var1` and `var2` are not compatible or meaningful. -#' Examples of meaningful combinations are winter vs summer values of the same climate var or minimum vs. -#' maximum temperatures. -#' -#' Downloads of GCM time series take a long time. The `plot_timeSeries_input()` function can take >1hr -#' to run for the first time it is called for a location. We are looking into ways to speed this up, but until then -#' we recommend users dedicate some time to run this function in background. Once the time series are cached, they -#' don't need to be downloaded again. +#' to make multiple calls to [`plot_timeSeries()`] without needing to generate the inputs each time. #' -#' @param X A `data.table` object produced using the function [`plot_timeSeries_input()`]. This -#' table can include more models, scenarios, and variables than are used in individual calls to +#' Some combinations of `var1` and `var2` are not compatible or meaningful. +#' Examples of meaningful combinations are winter vs summer values of the same climate var or minimum vs. +#' maximum temperatures. +#' +#' Downloads of GCM time series take a long time. The `plot_timeSeries_input()` function can take >1hr +#' to run for the first time it is called for a location. We are looking into ways to speed this up, but until then +#' we recommend users dedicate some time to run this function in background. Once the time series are cached, they +#' don't need to be downloaded again. +#' +#' @param X A `data.table` object produced using the function [`plot_timeSeries_input()`]. This +#' table can include more models, scenarios, and variables than are used in individual calls to #' [`plot_timeSeries()`]. #' @inheritParams downscale #' @param var1 character. A climate var. options are [`list_vars()`]. -#' @param var2 character. A second climate var to plot in combination with `var1`. +#' @param var2 character. A second climate var to plot in combination with `var1`. #' options are [`list_vars()`]. -#' @param showObserved logical. Plot a time series of observed climate. -#' @param showrange logical. Plot a shaded region indicating the minimum and maximum of the -#' selected ensemble of GCM simulations for each selected scenario. -#' @param yfit logical. Set the range of the y axis to the range of the visible data. If `FALSE` -#' the y axis is the range of all values of `var1` (and `var2` if applicable) in the -#' input table defined by `X`. +#' @param showObserved logical. Plot a time series of observed climate. +#' @param showrange logical. Plot a shaded region indicating the minimum and maximum of the +#' selected ensemble of GCM simulations for each selected scenario. +#' @param yfit logical. Set the range of the y axis to the range of the visible data. If `FALSE` +#' the y axis is the range of all values of `var1` (and `var2` if applicable) in the +#' input table defined by `X`. #' @param cex Numeric. The magnification factor for text size. Default is 1. -#' @param mar A numerical vector of length 4, giving the margin sizes in number of lines of text: c(bottom, left, +#' @param mar A numerical vector of length 4, giving the margin sizes in number of lines of text: c(bottom, left, #' top, right). The default is c(3,3,0.1,4). -#' @param showmean logical. Plot the ensemble mean time series. Multi-model ensemble means are -#' calculated from the mean of simulations for each model. -#' @param compile logical. Compile multiple global climate models into a multi-model ensemble. -#' If `FALSE` the single-model ensembles are plotted individually. -#' @param simplify logical. Simplify the ensemble range and mean using a smoothing spline. +#' @param showmean logical. Plot the ensemble mean time series. Multi-model ensemble means are +#' calculated from the mean of simulations for each model. +#' @param compile logical. Compile multiple global climate models into a multi-model ensemble. +#' If `FALSE` the single-model ensembles are plotted individually. +#' @param simplify logical. Simplify the ensemble range and mean using a smoothing spline. #' @param refline logical. Plot the 1961-1990 reference period mean for the selected var -#' and extend this line to the year 2100 as a visual reference. -#' @param refline.obs logical. Plot the 1961-1990 reference period mean for the observational data. -#' This should be the same as the reference line for the GCM time series. -#' @param pal character. color palette. Options are "scenario", for use when comparing scenarios, -#' and "gcms", for use when comparing GCMs. -#' @param label.endyear logical. Add a label of the final year of the observational time series. -#' @param endlabel character. Add a label to the end of each simulated time series. Options +#' and extend this line to the year 2100 as a visual reference. +#' @param refline.obs logical. Plot the 1961-1990 reference period mean for the observational data. +#' This should be the same as the reference line for the GCM time series. +#' @param pal character. color palette. Options are "scenario", for use when comparing scenarios, +#' and "gcms", for use when comparing GCMs. +#' @param label.endyear logical. Add a label of the final year of the observational time series. +#' @param endlabel character. Add a label to the end of each simulated time series. Options #' are "change", to indicate the change in year 2100 relative to the 1961-1990 baseline, or "gcms" -#' to indicate the global climate model. -#' @param yearmarkers logical. Add white points to the observational time series as a visual aid. +#' to indicate the global climate model. +#' @param yearmarkers logical. Add white points to the observational time series as a visual aid. #' @param yearlines logical. Add vertical lines on every fifth year as a visual reference -#' @param legend_pos character. Position of the legend. Viable options are c("bottomright", -#' "bottomleft", "topleft", and "topright"). +#' @param legend_pos character. Position of the legend. Viable options are "bottomright", +#' "bottomleft", "topleft", and "topright". #' #' @return NULL. Draws a plot in the active graphics device. #' #' @examples -#' if(FALSE){ -#' # data frame of arbitrary points -#' my_points <- data.frame(lon = c(-127.7300,-127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) +#' if (FALSE) { +#' # data frame of arbitrary points +#' my_points <- data.frame( +#' lon = c(-127.7300, -127.7500), +#' lat = c(55.34114, 55.25), +#' elev = c(711, 500), +#' id = 1:2 +#' ) +#' +#' # generate the input data +#' my_data <- plot_timeSeries_input(my_points) #' -#' # generate the input data -#' my_data <- plot_timeSeries_input(my_points) +#' # use the input to create a plot +#' plot_timeSeries(my_data, var1 = "Tmin_sm") #' -#' # use the input to create a plot -#' plot_timeSeries(my_data, var1 = "Tmin_sm") +#' # compare observational time series +#' plot_timeSeries(my_data, var1 = "Tmin_sm", obs_ts_dataset = c("cru.gpcc", "climatena")) #' -#' # compare observational time series -#' plot_timeSeries(my_data, var1 = "Tmin_sm", obs_ts_dataset = c("cru.gpcc", "climatena")) -#' -#' # compare mean daily minimum and maximum temperatures -#' plot_timeSeries(my_data, var1 = "Tmin_sm", var2 = "Tmax_sm") +#' # compare mean daily minimum and maximum temperatures +#' plot_timeSeries(my_data, var1 = "Tmin_sm", var2 = "Tmax_sm") #' -#' # compare summer and winter temperatures (without simplifying the ensemble range) -#' plot_timeSeries(my_data, var1 = "Tmax_sm", var2 = "Tmax_wt", simplify = FALSE) +#' # compare summer and winter temperatures (without simplifying the ensemble range) +#' plot_timeSeries(my_data, var1 = "Tmax_sm", var2 = "Tmax_wt", simplify = FALSE) #' -#' # compare global climate models -#' plot_timeSeries(my_data, gcms = list_gcms()[c(7,13)], pal = "gcms", ssps = list_ssps()[2], showmean = FALSE, compile = FALSE, simplify = FALSE, endlabel = "gcms", mar=c(3,3,0.1,6), showObserved = FALSE) +#' # compare global climate models +#' plot_timeSeries(my_data, gcms = list_gcms()[c(7, 13)], pal = "gcms", +#' ssps = list_ssps()[2], showmean = FALSE, compile = FALSE, +#' simplify = FALSE, endlabel = "gcms", mar = c(3, 3, 0.1, 6), +#' showObserved = FALSE) #' -#' # export plot to a temporary directory, including a title -#' figDir <- tempdir() -#' png( -#' filename = file.path(figDir, "plot_test.png"), type = "cairo", units = "in", -#' width = 6, height = 5, pointsize = 10, res = 300 -#' ) -#' plot_timeSeries(my_data, var1 = "Tmin_sm", mar=c(3,3,2,4)) -#' title("Historical and projected summer night-time warming in the Bulkley Valley, BC") -#' dev.off() +#' # export plot to a temporary directory, including a title +#' figDir <- tempdir() +#' png( +#' filename = file.path(figDir, "plot_test.png"), type = "cairo", units = "in", +#' width = 6, height = 5, pointsize = 10, res = 300 +#' ) +#' plot_timeSeries(my_data, var1 = "Tmin_sm", mar = c(3, 3, 2, 4)) +#' title("Historical and projected summer night-time warming in the Bulkley Valley, BC") +#' dev.off() #' } -#' #' #' @importFrom scales alpha #' @importFrom stinepack stinterp @@ -109,299 +116,303 @@ plot_timeSeries <- function( var2 = NULL, showObserved = TRUE, obs_ts_dataset = "climatena", - gcms = list_gcms()[c(1,4,5,6,7,10,11,12)], + gcms = list_gcms()[c(1, 4, 5, 6, 7, 10, 11, 12)], ssps = list_ssps()[1:3], - showrange = TRUE, + showrange = TRUE, yfit = TRUE, cex = 1, - mar=c(3,3,0.1,4), + mar = c(3, 3, 0.1, 4), showmean = TRUE, compile = TRUE, simplify = TRUE, refline = FALSE, refline.obs = TRUE, pal = "scenario", - label.endyear = FALSE, - endlabel = "change", + label.endyear = FALSE, + endlabel = "change", yearmarkers = TRUE, yearlines = FALSE, legend_pos = "topleft") { if (!requireNamespace("scales", quietly = TRUE)) { stop("package scales must be installed to use this function") - } else if (!requireNamespace("stinepack", quietly = TRUE)) { + } else if (!requireNamespace("stinepack", quietly = TRUE)) { stop("package stinepack must be installed to use this function") } else { - data("variables", envir = environment()) + data("variables", envir = environment()) - # Scenario definitions - scenarios.selected <- c("historical", ssps) - scenarios <- c("historical", list_ssps()) - scenario.names <- c("Historical simulations", "SSP1-2.6", "SSP2-4.5", "SSP3-7.0", "SSP5-8.5") - pal.scenario <- c("gray60", "dodgerblue4", "seagreen", "darkorange3", "darkred") # these roughly match the IPCC standard scenario colours. - - # GCM color palette (from https://mk.bcgsc.ca/colorblind/) - pal.gcms <- c("#004949","#009292","#ff6db6","#ffb6db", "#490092","#006ddb","#b66dff","#6db6ff","#b6dbff", "#920000","#924900","#db6d00","#24ff24") - - # function for specifying the color - colSelect <- function(scenario, gcms){ - if(missing(gcms)){ - col = pal.scenario[which(scenarios==scenario)] - } else { - col = if(pal=="gcms") pal.gcms[which(list_gcms()==gcms)] else pal.scenario[which(scenarios==scenario)] - } - return(col) - } - - # yeartime definitions - monthcodes <- c("01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12") - seasonmonth.mat <- matrix(monthcodes[c(12, 1:11)],4, byrow=TRUE) - seasons <- c("wt", "sp", "sm", "at") - season.names <- c("Winter", "Spring", "Summer", "Autumn") - yeartimes <- c(seasons, monthcodes) - yeartime.names <- c(season.names, month.name) - - # ensemble statistics definitions - ensstats <- c("ensmin", "ensmax", "ensmean") - - ## Assemble the data that will be used in the plot (for setting the ylim) - alldata <- X[, if(is.null(var2)) get(var1) else c(get(var1), get(var2))] # a vector of all potential data on the plot for setting the ylim (y axis range) - visibledata <- X[GCM%in%c(NA, gcms) & SSP%in%c(NA, ssps), (if(is.null(var2)) get(var1) else c(get(var1), get(var2)))] # a vector of all visible data on the plot for setting the ylim (y axis range) - - # components of the var - nums <- if(is.null(var2)) 1 else 1:2 #nums is the set of var numbers (var1 or var2) (this is used later on as well) - for(num in nums){ - var <- get(paste("var",num,sep="")) - assign(paste0("yeartime", num), as.character(variables[Code == var, "Code_Time"]) ) - assign(paste0("element", num), as.character(variables[Code == var, "Code_Element"]) ) + # Scenario definitions + scenarios.selected <- c("historical", ssps) + scenarios <- c("historical", list_ssps()) + scenario.names <- c("Historical simulations", "SSP1-2.6", "SSP2-4.5", "SSP3-7.0", "SSP5-8.5") + pal.scenario <- c("gray60", "dodgerblue4", "seagreen", "darkorange3", "darkred") # these roughly match the IPCC standard scenario colours. + + # GCM color palette (from https://mk.bcgsc.ca/colorblind/) + pal.gcms <- c("#004949", "#009292", "#ff6db6", "#ffb6db", "#490092", "#006ddb", "#b66dff", "#6db6ff", "#b6dbff", "#920000", "#924900", "#db6d00", "#24ff24") + + # function for specifying the color + colSelect <- function(scenario, gcms) { + if (missing(gcms)) { + col <- pal.scenario[which(scenarios == scenario)] + } else { + col <- if (pal == "gcms") pal.gcms[which(list_gcms() == gcms)] else pal.scenario[which(scenarios == scenario)] } - - # PLOT - par(mfrow=c(1,1), mar=mar, mgp=c(1.75, 0.25, 0), cex=cex) - # y axis title. - if(is.null(var2)){ #if there is no second var - ylab <- paste(yeartime.names[which(yeartimes==yeartime1)], variables[Code==var1, "Element"]) - } else - if(element1==element2){ #if both variables have the same element - ylab <- variables[Code==var1, "Element"] - } else - if(yeartime1==yeartime2){ #if both variables have the same yeartime - ylab <- paste(yeartime.names[which(yeartimes==yeartime1)], variables[Code==var1, "Element"], "or", variables[Code==var2, "Element"]) - } else { #if variables1 and 2 have different elements and yeartimes - ylab <- paste(yeartime.names[which(yeartimes==yeartime1)], variables[Code==var1, "Element"], "or", variables[Code==var2, "Element"]) - } - plot(0, col="white", xlim=c(1900, 2100), ylim=range(if(yfit) visibledata else alldata, na.rm = TRUE), xaxs="i", xaxt="n", tck=0, xlab="", ylab=ylab) - axis(1, at=seq(1850,2100,25), labels = seq(1850,2100,25), tck=0) - - num <- 1 - for(num in nums){ - yeartime <- get(paste("yeartime",num,sep="")) - element <- get(paste("element",num,sep="")) - var <- get(paste("var",num,sep="")) - - # function for plotting time series for gcms or compiled ensemble - plot.ensemble <- function(x) { - # scenario <- scenarios.selected[1] - temp.historical <- x[is.na(SSP), c("PERIOD", "RUN", var), with=FALSE] - x.historical <- as.numeric(temp.historical[, .(min = min(get(var))), by = PERIOD][["PERIOD"]]) - ensmin.historical <- temp.historical[RUN!="ensembleMean", .(min = min(get(var), na.rm=TRUE)), by = PERIOD][["min"]] - ensmax.historical <- temp.historical[RUN!="ensembleMean", .(max = max(get(var), na.rm=TRUE)), by = PERIOD][["max"]] - ensmean.historical <- temp.historical[RUN=="ensembleMean", .(mean = mean(get(var), na.rm=TRUE)), by = PERIOD][["mean"]] #calculate mean using the single-model ensembleMean, to ensure that the mean isn't biased towards models with more runs - - for(scenario in scenarios.selected[order(c(1,4,5,3,2)[which(scenarios%in%scenarios.selected)])][-1]){ - temp <- x[SSP==scenario, c("PERIOD", "RUN", var), with=FALSE] - x.temp <- as.numeric(temp[, .(min = min(get(var))), by = PERIOD][["PERIOD"]]) - ensmin.temp <- temp[RUN!="ensembleMean", .(min = min(get(var), na.rm=TRUE)), by = PERIOD][["min"]] - ensmax.temp <- temp[RUN!="ensembleMean", .(max = max(get(var), na.rm=TRUE)), by = PERIOD][["max"]] - ensmean.temp <- temp[RUN=="ensembleMean", .(mean = mean(get(var), na.rm=TRUE)), by = PERIOD][["mean"]] #calculate mean using the single-model ensembleMean, to ensure that the mean isn't biased towards models with more runs - assign(paste("x", scenario, sep="."), c(x.historical[length(x.historical)], x.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. - assign(paste("ensmin", scenario, sep="."), c(ensmin.historical[length(ensmin.historical)], ensmin.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. - assign(paste("ensmax", scenario, sep="."), c(ensmax.historical[length(ensmax.historical)], ensmax.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. - assign(paste("ensmean", scenario, sep="."), c(ensmean.historical[length(ensmean.historical)], ensmean.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. - } - - if(showrange) { - if(simplify==FALSE){ - for(scenario in scenarios.selected[order(c(1,4,5,3,2)[which(scenarios%in%scenarios.selected)])]){ - x3 <- get(paste("x", scenario, sep=".")) - polygon(c(x3, rev(x3)), c(get(paste("ensmin", scenario, sep=".")), rev(get(paste("ensmax", scenario, sep=".")))), col=alpha(colSelect(scenario, gcms), 0.35), border=colSelect(scenario, gcms)) - } - } else { - scenarios.select <- scenarios.selected[order(c(1,4,5,3,2)[which(scenarios%in%scenarios.selected)])][-1] - for(scenario in scenarios.select){ - - if(scenario==scenarios.select[1]){ # we need to run spline through the historical/projected transition - x4 <- c(x.historical, get(paste("x", scenario, sep="."))[-1]) - y.ensmin <- c(ensmin.historical, get(paste("ensmin", scenario, sep="."))[-1]) - y.ensmax <- c(ensmax.historical, get(paste("ensmax", scenario, sep="."))[-1]) - s.ensmin <- smooth.spline(x4,y.ensmin, df=8) - s.ensmax <- smooth.spline(x4,y.ensmax, df=8) - subset.hist <- which(x4%in%x.historical) - subset.proj <- which(x4%in%get(paste("x", scenario, sep="."))) - polygon(c(s.ensmin$x[subset.hist], rev(s.ensmax$x[subset.hist])), c(s.ensmin$y[subset.hist], rev(s.ensmax$y[subset.hist])), col=alpha(if(pal=="gcms") colSelect(scenario, gcms) else pal.scenario[which(scenarios=="historical")], 0.35), border=NA) - polygon(c(s.ensmin$x[subset.proj], rev(s.ensmax$x[subset.proj])), c(s.ensmin$y[subset.proj], rev(s.ensmax$y[subset.proj])), col=alpha(colSelect(scenario, gcms), 0.35), border=NA) - } else { # this second routine uses interpolation splines so that the starting point for all scenarios is the same - x5 <- c(x.historical, get(paste("x", scenario, sep="."))[-1]) - y.ensmin2 <- c(ensmin.historical, get(paste("ensmin", scenario, sep="."))[-1]) - y.ensmax2 <- c(ensmax.historical, get(paste("ensmax", scenario, sep="."))[-1]) - s.ensmin2 <- smooth.spline(x5,y.ensmin2, df=8) - s.ensmax2 <- smooth.spline(x5,y.ensmax2, df=8) - knots.hist <- which(x5%in%c(seq(1860, 2000, 20), 2014)) - knots.proj <- which(x5%in%c(seq(2030, 2090, 20), 2100)) - s.ensmin3 <- stinterp(x5[c(knots.hist, knots.proj)],c(s.ensmin$y[knots.hist], s.ensmin2$y[knots.proj]), x5) - s.ensmax3 <- stinterp(x5[c(knots.hist, knots.proj)],c(s.ensmax$y[knots.hist], s.ensmax2$y[knots.proj]), x5) - polygon(c(s.ensmin3$x[subset.proj], rev(s.ensmax3$x[subset.proj])), c(s.ensmin3$y[subset.proj], rev(s.ensmax3$y[subset.proj])), col=alpha(colSelect(scenario, gcms), 0.35), border=NA) - } + return(col) + } + + # yeartime definitions + monthcodes <- c("01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12") + seasonmonth.mat <- matrix(monthcodes[c(12, 1:11)], 4, byrow = TRUE) + seasons <- c("wt", "sp", "sm", "at") + season.names <- c("Winter", "Spring", "Summer", "Autumn") + yeartimes <- c(seasons, monthcodes) + yeartime.names <- c(season.names, month.name) + + # ensemble statistics definitions + ensstats <- c("ensmin", "ensmax", "ensmean") + + ## Assemble the data that will be used in the plot (for setting the ylim) + alldata <- X[, if (is.null(var2)) get(var1) else c(get(var1), get(var2))] # a vector of all potential data on the plot for setting the ylim (y axis range) + visibledata <- X[GCM %in% c(NA, gcms) & SSP %in% c(NA, ssps), (if (is.null(var2)) get(var1) else c(get(var1), get(var2)))] # a vector of all visible data on the plot for setting the ylim (y axis range) + + # components of the var + nums <- if (is.null(var2)) 1 else 1:2 # nums is the set of var numbers (var1 or var2) (this is used later on as well) + for (num in nums) { + var <- get(paste("var", num, sep = "")) + assign(paste0("yeartime", num), as.character(variables[Code == var, "Code_Time"])) + assign(paste0("element", num), as.character(variables[Code == var, "Code_Element"])) + } + + # PLOT + par(mfrow = c(1, 1), mar = mar, mgp = c(1.75, 0.25, 0), cex = cex) + # y axis title. + if (is.null(var2)) { # if there is no second var + ylab <- paste(yeartime.names[which(yeartimes == yeartime1)], variables[Code == var1, "Element"]) + } else if (element1 == element2) { # if both variables have the same element + ylab <- variables[Code == var1, "Element"] + } else if (yeartime1 == yeartime2) { # if both variables have the same yeartime + ylab <- paste(yeartime.names[which(yeartimes == yeartime1)], variables[Code == var1, "Element"], "or", variables[Code == var2, "Element"]) + } else { # if variables1 and 2 have different elements and yeartimes + ylab <- paste(yeartime.names[which(yeartimes == yeartime1)], variables[Code == var1, "Element"], "or", variables[Code == var2, "Element"]) + } + plot(0, col = "white", xlim = c(1900, 2100), ylim = range(if (yfit) visibledata else alldata, na.rm = TRUE), xaxs = "i", xaxt = "n", tck = 0, xlab = "", ylab = ylab) + axis(1, at = seq(1850, 2100, 25), labels = seq(1850, 2100, 25), tck = 0) + + num <- 1 + for (num in nums) { + yeartime <- get(paste("yeartime", num, sep = "")) + element <- get(paste("element", num, sep = "")) + var <- get(paste("var", num, sep = "")) + + # function for plotting time series for gcms or compiled ensemble + plot.ensemble <- function(x) { + # scenario <- scenarios.selected[1] + temp.historical <- x[is.na(SSP), c("PERIOD", "RUN", var), with = FALSE] + x.historical <- as.numeric(temp.historical[, .(min = min(get(var))), by = PERIOD][["PERIOD"]]) + ensmin.historical <- temp.historical[RUN != "ensembleMean", .(min = min(get(var), na.rm = TRUE)), by = PERIOD][["min"]] + ensmax.historical <- temp.historical[RUN != "ensembleMean", .(max = max(get(var), na.rm = TRUE)), by = PERIOD][["max"]] + ensmean.historical <- temp.historical[RUN == "ensembleMean", .(mean = mean(get(var), na.rm = TRUE)), by = PERIOD][["mean"]] # calculate mean using the single-model ensembleMean, to ensure that the mean isn't biased towards models with more runs + + for (scenario in scenarios.selected[order(c(1, 4, 5, 3, 2)[which(scenarios %in% scenarios.selected)])][-1]) { + temp <- x[SSP == scenario, c("PERIOD", "RUN", var), with = FALSE] + x.temp <- as.numeric(temp[, .(min = min(get(var))), by = PERIOD][["PERIOD"]]) + ensmin.temp <- temp[RUN != "ensembleMean", .(min = min(get(var), na.rm = TRUE)), by = PERIOD][["min"]] + ensmax.temp <- temp[RUN != "ensembleMean", .(max = max(get(var), na.rm = TRUE)), by = PERIOD][["max"]] + ensmean.temp <- temp[RUN == "ensembleMean", .(mean = mean(get(var), na.rm = TRUE)), by = PERIOD][["mean"]] # calculate mean using the single-model ensembleMean, to ensure that the mean isn't biased towards models with more runs + assign(paste("x", scenario, sep = "."), c(x.historical[length(x.historical)], x.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. + assign(paste("ensmin", scenario, sep = "."), c(ensmin.historical[length(ensmin.historical)], ensmin.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. + assign(paste("ensmax", scenario, sep = "."), c(ensmax.historical[length(ensmax.historical)], ensmax.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. + assign(paste("ensmean", scenario, sep = "."), c(ensmean.historical[length(ensmean.historical)], ensmean.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. + } + + if (showrange) { + if (simplify == FALSE) { + for (scenario in scenarios.selected[order(c(1, 4, 5, 3, 2)[which(scenarios %in% scenarios.selected)])]) { + x3 <- get(paste("x", scenario, sep = ".")) + polygon(c(x3, rev(x3)), c(get(paste("ensmin", scenario, sep = ".")), rev(get(paste("ensmax", scenario, sep = ".")))), col = alpha(colSelect(scenario, gcms), 0.35), border = colSelect(scenario, gcms)) + } + } else { + scenarios.select <- scenarios.selected[order(c(1, 4, 5, 3, 2)[which(scenarios %in% scenarios.selected)])][-1] + for (scenario in scenarios.select) { + if (scenario == scenarios.select[1]) { # we need to run spline through the historical/projected transition + x4 <- c(x.historical, get(paste("x", scenario, sep = "."))[-1]) + y.ensmin <- c(ensmin.historical, get(paste("ensmin", scenario, sep = "."))[-1]) + y.ensmax <- c(ensmax.historical, get(paste("ensmax", scenario, sep = "."))[-1]) + s.ensmin <- smooth.spline(x4, y.ensmin, df = 8) + s.ensmax <- smooth.spline(x4, y.ensmax, df = 8) + subset.hist <- which(x4 %in% x.historical) + subset.proj <- which(x4 %in% get(paste("x", scenario, sep = "."))) + polygon(c(s.ensmin$x[subset.hist], rev(s.ensmax$x[subset.hist])), c(s.ensmin$y[subset.hist], rev(s.ensmax$y[subset.hist])), col = alpha(if (pal == "gcms") colSelect(scenario, gcms) else pal.scenario[which(scenarios == "historical")], 0.35), border = NA) + polygon(c(s.ensmin$x[subset.proj], rev(s.ensmax$x[subset.proj])), c(s.ensmin$y[subset.proj], rev(s.ensmax$y[subset.proj])), col = alpha(colSelect(scenario, gcms), 0.35), border = NA) + } else { # this second routine uses interpolation splines so that the starting point for all scenarios is the same + x5 <- c(x.historical, get(paste("x", scenario, sep = "."))[-1]) + y.ensmin2 <- c(ensmin.historical, get(paste("ensmin", scenario, sep = "."))[-1]) + y.ensmax2 <- c(ensmax.historical, get(paste("ensmax", scenario, sep = "."))[-1]) + s.ensmin2 <- smooth.spline(x5, y.ensmin2, df = 8) + s.ensmax2 <- smooth.spline(x5, y.ensmax2, df = 8) + knots.hist <- which(x5 %in% c(seq(1860, 2000, 20), 2014)) + knots.proj <- which(x5 %in% c(seq(2030, 2090, 20), 2100)) + s.ensmin3 <- stinterp(x5[c(knots.hist, knots.proj)], c(s.ensmin$y[knots.hist], s.ensmin2$y[knots.proj]), x5) + s.ensmax3 <- stinterp(x5[c(knots.hist, knots.proj)], c(s.ensmax$y[knots.hist], s.ensmax2$y[knots.proj]), x5) + polygon(c(s.ensmin3$x[subset.proj], rev(s.ensmax3$x[subset.proj])), c(s.ensmin3$y[subset.proj], rev(s.ensmax3$y[subset.proj])), col = alpha(colSelect(scenario, gcms), 0.35), border = NA) } } } - - if(refline){ - ref.temp <- mean(ensmean.historical[which(x.historical%in%1961:1990)]) - lines(1961:1990, rep(ref.temp, 30), lwd=2) - lines(c(1990,2100), rep(ref.temp, 2), lty=2) + } + + if (refline) { + ref.temp <- mean(ensmean.historical[which(x.historical %in% 1961:1990)]) + lines(1961:1990, rep(ref.temp, 30), lwd = 2) + lines(c(1990, 2100), rep(ref.temp, 2), lty = 2) + } + + # overlay the ensemble mean lines on top of all polygons + for (scenario in scenarios.selected[order(c(1, 4, 5, 3, 2)[which(scenarios %in% scenarios.selected)])]) { + # calculate a spline through the time series (used for plotting and the text warming value) + if (scenario == "historical") { # need to run spline through the historical/projected transition + x4 <- c(x.historical, get(paste("x", scenarios.selected[2], sep = "."))) + y4 <- c(ensmean.historical, get(paste("ensmean", scenarios.selected[2], sep = "."))) + } else { + x4 <- c(x.historical, get(paste("x", scenario, sep = "."))) + y4 <- c(ensmean.historical, get(paste("ensmean", scenario, sep = "."))) } - - # overlay the ensemble mean lines on top of all polygons - for(scenario in scenarios.selected[order(c(1,4,5,3,2)[which(scenarios%in%scenarios.selected)])]){ - - # calculate a spline through the time series (used for plotting and the text warming value) - if(scenario=="historical"){ # need to run spline through the historical/projected transition - x4 <- c(x.historical, get(paste("x", scenarios.selected[2], sep="."))) - y4 <- c(ensmean.historical, get(paste("ensmean", scenarios.selected[2], sep="."))) + s4 <- smooth.spline(x4, y4, df = 10) + subset <- which(x4 %in% get(paste("x", scenario, sep = "."))) + + # plot the ensemble mean + if (showmean) { + if (simplify) { + lines(x = s4$x[subset], y = s4$y[subset], col = colSelect(scenario, gcms), lwd = 2) } else { - x4 <- c(x.historical, get(paste("x", scenario, sep="."))) - y4 <- c(ensmean.historical, get(paste("ensmean", scenario, sep="."))) - } - s4 <- smooth.spline(x4,y4, df=10) - subset <- which(x4%in%get(paste("x", scenario, sep="."))) - - # plot the ensemble mean - if(showmean){ - if(simplify){ - lines(x=s4$x[subset], y=s4$y[subset], col=colSelect(scenario, gcms), lwd=2) - } else lines(x=get(paste("x", scenario, sep=".")), y=get(paste("ensmean", scenario, sep=".")), col=colSelect(scenario, gcms), lwd=2) + lines(x = get(paste("x", scenario, sep = ".")), y = get(paste("ensmean", scenario, sep = ".")), col = colSelect(scenario, gcms), lwd = 2) } - - # end labels - if(is.null(endlabel)==FALSE){ - if(scenario != "historical"){ - par(xpd=TRUE) - baseline <- mean(ensmean.historical[which(x.historical%in%1961:1990)]) - projected <- s4$y[length(s4$y)] - if(endlabel=="change"){ - if(element%in%c("PPT", "PAS", "CMD", "MAP", "MSP", "DDsub18", "DD18", "DDsub0", "DD5", "Eref")){ - change <- round(projected/baseline-1,2) - if(is.na(change)==FALSE) text(2099,projected, if(change>0) paste("+",change*100,"%", sep="") else paste(change*100,"%", sep=""), col=colSelect(scenario, gcms), pos=4, font=2, cex=1) - } else if(element%in%c("Tave", "Tmin", "Tmax", "MCMT", "MWMT", "EXT", "EMT", "MAT")){ - change <- round(projected-baseline,1) - if(is.na(change)==FALSE) text(2099,projected, if(change>0) bquote("+" * .(change) * degree * C) else bquote(.(change) * degree * C), col=colSelect(scenario, gcms), pos=4, font=2, cex=1) - } else if(element%in%c("RH")){ - change <- round(projected-baseline,1) - if(is.na(change)==FALSE) text(2099,projected, if(change>0) paste("+",change,"%", sep="") else paste(change,"%", sep=""), col=colSelect(scenario, gcms), pos=4, font=2, cex=1) - } else { - change <- round(projected-baseline,1) - if(is.na(change)==FALSE) text(2099,projected, if(change>0) paste("+",change, sep="") else paste(change, sep=""), col=colSelect(scenario, gcms), pos=4, font=2, cex=1) - } - } - if(endlabel=="gcms"){ - text(2099,projected, if(compile) "ensemble" else gcms, col=colSelect(scenario, gcms), pos=4, font=1, cex=1) + } + + # end labels + if (is.null(endlabel) == FALSE) { + if (scenario != "historical") { + par(xpd = TRUE) + baseline <- mean(ensmean.historical[which(x.historical %in% 1961:1990)]) + projected <- s4$y[length(s4$y)] + if (endlabel == "change") { + if (element %in% c("PPT", "PAS", "CMD", "MAP", "MSP", "DDsub18", "DD18", "DDsub0", "DD5", "Eref")) { + change <- round(projected / baseline - 1, 2) + if (is.na(change) == FALSE) text(2099, projected, if (change > 0) paste("+", change * 100, "%", sep = "") else paste(change * 100, "%", sep = ""), col = colSelect(scenario, gcms), pos = 4, font = 2, cex = 1) + } else if (element %in% c("Tave", "Tmin", "Tmax", "MCMT", "MWMT", "EXT", "EMT", "MAT")) { + change <- round(projected - baseline, 1) + if (is.na(change) == FALSE) text(2099, projected, if (change > 0) bquote("+" * .(change) * degree * C) else bquote(.(change) * degree * C), col = colSelect(scenario, gcms), pos = 4, font = 2, cex = 1) + } else if (element %in% c("RH")) { + change <- round(projected - baseline, 1) + if (is.na(change) == FALSE) text(2099, projected, if (change > 0) paste("+", change, "%", sep = "") else paste(change, "%", sep = ""), col = colSelect(scenario, gcms), pos = 4, font = 2, cex = 1) + } else { + change <- round(projected - baseline, 1) + if (is.na(change) == FALSE) text(2099, projected, if (change > 0) paste("+", change, sep = "") else paste(change, sep = ""), col = colSelect(scenario, gcms), pos = 4, font = 2, cex = 1) } - par(xpd=FALSE) } + if (endlabel == "gcms") { + text(2099, projected, if (compile) "ensemble" else gcms, col = colSelect(scenario, gcms), pos = 4, font = 1, cex = 1) + } + par(xpd = FALSE) } } - - # Text to identify the time of year - if(!is.null(var2)){ #if there is no second var - if(element1==element2){ - label <- yeartime.names[which(yeartimes==yeartime)] - } else { - label <- paste(yeartime.names[which(yeartimes==yeartime)], element) - } - temp <- get("ensmax.historical") - text(1925,mean(temp[10:40]), label, col="black", pos=3, font=2, cex=1) + } + + # Text to identify the time of year + if (!is.null(var2)) { # if there is no second var + if (element1 == element2) { + label <- yeartime.names[which(yeartimes == yeartime)] + } else { + label <- paste(yeartime.names[which(yeartimes == yeartime)], element) } + temp <- get("ensmax.historical") + text(1925, mean(temp[10:40]), label, col = "black", pos = 3, font = 2, cex = 1) } - - if(compile){ #this plots a single envelope for the ensemble as a whole - temp.data <- X[GCM%in%gcms, c("PERIOD", "SSP", "RUN", var), with=FALSE] - plot.ensemble(temp.data) - - } else for(gcm in gcms){ #this plots of individual GCM ensembles. - temp.data <- X[GCM==gcm, c("PERIOD", "SSP", "RUN", var), with=FALSE] + } + + if (compile) { # this plots a single envelope for the ensemble as a whole + temp.data <- X[GCM %in% gcms, c("PERIOD", "SSP", "RUN", var), with = FALSE] + plot.ensemble(temp.data) + } else { + for (gcm in gcms) { # this plots of individual GCM ensembles. + temp.data <- X[GCM == gcm, c("PERIOD", "SSP", "RUN", var), with = FALSE] plot.ensemble(temp.data) - + print(gcms) } - - # overlay the 5-year lines on top of all polygons - if(yearlines){ - for(n in seq(1905, 2095, 5)){ - lines(c(n, n), c(-9999, 9999), col="grey", lty=2) - } + } + + # overlay the 5-year lines on top of all polygons + if (yearlines) { + for (n in seq(1905, 2095, 5)) { + lines(c(n, n), c(-9999, 9999), col = "grey", lty = 2) } - - if(showObserved){ - # add in observations - obs.colors <- c("black", "blue", "red") - obs.options <- c("climatena", "cru.gpcc") ##, "era5" - for(obs.dataset in obs_ts_dataset){ #TODO update this code block once i know how the datasets are identified in the climr output - obs.color <- obs.colors[which(obs.options==obs.dataset)] - x.obs <- as.numeric(X[DATASET==obs.dataset & PERIOD%in%1900:2100, "PERIOD"][[1]]) - y.obs <- X[DATASET==obs.dataset & PERIOD%in%1900:2100, get(var)] - recent.obs <- mean(y.obs[which(x.obs%in%2014:2023)], na.rm=TRUE) - baseline.obs <- mean(y.obs[which(x.obs%in%1961:1990)], na.rm=TRUE) - end <- max(which(!is.na(y.obs))) - lines(x.obs[which(x.obs<1951)], y.obs[which(x.obs<1951)], lwd=3, lty=3, col=obs.color) - lines(x.obs[which(x.obs>1949)], y.obs[which(x.obs>1949)], lwd=4, col=obs.color) - if(yearmarkers){ - points(x.obs[which(x.obs>1949)], y.obs[which(x.obs>1949)], pch=21, bg="white", col=obs.color, cex=0.4) - points(x.obs[which(x.obs>1949)[seq(1,999,5)]], y.obs[which(x.obs>1949)[seq(1,999,5)]], pch=21, bg="white", col=obs.color, cex=0.7) - } - if(label.endyear){ - points(x.obs[end], y.obs[end], pch=16, cex=1, col=obs.color) - text(x.obs[end], y.obs[end], x.obs[end], pos= 4, offset = 0.25, col=obs.color, cex=1,) - } - if(refline.obs){ - lines(1961:1990, rep(baseline.obs, 30), lwd=1, col=obs.color) - lines(c(1990,2100), rep(baseline.obs, 2), lty=2, col=obs.color) - } + } + + if (showObserved) { + # add in observations + obs.colors <- c("black", "blue", "red") + obs.options <- c("climatena", "cru.gpcc") ## , "era5" + for (obs.dataset in obs_ts_dataset) { # TODO update this code block once i know how the datasets are identified in the climr output + obs.color <- obs.colors[which(obs.options == obs.dataset)] + x.obs <- as.numeric(X[DATASET == obs.dataset & PERIOD %in% 1900:2100, "PERIOD"][[1]]) + y.obs <- X[DATASET == obs.dataset & PERIOD %in% 1900:2100, get(var)] + recent.obs <- mean(y.obs[which(x.obs %in% 2014:2023)], na.rm = TRUE) + baseline.obs <- mean(y.obs[which(x.obs %in% 1961:1990)], na.rm = TRUE) + end <- max(which(!is.na(y.obs))) + lines(x.obs[which(x.obs < 1951)], y.obs[which(x.obs < 1951)], lwd = 3, lty = 3, col = obs.color) + lines(x.obs[which(x.obs > 1949)], y.obs[which(x.obs > 1949)], lwd = 4, col = obs.color) + if (yearmarkers) { + points(x.obs[which(x.obs > 1949)], y.obs[which(x.obs > 1949)], pch = 21, bg = "white", col = obs.color, cex = 0.4) + points(x.obs[which(x.obs > 1949)[seq(1, 999, 5)]], y.obs[which(x.obs > 1949)[seq(1, 999, 5)]], pch = 21, bg = "white", col = obs.color, cex = 0.7) + } + if (label.endyear) { + points(x.obs[end], y.obs[end], pch = 16, cex = 1, col = obs.color) + text(x.obs[end], y.obs[end], x.obs[end], pos = 4, offset = 0.25, col = obs.color, cex = 1, ) + } + if (refline.obs) { + lines(1961:1990, rep(baseline.obs, 30), lwd = 1, col = obs.color) + lines(c(1990, 2100), rep(baseline.obs, 2), lty = 2, col = obs.color) } } - print(num) - } - - if(showObserved){ - # Sources legend - a <- if("climatena"%in%obs_ts_dataset) 1 else NA - b <- if("cru.gpcc"%in%obs_ts_dataset) 2 else NA - c <- if("era5"%in%obs_ts_dataset) 3 else NA - d <- if(length(gcms>0)) 4 else NA - s <- !is.na(c(a,b,c,d)) - legend.GCM <- if(length(gcms)>1) paste("Simulated (", length(gcms), " GCMs)", sep="") else paste("Simulated (", gcms, ")", sep="") - legend(legend_pos, title = "", legend=c("Observed (ClimateNA)", "Observed (CRU/GPCC)", "Observed (ERA5)", legend.GCM)[s], bty="n", - lty=c(1,1,1,1)[s], - col=c(obs.colors, "gray")[s], - lwd=c(4,4,4,2)[s], - pch=c(NA,NA,NA,NA)[s], - pt.bg = c(NA, NA,NA,NA)[s], - pt.cex=c(NA,NA,NA,NA)[s]) } - - #Scenario legend - if(pal=="gcms"){ - s <- which(list_gcms()%in%gcms) - legend(c("top", "bottom")[if(length(grep("top", legend_pos))==1) 1 else 2], title = "GCMs", legend=gcms, bty="n", - col=pal.gcms[s], pch=22, pt.bg = alpha(pal.gcms[s], 0.35), pt.cex=2) - } else { - s <- rev(which(scenarios[-1]%in%scenarios.selected)) - legend(c("top", "bottom")[if(length(grep("top", legend_pos))==1) 1 else 2], title = "Scenarios", legend=c("Historical", scenario.names[-1][s]), bty="n", - lty=c(NA,NA,NA,NA,NA)[c(1,s+1)], col=pal.scenario[c(1,s+1)], lwd=c(NA,NA,NA,NA,NA)[c(1,s+1)], pch=c(22,22,22,22,22)[c(1,s+1)], pt.bg = alpha(pal.scenario[c(1,s+1)], 0.35), pt.cex=c(2,2,2,2,2)[c(1,s+1)]) - } - - # mtext(paste(" Created using climr (https://bcgov.github.io/climr/)"), side=1, line=1.5, adj=0.0, font=1, cex=1.1, col="gray") - - box() - + print(num) + } + + if (showObserved) { + # Sources legend + a <- if ("climatena" %in% obs_ts_dataset) 1 else NA + b <- if ("cru.gpcc" %in% obs_ts_dataset) 2 else NA + c <- if ("era5" %in% obs_ts_dataset) 3 else NA + d <- if (length(gcms > 0)) 4 else NA + s <- !is.na(c(a, b, c, d)) + legend.GCM <- if (length(gcms) > 1) paste("Simulated (", length(gcms), " GCMs)", sep = "") else paste("Simulated (", gcms, ")", sep = "") + legend(legend_pos, + title = "", legend = c("Observed (ClimateNA)", "Observed (CRU/GPCC)", "Observed (ERA5)", legend.GCM)[s], bty = "n", + lty = c(1, 1, 1, 1)[s], + col = c(obs.colors, "gray")[s], + lwd = c(4, 4, 4, 2)[s], + pch = c(NA, NA, NA, NA)[s], + pt.bg = c(NA, NA, NA, NA)[s], + pt.cex = c(NA, NA, NA, NA)[s] + ) } + + # Scenario legend + if (pal == "gcms") { + s <- which(list_gcms() %in% gcms) + legend(c("top", "bottom")[if (length(grep("top", legend_pos)) == 1) 1 else 2], + title = "GCMs", legend = gcms, bty = "n", + col = pal.gcms[s], pch = 22, pt.bg = alpha(pal.gcms[s], 0.35), pt.cex = 2 + ) + } else { + s <- rev(which(scenarios[-1] %in% scenarios.selected)) + legend(c("top", "bottom")[if (length(grep("top", legend_pos)) == 1) 1 else 2], + title = "Scenarios", legend = c("Historical", scenario.names[-1][s]), bty = "n", + lty = c(NA, NA, NA, NA, NA)[c(1, s + 1)], col = pal.scenario[c(1, s + 1)], lwd = c(NA, NA, NA, NA, NA)[c(1, s + 1)], pch = c(22, 22, 22, 22, 22)[c(1, s + 1)], pt.bg = alpha(pal.scenario[c(1, s + 1)], 0.35), pt.cex = c(2, 2, 2, 2, 2)[c(1, s + 1)] + ) + } + + # mtext(paste(" Created using climr (https://bcgov.github.io/climr/)"), side=1, line=1.5, adj=0.0, font=1, cex=1.1, col="gray") + + box() + } } diff --git a/man/plot_timeSeries.Rd b/man/plot_timeSeries.Rd index 83b4d13..299371e 100644 --- a/man/plot_timeSeries.Rd +++ b/man/plot_timeSeries.Rd @@ -91,8 +91,8 @@ to indicate the global climate model.} \item{yearlines}{logical. Add vertical lines on every fifth year as a visual reference} -\item{legend_pos}{character. Position of the legend. Viable options are c("bottomright", -"bottomleft", "topleft", and "topright").} +\item{legend_pos}{character. Position of the legend. Viable options are "bottomright", +"bottomleft", "topleft", and "topright".} } \value{ NULL. Draws a plot in the active graphics device. @@ -123,38 +123,45 @@ we recommend users dedicate some time to run this function in background. Once t don't need to be downloaded again. } \examples{ -if(FALSE){ -# data frame of arbitrary points -my_points <- data.frame(lon = c(-127.7300,-127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) - -# generate the input data -my_data <- plot_timeSeries_input(my_points) - -# use the input to create a plot -plot_timeSeries(my_data, var1 = "Tmin_sm") - -# compare observational time series -plot_timeSeries(my_data, var1 = "Tmin_sm", obs_ts_dataset = c("cru.gpcc", "climatena")) - -# compare mean daily minimum and maximum temperatures -plot_timeSeries(my_data, var1 = "Tmin_sm", var2 = "Tmax_sm") - -# compare summer and winter temperatures (without simplifying the ensemble range) -plot_timeSeries(my_data, var1 = "Tmax_sm", var2 = "Tmax_wt", simplify = FALSE) - -# compare global climate models -plot_timeSeries(my_data, gcms = list_gcms()[c(7,13)], pal = "gcms", ssps = list_ssps()[2], showmean = FALSE, compile = FALSE, simplify = FALSE, endlabel = "gcms", mar=c(3,3,0.1,6), showObserved = FALSE) - -# export plot to a temporary directory, including a title -figDir <- tempdir() -png( - filename = file.path(figDir, "plot_test.png"), type = "cairo", units = "in", - width = 6, height = 5, pointsize = 10, res = 300 -) -plot_timeSeries(my_data, var1 = "Tmin_sm", mar=c(3,3,2,4)) -title("Historical and projected summer night-time warming in the Bulkley Valley, BC") -dev.off() +if (FALSE) { + # data frame of arbitrary points + my_points <- data.frame( + lon = c(-127.7300, -127.7500), + lat = c(55.34114, 55.25), + elev = c(711, 500), + id = 1:2 + ) + + # generate the input data + my_data <- plot_timeSeries_input(my_points) + + # use the input to create a plot + plot_timeSeries(my_data, var1 = "Tmin_sm") + + # compare observational time series + plot_timeSeries(my_data, var1 = "Tmin_sm", obs_ts_dataset = c("cru.gpcc", "climatena")) + + # compare mean daily minimum and maximum temperatures + plot_timeSeries(my_data, var1 = "Tmin_sm", var2 = "Tmax_sm") + + # compare summer and winter temperatures (without simplifying the ensemble range) + plot_timeSeries(my_data, var1 = "Tmax_sm", var2 = "Tmax_wt", simplify = FALSE) + + # compare global climate models + plot_timeSeries(my_data, gcms = list_gcms()[c(7, 13)], pal = "gcms", + ssps = list_ssps()[2], showmean = FALSE, compile = FALSE, + simplify = FALSE, endlabel = "gcms", mar = c(3, 3, 0.1, 6), + showObserved = FALSE) + + # export plot to a temporary directory, including a title + figDir <- tempdir() + png( + filename = file.path(figDir, "plot_test.png"), type = "cairo", units = "in", + width = 6, height = 5, pointsize = 10, res = 300 + ) + plot_timeSeries(my_data, var1 = "Tmin_sm", mar = c(3, 3, 2, 4)) + title("Historical and projected summer night-time warming in the Bulkley Valley, BC") + dev.off() } - }