mapping_script.R

#################################
##                             ##
##  EU-Norway VMS & Logbook    ##
##   Cod Data Processing       ##
##                             ##
#################################

## This script takes the standardised EFLALO and TACSAT files which are held
## at national level, and processes them to produce an aggregated data product
## which can be combined with those generated by other countries to visualise
## the distribution of cod landings from the North Sea.

## The assumed input data to the script are the processed and quality-assured 
## files used to fulfill the ICES VMS data call


#- Clear workspace
rm(list=ls())

## set location where your EFLALO and TACSAT files are
folder.path <- "C:/Work/VMS/Results/"

## required libraries
library(raster)
library(vmstools) #- download from www.vmstools.org
library(sp)

country <- "UK"

years <- c(2014, 2015, 2016, 2017, 2018)
## we are more interested in recent years, so lets see if we can use 2014 onwards
mesh.sizes <- c("below_100mm", "100mm_plus")
## and we are going to have two mesh size categories, above and below 100mm

data(europa)
pdf(file = paste(folder.path, "Cod_Maps.pdf", sep = ""))


for(i in 1:length(years)){
for(j in 1:length(mesh.sizes)){
  
load(paste(folder.path, "cleanEflalo", years[i], ".Rdata", sep = ""))
load(paste(folder.path, "cleanTacsat", years[i], ".Rdata", sep = ""))

# in our data we have
## "OTB" "OTT" "GNS" "FPO" "OTM" "PTM" "TBN" "PTB" "LX"  "TB"  "GNC"
## "GN"  "LHP" "GTR" "GND" "LHM" "DRB" "SSC" "SDN" "LL"  "TBB" "LLS"

# we want 
## "OTB", "OTT", "PTB", "SSC", "SDN"
## at mesh sizes >= 100mm for TR1, and <100mm for TR2

## we are only interested in trips which record a landing of cod, coming from the North Sea stock area (3a, 4 and 7d)
## we aren't interested in other species, or in what quantities...

## contstrain eflalo to be trips landing cod using selected gear types, in correct area
eflalo <- eflalo[eflalo$LE_KG_COD > 0 & eflalo$LE_GEAR %in% c("OTB", "OTT", "PTB", "SSC", "SDN") & eflalo$LE_DIV %in% c("27.3.a", "27.4.a", "27.4.b", "27.4.c", "27.7.d"),]

if(j == 1){eflalo <- eflalo[eflalo$LE_MSZ < 100,]}
if(j == 2){eflalo <- eflalo[eflalo$LE_MSZ >= 100,]}

## create a unique ID field for each vessel and date, for the eflalo and tacsat data
tacsat$UID <- paste(tacsat$VE_REF, tacsat$SI_DATE, sep="-")
eflalo$UID <- paste(eflalo$VE_REF, eflalo$LE_CDAT, sep="-")


## assigns fishing trip reference and gear values to tacsat data, where the unique ID field matches
tacsat$FT_REF  <- paste(eflalo$FT_REF)[match(tacsat$UID, eflalo$UID)]
tacsat$LE_GEAR <- paste(eflalo$LE_GEAR)[match(tacsat$UID, eflalo$UID)]
tacsat$LE_MSZ  <- paste(eflalo$LE_MSZ)[match(tacsat$UID, eflalo$UID)]


## creates a null field representing fishing activity in the tacsat data
## (we are treating multiple gear types in a single step here, which is why I did it this way)
tacsat$ACTIVITY <- 0

 ## assigns a value of 1 to pings at fishing speeds - the range used is 3-5 knots for otter trawls and 0.5 to 5 for seines
## other fleets may have their own more appropriate ranges 
tacsat$ACTIVITY[tacsat$LE_GEAR %in% c("OTB", "OTT", "PTB") & tacsat$SI_SP >=3 & tacsat$SI_SP <=5] <- 1
tacsat$ACTIVITY[tacsat$LE_GEAR %in% c("SSC", "SDN") & tacsat$SI_SP >=0.5 & tacsat$SI_SP <=5] <- 1

## removes tacsat data which have not been correctly assigned to a fishing trip record
tacsat <- tacsat[!is.na(tacsat$FT_REF),]

## removes tacsat data which is not at fishing speeds
tacsat <- tacsat[tacsat$ACTIVITY == 1,]

## count the number of pings at fishing speeds per each vessel 
day.pings <- tapply(tacsat$ACTIVITY, tacsat$UID, sum)

sum(day.pings == 0) ## you shouldn't have an zeroes

# aggregate cod catches to UID level first
eflalo_short<-eflalo[,c("UID", "LE_KG_COD")]
eflalo_agg<-aggregate(eflalo_short$LE_KG_COD, by=list(eflalo_short$UID), FUN=sum,na.rm=T)
names(eflalo_agg)<-c("UID", "LE_KG_COD")
    
## calculate the catch per ping for each vessel/day
#cod.per.ping <- eflalo$LE_KG_COD/day.pings[match(eflalo$UID,names(day.pings))]
cod.per.ping  <- eflalo_agg$LE_KG_COD/day.pings[match(eflalo_agg$UID,names(day.pings))]
  
tacsat$COD <- cod.per.ping[match(tacsat$UID, names(cod.per.ping))]
## assign a value of cod catch to each ping


sp.cod <- SpatialPointsDataFrame(data.frame(tacsat$SI_LONG, tacsat$SI_LATI), data.frame(tacsat$COD/1000))
## creates a spatial point data frame for cod landings associated with each ping

cod.grid <- raster(ncol = 128, nrow = 112, xmn = -4, xmx = 12, ymn = 48, ymx = 62)
## grids the data at approx. 7.5nm scale (minimum used so far for cod spatial closures)


cod.raster <- rasterize(sp.cod, cod.grid, field = "tacsat.COD.1000", fun= sum)
## sums the cod point data over the raster

plot(europa, ylim=c(48, 62), xlim=c(-4, 12), col="grey70", asp=1.5)
box()
plot(cod.raster, add = T)
title(main = years[i], sub = mesh.sizes[j])


writeRaster(cod.raster, paste(folder.path, "Cod_", country, "_", years[i], "_", mesh.sizes[j],sep=""), overwrite = TRUE, format = "raster")
## writes the output
}}

dev.off()
## closes the pdf file