Nothing
#
# compare_two_runs_catch.R
#
#' Create a tornado plot of the differences in landings and discards in two model runs.
#'
#' Reads final year annual landings and discards results generated by two e2e_run() model runs with different inputs, and plot the difference
#' between the two as tornado diagrams.
#'
#' Inshore and offshore annual landings and discards 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_landingcomposition_by_gear-*.csv and /results/Modelname/Variantname/ZONE_discardcomposition_by_gear-*.csv
#' where ZONE is INSHORE or OFFSHORE 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) for landings (black for discards). Conversely scenario values < baseline values
#' are shown by red bars to the left for landings (grey for discards).
#'
#' 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 (changeland and changedisc). The first if the proportional difference in landings between the scenario run and the baseline
#' the second is the proportional difference in discards. 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: May 2020 |
# | |
# ---------------------------------------------------------------------
compare_two_runs_catch <- 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 landings data
print("Using baseline data held in memory from an existing model run")
final.year.outputs1 <- elt(results1, "final.year.outputs")
if(zone=="I") baselinelanddata <- elt(final.year.outputs1,"inshore_landmat")
if(zone=="O") baselinelanddata <- elt(final.year.outputs1,"offshore_landmat")
if(zone=="W") baselinelanddata <- elt(final.year.outputs1,"inshore_landmat") + elt(final.year.outputs1,"offshore_landmat")
if(zone=="I") baselinediscdata <- elt(final.year.outputs1,"inshore_discmat")
if(zone=="O") baselinediscdata <- elt(final.year.outputs1,"offshore_discmat")
if(zone=="W") baselinediscdata <- elt(final.year.outputs1,"inshore_discmat") + elt(final.year.outputs1,"offshore_discmat")
}
if(from.csv2==FALSE) {
#scenario run landings data
print("Using scenario data held in memory from an existing model run")
final.year.outputs2 <- elt(results2, "final.year.outputs")
if(zone=="I") scenariolanddata <- elt(final.year.outputs2,"inshore_landmat")
if(zone=="O") scenariolanddata <- elt(final.year.outputs2,"offshore_landmat")
if(zone=="W") scenariolanddata <- elt(final.year.outputs2,"inshore_landmat") + elt(final.year.outputs2,"offshore_landmat")
if(zone=="I") scenariodiscdata <- elt(final.year.outputs2,"inshore_discmat")
if(zone=="O") scenariodiscdata <- elt(final.year.outputs2,"offshore_discmat")
if(zone=="W") scenariodiscdata <- elt(final.year.outputs2,"inshore_discmat") + elt(final.year.outputs2,"offshore_discmat")
}
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")
baselandifile <- csvname(resultsdir1, "INSHORE_landingcomposition_by_gear", model.ident1)
baselandofile <- csvname(resultsdir1, "OFFSHORE_landingcomposition_by_gear", model.ident1)
basediscifile <- csvname(resultsdir1, "INSHORE_discardcomposition_by_gear", model.ident1)
basediscofile <- csvname(resultsdir1, "OFFSHORE_discardcomposition_by_gear", model.ident1)
check.exists(baselandifile)
check.exists(baselandofile)
check.exists(basediscifile)
check.exists(basediscofile)
baselinelandidata<-readcsv(baselandifile,row.names=1)
baselinelandodata<-readcsv(baselandofile,row.names=1)
baselinediscidata<-readcsv(basediscifile,row.names=1)
baselinediscodata<-readcsv(basediscofile,row.names=1)
if(zone=="I") {
baselinelanddata <- baselinelandidata
baselinediscdata <- baselinediscidata
}
if(zone=="O") {
baselinelanddata <- baselinelandodata
baselinediscdata <- baselinediscodata
}
if(zone=="W") {
baselinelanddata <- baselinelandidata + baselinelandodata
baselinediscdata <- baselinediscidata + baselinediscodata
}
}
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")
scenlandifile <- csvname(resultsdir2, "INSHORE_landingcomposition_by_gear", model.ident2)
scenlandofile <- csvname(resultsdir2, "OFFSHORE_landingcomposition_by_gear", model.ident2)
scendiscifile <- csvname(resultsdir2, "INSHORE_discardcomposition_by_gear", model.ident2)
scendiscofile <- csvname(resultsdir2, "OFFSHORE_discardcomposition_by_gear", model.ident2)
check.exists(scenlandifile)
check.exists(scenlandofile)
check.exists(scendiscifile)
check.exists(scendiscofile)
scenariolandidata<-readcsv(scenlandifile,row.names=1)
scenariolandodata<-readcsv(scenlandofile,row.names=1)
scenariodiscidata<-readcsv(scendiscifile,row.names=1)
scenariodiscodata<-readcsv(scendiscofile,row.names=1)
if(zone=="I") {
scenariolanddata <- scenariolandidata
scenariodiscdata <- scenariodiscidata
}
if(zone=="O") {
scenariolanddata <- scenariolandodata
scenariodiscdata <- scenariodiscodata
}
if(zone=="W") {
scenariolanddata <- scenariolandidata+scenariolandodata
scenariodiscdata <- scenariodiscidata+scenariodiscodata
}
}
baselineland <- data.frame(rep(0,11))
scenarioland <- data.frame(rep(0,11))
baselinedisc <- data.frame(rep(0,11))
scenariodisc <- data.frame(rep(0,11))
baselineland[1,1]<-sum(baselinelanddata[11,],na.rm=TRUE) #Macrophytes
baselineland[2,1]<-sum(baselinelanddata[5,],na.rm=TRUE) #susp/dep benthos
baselineland[3,1]<-sum(baselinelanddata[6,],na.rm=TRUE) #carn/scav benthos
baselineland[4,1]<-sum(baselinelanddata[7,],na.rm=TRUE) #pel inver
baselineland[5,1]<-sum(baselinelanddata[1,],na.rm=TRUE) #planktiv fish
baselineland[6,1]<-sum(baselinelanddata[4,],na.rm=TRUE) #Mig fish
baselineland[7,1]<-sum(baselinelanddata[3,],na.rm=TRUE) #NQ demersal fish
baselineland[8,1]<-sum(baselinelanddata[2,],na.rm=TRUE) #QL demersal fish
baselineland[9,1]<-sum(baselinelanddata[8,],na.rm=TRUE) #Birds
baselineland[10,1]<-sum(baselinelanddata[9,],na.rm=TRUE) #Pinnipeds
baselineland[11,1]<-sum(baselinelanddata[10,],na.rm=TRUE) #Cetaceans
scenarioland[1,1]<-sum(scenariolanddata[11,],na.rm=TRUE) #Macrophytes
scenarioland[2,1]<-sum(scenariolanddata[5,],na.rm=TRUE) #susp/dep benthos
scenarioland[3,1]<-sum(scenariolanddata[6,],na.rm=TRUE) #carn/scav benthos
scenarioland[4,1]<-sum(scenariolanddata[7,],na.rm=TRUE) #pel inver
scenarioland[5,1]<-sum(scenariolanddata[1,],na.rm=TRUE) #planktiv fish
scenarioland[6,1]<-sum(scenariolanddata[4,],na.rm=TRUE) #Mig fish
scenarioland[7,1]<-sum(scenariolanddata[3,],na.rm=TRUE) #NQ demersal fish
scenarioland[8,1]<-sum(scenariolanddata[2,],na.rm=TRUE) #QL demersal fish
scenarioland[9,1]<-sum(scenariolanddata[8,],na.rm=TRUE) #Birds
scenarioland[10,1]<-sum(scenariolanddata[9,],na.rm=TRUE) #Pinnipeds
scenarioland[11,1]<-sum(scenariolanddata[10,],na.rm=TRUE) #Cetaceans
baselinedisc[1,1]<-sum(baselinediscdata[11,],na.rm=TRUE) #Macrophytes
baselinedisc[2,1]<-sum(baselinediscdata[5,],na.rm=TRUE) #susp/dep benthos
baselinedisc[3,1]<-sum(baselinediscdata[6,],na.rm=TRUE) #carn/scav benthos
baselinedisc[4,1]<-sum(baselinediscdata[7,],na.rm=TRUE) #pel inver
baselinedisc[5,1]<-sum(baselinediscdata[1,],na.rm=TRUE) #planktiv fish
baselinedisc[6,1]<-sum(baselinediscdata[4,],na.rm=TRUE) #Mig fish
baselinedisc[7,1]<-sum(baselinediscdata[3,],na.rm=TRUE) #NQ demersal fish
baselinedisc[8,1]<-sum(baselinediscdata[2,],na.rm=TRUE) #QL demersal fish
baselinedisc[9,1]<-sum(baselinediscdata[8,],na.rm=TRUE) #Birds
baselinedisc[10,1]<-sum(baselinediscdata[9,],na.rm=TRUE) #Pinnipeds
baselinedisc[11,1]<-sum(baselinediscdata[10,],na.rm=TRUE) #Cetaceans
scenariodisc[1,1]<-sum(scenariodiscdata[11,],na.rm=TRUE) #Macrophytes
scenariodisc[2,1]<-sum(scenariodiscdata[5,],na.rm=TRUE) #susp/dep benthos
scenariodisc[3,1]<-sum(scenariodiscdata[6,],na.rm=TRUE) #carn/scav benthos
scenariodisc[4,1]<-sum(scenariodiscdata[7,],na.rm=TRUE) #pel inver
scenariodisc[5,1]<-sum(scenariodiscdata[1,],na.rm=TRUE) #planktiv fish
scenariodisc[6,1]<-sum(scenariodiscdata[4,],na.rm=TRUE) #Mig fish
scenariodisc[7,1]<-sum(scenariodiscdata[3,],na.rm=TRUE) #NQ demersal fish
scenariodisc[8,1]<-sum(scenariodiscdata[2,],na.rm=TRUE) #QL demersal fish
scenariodisc[9,1]<-sum(scenariodiscdata[8,],na.rm=TRUE) #Birds
scenariodisc[10,1]<-sum(scenariodiscdata[9,],na.rm=TRUE) #Pinnipeds
scenariodisc[11,1]<-sum(scenariodiscdata[10,],na.rm=TRUE) #Cetaceans
rownames(baselineland) <-c("Macrophytes",
"Susp/dep feeding benthos",
"Carn/scav feeding benthos",
"Pelagic invertebrates",
"Planktivorous fish",
"Migratory fish",
"Non-quota demersal fish",
"Quota-limited demersal fish",
"Birds",
"Pinnipeds",
"LANDINGS Cetaceans")
rownames(baselinedisc) <-c("Macrophytes",
"Susp/dep feeding benthos",
"Carn/scav feeding benthos",
"Pelagic invertebrates",
"Planktivorous fish",
"Migratory fish",
"Non-quota demersal fish",
"Quota-limited demersal fish",
"Birds",
"Pinnipeds",
"DISCARDS Cetaceans")
rownames(baselinedisc)<-rownames(baselinedisc)
rownames(scenarioland)<-rownames(baselineland)
rownames(scenariodisc)<-rownames(baselinedisc)
changeland<-baselineland
changeland[,1]<-NA
changeland[,2]<-NA
colnames(changeland)<-c("LG","PC")
changedisc<-baselinedisc
changedisc[,1]<-NA
changedisc[,2]<-NA
colnames(changedisc)<-c("LG","PC")
for(zz in 1:nrow(changeland)) {
if( (is.na(scenarioland[zz,1])==FALSE) &
(is.na(baselineland[zz,1])==FALSE) &
(baselineland[zz,1]>0) ) {
if(scenarioland[zz,1]>0) changeland[zz,1]<-log10(scenarioland[zz,1]/baselineland[zz,1])
if(scenarioland[zz,1]==0) changeland[zz,1]<-log10((scenarioland[zz,1]+1E-20)/baselineland[zz,1])
changeland[zz,2]<-((scenarioland[zz,1]/baselineland[zz,1])-1)*100
}
}
for(zz in 1:nrow(changedisc)) {
if( (is.na(scenariodisc[zz,1])==FALSE) &
(is.na(baselinedisc[zz,1])==FALSE) &
(baselinedisc[zz,1]>0) ) {
if(scenariodisc[zz,1]>0) changedisc[zz,1]<-log10(scenariodisc[zz,1]/baselinedisc[zz,1])
if(scenariodisc[zz,1]==0) changedisc[zz,1]<-log10((scenariodisc[zz,1]+1E-20)/baselinedisc[zz,1])
changedisc[zz,2]<-((scenariodisc[zz,1]/baselinedisc[zz,1])-1)*100
}
}
if(log.pc=="LG"){
changeland2p<-as.data.frame(changeland[,1])
changedisc2p<-as.data.frame(changedisc[,1])
}
if(log.pc=="PC"){
changeland2p<-as.data.frame(changeland[,2])
changedisc2p<-as.data.frame(changedisc[,2])
}
rownames(changeland2p)<-rownames(changeland)
rownames(changedisc2p)<-rownames(changedisc)
overlabelland<-rep("",nrow(changeland2p))
overlabeldisc<-rep("",nrow(changedisc2p))
for(zz in 1:nrow(changeland2p)){
if(is.na(changeland2p[zz,1])==FALSE){
if(changeland2p[zz,1]<(bpmin)){
overlabelland[zz]<-as.character( (floor(100*changeland2p[zz,1]))/100 )
changeland2p[zz,1]<-bpmin
}
}
}
for(zz in 1:nrow(changedisc2p)){
if(is.na(changedisc2p[zz,1])==FALSE){
if(changedisc2p[zz,1]<(bpmin)){
overlabeldisc[zz]<-as.character( (floor(100*changedisc2p[zz,1]))/100 )
changedisc2p[zz,1]<-bpmin
}
}
}
wcolvec<-rep("green",length(changeland2p[,1]))
wcolvec[which(changeland2p[,1]<0)]<-"red"
scolvec<-rep("black",length(changedisc2p[,1]))
scolvec[which(changedisc2p[,1]<0)]<-"grey"
par(mfrow=c(2,1))
par(mar=c(1.5,13,2,1))
barplot(changeland2p[,1],horiz=TRUE,col=wcolvec,names.arg=rownames(changeland2p),las=1,xlim=c(bpmin,bpmax),cex.axis=0.9,cex.names=0.9,main=maintitle,cex.main=0.9)
for(zz in 1:length(overlabelland)){
text(0.7*bpmin,((zz-0.4)*1.2),overlabelland[zz],cex=0.8)
}
abline(v=0)
par(mar=c(3.5,13,1,1))
barplot(changedisc2p[,1],horiz=TRUE,col=scolvec,names.arg=rownames(changedisc2p),las=1,xlim=c(bpmin,bpmax),cex.axis=0.9,cex.names=0.9)
for(zz in 1:length(overlabeldisc)){
text(0.7*bpmin,((zz-0.4)*1.2),overlabeldisc[zz],cex=0.8)
}
if(log.pc=="PC") mtext("Percent change in landings & discards",cex=1,side=1,line=2)
if(log.pc=="LG") mtext("Log10 change in landings & discards",cex=1,side=1,line=2)
abline(v=0)
##################################################################
#Return the two dataframes for change in landings and discards
#Replace the rownames to remove the plaot labels LANDINGS and DISCARDS
rownames(changeland) <-c("Macrophytes",
"Susp/dep feeding benthos",
"Carn/scav feeding benthos",
"Pelagic invertebrates",
"Planktivorous fish",
"Migratory fish",
"Non-quota demersal fish",
"Quota-limited demersal fish",
"Birds",
"Pinnipeds",
"Cetaceans")
rownames(changedisc) <-c("Macrophytes",
"Susp/dep feeding benthos",
"Carn/scav feeding benthos",
"Pelagic invertebrates",
"Planktivorous fish",
"Migratory fish",
"Non-quota demersal fish",
"Quota-limited demersal fish",
"Birds",
"Pinnipeds",
"Cetaceans")
list(changeland=changeland,
changedisc=changedisc)
}
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.