R/plot_inshore_vs_offshore_anavmass_with_ci.R

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

#
# plot_inshore_vs_offshore_anavmass_with_ci.R
#
#' Plots of annual average inshore and offshore mass densities generated a Monte Carlo simulation
#'
#' Create a plot showing the distribution of credible values of annual average inshore and offshore mass densities over the final year of a simulation generated by the Monte Carlo e2e_run_mc() function.
#'
#' Daily interval post-processed data from the e2e_run_mc() function are stored in the files
#' /results/Modelname/Variantname/CredInt/CredInt_processed_AAMresults_offshore-*.csv and CredInt_processed_AAMresults_inshore-*.csv,
#' where * 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 an example data set for one of the two North Sea model variants supplied with the package.
#'
#' The format of the output is a single panel with rows of horizintally orientated box-and-whisker plots. Each row represents a different
#' state variable in the model (bacteria and detritus, phytoplankton etc etc). Each row consists of two box-and-whiskers, one (red) fof the offshore
#' zone and the other (blue) for the inshore zone. The data represented are distributions of credible values of the mass densities (mMn/m2) of each variable, averaged
#' over the final year of a model run. The distribution of credible values is derived from
#' a Monte Carlo simulation scheme performed by the e2e_run_mc() function.
#'
#' In each plot the central tick mark of the box-and-whisker represents the median mass density, the box represents the 50% credible interval (spanning quartiles of the cumulative likelihood
#' of simulated values) and the whiskers represent the 99% credible intervals.
#'
#' For details of how the distribution of credible output values from StrathE2E are calculated see help(e2e_run_mc).
#'
#' @param model R-list object defining the model configuration compiled by the e2e_read() function
#' @param use.example Logical. If TRUE use pre-computed example data from the internal North Sea model rather than user-generated data, (default=FALSE).
#'
#' @return Graphical display in a new graphics window
#'
#' @importFrom graphics text
#'
#' @noRd
#
# ------------------------------------------------------------------------------

plot_inshore_vs_offshore_anavmass_with_ci <- function(model, use.example=FALSE) {

start_par = par()$mfrow
on.exit(par(mfrow = start_par))

	resultsdir	<- elt(model, "setup", "resultsdir")
	model.ident	<- elt(model, "setup", "model.ident")
	model.name 	<- elt(model, "setup", "model.name")
	model.variant 	<- elt(model, "setup", "model.variant")

corefilename_o<-"CredInt_processed_AAMresults_offshore"
corefilename_i<-"CredInt_processed_AAMresults_inshore"

if(use.example==TRUE){
	O_results <- get.example.results(model.name, model.variant, corefilename_o, CREDINT_DIR)
	I_results <- get.example.results(model.name, model.variant, corefilename_i, CREDINT_DIR)
}

if(use.example==FALSE){
	credpath	<- makepath(resultsdir, CREDINT_DIR)
	credfile_o	<- csvname(credpath, corefilename_o, model.ident)
	credfile_i	<- csvname(credpath, corefilename_i, model.ident)
	message("Reading credible interval processed data from '", credfile_o, "'")
	message("Reading credible interval processed data from '", credfile_i, "'")

	O_results<-readcsv(credfile_o, row.names=1)
	I_results<-readcsv(credfile_i, row.names=1)
}

	Nind<-13

	#Plot the inshore and offshore results

	#plotlabels<-c("bact.&.detritus","phyt","s/d.benth","omniv.zoo","c/s.benth","carn.zoo","plankt.fish","mig.fish","dem.fish","bird.&.mam")
	plotlabels<-c("Bact.&.detritus","Macrophytes","Phytoplankton","S/d.benthos","Omniv.zoo","C/s.benthos","Carn.zoo","Plankt.fish","Mig.fish","Dem.fish","Birds","Pinnipeds","Cetaceans")

	modeldata2plot<-log10(O_results[2:6,1])
	for(jj in 2:Nind) { modeldata2plot<-c(modeldata2plot,log10(O_results[2:6,jj]))}
	array2plot<- array(dim=c(5,Nind),modeldata2plot)
	colnames(array2plot)<-plotlabels
	bxpdata_O<-list(stats=array2plot,n=rep(100,Nind),conf=NULL,out=numeric(length=0),names=plotlabels)

	modeldata2plot<-log10(I_results[2:6,1])
	for(jj in 2:Nind) { modeldata2plot<-c(modeldata2plot,log10(I_results[2:6,jj]))}
	array2plot<- array(dim=c(5,Nind),modeldata2plot)
	colnames(array2plot)<-plotlabels
	bxpdata_I<-list(stats=array2plot,n=rep(100,Nind),conf=NULL,out=numeric(length=0),names=plotlabels)

        combdata<-rbind(bxpdata_O$stats,bxpdata_I$stats)
        combdata[1:5,2]<-NA  # reset offshore macrophytes to NA instead of -Inf
	ulm<-max(combdata[c(5,10),],na.rm=TRUE)
	llm<-min(combdata[c(2,7),],na.rm=TRUE)

	ymin<-floor(llm)
	ymax<-ceiling(ulm)

	par(mar=c(5,8,1,1))
	bxp(bxpdata_O,boxwex=0.35,at=seq(1,Nind),ylim=c(ymin,ymax),cex.axis=1.2,show.names=TRUE,horizontal=TRUE,las=1,boxcol="red",whiskcol="red",whisklty="solid",medcol="red",staplecol="red")
	#axis(side=1,las=1,cex.axis=1.2)
	#axis(side=1,las=1,cex.axis=1.2,at=c(-1,0,-1,2))
#	mtext("Log(annual average biomass density)",cex=1.2,side=1,line=3)
	mtext(bquote("Annual average biomass log"[10] ~ "mMN.m"^-2),cex=1.2,side=1,line=3)
	bxp(bxpdata_I,add=TRUE,boxwex=0.35,at=seq(1,Nind)+0.35,yaxt="n",xaxt="n",horizontal=TRUE,boxcol="blue",whiskcol="blue",whisklty="solid",medcol="blue",staplecol="blue")

	legend("bottomleft",box.lty=0,bg="transparent",c("Inshore","Offshore"),col=c("blue","red"),lty=c(1,1),lwd=c(2,2),pt.cex=c(1,1),cex=c(0.9,0.9))

	abline(h=1.7,lty="dashed")
	abline(h=2.7,lty="dashed")
	abline(h=3.7,lty="dashed")
	abline(h=4.7,lty="dashed")
	abline(h=5.7,lty="dashed")
	abline(h=6.7,lty="dashed")
	abline(h=7.7,lty="dashed")
	abline(h=8.7,lty="dashed")
	abline(h=9.7,lty="dashed")
	abline(h=10.7,lty="dashed")
	abline(h=11.7,lty="dashed")
	abline(h=12.7,lty="dashed")

}