Nothing
#
# compare_two_runs_aam.R
#
#' Create a tornado plot of the differences in annual averaged masses of variable in two model runs.
#'
#' Reads final year annual average mass results generated by two e2e_run() model runs with different inputs, and plot the difference
#' between the two as tornado diagrams for the water column and seabed components separately.
#'
#' Annual average mass results from the final year of runs of the e2e_run() function are held in memory and also stored as CSV files as
#' /results/Modelname/Variantname/ZONE_model_anav_biomass-*.csv
#' where ZONE is INSHORE, OFFSHORE or WHOLEDOMAIN and * represents the model run identifier (model.ident) text embedded in the R-list object created by the e2e_read() function.
#'
#' Optionally, the function can read these data from files created by past model runs, or use data still held as results objects from recent runs.
#'
#' The tornado plots show differences between the runs as horizontal barplots. If the first set of model results is regarded as
#' a baseline run, and the second as a scenario run (e.g. with different driving data), then scenario values > baseline values are
#' shown by green bars to the right of 'no difference' (i.e scenario=baseline). Conversely scenario values < baseline values
#' are shown by red bars to the left.
#'
#' The bars can be plotted on a log10 scale or a percentage scale governed by argument settings.
#'
#' The function returns a list object comprising two dataframes (changewater and changeseabed). The first if the proportional difference in water column variables between the scenario run and the baseline
#' the second is the proportional difference in seabed variables. The first column in each dataframe is the proportional difference expressed on a log10 scale, the second column as a percentage.
#'
#' @param model1 R-list object defining the model configuration compiled by the e2e_read() function for the first (baseline) run. Set to NA if argument 'from.csv1'=FALSE. Default = NA
#' @param from.csv1 TRUE = read baseline data from csv files created by a past model run. FALSE = use baseline data held in memory. Default = FALSE.
#' @param results1 R-list object containing the results from the first (baseline) run of the e2e_run() function. Set to NA if argument 'from.csv1'=TRUE. Default = value expected.
#' @param model2 R-list object defining the model configuration compiled by the e2e_read() function for the second (scenario) run. Set to NA if argument 'from.csv2'=FALSE. Default = NA
#' @param from.csv2 TRUE = read scenario data from csv files created by a past model run. FALSE = use scenario data held in memory. Default = FALSE.
#' @param results2 R-list object containing the results from the second (scenario) run of the e2e_run() function. Set to NA if argument 'from.csv1'=TRUE. Default = value expected.
#' @param log.pc Value="LG" for data to be plotted on a log10 scale, value = "PC" for data to be plotted on a percentage difference scale. Default = "PC".
#' @param zone Value = "O" for offshore, "I" for inshore, or "W" for whole model domain (all upper case). Default = "W".
#' @param bpmin Axis minimum for plot - i.e. the maximum NEGATIVE value of (scenario-baseline). Default = -50, given log.pc="PC" (percentage differences). Needs to be reset to e.g. -0.3 if log.pc="LG" (log scale).
#' @param bpmax Axis maximum for plot - i.e. the maximum POSITIVE value of (scenario-baseline). Default = +50, given log.pc="PC" (percentage differences). Needs to be reset to e.g. +0.3 if log.pc="LG" (log scale).
#' @param maintitle A descriptive text field (in quotes) to be added above the plot. Keep to 45 characters including spaces. Default="".
#'
#' @return List object comprising two dataframes of the displayed data, graphical display in a new graphics window.
#'
#' @importFrom graphics text
#'
#' @noRd
#
# ---------------------------------------------------------------------
# | |
# | Authors: Mike Heath |
# | Department of Mathematics and Statistics |
# | University of Strathclyde, Glasgow |
# | |
# | Date of this version: December 2019 |
# | |
# ---------------------------------------------------------------------
compare_two_runs_aam <- function(model1=NA,from.csv1=FALSE,results1,
model2=NA,from.csv2=FALSE,results2,
log.pc="PC",
zone="W",
bpmin=(-50),bpmax=(+50),
maintitle="" ) {
start_par = par()$mfrow
on.exit(par(mfrow = start_par))
if((zone %in% c("O","I","W"))==FALSE) {
stop("Retry with argument zone = O, I or W (upper case in double quotes)! \n")
}
if((log.pc %in% c("LG","PC"))==FALSE) {
stop("Retry with argument log.pc = LG or PC (upper case in double quotes)! \n")
}
if(bpmin>0) {
stop("bpmin argument should be a negative value. Try again! \n")
}
if(bpmax<0) {
stop("bpmax argument should be a positive value. Try again! \n")
}
if(log.pc=="LG" & (bpmin<(-2) | bpmax > 2)==TRUE) {
message("bpmin and/or bpmax arguments seem very large given log.pc = LG. Are you sure!")
}
if(log.pc=="PC" & (bpmin>(-1) | bpmax<1)==TRUE) {
message("bpmin and/or bpmax arguments seem very small given log.pc = PC. Are you sure!")
}
if(from.csv1==FALSE) {
#baseline run data
print("Using baseline data held in memory from an existing model run")
final.year.outputs1 <- elt(results1, "final.year.outputs")
if(zone=="I") baselinedata <- elt(final.year.outputs1,"mass_results_inshore")
if(zone=="O") baselinedata <- elt(final.year.outputs1,"mass_results_offshore")
if(zone=="W") baselinedata <- elt(final.year.outputs1,"mass_results_wholedomain")
}
if(from.csv2==FALSE) {
#scenario run data
print("Using scenario data held in memory from an existing model run")
final.year.outputs2 <- elt(results2, "final.year.outputs")
if(zone=="I") scenariodata <- elt(final.year.outputs2,"mass_results_inshore")
if(zone=="O") scenariodata <- elt(final.year.outputs2,"mass_results_offshore")
if(zone=="W") scenariodata <- elt(final.year.outputs2,"mass_results_wholedomain")
}
if(from.csv1==TRUE) {
#baseline run data
print("Using baseline data held in a csv files from a past model run")
resultsdir1 <- elt(model1, "setup", "resultsdir")
model.ident1 <- elt(model1, "setup", "model.ident")
# model.name1 <- elt(model1, "setup", "model.name")
# model.variant1 <- elt(model1, "setup", "model.variant")
if(zone=="I") basefile <- csvname(resultsdir1, "INSHORE_model_anav_biomass", model.ident1)
if(zone=="O") basefile <- csvname(resultsdir1, "OFFSHORE_model_anav_biomass", model.ident1)
if(zone=="W") basefile <- csvname(resultsdir1, "WHOLEDOMAIN_model_anav_biomass", model.ident1)
check.exists(basefile)
baselinedata<-readcsv(basefile)
}
if(from.csv2==TRUE) {
#scenario run data
print("Using scenario data held in a csv files from a past model run")
resultsdir2 <- elt(model2, "setup", "resultsdir")
model.ident2 <- elt(model2, "setup", "model.ident")
# model.name2 <- elt(model2, "setup", "model.name")
# model.variant2 <- elt(model2, "setup", "model.variant")
if(zone=="I") scenfile <- csvname(resultsdir2, "INSHORE_model_anav_biomass", model.ident2)
if(zone=="O") scenfile <- csvname(resultsdir2, "OFFSHORE_model_anav_biomass", model.ident2)
if(zone=="W") scenfile <- csvname(resultsdir2, "WHOLEDOMAIN_model_anav_biomass", model.ident2)
check.exists(scenfile)
scenariodata<-readcsv(scenfile)
}
baselinewater <- data.frame(rep(0,15))
scenariowater <- data.frame(rep(0,15))
baselineseabed <- data.frame(rep(0,8))
scenarioseabed <- data.frame(rep(0,8))
baselinewater[1,1]<-sum(baselinedata[1:2,1],na.rm=TRUE) #wc detritus
baselinewater[2,1]<-sum(baselinedata[8:9,1],na.rm=TRUE) #wc ammonia
baselinewater[3,1]<-sum(baselinedata[11:12,1],na.rm=TRUE) #wc nitrate
baselinewater[4,1]<-(baselinedata[14,1]) #macrophytes
baselinewater[5,1]<-sum(baselinedata[15:16,1],na.rm=TRUE) #phyt
baselinewater[6,1]<-(baselinedata[19,1]) #sdb larvae
baselinewater[7,1]<-(baselinedata[17,1] ) #omnizoo
baselinewater[8,1]<-(baselinedata[21,1]) #csb larvae
baselinewater[9,1]<-(baselinedata[18,1]) #carnzoo
baselinewater[10,1]<-sum(baselinedata[23:24,1],na.rm=TRUE) #pelfish settled + larvae
baselinewater[11,1]<-(baselinedata[25,1]) #migfish
baselinewater[12,1]<-sum(baselinedata[26:27,1],na.rm=TRUE) #demfish settled + larvae
baselinewater[13,1]<-(baselinedata[28,1] ) #birds
baselinewater[14,1]<-(baselinedata[29,1] ) #pinnipeds
baselinewater[15,1]<-(baselinedata[30,1] ) #cetaceans
baselineseabed[1,1]<-(baselinedata[5,1]) #discards
baselineseabed[2,1]<-(baselinedata[6,1]) #corpses
baselineseabed[3,1]<-(baselinedata[3,1]) #sed labile + refractory det
baselineseabed[4,1]<-(baselinedata[7,1]) #acrophyte debris
baselineseabed[5,1]<-(baselinedata[10,1]) #sed ammonia
baselineseabed[6,1]<-(baselinedata[13,1]) #sed nitrate
baselineseabed[7,1]<-(baselinedata[20,1]) #sdb settled
baselineseabed[8,1]<-(baselinedata[22,1]) #csb settled
scenariowater[1,1]<-sum(scenariodata[1:2,1],na.rm=TRUE) #wc detritus
scenariowater[2,1]<-sum(scenariodata[8:9,1],na.rm=TRUE) #wc ammonia
scenariowater[3,1]<-sum(scenariodata[11:12,1],na.rm=TRUE) #wc nitrate
scenariowater[4,1]<-(scenariodata[14,1]) #macrophytes
scenariowater[5,1]<-sum(scenariodata[15:16,1],na.rm=TRUE) #phyt
scenariowater[6,1]<-(scenariodata[19,1]) #sdb larvae
scenariowater[7,1]<-(scenariodata[17,1] ) #omnizoo
scenariowater[8,1]<-(scenariodata[21,1]) #csb larvae
scenariowater[9,1]<-(scenariodata[18,1]) #carnzoo
scenariowater[10,1]<-sum(scenariodata[23:24,1],na.rm=TRUE) #pelfish settled + larvae
scenariowater[11,1]<-(scenariodata[25,1]) #migfish
scenariowater[12,1]<-sum(scenariodata[26:27,1],na.rm=TRUE) #demfish settled + larvae
scenariowater[13,1]<-(scenariodata[28,1] ) #birds
scenariowater[14,1]<-(scenariodata[29,1] ) #pinnipeds
scenariowater[15,1]<-(scenariodata[30,1] ) #cetaceans
scenarioseabed[1,1]<-(scenariodata[5,1]) #discards
scenarioseabed[2,1]<-(scenariodata[6,1]) #corpses
scenarioseabed[3,1]<-(scenariodata[3,1]) #sed labile + refractory det
scenarioseabed[4,1]<-(scenariodata[7,1]) #macrophyte debris
scenarioseabed[5,1]<-(scenariodata[10,1]) #sed ammonia
scenarioseabed[6,1]<-(scenariodata[13,1]) #sed nitrate
scenarioseabed[7,1]<-(scenariodata[20,1]) #sdb settled
scenarioseabed[8,1]<-(scenariodata[22,1]) #csb settled
rownames(baselinewater)<-c("Suspended bacteria & detritus",
"Water column ammonia","Water column nitrate",
"Macrophytes",
"Phytoplankton",
"Benthic susp/dep feeder larvae",
"Omnivorous zooplankton",
"Benthic carn/scav feeder larvae",
"Carnivorous zooplankton",
"Plantiv. fish adults+larvae",
"Migratory fish",
"Demersal fish adults+larvae",
"Birds","Pinnipeds","Cetaceans")
rownames(baselineseabed)<-c("Fishery discards & offal",
"Corpses",
"Sediment bacteria & detritus",
"Macrophyte debris",
"Porewater ammonia",
"Porewater nitrate",
"Benthic susp/dep feeders",
"Benthic carn/scav feeders")
rownames(scenariowater)<-rownames(baselinewater)
rownames(scenarioseabed)<-rownames(baselineseabed)
changewater<-baselinewater
changewater[,1]<-NA
changewater[,2]<-NA
colnames(changewater)<-c("LG","PC")
changeseabed<-baselineseabed
changeseabed[,1]<-NA
changeseabed[,2]<-NA
colnames(changeseabed)<-c("LG","PC")
for(zz in 1:nrow(changewater)) {
if( (is.na(scenariowater[zz,1])==FALSE) &
(is.na(baselinewater[zz,1])==FALSE) &
(baselinewater[zz,1]>0) ) {
if(scenariowater[zz,1]>0) changewater[zz,1]<-log10(scenariowater[zz,1]/baselinewater[zz,1])
if(scenariowater[zz,1]==0) changewater[zz,1]<-log10((scenariowater[zz,1]+1E-20)/baselinewater[zz,1])
changewater[zz,2]<-((scenariowater[zz,1]/baselinewater[zz,1])-1)*100
}
}
for(zz in 1:nrow(changeseabed)) {
if( (is.na(scenarioseabed[zz,1])==FALSE) &
(is.na(baselineseabed[zz,1])==FALSE) &
(baselineseabed[zz,1]>0) ) {
if(scenarioseabed[zz,1]>0) changeseabed[zz,1]<-log10(scenarioseabed[zz,1]/baselineseabed[zz,1])
if(scenarioseabed[zz,1]==0) changeseabed[zz,1]<-log10((scenarioseabed[zz,1]+1E-20)/baselineseabed[zz,1])
changeseabed[zz,2]<-((scenarioseabed[zz,1]/baselineseabed[zz,1])-1)*100
}
}
if(log.pc=="LG"){
changewater2p<-as.data.frame(changewater[,1])
changeseabed2p<-as.data.frame(changeseabed[,1])
}
if(log.pc=="PC"){
changewater2p<-as.data.frame(changewater[,2])
changeseabed2p<-as.data.frame(changeseabed[,2])
}
rownames(changewater2p)<-rownames(changewater)
rownames(changeseabed2p)<-rownames(changeseabed)
overlabelwater<-rep("",nrow(changewater2p))
overlabelseabed<-rep("",nrow(changeseabed2p))
for(zz in 1:nrow(changewater2p)){
if(is.na(changewater2p[zz,1])==FALSE){
if(changewater2p[zz,1]<(bpmin)){
overlabelwater[zz]<-as.character( (floor(100*changewater2p[zz,1]))/100 )
changewater2p[zz,1]<-bpmin
}
}
}
for(zz in 1:nrow(changeseabed2p)){
if(is.na(changeseabed2p[zz,1])==FALSE){
if(changeseabed2p[zz,1]<(bpmin)){
overlabelseabed[zz]<-as.character( (floor(100*changeseabed2p[zz,1]))/100 )
changeseabed2p[zz,1]<-bpmin
}
}
}
wcolvec<-rep("green",length(changewater2p[,1]))
wcolvec[which(changewater2p[,1]<0)]<-"red"
scolvec<-rep("green",length(changeseabed2p[,1]))
scolvec[which(changeseabed2p[,1]<0)]<-"red"
par(mfrow=c(2,1))
par(mar=c(1.5,13,2,1))
barplot(changewater2p[,1],horiz=TRUE,col=wcolvec,names.arg=rownames(changewater2p),las=1,xlim=c(bpmin,bpmax),cex.axis=0.9,cex.names=0.75,main=maintitle,cex.main=0.9)
for(zz in 1:length(overlabelwater)){
text(0.7*bpmin,((zz-0.4)*1.2),overlabelwater[zz],cex=0.7)
}
abline(v=0)
par(mar=c(3.5,13,1,1))
barplot(changeseabed2p[,1],horiz=TRUE,col=scolvec,names.arg=rownames(changeseabed2p),las=1,xlim=c(bpmin,bpmax),cex.axis=0.8,cex.names=0.75)
for(zz in 1:length(overlabelseabed)){
text(0.7*bpmin,((zz-0.3)*1.17),overlabelseabed[zz],cex=0.8)
}
if(log.pc=="PC") mtext("Percent change in annual average mass",cex=1,side=1,line=2)
if(log.pc=="LG") mtext("Log10 change in annual average mass",cex=1,side=1,line=2)
abline(v=0)
##################################################################
#Return the two dataframes for change in water column and seabed variables
list(changewater=changewater,
changeseabed=changeseabed)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.