From dca72141a5c898d51081216645ef90a148c3a23d Mon Sep 17 00:00:00 2001 From: cmahony Date: Sun, 11 Feb 2024 14:30:27 -0800 Subject: [PATCH 01/17] first draft of plot_bivariate() --- R/plot_bivariate.R | 131 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 131 insertions(+) create mode 100644 R/plot_bivariate.R diff --git a/R/plot_bivariate.R b/R/plot_bivariate.R new file mode 100644 index 00000000..d53cf1e9 --- /dev/null +++ b/R/plot_bivariate.R @@ -0,0 +1,131 @@ +#' Bivariate plots of 21st century climate change at user-selected locations +#' +#' @description +#' +#' +#' @details +#' +#' +#' @template xyz +#' +#' @return `data.frame` of downscaled climate variables for each location. +#' All outputs are returned in one table. +#' +#' @importFrom data.table getDTthreads setDTthreads rbindlist setkey +#' @importFrom stinepack stinterp TODO: add this package to dependencies +#' +#' @examples { +#' library(data.table) +#' +#' } +#' @export + +Test_Locations_VanIsl <- read.csv("data-raw/Test_Locations_VanIsl.csv") +xyz <- Test_Locations_VanIsl +names(xyz) <- c("lat", "lon", "elev") +xyz <- data.frame(Test_Locations_VanIsl, id=1:dim(Test_Locations_VanIsl)[1]) +names(xyz) <- c("lat", "lon", "elev", "id") +plot_bivariate <- function( + xyz, + xvar = "Tave_sm", + yvar = "PPT_sm", + percent_x = NULL, + percent_y = NULL, + period_focal = list_gcm_period()[1], + gcm_models = list_gcm()[c(1,4,5,6,7,10,11,12)], ## TODO: check that i am passing this to climr_downscale properly + ssp = list_ssp()[2], ## TODO: check that i am passing this to climr_downscale properly + legend_pos = "topright", + show_runs = T, + show_ensMean = T, + show_observed = T, + show_trajectories = T, + interactive = F, + ... +) { + + variables_lookup <- read.csv("data-raw/variables_lookup.csv") + + # variable types for default scaling (percent or absolute) + xvar_type <- variables_lookup$Type[which(variables_lookup$Code==xvar)] + yvar_type <- variables_lookup$Type[which(variables_lookup$Code==yvar)] + + colors = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)][-1] + set.seed(2) + ColScheme <- sample(colors,length(gcm_models)) # TODO select better colors + + # generate the climate data + data <- climr_downscale(xyz, + which_normal = "auto", + historic_period = list_historic()[1], + gcm_models = list_gcm()[c(1,4,5,6,7,10,11,12)], ## TODO: this is wrong. need to recieve this as passed from plot_bivariate() + ssp = list_ssp()[2], ## TODO: this is wrong. need to receive this as passed from plot_bivariate() + gcm_period = list_gcm_period(), + max_run = 10, + vars = c(xvar, yvar) + ) + + # convert absolute values to anomalies + data[, xanom := if(xvar_type=="ratio") (get(xvar)/get(xvar)[1]-1) else (get(xvar) - get(xvar)[1]), by = id] + data[, yanom := if(yvar_type=="ratio") (get(yvar)/get(yvar)[1]-1) else (get(yvar) - get(yvar)[1]), by = id] + + # collapse the points down to a mean anomaly + data <- data[, .(xanom = mean(xanom), yanom = mean(yanom)), by = .(GCM, SSP, RUN, PERIOD)] + + # initiate the plot + par(mar=c(3,4,0,1), mgp=c(1.25, 0.25,0), cex=1.5) + plot(data_mean$xanom,data_mean$yanom,col="white", tck=0, xaxt="n", yaxt="n", ylab="", + xlab=paste("Change in", variables_lookup$Variable[which(variables_lookup$Code==xvar)]) + ) + par(mgp=c(2.5,0.25, 0)) + title(ylab=paste("Change in", variables_lookup$Variable[which(variables_lookup$Code==yvar)])) + lines(c(0,0), c(-99,99), lty=2, col="gray") + lines(c(-99,99), c(0,0), lty=2, col="gray") + + # ensemble mean for the selected period + ensMean <- data[!is.na(GCM) & RUN == "ensembleMean" & PERIOD == period_focal, .(xanom = mean(xanom), yanom = mean(yanom)), ] + points(ensMean$xanom,ensMean$yanom, pch=43, col="gray", cex=3) + + # observed climate + obs <- data[is.na(GCM) & PERIOD == period_focal] + points(obs$xanom ,obs$yanom, pch=22, bg="gray", cex=2.5) + # text(x1,y1, "1991-2019", cex=1.15, font=2, pos=4, col="gray", offset=0.9) + + # GCM projections + gcm=list_gcm()[1] + for(gcm in gcm_models){ + i=which(gcm_models==gcm) + x <- c(0, data[GCM==gcm & RUN == "ensembleMean", xanom]) + y <- c(0, data[GCM==gcm & RUN == "ensembleMean", yanom]) + x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] + y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] + if(length(unique(sign(diff(x))))==1){ + x3 <- if(unique(sign(diff(x)))==-1) rev(x) else x + y3 <- if(unique(sign(diff(x)))==-1) rev(y) else y + s <- stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data$xanom))/500)) # way better than interpSpline, not prone to oscillations + lines(s, col=ColScheme[i], lwd=2) + } else lines(x, y, col=ColScheme[i], lwd=2) + j=which(list_gcm_period()==period_focal)+1 + points(x,y, pch=16, col=ColScheme[i], cex=0.5) + points(x[j],y[j], pch=21, bg=ColScheme[i], cex=2.5) + points(x.runs,y.runs, pch=21, bg=ColScheme[i], cex=1) + text(x[j],y[j], substr(gcm_models, 1, 2)[i], cex=0.5, font=2) + } + + # axis labels + axis(1, at=pretty(data_mean$xanom), labels=if(xvar_type=="ratio") paste(pretty(data_mean$xanom)*100, "%", sep="") else pretty(data_mean$xanom), tck=0) + axis(2, at=pretty(data_mean$yanom), labels=if(yvar_type=="ratio") paste(pretty(data_mean$yanom)*100, "%", sep="") else pretty(data_mean$yanom), las=2, tck=0) + + # Legend TODO: modify for user selection of plot elements. + legend(legend_pos, + legend = c("Observed historical climate", "individual GCM simulation", "GCM mean", "GCM mean trajectory", "Ensemble mean"), + pch = c(22, 21, 21, 16, 43), + pt.cex = c(2.5, 1, 2.5, 0.5, 3), + pt.bg = c("gray", "gray", "gray", NA, NA), + col = c(1,1,1,"gray","gray"), + lty = c(NA, NA, NA, 1, NA), + lwd = c(NA, NA, NA, 2, NA), + bty = "n", cex=0.8 + ) + +} + From 00b20acad97c26b0ab5d60ede2e7f08dffddb43e Mon Sep 17 00:00:00 2001 From: cmahony Date: Sun, 11 Feb 2024 22:16:28 -0800 Subject: [PATCH 02/17] refine base plot and add documentation --- R/plot_bivariate.R | 146 +++++++++++++++++++++++++-------------------- 1 file changed, 80 insertions(+), 66 deletions(-) diff --git a/R/plot_bivariate.R b/R/plot_bivariate.R index d53cf1e9..79134d8e 100644 --- a/R/plot_bivariate.R +++ b/R/plot_bivariate.R @@ -1,36 +1,44 @@ -#' Bivariate plots of 21st century climate change at user-selected locations +#' Bivariate climate change plots #' #' @description -#' +#' Bivariate plots of 21st century climate change for user-selected locations and climate variables. #' #' @details +#' The input table `xyz` can be a single location or multiple locations. If multiple locations, the plot provides the mean of the anomalies for these locations. #' +#' The climate change trajectories provided by `show_trajectories` are connected with an interpolation spline when the x variable is monotonic; otherwise the trajectory points are connected by straight lines. #' #' @template xyz +#' @param xvar character. x-axis variable. options are `list_variables()`. +#' @param yvar character. y-axis variable. options are `list_variables()`. +#' @param period_focal character. The 20-year period for which to plot the ensemble detail. options are `list_gcm_period()`. +#' @param legend_pos character. Position of the legend. Options are c("bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right", "center"). +#' @param show_runs logical. If TRUE, the individual runs of the model are plotted (for `period_focal` only) in addition to the single-model ensemble mean. +#' @param show_ensMean logical. If TRUE, the multi-model ensemble mean is plotted (for `period_focal` only). +#' @param show_observed logical. If TRUE, the 2001-2020 observed climate is plotted. +#' @param show_trajectories logical. If TRUE, the values of the single-model ensemble mean are plotted for all 20-year periods in `list_gcm_period()`, connected by an interpolation spline. +#' @param interactive logical. If TRUE, an interactive plot is generated using plotly(). If FALSE, a plot is generated using base graphics. +#' +#' @return TODO #' -#' @return `data.frame` of downscaled climate variables for each location. -#' All outputs are returned in one table. -#' -#' @importFrom data.table getDTthreads setDTthreads rbindlist setkey +#' @importFrom data.table TODO: find out what i need to do to get this command right. #' @importFrom stinepack stinterp TODO: add this package to dependencies #' #' @examples { #' library(data.table) -#' +#' TODO: add example usage #' } #' @export -Test_Locations_VanIsl <- read.csv("data-raw/Test_Locations_VanIsl.csv") -xyz <- Test_Locations_VanIsl -names(xyz) <- c("lat", "lon", "elev") -xyz <- data.frame(Test_Locations_VanIsl, id=1:dim(Test_Locations_VanIsl)[1]) -names(xyz) <- c("lat", "lon", "elev", "id") +# Test_Locations_VanIsl <- read.csv("data-raw/Test_Locations_VanIsl.csv") +# xyz <- data.frame(Test_Locations_VanIsl, id=1:dim(Test_Locations_VanIsl)[1]) +# names(xyz) <- c("lat", "lon", "elev", "id") plot_bivariate <- function( xyz, xvar = "Tave_sm", yvar = "PPT_sm", - percent_x = NULL, - percent_y = NULL, + # percent_x = NULL, TODO: set up an override for ratio variables being expressed as percent anomalies + # percent_y = NULL, TODO: set up an override for ratio variables being expressed as percent anomalies period_focal = list_gcm_period()[1], gcm_models = list_gcm()[c(1,4,5,6,7,10,11,12)], ## TODO: check that i am passing this to climr_downscale properly ssp = list_ssp()[2], ## TODO: check that i am passing this to climr_downscale properly @@ -39,7 +47,7 @@ plot_bivariate <- function( show_ensMean = T, show_observed = T, show_trajectories = T, - interactive = F, + interactive = F, ## TODO: add a plotly version of the plot ... ) { @@ -70,62 +78,68 @@ plot_bivariate <- function( # collapse the points down to a mean anomaly data <- data[, .(xanom = mean(xanom), yanom = mean(yanom)), by = .(GCM, SSP, RUN, PERIOD)] - - # initiate the plot - par(mar=c(3,4,0,1), mgp=c(1.25, 0.25,0), cex=1.5) - plot(data_mean$xanom,data_mean$yanom,col="white", tck=0, xaxt="n", yaxt="n", ylab="", - xlab=paste("Change in", variables_lookup$Variable[which(variables_lookup$Code==xvar)]) - ) - par(mgp=c(2.5,0.25, 0)) - title(ylab=paste("Change in", variables_lookup$Variable[which(variables_lookup$Code==yvar)])) - lines(c(0,0), c(-99,99), lty=2, col="gray") - lines(c(-99,99), c(0,0), lty=2, col="gray") - # ensemble mean for the selected period ensMean <- data[!is.na(GCM) & RUN == "ensembleMean" & PERIOD == period_focal, .(xanom = mean(xanom), yanom = mean(yanom)), ] - points(ensMean$xanom,ensMean$yanom, pch=43, col="gray", cex=3) - # observed climate obs <- data[is.na(GCM) & PERIOD == period_focal] - points(obs$xanom ,obs$yanom, pch=22, bg="gray", cex=2.5) - # text(x1,y1, "1991-2019", cex=1.15, font=2, pos=4, col="gray", offset=0.9) - # GCM projections - gcm=list_gcm()[1] - for(gcm in gcm_models){ - i=which(gcm_models==gcm) - x <- c(0, data[GCM==gcm & RUN == "ensembleMean", xanom]) - y <- c(0, data[GCM==gcm & RUN == "ensembleMean", yanom]) - x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] - y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] - if(length(unique(sign(diff(x))))==1){ - x3 <- if(unique(sign(diff(x)))==-1) rev(x) else x - y3 <- if(unique(sign(diff(x)))==-1) rev(y) else y - s <- stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data$xanom))/500)) # way better than interpSpline, not prone to oscillations - lines(s, col=ColScheme[i], lwd=2) - } else lines(x, y, col=ColScheme[i], lwd=2) - j=which(list_gcm_period()==period_focal)+1 - points(x,y, pch=16, col=ColScheme[i], cex=0.5) - points(x[j],y[j], pch=21, bg=ColScheme[i], cex=2.5) - points(x.runs,y.runs, pch=21, bg=ColScheme[i], cex=1) - text(x[j],y[j], substr(gcm_models, 1, 2)[i], cex=0.5, font=2) + if(interactive == F){ + # BASE PLOT + # initiate the plot + par(mar=c(3,4,0,1), mgp=c(1.25, 0.25,0), cex=1.5) + plot(data_mean$xanom,data_mean$yanom,col="white", tck=0, xaxt="n", yaxt="n", ylab="", + xlab=paste("Change in", variables_lookup$Variable[which(variables_lookup$Code==xvar)]) + ) + par(mgp=c(2.5,0.25, 0)) + title(ylab=paste("Change in", variables_lookup$Variable[which(variables_lookup$Code==yvar)])) + lines(c(0,0), c(-99,99), lty=2, col="gray") + lines(c(-99,99), c(0,0), lty=2, col="gray") + + if(show_ensMean) points(ensMean$xanom,ensMean$yanom, pch=43, col="gray", cex=3) + if(show_observed) points(obs$xanom ,obs$yanom, pch=22, bg="gray", cex=2.5) + # text(x1,y1, "2001-2020", cex=1.15, font=2, pos=4, col="gray", offset=0.9) + + # GCM projections + gcm=list_gcm()[1] + for(gcm in gcm_models){ + i=which(gcm_models==gcm) + x2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", xanom]) + y2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", yanom]) + x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] + y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] + if(show_trajectories){ + if(length(unique(sign(diff(x2))))==1){ + x3 <- if(unique(sign(diff(x2)))==-1) rev(x2) else x2 + y3 <- if(unique(sign(diff(x2)))==-1) rev(y2) else y2 + s <- stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data$xanom))/500)) # way better than interpSpline, not prone to oscillations + lines(s, col=ColScheme[i], lwd=2) + } else lines(x2, y2, col=ColScheme[i], lwd=2) + points(x2,y2, pch=16, col=ColScheme[i], cex=0.5) + } + j=which(list_gcm_period()==period_focal)+1 + points(x2[j],y2[j], pch=21, bg=ColScheme[i], cex=2.5) + if(show_runs) points(x.runs,y.runs, pch=21, bg=ColScheme[i], cex=1) + text(x2[j],y2[j], substr(gcm_models, 1, 2)[i], cex=0.5, font=2) + } + + # axis labels + axis(1, at=pretty(data_mean$xanom), labels=if(xvar_type=="ratio") paste(pretty(data_mean$xanom)*100, "%", sep="") else pretty(data_mean$xanom), tck=0) + axis(2, at=pretty(data_mean$yanom), labels=if(yvar_type=="ratio") paste(pretty(data_mean$yanom)*100, "%", sep="") else pretty(data_mean$yanom), las=2, tck=0) + + # Legend + s <- c(show_observed, show_runs, T, show_trajectories, show_ensMean) + legend(legend_pos, + legend = c("Observed climate (2001-2020)", "individual GCM simulation", "GCM mean", "GCM mean trajectory", "Ensemble mean")[s], + pch = c(22, 21, 21, 16, 43)[s], + pt.cex = c(2.5, 1, 2.5, 0.5, 3)[s], + pt.bg = c("gray", "gray", "gray", NA, NA)[s], + col = c(1,1,1,"gray","gray")[s], + lty = c(NA, NA, NA, 1, NA)[s], + lwd = c(NA, NA, NA, 2, NA)[s], + bty = "n", cex=0.8 + ) + } else { + # PLOTLY PLOT } - - # axis labels - axis(1, at=pretty(data_mean$xanom), labels=if(xvar_type=="ratio") paste(pretty(data_mean$xanom)*100, "%", sep="") else pretty(data_mean$xanom), tck=0) - axis(2, at=pretty(data_mean$yanom), labels=if(yvar_type=="ratio") paste(pretty(data_mean$yanom)*100, "%", sep="") else pretty(data_mean$yanom), las=2, tck=0) - - # Legend TODO: modify for user selection of plot elements. - legend(legend_pos, - legend = c("Observed historical climate", "individual GCM simulation", "GCM mean", "GCM mean trajectory", "Ensemble mean"), - pch = c(22, 21, 21, 16, 43), - pt.cex = c(2.5, 1, 2.5, 0.5, 3), - pt.bg = c("gray", "gray", "gray", NA, NA), - col = c(1,1,1,"gray","gray"), - lty = c(NA, NA, NA, 1, NA), - lwd = c(NA, NA, NA, 2, NA), - bty = "n", cex=0.8 - ) - } From 32844c188f1ab26f200940bb08ada0724d850591 Mon Sep 17 00:00:00 2001 From: cmahony Date: Mon, 12 Feb 2024 07:16:52 -0800 Subject: [PATCH 03/17] minor revisions --- R/plot_bivariate.R | 66 ++++++++++++++++++++++++++++++---------------- 1 file changed, 44 insertions(+), 22 deletions(-) diff --git a/R/plot_bivariate.R b/R/plot_bivariate.R index 79134d8e..5d7a4589 100644 --- a/R/plot_bivariate.R +++ b/R/plot_bivariate.R @@ -8,6 +8,8 @@ #' #' The climate change trajectories provided by `show_trajectories` are connected with an interpolation spline when the x variable is monotonic; otherwise the trajectory points are connected by straight lines. #' +#' This plot is designed to be used with a single SSP scenario. If multiple scenarios are passed to the plot, the GCM means and ensemble mean are averaged across the scenarios, but the individual runs for all scenarios are plotted separately. +#' #' @template xyz #' @param xvar character. x-axis variable. options are `list_variables()`. #' @param yvar character. y-axis variable. options are `list_variables()`. @@ -25,14 +27,25 @@ #' @importFrom stinepack stinterp TODO: add this package to dependencies #' #' @examples { -#' library(data.table) -#' TODO: add example usage +#' library(data.table) # TODO: what do i need to do with library(plotly) and library(stinepack)? +#' +#' # data frame of arbitrary points on vancouver island +#' my_points <- data.frame(lon = c(-123.4404, -123.5064, -124.2317), +#' lat = c(48.52631, 48.46807, 49.21999), +#' elev = c(52, 103, 357), +#' id = LETTERS[1:3] +#' ) +#' +#' # plot without export +#' plot_bivariate(my_points) +#' +#' # export plot to the working directory +#' png(filename="plot_test.png", type="cairo", units="in", width=6, height=5, pointsize=10, res=300) +#' plot_bivariate(my_points) +#' dev.off() #' } #' @export -# Test_Locations_VanIsl <- read.csv("data-raw/Test_Locations_VanIsl.csv") -# xyz <- data.frame(Test_Locations_VanIsl, id=1:dim(Test_Locations_VanIsl)[1]) -# names(xyz) <- c("lat", "lon", "elev", "id") plot_bivariate <- function( xyz, xvar = "Tave_sm", @@ -42,7 +55,7 @@ plot_bivariate <- function( period_focal = list_gcm_period()[1], gcm_models = list_gcm()[c(1,4,5,6,7,10,11,12)], ## TODO: check that i am passing this to climr_downscale properly ssp = list_ssp()[2], ## TODO: check that i am passing this to climr_downscale properly - legend_pos = "topright", + legend_pos = "bottomleft", show_runs = T, show_ensMean = T, show_observed = T, @@ -51,11 +64,11 @@ plot_bivariate <- function( ... ) { - variables_lookup <- read.csv("data-raw/variables_lookup.csv") + data("variables") # variable types for default scaling (percent or absolute) - xvar_type <- variables_lookup$Type[which(variables_lookup$Code==xvar)] - yvar_type <- variables_lookup$Type[which(variables_lookup$Code==yvar)] + xvar_type <- variables$Scale[which(variables$Code==xvar)] + yvar_type <- variables$Scale[which(variables$Code==yvar)] colors = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)][-1] set.seed(2) @@ -73,11 +86,11 @@ plot_bivariate <- function( ) # convert absolute values to anomalies - data[, xanom := if(xvar_type=="ratio") (get(xvar)/get(xvar)[1]-1) else (get(xvar) - get(xvar)[1]), by = id] - data[, yanom := if(yvar_type=="ratio") (get(yvar)/get(yvar)[1]-1) else (get(yvar) - get(yvar)[1]), by = id] + data[, xanom := if(xvar_type=="Log") (get(xvar)/get(xvar)[1]-1) else (get(xvar) - get(xvar)[1]), by = id] + data[, yanom := if(yvar_type=="Log") (get(yvar)/get(yvar)[1]-1) else (get(yvar) - get(yvar)[1]), by = id] # collapse the points down to a mean anomaly - data <- data[, .(xanom = mean(xanom), yanom = mean(yanom)), by = .(GCM, SSP, RUN, PERIOD)] + data <- data[, .(xanom = mean(xanom), yanom = mean(yanom)), by = .(GCM, SSP, RUN, PERIOD)] # TODO: need to give special treatment to SSPs; perhaps by averaging the ensemble mean but retaining the individual runs # ensemble mean for the selected period ensMean <- data[!is.na(GCM) & RUN == "ensembleMean" & PERIOD == period_focal, .(xanom = mean(xanom), yanom = mean(yanom)), ] # observed climate @@ -88,10 +101,10 @@ plot_bivariate <- function( # initiate the plot par(mar=c(3,4,0,1), mgp=c(1.25, 0.25,0), cex=1.5) plot(data_mean$xanom,data_mean$yanom,col="white", tck=0, xaxt="n", yaxt="n", ylab="", - xlab=paste("Change in", variables_lookup$Variable[which(variables_lookup$Code==xvar)]) + xlab=paste("Change in", variables$Variable[which(variables$Code==xvar)]) ) par(mgp=c(2.5,0.25, 0)) - title(ylab=paste("Change in", variables_lookup$Variable[which(variables_lookup$Code==yvar)])) + title(ylab=paste("Change in", variables$Variable[which(variables$Code==yvar)])) lines(c(0,0), c(-99,99), lty=2, col="gray") lines(c(-99,99), c(0,0), lty=2, col="gray") @@ -99,14 +112,21 @@ plot_bivariate <- function( if(show_observed) points(obs$xanom ,obs$yanom, pch=22, bg="gray", cex=2.5) # text(x1,y1, "2001-2020", cex=1.15, font=2, pos=4, col="gray", offset=0.9) - # GCM projections - gcm=list_gcm()[1] + # plot individual runs + if(show_runs){ + for(gcm in gcm_models){ + i=which(gcm_models==gcm) + x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] + y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] + points(x.runs,y.runs, pch=21, bg=ColScheme[i], cex=1) + } + } + + # plot model means and trajectories for(gcm in gcm_models){ i=which(gcm_models==gcm) x2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", xanom]) y2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", yanom]) - x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] - y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] if(show_trajectories){ if(length(unique(sign(diff(x2))))==1){ x3 <- if(unique(sign(diff(x2)))==-1) rev(x2) else x2 @@ -118,20 +138,19 @@ plot_bivariate <- function( } j=which(list_gcm_period()==period_focal)+1 points(x2[j],y2[j], pch=21, bg=ColScheme[i], cex=2.5) - if(show_runs) points(x.runs,y.runs, pch=21, bg=ColScheme[i], cex=1) text(x2[j],y2[j], substr(gcm_models, 1, 2)[i], cex=0.5, font=2) } # axis labels - axis(1, at=pretty(data_mean$xanom), labels=if(xvar_type=="ratio") paste(pretty(data_mean$xanom)*100, "%", sep="") else pretty(data_mean$xanom), tck=0) - axis(2, at=pretty(data_mean$yanom), labels=if(yvar_type=="ratio") paste(pretty(data_mean$yanom)*100, "%", sep="") else pretty(data_mean$yanom), las=2, tck=0) + axis(1, at=pretty(data_mean$xanom), labels=if(xvar_type=="Log") paste(pretty(data_mean$xanom)*100, "%", sep="") else pretty(data_mean$xanom), tck=0) + axis(2, at=pretty(data_mean$yanom), labels=if(yvar_type=="Log") paste(pretty(data_mean$yanom)*100, "%", sep="") else pretty(data_mean$yanom), las=2, tck=0) # Legend s <- c(show_observed, show_runs, T, show_trajectories, show_ensMean) legend(legend_pos, legend = c("Observed climate (2001-2020)", "individual GCM simulation", "GCM mean", "GCM mean trajectory", "Ensemble mean")[s], pch = c(22, 21, 21, 16, 43)[s], - pt.cex = c(2.5, 1, 2.5, 0.5, 3)[s], + pt.cex = c(2, .8, 2, 0.5, 2.2)[s], pt.bg = c("gray", "gray", "gray", NA, NA)[s], col = c(1,1,1,"gray","gray")[s], lty = c(NA, NA, NA, 1, NA)[s], @@ -143,3 +162,6 @@ plot_bivariate <- function( } } + + + From c33b628de02c99cce917c1e81d60d15c1a7c91b4 Mon Sep 17 00:00:00 2001 From: cmahony Date: Mon, 12 Feb 2024 14:24:02 -0800 Subject: [PATCH 04/17] add the plotly plot --- R/plot_bivariate.R | 98 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 83 insertions(+), 15 deletions(-) diff --git a/R/plot_bivariate.R b/R/plot_bivariate.R index 5d7a4589..c71ce1f5 100644 --- a/R/plot_bivariate.R +++ b/R/plot_bivariate.R @@ -6,7 +6,7 @@ #' @details #' The input table `xyz` can be a single location or multiple locations. If multiple locations, the plot provides the mean of the anomalies for these locations. #' -#' The climate change trajectories provided by `show_trajectories` are connected with an interpolation spline when the x variable is monotonic; otherwise the trajectory points are connected by straight lines. +#' The climate change trajectories provided by `show_trajectories` are points for each of the five 20-year periods specified by `list_gcm_period()`. These points are connected with an interpolation spline when the x variable is monotonic; otherwise the trajectory points are connected by straight lines. #' #' This plot is designed to be used with a single SSP scenario. If multiple scenarios are passed to the plot, the GCM means and ensemble mean are averaged across the scenarios, but the individual runs for all scenarios are plotted separately. #' @@ -19,22 +19,22 @@ #' @param show_ensMean logical. If TRUE, the multi-model ensemble mean is plotted (for `period_focal` only). #' @param show_observed logical. If TRUE, the 2001-2020 observed climate is plotted. #' @param show_trajectories logical. If TRUE, the values of the single-model ensemble mean are plotted for all 20-year periods in `list_gcm_period()`, connected by an interpolation spline. -#' @param interactive logical. If TRUE, an interactive plot is generated using plotly(). If FALSE, a plot is generated using base graphics. +#' @param interactive logical. If TRUE, an interactive plot is generated using `{plotly}`. If FALSE, a plot is generated using base graphics. #' #' @return TODO #' #' @importFrom data.table TODO: find out what i need to do to get this command right. -#' @importFrom stinepack stinterp TODO: add this package to dependencies +#' @importFrom stinepack stinterp TODO: add this package to dependencies (as a suggest) #' #' @examples { #' library(data.table) # TODO: what do i need to do with library(plotly) and library(stinepack)? #' #' # data frame of arbitrary points on vancouver island -#' my_points <- data.frame(lon = c(-123.4404, -123.5064, -124.2317), -#' lat = c(48.52631, 48.46807, 49.21999), -#' elev = c(52, 103, 357), -#' id = LETTERS[1:3] -#' ) +# my_points <- data.frame(lon = c(-123.4404, -123.5064, -124.2317), +# lat = c(48.52631, 48.46807, 49.21999), +# elev = c(52, 103, 357), +# id = LETTERS[1:3] +# ) #' #' # plot without export #' plot_bivariate(my_points) @@ -46,6 +46,7 @@ #' } #' @export +## TODO: structure the climr_downscale calls so that a user can't break the plot. plot_bivariate <- function( xyz, xvar = "Tave_sm", @@ -60,7 +61,7 @@ plot_bivariate <- function( show_ensMean = T, show_observed = T, show_trajectories = T, - interactive = F, ## TODO: add a plotly version of the plot + interactive = F, ... ) { @@ -78,8 +79,8 @@ plot_bivariate <- function( data <- climr_downscale(xyz, which_normal = "auto", historic_period = list_historic()[1], - gcm_models = list_gcm()[c(1,4,5,6,7,10,11,12)], ## TODO: this is wrong. need to recieve this as passed from plot_bivariate() - ssp = list_ssp()[2], ## TODO: this is wrong. need to receive this as passed from plot_bivariate() + gcm_models = gcm_models, ## TODO: this is wrong. need to recieve this as passed from plot_bivariate() + ssp = ssp, ## TODO: this is wrong. need to receive this as passed from plot_bivariate() gcm_period = list_gcm_period(), max_run = 10, vars = c(xvar, yvar) @@ -97,10 +98,11 @@ plot_bivariate <- function( obs <- data[is.na(GCM) & PERIOD == period_focal] if(interactive == F){ + # BASE PLOT # initiate the plot par(mar=c(3,4,0,1), mgp=c(1.25, 0.25,0), cex=1.5) - plot(data_mean$xanom,data_mean$yanom,col="white", tck=0, xaxt="n", yaxt="n", ylab="", + plot(data$xanom,data$yanom,col="white", tck=0, xaxt="n", yaxt="n", ylab="", xlab=paste("Change in", variables$Variable[which(variables$Code==xvar)]) ) par(mgp=c(2.5,0.25, 0)) @@ -142,8 +144,8 @@ plot_bivariate <- function( } # axis labels - axis(1, at=pretty(data_mean$xanom), labels=if(xvar_type=="Log") paste(pretty(data_mean$xanom)*100, "%", sep="") else pretty(data_mean$xanom), tck=0) - axis(2, at=pretty(data_mean$yanom), labels=if(yvar_type=="Log") paste(pretty(data_mean$yanom)*100, "%", sep="") else pretty(data_mean$yanom), las=2, tck=0) + axis(1, at=pretty(data$xanom), labels=if(xvar_type=="Log") paste(pretty(data$xanom)*100, "%", sep="") else pretty(data$xanom), tck=0) + axis(2, at=pretty(data$yanom), labels=if(yvar_type=="Log") paste(pretty(data$yanom)*100, "%", sep="") else pretty(data$yanom), las=2, tck=0) # Legend s <- c(show_observed, show_runs, T, show_trajectories, show_ensMean) @@ -157,11 +159,77 @@ plot_bivariate <- function( lwd = c(NA, NA, NA, 2, NA)[s], bty = "n", cex=0.8 ) + } else { + # PLOTLY PLOT + library(plotly) # TODO: figure out how to manage dependencies + # TODO: the colors aren't plotting correctly. need to fix this. + + #initiate the plot + fig <- plot_ly(x=data$xanom,y=data$yanom, type = 'scatter', mode = 'markers', marker = list(color ="lightgray", size=5), hoverinfo="none", color="All models/scenarios/runs/periods") + + # axis titles + fig <- fig %>% layout(xaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==xvar)]), range=range(data$xanom)), + yaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==yvar)]), range=range(data$yanom)) + ) + + # observed climate + fig <- fig %>% add_markers(obs$xanom ,obs$yanom, name="Observed Climate (2001-2020)", text="observed\n(2001-2020)", hoverinfo="text", + marker = list(size = 25, color = "grey"), symbol = 43) + + # ensemble mean + fig <- fig %>% add_markers(ensMean$xanom,ensMean$yanom, name="Ensemble mean", text="Ensemble mean", hoverinfo="text", + marker = list(size = 20, color = "grey", symbol = 3)) + + # plot individual runs + if(show_runs){ + for(gcm in gcm_models){ + i=which(gcm_models==gcm) + x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] + y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] + runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, RUN] + fig <- fig %>% add_markers(x=x.runs,y=y.runs, color = ColScheme[i], name="Individual GCM runs", text=paste(gcm_models[i], runs), hoverinfo="text", showlegend = F, + marker = list(size = 7, color = ColScheme[i], line = list(color = "black", width = 1)), legendgroup=paste("group", i, sep="")) + } + } + + # GCM mean trajectories + # plot model means and trajectories + gcm = gcm_models[2] + for(gcm in gcm_models){ + i=which(gcm_models==gcm) + x2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", xanom]) + y2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", yanom]) + if(show_trajectories){ + if(length(unique(sign(diff(x2))))==1){ + x3 <- if(unique(sign(diff(x2)))==-1) rev(x2) else x2 + y3 <- if(unique(sign(diff(x2)))==-1) rev(y2) else y2 + s <- stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data$xanom))/500)) # way better than interpSpline, not prone to oscillations + fig <- fig %>% add_trace(x=s$x, y=s$y, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) + } else { + fig <- fig %>% add_trace(x=x2, y=y2, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) + } + fig <- fig %>% add_markers(x=x2,y=y2, color = ColScheme[i], text=gcm_models[i], hoverinfo="text", + marker = list(size = 8, color = ColScheme[i]), legendgroup=paste("group", i, sep=""), showlegend = FALSE) + } + j=which(list_gcm_period()==period_focal)+1 + fig <- fig %>% add_markers(x2[j],y2[j], color = gcm_models[i], colors=ColScheme[i], text=gcm_models[i], + marker = list(size = 20, color = ColScheme[i], line = list(color = "black", width = 1)), + legendgroup=paste("group", i, sep="")) + + fig <- fig %>% add_annotations(x=x2[j],y=y2[j], text = sprintf("%s", substr(gcm_models, 1, 2)[i]), xanchor = 'center', yanchor = 'center', showarrow = F, + legendgroup=paste("group", i, sep="")) + } + + if(xvar_type=="Log") fig <- fig %>% layout(xaxis = list(tickformat = "%")) + if(yvar_type=="Log") fig <- fig %>% layout(yaxis = list(tickformat = "%")) + + fig + } } - +plot_bivariate(my_points, xvar = "PPT_sm", yvar = "Tave_sm") From 10b8ac2514931b60f3ec8ff35a4eabcb42b3b76e Mon Sep 17 00:00:00 2001 From: cmahony Date: Mon, 12 Feb 2024 15:12:29 -0800 Subject: [PATCH 05/17] redoc --- NAMESPACE | 1 + R/plot_bivariate.R | 358 +++++++++++++++++++++--------------------- man/get_bb.Rd | 5 +- man/plot_bivariate.Rd | 95 +++++++++++ 4 files changed, 280 insertions(+), 179 deletions(-) create mode 100644 man/plot_bivariate.Rd diff --git a/NAMESPACE b/NAMESPACE index f5145e0b..7671130c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -21,6 +21,7 @@ export(list_ssp) export(list_variables) export(normal_input) export(pgGetTerra) +export(plot_bivariate) import(data.table) import(uuid) importFrom(DBI,dbGetQuery) diff --git a/R/plot_bivariate.R b/R/plot_bivariate.R index c71ce1f5..9eca4eed 100644 --- a/R/plot_bivariate.R +++ b/R/plot_bivariate.R @@ -5,231 +5,235 @@ #' #' @details #' The input table `xyz` can be a single location or multiple locations. If multiple locations, the plot provides the mean of the anomalies for these locations. -#' #' The climate change trajectories provided by `show_trajectories` are points for each of the five 20-year periods specified by `list_gcm_period()`. These points are connected with an interpolation spline when the x variable is monotonic; otherwise the trajectory points are connected by straight lines. -#' #' This plot is designed to be used with a single SSP scenario. If multiple scenarios are passed to the plot, the GCM means and ensemble mean are averaged across the scenarios, but the individual runs for all scenarios are plotted separately. #' #' @template xyz #' @param xvar character. x-axis variable. options are `list_variables()`. #' @param yvar character. y-axis variable. options are `list_variables()`. #' @param period_focal character. The 20-year period for which to plot the ensemble detail. options are `list_gcm_period()`. -#' @param legend_pos character. Position of the legend. Options are c("bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right", "center"). +#' @inheritParams climr_downscale +#' @param legend_pos character. Position of the legend. Options are `c("bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right", "center")`. #' @param show_runs logical. If TRUE, the individual runs of the model are plotted (for `period_focal` only) in addition to the single-model ensemble mean. #' @param show_ensMean logical. If TRUE, the multi-model ensemble mean is plotted (for `period_focal` only). #' @param show_observed logical. If TRUE, the 2001-2020 observed climate is plotted. #' @param show_trajectories logical. If TRUE, the values of the single-model ensemble mean are plotted for all 20-year periods in `list_gcm_period()`, connected by an interpolation spline. #' @param interactive logical. If TRUE, an interactive plot is generated using `{plotly}`. If FALSE, a plot is generated using base graphics. #' -#' @return TODO +#' @return NULL. Draws a plot in the active graphics device. #' -#' @importFrom data.table TODO: find out what i need to do to get this command right. -#' @importFrom stinepack stinterp TODO: add this package to dependencies (as a suggest) -#' -#' @examples { -#' library(data.table) # TODO: what do i need to do with library(plotly) and library(stinepack)? -#' +#' @examples #' # data frame of arbitrary points on vancouver island -# my_points <- data.frame(lon = c(-123.4404, -123.5064, -124.2317), -# lat = c(48.52631, 48.46807, 49.21999), -# elev = c(52, 103, 357), -# id = LETTERS[1:3] -# ) +#' my_points <- data.frame(lon = c(-123.4404, -123.5064, -124.2317), +#' lat = c(48.52631, 48.46807, 49.21999), +#' elev = c(52, 103, 357), +#' id = LETTERS[1:3] +#' ) #' -#' # plot without export +#' # draw the plot #' plot_bivariate(my_points) #' #' # export plot to the working directory #' png(filename="plot_test.png", type="cairo", units="in", width=6, height=5, pointsize=10, res=300) #' plot_bivariate(my_points) #' dev.off() -#' } +#' #' @export -## TODO: structure the climr_downscale calls so that a user can't break the plot. plot_bivariate <- function( xyz, xvar = "Tave_sm", yvar = "PPT_sm", # percent_x = NULL, TODO: set up an override for ratio variables being expressed as percent anomalies # percent_y = NULL, TODO: set up an override for ratio variables being expressed as percent anomalies + which_normal = "auto", # TODO: remove once the bug is fixed. period_focal = list_gcm_period()[1], - gcm_models = list_gcm()[c(1,4,5,6,7,10,11,12)], ## TODO: check that i am passing this to climr_downscale properly - ssp = list_ssp()[2], ## TODO: check that i am passing this to climr_downscale properly + gcm_models = list_gcm()[c(1,4,5,6,7,10,11,12)], + ssp = list_ssp()[2], + historic_period = list_historic()[1], + gcm_period = list_gcm_period(), + max_run = 10, legend_pos = "bottomleft", - show_runs = T, - show_ensMean = T, - show_observed = T, - show_trajectories = T, - interactive = F, - ... + show_runs = TRUE, + show_ensMean = TRUE, + show_observed = TRUE, + show_trajectories = TRUE, + interactive = FALSE, + cache = TRUE ) { - data("variables") - - # variable types for default scaling (percent or absolute) - xvar_type <- variables$Scale[which(variables$Code==xvar)] - yvar_type <- variables$Scale[which(variables$Code==yvar)] - - colors = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)][-1] - set.seed(2) - ColScheme <- sample(colors,length(gcm_models)) # TODO select better colors - - # generate the climate data - data <- climr_downscale(xyz, - which_normal = "auto", - historic_period = list_historic()[1], - gcm_models = gcm_models, ## TODO: this is wrong. need to recieve this as passed from plot_bivariate() - ssp = ssp, ## TODO: this is wrong. need to receive this as passed from plot_bivariate() - gcm_period = list_gcm_period(), - max_run = 10, - vars = c(xvar, yvar) - ) - - # convert absolute values to anomalies - data[, xanom := if(xvar_type=="Log") (get(xvar)/get(xvar)[1]-1) else (get(xvar) - get(xvar)[1]), by = id] - data[, yanom := if(yvar_type=="Log") (get(yvar)/get(yvar)[1]-1) else (get(yvar) - get(yvar)[1]), by = id] - - # collapse the points down to a mean anomaly - data <- data[, .(xanom = mean(xanom), yanom = mean(yanom)), by = .(GCM, SSP, RUN, PERIOD)] # TODO: need to give special treatment to SSPs; perhaps by averaging the ensemble mean but retaining the individual runs - # ensemble mean for the selected period - ensMean <- data[!is.na(GCM) & RUN == "ensembleMean" & PERIOD == period_focal, .(xanom = mean(xanom), yanom = mean(yanom)), ] - # observed climate - obs <- data[is.na(GCM) & PERIOD == period_focal] - - if(interactive == F){ - - # BASE PLOT - # initiate the plot - par(mar=c(3,4,0,1), mgp=c(1.25, 0.25,0), cex=1.5) - plot(data$xanom,data$yanom,col="white", tck=0, xaxt="n", yaxt="n", ylab="", - xlab=paste("Change in", variables$Variable[which(variables$Code==xvar)]) - ) - par(mgp=c(2.5,0.25, 0)) - title(ylab=paste("Change in", variables$Variable[which(variables$Code==yvar)])) - lines(c(0,0), c(-99,99), lty=2, col="gray") - lines(c(-99,99), c(0,0), lty=2, col="gray") - - if(show_ensMean) points(ensMean$xanom,ensMean$yanom, pch=43, col="gray", cex=3) - if(show_observed) points(obs$xanom ,obs$yanom, pch=22, bg="gray", cex=2.5) - # text(x1,y1, "2001-2020", cex=1.15, font=2, pos=4, col="gray", offset=0.9) - - # plot individual runs - if(show_runs){ - for(gcm in gcm_models){ - i=which(gcm_models==gcm) - x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] - y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] - points(x.runs,y.runs, pch=21, bg=ColScheme[i], cex=1) - } - } - - # plot model means and trajectories - for(gcm in gcm_models){ - i=which(gcm_models==gcm) - x2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", xanom]) - y2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", yanom]) - if(show_trajectories){ - if(length(unique(sign(diff(x2))))==1){ - x3 <- if(unique(sign(diff(x2)))==-1) rev(x2) else x2 - y3 <- if(unique(sign(diff(x2)))==-1) rev(y2) else y2 - s <- stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data$xanom))/500)) # way better than interpSpline, not prone to oscillations - lines(s, col=ColScheme[i], lwd=2) - } else lines(x2, y2, col=ColScheme[i], lwd=2) - points(x2,y2, pch=16, col=ColScheme[i], cex=0.5) - } - j=which(list_gcm_period()==period_focal)+1 - points(x2[j],y2[j], pch=21, bg=ColScheme[i], cex=2.5) - text(x2[j],y2[j], substr(gcm_models, 1, 2)[i], cex=0.5, font=2) - } - - # axis labels - axis(1, at=pretty(data$xanom), labels=if(xvar_type=="Log") paste(pretty(data$xanom)*100, "%", sep="") else pretty(data$xanom), tck=0) - axis(2, at=pretty(data$yanom), labels=if(yvar_type=="Log") paste(pretty(data$yanom)*100, "%", sep="") else pretty(data$yanom), las=2, tck=0) - - # Legend - s <- c(show_observed, show_runs, T, show_trajectories, show_ensMean) - legend(legend_pos, - legend = c("Observed climate (2001-2020)", "individual GCM simulation", "GCM mean", "GCM mean trajectory", "Ensemble mean")[s], - pch = c(22, 21, 21, 16, 43)[s], - pt.cex = c(2, .8, 2, 0.5, 2.2)[s], - pt.bg = c("gray", "gray", "gray", NA, NA)[s], - col = c(1,1,1,"gray","gray")[s], - lty = c(NA, NA, NA, 1, NA)[s], - lwd = c(NA, NA, NA, 2, NA)[s], - bty = "n", cex=0.8 - ) - + if(!requireNamespace("stinepack")){ + stop("package stinepack must be installed to use this function") } else { - # PLOTLY PLOT - library(plotly) # TODO: figure out how to manage dependencies - # TODO: the colors aren't plotting correctly. need to fix this. - - #initiate the plot - fig <- plot_ly(x=data$xanom,y=data$yanom, type = 'scatter', mode = 'markers', marker = list(color ="lightgray", size=5), hoverinfo="none", color="All models/scenarios/runs/periods") - - # axis titles - fig <- fig %>% layout(xaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==xvar)]), range=range(data$xanom)), - yaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==yvar)]), range=range(data$yanom)) + data("variables") + + # variable types for default scaling (percent or absolute) + xvar_type <- variables$Scale[which(variables$Code==xvar)] + yvar_type <- variables$Scale[which(variables$Code==yvar)] + + colors = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = TRUE)][-1] + set.seed(2) + ColScheme <- sample(colors,length(gcm_models)) # TODO select better colors + + # generate the climate data + data <- climr_downscale(xyz, + which_normal = which_normal, # TODO: remove once the bug is fixed. + historic_period = historic_period, + gcm_models = gcm_models, + ssp = ssp, + gcm_period = gcm_period, + max_run = max_run, + vars = c(xvar, yvar), + cache = cache ) - # observed climate - fig <- fig %>% add_markers(obs$xanom ,obs$yanom, name="Observed Climate (2001-2020)", text="observed\n(2001-2020)", hoverinfo="text", - marker = list(size = 25, color = "grey"), symbol = 43) + # convert absolute values to anomalies + data[, xanom := if(xvar_type=="Log") (get(xvar)/get(xvar)[1]-1) else (get(xvar) - get(xvar)[1]), by = id] + data[, yanom := if(yvar_type=="Log") (get(yvar)/get(yvar)[1]-1) else (get(yvar) - get(yvar)[1]), by = id] - # ensemble mean - fig <- fig %>% add_markers(ensMean$xanom,ensMean$yanom, name="Ensemble mean", text="Ensemble mean", hoverinfo="text", - marker = list(size = 20, color = "grey", symbol = 3)) + # collapse the points down to a mean anomaly + data <- data[, .(xanom = mean(xanom), yanom = mean(yanom)), by = .(GCM, SSP, RUN, PERIOD)] # TODO: need to give special treatment to SSPs; perhaps by averaging the ensemble mean but retaining the individual runs + # ensemble mean for the selected period + ensMean <- data[!is.na(GCM) & RUN == "ensembleMean" & PERIOD == period_focal, .(xanom = mean(xanom), yanom = mean(yanom)), ] + # observed climate + obs <- data[is.na(GCM) & PERIOD == period_focal] - # plot individual runs - if(show_runs){ + if(interactive == FALSE){ + + # BASE PLOT + # initiate the plot + par(mar=c(3,4,0,1), mgp=c(1.25, 0.25,0), cex=1.5) + plot(data$xanom,data$yanom,col="white", tck=0, xaxt="n", yaxt="n", ylab="", + xlab=paste("Change in", variables$Variable[which(variables$Code==xvar)]) + ) + par(mgp=c(2.5,0.25, 0)) + title(ylab=paste("Change in", variables$Variable[which(variables$Code==yvar)])) + lines(c(0,0), c(-99,99), lty=2, col="gray") + lines(c(-99,99), c(0,0), lty=2, col="gray") + + if(show_ensMean) points(ensMean$xanom,ensMean$yanom, pch=43, col="gray", cex=3) + if(show_observed) points(obs$xanom ,obs$yanom, pch=22, bg="gray", cex=2.5) + # text(x1,y1, "2001-2020", cex=1.15, font=2, pos=4, col="gray", offset=0.9) + + # plot individual runs + if(show_runs){ + for(gcm in gcm_models){ + i=which(gcm_models==gcm) + x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] + y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] + points(x.runs,y.runs, pch=21, bg=ColScheme[i], cex=1) + } + } + + # plot model means and trajectories for(gcm in gcm_models){ i=which(gcm_models==gcm) - x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] - y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] - runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, RUN] - fig <- fig %>% add_markers(x=x.runs,y=y.runs, color = ColScheme[i], name="Individual GCM runs", text=paste(gcm_models[i], runs), hoverinfo="text", showlegend = F, - marker = list(size = 7, color = ColScheme[i], line = list(color = "black", width = 1)), legendgroup=paste("group", i, sep="")) - } - } - - # GCM mean trajectories - # plot model means and trajectories - gcm = gcm_models[2] - for(gcm in gcm_models){ - i=which(gcm_models==gcm) - x2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", xanom]) - y2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", yanom]) - if(show_trajectories){ - if(length(unique(sign(diff(x2))))==1){ - x3 <- if(unique(sign(diff(x2)))==-1) rev(x2) else x2 - y3 <- if(unique(sign(diff(x2)))==-1) rev(y2) else y2 - s <- stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data$xanom))/500)) # way better than interpSpline, not prone to oscillations - fig <- fig %>% add_trace(x=s$x, y=s$y, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) - } else { - fig <- fig %>% add_trace(x=x2, y=y2, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) + x2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", xanom]) + y2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", yanom]) + if(show_trajectories){ + if(length(unique(sign(diff(x2))))==1){ + x3 <- if(unique(sign(diff(x2)))==-1) rev(x2) else x2 + y3 <- if(unique(sign(diff(x2)))==-1) rev(y2) else y2 + s <- stinepack::stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data$xanom))/500)) + lines(s, col=ColScheme[i], lwd=2) + } else lines(x2, y2, col=ColScheme[i], lwd=2) + points(x2,y2, pch=16, col=ColScheme[i], cex=0.5) } - fig <- fig %>% add_markers(x=x2,y=y2, color = ColScheme[i], text=gcm_models[i], hoverinfo="text", - marker = list(size = 8, color = ColScheme[i]), legendgroup=paste("group", i, sep=""), showlegend = FALSE) + j=which(list_gcm_period()==period_focal)+1 + points(x2[j],y2[j], pch=21, bg=ColScheme[i], cex=2.5) + text(x2[j],y2[j], substr(gcm_models, 1, 2)[i], cex=0.5, font=2) } - j=which(list_gcm_period()==period_focal)+1 - fig <- fig %>% add_markers(x2[j],y2[j], color = gcm_models[i], colors=ColScheme[i], text=gcm_models[i], - marker = list(size = 20, color = ColScheme[i], line = list(color = "black", width = 1)), - legendgroup=paste("group", i, sep="")) - fig <- fig %>% add_annotations(x=x2[j],y=y2[j], text = sprintf("%s", substr(gcm_models, 1, 2)[i]), xanchor = 'center', yanchor = 'center', showarrow = F, + # axis labels + axis(1, at=pretty(data$xanom), labels=if(xvar_type=="Log") paste(pretty(data$xanom)*100, "%", sep="") else pretty(data$xanom), tck=0) + axis(2, at=pretty(data$yanom), labels=if(yvar_type=="Log") paste(pretty(data$yanom)*100, "%", sep="") else pretty(data$yanom), las=2, tck=0) + + # Legend + s <- c(show_observed, show_runs, TRUE, show_trajectories, show_ensMean) + legend(legend_pos, + legend = c("Observed climate (2001-2020)", "individual GCM simulation", "GCM mean", "GCM mean trajectory", "Ensemble mean")[s], + pch = c(22, 21, 21, 16, 43)[s], + pt.cex = c(2, .8, 2, 0.5, 2.2)[s], + pt.bg = c("gray", "gray", "gray", NA, NA)[s], + col = c(1,1,1,"gray","gray")[s], + lty = c(NA, NA, NA, 1, NA)[s], + lwd = c(NA, NA, NA, 2, NA)[s], + bty = "n", cex=0.8 + ) + + } else { + + if(!require(plotly)){ + stop("package 'plotly' must be installed when 'interactive==TRUE'") + } else { + + # PLOTLY PLOT + # TODO: the colors aren't plotting correctly. need to fix this. + + #initiate the plot + fig <- plot_ly(x=data$xanom,y=data$yanom, type = 'scatter', mode = 'markers', marker = list(color ="lightgray", size=5), hoverinfo="none", color="All models/scenarios/runs/periods") + + # axis titles + fig <- fig %>% layout(xaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==xvar)]), range=range(data$xanom)), + yaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==yvar)]), range=range(data$yanom)) + ) + + # observed climate + fig <- fig %>% add_markers(obs$xanom ,obs$yanom, name="Observed Climate (2001-2020)", text="observed\n(2001-2020)", hoverinfo="text", + marker = list(size = 25, color = "grey"), symbol = 43) + + # ensemble mean + fig <- fig %>% add_markers(ensMean$xanom,ensMean$yanom, name="Ensemble mean", text="Ensemble mean", hoverinfo="text", + marker = list(size = 20, color = "grey", symbol = 3)) + + # plot individual runs + if(show_runs){ + for(gcm in gcm_models){ + i=which(gcm_models==gcm) + x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] + y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] + runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, RUN] + fig <- fig %>% add_markers(x=x.runs,y=y.runs, color = ColScheme[i], name="Individual GCM runs", text=paste(gcm_models[i], runs), hoverinfo="text", showlegend = FALSE, + marker = list(size = 7, color = ColScheme[i], line = list(color = "black", width = 1)), legendgroup=paste("group", i, sep="")) + } + } + + # GCM mean trajectories + # plot model means and trajectories + gcm = gcm_models[2] + for(gcm in gcm_models){ + i=which(gcm_models==gcm) + x2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", xanom]) + y2 <- c(0, data[GCM==gcm & RUN == "ensembleMean", yanom]) + if(show_trajectories){ + if(length(unique(sign(diff(x2))))==1){ + x3 <- if(unique(sign(diff(x2)))==-1) rev(x2) else x2 + y3 <- if(unique(sign(diff(x2)))==-1) rev(y2) else y2 + s <- stinepack::stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data$xanom))/500)) # way better than interpSpline, not prone to oscillations + fig <- fig %>% add_trace(x=s$x, y=s$y, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) + } else { + fig <- fig %>% add_trace(x=x2, y=y2, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) + } + fig <- fig %>% add_markers(x=x2,y=y2, color = ColScheme[i], text=gcm_models[i], hoverinfo="text", + marker = list(size = 8, color = ColScheme[i]), legendgroup=paste("group", i, sep=""), showlegend = FALSE) + } + j=which(list_gcm_period()==period_focal)+1 + fig <- fig %>% add_markers(x2[j],y2[j], color = gcm_models[i], colors=ColScheme[i], text=gcm_models[i], + marker = list(size = 20, color = ColScheme[i], line = list(color = "black", width = 1)), legendgroup=paste("group", i, sep="")) + + fig <- fig %>% add_annotations(x=x2[j],y=y2[j], text = sprintf("%s", substr(gcm_models, 1, 2)[i]), xanchor = 'center', yanchor = 'center', showarrow = FALSE, + legendgroup=paste("group", i, sep="")) + } + + if(xvar_type=="Log") fig <- fig %>% layout(xaxis = list(tickformat = "%")) + if(yvar_type=="Log") fig <- fig %>% layout(yaxis = list(tickformat = "%")) + + fig + } } - - if(xvar_type=="Log") fig <- fig %>% layout(xaxis = list(tickformat = "%")) - if(yvar_type=="Log") fig <- fig %>% layout(yaxis = list(tickformat = "%")) - - fig - } } -plot_bivariate(my_points, xvar = "PPT_sm", yvar = "Tave_sm") - diff --git a/man/get_bb.Rd b/man/get_bb.Rd index a0573490..e244d552 100644 --- a/man/get_bb.Rd +++ b/man/get_bb.Rd @@ -7,10 +7,11 @@ get_bb(in_xyz) } \arguments{ -\item{in_xyz}{\code{data.table} (or \code{data.frame}) of points to downscale} +\item{in_xyz}{\code{data.table} (or \code{data.frame}) of points to downscale +with columns "lon", "lat", "elev" and "id"} } \value{ -integer. A bounding box (e.g. \code{c(51,50,-121,-122)}) +integer. A bounding box (e.g. \code{c(51, 50, -121, -122)}) } \description{ Find bounding box of data diff --git a/man/plot_bivariate.Rd b/man/plot_bivariate.Rd new file mode 100644 index 00000000..9d402f80 --- /dev/null +++ b/man/plot_bivariate.Rd @@ -0,0 +1,95 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_bivariate.R +\name{plot_bivariate} +\alias{plot_bivariate} +\title{Bivariate climate change plots} +\usage{ +plot_bivariate( + xyz, + xvar = "Tave_sm", + yvar = "PPT_sm", + which_normal = "auto", + period_focal = list_gcm_period()[1], + gcm_models = list_gcm()[c(1, 4, 5, 6, 7, 10, 11, 12)], + ssp = list_ssp()[2], + historic_period = list_historic()[1], + gcm_period = list_gcm_period(), + max_run = 10, + legend_pos = "bottomleft", + show_runs = TRUE, + show_ensMean = TRUE, + show_observed = TRUE, + show_trajectories = TRUE, + interactive = FALSE, + cache = TRUE +) +} +\arguments{ +\item{xyz}{a \code{data.frame} with the following columns "long", "lat", "elev", +and a unique "id". Any extra columns will be ignored and not output.} + +\item{xvar}{character. x-axis variable. options are \code{list_variables()}.} + +\item{yvar}{character. y-axis variable. options are \code{list_variables()}.} + +\item{which_normal}{character. Which normal layer to use. +Default is "auto", which selects the highest resolution normal for each point. +Other options are one of \code{\link[=list_normal]{list_normal()}}} + +\item{period_focal}{character. The 20-year period for which to plot the ensemble detail. options are \code{list_gcm_period()}.} + +\item{gcm_models}{character. Vector of GCM names. Options are \code{\link[=list_gcm]{list_gcm()}}. Used for gcm periods, gcm timeseries, and historic timeseries. Default \code{NULL}} + +\item{ssp}{character. Vector of labels of the shared socioeconomic pathways to use. +See available options with \code{\link[=list_ssp]{list_ssp()}}. Defaults to all SSPs available.} + +\item{historic_period}{character. Which historic period. Default \code{NULL}} + +\item{gcm_period}{character. Requested future periods. Options are \code{\link[=list_gcm_period]{list_gcm_period()}}} + +\item{max_run}{integer. Maximum number of model runs to include. +A value of 0 returns the \code{ensembleMean} only. Runs are included in the order they +are found in the models data until \code{max_run} is reached. Defaults to 0L.} + +\item{legend_pos}{character. Position of the legend. Options are \code{c("bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right", "center")}.} + +\item{show_runs}{logical. If TRUE, the individual runs of the model are plotted (for \code{period_focal} only) in addition to the single-model ensemble mean.} + +\item{show_ensMean}{logical. If TRUE, the multi-model ensemble mean is plotted (for \code{period_focal} only).} + +\item{show_observed}{logical. If TRUE, the 2001-2020 observed climate is plotted.} + +\item{show_trajectories}{logical. If TRUE, the values of the single-model ensemble mean are plotted for all 20-year periods in \code{list_gcm_period()}, connected by an interpolation spline.} + +\item{interactive}{logical. If TRUE, an interactive plot is generated using \code{{plotly}}. If FALSE, a plot is generated using base graphics.} + +\item{cache}{logical. Cache data locally? Default \code{TRUE}} +} +\value{ +NULL. Draws a plot in the active graphics device. +} +\description{ +Bivariate plots of 21st century climate change for user-selected locations and climate variables. +} +\details{ +The input table \code{xyz} can be a single location or multiple locations. If multiple locations, the plot provides the mean of the anomalies for these locations. +The climate change trajectories provided by \code{show_trajectories} are points for each of the five 20-year periods specified by \code{list_gcm_period()}. These points are connected with an interpolation spline when the x variable is monotonic; otherwise the trajectory points are connected by straight lines. +This plot is designed to be used with a single SSP scenario. If multiple scenarios are passed to the plot, the GCM means and ensemble mean are averaged across the scenarios, but the individual runs for all scenarios are plotted separately. +} +\examples{ +# data frame of arbitrary points on vancouver island +my_points <- data.frame(lon = c(-123.4404, -123.5064, -124.2317), + lat = c(48.52631, 48.46807, 49.21999), + elev = c(52, 103, 357), + id = LETTERS[1:3] +) + +# draw the plot +plot_bivariate(my_points) + +# export plot to the working directory +png(filename="plot_test.png", type="cairo", units="in", width=6, height=5, pointsize=10, res=300) +plot_bivariate(my_points) +dev.off() + +} From ba3ac023a310df7309d64e9322c6ffc84c5fc196 Mon Sep 17 00:00:00 2001 From: cmahony Date: Mon, 12 Feb 2024 15:16:31 -0800 Subject: [PATCH 06/17] add plotly and stinepack to DESCRIPTION --- DESCRIPTION | 2 ++ 1 file changed, 2 insertions(+) diff --git a/DESCRIPTION b/DESCRIPTION index 1391680a..1edaf54c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -37,8 +37,10 @@ Imports: uuid Suggests: parallel, + plotly, rmarkdown, remotes, + stinepack, testthat (>= 3.0.0), withr Depends: From ad9057233be29ff6780badbcced2c776d957def2 Mon Sep 17 00:00:00 2001 From: cmahony Date: Mon, 12 Feb 2024 15:41:33 -0800 Subject: [PATCH 07/17] bump and news --- DESCRIPTION | 4 ++-- NEWS.md | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 1edaf54c..f288d55b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: climr Title: Downscaling of global climate data -Version: 0.0.1.9003 -Date: 05-02-2024 +Version: 0.0.1.9004 +Date: 12-02-2024 Authors@R: c( person("Kiri","Daust", email = "kiri.daust@gov.bc.ca", role = "aut"), person("Colin", "Mahony", email = "Colin.Mahony@gov.bc.ca", role = c("aut", "cre"), diff --git a/NEWS.md b/NEWS.md index 1f783651..352d7a2d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,7 @@ * code further streamlined * new messages warn user about meaningless `downscale`/`climr_downscale` argument combinations * argument options in `climr_downscale(..., which_normal)` now match the options of `normal_input(..., normal)` +* add `plot_bivariate()` function to generate plots showing climate model ensemble variation in recent and future climate change. ## Behaviour changes * `xyz` (argument to `climr_downscale` and `downscale`) and `in_xyz` (argument to `get_bb`), must now be a 4 column `data.table` (or coercible class) with `lon`, `lat`, `elev` and `id` columns. All other columns are ignored and NOT returned. Column order no longer matters. From 233d98fd1cce1eab712bfa4ad062fa2cc63b0613 Mon Sep 17 00:00:00 2001 From: cmahony Date: Mon, 12 Feb 2024 17:26:18 -0800 Subject: [PATCH 08/17] final edits --- R/plot_bivariate.R | 60 +++++++++++++++++++++++-------------------- man/downscaling.Rd | 4 +-- man/plot_bivariate.Rd | 14 +++++----- 3 files changed, 42 insertions(+), 36 deletions(-) diff --git a/R/plot_bivariate.R b/R/plot_bivariate.R index 9eca4eed..6252ccd2 100644 --- a/R/plot_bivariate.R +++ b/R/plot_bivariate.R @@ -1,7 +1,14 @@ #' Bivariate climate change plots #' #' @description -#' Bivariate plots of 21st century climate change for user-selected locations and climate variables. +#' Bivariate plots of 21st century climate change for user-selected locations and climate variables. +#' The purposes of the plot are to +#' \enumerate{ +#' \item show differences in climate change trends among global climate models (GCMs), +#' \item show the differences between multiple simulations of each model, and +#' \item compare simulated climate change to observed climate change in the 2001-2020 period. +#' } +#' All climate changes are relative to the 1961-1990 reference period normals. #' #' @details #' The input table `xyz` can be a single location or multiple locations. If multiple locations, the plot provides the mean of the anomalies for these locations. @@ -23,7 +30,7 @@ #' @return NULL. Draws a plot in the active graphics device. #' #' @examples -#' # data frame of arbitrary points on vancouver island +#' # data frame of arbitrary points on Vancouver Island #' my_points <- data.frame(lon = c(-123.4404, -123.5064, -124.2317), #' lat = c(48.52631, 48.46807, 49.21999), #' elev = c(52, 103, 357), @@ -46,7 +53,6 @@ plot_bivariate <- function( yvar = "PPT_sm", # percent_x = NULL, TODO: set up an override for ratio variables being expressed as percent anomalies # percent_y = NULL, TODO: set up an override for ratio variables being expressed as percent anomalies - which_normal = "auto", # TODO: remove once the bug is fixed. period_focal = list_gcm_period()[1], gcm_models = list_gcm()[c(1,4,5,6,7,10,11,12)], ssp = list_ssp()[2], @@ -72,13 +78,11 @@ plot_bivariate <- function( xvar_type <- variables$Scale[which(variables$Code==xvar)] yvar_type <- variables$Scale[which(variables$Code==yvar)] - colors = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = TRUE)][-1] - set.seed(2) - ColScheme <- sample(colors,length(gcm_models)) # TODO select better colors + colors <- c("#A6CEE3", "#1F78B4", "#B2DF8A", "#33A02C", "#FB9A99", "#E31A1C", "#FDBF6F", "#FF7F00", "#CAB2D6", "#6A3D9A", "#1e90ff", "#B15928", "#FFFF99") + ColScheme <- colors[1:length(gcm_models)] # generate the climate data data <- climr_downscale(xyz, - which_normal = which_normal, # TODO: remove once the bug is fixed. historic_period = historic_period, gcm_models = gcm_models, ssp = ssp, @@ -93,7 +97,9 @@ plot_bivariate <- function( data[, yanom := if(yvar_type=="Log") (get(yvar)/get(yvar)[1]-1) else (get(yvar) - get(yvar)[1]), by = id] # collapse the points down to a mean anomaly - data <- data[, .(xanom = mean(xanom), yanom = mean(yanom)), by = .(GCM, SSP, RUN, PERIOD)] # TODO: need to give special treatment to SSPs; perhaps by averaging the ensemble mean but retaining the individual runs + data.all <- copy(data[, .(xanom = mean(xanom), yanom = mean(yanom)), by = .(GCM, SSP, RUN, PERIOD)]) + # collapse the SSP field to calculate a single-model ensemble mean + data <- copy(data.all[, .(xanom = mean(xanom), yanom = mean(yanom)), by = .(GCM, RUN, PERIOD)]) # ensemble mean for the selected period ensMean <- data[!is.na(GCM) & RUN == "ensembleMean" & PERIOD == period_focal, .(xanom = mean(xanom), yanom = mean(yanom)), ] # observed climate @@ -104,7 +110,7 @@ plot_bivariate <- function( # BASE PLOT # initiate the plot par(mar=c(3,4,0,1), mgp=c(1.25, 0.25,0), cex=1.5) - plot(data$xanom,data$yanom,col="white", tck=0, xaxt="n", yaxt="n", ylab="", + plot(data.all$xanom,data.all$yanom,col="white", tck=0, xaxt="n", yaxt="n", ylab="", xlab=paste("Change in", variables$Variable[which(variables$Code==xvar)]) ) par(mgp=c(2.5,0.25, 0)) @@ -120,8 +126,8 @@ plot_bivariate <- function( if(show_runs){ for(gcm in gcm_models){ i=which(gcm_models==gcm) - x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] - y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] + x.runs <- data.all[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] + y.runs <- data.all[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] points(x.runs,y.runs, pch=21, bg=ColScheme[i], cex=1) } } @@ -135,7 +141,7 @@ plot_bivariate <- function( if(length(unique(sign(diff(x2))))==1){ x3 <- if(unique(sign(diff(x2)))==-1) rev(x2) else x2 y3 <- if(unique(sign(diff(x2)))==-1) rev(y2) else y2 - s <- stinepack::stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data$xanom))/500)) + s <- stinepack::stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data.all$xanom))/500)) lines(s, col=ColScheme[i], lwd=2) } else lines(x2, y2, col=ColScheme[i], lwd=2) points(x2,y2, pch=16, col=ColScheme[i], cex=0.5) @@ -146,8 +152,8 @@ plot_bivariate <- function( } # axis labels - axis(1, at=pretty(data$xanom), labels=if(xvar_type=="Log") paste(pretty(data$xanom)*100, "%", sep="") else pretty(data$xanom), tck=0) - axis(2, at=pretty(data$yanom), labels=if(yvar_type=="Log") paste(pretty(data$yanom)*100, "%", sep="") else pretty(data$yanom), las=2, tck=0) + axis(1, at=pretty(data.all$xanom), labels=if(xvar_type=="Log") paste(pretty(data.all$xanom)*100, "%", sep="") else pretty(data.all$xanom), tck=0) + axis(2, at=pretty(data.all$yanom), labels=if(yvar_type=="Log") paste(pretty(data.all$yanom)*100, "%", sep="") else pretty(data.all$yanom), las=2, tck=0) # Legend s <- c(show_observed, show_runs, TRUE, show_trajectories, show_ensMean) @@ -169,19 +175,18 @@ plot_bivariate <- function( } else { # PLOTLY PLOT - # TODO: the colors aren't plotting correctly. need to fix this. - + #initiate the plot - fig <- plot_ly(x=data$xanom,y=data$yanom, type = 'scatter', mode = 'markers', marker = list(color ="lightgray", size=5), hoverinfo="none", color="All models/scenarios/runs/periods") + fig <- plot_ly(x=data.all$xanom,y=data.all$yanom, type = 'scatter', mode = 'markers', marker = list(color ="lightgray", size=5), hoverinfo="none", color="All models/scenarios/runs/periods") # axis titles - fig <- fig %>% layout(xaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==xvar)]), range=range(data$xanom)), - yaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==yvar)]), range=range(data$yanom)) + fig <- fig %>% layout(xaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==xvar)]), range=range(data.all$xanom)), + yaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==yvar)]), range=range(data.all$yanom)) ) # observed climate fig <- fig %>% add_markers(obs$xanom ,obs$yanom, name="Observed Climate (2001-2020)", text="observed\n(2001-2020)", hoverinfo="text", - marker = list(size = 25, color = "grey"), symbol = 43) + marker = list(size = 25, color = "grey", symbol = 1)) # ensemble mean fig <- fig %>% add_markers(ensMean$xanom,ensMean$yanom, name="Ensemble mean", text="Ensemble mean", hoverinfo="text", @@ -191,10 +196,11 @@ plot_bivariate <- function( if(show_runs){ for(gcm in gcm_models){ i=which(gcm_models==gcm) - x.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] - y.runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] - runs <- data[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, RUN] - fig <- fig %>% add_markers(x=x.runs,y=y.runs, color = ColScheme[i], name="Individual GCM runs", text=paste(gcm_models[i], runs), hoverinfo="text", showlegend = FALSE, + x.runs <- data.all[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, xanom] + y.runs <- data.all[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] + runs <- data.all[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, RUN] + ssps <- data.all[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, SSP] + fig <- fig %>% add_markers(x=x.runs,y=y.runs, color = ColScheme[i], name="Individual GCM runs", text=paste(gcm_models[i], ssps, runs), hoverinfo="text", showlegend = if(i==1) TRUE else FALSE, marker = list(size = 7, color = ColScheme[i], line = list(color = "black", width = 1)), legendgroup=paste("group", i, sep="")) } } @@ -210,7 +216,7 @@ plot_bivariate <- function( if(length(unique(sign(diff(x2))))==1){ x3 <- if(unique(sign(diff(x2)))==-1) rev(x2) else x2 y3 <- if(unique(sign(diff(x2)))==-1) rev(y2) else y2 - s <- stinepack::stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data$xanom))/500)) # way better than interpSpline, not prone to oscillations + s <- stinepack::stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data.all$xanom))/500)) # way better than interpSpline, not prone to oscillations fig <- fig %>% add_trace(x=s$x, y=s$y, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) } else { fig <- fig %>% add_trace(x=x2, y=y2, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) @@ -219,7 +225,7 @@ plot_bivariate <- function( marker = list(size = 8, color = ColScheme[i]), legendgroup=paste("group", i, sep=""), showlegend = FALSE) } j=which(list_gcm_period()==period_focal)+1 - fig <- fig %>% add_markers(x2[j],y2[j], color = gcm_models[i], colors=ColScheme[i], text=gcm_models[i], + fig <- fig %>% add_markers(x2[j],y2[j], color = gcm_models[i], colors=ColScheme[i], marker = list(size = 20, color = ColScheme[i], line = list(color = "black", width = 1)), legendgroup=paste("group", i, sep="")) @@ -235,5 +241,3 @@ plot_bivariate <- function( } } } - - diff --git a/man/downscaling.Rd b/man/downscaling.Rd index b78b0da5..8203e8c1 100644 --- a/man/downscaling.Rd +++ b/man/downscaling.Rd @@ -7,7 +7,7 @@ \usage{ climr_downscale( xyz, - which_normal = c("auto", list_normal()), + which_normal = "auto", historic_period = NULL, historic_ts = NULL, gcm_models = NULL, @@ -45,7 +45,7 @@ and a unique "id". Any extra columns will be ignored and not output.} \item{which_normal}{character. Which normal layer to use. Default is "auto", which selects the highest resolution normal for each point. -Other options are one of \code{\link[=list_normal]{list_normal()}}} +Other options are one of \code{\link[=list_normal]{list_normal()}}.} \item{historic_period}{character. Which historic period. Default \code{NULL}} diff --git a/man/plot_bivariate.Rd b/man/plot_bivariate.Rd index 9d402f80..e7444bf8 100644 --- a/man/plot_bivariate.Rd +++ b/man/plot_bivariate.Rd @@ -8,7 +8,6 @@ plot_bivariate( xyz, xvar = "Tave_sm", yvar = "PPT_sm", - which_normal = "auto", period_focal = list_gcm_period()[1], gcm_models = list_gcm()[c(1, 4, 5, 6, 7, 10, 11, 12)], ssp = list_ssp()[2], @@ -32,10 +31,6 @@ and a unique "id". Any extra columns will be ignored and not output.} \item{yvar}{character. y-axis variable. options are \code{list_variables()}.} -\item{which_normal}{character. Which normal layer to use. -Default is "auto", which selects the highest resolution normal for each point. -Other options are one of \code{\link[=list_normal]{list_normal()}}} - \item{period_focal}{character. The 20-year period for which to plot the ensemble detail. options are \code{list_gcm_period()}.} \item{gcm_models}{character. Vector of GCM names. Options are \code{\link[=list_gcm]{list_gcm()}}. Used for gcm periods, gcm timeseries, and historic timeseries. Default \code{NULL}} @@ -70,6 +65,13 @@ NULL. Draws a plot in the active graphics device. } \description{ Bivariate plots of 21st century climate change for user-selected locations and climate variables. +The purposes of the plot are to +\enumerate{ +\item show differences in climate change trends among global climate models (GCMs), +\item show the differences between multiple simulations of each model, and +\item compare simulated climate change to observed climate change in the 2001-2020 period. +} +All climate changes are relative to the 1961-1990 reference period normals. } \details{ The input table \code{xyz} can be a single location or multiple locations. If multiple locations, the plot provides the mean of the anomalies for these locations. @@ -77,7 +79,7 @@ The climate change trajectories provided by \code{show_trajectories} are points This plot is designed to be used with a single SSP scenario. If multiple scenarios are passed to the plot, the GCM means and ensemble mean are averaged across the scenarios, but the individual runs for all scenarios are plotted separately. } \examples{ -# data frame of arbitrary points on vancouver island +# data frame of arbitrary points on Vancouver Island my_points <- data.frame(lon = c(-123.4404, -123.5064, -124.2317), lat = c(48.52631, 48.46807, 49.21999), elev = c(52, 103, 357), From abe32e166d4f7a7b4b6199ab48ed3411868b5682 Mon Sep 17 00:00:00 2001 From: CeresBarros Date: Tue, 13 Feb 2024 15:34:45 -0800 Subject: [PATCH 09/17] fix bad merge --- DESCRIPTION | 2 -- 1 file changed, 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 58cbf6e2..429712c3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,8 +2,6 @@ Package: climr Title: Downscaling of global climate data Version: 0.0.2.9004 Date: 12-02-2024 -Version: 0.0.2.9004 -Date: 12-02-2024 Authors@R: c( person("Kiri","Daust", email = "kiri.daust@gov.bc.ca", role = "aut"), person("Colin", "Mahony", email = "Colin.Mahony@gov.bc.ca", role = c("aut", "cre"), From 9d94469325e826592c724e8bcdabb6c93ed51b59 Mon Sep 17 00:00:00 2001 From: CeresBarros Date: Tue, 13 Feb 2024 15:41:18 -0800 Subject: [PATCH 10/17] install XQuartx in macOS GHA nodes --- .github/workflows/R-CMD-check.yaml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 5b5b2c55..53f7f32f 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -39,9 +39,15 @@ jobs: R_REPRODUCIBLE_RUN_ALL_TESTS: true steps: + - name: Install X11 dependencies on MacOS + if: runner.os == 'macOS' + run: | + brew install --cask xquartz + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 + - uses: r-lib/actions/setup-r@v2 with: From 7de92a89de4d545788366a3e5737502cdaa0a3a2 Mon Sep 17 00:00:00 2001 From: CeresBarros Date: Tue, 13 Feb 2024 16:03:02 -0800 Subject: [PATCH 11/17] add missing global vars --- R/globalVars.R | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/R/globalVars.R b/R/globalVars.R index f5e517bc..8a70f99d 100644 --- a/R/globalVars.R +++ b/R/globalVars.R @@ -1,6 +1,8 @@ utils::globalVariables(c( - "..cols", "..expectedCols", "fullnm", "id", "id_orig", "lat", "lon", + "..cols", "..expectedCols", "fullnm", "GCM", "id", "id_orig", "lat", "lon", "laynum", "mod", "numlay", "period", - "Period", "Period1", "Period2", "run", "Run", - "scenario", "Scenario", "skip_if_not_installed", "var", "Year", "." + "Period", "PERIOD", "Period1", "Period2", "run", "Run", "RUN", + "scenario", "Scenario", "skip_if_not_installed", "SSP", + "xanom", "var", + "variables", "yanom", "Year", ".", "%>%" )) From 3412761654da0c2ff1a83be69b3ddfedc9bc13b0 Mon Sep 17 00:00:00 2001 From: CeresBarros Date: Tue, 13 Feb 2024 16:03:18 -0800 Subject: [PATCH 12/17] added missing imports --- R/plot_bivariate.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/plot_bivariate.R b/R/plot_bivariate.R index 6252ccd2..642c8cbf 100644 --- a/R/plot_bivariate.R +++ b/R/plot_bivariate.R @@ -28,6 +28,8 @@ #' @param interactive logical. If TRUE, an interactive plot is generated using `{plotly}`. If FALSE, a plot is generated using base graphics. #' #' @return NULL. Draws a plot in the active graphics device. +#' +#' @importFrom graphics axis legend lines par points text title #' #' @examples #' # data frame of arbitrary points on Vancouver Island From 267d2ac9c85aa022d32b42f9f7e2e28e9c92c77e Mon Sep 17 00:00:00 2001 From: CeresBarros Date: Tue, 13 Feb 2024 16:03:55 -0800 Subject: [PATCH 13/17] follow CRAN convention: don't load plotly --- R/plot_bivariate.R | 44 ++++++++++++++++++++++---------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/R/plot_bivariate.R b/R/plot_bivariate.R index 642c8cbf..41074120 100644 --- a/R/plot_bivariate.R +++ b/R/plot_bivariate.R @@ -172,27 +172,27 @@ plot_bivariate <- function( } else { - if(!require(plotly)){ + if(!requireNamespace("plotly")){ stop("package 'plotly' must be installed when 'interactive==TRUE'") } else { # PLOTLY PLOT - + #initiate the plot - fig <- plot_ly(x=data.all$xanom,y=data.all$yanom, type = 'scatter', mode = 'markers', marker = list(color ="lightgray", size=5), hoverinfo="none", color="All models/scenarios/runs/periods") + fig <- plotly::plot_ly(x=data.all$xanom,y=data.all$yanom, type = 'scatter', mode = 'markers', marker = list(color ="lightgray", size=5), hoverinfo="none", color="All models/scenarios/runs/periods") # axis titles - fig <- fig %>% layout(xaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==xvar)]), range=range(data.all$xanom)), - yaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==yvar)]), range=range(data.all$yanom)) + fig <- fig %>% plotly::layout(xaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==xvar)]), range=range(data.all$xanom)), + yaxis = list(title=paste("Change in", variables$Variable[which(variables$Code==yvar)]), range=range(data.all$yanom)) ) # observed climate - fig <- fig %>% add_markers(obs$xanom ,obs$yanom, name="Observed Climate (2001-2020)", text="observed\n(2001-2020)", hoverinfo="text", - marker = list(size = 25, color = "grey", symbol = 1)) + fig <- fig %>% plotly::add_markers(obs$xanom ,obs$yanom, name="Observed Climate (2001-2020)", text="observed\n(2001-2020)", hoverinfo="text", + marker = list(size = 25, color = "grey", symbol = 1)) # ensemble mean - fig <- fig %>% add_markers(ensMean$xanom,ensMean$yanom, name="Ensemble mean", text="Ensemble mean", hoverinfo="text", - marker = list(size = 20, color = "grey", symbol = 3)) + fig <- fig %>% plotly::add_markers(ensMean$xanom,ensMean$yanom, name="Ensemble mean", text="Ensemble mean", hoverinfo="text", + marker = list(size = 20, color = "grey", symbol = 3)) # plot individual runs if(show_runs){ @@ -202,8 +202,8 @@ plot_bivariate <- function( y.runs <- data.all[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, yanom] runs <- data.all[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, RUN] ssps <- data.all[GCM==gcm & RUN != "ensembleMean" & PERIOD == period_focal, SSP] - fig <- fig %>% add_markers(x=x.runs,y=y.runs, color = ColScheme[i], name="Individual GCM runs", text=paste(gcm_models[i], ssps, runs), hoverinfo="text", showlegend = if(i==1) TRUE else FALSE, - marker = list(size = 7, color = ColScheme[i], line = list(color = "black", width = 1)), legendgroup=paste("group", i, sep="")) + fig <- fig %>% plotly::add_markers(x=x.runs,y=y.runs, color = ColScheme[i], name="Individual GCM runs", text=paste(gcm_models[i], ssps, runs), hoverinfo="text", showlegend = if(i==1) TRUE else FALSE, + marker = list(size = 7, color = ColScheme[i], line = list(color = "black", width = 1)), legendgroup=paste("group", i, sep="")) } } @@ -219,24 +219,24 @@ plot_bivariate <- function( x3 <- if(unique(sign(diff(x2)))==-1) rev(x2) else x2 y3 <- if(unique(sign(diff(x2)))==-1) rev(y2) else y2 s <- stinepack::stinterp(x3,y3, seq(min(x3),max(x3), diff(range(data.all$xanom))/500)) # way better than interpSpline, not prone to oscillations - fig <- fig %>% add_trace(x=s$x, y=s$y, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) + fig <- fig %>% plotly::add_trace(x=s$x, y=s$y, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) } else { - fig <- fig %>% add_trace(x=x2, y=y2, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) + fig <- fig %>% plotly::add_trace(x=x2, y=y2, color = ColScheme[i], type = 'scatter', mode = 'lines', line = list(color=ColScheme[i], width = 2), marker=NULL, legendgroup=paste("group", i, sep=""), showlegend = FALSE) } - fig <- fig %>% add_markers(x=x2,y=y2, color = ColScheme[i], text=gcm_models[i], hoverinfo="text", - marker = list(size = 8, color = ColScheme[i]), legendgroup=paste("group", i, sep=""), showlegend = FALSE) + fig <- fig %>% plotly::add_markers(x=x2,y=y2, color = ColScheme[i], text=gcm_models[i], hoverinfo="text", + marker = list(size = 8, color = ColScheme[i]), legendgroup=paste("group", i, sep=""), showlegend = FALSE) } j=which(list_gcm_period()==period_focal)+1 - fig <- fig %>% add_markers(x2[j],y2[j], color = gcm_models[i], colors=ColScheme[i], - marker = list(size = 20, color = ColScheme[i], line = list(color = "black", width = 1)), - legendgroup=paste("group", i, sep="")) + fig <- fig %>% plotly::add_markers(x2[j],y2[j], color = gcm_models[i], colors=ColScheme[i], + marker = list(size = 20, color = ColScheme[i], line = list(color = "black", width = 1)), + legendgroup=paste("group", i, sep="")) - fig <- fig %>% add_annotations(x=x2[j],y=y2[j], text = sprintf("%s", substr(gcm_models, 1, 2)[i]), xanchor = 'center', yanchor = 'center', showarrow = FALSE, - legendgroup=paste("group", i, sep="")) + fig <- fig %>% plotly::add_annotations(x=x2[j],y=y2[j], text = sprintf("%s", substr(gcm_models, 1, 2)[i]), xanchor = 'center', yanchor = 'center', showarrow = FALSE, + legendgroup=paste("group", i, sep="")) } - if(xvar_type=="Log") fig <- fig %>% layout(xaxis = list(tickformat = "%")) - if(yvar_type=="Log") fig <- fig %>% layout(yaxis = list(tickformat = "%")) + if(xvar_type=="Log") fig <- fig %>% plotly::layout(xaxis = list(tickformat = "%")) + if(yvar_type=="Log") fig <- fig %>% plotly::layout(yaxis = list(tickformat = "%")) fig } From dd8727028689490d7eebb6cbc7aa51c3d27f612c Mon Sep 17 00:00:00 2001 From: CeresBarros Date: Tue, 13 Feb 2024 16:04:20 -0800 Subject: [PATCH 14/17] follow CRAN convention: use data() in current env --- R/plot_bivariate.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/plot_bivariate.R b/R/plot_bivariate.R index 41074120..c735aa04 100644 --- a/R/plot_bivariate.R +++ b/R/plot_bivariate.R @@ -74,7 +74,7 @@ plot_bivariate <- function( stop("package stinepack must be installed to use this function") } else { - data("variables") + data("variables", envir = environment()) # variable types for default scaling (percent or absolute) xvar_type <- variables$Scale[which(variables$Code==xvar)] From 02d257f6cd869b4c9cf3ba4e49522436153c153f Mon Sep 17 00:00:00 2001 From: CeresBarros Date: Tue, 13 Feb 2024 16:04:32 -0800 Subject: [PATCH 15/17] re-doc --- NAMESPACE | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/NAMESPACE b/NAMESPACE index 7671130c..fd1be1f0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -45,6 +45,13 @@ importFrom(data.table,setkey) importFrom(data.table,setorder) importFrom(grDevices,hcl.colors) importFrom(grDevices,palette) +importFrom(graphics,axis) +importFrom(graphics,legend) +importFrom(graphics,lines) +importFrom(graphics,par) +importFrom(graphics,points) +importFrom(graphics,text) +importFrom(graphics,title) importFrom(methods,is) importFrom(pool,dbPool) importFrom(pool,poolClose) From 3b9ed56ad220ced47a9415bca7d88b2db94379b4 Mon Sep 17 00:00:00 2001 From: CeresBarros Date: Tue, 13 Feb 2024 16:04:56 -0800 Subject: [PATCH 16/17] added ploty/stinepack to GHA set up --- .github/workflows/R-CMD-check.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 53f7f32f..c8fafd96 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -58,9 +58,11 @@ jobs: - uses: r-lib/actions/setup-r-dependencies@v2 with: extra-packages: | + any::plotly any::rcmdcheck any::Rcpp any::remotes + any::stinepack any::withr - uses: r-lib/actions/check-r-package@v2 From 54f928e14fcc13a801d9a1bd88540693cf606665 Mon Sep 17 00:00:00 2001 From: CeresBarros Date: Tue, 13 Feb 2024 16:05:38 -0800 Subject: [PATCH 17/17] added vignette builder to DESCRIPTION --- DESCRIPTION | 2 ++ 1 file changed, 2 insertions(+) diff --git a/DESCRIPTION b/DESCRIPTION index 429712c3..8f71abf5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -42,8 +42,10 @@ Suggests: remotes, stinepack, testthat (>= 3.0.0), + utils, withr Depends: R (>= 4.0) Config/testthat/edition: 3 LazyData: true +VignetteBuilder: utils