inst/ExploitationEstimates.r

	options(stringsAsFactors=F)
	load_all('~/git/bio.eels')
	dadir = file.path(project.datadirectory('bio.eels'),'data')
	fidir = file.path(project.datadirectory('bio.eels'),'figures')

	fis = read.csv(file.path(dadir,'ElverFishery_revised_for_adam.csv'))
	ar = read.csv(file.path(dadir,'Framework_AreaTable_rgb_July_25.csv'))
	ind = read.csv(file.path(dadir,'ERCindexandfisheryto2018.csv'))
	#check naming
	fG = unique(fis$GAZETTED_NAME)
	aG = unique(ar$DRAINAGE_NAME)

	h = merge(fis,ar,by.x=c('GAZETTED_NAME'),by.y = 'DRAINAGE_NAME') #lose 195 records where GAZETTED_NAME is NA

	uH = unique(h$GAZETTED_NAME)
	out = list()
	m = 0

	com = aggregate(cbind(EFFORT_HOURS,LOG_WEIGHT)~GAZETTED_NAME+LICENCE+YYYY+DRAINAGE.AREA.KM2,data=h,FUN=sum)
ad = aggregate(DRAINAGE.AREA.KM2~YYYY,data=h,FUN=sum)
ad$TAC = ad$DRAINAGE.AREA.KM2*0.69


	comD = aggregate(cbind(EFFORT_HOURS,LOG_WEIGHT)~GAZETTED_NAME+YYYY+DRAINAGE.AREA.KM2,data=h,FUN=sum)

comL = aggregate(cbind(EFFORT_HOURS,LOG_WEIGHT,DRAINAGE.AREA.KM2)~YYYY+LICENCE,data=h,FUN=sum)

#Production by river km from EscapementModel
## east river chester = 31
 load(file.path(project.datadirectory('bio.eels'),'jagsModels', 'escapementModelCPUETYNO2002.rdata'))
 dadir = file.path(project.datadirectory('bio.eels'),'data')
		 inda = read.csv(file.path(dadir,'ERCindexandfisheryto2018.csv'))
		 inda = inda[7:nrow(inda),] 
		 w = which(inda$Year==2002)
		 inda$KgTot[w] <- 536 #Not a highly confident number but better than model likely
		 inda$EscKG[w] = inda$KgTot[w] - inda$ERC_L[w]
		 inda$URS = inda$KgTot / 136.59
	
	A = y$URS
	yrs = data.frame(Ind = 1:23, Yrs=1996:2018)
dN = split(comD,f=comD$GAZETTED_NAME)
meds=list()
m=0
for(i in 1:length(dN)){
	U = dN[[i]]
	dR = unique(U$GAZETTED_NAME)
	outs = list()
	yys = c()
	yy = c()
	for(j in 1:nrow(U)){
		m=m+1
		yR = U[j,'YYYY']
		uU = U[j,'DRAINAGE.AREA.KM2']
		L = U[j,'LOG_WEIGHT']
		P = subset(inda,Year==yR)$URS
		ii = yrs[which(yrs$Yrs==yR),'Ind']
		yys = c(yys,yR)
		outs[[j]] = L / (as.vector(A[ii,,])*as.numeric(uU))
		yy[[j]] = L/(P*as.numeric(uU))
		meds[[m]] = c(median(outs[[j]]),uU,dR,yy[[j]],yR)
		}

	tsHists(x=outs,xaxis=yys,xlab='Year',ylab='Exploitation',main=dR,ylim=c(0,1))
	points(yys,yy,col='blue',pch=16)
	title(dR)
	abline(h=0.69,lwd=2,col='blue')
	abline(h=0.49,lwd=2,lty=2,col='blue')
	savePlot(file.path(project.figuredirectory('bio.eels'),paste('Exploitation',dR,'NO2002.png',sep="")),type='png')
	}
	#savePlot(file.path(project.figuredirectory('bio.eels'),paste('Exploitation','EastRiverChesterNoViolins','NO2002.png',sep="")),type='png')

a=as.data.frame(do.call(rbind,meds))
a$V1=as.numeric(a$V1)
a$V2=as.numeric(a$V2)
a$V4 = as.numeric(a$V4)
a$V5 = as.numeric(a$V5)

g = subset(a,V1>0.69)
aggregate(V1~V3+V2,data=g,FUN=function(x) length(unique(x)))
a$cols = ifelse(a$V5>=2008,'blue','black')
a$pchs = ifelse(a$V5>=2008,8,16)
with(a,plot(log(V2),V1,xlab='log(DRAINAGE.AREA.KM2)',ylab='Exploitation',pch=pchs,col=cols))
abline(h=0.69,col='blue',lwd=2)
abline(h=0.49,col='blue',lwd=2,lty=2)

savePlot(file.path(project.figuredirectory('bio.eels'),paste('ExploitationbyDrainageAreaCols.png',sep="")),type='png')

with(a,plot(log(V2),V4,xlab='log(DRAINAGE.AREA.KM2)',ylab='Exploitation',pch=16))
abline(h=0.69,col='blue',lwd=2)
abline(h=0.49,col='blue',lwd=2)

savePlot(file.path(project.figuredirectory('bio.eels'),paste('ExploitationbyDrainageAreaNonModelKgperKm2.png',sep="")),type='png')


write.csv(a,file=file.path(project.datadirectory('bio.eels'),'data','MedianExploitationByRiver.csv'))


###do the exploitatoin raw with the data to come up the these relationships
median(y$URS[-7,,])*0.69
median(y$URS[,,])*0.69
median(y$URS[-7,,])*0.49 #FSPR50


#Median across full time series of modelled resutls
u30 = median(y$URS[,,])*0.69
u50 = median(y$URS[,,])*0.49

u30R = median(inda$URS,na.rm=T)*0.69
u50R= median(inda$URS,na.rm=T)*0.49

u30BS = 2.55*0.69
u50BS = 2.55*0.49

with(inda,plot(Year,URS,pch=16,ylab='Kg/km2 ERC'))
abline(h=2.339117,col='black')
text(1999,2.45,'Median=2.34')
abline(h=2.589257,col='black')
text(1999,2.65,'Mean=2.59')
abline(h=2.55,col='blue')
text(2006,2.65,'MedianBS=2.55',col='blue')

savePlot(file.path(project.figuredirectory('bio.eels'),paste('MedianMeanKgperKm2.png',sep="")),type='png')

km2 = 0:1000

plot((km2),u30*km2,col='blue',lwd=2,lty=1,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines(km2,u30R*km2,col='blue',lwd=2,lty=4,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines((km2),u50*km2,col='black',lwd=2,lty=1,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines((km2),u50R*km2,col='black',lwd=2,lty=4,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
abline(h=400,lwd=2,lty=3)
legend('topleft',c('TAC SPR50 Modelled','TAC SPR30 Modelled','TAC SPR50 Raw','TAC SPR30 Raw','Constant TAC'),lty=c(1,1,4,4,3),col=c('black','blue','black','blue','black'),bty='n',cex=0.75)
savePlot(file.path(project.figuredirectory('bio.eels'),paste('HabitatTACbySPRRawModelled.png',sep="")),type='png')


plot((km2),u30*km2,col='blue',lwd=2,lty=1,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC or Landings')
lines(km2,u30R*km2,col='blue',lwd=2,lty=4,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines((km2),u50*km2,col='black',lwd=2,lty=1,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines((km2),u50R*km2,col='black',lwd=2,lty=4,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
abline(h=400,lwd=2,lty=3)
legend('topleft',c('TAC SPR50 Modelled','TAC SPR30 Modelled','TAC SPR50 Raw','TAC SPR30 Raw','Constant TAC'),lty=c(1,1,4,4,3),col=c('black','blue','black','blue','black'),bty='n',cex=0.75)
with(comD,points((DRAINAGE.AREA.KM2),LOG_WEIGHT,pch=16,cex=EFFORT_HOURS/max(EFFORT_HOURS)*2))
savePlot(file.path(project.figuredirectory('bio.eels'),paste('HabitatTACbySPRRawModelledwithData.png',sep="")),type='png')


plot((km2),u30*km2,col='blue',lwd=2,lty=1,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC or Landings',xlim=c(0,250),ylim=c(0,425))
lines(km2,u30R*km2,col='blue',lwd=2,lty=4,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines((km2),u50*km2,col='black',lwd=2,lty=1,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines((km2),u50R*km2,col='black',lwd=2,lty=4,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
abline(h=400,lwd=2,lty=3)
legend('topleft',c('TAC SPR50 Modelled','TAC SPR30 Modelled','TAC SPR50 Raw','TAC SPR30 Raw','Constant TAC'),lty=c(1,1,4,4,3),col=c('black','blue','black','blue','black'),bty='n',cex=0.75)
with(comD,points((DRAINAGE.AREA.KM2),LOG_WEIGHT,pch=16,cex=EFFORT_HOURS/max(EFFORT_HOURS)*2))
savePlot(file.path(project.figuredirectory('bio.eels'),paste('HabitatTACbySPRRawModelledwithDataZoomed<250.png',sep="")),type='png')




###early years
km2 = 0:1000

plot((km2),u30*km2,col='blue',lwd=2,lty=1,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC or Landings',xlim=c(0,250),ylim=c(0,450))
lines(km2,u30R*km2,col='blue',lwd=2,lty=4,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines(km2,u30BS*km2,col='blue',lwd=2,lty=6,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')

lines((km2),u50*km2,col='black',lwd=2,lty=1,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines((km2),u50R*km2,col='black',lwd=2,lty=4,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines(km2,u50BS*km2,col='black',lwd=2,lty=6,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')

abline(h=400,lwd=2,lty=3)
legend('topleft',c('TAC SPR50 Modelled','TAC SPR30 Modelled','TAC SPR50 Raw','TAC SPR30 Raw','Constant TAC'),lty=c(1,1,4,4,3),col=c('black','blue','black','blue','black'),bty='n',cex=0.75)
comD$pchs=ifelse(comD$YYYY>=2008,10,16)
comD$cols = ifelse(comD$YYYY>=2008,'blue','black')

with(comD,points((DRAINAGE.AREA.KM2),LOG_WEIGHT,pch=pchs,col=cols,cex=EFFORT_HOURS/max(EFFORT_HOURS)*3))
savePlot(file.path(project.figuredirectory('bio.eels'),paste('ColsHabitatTACbySPRRawModelledwithDataZoomed<250.png',sep="")),type='png')



plot((km2),u30*km2,col='blue',lwd=2,lty=1,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC or Landings')
lines(km2,u30R*km2,col='blue',lwd=2,lty=4,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines(km2,u30BS*km2,col='blue',lwd=2,lty=6,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')

lines((km2),u50*km2,col='black',lwd=2,lty=1,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines((km2),u50R*km2,col='black',lwd=2,lty=4,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')
lines(km2,u50BS*km2,col='black',lwd=2,lty=6,type='l',xlab='Habtitat Area Km2',ylab='Habitat Based TAC')

abline(h=400,lwd=2,lty=3)
legend('topleft',c('TAC SPR50 Modelled','TAC SPR30 Modelled','TAC SPR50 Raw','TAC SPR30 Raw','Constant TAC'),lty=c(1,1,4,4,3),col=c('black','blue','black','blue','black'),bty='n',cex=0.75)
comD$pchs=ifelse(comD$YYYY>=2008,10,16)
comD$cols = ifelse(comD$YYYY>=2008,'blue','black')

with(comD,points((DRAINAGE.AREA.KM2),LOG_WEIGHT,pch=pchs,col=cols,cex=EFFORT_HOURS/max(EFFORT_HOURS)*3))
savePlot(file.path(project.figuredirectory('bio.eels'),paste('ColsHabitatTACbySPRRawModelledwithData.png',sep="")),type='png')

##
comD$HAB=comD$DRAINAGE.AREA.KM2*u30R

comD$Ind = ifelse(comD$HAB<comD$LOG_WEIGHT,1,0)
#Proportion of small rivers with over exploitatoin in some years

ccD = subset(comD,DRAINAGE.AREA.KM2<250)
sum(ccD$Ind)/nrow(ccD)
AMCOOK/bio.eels documentation built on May 20, 2019, 4:13 p.m.