R/boxplot_monthly_compare_observations.R

Defines functions boxplot_monthly_compare_observations

#
# boxplot_monthly_compare_observations.R
#
#' Box-and-whisker plots of observed monthly data on the state of the ecosystem with equivalent properties derived from a StrathE2E run.
#'
#' Multi-panel plot comparing annual cycles of monthly averaged observational data on nutrinet and plankton concentrations with equivalent data derived from the final year of a model run generated by the e2e_run() function.
#'
#' The function plots a multi-panel page of box-and-whisker plots showing the medians and variability ranges (quartiles as box-and-whisker) of observational data on nitrate, ammonia, chlorophyll, omnivorous and carnivorous zooplankton, and benthos larvae concentrations
#' alongside single-value data (in red) on equivalent measures derived from the final year of a e2e_run() fuunction call.
#'
#' The observational data to be plotted are loaded from the folder Modelname/Variantname/Target/monthly_observed_*.csv as part of a e2e_read() function call and are built into the R-list object generated by e2e_read()
#'
#' The corresponding measures derived from the final year of a model run generated by the e2e_run() function call are located in /results/Modelname/Variantname/model_target_monthlyresults-*.csv, and in the R-list object generated by StrathE2E()
#'
#' @param model R-list object defining the model configuration compiled by the e2e_read() function.
#' @param use.saved Logical. If TRUE use data from a prior user-defined run held as csv files data in the current results folder, (default=FALSE).
#' @param results R-list object of baseline 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_monthly_compare_observations <- function(model, use.saved=FALSE,results=NULL) {

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

if(use.saved==TRUE & is.list(results)==TRUE){
        stop("use.saved is TRUE but a dataframe object has also been specified ! \n")
}

#       Need to access raw monthly target data from the model setup in the event that the option to
#       plot data from a saved output file is selected. This is because the monthly target data are
#       not saved anywhere in the set of csv output files generated by a model run, though they are
#       in the results list object
	setup		<- elt(model, "setup")
	model.path	<- elt(setup, "model.path")

if(use.saved==TRUE) {
	resultsdir	<- elt(model, "setup", "resultsdir")
	model.ident	<- elt(model, "setup", "model.ident")
	datafile	<- csvname(resultsdir, "model_monthlyresults", model.ident)
	print(paste("Using data held in a file ",datafile," from a past model run"))
	check.exists(datafile)
	modelmonthly<-readcsv(datafile)
	obstargetdataset<-read_monthly_target_data(model.path)
}

if(use.saved==FALSE & is.list(results)==TRUE){
	final.year.outputs	<- elt(results, "final.year.outputs")
	modelmonthly		<- elt(final.year.outputs, "monthly.averages")
	obstargetdataset	<- elt(final.year.outputs, "monthly.target.data")
}

#Structure of the observed target data
#Format expected = 7 columns
#Month	Variable	median	lower_centile	upper_centile	Units	low_cent_value	upp_cent_value	Comments
#The variable names expected are:
#surface_nitrate
#deep_nitrate
#surface_ammonia
#deep_ammonia
#surface_chlorophyll
#omniv_zooplankton
#carniv_zooplankton
#larvae_susp_dep_benthos
#larvae_carn_scav_benthos

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

#Make an array to hold the synthesised observational monthly data
ntargobs<-100
monthtarget<-array(0,dim=c(ntargobs,12,10))
#dimensions observations,month, parameter, 

#Make an array to hold the modelled monthly data
nmodobs<-10
MODmonthtarget<-array(0,dim=c(nmodobs,12,10))
#dimensions observations,month, parameter, 
MODmonthtarget[1:nmodobs,1,1]<-modelmonthly$surfnitratemMm3[1]
MODmonthtarget[1:nmodobs,2,1]<-modelmonthly$surfnitratemMm3[2]
MODmonthtarget[1:nmodobs,3,1]<-modelmonthly$surfnitratemMm3[3]
MODmonthtarget[1:nmodobs,4,1]<-modelmonthly$surfnitratemMm3[4]
MODmonthtarget[1:nmodobs,5,1]<-modelmonthly$surfnitratemMm3[5]
MODmonthtarget[1:nmodobs,6,1]<-modelmonthly$surfnitratemMm3[6]
MODmonthtarget[1:nmodobs,7,1]<-modelmonthly$surfnitratemMm3[7]
MODmonthtarget[1:nmodobs,8,1]<-modelmonthly$surfnitratemMm3[8]
MODmonthtarget[1:nmodobs,9,1]<-modelmonthly$surfnitratemMm3[9]
MODmonthtarget[1:nmodobs,10,1]<-modelmonthly$surfnitratemMm3[10]
MODmonthtarget[1:nmodobs,11,1]<-modelmonthly$surfnitratemMm3[11]
MODmonthtarget[1:nmodobs,12,1]<-modelmonthly$surfnitratemMm3[12]

MODmonthtarget[1:nmodobs,1,2]<-modelmonthly$deepnitratemMm3[1]
MODmonthtarget[1:nmodobs,2,2]<-modelmonthly$deepnitratemMm3[2]
MODmonthtarget[1:nmodobs,3,2]<-modelmonthly$deepnitratemMm3[3]
MODmonthtarget[1:nmodobs,4,2]<-modelmonthly$deepnitratemMm3[4]
MODmonthtarget[1:nmodobs,5,2]<-modelmonthly$deepnitratemMm3[5]
MODmonthtarget[1:nmodobs,6,2]<-modelmonthly$deepnitratemMm3[6]
MODmonthtarget[1:nmodobs,7,2]<-modelmonthly$deepnitratemMm3[7]
MODmonthtarget[1:nmodobs,8,2]<-modelmonthly$deepnitratemMm3[8]
MODmonthtarget[1:nmodobs,9,2]<-modelmonthly$deepnitratemMm3[9]
MODmonthtarget[1:nmodobs,10,2]<-modelmonthly$deepnitratemMm3[10]
MODmonthtarget[1:nmodobs,11,2]<-modelmonthly$deepnitratemMm3[11]
MODmonthtarget[1:nmodobs,12,2]<-modelmonthly$deepnitratemMm3[12]

MODmonthtarget[1:nmodobs,1,3]<-modelmonthly$surfammoniamMm3[1]
MODmonthtarget[1:nmodobs,2,3]<-modelmonthly$surfammoniamMm3[2]
MODmonthtarget[1:nmodobs,3,3]<-modelmonthly$surfammoniamMm3[3]
MODmonthtarget[1:nmodobs,4,3]<-modelmonthly$surfammoniamMm3[4]
MODmonthtarget[1:nmodobs,5,3]<-modelmonthly$surfammoniamMm3[5]
MODmonthtarget[1:nmodobs,6,3]<-modelmonthly$surfammoniamMm3[6]
MODmonthtarget[1:nmodobs,7,3]<-modelmonthly$surfammoniamMm3[7]
MODmonthtarget[1:nmodobs,8,3]<-modelmonthly$surfammoniamMm3[8]
MODmonthtarget[1:nmodobs,9,3]<-modelmonthly$surfammoniamMm3[9]
MODmonthtarget[1:nmodobs,10,3]<-modelmonthly$surfammoniamMm3[10]
MODmonthtarget[1:nmodobs,11,3]<-modelmonthly$surfammoniamMm3[11]
MODmonthtarget[1:nmodobs,12,3]<-modelmonthly$surfammoniamMm3[12]

MODmonthtarget[1:nmodobs,1,4]<-modelmonthly$deepammoniamMm3[1]
MODmonthtarget[1:nmodobs,2,4]<-modelmonthly$deepammoniamMm3[2]
MODmonthtarget[1:nmodobs,3,4]<-modelmonthly$deepammoniamMm3[3]
MODmonthtarget[1:nmodobs,4,4]<-modelmonthly$deepammoniamMm3[4]
MODmonthtarget[1:nmodobs,5,4]<-modelmonthly$deepammoniamMm3[5]
MODmonthtarget[1:nmodobs,6,4]<-modelmonthly$deepammoniamMm3[6]
MODmonthtarget[1:nmodobs,7,4]<-modelmonthly$deepammoniamMm3[7]
MODmonthtarget[1:nmodobs,8,4]<-modelmonthly$deepammoniamMm3[8]
MODmonthtarget[1:nmodobs,9,4]<-modelmonthly$deepammoniamMm3[9]
MODmonthtarget[1:nmodobs,10,4]<-modelmonthly$deepammoniamMm3[10]
MODmonthtarget[1:nmodobs,11,4]<-modelmonthly$deepammoniamMm3[11]
MODmonthtarget[1:nmodobs,12,4]<-modelmonthly$deepammoniamMm3[12]

MODmonthtarget[1:nmodobs,1,5]<-modelmonthly$surfchlmgm3[1]
MODmonthtarget[1:nmodobs,2,5]<-modelmonthly$surfchlmgm3[2]
MODmonthtarget[1:nmodobs,3,5]<-modelmonthly$surfchlmgm3[3]
MODmonthtarget[1:nmodobs,4,5]<-modelmonthly$surfchlmgm3[4]
MODmonthtarget[1:nmodobs,5,5]<-modelmonthly$surfchlmgm3[5]
MODmonthtarget[1:nmodobs,6,5]<-modelmonthly$surfchlmgm3[6]
MODmonthtarget[1:nmodobs,7,5]<-modelmonthly$surfchlmgm3[7]
MODmonthtarget[1:nmodobs,8,5]<-modelmonthly$surfchlmgm3[8]
MODmonthtarget[1:nmodobs,9,5]<-modelmonthly$surfchlmgm3[9]
MODmonthtarget[1:nmodobs,10,5]<-modelmonthly$surfchlmgm3[10]
MODmonthtarget[1:nmodobs,11,5]<-modelmonthly$surfchlmgm3[11]
MODmonthtarget[1:nmodobs,12,5]<-modelmonthly$surfchlmgm3[12]

MODmonthtarget[1:nmodobs,1,6]<-modelmonthly$omnizoomMNm3[1]
MODmonthtarget[1:nmodobs,2,6]<-modelmonthly$omnizoomMNm3[2]
MODmonthtarget[1:nmodobs,3,6]<-modelmonthly$omnizoomMNm3[3]
MODmonthtarget[1:nmodobs,4,6]<-modelmonthly$omnizoomMNm3[4]
MODmonthtarget[1:nmodobs,5,6]<-modelmonthly$omnizoomMNm3[5]
MODmonthtarget[1:nmodobs,6,6]<-modelmonthly$omnizoomMNm3[6]
MODmonthtarget[1:nmodobs,7,6]<-modelmonthly$omnizoomMNm3[7]
MODmonthtarget[1:nmodobs,8,6]<-modelmonthly$omnizoomMNm3[8]
MODmonthtarget[1:nmodobs,9,6]<-modelmonthly$omnizoomMNm3[9]
MODmonthtarget[1:nmodobs,10,6]<-modelmonthly$omnizoomMNm3[10]
MODmonthtarget[1:nmodobs,11,6]<-modelmonthly$omnizoomMNm3[11]
MODmonthtarget[1:nmodobs,12,6]<-modelmonthly$omnizoomMNm3[12]

MODmonthtarget[1:nmodobs,1,7]<-modelmonthly$carnzoomMNm3[1]
MODmonthtarget[1:nmodobs,2,7]<-modelmonthly$carnzoomMNm3[2]
MODmonthtarget[1:nmodobs,3,7]<-modelmonthly$carnzoomMNm3[3]
MODmonthtarget[1:nmodobs,4,7]<-modelmonthly$carnzoomMNm3[4]
MODmonthtarget[1:nmodobs,5,7]<-modelmonthly$carnzoomMNm3[5]
MODmonthtarget[1:nmodobs,6,7]<-modelmonthly$carnzoomMNm3[6]
MODmonthtarget[1:nmodobs,7,7]<-modelmonthly$carnzoomMNm3[7]
MODmonthtarget[1:nmodobs,8,7]<-modelmonthly$carnzoomMNm3[8]
MODmonthtarget[1:nmodobs,9,7]<-modelmonthly$carnzoomMNm3[9]
MODmonthtarget[1:nmodobs,10,7]<-modelmonthly$carnzoomMNm3[10]
MODmonthtarget[1:nmodobs,11,7]<-modelmonthly$carnzoomMNm3[11]
MODmonthtarget[1:nmodobs,12,7]<-modelmonthly$carnzoomMNm3[12]

MODmonthtarget[1:nmodobs,1,8]<-modelmonthly$benthslarmMNm3[1]
MODmonthtarget[1:nmodobs,2,8]<-modelmonthly$benthslarmMNm3[2]
MODmonthtarget[1:nmodobs,3,8]<-modelmonthly$benthslarmMNm3[3]
MODmonthtarget[1:nmodobs,4,8]<-modelmonthly$benthslarmMNm3[4]
MODmonthtarget[1:nmodobs,5,8]<-modelmonthly$benthslarmMNm3[5]
MODmonthtarget[1:nmodobs,6,8]<-modelmonthly$benthslarmMNm3[6]
MODmonthtarget[1:nmodobs,7,8]<-modelmonthly$benthslarmMNm3[7]
MODmonthtarget[1:nmodobs,8,8]<-modelmonthly$benthslarmMNm3[8]
MODmonthtarget[1:nmodobs,9,8]<-modelmonthly$benthslarmMNm3[9]
MODmonthtarget[1:nmodobs,10,8]<-modelmonthly$benthslarmMNm3[10]
MODmonthtarget[1:nmodobs,11,8]<-modelmonthly$benthslarmMNm3[11]
MODmonthtarget[1:nmodobs,12,8]<-modelmonthly$benthslarmMNm3[12]

MODmonthtarget[1:nmodobs,1,9]<-modelmonthly$benthclarmMNm3[1]
MODmonthtarget[1:nmodobs,2,9]<-modelmonthly$benthclarmMNm3[2]
MODmonthtarget[1:nmodobs,3,9]<-modelmonthly$benthclarmMNm3[3]
MODmonthtarget[1:nmodobs,4,9]<-modelmonthly$benthclarmMNm3[4]
MODmonthtarget[1:nmodobs,5,9]<-modelmonthly$benthclarmMNm3[5]
MODmonthtarget[1:nmodobs,6,9]<-modelmonthly$benthclarmMNm3[6]
MODmonthtarget[1:nmodobs,7,9]<-modelmonthly$benthclarmMNm3[7]
MODmonthtarget[1:nmodobs,8,9]<-modelmonthly$benthclarmMNm3[8]
MODmonthtarget[1:nmodobs,9,9]<-modelmonthly$benthclarmMNm3[9]
MODmonthtarget[1:nmodobs,10,9]<-modelmonthly$benthclarmMNm3[10]
MODmonthtarget[1:nmodobs,11,9]<-modelmonthly$benthclarmMNm3[11]
MODmonthtarget[1:nmodobs,12,9]<-modelmonthly$benthclarmMNm3[12]

MODmonthtarget[1:nmodobs,1:12,10]<-MODmonthtarget[1:nmodobs,1:12,8]+MODmonthtarget[1:nmodobs,1:12,9]
#Combine the s/d and c/s benthos




monlab<-c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")

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

	installtargdata<-function(qr,obspar,monthtarget){
		#Make sure that the selected data are in the correct order
		data<-data.frame("Month"=monlab,"median"=rep(NA,12),"17centile"=rep(NA,12),"83centile"=rep(NA,12))
		for(mmm in 1:12){
			data[mmm,2:4]<-seldata[(which(seldata$Month==data$Month[mmm])),3:5]
		}
		for(mmm in 1:12){
			sdpos<-(data[mmm,4]-data[mmm,2])/(qr/2)
			sdneg<-(data[mmm,2]-data[mmm,3])/(qr/2)
			for(kkk in 1:ntargobs){
				rand<-rnorm(1,0,1)
				if(rand<0) dev<-rand*sdneg
				if(rand>=0) dev<-rand*sdpos
				monthtarget[kkk,mmm,obspar]<-data[mmm,2]+dev
			}
		}

		return(monthtarget)
	}
	#~~~~~~~~~~~~~




#~~~~~~~~~~~~~
#Read in the monthly surface nitrate target data file
#qr=number of sd in the observed quantile range
#obspar=position of the data in the data array
seldata<-subset(obstargetdataset,obstargetdataset$Variable=="surface_nitrate")
obspar<-1
qr<-3.5
if(seldata[1,7]==17 & seldata[1,8]==83) qr<-3.5
if(seldata[1,7]==5 & seldata[1,8]==95) qr<-7
monthtarget<-installtargdata(qr,obspar,monthtarget)


#~~~~~~~~~~~~~
#Read in the monthly deep nitrate target data file
#qr=number of sd in the observed quantile range
#obspar=position of the data in the data array
seldata<-subset(obstargetdataset,obstargetdataset$Variable=="deep_nitrate")
obspar<-2
qr<-3.5
if(seldata[1,7]==17 & seldata[1,8]==83) qr<-3.5
if(seldata[1,7]==5 & seldata[1,8]==95) qr<-7
monthtarget<-installtargdata(qr,obspar,monthtarget)


#~~~~~~~~~~~~~
#Read in the monthly surface ammonia target data file
#qr=number of sd in the observed quantile range
#obspar=position of the data in the data array
seldata<-subset(obstargetdataset,obstargetdataset$Variable=="surface_ammonia")
obspar<-3
qr<-3.5
if(seldata[1,7]==17 & seldata[1,8]==83) qr<-3.5
if(seldata[1,7]==5 & seldata[1,8]==95) qr<-7
monthtarget<-installtargdata(qr,obspar,monthtarget)


#~~~~~~~~~~~~~
#Read in the monthly deep ammonia target data file
#qr=number of sd in the observed quantile range
#obspar=position of the data in the data array
seldata<-subset(obstargetdataset,obstargetdataset$Variable=="deep_ammonia")
obspar<-4
qr<-3.5
if(seldata[1,7]==17 & seldata[1,8]==83) qr<-3.5
if(seldata[1,7]==5 & seldata[1,8]==95) qr<-7
monthtarget<-installtargdata(qr,obspar,monthtarget)


#~~~~~~~~~~~~~
#Read in the monthly surface chlorophyll target data file
#qr=number of sd in the observed quantile range
#obspar=position of the data in the data array
seldata<-subset(obstargetdataset,obstargetdataset$Variable=="surface_chlorophyll")
obspar<-5
qr<-3.5
if(seldata[1,7]==17 & seldata[1,8]==83) qr<-3.5
if(seldata[1,7]==5 & seldata[1,8]==95) qr<-7
monthtarget<-installtargdata(qr,obspar,monthtarget)


#~~~~~~~~~~~~~
#Read in the monthly omniv zooplankton target data file
#qr=number of sd in the observed quantile range
#obspar=position of the data in the data array
seldata<-subset(obstargetdataset,obstargetdataset$Variable=="omniv_zooplankton")
obspar<-6
qr<-3.5
if(seldata[1,7]==17 & seldata[1,8]==83) qr<-3.5
if(seldata[1,7]==5 & seldata[1,8]==95) qr<-7
monthtarget<-installtargdata(qr,obspar,monthtarget)


#~~~~~~~~~~~~~
#Read in the monthly carniv zooplankton target data file
#qr=number of sd in the observed quantile range
#obspar=position of the data in the data array
seldata<-subset(obstargetdataset,obstargetdataset$Variable=="carniv_zooplankton")
obspar<-7
qr<-3.5
if(seldata[1,7]==17 & seldata[1,8]==83) qr<-3.5
if(seldata[1,7]==5 & seldata[1,8]==95) qr<-7
monthtarget<-installtargdata(qr,obspar,monthtarget)


#~~~~~~~~~~~~~
#Read in the monthly larvae of susp/dep feeding benthos target data file
#qr=number of sd in the observed quantile range
#obspar=position of the data in the data array
seldata<-subset(obstargetdataset,obstargetdataset$Variable=="larvae_susp_dep_benthos")
obspar<-8
qr<-3.5
if(seldata[1,7]==17 & seldata[1,8]==83) qr<-3.5
if(seldata[1,7]==5 & seldata[1,8]==95) qr<-7
monthtarget<-installtargdata(qr,obspar,monthtarget)


#~~~~~~~~~~~~~
#Read in the monthly larvae of carn/scav feeding benthos target data file
#qr=number of sd in the observed quantile range
#obspar=position of the data in the data array
seldata<-subset(obstargetdataset,obstargetdataset$Variable=="larvae_carn_scav_benthos")
obspar<-9
qr<-3.5
if(seldata[1,7]==17 & seldata[1,8]==83) qr<-3.5
if(seldata[1,7]==5 & seldata[1,8]==95) qr<-7
monthtarget<-installtargdata(qr,obspar,monthtarget)



#~~~~~~~~~~~~~
#COMBINE THE TWO TYPES OF BENTHIC LARVAE
#Read in the monthly susp benlar target data file
seldata<-subset(obstargetdataset,obstargetdataset$Variable=="larvae_susp_dep_benthos")
#Read in the monthly carn benlar target data file
seldata2<-subset(obstargetdataset,obstargetdataset$Variable=="larvae_carn_scav_benthos")
seldata[,3]<-seldata[,3]+seldata2[,3]
seldata[,4]<-seldata[,4]+seldata2[,4]
seldata[,5]<-seldata[,5]+seldata2[,5]
obspar<-10
qr<-3.5
if(seldata[1,7]==17 & seldata[1,8]==83) qr<-3.5
if(seldata[1,7]==5 & seldata[1,8]==95) qr<-7
monthtarget<-installtargdata(qr,obspar,monthtarget)



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	plotdata_mso<-function(monthtarget, MODmonthtarget, obspar, monlab){

	obsplot<-as.data.frame(monthtarget[,,obspar])
	names(obsplot)<-monlab
	modplot<-as.data.frame(MODmonthtarget[,,obspar])



	# Tricky bit - need to find global maximum across obs and model when there may or may not be some NAs in the obs

	#Is there a mixture of data and NA in the observations?
	#Check if sum with na.rm=TRUE is non zero
	obschecksum_0<-sum(obsplot,na.rm=TRUE)
	obschecksum1_0<-sum(monthtarget[,,1],na.rm=TRUE)
	obschecksum2_0<-sum(monthtarget[,,2],na.rm=TRUE)
	obschecksum3_0<-sum(monthtarget[,,3],na.rm=TRUE)
	obschecksum4_0<-sum(monthtarget[,,4],na.rm=TRUE)
	obschecksum_1<-sum(obsplot)
	obschecksum1_1<-sum(monthtarget[,,1])
	obschecksum2_1<-sum(monthtarget[,,2])
	obschecksum3_1<-sum(monthtarget[,,3])
	obschecksum4_1<-sum(monthtarget[,,4])

	MIXFLAG<-0
	MIXFLAG1<-0
	MIXFLAG2<-0
	MIXFLAG3<-0
	MIXFLAG4<-0

	if(obschecksum_0>0 & is.na(obschecksum_1)==TRUE) MIXFLAG<-1    # means there is a mixture of data and NAs in the data
	if(obschecksum1_0>0 & is.na(obschecksum1_1)==TRUE) MIXFLAG1<-1
	if(obschecksum2_0>0 & is.na(obschecksum2_1)==TRUE) MIXFLAG2<-1
	if(obschecksum3_0>0 & is.na(obschecksum3_1)==TRUE) MIXFLAG3<-1
	if(obschecksum4_0>0 & is.na(obschecksum4_1)==TRUE) MIXFLAG4<-1

	if(obspar==1 | obspar==2){
	if(MIXFLAG1==1 & MIXFLAG2==1){
		ymax<- max(0, max(as.data.frame(monthtarget[,,1]),na.rm=TRUE), max(as.data.frame(monthtarget[,,2]),na.rm=TRUE), max(as.data.frame(MODmonthtarget[,,1])), max(as.data.frame(MODmonthtarget[,,2])),na.rm=TRUE )
	} else if(MIXFLAG1==1 & MIXFLAG2==0){
		ymax<- max(0, max(as.data.frame(monthtarget[,,1]),na.rm=TRUE), max(as.data.frame(monthtarget[,,2])), max(as.data.frame(MODmonthtarget[,,1])), max(as.data.frame(MODmonthtarget[,,2])),na.rm=TRUE )
	} else if(MIXFLAG1==0 & MIXFLAG2==1){
		ymax<- max(0, max(as.data.frame(monthtarget[,,1])), max(as.data.frame(monthtarget[,,2]),na.rm=TRUE), max(as.data.frame(MODmonthtarget[,,1])), max(as.data.frame(MODmonthtarget[,,2])),na.rm=TRUE )
	} else if(MIXFLAG1==0 & MIXFLAG2==0){
		ymax<- max(0, max(as.data.frame(monthtarget[,,1])), max(as.data.frame(monthtarget[,,2])), max(as.data.frame(MODmonthtarget[,,1])), max(as.data.frame(MODmonthtarget[,,2])),na.rm=TRUE )
	}
	}

	if(obspar==3 | obspar==4){
	if(MIXFLAG3==1 & MIXFLAG4==1){
		ymax<- max(0, max(as.data.frame(monthtarget[,,3]),na.rm=TRUE), max(as.data.frame(monthtarget[,,4]),na.rm=TRUE), max(as.data.frame(MODmonthtarget[,,3])), max(as.data.frame(MODmonthtarget[,,4])),na.rm=TRUE )
	} else if(MIXFLAG3==1 & MIXFLAG4==0){
		ymax<- max(0, max(as.data.frame(monthtarget[,,3]),na.rm=TRUE), max(as.data.frame(monthtarget[,,4])), max(as.data.frame(MODmonthtarget[,,3])), max(as.data.frame(MODmonthtarget[,,4])),na.rm=TRUE )
	} else if(MIXFLAG3==0 & MIXFLAG4==1){
		ymax<- max(0, max(as.data.frame(monthtarget[,,3])), max(as.data.frame(monthtarget[,,4]),na.rm=TRUE), max(as.data.frame(MODmonthtarget[,,3])), max(as.data.frame(MODmonthtarget[,,4])),na.rm=TRUE )
	} else if(MIXFLAG3==0 & MIXFLAG4==0){
		ymax<- max(0, max(as.data.frame(monthtarget[,,3])), max(as.data.frame(monthtarget[,,4])), max(as.data.frame(MODmonthtarget[,,3])), max(as.data.frame(MODmonthtarget[,,4])),na.rm=TRUE )
	}
	}

	if(obspar>4){
	if(MIXFLAG==1){
		ymax<- max(0, max(obsplot,na.rm=TRUE),max(as.data.frame(modplot)),na.rm=TRUE )
	} else if(MIXFLAG==0){
		ymax<- max(0, max(obsplot),max(as.data.frame(modplot)),na.rm=TRUE )
	}
	}

	if(ymax==0 | is.na(ymax)==TRUE) ymax<-0.1

		boxplot(obsplot,range=0,boxwex=0.25,ylim=c(0,ymax*1.1),las=1,cex.axis=1.1,yaxt="n",xaxt="n")

		axis(labels=monlab, at=seq(1,12),side=1,las=1,cex.axis=1.1,padj=-0.55)

		if(obspar==1){
			axis(side=2,cex.lab=1.0,las=1)
			mtext("Surf.nitrate",cex=0.8,side=2,line=4)
			mtext(bquote(mMN.m^-3),cex=0.6,side=2,line=2.7)
		}
		
		if(obspar==2){
			axis(side=2,cex.lab=1.0,las=1)
			mtext("Deep nitrate",cex=0.8,side=2,line=4)
			mtext(bquote(mMN.m^-3),cex=0.6,side=2,line=2.7)
		}
		
		if(obspar==3){
			axis(side=2,cex.lab=1.0,las=1)
			mtext("Surf.ammonia",cex=0.8,side=2,line=4)
			mtext(bquote(mMN.m^-3),cex=0.6,side=2,line=2.7)
		}
		
		if(obspar==4){
			axis(side=2,cex.lab=1.0,las=1)
			mtext("Deep ammonia",cex=0.8,side=2,line=4)
			mtext(bquote(mMN.m^-3),cex=0.6,side=2,line=2.7)
		}
		
		if(obspar==5){
			axis(side=2,cex.lab=1.0,las=1)
			mtext("Chlorophyll",cex=0.8,side=2,line=4)
			mtext(bquote(mg.m^-3),cex=0.6,side=2,line=2.7)
		}
		
		if(obspar==6){
			axis(side=2,cex.lab=1.0,las=1)
			mtext("Omniv.zoo",cex=0.8,side=2,line=4)
			mtext(bquote(mMN.m^-3),cex=0.6,side=2,line=2.7)
		}
		
		if(obspar==7){
			axis(side=2,cex.lab=1.0,las=1)
			mtext("Carniv.zoo",cex=0.8,side=2,line=4)
			mtext(bquote(mMN.m^-3),cex=0.6,side=2,line=2.7)
		}
		
		if(obspar==8){
			axis(side=2,cex.lab=1.0,las=1)
			mtext("Larv.s/d.benth.",cex=0.8,side=2,line=4)
			mtext(bquote(mMN.m^-3),cex=0.6,side=2,line=2.7)
		}
		
		if(obspar==9){
			axis(side=2,cex.lab=1.0,las=1)
			mtext("Larv.c/s.benth.",cex=0.8,side=2,line=4)
			mtext(bquote(mMN.m^-3),cex=0.6,side=2,line=2.7)
		}
		
		if(obspar==10){
			axis(side=2,cex.lab=1.0,las=1)
			mtext("Benthos larvae (all)",cex=0.8,side=2,line=4)
			mtext(bquote(mMN.m^-3),cex=0.6,side=2,line=2.7)
		}


		boxplot(modplot,add=TRUE,range=0,boxwex=0.25,at=1:12+0.3,par(lty=1),yaxt="n",xaxt="n",
		boxcol="red",whiskcol="red",whisklty="solid",medcol="red",staplecol="red")

	}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



par(mfrow=c(4,2))

par(mar=c(3,6,0.6,0.5))


plotdata_mso(monthtarget, MODmonthtarget, 1, monlab)

plotdata_mso(monthtarget, MODmonthtarget, 2, monlab)

plotdata_mso(monthtarget, MODmonthtarget, 3, monlab)

plotdata_mso(monthtarget, MODmonthtarget, 4, monlab)

plotdata_mso(monthtarget, MODmonthtarget, 5, monlab)

plotdata_mso(monthtarget, MODmonthtarget, 6, monlab)

plotdata_mso(monthtarget, MODmonthtarget, 7, monlab)

#plotdata_mso(monthtarget, MODmonthtarget, 8, monlab)

#plotdata_mso(monthtarget, MODmonthtarget, 9, monlab)

plotdata_mso(monthtarget, MODmonthtarget, 10, monlab)


	legend(grconvertX(0.425, "ndc", "user"), grconvertY(0.045, "ndc", "user"),
	c("observations","model"), fill = c("black","red"), ncol=2, 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.