Nothing
R/boxplot_monthly_compare_observations.R
In StrathE2E2: End-to-End Marine Food Web Model
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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.