R/compare_two_runs_aam.R

In StrathE2E2: End-to-End Marine Food Web Model

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



}