R/plot_gear_activity_distribution.R

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

#
# plot_gear_activity_distribution.R
#
#' Plot the distribution of fishing activity rate
#'
#' Plots a colour-scaled matrix of the distribiution of fishing gear activity density across seabed habitats.
#'
#' The function plots a matrix of seabed habitats vs fishing gears with each cell colour-shaded to indicate activity density (white = 0, purple = high).
#'
#' The spatial distribution of fishing gear activity in the model is defined by two input data sets:
#' 
#' \itemize{
#'   \item Vector of whole domain activity density of each fishing gear (s/d per m2 of whole model domain)
#'   \item Matrix of the proportional distribution of whole domain activity density of each gear (rows) across seabed habitats (columns)
#' }
#'
#' The activity density of each gear in each habitat is then obtained by multiplying the whole domain activity density into the proportional
#' distribution matrix, and dividing by the area-proportions of habitats in the domain. The units of habitat-specific activity density
#' are then (s/d/m2).
#'
#' The vector of area-proportions of each habitat in the model domain is part of the moldel configuration parameter set.
#' 
#' The calculation take account of any activity multipliers set in the csv inputs to be applied to the activity density.
#'
#' When interpreting patterns of activity density, bear in mind that the impact (harvest ratio, seabed abrasion) per unit activity can vary considerably between gears.
#'
#' @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_activity_distribution <- function(model) {

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

data                <- elt(model, "data")
physical.parameters <- elt(data,"physical.parameters")
habareas       <- elt(physical.parameters,"habitat_areas")

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")
FGSdata <- elt(fleet.model,"gear_habitat_activity")
gear_mult <- elt(fleet.model,"gear_mult")

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

#Get the proportional spatial distribution data
PSPACE<-(as.matrix(FGSdata))
rownames(PSPACE)<-gear_labels
colnames(PSPACE)<-c("S0","S1","S2","S3","D0","D1","D2","D3")

habitat_labels_l <-c("Inshore rock ",
                     "Inshore fine ",
                     "Inshore medium ",
                     "Inshore coarse ",
                     "Offshore rock ",
                     "Offshore fine ",
                     "Offshore medium ",
                     "Offshore coarse ")

#Make a matrix of the activity density by habitat
ACTSPACE<-PSPACE
ACTSPACE[,]<-NA

for(j in 1:ncol(PSPACE)){
if(habareas[j]>0) ACTSPACE[,j]<-PSPACE[,j]*ACTdens/habareas[j]
qq<-which(ACTSPACE[,j]==0)
ACTSPACE[qq,j]<-NA
}

gear_labels_l <- gear_labels[12:1]

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

plotmat<-ACTSPACE
#Rows are gears columns are gears - 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,8,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:8, par("usr")[1] - 0, srt = 45, adj = 1,labels = habitat_labels_l, xpd = TRUE,cex=1)
    mtext(bquote("Activity density (log"[e] ~ "s" ~ ".m"^-2 ~ ".d"^-1 ~ ")"),cex=1.2,side=3,line=0.5)
    abline(v=seq(0.5,7.5,by=1))
    abline(h=seq(1.5,12.5,by=1))
    legend(grconvertX(0.38, "ndc", "user"), grconvertY(0.09, "ndc", "user"),
    c( "no activity",
      ((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,habitats=habitat_labels_l))

return(retdata)

}