#' A plotting suite for Hydra.
#'
#'These functions comprise plotting routines for each component of the data
#'
#' #############################################################################
#' A complete rework is required. Code is clunky and prone to issues when expand
#'##############################################################################
#'
#'This isnt exported with the package
#'@param inputData data passed
#'
# plots Intake, stomach content,
binWidths_Calc <- function(inputData){
#binWidths <- read.table(paste0(inPath,"length_sizebins.csv"),sep=",",header=TRUE)
binWidths <- inputData$binwidth
# finds midpoints of bins for biomass plot
#binWidths <- t(head(t(binWidths),-1))
#class(binWidths) <- "numeric"
cumSumBinWidths <- t(apply(binWidths,1,cumsum))
v <- rep(0,dim(binWidths)[1])
cumSumBinWidths <- cbind(v,cumSumBinWidths)
midPoints <- cumSumBinWidths + (0.5*cbind(binWidths,v))
midPoints <- t(head(t(midPoints),-1))
nClass <- dim(binWidths)[2]
return(list(midPoints=midPoints,cumSumBinWidths=cumSumBinWidths,nClass=nClass))
}
# plots the rates of survival and discrad for species where entries are no zero
plot_Discards_Survival <- function(outPath,inputData){
#singles <- read.table(paste0(inPath,"singletons.csv"),sep=",",row.names=1,header=FALSE)
#discards <- read.csv(paste0(inPath,"fishing_discards.csv"),sep=",",skip=3,header=TRUE,row.names=1)
#survival <- read.csv(paste0(inPath,"fishing_discardsSurvival.csv"),sep=",",skip=3,header=TRUE,row.names=1)
#singletons <- simplify2array(singles)
#names(singletons) <- row.names(singles)
nFleets <- inputData$Nfleets
nSpecies <- inputData$Nspecies
nClass <- inputData$Nsizebins
discards <- inputData$discardCoef
survival <- inputData$discardSurvival
for (isp in 1:nSpecies) {
# take in chunks of nFleets
istart <- (isp-1)*nFleets + 1
iend <- istart+nFleets-1
speciesDataD <- discards[istart:iend,]
if(sum(as.numeric(as.matrix(speciesDataD))) > 0) {
speciesGearName <- row.names(speciesDataD)[1]
numLetters <- nchar(speciesGearName)
speciesName <- substr(speciesGearName,1,numLetters-2)
gearType <- substr(speciesGearName,numLetters-1,numLetters)
matplot(t(speciesDataD),type="l")
}
}
}
plot_observed_Biomass_Catch <- function(outPath,inputData){
#biomass <- read.table(paste0(inPath,"observation_biomass.csv"),sep=",",header=TRUE)
biomass <- t(inputData$observedBiomass)
#catch <- read.table(paste0(inPath,"observation_catch.csv"),sep=",",header=TRUE)
catch <- t(inputData$observedCatch)
nCols <- dim(biomass)[2]
png(paste0(outPath,"/Hydra_Biomass.png"),width=11.5,height=8,units="in",res=300)
par(oma=c(1,4,1,2))
matplot(biomass[,"year"],biomass[,c(2:nCols)], type="l",xlab="Year",ylab = "",main="Observed Biomass")
legend("topright",legend=colnames(biomass)[2:nCols],col=seq_len(nCols-1),lty=seq_len(nCols-1))
graphics.off()
png(paste0(outPath,"/Hydra_Catch.png"),width=11.5,height=8,units="in",res=300)
par(oma=c(1,4,1,2))
matplot(catch[,"year"],catch[,c(2:nCols)], type="l",xlab="Year",ylab = "",main="Observed Catch")
legend("topright",legend=colnames(catch)[2:nCols],col=seq_len(nCols-1),lty=seq_len(nCols-1))
graphics.off()
# png(paste0(outPath,"/Hydra_biomass.png"),width=11.5,height=8,units="in",res=300)
}
# plots all covariate data
plot_CovariateData <- function(outPath,inputData) {
#growthCov <- read.table(paste0(inPath,"growth_covariates.csv"),sep=",",header=TRUE)
#maturityCov <- read.table(paste0(inPath,"maturity_covariates.csv"),sep=",",header=TRUE)
#recruitmentCov <- read.table(paste0(inPath,"recruitment_covariates.csv"),sep=",",header=TRUE)
growthEffects <- inputData$growthCovEffects #read.table(paste0(inPath,"growth_covariateEffects.csv"),sep=",",header=TRUE)
maturityEffects <- inputData$maturityCovEffects #read.table(paste0(inPath,"maturity_covariateEffects.csv"),sep=",",header=TRUE)
recruitmentEffects <- inputData$recruitCovEffects #read.table(paste0(inPath,"recruitment_covariateEffects.csv"),sep=",",header=TRUE)
nSpecies <- dim(growthEffects)[1]
effects <- cbind(growthEffects,maturityEffects,recruitmentEffects)
covariates <- t(rbind(inputData$growthCov,inputData$maturityCov,inputData$recruitmentCov))
names(covariates) <- c("temperature(Growth)","maturity","recruitment")
nCovs <- dim(covariates)[2]
png(paste0(outPath,"/Hydra_Covariates.png"),width=8,height=11.5,units="in",res=300)
mat <- matrix(c(1,2,3,4,5,6),3,2)
graphics::layout(mat,widths=c(6,2))
par(mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+0.5)
for (icov in 1:nCovs) {
plot(covariates[,icov],type="l",cex=1.6,lwd=2)
legend("topleft",legend=names(covariates)[icov],bty="n",cex=1.5)
}
for (icov in 1:nCovs) {
plot(effects[,icov],c(1:nSpecies),type="l",cex=1.6,lwd=2)
mtext("Species Effect",line=1.9,side=2,cex=1)
}
graphics.off()
}
# size preferenmce * weight ratio
plot_SizePreference_weightRatio <- function(outPath,inputData) {
#sizePref <- as.data.frame(t(read.table(paste0(inPath,"mortality_M2sizePreference.csv"),sep=",",row.names=1,header=TRUE)))
speciesNames <- colnames(inputData$M2sizePrefMu)
speciesNames[speciesNames=="goosefish"] <- "monkfish"
nSpecies <- length(speciesNames)
# weight ratio
out <- binWidths_Calc(inputData)
midPoints <- out$midPoints
nClass <- out$nClass
png(paste0(outPath,"/Hydra_SizePreference.png"),width=8,height= 11.5,units="in",res=300)
par(mfrow=c(nSpecies/2,2),mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+.5)
weightRatio <- seq(0.0001,10,.0001)
meanSizePreference <- vector(mode="numeric",length = nSpecies)
for (isp in 1:nSpecies) {
meanlog <- inputData$M2sizePrefMu[isp]
sdlog <- inputData$M2sizePrefSigma[isp]
pdfLogN <- (1/(weightRatio*sdlog*sqrt(2*pi))) *exp(-((log(weightRatio)-meanlog)^2)/(2*sdlog^2))
meanSizePreference[isp] <- weightRatio[which(max(pdfLogN)==pdfLogN)]
plot(log(weightRatio),pdfLogN,type="l",xlab="",ylab="",yaxt="n",xaxt="n")
axis(1, at=c(log(c(.0001,.001, 0.01, .1,1,10))), labels=c(.0001,.001,.01,.1,1,10))
legend("topleft",legend=speciesNames[isp])
}
mtext("Weight Ratio",line = 1.7,outer=T, side=1,cex = 1.5)
mtext("PDF",line = -.5,outer=T, side=2,cex = 1.5)
graphics.off()
# weight ratio is calculated for each pair of species. However some species dont interact .We only print ratios of predators
#foodweb <- as.matrix(read.table(paste0(inPath,"foodweb.csv"),sep=",",header=TRUE,row.names=1))
predators <- which(colSums(inputData$foodweb)>0)
nPredators <- length(predators)
# calculate weight at each size class for each species
#lengthWeight <- read.table(paste0(inPath,"lengthweight_species.csv"),sep=",",header=TRUE,row.names=1)
weight <- vector(mode = "list", length = nSpecies)
for (isp in 1:nSpecies) {
a <- inputData$lenwta[isp]
b <- inputData$lenwtb[isp]
weight[[isp]] <- a*midPoints[isp,]^b
}
# for each predator we find ratios for each prey and plot
for (isp in 1:nPredators) {
ipred <- predators[isp]
prey <- which(as.logical(inputData$foodweb[,ipred]))
preyNames <- speciesNames[prey]
nPrey <- length(prey)
png(paste0(outPath,"/Hydra_weightRatio_",speciesNames[ipred],".png"),width=8,height= 11.5,units="in",res=300)
par(mfrow=c(nSpecies/2,2),mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+.5)
for (iprey in 1:nPrey) {
preyW <- matrix(rep(weight[[prey[iprey]]],nClass),nrow=nClass)
predW <- matrix(rep(weight[[ipred]],nClass),ncol=nClass,byrow=TRUE)
wR <- (preyW/predW)
wR[lower.tri(wR,diag=TRUE)] <- NA
matplot(t(wR),type="b",ylim=c(0,max(wR,na.rm = TRUE)),xlab="",ylab="")
abline(meanSizePreference[isp],0,col="black")
legend("topleft",legend=preyNames[iprey])
#invisible(readline("Press [Enter] ..."))
}
mtext(paste0("Predator: ",speciesNames[ipred]),line = -.5,outer=T, side=3,cex = 1.5)
mtext("Predator Size Class",line = 1.7,outer=T, side=1,cex = 1.5)
graphics.off()
}
}
# plots Fecundity
plot_Maturity <- function(outPath,inputData) {
maturity_cov <- inputData$maturityCov #(lazy data)
maturity_effects <- inputData$maturityCovEffects #(lazy data)
speciesNames <- names(inputData$maturityNu)
speciesNames[speciesNames=="goosefish"] <- "monkfish"
nSpecies <- length(speciesNames)
out <- binWidths_Calc(inputData)
midPoints <- out$midPoints
nClass <- out$nClass
png(paste0(outPath,"/Hydra_Maturity.png"),width=8,height= 11.5,units="in",res=300)
par(mfrow=c(nSpecies/2,2),mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+.5)
for (isp in 1:nSpecies) {
nu <- inputData$maturityNu[isp]
omega <- inputData$maturityOmega[isp]
maturityCalc <- 1/(1+exp(-(nu + omega*midPoints[isp,])))
if (speciesNames[isp]=="spinydog") {maturityCalc=c(rep(0,nClass-1),1)} # hardcoded in ADMB
plot(midPoints[isp,],maturityCalc,type="b")
legend("topleft",legend=speciesNames[isp])
}
mtext("Length (cm)",line = 1.7,outer=T, side=1,cex = 1.5)
mtext("Proportion Mature",line = -.5,outer=T, side=2,cex = 1.5)
graphics.off()
}
# intake
plot_Intake <- function(outPath,inputData) {
speciesNames <- names(inputData$intakeAlpha)
nSpecies <- length(speciesNames)
nClasses <- dim(inputData$intakeStomach)[2]
png(paste0(outPath,"/Hydra_Intake.png"),width=8,height= 11.5,units="in",res=300)
par(mfrow=c(nSpecies/2,2),mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+.5)
temperature <- seq(from=7.97, to=9.97, by=.01)
for (isp in 1:nSpecies) {
alpha <- inputData$intakeAlpha[isp]
beta <- inputData$intakeBeta[isp]
StomachContent <- inputData$intakeStomach[isp,]
IntakeCalc <- 24*alpha*exp(beta*temperature)
minI <- min(IntakeCalc)*min(StomachContent)
maxI <- max(IntakeCalc)*max(StomachContent)
if (speciesNames[isp] == "goosefish") {speciesNames[isp] <- "monkfish"}
for (isize in 1:nClasses) {
if (isize == 1) {
# plot(temperature,IntakeCalc,type="l")
plot(temperature,IntakeCalc*as.numeric(StomachContent[isize]),type="l",ylim=c(minI,maxI))
} else {
lines(temperature,IntakeCalc*as.numeric(StomachContent[isize]))
}
}
legend("topleft",legend=speciesNames[isp])
}
mtext("Temperature ('C)",line = 1.7,outer=T, side=1,cex = 1.5)
mtext("Daily Intake (g)",line = -.5,outer=T, side=2,cex = 1.5)
graphics.off()
}
# plots Fecundity
plot_Fecundity <- function(outPath,inputData) {
#fecundity <- as.data.frame(t(read.table(paste0(inPath,"fecundity_species.csv"),sep=",",row.names=1,header=TRUE)))
fecundity_theta <- inputData$fecundityTheta # (lazy data)
out <- binWidths_Calc(inputData)
midPoints <- out$midPoints
speciesNames <- names(inputData$fecundityd)
#fecundity_theta <- as.numeric(fecundity_theta[,-1])
speciesNames[speciesNames=="goosefish"] <- "monkfish"
nSpecies <- length(speciesNames)
png(paste0(outPath,"/Hydra_Fecundity.png"),width=8,height= 11.5,units="in",res=300)
par(mfrow=c(nSpecies/2,2),mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+.5)
for (isp in 1:nSpecies) {
fecund <- inputData$fecundityd[isp]*as.numeric(fecundity_theta[isp,])*(midPoints[isp,]^inputData$fecundityh[isp])
plot(midPoints[isp,],fecund,main="",type = "b")
legend("topleft",legend=speciesNames[isp])
}
mtext("Length (cm)",line = 1.7,outer=T, side=1,cex = 1.5)
mtext("Fecundity",line = -.5,outer=T, side=2,cex = 1.5)
graphics.off()
}
# plots all paramaterized versions of SR. Shepher, BH, ricker, segmented etc
plot_StockRecruitment <- function(outPath,inputData) {
##########################################################################3
# plot stock recruitment
#SR <- read.table(paste0(inPath,"recruitment_species.csv"),sep=",",row.names=1,header=TRUE)
recTypes <- inputData$recType
speciesNames <- inputData$speciesList
speciesNames[speciesNames=="goosefish"] <- "monkfish"
nSpecies <- length(speciesNames)
SR_types <- c("EggRicker","DS","Gamma","Ricker","BH","Shepherd","Hockey","Segmented")
png(paste0(outPath,"/Hydra_StockRecruitment.png"),width=8,height= 11.5,units="in",res=300)
par(mfrow=c(nSpecies/2,2),mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+.5)
attach(inputData)
for (isp in 1:nSpecies) {
# currentType <- SR_types[SR["type",isp]]
currentType <- SR_types[recTypes[isp]]
assign("tempAlpha",paste0("alpha",currentType))
alpha <- eval(as.name(tempAlpha))[isp]
assign("tempBeta",paste0("beta",currentType))
beta <- eval(as.name(tempBeta))[isp]
assign("tempShape",paste0("shape",currentType))
shape <- eval(as.name(tempShape))[isp]
if (recTypes[isp] == 1) { # Egg production : ricker
SSB <- seq(0,100000,5000)
R <- alpha*SSB*exp(-beta*SSB)
} else if (recTypes[isp] == 2) {# Deriso-Schnute
SSB <- seq(0,100000,5000)
R <- alpha*SSB*(1-beta*shape*SSB)^(1/shape)
} else if (recTypes[isp] == 3){ # Gamma
SSB <- seq(0,100000,5000)
R <- alpha*(SSB^shape)*exp(-beta*SSB)
} else if (recTypes[isp] == 4) { # Ricker
SSB <- seq(0,100000,5000)
R <- alpha*SSB*exp(-beta*SSB)
} else if (recTypes[isp] == 5) { # Beverton Holt
SSB <- seq(0,100000,5000)
R <- alpha*SSB/(1+ beta*SSB)
} else if (recTypes[isp] == 6) { # Shepherd
SSB <- seq(0,100000,5000)
R <- alpha*SSB/(1+ (SSB/beta)^shape)
} else if (recTypes[isp] == 7) { # Hockey Stick
breakpoint <- shape
SSB <- seq(0,2.5*breakpoint,breakpoint)
ind <- SSB > breakpoint
R <- alpha*SSB + (beta*(SSB-breakpoint))*ind
} else if (recTypes[isp] == 8) { # Segmented
breakpoint <- shape
SSB <- seq(0,2.5*breakpoint,breakpoint)
ind <- SSB > breakpoint
R <- alpha*SSB + (beta*(SSB-breakpoint))*ind
}
plot(SSB,R,main="",type = "l",ylim=c(0,max(R)))
legend("topleft",legend=speciesNames[isp],bty="n",cex=1.5)
legend("bottomright",legend=currentType,bty="n",cex=1)
}
detach(inputData)
mtext("SSB",line = 1.7,outer=T, side=1,cex = 1.5)
mtext("Recruitment",line = -.5,outer=T, side=2,cex = 1.5)
graphics.off()
}
# plots observed temperature
plot_Temperature <- function(outPath,inputData) {
temperature <- inputData$observedTemperature #read.table(paste0(inPath,"observation_temperature.csv"),sep=",",header=TRUE)
png(paste0(outPath,"/Hydra_Temperature.png"),width=11.5,height=8,units="in",res=300)
par(oma=c(1,4,1,2))
plot(temperature["year",],temperature["temperature",], type="b",xlab="Year",ylab = "Temperature ('C)")
graphics.off()
}
# plots fooweb structure in matrix and graph format
plot_Foodweb <- function(outPath,inputData) {
#foodweb <- as.matrix(read.table(paste0(inPath,"foodweb.csv"),sep=",",header=TRUE,row.names=1))
speciesnames <- inputData$speciesList #read.table(paste0(inPath,"speciesList.csv"),sep=",",header=TRUE)
guilds <- unique(as.data.frame(cbind(inputData$guildNames,inputData$guildMembership))) # (lazy data)
names(guilds) <- c("guildNames","guildMembership")
guilds$guildMembership <- as.integer(guilds$guildMembership)
guilds <- guilds[with(guilds,order(guildMembership)),]
nGuilds <- dim(guilds)[1]
colorVec <- viridis::viridis(nGuilds)
foodweb <- inputData$foodweb
web <- igraph::graph_from_adjacency_matrix(foodweb,mode="directed") #(lazy data)
# average-linkage clustering method
#jaccard <- similarity(web,vids=V(web),mode="all",loops = TRUE,method="jaccard")
#cc = hclust(1-as.dist(jaccard), method = "average")
#plot(cc)
#return(jaccard)
# make igraph object
# add attributed for plottinf
#webAttributes <- data.frame(guild=NULL,guildMember=NULL,color=NULL)
nSpecies <- dim(foodweb)[1]
for (isp in 1:nSpecies) {
#ind <- which(names(igraph::V(web)[isp]) == speciesList$species)
ind <- which(names(igraph::V(web)[isp]) == speciesnames)
igraph::V(web)$guildmember[isp] <- inputData$guildMembership[ind]
igraph::V(web)$guild[isp] <- as.character(inputData$guildNames[ind])
igraph::V(web)$color[isp] <- colorVec[inputData$guildMembership[isp]]
igraph::V(web)$predPrey[isp] <- inputData$predOrPrey[isp]
}
l <- igraph::layout_with_lgl(web)
tl <- NetIndices::TrophInd(foodweb) # determins tropic level. netIndices package
l[,2] <- tl$TL # replace the y location column with tropich level
png(paste0(outPath,"/Hydra_foodWeb_matrix.png"),width=4,height=4,units="in",res=300)
par(oma=c(1,0,1,1))
bipartite::visweb(web = foodweb, prednames = TRUE, preynames = TRUE, labsize = .5,
type = "none", clear = FALSE)
graphics.off()
png(paste0(outPath,"/Hydra_foodWeb_network.png"),width=6,height=6,units="in",res=300)
par(oma=c(1,1,1,1))
igraph::plot.igraph(web,edge.arrow.size=.5,vertex.color=igraph::V(web)$color,vertex.size= 20* 2^(igraph::V(web)$predPrey),
layout=l)
legend("topright",legend=guilds$guildNames,pch=21,
pt.bg=colorVec, pt.cex=2, cex=.8, bty="n", ncol=1)
graphics.off()
}
### plots Y1N and corresponding starting biomass
### also plots length weight relationsships and bin widths
# all connected
plot_Y1N_Biomass_lengthweight_bins <- function(outPath,inputData) {
cex <- 0.8
#Y1N <- read.table(paste0(inPath,"observation_Y1N.csv"),sep=",",header=TRUE,row.names=1)
out <- binWidths_Calc(inputData)
midPoints <- out$midPoints
cumSumBinWidths <- out$cumSumBinWidths
#singletons <- read.table(paste0(inPath,"singletons.csv"),sep=",",header=TRUE,row.names=1)
weightConversion <- inputData$wtconv # from grams
#lengthWeight <- read.table(paste0(inPath,"lengthweight_species.csv"),sep=",",header=TRUE,row.names=1)
nSpecies <- length(row.names(inputData$Y1N))
speciesNames <- row.names(inputData$Y1N)
nBins <- dim(inputData$Y1N)[2]
## binwidths
#png(system.file("rmd","Hydra_binWidths.png",package="hydradata"),width=11.5,height=8,units="in",res=300)
png(paste0(outPath,"/Hydra_binWidths.png"),width=11.5,height=8,units="in",res=300)
par(oma=c(1,4,1,2))
for (isp in 1:nSpecies) {
if (isp == 1) {
plot(cumSumBinWidths[isp,],rep(isp,nBins+1),type="b",xlim=c(min(cumSumBinWidths),max(cumSumBinWidths))
,ylim=c(0,nSpecies+1),xlab="length (cm)",ylab=NA,yaxt="n",xaxt="n")
lines(cumSumBinWidths[isp,],rep(isp,nBins+1))
} else {
lines(cumSumBinWidths[isp,],rep(isp,nBins+1))
points(cumSumBinWidths[isp,],rep(isp,nBins+1))
}
}
axis(side=2, at=c(1:10),labels=speciesNames,las=1)
axis(side=1, at=seq(0,140,20),labels=seq(0,140,20))
graphics.off()
## length-weight relationship
png(paste0(outPath,"/Hydra_lengthWeight.png"),width=8,height=11.5,units="in",res=300)
par(mfrow=c(nSpecies/2,2),mar=c(2,4,0,4)+.0,oma=c(3,1,1,1)+0.0)
for (isp in 1:nSpecies) {
lengthOfFish <- seq(0,cumSumBinWidths[isp,nBins+1],1)
a <- inputData$lenwta[isp]
b <- inputData$lenwtb[isp]
plot(lengthOfFish,a*lengthOfFish^b,type="l",ylab="",xlab="")
if ((isp %% 2) == 1) {mtext(side=2,line=2.5,"weight (g)",cex=cex)}
if ((isp==nSpecies) | (isp ==(nSpecies-1))) {mtext(side=1,line=2.5,cex=cex,"length (cm)")}
legend("topleft",legend=speciesNames[isp])
}
graphics.off()
## Y1N and biomass
png(paste0(outPath,"/Hydra_Y1N_Biomass.png"),width=8,height=11.5,units="in",res=300)
par(mfrow=c(nSpecies/2,2),mar=c(2,4,0,4)+.0,oma=c(3,1,1,1)+0.0)
for (isp in 1:nSpecies) {
a <- inputData$lenwta[isp] # (lazy data)
b <- inputData$lenwtb[isp] # (lazy data)
plot(1:nBins,inputData$Y1N[isp,],type="b",ylab="",xlab="")
if ((isp %% 2) == 1) {mtext(side=2,line=2,"# individuals (mil)",cex=cex)}
if ((isp==nSpecies) | (isp ==(nSpecies-1))) {mtext(side=1,line=2,cex=cex,"size class")}
par(new=T)
# biomass part
weightG <- a*midPoints[isp,]^b
plot(1:nBins,weightConversion*inputData$Y1N[isp,]*weightG,axes=F,xlab=NA,ylab=NA,type="l",lty=2)
axis(side=4)
if ((isp %% 2) == 0) {mtext(side=4,line=2.2,cex=cex,"Biomass (g)")}
legend("topleft",legend=c(speciesNames[isp], "Y1N","Biomass"),
lty=c(0,1,2),pch=c(NA,NA,NA), col=c("black", "black","black"))
}
graphics.off()
}
# fishing selectivities
plot_Selectivities <- function(outPath,inputData) {
# selectivity
# requires selectivity_c, selectivity_d and size bins
png(paste0(outPath,"/Hydra_Selectivities.png"),width=8,height=11.5,units="in",res=300)
sel_c <- inputData$fisherySelectivityc # lazy data # read.table(paste0(inPath,"fishing_selectivityc.csv"),sep=",",header=TRUE,row.names=1)
sel_d <- inputData$fisherySelectivityd # lazy data # read.table(paste0(inPath,"fishing_selectivityd.csv"),sep=",",header=TRUE,row.names=1)
bins <- inputData$binwidth # lazy data #read.table(paste0(inPath,"length_sizebins.csv"),sep=",",header=TRUE)
bins <- bins[,1:ncol(bins)-1]
nFleets <- ncol(sel_c)
nSpecies <- nrow(bins)
nClasses <- ncol(bins)
selectivity <- matrix(list(data=NA),nSpecies,nFleets)
selectivityMid <- matrix(list(data=NA),nSpecies,nFleets)
par(mfrow=c(nSpecies/2,2),mar=c(0,2,1,1)+.75,oma=c(3,1,1,1)+0.5)
#colVec <- c("black","green","red","blue")
colVec <- c("black","blue","red")
lineTypeVec <- c(1,2,3,4)
for (isp in 1:nSpecies) {
binsInterval <- unlist(bins[isp,])
cum <- c(0,cumsum(binsInterval))
midBins <- cum[1:nClasses]+(diff(cum)/2)
for (jfleet in 1:(nFleets-1)) {
x <- seq(0,tail(cum,1),1)
#print(1/(1 + exp(-sel_c[isp,jfleet]-sel_d[isp,jfleet]*midBins)))
selectivity[[isp,jfleet]] <- 1/(1 + exp(-sel_c[isp,jfleet]-sel_d[isp,jfleet]*x))
selectivityMid[[isp,jfleet]] <- 1/(1 + exp(-sel_c[isp,jfleet]-sel_d[isp,jfleet]*midBins))
if (jfleet == 1) {
plot(x,selectivity[[isp,jfleet]],type="l",lwd=2,ylab="selectivity",xlab="length(cm)",xlim=c(0,max(midBins)),ylim=c(0,1)) #main=rownames(sel_c)[isp]
#ext(midBins[1],.95,rownames(sel_c)[isp],pos=4,cex=1.5)
if (rownames(sel_c)[isp]=="goosefish") {legend("topleft",legend="monkfish",bty="n",cex=1.5)} else {
legend("topleft",legend=rownames(sel_c)[isp],bty="n",cex=1.5)
}
} else {
lines(x,selectivity[[isp,jfleet]],col=colVec[jfleet],lty=lineTypeVec[jfleet],lwd=2)
points(midBins,selectivityMid[[isp,jfleet]],col=colVec[jfleet],pch=19)
}
}
}
mtext("Length (cm)",line = 1.9,outer=T, side=1,cex = 1.5)
mtext("Selectivity",line = -.3,outer=T, side=2,cex = 1.5)
# add text as a legend to plot in separate panel
#plot(1, type="n", axes=FALSE, xlab="", ylab="")
#legend(.5,1.5,legend=colnames(sel_c)[1:(nFleets-1)],fill=colVec,cex=1.8)
graphics.off()
}
# effort by fleet
plot_FishingEffort <- function(outPath,inputData) {
#plot effort
obs_effort <- t(inputData$observedEffort) #lazy data # read.table(paste0(inPath,"observation_effort.csv"),sep=",",header=TRUE)
nFleets <- dim(obs_effort)[2]-1
fleetNames <- colnames(obs_effort)
png(paste0(outPath,"/Hydra_EffortAnomoly.png"),width=8,height=11.5,units="in",res=300)
par(mfrow=c(nFleets,1),mar=c(0,2,1,1)+0.75,oma=c(2,1,1,1)+0.5)
for (ifleet in 1:nFleets) {
#plot(obs_effort[,1],obs_effort[,ifleet+1],type="l",main=colnames(obs_effort)[ifleet+1],xlab="",ylab="Effort")
anomoly <- (obs_effort[,ifleet+1] - mean(obs_effort[,ifleet+1]))/sd(obs_effort[,ifleet+1])
if (any(is.nan(anomoly))) {
# fllet not fully parameterized
plot(0,0)
next
}
plot(obs_effort[,1],anomoly,type="l", cex=1.6,lwd=2)
abline(h=0)
if (fleetNames[ifleet+1] =="benthic") {fleetNames[ifleet+1] <- "Demersal Trawl"}
if (fleetNames[ifleet+1] =="pelagic") {fleetNames[ifleet+1] <- "Pelagic Trawl"}
if (fleetNames[ifleet+1] =="longline") {fleetNames[ifleet+1] <- "Fixed Gear"}
legend("topleft",legend=fleetNames[ifleet+1],bty="n",cex=1.5)
}
mtext("Anomoly",line=-.6,outer=T,side=2,cex=1.5)
graphics.off()
png(paste0(outPath,"/Hydra_Effort.png"),width=8,height=11.5,units="in",res=300)
par(mfrow=c(nFleets,1),mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+0.5)
for (ifleet in 1:nFleets) {
plot(obs_effort[,1],obs_effort[,ifleet+1],type="l",cex=1.6,lwd=2)
legend("topleft",legend=fleetNames[ifleet+1],bty="n",cex=1.5)
}
mtext("Effort (Hours Fished)",line=-.6,outer=T,side=2,cex=1.5)
graphics.off()
############################################################################
# plot fishing _q's
#fishing_q <- read.table(paste0(inPath,"fishing_q.csv"),sep=",",row.names=1,header=TRUE)
fishing_q <- inputData$fisheryq # lazy data
fleetNames <- colnames(fishing_q)
nSpecies <- dim(fishing_q)[1]
png(paste0(outPath,"/Hydra_FishingQ.png"),width=8,height=11.5,units="in",res=300)
#png("HydraFishingQ.png",width=11.5,height=8,units="in",res=300)
par(mfrow=c(nFleets,1),mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+0.5)
#par(mfrow=c(1,nFleets),mar=c(0,2,5,1)+0.75,oma=c(3,1,5,1)+0.5)
for (ifleet in 1:nFleets) {
#plot(fishing_q[,ifleet],ylab="q",xaxt="n",main=colnames(fishing_q)[ifleet])
barplot(fishing_q[,ifleet],ylab="q",names.arg = rownames(fishing_q),cex.names = 1.1)
# legend("topleft",legend=colnames(fishing_q)[ifleet],bty="n",cex=1.7)
if (fleetNames[ifleet] =="benthic") {fleetNames[ifleet] <- "Demersal Trawl"}
if (fleetNames[ifleet] =="pelagic") {fleetNames[ifleet] <- "Pelagic Trawl"}
if (fleetNames[ifleet] =="longline") {fleetNames[ifleet] <- "Fixed Gear"}
legend("topleft",legend=fleetNames[ifleet],bty="n",cex=1.7)
#axis(side=1, at=c(1:10), labels=rownames(fishing_q))
}
#mtext("Fishery Q",line=-.6,outer=T,side=2,cex=1.5)
graphics.off()
###########################################################################
# plot exploitation rates as a histogram for each fleet
for (ifleet in 1:nFleets) {
png(paste0(paste0(outPath,"/Hydra_ExploitationRate_"),names(fishing_q)[ifleet],".png"),width=8,height=11.5,units="in",res=300)
par(mfrow=c(nSpecies/2,2),mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+0.5)
effort <- obs_effort[,ifleet+1]
for (isp in 1:nSpecies) {
hist(fishing_q[isp,ifleet]*effort,main="",breaks=10,xlab="Exploitation")
if (rownames(fishing_q)[isp]=="goosefish") {legend("topleft",legend="monkfish",bty="n",cex=1.5)} else {
legend("topright",legend=rownames(fishing_q)[isp],bty="n",cex=1.5)
}
}
mtext("Exploitation",line = 1.7,outer=T, side=1,cex = 1.5)
mtext("Frequency (of yrs)",line = -.5,outer=T, side=2,cex = 1.5)
title(names(fishing_q)[ifleet],outer=TRUE,line=-1,cex.main=2)
graphics.off()
}
}
# species growth
plot_Growth <- function(outPath,inputData) {
##########################################################################3
# plot growth curves
#growth <- read.table(paste0(inPath,"growth_species.csv"),sep=",",row.names=1,header=TRUE)
growthBins <- inputData$binwidth # lazy data read.table(paste0(inPath,"length_sizebins.csv"),sep=",",header=TRUE)
speciesNames <- row.names(growthBins)
nSpecies <- length(speciesNames)
speciesNames[speciesNames=="goosefish"] <- "monkfish"
png(paste0(outPath,"/Hydra_Growth.png"),width=8,height=11.5,units="in",res=300)
par(mfrow=c(nSpecies/2,2),mar=c(0,2,1,1)+0.75,oma=c(3,1,1,1)+.5)
growthFunctions <- c("na","exponential","na","von Bertalanffy")
maxLinf <- max(inputData$growthLinf)
t <- seq(0,40,1)
for (isp in 1:nSpecies) {
growthThroughInterval1 <- growthBins[isp,1]
if(inputData$growthType[isp] == 2) {# exponential
psi <- inputData$growthPsi[isp]
kappa_r <- inputData$growthKappa[isp]
growthFunc <- psi*t^kappa_r
delta_t <- (growthThroughInterval1/psi)^(1/kappa_r)
func <- growthFunctions[2]
} else if (inputData$growthType[isp] == 4) {# von Bert
k <- inputData$growthK[isp]
linf <- inputData$growthLinf[isp]
delta_t <- (1/k)*log(linf/(linf-growthThroughInterval1))
growthFunc <- linf*(1-exp(-k*t))
func <- growthFunctions[4]
}
plot(t,growthFunc,type="l",xlab="age",ylab="length",ylim=c(0,maxLinf))
lines(rep(delta_t,2),c(0,linf),col="red",lty=2)
legend("topleft",legend=speciesNames[isp],bty="n",cex=1.5)
legend("bottomright",legend=func,bty="n",cex=1)
}
mtext("Age (yrs)",line = 1.7,outer=T, side=1,cex = 1.5)
mtext("Length (cm)",line = -.5,outer=T, side=2,cex = 1.5)
graphics.off()
}
# step/ramp for assessment module
plot_Assessment_Step <- function(outPath,inputData) {
# plot assessment thresholds - step function
png(paste0(outPath,"/Hydra_Assessment_Thresholds.png"),width=8,height= 11.5,units="in",res=300)
#thresholds <- read.table(paste0(inPath,"assessmentThresholds.csv"),sep=",",row.names=1,header=TRUE)
nExperiments <- dim(inputData$exploitationOptions)[2] # lazy data
nSteps <- dim(inputData$thresholds)[1]
par(mfrow=c(nExperiments/2,2),mar=c(4,4,2,1)+.5,oma=c(1,1,1,1)+.0)
actionLevels <- inputData$thresholds #as.numeric(rownames(thresholds))
actionLevels <- head(actionLevels,-1) # removes last value
for (iExp in 1:nExperiments) {
stepObj <- stepfun(actionLevels,inputData$exploitationOptions[,iExp])
plot.stepfun(stepObj,ylim=c(0,max(inputData$exploitationOptions)),main="" ,ylab = "Exploitation Rate",xlab="B/B0")
abline(inputData$exploitationOptions[2,iExp],0,col="red")
legend("topleft",legend=paste0("Max rate = ",max(inputData$exploitationOptions[,iExp])))
}
graphics.off()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.