Nothing
R/boxplot_monthly_compare_observations_with_ci.R
In StrathE2E2: End-to-End Marine Food Web Model
Defines functions
boxplot_monthly_compare_observations_with_ci
#
# boxplot_monthly_compare_observations_with_ci.R
#
#' Box-and-whisker plots of observed monthly data on the state of the ecosystem with equivalent properties derived from a Monte Carlo StrathE2E run.
#'
#' Creates a multi-panel plot comparing annual cycles of monthly averaged observational data on nutrinet and plankton concentrations together with the distribution of credible values of equivalent data derived from the final year of a model run generated by the e2e_run_mc() 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 observational data on nitrate, ammonia, chlorophyll, omnivorous and carnivorous zooplankton, and benthos larvae concentrations
#' (shown in black), alongside comparable box-and-whisker plots (shown in red) of equivalent measures derived from the final years of an ensemble of model runs generted by a Monte Carlo methodology (e2e_run_mc() function).
#'
#' 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()
#'
#' Optionally the function can read an example data set for one of the two North Sea model variants supplied with the package.
#'
#' The corresponding measures derived from the final year of a model run generated by the e2e_run_mc() function call are located in /results/Modelname/Variantname/CredInt_processed_monthly_mass-*.csv, and in the R-list object generated by e2e_run_mc()
#'
#' @param model R-list object defining the model configuration compiled by the e2e_read() function.
#' @param use.example Logical. If TRUE use pre-computed example data from the internal North Sea model rather than user-generated data, (default=FALSE).
#'
#' @return Graphical display in a new graphics window.
#'
#' @importFrom graphics 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_with_ci <- function(model, use.example=FALSE) {
start_par = par()$mfrow
on.exit(par(mfrow = start_par))
resultsdir <- elt(model, "setup", "resultsdir")
model.path <- elt(model, "setup", "model.path")
model.ident <- elt(model, "setup", "model.ident")
model.name <- elt(model, "setup", "model.name")
model.variant <- elt(model, "setup", "model.variant")
#Read the observed data file
#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
obstargetdataset <- get.model.file(model.path, TARGET_DATA_DIR, file.pattern=MONTHLY_TARGET_DATA)
corefilename<-"CredInt_processed_monthly_mass"
if(use.example==TRUE){
credintervaldata <- get.example.results(model.name, model.variant, corefilename, CREDINT_DIR)
}
if(use.example==FALSE){
credpath <- makepath(resultsdir, CREDINT_DIR)
credfile <- csvname(credpath, corefilename, model.ident)
message("Reading credible interval processed data from '", credfile, "'")
if (! file.exists(credfile)) {
message("Error: cannot find credible interval output file: ", credfile)
stop("Please run the Monte-Carlo function!\n")
}
credintervaldata <- readcsv(credfile, row.names=1) # first column is row names
}
#Make an array to hold the synthesised observational monthly data
ntargobs<-100
monthtarget<-array(0,dim=c(ntargobs,12,10)) #dimensions observations,month, parameter,
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)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Generate the list objects needed by the bxp plotting function
for(iii in 1:9){
credrows<- seq( ((iii-1)*(5+1))+2,((iii-1)*(5+1))+(5+1) )
modeldata2plot<-(credintervaldata[credrows,1])
for(jj in 2:12) { modeldata2plot<-c(modeldata2plot,(credintervaldata[credrows,jj]))}
array2plot<- array(dim=c(5,12),modeldata2plot)
bxpdata<-list(stats=array2plot,n=rep(100,12),conf=NULL,out=numeric(length=0))
# bxp(bxpdata,boxwex=0.3,at=seq(1,12)+0.35,xlim=c(0,13),ylim=c(0,max(modeldata2plot)*1.1))
if(iii==1) bxpdata1<-bxpdata
if(iii==2) bxpdata2<-bxpdata
if(iii==3) bxpdata3<-bxpdata
if(iii==4) bxpdata4<-bxpdata
if(iii==5) bxpdata5<-bxpdata
if(iii==6) bxpdata6<-bxpdata
if(iii==7) bxpdata7<-bxpdata
if(iii==8) bxpdata8<-bxpdata
if(iii==9) bxpdata9<-bxpdata
}
bxpdata10<-list(stats=(bxpdata9$stats+bxpdata8$stats),n=rep(100,12),conf=NULL,out=numeric(length=0))
#Combines the two types of benthic larvae
#Package all the bxpdata objects up into a list t pass into the plotting function
bxpdata.all<-list(bxpdata1=bxpdata1,
bxpdata2=bxpdata2,
bxpdata3=bxpdata3,
bxpdata4=bxpdata4,
bxpdata5=bxpdata5,
bxpdata6=bxpdata6,
bxpdata7=bxpdata7,
bxpdata8=bxpdata8,
bxpdata9=bxpdata9,
bxpdata10=bxpdata10)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
plotdata_mo<-function(monthtarget, bxpdata.all, obspar, monlab){
obsplot<-as.data.frame(monthtarget[,,obspar])
names(obsplot)<-monlab
if(obspar==1) bxpdata<-bxpdata.all$bxpdata1
if(obspar==2) bxpdata<-bxpdata.all$bxpdata2
if(obspar==3) bxpdata<-bxpdata.all$bxpdata3
if(obspar==4) bxpdata<-bxpdata.all$bxpdata4
if(obspar==5) bxpdata<-bxpdata.all$bxpdata5
if(obspar==6) bxpdata<-bxpdata.all$bxpdata6
if(obspar==7) bxpdata<-bxpdata.all$bxpdata7
if(obspar==8) bxpdata<-bxpdata.all$bxpdata8
if(obspar==9) bxpdata<-bxpdata.all$bxpdata9
if(obspar==10) bxpdata<-bxpdata.all$bxpdata10
modplot<-bxpdata$stats
# 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(bxpdata.all$bxpdata1$stats)), max(as.data.frame(bxpdata.all$bxpdata2$stats)),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(bxpdata.all$bxpdata1$stats)), max(as.data.frame(bxpdata.all$bxpdata2$stats)),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(bxpdata.all$bxpdata1$stats)), max(as.data.frame(bxpdata.all$bxpdata2$stats)),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(bxpdata.all$bxpdata1$stats)), max(as.data.frame(bxpdata.all$bxpdata2$stats)),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(bxpdata.all$bxpdata3$stats)), max(as.data.frame(bxpdata.all$bxpdata4$stats)),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(bxpdata.all$bxpdata3$stats)), max(as.data.frame(bxpdata.all$bxpdata4$stats)),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(bxpdata.all$bxpdata3$stats)), max(as.data.frame(bxpdata.all$bxpdata4$stats)),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(bxpdata.all$bxpdata3$stats)), max(as.data.frame(bxpdata.all$bxpdata4$stats)),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)
}
bxp(bxpdata,add=TRUE,boxwex=0.25,at=1:12+0.35,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_mo(monthtarget, bxpdata.all, 1, monlab)
plotdata_mo(monthtarget, bxpdata.all, 2, monlab)
plotdata_mo(monthtarget, bxpdata.all, 3, monlab)
plotdata_mo(monthtarget, bxpdata.all, 4, monlab)
plotdata_mo(monthtarget, bxpdata.all, 5, monlab)
plotdata_mo(monthtarget, bxpdata.all, 6, monlab)
plotdata_mo(monthtarget, bxpdata.all, 7, monlab)
#plotdata_mo(monthtarget, bxpdata.all, 8, monlab)
#plotdata_mo(monthtarget, bxpdata.all, 9, monlab)
plotdata_mo(monthtarget, bxpdata.all, 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_with_ci
#
# boxplot_monthly_compare_observations_with_ci.R
#
#' Box-and-whisker plots of observed monthly data on the state of the ecosystem with equivalent properties derived from a Monte Carlo StrathE2E run.
#'
#' Creates a multi-panel plot comparing annual cycles of monthly averaged observational data on nutrinet and plankton concentrations together with the distribution of credible values of equivalent data derived from the final year of a model run generated by the e2e_run_mc() 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 observational data on nitrate, ammonia, chlorophyll, omnivorous and carnivorous zooplankton, and benthos larvae concentrations
#' (shown in black), alongside comparable box-and-whisker plots (shown in red) of equivalent measures derived from the final years of an ensemble of model runs generted by a Monte Carlo methodology (e2e_run_mc() function).
#'
#' 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()
#'
#' Optionally the function can read an example data set for one of the two North Sea model variants supplied with the package.
#'
#' The corresponding measures derived from the final year of a model run generated by the e2e_run_mc() function call are located in /results/Modelname/Variantname/CredInt_processed_monthly_mass-*.csv, and in the R-list object generated by e2e_run_mc()
#'
#' @param model R-list object defining the model configuration compiled by the e2e_read() function.
#' @param use.example Logical. If TRUE use pre-computed example data from the internal North Sea model rather than user-generated data, (default=FALSE).
#'
#' @return Graphical display in a new graphics window.
#'
#' @importFrom graphics 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_with_ci <- function(model, use.example=FALSE) {
start_par = par()$mfrow
on.exit(par(mfrow = start_par))
resultsdir <- elt(model, "setup", "resultsdir")
model.path <- elt(model, "setup", "model.path")
model.ident <- elt(model, "setup", "model.ident")
model.name <- elt(model, "setup", "model.name")
model.variant <- elt(model, "setup", "model.variant")
#Read the observed data file
#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
obstargetdataset <- get.model.file(model.path, TARGET_DATA_DIR, file.pattern=MONTHLY_TARGET_DATA)
corefilename<-"CredInt_processed_monthly_mass"
if(use.example==TRUE){
credintervaldata <- get.example.results(model.name, model.variant, corefilename, CREDINT_DIR)
}
if(use.example==FALSE){
credpath <- makepath(resultsdir, CREDINT_DIR)
credfile <- csvname(credpath, corefilename, model.ident)
message("Reading credible interval processed data from '", credfile, "'")
if (! file.exists(credfile)) {
message("Error: cannot find credible interval output file: ", credfile)
stop("Please run the Monte-Carlo function!\n")
}
credintervaldata <- readcsv(credfile, row.names=1) # first column is row names
}
#Make an array to hold the synthesised observational monthly data
ntargobs<-100
monthtarget<-array(0,dim=c(ntargobs,12,10)) #dimensions observations,month, parameter,
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)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Generate the list objects needed by the bxp plotting function
for(iii in 1:9){
credrows<- seq( ((iii-1)*(5+1))+2,((iii-1)*(5+1))+(5+1) )
modeldata2plot<-(credintervaldata[credrows,1])
for(jj in 2:12) { modeldata2plot<-c(modeldata2plot,(credintervaldata[credrows,jj]))}
array2plot<- array(dim=c(5,12),modeldata2plot)
bxpdata<-list(stats=array2plot,n=rep(100,12),conf=NULL,out=numeric(length=0))
# bxp(bxpdata,boxwex=0.3,at=seq(1,12)+0.35,xlim=c(0,13),ylim=c(0,max(modeldata2plot)*1.1))
if(iii==1) bxpdata1<-bxpdata
if(iii==2) bxpdata2<-bxpdata
if(iii==3) bxpdata3<-bxpdata
if(iii==4) bxpdata4<-bxpdata
if(iii==5) bxpdata5<-bxpdata
if(iii==6) bxpdata6<-bxpdata
if(iii==7) bxpdata7<-bxpdata
if(iii==8) bxpdata8<-bxpdata
if(iii==9) bxpdata9<-bxpdata
}
bxpdata10<-list(stats=(bxpdata9$stats+bxpdata8$stats),n=rep(100,12),conf=NULL,out=numeric(length=0))
#Combines the two types of benthic larvae
#Package all the bxpdata objects up into a list t pass into the plotting function
bxpdata.all<-list(bxpdata1=bxpdata1,
bxpdata2=bxpdata2,
bxpdata3=bxpdata3,
bxpdata4=bxpdata4,
bxpdata5=bxpdata5,
bxpdata6=bxpdata6,
bxpdata7=bxpdata7,
bxpdata8=bxpdata8,
bxpdata9=bxpdata9,
bxpdata10=bxpdata10)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
plotdata_mo<-function(monthtarget, bxpdata.all, obspar, monlab){
obsplot<-as.data.frame(monthtarget[,,obspar])
names(obsplot)<-monlab
if(obspar==1) bxpdata<-bxpdata.all$bxpdata1
if(obspar==2) bxpdata<-bxpdata.all$bxpdata2
if(obspar==3) bxpdata<-bxpdata.all$bxpdata3
if(obspar==4) bxpdata<-bxpdata.all$bxpdata4
if(obspar==5) bxpdata<-bxpdata.all$bxpdata5
if(obspar==6) bxpdata<-bxpdata.all$bxpdata6
if(obspar==7) bxpdata<-bxpdata.all$bxpdata7
if(obspar==8) bxpdata<-bxpdata.all$bxpdata8
if(obspar==9) bxpdata<-bxpdata.all$bxpdata9
if(obspar==10) bxpdata<-bxpdata.all$bxpdata10
modplot<-bxpdata$stats
# 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(bxpdata.all$bxpdata1$stats)), max(as.data.frame(bxpdata.all$bxpdata2$stats)),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(bxpdata.all$bxpdata1$stats)), max(as.data.frame(bxpdata.all$bxpdata2$stats)),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(bxpdata.all$bxpdata1$stats)), max(as.data.frame(bxpdata.all$bxpdata2$stats)),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(bxpdata.all$bxpdata1$stats)), max(as.data.frame(bxpdata.all$bxpdata2$stats)),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(bxpdata.all$bxpdata3$stats)), max(as.data.frame(bxpdata.all$bxpdata4$stats)),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(bxpdata.all$bxpdata3$stats)), max(as.data.frame(bxpdata.all$bxpdata4$stats)),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(bxpdata.all$bxpdata3$stats)), max(as.data.frame(bxpdata.all$bxpdata4$stats)),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(bxpdata.all$bxpdata3$stats)), max(as.data.frame(bxpdata.all$bxpdata4$stats)),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)
}
bxp(bxpdata,add=TRUE,boxwex=0.25,at=1:12+0.35,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_mo(monthtarget, bxpdata.all, 1, monlab)
plotdata_mo(monthtarget, bxpdata.all, 2, monlab)
plotdata_mo(monthtarget, bxpdata.all, 3, monlab)
plotdata_mo(monthtarget, bxpdata.all, 4, monlab)
plotdata_mo(monthtarget, bxpdata.all, 5, monlab)
plotdata_mo(monthtarget, bxpdata.all, 6, monlab)
plotdata_mo(monthtarget, bxpdata.all, 7, monlab)
#plotdata_mo(monthtarget, bxpdata.all, 8, monlab)
#plotdata_mo(monthtarget, bxpdata.all, 9, monlab)
plotdata_mo(monthtarget, bxpdata.all, 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.