R/plot_gear_guild_harvest_ratios.R

In StrathE2E2: End-to-End Marine Food Web Model

#
# plot_gear_guild_harvest_ratios.R
#
#' Plot the distribution of harvest ratios across guilds.
#'
#' Plots a colour-scaled matrix of the distribiution of harvest ratios due to each gear across across the range of guilds in the model.
#'
#' The function plots a matrix of the values of harvest ratio (proportion of exploitable biomass of each guild capatured per day) for each guild in the ecology model (columns) arising from the activity of each gear (rows). 
#' Hence each row represents the selectivity pattern of each gear with respect to guilds.
#' Cells in the matrix are colour-shaded to indicate log(e) transformed harvest ratio (colour gradient: white = 0, purple = high).
#'
#' Harvest ratio is calculated by the fishing fleet model and piped into the ecology model. The fleet model takes inputs of activity density by each gear,
#' and a matrix of fishing power, to calculate effort (effort = activity * power). Harvest ratio is assumed to be lineally related to effort by a scaling
#' coefficient which is input as a parameter to the fleet model.
#' 
#' Note that the visualization of harvest ratio generated by this function is based on the csv input data to the fleet model (including any multipliers to be applied to
#' either fishing activity or harvest ratio scaling parameters). However, it does not reflect any effects on harvest ratio arising from discarding scenarios since
#' these are dynamic and generated at run-time within the ecology model.
#'
#' @param model R-list object defining the model configuration compiled by the e2e_read() function
#'
#' @return Graphical display in a new graphics window and a list object comprising the plotted data and axis labels.
#'
#' @importFrom graphics grconvertX grconvertY image
#' @importFrom grDevices colorRampPalette
#'
#' @noRd
#
# ------------------------------------------------------------------------------

plot_gear_guild_harvest_ratios <- function(model) {

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

data                <- elt(model, "data")
fleet.model<-elt(data, "fleet.model")

gear_labels <- elt(fleet.model,"gear_labels")
gear_codes <- elt(fleet.model,"gear_codes")
ACTdens <- elt(fleet.model,"gear_activity")
FGSdata <- elt(fleet.model,"gear_habitat_activity")
gear_mult <- elt(fleet.model,"gear_mult")
HRscale_vector <- elt(fleet.model,"HRscale_vector")
HRscale_vector_multiplier <- elt(fleet.model,"HRscale_vector_multiplier")
power  <-  elt(fleet.model,"gear_group_rel_power")

ACTdens<-ACTdens*gear_mult   # Usually gear_mult will be 1.0 for all gears but could be otherwise
HRscale_vector<-HRscale_vector*HRscale_vector_multiplier # Usually gear_multHRscale_vector_multiplier will be 1.0 for all groups but could be otherwise

effort<-power
effort[]<-0
#Multiply the activity density into the relative power matrix to get effort
for(jj in 1:(ncol(effort))){
effort[,jj] <- power[,jj] * ACTdens
}

#Multiply the HRscale vector into the effort matrix to get HR
HRmatrix<-effort
HRmatrix[]<-0
for(jj in 1:(nrow(HRmatrix))){
HRmatrix[jj,]<-effort[jj,] * HRscale_vector
}

PSPACE<-(as.matrix(HRmatrix))
rownames(PSPACE)<-gear_labels
colnames(PSPACE)<-colnames(power)

group_labels_l <-c("Planktiv.fish ",
                     "Demersal fish ",
                     "Migratory fish ",
                     "Susp/dep.benthos ",
                     "Carn/scav.benthos ",
                     "Carniv.zoo ",
                     "Birds ",
                     "Pinnipeds ",
		     "Cetaceans ",
		     "Macrophytes ")

for(j in 1:ncol(PSPACE)){
for(k in 1:nrow(PSPACE)){
if(PSPACE[k,j]==0) PSPACE[k,j]<-NA
}
}

gear_labels_l <- gear_labels[12:1]

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

plotmat<-PSPACE
#Rows are gears columns are groups - gear1 group 1 in top left
#So we need to invert the matrix top to bottom
plotmat_i<-(plotmat)
plotmat_i[,]<-0
for(i in 1:nrow(plotmat_i)){
plotmat_i[i,]<-plotmat[(nrow(plotmat_i)-(i-1)),]
}
rownames(plotmat_i)<-rownames(plotmat)[seq(12,1,by=-1)]

#Now we need to transpose the plotmat_i matrix to get bottom-left to top-right
plotmat_it<-t(plotmat_i)

plotmat_it<-log(plotmat_it) # log transform the data

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

par(mfrow=c(1,1))
par(mar=c(8,16,2,1))
    nsamples <- 20
    gridmax<-max(plotmat_it,na.rm=TRUE)
    gridmin<-min(plotmat_it,na.rm=TRUE)
    colMap <- colorRampPalette(c("white","purple" ))(nsamples)    #normal colour scale
    image(seq(1,10,by=1),seq(1,12,by=1),(plotmat_it[,1:12]),col=colMap,xaxt="n",yaxt="n",ann=FALSE)
    axis(1, labels = FALSE,tck=0)
    axis(side=2,at=seq(1,12),labels=gear_labels_l,las=1,cex.axis=0.9)
    text(1:10, par("usr")[1] - 0, srt = 45, adj = 1,labels = group_labels_l, xpd = TRUE,cex=1)
    mtext(bquote("Harvest ratio (log"[e] ~ "d"^-1 ~ ")"),cex=1.2,side=3,line=0.5)
    abline(v=seq(0.5,9.5,by=1))
    abline(h=seq(1.5,12.5,by=1))
    legend(grconvertX(0.45, "ndc", "user"), grconvertY(0.09, "ndc", "user"),
    c( "no catch",
      ((floor((gridmin+((gridmax-gridmin)/2))*100))/100),
      ((floor((gridmax)*100))/100) ),
    fill = colMap[c(1, (nsamples/2), nsamples)], ncol=3, bty="n", xpd = NA)

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

retdata<-list(data=plotmat_it,
	      labels=list(gears=gear_labels_l,guilds=group_labels_l))

return(retdata)

}