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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.