data-raw/firsttrials.r

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

# library declarations -------------------------------------------------
#library(r4cpue)
#library(r4maps)
library(rutilsMH)
library(diagrams)
#library(MQMF)
library(sp)
library(rgdal)
library(abspatial)
library(knitr)

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

# blacklip 1 Ha grid locations
ab <- readRDS(paste0(datadir,"G1Ha_2019_05_13.rds"))

#ab <- readRDS(paste0(datadir,"AllKud90_2018_10_17.rds"))

#ab <- load(paste0(datadir,"fullspatial_2018-10-23.RData"))

str(ab,max.level = 2)


abdat <- getlatlong(ab)

str(abdat,max.level=1)

# plot actaeons -----------------------------------------------------------
bruny <- c(147.075,147.45,-43.05,-43.55) # left, right, up and down
blk13cde <- c(146.9,147.025,-43.5,-43.625)
actaeons <- c(146.98,147.02,-43.51,-43.59)
invar <- "blkg"


label <- c(paste0("blkg",c(2012:2018,"total")))
plotprep(width=6,height=10, newdev=FALSE)
plotgrid(actaeons,abdat,plots=c(2,2),label,index=1:4)



addat <- getpoints(actaeons,abdat)
ans <- summarygrid(addat,invar)
nactive <- dim(addat)
whchvar <- "blkg"

label <- paste0(whchvar,c(2012:2018))

lab1 <- label[1]
for (i in 2:7) { # i = 2
  lab2 <- label[i]
  
  pick1 <- which(addat[,lab1] > 0)
  pick2 <- which(addat[,lab2] > 0)
  overlapoid <- which(addat[pick1,"oid"] %in% addat[pick2,"oid"])
  N1 <- length(pick1)
  N2 <- length(pick2)
  No <- length(overlapoid)
  cat(label[i],N1-No,No,N2-No,"\n")
}

# cumulative catch --------------------------------------------------------
invar <- "blkg2018"
pick <- which(addat[,invar] > 0)
nall <- length(pick)

bl12 <- sort(addat[pick,invar],decreasing=TRUE)
cbl12 <- cumsum(bl12)/sum(bl12)
plotprep(newdev=TRUE)
plot((1:nall)/nall,cbl12,type="l",lwd=2,panel.first=grid(),
     xlab="Proportion of Records",ylab="Proportion of Total Catch")
abline(h=c(0.2,0.5,0.75,0.9),col=3)
abline(v=c(seq(0.1,0.9,0.1)),col="grey",lty=2)




pick <- which(addat[,"blkg2012"] < 100)
FMR::plot2(newdev=FALSE)
hist(addat[pick,"blkg2012"],breaks=30,main="",col=2)
outH <- hist(log(addat[pick,"blkg2012"]),breaks=30,main="",col=2)


mids <- outH$mids
cw <- mids[2]-mids[1]
counts <- outH$counts
# Guess some normal parameters and plots those curves on histogram
tot <- sum(counts,na.rm=TRUE)
av <- c(3.9,5.7)   # the means
stdev <- c(0.75,0.6)  # the standardi deviations
prop1 <- 0.45         # the proportion of observations in cohort 1
cohort1 <- (tot*prop1*cw)*dnorm(mids,av[1],stdev[1]) # expected counts for
cohort2 <- (tot*(1-prop1)*cw)*dnorm(mids,av[2],stdev[2]) # each cohort
# the (tot*prop1*cw) scales the likelihoods to suit the 2mm class width
lines(mids,cohort1,lwd=2,col=1)
lines(mids,cohort2,lwd=2,col=4)

mnnegLL <- function(pars,obs,vals,midval=TRUE) { # when using mid-points of class widths 
  tot <- sum(obs)  #  tot and cw are for scaling
  cw <- vals[2] - vals[1]
  n <- (length(pars)+ 1)/3 # number of cohorts
  props <- tail(pars,(n-1)) # input n-1 proportions, 
  props <- c(props,(1-sum(props))) # the last is the difference
  nval <- length(vals)
  if (midval) {
    freqs <- numeric(nval) # make a vector of zeros to hold freqs
    for (i in 1:n) # step through the cohorts
      freqs <- freqs + (tot*props[i]*cw)*dnorm(vals,pars[i],pars[i+n])
  } else {
    freqs <- numeric(nval-1) # make a vector of zeros to hold freqs
    for (i in 1:n) { # step through the cohorts
      tmpcump <- (tot*props[i])*pnorm(vals,pars[i],pars[i+n])
      freqs <- freqs + (tmpcump[2:nval] - tmpcump[1:(nval-1)]) # lower from upper
    } 
  } # end of if statement
  # now return -veLL for multinomial; does not use dmultinom
  return(-sum(obs * log(freqs/sum(freqs,na.rm=TRUE)),na.rm=TRUE))
} # end of mnnegLL

av <- c(3.9,5.7)   # the means
stdev <- c(0.75,0.6)  # the standardi deviations
prop1 <- 0.45         # the proportion of observations in cohort 1
pars <- (c(av,stdev,prop1))
# use nlm to find optimum parameters when using size class mid and bounds
bestmod <- nlm(f=mnnegLL,p=pars,obs=counts,vals=mids, 
               typsize=MQMF::magnitude(pars),iterlim=1000)
outfit(bestmod); cat("\n")
round((bestmod$estimate),6)

pick <- which(addat[,"blkg2012"] > 0)
FMR::plot2(newdev=FALSE)
hist(addat[pick,"blkg2012"],breaks=30,main="",col=2)
outH <- hist(log(addat[pick,"blkg2012"]),breaks=30,main="",col=2)
mids <- outH$mids
cw <- mids[2]-mids[1]
counts <- outH$counts
# Guess some normal parameters and plots those curves on histogram
tot <- sum(counts,na.rm=TRUE)

pars <- bestmod$estimate # best estimate using mid-points
cohort1 <- (tot*pars[5]*cw)*dnorm(mids,pars[1],pars[3]) 
cohort2 <- (tot*(1-pars[5])*cw)*dnorm(mids,pars[2],pars[4]) 
lines(mids,cohort1,lwd=2,col=1) # plot the mid-point estimates
lines(mids,cohort2,lwd=2,col=4)
# Plot Hexagon Catches -----------------------------------------------------
# Maria Island ---------------------------------------------------

plotprep(width=5,height=5.0,newdev = FALSE)
maria <- c(148,148.1,-42.675,-42.75)
maptas(view=maria,gridon=0.1) 
b12 <- mapgrid(view=maria,abdat,x="blkg2012",mult=3,col=2) 
b13 <- mapgrid(view=maria,abdat,x="blkg2013",mult=3,col=4,loconly=TRUE) 

plotprep(width=5,height=5.0,newdev = FALSE)
maria <- c(148,148.1,-42.675,-42.75)
maptas(view=maria,gridon=0.1) 
addat <- getpoints(maria,abdat)
mappoints(addat)
b12 <- mapgrid(view=maria,abdat,x="blkg2012",mult=3,col=2) 
b13 <- mapgrid(view=maria,abdat,x="blkg2013",mult=3,col=4) 
btot <- mapgrid(view=maria,abdat,x="blkgtotal",mult=3,col=3) 

getvar <- "blkg"
label <- c(paste0(getvar,c(2012:2018,"total")))
num <- length(label)
plotprep(width=6,height=10, newdev=FALSE)
par(mfcol=c(4,2),mai=c(0.3,0.3,0.05,0.05),oma=c(1.0,1.0,0.0,0.0)) 
par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)  
maria <- c(148,148.1,-42.675,-42.75)
for (i in 1:num) {
  maptas(view=maria,gridon=0.1,defineplot=FALSE,maintitle = label[i]) 
  mapgrid(view=maria,abdat,inpch=21,x=label[i],mult=3,col=1,infill=2) 
}
mtext("Longitude",side=1,outer=TRUE,line=0.0,cex=1.1,font=7)
mtext("Latitude",side=2,outer=TRUE,line=0.0,cex=1.1,font=7)
addat <- getpoints(maria,abdat)



summarygrid(addat,getvar)


# Growth in area fished----------------------------------------------------
addat <- getsubblock(c("13C","13D","13E"),abdat)

rown <- dim(addat)[1]
addat[,"oid"] <- as.numeric(addat[,"oid"])
pick <- grep("blkg",colnames(addat))
yrs <- 2012:2018
nyr <- length(yrs)
columns <- c("Year","Total","P(tot)","Unfound","New")
knownarea <- matrix(0,nrow=nyr,ncol=length(columns),dimnames=list(yrs,columns))

pickr <- which(addat[,pick[1]] > 0)
found <- addat[pickr,"oid"]
unfound <- addat[-pickr,"oid"]
knownarea[1,] <- c(yrs[1],length(found),length(found)/rown,length(unfound),
                   length(found))
for (yr in 2:nyr) {  # yr=2
  pickr <- which(addat[,pick[yr]] > 0)
  yrdat <- addat[pickr,]
  newoid <- which(yrdat[,"oid"] %in% unfound)
  found <- c(found,yrdat[newoid,"oid"])
  unfound <- unfound[-which(unfound %in% yrdat[,"oid"])]
  knownarea[yr,] <- c(yrs[yr],length(found),length(found)/rown,
                      length(unfound),length(newoid))
}
kable(knownarea)

addat <- getsubblock(c("13C","13D","13E"),abdat)

leftlong <- 146.86;  rightlong <- 147.03
uplat <- -43.51;  downlat <- -43.64
E13 <- c(leftlong,rightlong,uplat,downlat)
plotprep(width=7,height=9, newdev=FALSE)
maptas(view=E13,gridon=0.1) 
mapgrid(view=NA,addat,inpch=21,x="blkgtotal",mult=3,col=1,infill=2) 

addat <- getsubblock(c("13A","13B"),abdat)

leftlong <- 146.55;  rightlong <- 146.9
uplat <- -43.54;  downlat <- -43.65
W13 <- c(leftlong,rightlong,uplat,downlat)

plotprep(width=8,height=5, newdev=FALSE)
maptas(view=W13,gridon=0.1) 
mapgrid(view=NA,addat,inpch=21,x="blkgtotal",mult=3,col=1,infill=2) 


addat <- getblock(c(24),abdat)

leftlong <- 147.89;  rightlong <- 148.175
uplat <- -42.5;  downlat <- -42.765
E24 <- c(leftlong,rightlong,uplat,downlat)

plotprep(width=6,height=7, newdev=FALSE)
maptas(view=E24,gridon=0.1) 
mapgrid(view=NA,addat,inpch=21,x="blkgtotal",mult=3,col=1,infill=2) 

# overlap analyses ---------------------------------------------------------

data(abdat)
addat <- getsubblock(c("13C","13D","13E"),abdat)
addat[,"oid"] <- as.numeric(addat[,"oid"])
whchvar <- "blkg"
yrs <- 2012:2018

addat <- getsubblock(c("13A","13B"),abdat)
addat[,"oid"] <- as.numeric(addat[,"oid"])
whchvar <- "blkg"
yrs <- 2012:2018



answer <- overlap(addat,whichvar=whchvar,yrs=yrs)
answer

pickdat <-  which((addat[,12] > 60.0))
answer2 <- overlap(addat[pickdat,],whichvar=whchvar,yrs=yrs)
answer2

repeats <- apply(addat[,5:11],1,countgtzero)
counts <- as.matrix(table(repeats))
cbind(counts,counts/sum(counts))

addat <- getblock(12,abdat)
addat[,"oid"] <- as.numeric(addat[,"oid"])

answer <- overlap(addat)
answer
getrepeats(addat)


for (blk in 14:30) { # blk <- 6
  addat <- getblock(blk,abdat)
  print(blk)
  print(getrepeats(addat))
  
}



label <- paste0(whchvar,yrs)
pickcol <- grep(whchvar,colnames(addat)) # which columns contain year-catches
nyr <- length(yrs)

# 

leftlong <- 146.86;  rightlong <- 147.03
uplat <- -43.51;  downlat <- -43.64
E13 <- c(leftlong,rightlong,uplat,downlat)



# check there are no empty records in the available data 
totcatch <- rowSums(addat[,pickcol],na.rm=TRUE)
length(totcatch) # Number of records available for selected area
pick <- which(totcatch > 0.0)
length(pick)   # number of positive records available for area. 

whchvar <- "blkg"
yrs <- 2012:2018
label <- paste0(whchvar,yrs)



answer <- overlap(addat,whichvar=whchvar,yrs=yrs)
answer

pickdat <-  which((addat[,12] > 60.0))
answer2 <- overlap(addat[pickdat,],whichvar=whchvar,yrs=yrs)
answer2


dat <- as.numeric(as.matrix(addat[,5:11]))

pickdat <- which((addat[,12] <= 200) & (addat[,12] > 0.0))

                  
plotprep(width=7,height=4,newdev=FALSE)
hist(addat[pickdat,12],col=2,breaks=40)




  



  

  
  dat <- p1dat
  year <- yr+1
  plotprep(width=7,height=9, newdev=FALSE)
  maptas(view=E13,gridon=0.1) 
  mappoints(dat,inpch=21,incol=2,infill=3) 
  mappoints(dat[dat[,label[year]]>0,],inpch=21,incol=2,infill=2,incex=0.75) 
  mappoints(ydat[ydat[,label[yr]]>0,],inpch=21,incol=3,infill=3,incex=0.5) 
# getarea functions -------------------------------------------------------

#' Title
#'
#' @param inzone 
#' @param indat 
#'
#' @return
#' @export
#'
#' @examples
getzone <- function(inzone,indat) {
  if (any(inzone %ni% c("BS","E","N","W")))
    stop(cat("All members of ",inzone,"must be in c('BS','E','N','W') \n"))
  pick <- which(indat[,"zone"] %in% inzone)
  addat <- droplevels(indat[pick,])
  return(invisible(addat))
} # end of getzone


addat <- getzone(c("W","E"),abdat)
dim(addat)

addat <- Mgetblock(13,abdat)


# greenlip 1 Ha grid locations---------------------------------------------
abg <- readRDS(paste0(datadir,"G1HaGL_2018_10_23.rds"))

abgat <- getlatlong(abg) 

pick <- which(abgat$zone %in% c("BS","N"))
abgat <- abgat[pick,]
dim(abgat)

plotprep(width=8,height=4)
plot(abgat[,"long"],abgat[,"lat"],type="p",pch=16,cex=0.2)

plotprep(width=8,height=4)
leftlong <- 143;  rightlong <- 149 
uplat <- -39;  downlat <- -41.25 
maptas(c(leftlong,rightlong,uplat,downlat),gridon=1.0) 
mappoints(abgat,Long="long",Lat="lat") 


# Mixed Species CPUE -------------------------------------------------------

unique(abdat$zone)

pick <- which(abdat$fishyear > 2011)
ab1 <- droplevels(abdat[pick,])

cgb <- tapply(ab1$catchglkg,list(ab1$blockno,ab1$fishyear),sum,na.rm=TRUE)/1000
tot <- rowSums(cgb,na.rm=TRUE)
pickr <- which(tot > 0)
round(cbind(cgb[pickr,],tot[pickr]),3)

colSums(cgb,na.rm=TRUE)

cbb <- tapply(ab1$catchblkg,list(ab1$blockno,ab1$fishyear),sum,na.rm=TRUE)/1000
totb <- rowSums(cbb,na.rm=TRUE)

round(cbind(tot[pickr],totb[pickr]),3)

plotprep()
plot1(x=tot[pickr],y=totb[pickr],type="p",cex=1.2)

pickgblk <- c(1,2,3,4,5,31,32,33,35,36,38,39,48,49)



 dev.new(height=6.0,width=6,noRStudioGD = TRUE)
 view=c(143.5,149,-39,-44.0)
 maptas(view=view)
 long <- c(147.4,147.8,148,148.2)
 lat <- c(-43.6,-43.6,-43.3,-43.2)
 indata <- cbind(long,lat)
 mappoints(view,indata,incex=1.5)



# count up total oids in a region

plotprep()
outh <- inthist(addat[,"repeat"],col=2,width=0.8,border=3,lwd=2)



plotprep(width=6 height=8)


pick <- which(abdat$oid > 0)
length(pick)