p = bio.lobster::load.environment()
RLibrary("plyr","lattice","glmmADMB","ggplot2")
#LATEST DATA EXPORT FROM FSRS DATABASE:
#lobster.db("fsrs.redo")
lobster.db("fsrs")
#recruitment.trap.db('raw.redo',p=p)
FSRS.dat<-fsrs
FSRS.dat$VES_DATE<-paste(FSRS.dat$VESSEL_CD,FSRS.dat$HAUL_DATE,sep='.')
FSRS.dat$SYEAR<-FSRS.dat$HAUL_YEAR
FSRS.dat$HAUL_DATE<-as.Date(FSRS.dat$HAUL_DATE)
FSRS.dat$SYEAR[FSRS.dat$LFA%in%c("33","34")]<-as.numeric(substr(FSRS.dat$S_LABEL[FSRS.dat$LFA%in%c("33","34")],6,9))
FSRS.dat<-subset(FSRS.dat,SOAK_DAYS<6) # Remove soak days greater than 5, do not iclude berried females
FSRS.dat$HAUL_DATE<-as.Date(FSRS.dat$HAUL_DATE)
# this section is to deal with the fact that there are uneven binning going on for the different size categories
# it creates a pseudo CL (the mid point of each size category)
scd<-read.csv(file.path( project.datadirectory("lobster"), "data","inputs","FSRS_SIZE_CODES.csv"))
scd$LENGTH<-rowMeans(scd[c("MIN_S","MAX_S")])
FSRS.dat<-merge(FSRS.dat,scd[c("SIZE_CD","LENGTH")])
wa<-c(0.000608, 0.001413, 0.00482)
wb<-c(3.058, 2.875, 2.638)
FSRS.dat$WEIGHT<-NA
for(i in 1:3){
FSRS.dat$WEIGHT[FSRS.dat$SEX==i]<-FSRS.dat$LENGTH[FSRS.dat$SEX==i]^wb[i]*wa[i]
}
## Aggregate by unique vessal and day, summerizing total traps, legals and shorts
trap.lst<-lapply(with(FSRS.dat,tapply(TRAP_NO,VES_DATE,unique)),length)
trap.dat<-data.frame(VES_DATE=names(trap.lst),TOTAL_TRAPS=as.vector(unlist(trap.lst)))
short.lst<-with(subset(FSRS.dat,SHORT==1),tapply(TRAP_NO,VES_DATE,length))
short.dat<-data.frame(VES_DATE=names(short.lst),SHORTS=short.lst)
legal.lst<-with(subset(FSRS.dat,SHORT==0),tapply(TRAP_NO,VES_DATE,length))
legal.dat<-data.frame(VES_DATE=names(legal.lst),LEGALS=legal.lst)
recruit.lst<-with(subset(FSRS.dat,SHORT==1&SIZE_CD>7),tapply(WEIGHT,VES_DATE,sum))
recruit.dat<-data.frame(VES_DATE=names(recruit.lst),RECRUITS=recruit.lst)
legalbm.lst<-with(subset(FSRS.dat,SHORT==0),tapply(WEIGHT,VES_DATE,sum))
legalbm.dat<-data.frame(VES_DATE=names(legalbm.lst),BIOMASS=legalbm.lst)
FSRS_1.dat <- aggregate(cbind(VESSEL_CD, HAUL_DATE, DEPTH, LFA, LFA_GRID, TEMP, LAT_DD, LONG_DD, HAUL_YEAR, SYEAR)~VES_DATE,data=FSRS.dat,mean,na.rm=T)
FSRS_2.dat<-merge(trap.dat,merge(recruit.dat,merge(legalbm.dat,merge(short.dat,legal.dat,all=T),all=T),all=T),all=T)
FSRS_2.dat$SHORTS[is.na(FSRS_2.dat$SHORTS)]<-0
FSRS_2.dat$LEGALS[is.na(FSRS_2.dat$LEGALS)]<-0
FSRSvesday<-merge(FSRS_1.dat,FSRS_2.dat,all.x=T)
# Create column for week and day of season (WOS, DOS)
lfas<-unique(FSRSvesday$LFA[!is.na(FSRSvesday$LFA)])
FSRSvesday$WOS<-NA
FSRSvesday$DOS<-NA
for(a in 1:length(lfas)){
season<-sort(unique(FSRSvesday$SYEAR[FSRSvesday$LFA==lfas[a]]))
for(i in 1:length(season)){
FSRSvesday$WOS[FSRSvesday$SYEAR==season[i]&FSRSvesday$LFA==lfas[a]]<-floor((FSRSvesday$HAUL_DATE[FSRSvesday$SYEAR==season[i]&FSRSvesday$LFA==lfas[a]]-min(FSRSvesday$HAUL_DATE[FSRSvesday$SYEAR==season[i]&FSRSvesday$LFA==lfas[a]]))/7)+1
FSRSvesday$DOS[FSRSvesday$SYEAR==season[i]&FSRSvesday$LFA==lfas[a]]<-FSRSvesday$HAUL_DATE[FSRSvesday$SYEAR==season[i]&FSRSvesday$LFA==lfas[a]]-min(FSRSvesday$HAUL_DATE[FSRSvesday$SYEAR==season[i]&FSRSvesday$LFA==lfas[a]])+1
}
}
str(FSRSvesday)
subareas<-read.csv(file.path( project.datadirectory("bio.lobster"), "data","inputs","LFA2733subarea.csv"))
FSRSvesday<-merge(FSRSvesday,subareas,all.x=T)
write.csv(FSRSvesday,file.path( project.datadirectory("bio.lobster"), "data","products","FSRSrectraps.csv"),row.names=F)
FSRSvesday<-read.csv(file.path( project.datadirectory("bio.lobster"), "data","products","FSRSrectraps.csv"))
#plot raw data on shorts
b = aggregate(cbind(SHORTS,TOTAL_TRAPS)~LFA+SYEAR,data=FSRSvesday,FUN=sum)
b$CPUE= b$SHORTS / b$TOTAL_TRAPS
p <- ggplot()
p <- p + geom_point(data = b,
aes(y = CPUE, x = SYEAR),
shape = 16,
size = 3) + xlab("Year") + ylab("Lobsters / Trap") + theme(text = element_text(size=15)) +
theme_bw() + geom_line(data = b, aes(x = SYEAR, y = CPUE), colour = "black")+
facet_wrap( ~LFA, ncol=2,scales = "fixed")
p
#plot temperature across years
FSRSvesday$TEMP[FSRSvesday$TEMP==-99] <- NA
b = aggregate(TEMP~LFA+SYEAR+DOS,data=FSRSvesday,FUN=mean)
lf = unique(b$LFA)
for(i in lf) {
j = subset(b,LFA==i)
plot(1,1,type='n',xlim=c(min(j$DOS),max(j$DOS)),ylim=c(min(j$TEMP),max(j$TEMP)),main=i)
jj = unique(j$SYEAR)
m=0
gg = rev(shadesOfGrey(length(jj)))
for(l in jj) {
m=m+1
with(subset(j,SYEAR==l),lines(DOS,TEMP,col=gg[m]))
}
x11()
}
#FSRSvesday$HAUL_DATE<-as.Date(FSRSvesday$HAUL_DATE)
### Checking frequency of missing temp data
#missing.temps.percent<-lapply(lapply(sapply(2001:2014,function(x){1-with(subset(FSRSvesday,SYEAR==x&TEMP>(-5)),tapply(SYEAR,LFA,length))/with(subset(FSRSvesday,SYEAR==x),tapply(SYEAR,LFA,length))}),'*',100),round)
#names(missing.temps.percent)<-2001:2014
#missing.temps.percent
#-----------------------------------LFA27
#------------------------------North
LFA27north<-subset(FSRSvesday,subarea=='27 North')
LFA27n.sm<-FSRSmodel(LFA27north, response="SHORTS")
#LFA27n.lm<-FSRSmodel(LFA27north, response="LEGALS")
#------------------------------South
LFA27south<-subset(FSRSvesday,subarea=='27 South')
LFA27s.sm<-FSRSmodel(LFA27south, response="SHORTS")
#LFA27s.lm<-FSRSmodel(LFA27south, response="LEGALS")
#-------------------------------LFA28
LFA28<-subset(FSRSvesday,LFA==28)
LFA28.sm<-FSRSmodel(LFA28, response="SHORTS")
#LFA28.lm<-FSRSmodel(LFA28, response="LEGALS")
#-------------------------------LFA29
LFA29<-subset(FSRSvesday,LFA==29)
LFA29.sm<-FSRSmodel(LFA29, response="SHORTS")
#LFA29.lm<-FSRSmodel(LFA29, response="LEGALS")
#-------------------------------LFA30
LFA30<-subset(FSRSvesday,LFA==30)
LFA30.sm<-FSRSmodel(LFA30, response="SHORTS")
#LFA30.lm<-FSRSmodel(LFA30, response="LEGALS")
#-------------------------------LFA31A
LFA31A<-subset(FSRSvesday,LFA==31.1)
LFA31A.sm<-FSRSmodel(LFA31A, response="SHORTS")
#LFA31A.lm<-FSRSmodel(LFA31A, response="LEGALS")
#-------------------------------LFA31B
LFA31B<-subset(FSRSvesday,LFA==31.2)
LFA31B.sm<-FSRSmodel(LFA31B, response="SHORTS")
#LFA31B.lm<-FSRSmodel(LFA31B, response="LEGALS")
#-------------------------------LFA32
LFA32<-subset(FSRSvesday,LFA==32)
LFA32.sm<-FSRSmodel(LFA32, response="SHORTS")
#LFA32.lm<-FSRSmodel(LFA32, response="LEGALS")
#-----------------------------------LFA33
#------------------------------East
LFA33east<-subset(FSRSvesday,subarea=='33 East')
LFA33e.sm<-FSRSmodel(LFA33east, response="SHORTS",interaction=T)
#LFA33e.lm<-FSRSmodel(LFA33east, response="LEGALS",interaction=T)
#------------------------------West
LFA33west<-subset(FSRSvesday,subarea=='33 West')
LFA33w.sm<-FSRSmodel(LFA33west, response="SHORTS",interaction=T)
#LFA33w.lm<-FSRSmodel(LFA33west, response="LEGALS",interaction=T)
#-------------------------------LFA34
LFA34<-subset(FSRSvesday,LFA==34)
LFA34.sm<-FSRSmodel(LFA34, response="SHORTS",interaction=T)
#LFA34.lm<-FSRSmodel(LFA34, response="LEGALS",interaction=T)
#________________________
#------------------------------load previous runs
# FSRSvesday<-read.csv(file.path( project.datadirectory("bio.lobster"), "data","products","FSRSrectraps.csv"))
LFA27north<-subset(FSRSvesday,subarea=='27 North')
LFA27south<-subset(FSRSvesday,subarea=='27 South')
LFA28<-subset(FSRSvesday,LFA==28)
LFA29<-subset(FSRSvesday,LFA==29)
LFA30<-subset(FSRSvesday,LFA==30)
LFA31A<-subset(FSRSvesday,LFA==31.1)
LFA31B<-subset(FSRSvesday,LFA==31.2)
LFA32<-subset(FSRSvesday,LFA==32)
LFA33east<-subset(FSRSvesday,subarea=='33 East')
LFA33west<-subset(FSRSvesday,subarea=='33 West')
LFA27n.sm<-FSRSmodel(LFA27north, response="SHORTS",redo=F)
#LFA27n.lm<-FSRSmodel(LFA27north, response="LEGALS",redo=F)
LFA27s.sm<-FSRSmodel(LFA27south, response="SHORTS",redo=F)
#LFA27s.lm<-FSRSmodel(LFA27south, response="LEGALS",redo=F)
LFA28.sm<-FSRSmodel(LFA28, response="SHORTS",redo=F)
#LFA28.lm<-FSRSmodel(LFA28, response="LEGALS",redo=F)
LFA29.sm<-FSRSmodel(LFA29, response="SHORTS",redo=F)
#LFA29.lm<-FSRSmodel(LFA29, response="LEGALS",redo=F)
LFA30.sm<-FSRSmodel(LFA30, response="SHORTS",redo=F)
#LFA30.lm<-FSRSmodel(LFA30, response="LEGALS",redo=F)
LFA31A.sm<-FSRSmodel(LFA31A, response="SHORTS",redo=F)
#LFA31A.lm<-FSRSmodel(LFA31A, response="LEGALS",redo=F)
LFA31B.sm<-FSRSmodel(LFA31B, response="SHORTS",redo=F)
#LFA31B.lm<-FSRSmodel(LFA31B, response="LEGALS",redo=F)
LFA32.sm<-FSRSmodel(LFA32, response="SHORTS",redo=F)
#LFA32.lm<-FSRSmodel(LFA32, response="LEGALS",redo=F)
LFA33e.sm<-FSRSmodel(LFA33east, response="SHORTS",redo=F)
#LFA33e.lm<-FSRSmodel(LFA33east, response="LEGALS",redo=F)
LFA33w.sm<-FSRSmodel(LFA33west, response="SHORTS",redo=F)
#LFA33w.lm<-FSRSmodel(LFA33west, response="LEGALS",redo=F)
#LFA34.sm<-FSRSmodel(LFA34, response="SHORTS",redo=F)
#LFA34.lm<-FSRSmodel(LFA34, response="LEGALS",redo=F)
# plot model fit in most recent year
FSRSmodel.3dplot(LFA27n.sm$model,response="SHORTS")
FSRSmodel.3dplot(LFA27s.sm$model,response="SHORTS")
FSRSmodel.3dplot(LFA28.sm$model,response="SHORTS")
FSRSmodel.3dplot(LFA29.sm$model,response="SHORTS")
FSRSmodel.3dplot(LFA30.sm$model,response="SHORTS")
FSRSmodel.3dplot(LFA31A.sm$model,response="SHORTS")
FSRSmodel.3dplot(LFA31B.sm$model,response="SHORTS")
FSRSmodel.3dplot(LFA32.sm$model,response="SHORTS")
FSRSmodel.3dplot(LFA33e.sm$model,response="SHORTS")
FSRSmodel.3dplot(LFA33w.sm$model,response="SHORTS")
# compile results
LFA27n.sm$pData$Area<-'27N'
LFA27s.sm$pData$Area<-'27S'
LFA28.sm$pData$Area<-'28'
LFA29.sm$pData$Area<-'29'
LFA30.sm$pData$Area<-'30'
LFA31A.sm$pData$Area<-'31A'
LFA31B.sm$pData$Area<-'31B'
LFA32.sm$pData$Area<-'32'
LFA33e.sm$pData$Area<-'33E'
LFA33w.sm$pData$Area<-'33W'
shorts<-rbind(LFA27n.sm$pData,LFA27s.sm$pData,LFA28.sm$pData,LFA29.sm$pData,LFA30.sm$pData,LFA31A.sm$pData,LFA31B.sm$pData,LFA32.sm$pData,LFA33e.sm$pData,LFA33w.sm$pData)
write.csv(shorts,file.path( project.datadirectory("bio.lobster"), "data","products","FSRSmodelresultsSHORT.csv"),row.names=F)
shorts<-read.csv(file.path( project.datadirectory("bio.lobster"), "data","products","FSRSmodelresultsSHORT.csv"))
shorts.2829<-subset(shorts,Area %in% c("28","29"))
# shorts
pdf(file.path(project.figuredirectory('bio.lobster'),'FSRS.LFA27-33.2017.pdf'),width=8,height=10)
p <- ggplot()
p <- p + geom_point(data = shorts.2829,
aes(y = mu, x = YEAR),
shape = 16,
size = 3)
p <- p + xlab("Year") + ylab("Lobsters / Trap")
p <- p + theme(text = element_text(size=15)) + theme_bw()
p <- p + geom_line(data = shorts.2829,
aes(x = YEAR, y = mu),
colour = "black")
p <- p + geom_ribbon(data = shorts.2829,
aes(x = YEAR,
ymax = ub,
ymin = lb ),
alpha = 0.5)
p <- p + facet_wrap( ~Area, ncol=1,scales = "fixed")
print(p)
dev.off()
####33E and 33W
p33 = subset(shorts,Area %in% c("33W","33E"))
pdf(file.path(project.figuredirectory('bio.lobster'),'FSRS.LFA33.2017.pdf'),width=8,height=10)
p <- ggplot()
p <- p + geom_point(data = p33,
aes(y = mu, x = YEAR),
shape = 16,
size = 3)
p <- p + xlab("Year") + ylab("Lobsters / Trap")
p <- p + theme(text = element_text(size=15)) + theme_bw()
p <- p + geom_line(data = p33,
aes(x = YEAR, y = mu),
colour = "black")
p <- p + geom_ribbon(data = p33,
aes(x = YEAR,
ymax = ub,
ymin = lb ),
alpha = 0.5)
p <- p + facet_wrap( ~Area, ncol=1,scales = "fixed")
print(p)
savePlot(file.path(project.figuredirectory('bio.lobster'),'FSRS.LFA33.2017.png'),type='png')
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.