StrathE2E2
End-to-End Marine Food Web Model

Package index
Search the StrathE2E2 package
Vignettes
Functions
187
Source code
138
Man pages
32
Home
/
CRAN
/
StrathE2E2
/
R/boxplot_annual_compare_runs.R

R/boxplot_annual_compare_runs.R
In StrathE2E2: End-to-End Marine Food Web Model

Defines functions boxplot_annual_compare_runs 

#
# boxplot_annual_compare_runs.R
#
#' Box-and-whisker plots of baseline and scenario simulated annual target data on the state of the ecosystem with credible intervals derived from Monte Carlo StrathE2E runs.
#'
#' Creates a multi-panel plot comparing a range of simulated annual data on the state of the ecosystem from a baseline and a scenario mode run, with the distribution of credible values
#' generated by the e2e_run_mc() Monte Carlo function.
#'
#' For details of how the distribution of credible output values from StrathE2E are calculated see ?e2e_run_mc
#'
#' The function plots a multi-panel page of box-and-whisker plots showing the medians and variability ranges (quartiles as box-and-whisker) of a range of indices of ecosystem state from the final year of a baseline run
#' (shown in black), alongside comparable box-and-whisker plots (shown in red) of the same measures derived from the final year of a scenario run. The credible intervals for each case need ot be generted by the Monte Carlo methodology (e2e_run_mc() function).
#'
#' Individual panels of the plot represent groups of similar or related ecosystem properties - annual productions, annual fishery landings, annual P:B ratios, annual food consumtion and diet compositions, nutrient concentrations, and inshore : offshore abundanace ratios.
#'
#' Optionally the function can read an example data set for one of the two North Sea model variants supplied with the package.
#'
#' @param model1 R-list object defining the baseline model configuration compiled by the e2e_read() function.
#' @param ci.data1 Logical. If TRUE plot credible intervals around model results based on Monte Carlo simulation withteh e2e_run_mc() function, (default=FALSE).
#' @param use.saved1 Logical. If TRUE use data from a prior user-defined run held as csv files data in the current results folder, (default=FALSE).
#' @param use.example1 Logical. If TRUE use pre-computed example data from the internal North Sea model as the baseline rather than user-generated data, (default=FALSE).
#' @param results1 R-list object of baseline model output generated by the e2e_run(), (default=NULL).
#' @param model2 R-list object defining the baseline model configuration compiled by the e2e_read() function.
#' @param ci.data2 Logical. If TRUE plot credible intervals around model results based on Monte Carlo simulation withteh e2e_run_mc() function, (default=FALSE).
#' @param use.saved2 Logical. If TRUE use data from a prior user-defined run held as csv files data in the current results folder, (default=FALSE).
#' @param use.example2 Logical. If TRUE use pre-computed example data from the internal North Sea model as the scenario rather than user-generated data, (default=FALSE).
#' @param results2 R-list object of scenario model output generated by the e2e_run(), (default=NULL).
#'
#' @return Graphical display in a new graphics window.
#'
#' @importFrom graphics grconvertX grconvertY
#'
#' @noRd
#
# ---------------------------------------------------------------------
# |                                                                   |
# | Authors: Mike Heath, Ian Thurlbeck                                |
# | Department of Mathematics and Statistics                          |
# | University of Strathclyde, Glasgow                                |
# |                                                                   |
# | Date of this version: May 2020                                    |
# |                                                                   |
# ---------------------------------------------------------------------

boxplot_annual_compare_runs <- function(model1, ci.data1=FALSE, use.saved1=FALSE, use.example1=FALSE, results1=NULL,
                                        model2, ci.data2=FALSE, use.saved2=FALSE, use.example2=FALSE, results2=NULL) {

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

	resultsdir1		<- elt(model1, "setup", "resultsdir")
	model.ident1		<- elt(model1, "setup", "model.ident")
	model.path1		<- elt(model1, "setup", "model.path")
	model.name1		<- elt(model1, "setup", "model.name")
	model.variant1 		<- elt(model1, "setup", "model.variant")

	resultsdir2		<- elt(model2, "setup", "resultsdir")
	model.ident2		<- elt(model2, "setup", "model.ident")
	model.path2		<- elt(model2, "setup", "model.path")
	model.name2		<- elt(model2, "setup", "model.name")
	model.variant2 		<- elt(model2, "setup", "model.variant")

	annualtargetdata	<- read_annual_target_data(model.path1)

	corefilename<-"CredInt_processed_targetresults"

if(ci.data1==TRUE){

if(use.example1==TRUE){
	credintervaldata1 <- get.example.results(model.name1, model.variant1, corefilename, CREDINT_DIR)
}

if(use.example1==FALSE){
	credpath1	<- makepath(resultsdir1, CREDINT_DIR)
	credfile1	<- csvname(credpath1, corefilename, model.ident1)
	if (! file.exists(credfile1)) {
		message("Error: cannot find credible interval output file: ", credfile1)
		stop("Please run the Monte Carlo function!\n")
	}
	message("Reading credible interval processed data from '", credfile1, "'")
	credintervaldata1	<- readcsv(credfile1, row.names=1)	# first column is row names
}


}


if(ci.data2==TRUE){

if(use.example2==TRUE){
	credintervaldata2 <- get.example.results(model.name2, model.variant2, corefilename, CREDINT_DIR)
}

if(use.example2==FALSE){
	credpath2	<- makepath(resultsdir2, CREDINT_DIR)
	credfile2	<- csvname(credpath2, corefilename, model.ident2)
	if (! file.exists(credfile2)) {
		message("Error: cannot find credible interval output file: ", credfile2)
		stop("Please run the Monte Carlo function!\n")
	}
	message("Reading credible interval processed data from '", credfile2, "'")
	credintervaldata2	<- readcsv(credfile2, row.names=1)	# first column is row names
}

}


if(ci.data1==FALSE){

if(use.saved1==TRUE){
	datafile1	<- csvname(resultsdir1, "model_target_annualresults_plus_chi", model.ident1)
	print(paste("Using baseline data held in a file ",datafile1," from a past model run"))
	check.exists(datafile1)
	opt_results1<-readcsv(datafile1)
}

if(use.saved1==FALSE){
	opt_results1	<- elt(results1, "final.year.outputs", "opt_results")
}

}


if(ci.data2==FALSE){

if(use.saved2==TRUE){
	datafile2	<- csvname(resultsdir2, "model_target_annualresults_plus_chi", model.ident2)
	print(paste("Using scenario data held in a file ",datafile2," from a past model run"))
	check.exists(datafile2)
	opt_results2<-readcsv(datafile2)
}

if(use.saved2==FALSE){
	opt_results2	<- elt(results2, "final.year.outputs", "opt_results")
}

}

# --------------------------------------------------------------------------


if(ci.data1==FALSE){
    #convert optresults into credintervaldata format
	credintervaldata1 <- data.frame(rep(rep(NA,6)))
	for(jj in 2:nrow(opt_results1)){
		credintervaldata1[,jj]<-rep(rep(NA,6))
	}
	row.names(credintervaldata1)<-c("maxlik","lowlimitit","lowquart","median","uppquart","upplimit")
	for(jj in 1:6){
		credintervaldata1[jj,]<-opt_results1[,3]
	}
}

if(ci.data2==FALSE){
    #convert optresults into credinterval data format
	credintervaldata2 <- data.frame(rep(rep(NA,6)))
	for(jj in 2:nrow(opt_results2)){
		credintervaldata2[,jj]<-rep(rep(NA,6))
	}
	row.names(credintervaldata2)<-c("maxlik","lowlimitit","lowquart","median","uppquart","upplimit")
	for(jj in 1:6){
		credintervaldata2[jj,]<-opt_results2[,3]
	}
}


# -----------------------------------------------------------------------------------------

	#Here we are reading the database of annual target data file associated with model1 to gain access to variable names
	#so as to group them together into families for plotting and assign text strings to appear on plot axes

	annualtargetnames<-rep("xx",nrow(annualtargetdata))

	annualtargetnames[which(annualtargetdata[,4]=="Obs_TAPP")] <- "Total phyt."
	annualtargetnames[which(annualtargetdata[,4]=="Obs_NP")] <- "New primary phyt."
	annualtargetnames[which(annualtargetdata[,4]=="Obs_KelpP")] <- "Macrophyte carbon"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_OmnizooP")] <- "Omniv.zooplankton"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_CarnzooP")] <- "Carniv.zooplankton"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_PFishP")] <- "Planktiv.fish"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_DFishP")] <- "Demersal fish"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_BensuspP")] <- "Susp/dep.benthos"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_BencarnP")] <- "Carn/scav.benthos"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_birdP")] <- "Seabirds"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_sealP")] <- "Pinnipeds"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_cetaP")] <- "Cetaceans"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_maxbenthslar")] <- "Susp/dep.benthos larv"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_maxbenthclar")] <- "Carn/scav.benthos larv"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Conpfishfish")] <- "Plank.fish by fish"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Condfishfish")] <- "Dem.fish by fish"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Conzoofish")] <- "Zooplankton by fish"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Conzoocarnz")] <- "Omniv.zoo by carniv.zoo."
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Conbenfish")] <- "Benthos by fish"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Contotal_bird")] <- "Total by birds"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Proppfishbird")] <- "Plank.fish in bird diet"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Propdfishbird")] <- "Dem.fish in bird diet"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Propmfishbird")] <- "Mig.fish in bird diet"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Propdiscbird")] <- "Disc. in bird diet"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Contotal_seal")] <- "Total by pinnipeds"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Proppfishseal")] <- "Plank.fish in pinn. diet"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Propdfishseal")] <- "Dem.fish in pinn. diet"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Propmfishseal")] <- "Mig.fish in pinn. diet"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Contotal_ceta")] <- "Total by cetaceans"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Proppfishceta")] <- "Plank.fish in cet. diet"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Propdfishceta")] <- "Dem.fish in cet. diet"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Propmfishceta")] <- "Mig.fish in cet. diet"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Pland_livewt")] <- "Plank.fish landings"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Dland_livewt")] <- "Dem.fish landings"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Mland_livewt")] <- "Mig.fish landings"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Bsland_livewt")] <- "Susp/dep.benthos landings" 
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Bcland_livewt")] <- "Carn/scav.benthos landings"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Zcland_livewt")] <- "Pel.invert. landings"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Kland_livewt")] <- "Macrop. harvest"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_kelp_pb")] <- "Macrop. P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_benslar_pb")] <- "Susp/dep.benthos larv. P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_benclar_pb")] <- "Carn/scav.benthos larv. P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_bens_pb")] <- "Susp/dep.benthos P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_benc_pb")] <- "Carn/scav.benthos P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_omni_pb")] <- "Omniv.zooplankton P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_carn_pb")] <- "Carniv.zooplankton P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_fishplar_pb")] <- "Plank.fish larvae P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_fishdlar_pb")] <- "Dem.fish larvae P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_fishp_pb")] <- "Plank.fish P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_fishd_pb")] <- "Dem.fish P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_fishm_pb")] <- "Mig.fish P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_bird_pb")] <- "Bird P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_seal_pb")] <- "Pinniped P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_ceta_pb")] <- "Cetacean P/B"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_exud_C_kelp")] <- "Prop. macrop. prod. exuded"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_kelp_NC")] <- "Macrop. N/C ratio"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Denitrif")] <- "Denitrification"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Dfdiscardp")] <- "Dem.fish discard/catch"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_s_x_ammonia")] <- "Sand porewater ammonia"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_d_x_ammonia")] <- "Mud porewater ammonia"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_s_x_nitrate")] <- "Sand porewater nitrate"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_d_x_nitrate")] <- "Mud porewater nitrate"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_s_x_TON")] <- "Sand %TON"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_d_x_TON")] <- "Mud %TON"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_NDJF_s_nitrate")] <- "Winter surf.nitrate"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_MJJA_s_nitrate")] <- "Summer surf.nitrate"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_NDJF_d_nitrate")] <- "Winter deep nitrate"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_MJJA_d_nitrate")] <- "Summer deep nitrate"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_NDJF_s_ammonia")] <- "Winter surf.ammonia"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_MJJA_s_ammonia")] <- "Summer surf.ammonia"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_NDJF_d_ammonia")] <- "Winter deep ammonia"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_MJJA_d_ammonia")] <- "Summer deep ammonia"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_carn_io_ratio")] <- "Carniv.zooplanton"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_omni_io_ratio")] <- "Omniv.zooplankton"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_phyt_io_ratio")] <- "Surf.phytoplankton"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_nit_io_ratio")] <- "Surf.nitrate"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_amm_io_ratio")] <- "Surf.ammonia"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_pfish_io_ratio")] <- "Plank.fish"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_dfish_io_ratio")] <- "Dem.fish"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_birddisc")] <- "Bird by-catch"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_sealdisc")] <- "Pinniped by-catch"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_cetadisc")] <- "Cetacean by-catch"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_kelp_beachcast")] <- "Macrop. beach-cast"
	annualtargetnames[which(annualtargetdata[,4]=="Obs_Ctland_livewt")] <- "Cetacean landings"

	#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


	# Assign the new sanitized names to the columns in the credint dataframes.
	names(credintervaldata1)<-annualtargetnames
	names(credintervaldata2)<-annualtargetnames

	#Here we assign the data to families

	#IDs for Annual production rates
	set2plot1.1<-c(
		which(annualtargetdata[,4]=="Obs_KelpP"),
		which(annualtargetdata[,4]=="Obs_TAPP"),
		which(annualtargetdata[,4]=="Obs_NP"),
		which(annualtargetdata[,4]=="Obs_Denitrif"),
		which(annualtargetdata[,4]=="Obs_OmnizooP"),
		which(annualtargetdata[,4]=="Obs_CarnzooP"),
		which(annualtargetdata[,4]=="Obs_BensusP"),
		which(annualtargetdata[,4]=="Obs_BencarnP"),
		which(annualtargetdata[,4]=="Obs_PFishP"),
		which(annualtargetdata[,4]=="Obs_DFishP"),
		which(annualtargetdata[,4]=="Obs_birdP"),
		which(annualtargetdata[,4]=="Obs_sealP"),
		which(annualtargetdata[,4]=="Obs_cetaP")
	)

	#IDs for Annual fishery landings and by-catch (live weights)
	set2plot1.2<-c(
		which(annualtargetdata[,4]=="Obs_Pland_livewt"),
		which(annualtargetdata[,4]=="Obs_Dland_livewt"),
		which(annualtargetdata[,4]=="Obs_Mland_livewt"),
		which(annualtargetdata[,4]=="Obs_Bsland_livewt"),
		which(annualtargetdata[,4]=="Obs_Bcland_livewt"),
		which(annualtargetdata[,4]=="Obs_Zcland_livewt"),
		which(annualtargetdata[,4]=="Obs_Ctland_livewt"),
		which(annualtargetdata[,4]=="Obs_Kland_livewt"),
		which(annualtargetdata[,4]=="Obs_cetadisc"),
		which(annualtargetdata[,4]=="Obs_sealdisc"),
		which(annualtargetdata[,4]=="Obs_birddisc"),
		which(annualtargetdata[,4]=="Obs_Dfdiscardp")
	)

	#IDs for annual consumption rates
	set2plot2<-c(
		which(annualtargetdata[,4]=="Obs_Conzoocarnz"),
		which(annualtargetdata[,4]=="Obs_Conzoofish"),
		which(annualtargetdata[,4]=="Obs_Conbenfish"),
		which(annualtargetdata[,4]=="Obs_Conpfishfish"),
		which(annualtargetdata[,4]=="Obs_Condfishfish"),
		which(annualtargetdata[,4]=="Obs_Contotal_bird"),
		which(annualtargetdata[,4]=="Obs_Contotal_seal"),
		which(annualtargetdata[,4]=="Obs_Contotal_ceta"),

		which(annualtargetdata[,4]=="Obs_Proppfishbird"),
		which(annualtargetdata[,4]=="Obs_Propdfishbird"),
		which(annualtargetdata[,4]=="Obs_Propmfishbird"),
		which(annualtargetdata[,4]=="Obs_Propdiscbird"),
		which(annualtargetdata[,4]=="Obs_Proppfishseal"),
		which(annualtargetdata[,4]=="Obs_Propdfishseal"),
		which(annualtargetdata[,4]=="Obs_Propmfishseal"),
		which(annualtargetdata[,4]=="Obs_Proppfishceta"),
		which(annualtargetdata[,4]=="Obs_Propdfishceta"),
		which(annualtargetdata[,4]=="Obs_Propmfishceta")
	)

	#Annual PB and other ratios
	set2plot3<-c(
		which(annualtargetdata[,4]=="Obs_kelp_pb"),
		which(annualtargetdata[,4]=="Obs_kelp_NC"),
		which(annualtargetdata[,4]=="Obs_exud_C_kelp"),
		which(annualtargetdata[,4]=="Obs_kelp_beachcast"),

		which(annualtargetdata[,4]=="Obs_omni_pb"),
		which(annualtargetdata[,4]=="Obs_benslar_pb"),
		which(annualtargetdata[,4]=="Obs_benclar_pb"),
		which(annualtargetdata[,4]=="Obs_fishplar_pb"),
		which(annualtargetdata[,4]=="Obs_fishdlar_pb"),
		which(annualtargetdata[,4]=="Obs_carn_pb"),
		which(annualtargetdata[,4]=="Obs_bens_pb"),
		which(annualtargetdata[,4]=="Obs_benc_pb"),
		which(annualtargetdata[,4]=="Obs_fishp_pb"),
		which(annualtargetdata[,4]=="Obs_fishd_pb"),
		which(annualtargetdata[,4]=="Obs_fishm_pb"),
		which(annualtargetdata[,4]=="Obs_bird_pb"),
		which(annualtargetdata[,4]=="Obs_seal_pb"),
		which(annualtargetdata[,4]=="Obs_ceta_pb")
	)

	#Average nutrient concentaryions in water and sediments
	set2plot4<-c(
		which(annualtargetdata[,4]=="Obs_s_x_ammonia"),
		which(annualtargetdata[,4]=="Obs_d_x_ammonia"),
		which(annualtargetdata[,4]=="Obs_s_x_nitrate"),
		which(annualtargetdata[,4]=="Obs_d_x_nitrate"),
		which(annualtargetdata[,4]=="Obs_s_x_TON"),
		which(annualtargetdata[,4]=="Obs_d_x_TON"),
		which(annualtargetdata[,4]=="Obs_NDJF_s_nitrate"),
		which(annualtargetdata[,4]=="Obs_MJJA_s_nitrate"),
		which(annualtargetdata[,4]=="Obs_NDJF_d_nitrate"),
		which(annualtargetdata[,4]=="Obs_MJJA_d_nitrate"),
		which(annualtargetdata[,4]=="Obs_NDJF_s_ammonia"),
		which(annualtargetdata[,4]=="Obs_MJJA_s_ammonia"),
		which(annualtargetdata[,4]=="Obs_NDJF_d_ammonia"),
		which(annualtargetdata[,4]=="Obs_MJJA_d_ammonia")
	)

	#Inshore:offshore ratios
	set2plot5<-c(
		which(annualtargetdata[,4]=="Obs_amm_io_ratio"),
		which(annualtargetdata[,4]=="Obs_nit_io_ratio"),
		which(annualtargetdata[,4]=="Obs_phyt_io_ratio"),
		which(annualtargetdata[,4]=="Obs_omni_io_ratio"),
		which(annualtargetdata[,4]=="Obs_carn_io_ratio"),
		which(annualtargetdata[,4]=="Obs_pfish_io_ratio"),
		which(annualtargetdata[,4]=="Obs_dfish_io_ratio"),
		which(annualtargetdata[,4]=="Obs_bird_io_ratio"),
		which(annualtargetdata[,4]=="Obs_seal_io_ratio"),
		which(annualtargetdata[,4]=="Obs_ceta_io_ratio")
	)


#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	#Function to make each plot panel

	makeplot_aco<-function(credintervaldata1,credintervaldata2,family,ymi,yma,axlabel,lgplot=TRUE){	

	#Set y axis limits
	
	#Barrier to lower limits
	lowbarrier<-ymi/100

	#Global minimum of both data sets
			min12<-min(c(min(credintervaldata1[2,family],na.rm=TRUE),min(credintervaldata2[6,family],na.rm=TRUE) ))
			if(min12>lowbarrier){
				yaxmin<-min(ymi,min12,na.rm=TRUE)
			} else {
				yaxmin<-lowbarrier
			}

	#Global maximum of both data sets
			max12<-max(c(max(credintervaldata1[6,family],na.rm=TRUE),max(credintervaldata2[6,family],na.rm=TRUE) ))
			yaxmax<-max(yma,max12,na.rm=TRUE)

	plot_ci_data_aco(credintervaldata1,
			yaxmin,yaxmax,family,lgplot=lgplot,overplot=FALSE,txtscale=1)

	mtext(axlabel,cex=0.75,side=1,line=1.7)

	plot_ci_data_aco(credintervaldata2,
			yaxmin,yaxmax,family,lgplot=lgplot,overplot=TRUE,txtscale=1)

	}

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

#Function to plot or over-plot the modelci data as box and whiskers
plot_ci_data_aco<-function(statdata, yaxmin,yaxmax,var.list,lgplot=TRUE,overplot=FALSE,txtscale){

   		tmpmodeldata2plot<-(statdata[2:6,var.list])

		#Reset any data less than or equal to zero to 0.1*preferred axis minimum
   		for(jj in 1:length(var.list)) {
		tmpmodeldata2plot[(which(tmpmodeldata2plot[,jj]<=0)),jj] <- yaxmin*0.1
		}
		
		modeldata2plot<-tmpmodeldata2plot[,1]
   		for(jj in 2:length(var.list)) {
			modeldata2plot<-c(modeldata2plot,tmpmodeldata2plot[,jj])
		}

   		if(lgplot==TRUE){
			 modeldata2plot<-log10(modeldata2plot)
		}
   		array2plot<- array(dim=c(5,length(var.list)),modeldata2plot)
		colnames(array2plot)<-colnames(statdata)[var.list]
   		bxpdata<-list(stats=array2plot,n=rep(100,length(var.list)),conf=NULL,out=numeric(length=0),names=colnames(statdata)[var.list])

	ymin<-yaxmin
	ymax<-yaxmax
	if(lgplot==TRUE){
		ymin<-log10(yaxmin)
		ymax<-log10(yaxmax)
	}
  
	if(overplot==FALSE){
		bxp(bxpdata,boxwex=0.3,at=seq(1,length(var.list)),ylim=c(ymin,ymax),
            		show.names=TRUE,horizontal=TRUE,las=1,boxcol="black",whiskcol="black",whisklty="solid",medcol="black",staplecol="black",
			yaxt="n",xaxt="n",cex=0.4,cex.lab=txtscale,cex.axis=txtscale)
			axis(1,padj=-0.75)
		nl<-length(colnames(statdata)[var.list])
		labs<-colnames(statdata)[var.list]
		axis(side=2, at=(seq(1,nl,2))+0.25, las=1,labels=labs[seq(1,nl,2)], cex.axis=0.9)
		axis(side=2, at=(seq(2,nl,2))+0.25, las=1,labels=labs[seq(2,nl,2)], cex.axis=0.9)
	}

	if(overplot==TRUE){
		bxp(bxpdata,add=TRUE,boxwex=0.35,at=seq(1,length(var.list))+0.35,yaxt="n",xaxt="n",horizontal=TRUE,boxcol="red",whiskcol="red",whisklty="solid",medcol="red",staplecol="red",
			cex=0.4,cex.lab=txtscale,cex.axis=txtscale)
	}
}

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~




#Now the actual plotting

	par(mfrow=c(3,2))
	par(mar=c(3,12,0.4,2.75))

#..................................................................


	if(length(set2plot1.1)>0){

	axlabel<-bquote("Production log"[10] ~ "mMN.m"^-2 * ".y"^-1)
	makeplot_aco(credintervaldata1, credintervaldata2, set2plot1.1, 0.001, 2000, axlabel,lgplot=TRUE)

	}
#..................................................................

	if(length(set2plot1.2)>0){

	axlabel<-bquote("Catch log"[10] ~ "mMN.m"^-2 * ".y"^-1)
	makeplot_aco(credintervaldata1, credintervaldata2, set2plot1.2, 0.00001, 100, axlabel,lgplot=TRUE)

	}
#..................................................................

	if(length(set2plot2)>0){

	axlabel<-bquote("Consumption log"[10] ~ "mMN.m"^-2 * ".y"^-1)
	makeplot_aco(credintervaldata1, credintervaldata2, set2plot2, 0.0001, 200, axlabel,lgplot=TRUE)

	}
#..................................................................

	if(length(set2plot3)>0){

	axlabel<-bquote("log"[10] ~ "Annual ratio")
	makeplot_aco(credintervaldata1, credintervaldata2, set2plot3, 0.01, 100, axlabel,lgplot=TRUE)

	}
#..................................................................

	if(length(set2plot4)>0){

	axlabel<-bquote("Conc. log"[10] ~ "mMN.m"^-3 ~ "or % by weight")
	makeplot_aco(credintervaldata1, credintervaldata2, set2plot4, 0.01, 1000, axlabel,lgplot=TRUE)

	}
#..................................................................


	if(length(set2plot5)>0){

	axlabel<-bquote("Inshore:offshore ratios")
	makeplot_aco(credintervaldata1, credintervaldata2, set2plot5, 0, 3, axlabel,lgplot=FALSE)

	}

	legend(grconvertX(0.45, "ndc", "user"), grconvertY(0.13, "ndc", "user"),
	"scenario", fill = "red", ncol=1, bty="n", xpd = NA)
	legend(grconvertX(0.45, "ndc", "user"), grconvertY(0.10, "ndc", "user"),
	"baseline", fill = "black", ncol=1, bty="n", xpd = NA)

#..................................................................


}

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.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close