R/ccir_timeseries_exploitation_plots.r

In AMCOOK/bio.ccir:

#' @export

ccir_timeseries_exploitation_plots <- function(x, xrange = NULL, running.median =T , fdir=file.path(project.figuredirectory(bio.project),'ccir'), running.length = 3,fname = NULL,sex.ratio = NULL,bio.project ='bio.lobster',by.sex=F,combined.LFA=F,landings=NULL) {
	
  #x is output from ccir_collapse_summary 
	# both sexes same plot
	lff = x$LFA[1]
	if(is.null(fname)) fname = paste('TS.exploitation',lff,x$Grid[1],'png',sep=".")	
	if(is.null(xrange)) xr = range(x$Yr)
	if(is.null(sex.ratio)) sr = 0.5
	yr = c(0,1)
	fdir = fdir
	if(!dir.exists(fdir)) dir.create(fdir,recursive=T, showWarnings=F)

if(!combined.LFA){

  if(grepl(1,x$Grid[1])) Main = 'LFA 35'
  if(grepl(81,x$Grid[1])) Main = 'LFA 34'
  if(grepl(351,x$Grid[1])) Main = 'LFA 27 South'
	if(grepl(356,x$Grid[1])) Main = 'LFA 27 North'
	if(grepl(341,x$Grid[1])) Main = 'LFA 29'
	if(grepl(345,x$Grid[1])) Main = 'LFA 30'
	if(grepl(337,x$Grid[1])) Main = 'LFA 31A'
	if(grepl(332,x$Grid[1])) Main = 'LFA 31B'
	if(grepl(323,x$Grid[1])) Main = 'LFA 32'
	if(grepl(301,x$Grid[1])) Main = 'LFA 33 West'
	if(grepl(313,x$Grid[1])) Main = 'LFA 33 East'
	}

if(combined.LFA){
	a = split(x,f=x$Grid)
	aa = names(a)
	h = grep(names(landings[2]),aa)
	if(grepl(351,aa[h])) Main = 'LFA 27 Combined'
	if(grepl(301,aa[h])) Main = 'LFA 33 Combined'
	
	d = merge(a[[h]],landings[,1:2])
	names(d)[ncol(d)] <- 'Landings'
	h = grep(names(landings[3]),aa)
	e = merge(a[[h]],landings[,c(1,3)])
	names(e)[ncol(e)] <- 'Landings'
	e = rbind(d,e)
	ll = aggregate(Landings~Yr,data=e,FUN=sum)
	names(ll)[2] <- 'TotLand'
	e = merge(e,ll)
	r = cbind(Yr = e$Yr,e$Landings * e[,c('ERfl','ERfm','ERfu','ERf75')],TotLand = e$TotLand)
	xx = aggregate(cbind(ERfl,ERfm,ERfu,ERf75)~Yr+TotLand,data=r,FUN=sum)
	x = cbind(Yr=xx$Yr,xx[,c('ERfl','ERfm','ERfu','ERf75')] / xx$TotLand)
	x = x[order(x$Yr),]
	fname = paste('TS.exploitation',lff,'combined','png',sep=".")	
}

	png(filename=file.path(fdir,fname),width=8, height=5.5, units = "in", res = 800)	
	
	xlab = 'Year'
	ylab = 'Exploitation Index'
if(by.sex) {
with(subset(x,Sex==1),{
	Yr = Yr-0.1
	plot(Yr,ERfm,xlab=xlab,ylab=ylab,xlim=xr,ylim=yr,type='p',pch=16,col=1)
	arrows(x0 = Yr,y0 = ERfm,y1 = ERfu,length=0,lwd=1,col=1)
	arrows(x0 = Yr,y0 = ERfm,y1 = ERfl,length=0,lwd=1,col=1)
  		})

with(subset(x,Sex==2),{
	Yr = Yr+0.1
	points(Yr,ERfm,pch=18,col=2)
	arrows(x0 = Yr,y0 = ERfm,y1 = ERfu,length=0,lwd=1,col=2,lty=2)
	arrows(x0 = Yr,y0 = ERfm,y1 = ERfl,length=0,lwd=1,col=2,lty=2)
		})
x$Wts = x$NRef + x$NExp
  b = as.data.frame(sapply(split(x,x$Yr),function(x) weighted.mean(x$ERfm,w=x$Wts)))
  names(b) <- 'WtMean'
  b$Yr = rownames(b)
  with(b,{
      	rmean = runmed(WtMean,k=running.length,endrule='median')
        rmean.yr = Yr[1:length(Yr)]
		  lines(rmean.yr,rmean,lty=1,lwd=2,col='blue')
		  })
	legend('bottomleft',lty=c(0,0,1),pch=c(16,16,0),legend=c('Male','Female','Weighted Running Median'),col=c('black','red','blue'))
	} else {

		with(x,{
	Yr = Yr
	plot(Yr,ERfm,xlab=xlab,ylab=ylab,xlim=xr,ylim=yr,type='p',pch=16,col=1)
	arrows(x0 = Yr,y0 = ERfm,y1 = ERfu,length=0,lwd=1,col=1)
	arrows(x0 = Yr,y0 = ERfm,y1 = ERfl,length=0,lwd=1,col=1)
      	rmean = runmed(ERfm,k=running.length,endrule='median')
        rmean.yr = Yr[1:length(Yr)]
		  lines(rmean.yr,rmean,lty=1,lwd=2,col='blue')
		  })
	}
	title(Main)
	#browser()
	#savePlot(file.path(fdir,fname),type='png')
	dev.off()
	x$LFA = Main
	return(x)

}