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

knitr::opts_chunk$set(
    echo = FALSE,
    message = FALSE,
    warning = FALSE)
  options(knitr.kable.NA = "",
          knitr.table.format = "pandoc")

datadir <- "C:/A_Mal/Rcode/AbGeo026/geostat/" # change this to suit
options("show.signif.stars"=FALSE,"stringsAsFactors"=FALSE,
        "max.print"=50000,"width"=240)

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

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

l

Introduction

The Tasmanian abalone GPS data-logger data (Mundy, 2011) provides geo-referenced information about the location of each dive tender every ten seconds for the duration of each dive event. In addition, each dive event also has an associated temperature profile. The raw data consists of giga-bytes of geo-referenced dive-track data and associated catches. As a preliminary synthesis of this raw data it has been summarized in two ways:

1) kernel density estimates have been calculated for the effort in each dive event (including the 25th, 50th, 75th and 90th percentiles of the total implied area of each dive), and

2) after placing a 1-hectare hexagonal grid around the coastline of Tasmania the amount of effort (in 10 second gps points) that has occurred in each one is estimated, and the amount of catch taken in proportion to the effort relative to in each dive event included.

The two types of information derived from the gps logger data can be used to answer different questions about the progress of the fishery and the state of the abalone stock in each region.

The first objective of FRDC project 2017-026 (Mundy et al., 2018) was to "Characterize the statistical properties, coherence, interpret-ability and assumptions of spatial and classic indicators of fishery performance." The intent here was to begin to address that objective. We will begin with the Hexagon grid data (hex) and then consider the kernel density data (KUD - Vernal Density Utilization Distribution).

The hex data is available for both the blacklip and the greenlip stocks around Tasmania. This is contained within three files:

\newline

r tab_nums("T1", caption= " The available data files containing the KUD and the hex data. The dataes within each filename relate to the date when the data eas extracted from the central database. Those files with Kud in their names relate to Kernal Densities, with the associated number depicting the contour, hence AllKud90 is the 90% kernal density for all records. The G1Ha files are the 1 hectare grids for blacklip and greenlip (GL). Finally, the CatchEffortData contains the docket book CPUE data for comparison with the GPS data.")

files <- dir(datadir)
kable(files[c(2:4,7,10:12)], col.names=c("Files"))

\newline

l

Hectare Grid Data

We will step through the three hexagon data sets and illustrate their properties. The blacklip abalone hexagon data set has information about 58808 hexagon across seven years from 2012 - 2018.

gwide <- readRDS(paste0(datadir,"g1hawide_2019_06_07.rds"))
glong <- readRDS(paste0(datadir,"g1halong_2019_06_07.rds"))

dim(gwide)

dim(glong)
pick <- which(glong$blcatch_kg > 0)
length(unique(glong$oid[pick])) # blacklip only

sum(glong$blcatch_kg, na.rm=TRUE)/1000 
sum(gwide$blkgtotal, na.rm=TRUE)/1000

# blacklip 1 Ha grid locations
filenameV <- "G1HaV_2019_05_13.rds"
filename <- "G1Ha_2019_05_13.rds"
ab <- readRDS(paste0(datadir,filename))#get spatial obj
abdat <- getlatlong(ab)     # convert to data.frame + lat long
abdat$occur <- apply(abdat[,5:11],1,countgtzero)

\newline

r tab_nums("T2", caption= paste0(" The contents of ",filename,". Prefix blkg relate to blacklip catch in kg, prefix days is the days of fishing effort in a hexagon, mins is the total effort in each hex, divers is how many divers were active in a hex, kgday is an approximate CPUE, oid is a unique identifier for each hex, long and lat are longitude and latitude of the centroid of each hex, and occur is a count of how many years there was activity in each hex."))

props <- properties(abdat)  # characterize each field in data.frame
noid <- dim(abdat)[1]
kable(props[,1:4],row.names=TRUE)          # tabulate the characteristics of abdat
label <- paste0(" A schematic map of the coastline of Tasmania (black line) combined with the locations of the centroids of each of the ", noid, " 1-hectare hexagonal grids (red dots) that reported blacklip effort and catch data sometime between 2012 - 2018.")

\newline

There are r noid individual hexagonal grid cells defined around the Tasmanian coast from which blacklip abalone have been taken at least in some of the years 2012 - 2018. The location of each of these r noid unique hexagonal grid cells (or oids) is detailed by statistical zone, blockno, and subblockno. In addition there are also the distance from the closest shoreline and the latitude and longitude of each oid's centroid. Then, for each of the years 2012 - 2018 there is listed the total blacklip abalone catch (prefix blkg) in kg, the number of separate days in which it was fished, (prefix day), the number of minutes of effort summed for each year (prefix mins), the number of divers to visit in a given year (prefix divers), and finally the kg/day in each year (prefix kgday). Finally there is a count of how many years in which there was fishing activity for each hex or oid (occur).

To gain an understanding of the spread of blacklip data we can now plot up the location of each of the active 1-hectare grids. Where the coastline can be seen reflects areas which were not fished over that seven year period; generally this implies the coastline was unsuitable for abalone or was closed to all fishing.

leftlong <- 143.5;  rightlong <- 148.75 
uplat <- -39.0;  downlat <- -44.0 
alltas <- c(leftlong,rightlong,uplat,downlat)
maptas(view=alltas,gridon=1.0) 
mappoints(view=alltas,abdat) 

r fig_nums("F1",caption= label)

\newline

The same variables present in the blacklip data (r tab_nums("T2", display="num")) are present in the greenlip data, only there are fewer records.

# blacklip 1 Ha grid locations
filename <- "G1HaGL_2019_03_06.rds"
abG <- readRDS(paste0(datadir,filename))#get spatial obj
abGdat <- getlatlong(abG)     # convert to data.frame + lat long
abGdat$occur <- apply(abGdat[,5:11],1,countgtzero)
propsG <- suppressWarnings(properties(abGdat))  # characterize each field in data.frame

noidG <- dim(abGdat)[1]
label <- paste0(" A schematic map of the coastline of Tasmania (black line) combined with the locations of the centroids of each of the ", noidG, " 1-hectare hexagonal grids (red dots) that reported greenlip effort and catch data sometime between 2012 - 2018.")

\newline

maptas(view=alltas,gridon=1.0) 
mappoints(view=alltas,abGdat) 

r fig_nums("F2",caption = label)

\newline

The commercial catch and effort data from the logbooks (dockets) is in two datasets, one for the blacklip data and the other for the greenlip data. This data has the following properties.

filenameCE <- "CatchEffortData_2019-05-13.RData"
load(paste0(datadir,filenameCE))
numcebl <- dim(abCEbl)[1]
numcgbl <- dim(abCEgl)[1]

\newline

r tab_nums("T3", caption= paste0(" Yearly statistics for total catch and number of records for both blacklip and greenlip data from ",filenameV,". The four columns prefixed CE, relate to the docket CPUE data. The catches are comparable but not the number of records."))

kable(properties(abCEbl))

\newline

The greenlip data has the same fields except there is no plaindate. In each case they both have data from 1974 to 2019, but for the analyses we need to select for the years 2012 - 2018. These data are useful for a comparison to determine how much of the annual catch and effort is geo-referenced. We will need to determine whether the proportion coverage is essentially the same across the State of Tasmania.

\newline

picky <- which((abCEbl$fishyear > 2007) & (abCEbl$fishyear < 2019))
cbybl <- tapply(abCEbl$blips[picky],abCEbl$fishyear[picky],sum,na.rm=TRUE)/1000
rbybl <- tapply(abCEbl$blips[picky],abCEbl$fishyear[picky],countgtzero)
pickG <- which((abCEgl$fishyear > 2007) & (abCEgl$fishyear < 2019))
cbygl <- tapply(abCEgl$glips[pickG],abCEgl$fishyear[pickG],sum,na.rm=TRUE)/1000
rbygl <- tapply(abCEgl$glips[pickG],abCEgl$fishyear[pickG],countgtzero)
outcomeCE <- cbind(cbybl,rbybl,cbygl,rbygl)
colnames(outcomeCE) <- c("BL_C_CE","BL_R_CE","GL_C_CE","GL_R_CE")

\newline

# blacklip 1 Ha grid locations
filenameV <- "G1HaV_2019_05_13.rds"
abV <- readRDS(paste0(datadir,filenameV))#get spatial obj
abVdat <- getlatlong(abV)     # convert to data.frame + lat long
pickyr <- which(abVdat$fishyear > 2007)
abVdat <- droplevels(abVdat[pickyr,])
propsV <- properties(abVdat)  # characterize each field in data.frame
noidV <- dim(abVdat)[1]

tab3

r tab_nums("T4", caption=" The properties of the record by record GPS data relating to the hexagon data.")

kable(propsV)

maptas(view=alltas,gridon=1.0) 
mappoints(view=alltas,abVdat) 

Comparison with LogBook Data

The data from the GPS data loggers has been collected in addition to the usual CPUE data collection on the paper dockets (logbooks). In every year there are instances of the loggers failing through the batteries not being charged, from mistakenly not being turned on, and through hardware failure. However, by comparing the legally required docket totals with those obtained in association with valid gps data it is posible to compare the proportion of the total catches captured within the gps data sets.

\newline

r tab_nums("T5", caption= paste0(" Yearly statistics for total catch and number of records for both blacklip and greenlip data from ",filenameV,". The four columns prefixed CE, relate to the docket CPUE data. The catches are comparable but not the number of records."))

blcby <- tapply(abVdat$blcatch_kg,abVdat$fishyear,sum,na.rm=TRUE)/1000
blrby <- tapply(abVdat$blcatch_kg,abVdat$fishyear,countgtzero)
glcby <- tapply(abVdat$glcatch_kg,abVdat$fishyear,sum,na.rm=TRUE)/1000
glrby <- tapply(abVdat$glcatch_kg,abVdat$fishyear,countgtzero)
outcome <- cbind(blcby,blrby,glcby,glrby)
colnames(outcome) <- c("BL_Catch","BL_Records","GL_Catch","GL_Records")
out <- cbind(outcome,outcomeCE)
kable(out,digits=c(3,0,3,0,3,0,3,0))

\newline

Data collection prior to 2012 was voluntary and exploratory so, while the early data has its uses, the main analyses of the fishery data will be restricted to the years 2012 - 2018. The mean proportion coverage across the years 2012 - 2018 was 86.4% for blacklip and 85.4% for greenlip.

\newline

r tab_nums("T6", caption= " Proportion of blacklip and reported greenlip abalone catches reported while successfully using the GPS data loggers. Note the proportional coverage only became useable from 2012 onwards.")

outP <- cbind(out[,1]/out[,5],out[,3]/out[,7])
colnames(outP) <-  c("Proportion of Blacklip","Proportion of Greenlip")
kable(outP,digits=c(4,4))

Comparison of Catch-per-year by block for Greenlip

Both the catch-rate data from log-books, and the GPS logger data have information on greenlip catches by year and block. These can be compared.

pick <- which((abCEgl$fishyear > 2011) & (abCEgl$fishyear < 2019))
ceGL <- droplevels(abCEgl[pick,])
cbyb <- tapply(ceGL$glips,list(ceGL$blockno,ceGL$fishyear),sum,na.rm=TRUE)/1000

picky <- which(abVdat$fishyear > 2011)
abVdat <- droplevels(abVdat[picky,])
cbybH <- tapply(abVdat$glcatch_kg,list(abVdat$blockno,abVdat$fishyear),
                sum,na.rm=TRUE)/1000
sumH <- rowSums(cbybH,na.rm=TRUE)
cbybH <- cbybH[sumH > 0,]
nce <- as.numeric(rownames(cbyb))
ngp <- as.numeric(rownames(cbybH))
blk <- sort(unique(c(nce,ngp)))
outGC <- matrix(0,nrow=length(blk),ncol=3,dimnames=list(blk,c("CE","GPS","Prop")))
outGC[match(nce,blk),"CE"] <- rowSums(cbyb,na.rm=TRUE)
outGC[match(ngp,blk),"GPS"] <- rowSums(cbybH,na.rm=TRUE)
outGC[,"Prop"] <- outGC[,"GPS"]/outGC[,"CE"]

\newline

r tab_nums("T7", caption= " Comparison of blacklip and greenlip proportional representation.")

round(outGC,3)

l

Hexagons by Block

First we can consider the number of blacklip hexagons with catches reported in the different statistical reporting blocks.

hbb <- tapply(abdat$oid,abdat$blockno,countgtzero)
catch <- round(tapply(abdat$blkgtotal,abdat$blockno,sum,na.rm=TRUE)/1000,3)
effort <- round(tapply(abdat$minstotal,abdat$blockno,sum,na.rm=TRUE)/60,3)
pick10 <- which(abdat$minstotal > 10)
hbb10 <-  tapply(abdat$oid[pick10],abdat$blockno[pick10],countgtzero)
hbb <- as.data.frame(cbind(as.numeric(names(hbb)),hbb,hbb10,catch,effort))
colnames(hbb) <- c("block","Hexagons","Effgt10","catch","minutes")
hbbo <- hbb[order(hbb$catch,decreasing = TRUE),]

\newline

<<<<<<< HEAD \newline ======= Now plot the centroids of each oid

5671a399ddce3ee49a1626ae97070e16c278264e

r tab_nums("T8", caption= " Number of hexagons per statistical block reported compared to number of hexagons with greater than 10 minutes data over seven years. Total catch over 7 years is also reported. All data sorted by maximum number of records")

leftlong <- 143.7;  rightlong <- 144.2 
uplat <- -39.5;  downlat <- -40.25 
KingIs <- c(leftlong,rightlong,uplat,downlat)
maptas(view=KingIs,gridon=0) 
pick <- which(ki$minstotal > 10)
mappoints(view=KingIs,ki,incex=0.5) 
mappoints(view=KingIs,ki[pick,],incol=4) 
abline(h=-39.885,col=3)
abline(v=143.95,col=3)

kable(halftable(hbbo,yearcol="block",subdiv=2))

<<<<<<< HEAD

\newline

par(mfrow=c(2,1),mai=c(0.45,0.45,0.05,0.05),oma=c(0.0,0,0.0,0.0)) 
par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)  
plot1(hbbo$Hexagons,hbbo$Effgt10,type="p",incol = 2,defpar = FALSE)
model <- lm(hbbo$Effgt10 ~ hbbo$Hexagons)
abline(model,col=4)
hist(hbbo$Effgt10,hbbo$catch,type="p",incol = 2,defpar = FALSE)

Now examine the King Island blocks 1 - 4, which report large numbers of hexagons.

pick <- which(abdat$blockno == 49)
ki <- droplevels(abdat[pick,])
tapply(ki$blkgtotal,ki$blockno,sum,na.rm=TRUE)/1000
tapply(ki$minstotal,ki$blockno,sum,na.rm=TRUE)/60
pick10 <- which(ki$minstotal > 20)
length(pick10)

\newline

Now plot the centroids of each oid

leftlong <- 143.7;  rightlong <- 144.2 
uplat <- -39.5;  downlat <- -40.25 
KingIs <- c(leftlong,rightlong,uplat,downlat)
maptas(view=KingIs,gridon=0) 
pick <- which(ki$minstotal > 20)
mappoints(view=KingIs,ki,incex=0.5) 
mappoints(view=KingIs,ki[pick,],incol=4) 
abline(h=-39.885,col=3)
abline(v=143.95,col=3)

=======

5671a399ddce3ee49a1626ae97070e16c278264e

pick1 <- which((ki$blockno == 1) &(ki$minstotal < 2))
top <- getmax(ki[pick1,"blkgtotal"])
par(mfrow=c(1,1),mai=c(0.45,0.45,0.05,0.05),oma=c(0.0,0,0.0,0.0)) 
par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)  
#hist(ki[pick1,"blkgtotal"],breaks=20,col=2,main="")
#hist(ki[pick1,"minstotal"],breaks=11,col=2,main="")
plot1(ki[pick1,"minstotal"],ki[pick1,"blkgtotal"],type="p",defpar=FALSE,
      xlabel="Minutes",ylabel="Catch Kg",inpch = 1,cex=1.1)
text(0,0.975*top,length(pick1),pos=4)
text(0,0.9*top,round(sum(ki[pick1,"blkgtotal"])/1000,3),pos=4)

sum(ki[pick1,"blkgtotal"])/1000
sum(ki[,"blkgtotal"])/1000

l

References

Mundy, C. (2011) Using GPS technology to improve fishery dependent data collection in abalone fisheries FRDC Project No: 2006/029, Fisheries Research and Development Corporation and University of Tasmania. 122p.

Mundy, C., Haddon, M., Dichmont, C.D., and B. Venables (2018) Can spatial fishery-dependent data be used to determine Stock Status in a spatially structured fishery? FRDC Funding Application, 2017-026. University of Tasmania. 11p.

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