R/compare_two_runs_catch.R

Defines functions compare_two_runs_catch

#
# 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)


}

Try the StrathE2E2 package in your browser

Any scripts or data that you put into this service are public.

StrathE2E2 documentation built on Jan. 23, 2021, 1:07 a.m.