In haddonm/abspatial: A set of functions to assist mapping Abalone Spatial data

knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE)

options(knitr.kable.NA = "",
        knitr.table.format = "pandoc")

wkdir <- "C:/A_Mal/Rcode/abspatialuse/"

datadir <- "C:/A_Mal/Rcode/AbGeo026/geostat/"

options("show.signif.stars"=FALSE,"stringsAsFactors"=FALSE,
        "max.print"=50000,"width"=240)

library(rutilsMH)
#library(MQMF)
library(r4cpue)
library(r4maps)
suppressPackageStartupMessages(library(knitr))
library(captioner)
library(sp)
library(rgdal)
library(abspatial)

tab_nums <- captioner(prefix = "Table")
fig_nums <- captioner(prefix = "Figure")
equ_nums <- captioner(prefix = "Equ.")

l

Introduction

First, obviously we need to setup the required library calls and identify the datafiles.

# ab <- readRDS(paste0(datadir,"AllKud90_2019_05_13.rds"))
# abdat <- getlatlong(ab)
# pick <- which(abdat$fishyear > 2011) # remove earlier less complete GPS data
# ab1 <- droplevels(abdat[pick,])
# save(ab1,file = paste0(datadir,"ab1.RData"))

load(paste0(datadir,"ab1.RData"))
str(ab1,max.level=1)

\newline

Unfortunately, within the zone variable there is no GL zone identified, as it doesn't really relate to particular blocks but we can search out those blocks where any greenlip have been taken over the 7 years from 2012 - 2018 and then find out the relative greenlip to blacklip catches by block. We will first examine catches and cpue by block before we begin to explore potentially finer spatial scales.

\newline

# get the total catch of greenlip by block and year
totg <- tapply(ab1$catchglkg,ab1$blockno,sum,na.rm=TRUE)/1000
# which blocks have reported greenlip catches
pickgblk <- as.numeric(names(which(totg > 0.0)))
pickrec <- which(ab1$blockno %in% pickgblk)
ab2 <- droplevels(ab1[pickrec,])
cat("All > 2011 records ",dim(ab1)[1],"Greenlip records ",dim(ab2)[1],
    " proportion of greenlip records ",round((dim(ab2)[1]/dim(ab1)[1]),4))

\newline

cbg <- tapply(ab2$catchglkg,list(ab2$blockno,ab2$fishyear),sum,na.rm=TRUE)/1000
cbb <- tapply(ab2$catchblkg,list(ab2$blockno,ab2$fishyear),sum,na.rm=TRUE)/1000
totg <- rowSums(cbg,na.rm=TRUE)
totb <- rowSums(cbb,na.rm=TRUE)
# Combined catches of each only in blocks with greenlip
round(cbind(cbg,cbb),1) # Greenlip then Blacklip

\newline Obviously, some blocks are more mixed than others. If we select only those with at least 10% greenlip in total that will capture a large proportion of the total greenlip catch and remove a large proportion of the records

# Identify those blocks where catch may be sufficient to characterize cpue
totals <- cbind(totg,totb,totg/(totg+totb)) # Make matrix of totals
rown <- as.numeric(rownames(totals))
totals <- cbind(totals,rown)  #  PropG is proportion of greenlip in block
colnames(totals) <- c("Greenlip","Blacklip","PropG","Block")
lowerlim <- 0.1     # greenlip must be > lowerlim of all catches in a block
rowpick <- which(totals[,"PropG"] > lowerlim) # rowpick is the row in totals
pickblk <- as.numeric(names(rowpick))  # need names to identify blocks
pickrec <- which(ab2$blockno %in% pickblk) # use only those selected blocks
ab3 <- droplevels(ab2[pickrec,])
totblk <- sum(ab2$catchglkg,na.rm=TRUE)/1000 # what has been removed
totgblk <- sum(ab3$catchglkg,na.rm=TRUE)/1000
mix <- totgblk/totblk
records <- dim(ab3)[1]/dim(ab2)[1]
cat(round(mix,4)," Catch Retained or ",round(totgblk,3)," of ",
    round(totblk,3)," using ",lowerlim,"\n")
cat(dim(ab3)[1]," or ",round(records,4)," of ",dim(ab2)[1]," records, reserved using ",lowerlim,"\n\n")
round(totals[rowpick,],3) # see all of totals to see blocks ignored

\newline

There are some outlying records where a very small catch of greenlip is retained because the catch of blacklip is tiny. These might be best removed from consideration when considering CPUE. Although for the moment no action has been taken. One option would be to ensure that the sum of the total catch within a block was greater than some given total (perhaps 10 tonnes), which would eliminate blocks 40 and 47.

plot1(x=ab3[,"glip_qd_cpue"],y=ab3[,"glip_gps_cpue"],type="p",maxy=0,
      xlabel = "Quota Docket CPUE",ylabel="GPS CPUE")
lines(c(0,500),c(0,500),lwd=2,col=2)

\newline

The distribution of the

pickab <- which(ab3[,"catchglkg"] > 0)
glce <- ab3[pickab,"catchglkg"]/(ab3[pickab,"duration"]/60)
pick <- which(glce < 400)  # constrain plot to < log(6) = 403
glcegps <- ab3[,"glip_gps_cpue"]
pick2 <- which(glcegps < 400)
glceqd <- ab3[,"glip_qd_cpue"]
pick3 <- which(glceqd < 400)

parset(plots=c(3,1))
hist(log(glce[pick]),breaks=seq(-1.75,6,0.25),col=2,
     xlab="log(catchglkg/(duration/60))",main="",panel.first=grid())
hist(log(glcegps[pick2]),breaks=seq(-1.75,6,0.25),col=2,
     xlab="log(GPS cpue)",main="",panel.first=grid())
hist(log(glceqd[pick3]),breaks=seq(-1.75,6,0.25),col=2,
     xlab="log(QD cpue)", main="",panel.first=grid())

\newline

parset()
pickt <- which(ab3$duration < 20)
hist(ab3$duration[pickt],breaks=25,col=2,main="")

calculate the proportion of greenlip.

pickz <- which(ab3[,"catchglkg"] > 0) # cannot take log of zero
glce <- ab3[pickz,"catchglkg"]/(ab3[pickz,"duration"]/60)
propg <- ab3[pickz,"catchglkg"]/(ab3[pickz,"catchglkg"] + ab3[pickz,"catchblkg"])
parset(plots=c(2,1))
outh <- hist(propg,breaks=25,col=2,main="",xlab="Proportion of Greenlip",
             panel.first=grid())

plot1(x=propg,y=log(glce),type="p",defpar=FALSE,cex=0.5,maxy=0,
      ylabel="Log(CPUE)",xlabel="Proportion of Greenlip")
abline(v=0.05,lwd=1,col=2)

Preliminary Standardization

Let's examine a preliminary standardization that tries a few options of data selection for the whole data set. Possible factors for inclusion would include:

1. code - which I take to be some form of diverid
2. avgdpth - which is the average depth from the depth loggers
3. lmphr - linear metres per hour - the average speed of fishing
4. blockno - the statistical block number
5. dist_to_coast - the distance to coast of the centroid of the KUD?
6. fishmonth - the month of fishing
7. fishyear -Finally, the factor which summarizes the index of interest

There are interactions possible but the simple models will do for a start. Without more details regarding the precise spatial structure of fishing, which one can imagine would drive many of the zero catches, we will start by omitting zero catches of greenlip.

# prepare data and check for trends in zeros
pickz <- which(ab3[,"catchglkg"] > 0)
ab4 <- droplevels(ab3[pickz,])
glce <- ab4[,"catchglkg"]/(ab4[,"duration"]/60)
propg <- ab4[,"catchglkg"]/(ab4[,"catchglkg"] + ab4[,"catchblkg"])
# split the propg into 0.001 - 0.2, 0.2 - 0.4, 0.4-0.6, 0.6-0.8, & 0.8 - 0.99, 1.0
propfact <- trunc((propg * 10)/2) # try as a factor first
ab4$LnCE <- NA
ab4$propfact <- NA
ab4$LnCE <- log(glce)
ab4$propfact <- propfact
ab4$propg <- NA
ab4$propg <- propg
ab4$temp <- NA
ab4$temp <- round(ab4$avgtemp)

picktemp <- which((ab4$avgtemp < 25) & (ab4$avgtemp > 8) )
parset()
hist(ab4$avgtemp[picktemp],breaks=20,col=2,main="",xlab="Average Temperature")

The proportion of records with no greenlip varied from about 4.4% in 2012 to a maximum of about 5.25% in 2016 down to a minimum of 2% in 2018.

N <- dim(ab3)[1]
tmp <- droplevels(ab3[-pickz,])
empty <- tapply(tmp$catchglkg,list(tmp$blockno,tmp$fishyear),length)
plot1(x=2012:2018,y=colSums(empty,na.rm=TRUE)/N,
      ylabel="Proportion No Greenlip")
kable(empty)

Now conduct a standardization

pickrec <- which((ab4$propg > 0.1) & (ab4$avgtemp < 25) & (ab4$avgtemp > 8))
labelM <- c("fishyear","propfact","code","blockno","temp")#,"fishmonth") # simple model to start
ab5 <- makecategorical(labelM,ab4[pickrec,])

model <- makemodels(labelM,dependent="LnCE")
#out <- lm(model,data=ab5)
#summary(out)
#anova(out)
out <- standLM(model,ab5)
plotstand(out,bars=TRUE)

qqdiag(out$optModel,plotrug=TRUE,bins=seq(-4,3,0.2))

diagnosticPlot(out,inmodel=5,ab5)

The proportional representation of lower proportions of greenlip in the catches can be seen by plotting the proportion of records found in the different propfact levels. Remembering that propfact 5 equals 100% greenlip and 0 = >0 < 0.2, etc.

pby <- tapply(ab4$propfact,list(ab4$propfact,ab4$fishyear),length)
pbyp <- prop.table(pby,margin=2)
parset()
plot(0:5,pbyp[,1],type="l",lwd=2,col=1,panel.first=grid())
for (i in 2:7) lines(0:5,pbyp[,i],lwd=2,col=i)
legend("topleft",legend=2012:2018,col=1:7,lwd=3,bty="n")

impactplot(out)

tempfact <- getfact(out,"temp")
picktC <- which((ab4$avgtemp > 8) & (ab4$avgtemp < 24))
parset()
plot(ab4$avgtemp[picktC],ab4$LnCE[picktC],type="p",pch=16,col=rgb(1,0,0,1/20))
tapply(ab4$avgtemp[picktC],list(ab4$blockno[picktC],ab4$fishyear[picktC]),mean,na.rm=TRUE)

pickp <- which(ab4$propg == 1.0)
pickp2 <- which(ab4$propg < 0.1)
parset()
#viewpt <- c(143.5,149,-39,-41.5)
viewpt <- c(147.65,148.8,-39.6,-40.6)
#viewpt <- c(147.6,148.4,-40.2,-40.4)
maptas(view=viewpt)
mappoints(ab4[pickp,],inpch=16,incex=0.6)
mappoints(ab4[pickp2,],incol=4,inpch=16,incex=0.6)
legend("topright",c("100% Greenlip","<10% Greenlip"),col=c(2,4),lwd=3,
                    bty="n",cex=1.25)

haddonm/abspatial documentation built on June 7, 2019, 9:54 a.m.