Practicals7i - nielshintzen/vmstools GitHub Wiki

Practical 7: The ICES VMS and logbook datacall

Outputting data to ICES is an important deliverable for member states that need to adhere to the ICES datacall. Every year, around January, ICES sends out a request to supply VMS and logbook data. This data goes into a database, fishframe, and requires a certain format. The workflow below links many steps from P1 to P4 to get to an appropriate product for the ICES datacall.

The steps we will go through are: (please note that this is a suggested workflow, as the ICES call may vary from year to year, this is not a perfect template at all!)

  1. Cleaning the tacsat and eflalo data
  2. Merging eflalo to tacsat to get some info on gears used assigned to tacsat
  3. Define activity
  4. Assign euro's and kilo's from eflalo to tacsat
  5. Assign some lenght categories and metiers
  6. Assign some temporal and spatial references
  7. Aggregate data at the correct detail level in fishframe format

What is not part of the ICES datacall itself is calculating Swept Area Ratios while everything in the datacall is designed such that these values can be calculated afterwards within the ICES expert groups. Almost annually ICES produces maps of fishing pressure (in Swept Area Ratios, SAR) for clients and its own ecosystem overviews. In this practical we also show you how you can get there yourself.

Now onto processing the data

# 1. Cleaning the tacsat data
# We've done this before, so we go a bit quick over the next steps (see practical 2 for more info on cleaning)

data(tacsat)

#Remove points that are not possible, not unique, on land, in harbour or have an interval rate < 5 minutes
idx 	<- which(abs(tacsat$SI_LATI) > 90 | abs(tacsat$SI_LONG) > 180) 
idx 	<- unique(c(idx,which(tacsat$SI_HE < 0 | tacsat$SI_HE >  360))) 
tacsat$SI_DATIM <- as.POSIXct(paste(tacsat$SI_DATE,tacsat$SI_TIME,sep=" "), 
                              tz="GMT", format="%d/%m/%Y  %H:%M") #create one date-time stamp
                              
#Check whether there are any values in idx
length(idx)
#In this instance there are not any but if there were you would need to remove those rows, e.g. tacsat <- tacsat[-idx,]

uniqueTacsat    <- paste(tacsat$VE_REF,tacsat$SI_LATI,tacsat$SI_LONG,tacsat$SI_DATIM) 
tacsat          <- tacsat[!duplicated(uniqueTacsat),] #get rid of the duplicates
data(europa) #a dataset already build into the package with the coastlines embedded
                                       				
                                       				
idx     	<- pointOnLand(tacsat,europa,st_crs=4326);

pol     	<- tacsat[which(idx==1),] #points on land
tacsat  	<- tacsat[which(idx==0),] #tacsat not on land

data(harbours)
idx 	<- pointInHarbour(tacsat$SI_LONG,tacsat$SI_LATI,harbours,saveHarbourList=F)
pih 	<- tacsat[which(idx==1),]
tacsat 	<- tacsat[which(idx==0),]
tacsat 	<- sortTacsat(tacsat) #sort the data by vessel and time
tacsatp	<- intervalTacsat(tacsat,level="vessel",fill.na=T)
intThres<- 5
tacsat  <- tacsatp[which(tacsatp$INTV > intThres),]

# 1. Cleaning the eflalo data
data(eflalo)
eflalop           <- eflalo
eflalop$FT_DDATIM <- as.POSIXct(paste(eflalo$FT_DDAT,  eflalo$FT_DTIME,   sep=" "), tz="GMT", format="%d/%m/%Y  %H:%M") #turn date and time into one date-time stamp
eflalop$FT_LDATIM <- as.POSIXct(paste(eflalo$FT_LDAT,  eflalo$FT_LTIME,   sep=" "), tz="GMT", format="%d/%m/%Y  %H:%M") #turn date and time into one date-time stamp
idx               <- which(eflalop$FT_LDATIM >= eflalop$FT_DDATIM)
eflalo            <- eflalo[idx,]

The cleaning should be a breeze by now. Note that there may be more points to check for, but that is really for you to decide as you should hopefully be familiar with the data Now let's move on to linking VMS and Logbooks and defining activity

# 2. Merging eflalo to tacsat (see practical 3 for more info)
#- These TBB vessels all have a gear width of 24m
eflalo$LE_WIDTH   <- 24
tacsatp		  <- mergeEflalo2Tacsat(eflalo,tacsat)
tacsatp$LE_GEAR   <- eflalo$LE_GEAR[ match(tacsatp$FT_REF,eflalo$FT_REF)]
tacsatp$LE_MSZ     <- eflalo$LE_MSZ[  match(tacsatp$FT_REF,eflalo$FT_REF)]
tacsatp$VE_LEN    <- eflalo$VE_LEN[  match(tacsatp$FT_REF,eflalo$FT_REF)]
tacsatp$VE_KW     <- eflalo$VE_KW[   match(tacsatp$FT_REF,eflalo$FT_REF)]
tacsatp$LE_RECT   <- eflalo$LE_RECT[ match(tacsatp$FT_REF,eflalo$FT_REF)]
tacsatp$LE_MET    <- eflalo$LE_MET[  match(tacsatp$FT_REF,eflalo$FT_REF)]
tacsatp$LE_WIDTH  <- eflalo$LE_WIDTH[match(tacsatp$FT_REF,eflalo$FT_REF)]
tacsatp$VE_FLT    <- eflalo$VE_FLT[  match(tacsatp$FT_REF,eflalo$FT_REF)]
tacsatp$LE_CDAT   <- eflalo$LE_CDAT[ match(tacsatp$FT_REF,eflalo$FT_REF)]
tacsatp           <- subset(tacsatp,FT_REF != 0)

# 3. Define activity (here for simplicity reasons, we only select beamtrawling, see practical 3 for more info)
tacsatp		  <- subset(tacsatp,LE_GEAR == "TBB")
eflalo		  <- subset(eflalo,LE_GEAR == "TBB")


storeScheme               <- expand.grid(years = 1800:1801, months = 0,weeks = 0,
                                         analyse.by = unique(tacsatp[,"LE_GEAR"]))
storeScheme$peaks         <- NA
storeScheme$means         <- NA
storeScheme$fixPeaks      <- FALSE
storeScheme$sigma0        <- 0.911
storeScheme$means[which(storeScheme$analyse.by == "TBB")]  <- c("-11.5 -6 0 6 11.5")
acTa                      <- activityTacsat(tacsatp,units="year",analyse.by="LE_GEAR",
                                              storeScheme=storeScheme,plot=FALSE,level="all")
tacsatp$SI_STATE          <- acTa
tacsatp$SI_STATE[which(tacsatp$SI_STATE != "f")] <- 0
tacsatp$SI_STATE[which(tacsatp$SI_STATE == "f")] <- 1

As said before, we only analysed a subset of gears. We created our own run schedule for that. That is a powerful way of setting up the scenario's you wish to evaluate. Now we will dispatch the kilo's and euro's from our logbook data to the VMS pings

# 4. Adding up kilo's and euro's by row
idx			  <- kgeur(colnames(eflalo)) #We get the columns that contain either LE_KG or LE_EURO values
eflalo$LE_KGS		  <- rowSums(eflalo[,grep("LE_KG_",colnames(eflalo))],na.rm=T)
eflalo$LE_EUROS		  <- rowSums(eflalo[,grep("LE_EURO_",colnames(eflalo))],na.rm=T)

Assign length classes, metiers and spatial references to our records.

Now it's time to start to use the ICES vocabulary. Nice enough has ICES made a package available that helps us with this. Use install.packages("icesVocab", repos = "https://ices-tools-prod.r-universe.dev") to install this package that allows you access to all the icesVocab components.

Length classes

library(icesVocab)
library(dplyr)
# Get the values accepted in this vocabulary dataset
vlen_ices <- getCodeList("VesselLengthClass") ### Get DATSU Vocabulary list for selected dataset

# Filter the vessel length categories required  by  ICES VMS& Logbook datacall 
vlen_icesc =  vlen_ices%>%
  filter ( Key %in% c("VL0006", "VL0608", "VL0810", "VL1012", "VL1215" ,"VL1518", "VL1824" ,"VL2440" ,"VL40XX"))%>%
  select(Key)%>%
  arrange(Key)

#We use the 'cut' function here which assigns the absolute length value to each of the length categories
eflalo$LENGTHCAT     <- eflalo$VE_LEN %>% cut(breaks=c(0, 6, 8, 10, 12, 15, 18, 24, 40, 'inf' ), 
                                             right = FALSE    ,include.lowest = TRUE,
                                             labels =  vlen_icesc$Key)

Metiers

The other part that we need to add is the metiers. The VMS and logbook datacall asks for level 4, 5 and 6 metiers. Level 4 should be part of your data already and is stored in the 'LE_GEAR' column. We need to check however if all the gears we have are also part of the ICES vocabulary

metierL4         <-  getCodeList("GearTypeL4")
eflalo$LE_MET_4 <- eflalo$LE_GEAR
which(!eflalo$LE_MET_4 %in% metierL4$Key)

That is good news, all the gears in my dataset are accepted by the ICES database. Now onto metier level 5. This is not usually part of the logbook data and it requires you to make an assignment yourself. For our example it would be good to look at species composition first before we assign a specific 'target assemblage'. As I have selected only the TBB gears, I can narrow myself to looking into species already rather than the species * gear interaction.

totalLandings   <- as.data.frame(sort(colSums(eflalo[,grep("LE_EURO",colnames(eflalo))])))
colnames(totalLandings) <- c("EURO")
totalLandings$Species  <- substr(rownames(totalLandings),9,11)
totalLandings   <- subset(totalLandings,EURO > 0)

ggplot(data=totalLandings,aes(x=Species,y=EURO)) + geom_bar(stat="identity") + theme_bw()

In this case, we see sole, plaice, turbot, brill, dab and flounder landings, all being flatfish being targeted by the same fishery. We see some nephrops and some shrimp as well, for the sake of this example we'll assign these a different category.


idxCSH <- which(eflalo$LE_KG_CSH > 0)
idxNEP <- which(eflalo$LE_KG_NEP > 0 & eflalo$LE_KG_CSH == 0)
idxDEF <- which((eflalo$LE_KG_SOL > 0 | eflalo$LE_KG_PLE > 0 | eflalo$LE_KG_DAB > 0) & eflalo$LE_KG_CSH == 0 & eflalo$LE_KG_NEP == 0)

#- Let's also make a bin for those records not beloning to these three categories
idxREM <- which(!1:nrow(eflalo) %in% c(idxCSH,idxNEP,idxDEF))
REMlandings   <- as.data.frame(sort(colSums(eflalo[idxREM,grep("LE_EURO",colnames(eflalo))])))
colnames(REMlandings) <- c("EURO")
REMlandings$Species  <- substr(rownames(REMlandings),9,11)
REMlandings   <- subset(REMlandings,EURO > 0)

#- And let's plot the species composition. Apparently there has been some small pelagic fishery for smelt?
ggplot(data=REMlandings,aes(x=Species,y=EURO)) + geom_bar(stat="identity") + theme_bw() 

eflalo$LE_MET_5 <- NA
eflalo$LE_MET_5[idxCSH] <- "CRU"
eflalo$LE_MET_5[idxNEP] <- "MCD"
eflalo$LE_MET_5[idxDEF] <- "DEF"
eflalo$LE_MET_5[idxREM] <- "SPF" #this category is not very accurate, as there are other demersal fish in there as well

metierL5        <- getCodeList("TargetAssemblage")
which(!eflalo$LE_MET_5 %in% metierL5$Key)

That already required quite a bit of an analyses. If you need to do this routinely, it may help to standardize an approach to assign target assemblage to your eflalo records, e.g. based on gear information, mesh sizes and your knowledge about the fishery.

Finally we need to assign a metier level 6 category, which includes mesh size. The logbooks usually report the mesh size used, which is simply 1 number. For the metiers, we have a mesh size range. Let's get all options and then assign these accordingly.

metierL6        <- getCodeList("Metier6_FishingActivity")
tail(metierL6)

We can now see that the 'Key' column has gear, target species, a mesh size range and some other values in one string. Let's find some of our main metiers that are associated with TBB and DEF.

metierL6$Key[grep("TBB_DEF",metierL6$Key)]

So it turns out we have 15 different options to choose from! Note that there are metiers that overlap each other, e.g. a metier that allows mesh size >=40 mm would also apply to many of the other metiers e.g. the metiers which list 40-54mm as a valid range. As there are other metier definitions by area (e.g. North Sea is different from the Baltic), there are many overlapping metiers listed here. There is a full list available at the ICES website with area, gear and allowed mesh size ranges.

table(subset(eflalo,LE_MET_5 == "DEF")$LE_MSZ)

While we have mainly mesh size 80 in our dataset but also some smaller and bigger ones. Let's automate the selection of the correct metier here but we take ones that are relevant for the North Sea only.

#- Get metier meshranges
metierSplit <- lapply(as.list(metierL6$Key[grep("TBB_DEF",metierL6$Key)][c(7,9,14,6)]),strsplit,"_")
metierMeshRange <- lapply(metierSplit,function(x){x[[1]][3]})
metierMeshRange <- do.call(rbind,lapply(lapply(metierMeshRange,strsplit,"-"),function(x){x[[1]]}))
metierMeshRange <- gsub(">","", metierMeshRange)
metierMeshRange <- gsub("<","", metierMeshRange)
metierMeshRange <- gsub("=","", metierMeshRange)

#- Get gear and target codes
metierGear <- "TBB"
metierTarget <- "DEF"

#- Area codes
metierArea <- 4

#- Metier
metierCode <- metierL6$Key[grep("TBB_DEF",metierL6$Key)][c(7,9,14,6)]

#- merge back to 1 table
metiers <- as.data.frame(cbind(unlist(metierArea), unlist(metierGear), unlist(metierTarget), unlist(metierMeshRange), unlist(metierCode)), stringsAsFactors = FALSE)
colnames(metiers) <- c("Area", "Gear", "Target", "Min_MSZ", "Max_MSZ", "Metier")

#- Remove metier NA
metiers <- metiers[-c(176),] # altijd 176?

#- Convert to numeric
metiers$Min_MSZ <- as.numeric(metiers$Min_MSZ)
metiers$Max_MSZ <- as.numeric(metiers$Max_MSZ)

#- Assign metier L6
eflalo$LE_MET_6 <- NA
for(i in 1:nrow(metiers)){
  idx <- which(eflalo$LE_MET_4 == metiers$Gear[i] & eflalo$LE_MET_5 == metiers$Target[i] & eflalo$LE_MSZ >= metiers$Min_MSZ[i] & eflalo$LE_MSZ <= metiers$Max_MSZ[i])
  if(length(idx)>0)
    eflalo$LE_MET_6[idx] <- metiers$Metier[i]
}

Exercise 1

  1. Is there a difference in species composition and spatial location for the TBB level 6 metiers?
Show answer

  1. To get an answer to this question we can split up our data into the different metier types we have. table(eflalo$LE_MET_6 which are apparently 3 different metiers. We have one metier with very small mesh sizes, one with intermediate mesh sizes (the dominant metier here) and one with large meshes. Let's first look at the spatial distribution.
metiersL6 <- names(table(eflalo$LE_MET_6))
p1 <- plotTools(subset(eflalo,LE_MET_6==metiersL6[1]),level="ICESrectangle",xlim=c(-4,10),ylim=c(50,60),log=T,returnPlot=T)

p2 <- plotTools(subset(eflalo,LE_MET_6==metiersL6[2]),level="ICESrectangle",xlim=c(-4,10),ylim=c(50,60),log=T,returnPlot=T)

p3 <- plotTools(subset(eflalo,LE_MET_6==metiersL6[3]),level="ICESrectangle",xlim=c(-4,10),ylim=c(50,60),log=T,returnPlot=T)

require(gridExtra)
grid.arrange(p1$plot,p2$plot,p3$plot,ncol=3)

So at least it is pretty clear that the larger mesh size are further north, the smaller mesh sizes are around the Dutch coastline and the intermediate one is all around the (southern) North Sea.

Now let's look at the species composition

require(tidyr)
specCols     <- colnames(eflalo)[grep("LE_KG_",colnames(eflalo))]
specCompMet6 <- eflalo %>% 
                #analyse by metier group
                group_by(LE_MET_6) %>% 
                #sum each of the species columns
                summarize(across(all_of(specCols),\(x) sum(x,na.rm=T))) %>%  
                #ungroup to convert from wide to relational format
                ungroup() %>% pivot_longer(cols=specCols,names_to="Species",values_to="KG") %>% 
                #only take catches > 0 and group again to sort and take the top records
                filter(KG>0) %>% group_by(LE_MET_6)  %>%
                arrange(LE_MET_6,KG) %>% top_n(7)

ggplot(data=specCompMet6,aes(x=Species,y=KG)) + geom_bar(stat="identity") + theme_bw() + facet_wrap("LE_MET_6",scales="free")

There clearly is a difference. The bigger mesh size metier only seems to catch plaice while the intermediate mesh size has a clear mixture of sole and place. I was not expecting to see the species as listed for the smallest mesh size as this is usually a shrimp fishery! Something must be wrong with the data in this example dataset!

So far it's been quite a hassle, but we need this to assign a 'benthis metier' to our data as this will allow us to calculate swept area rations.

A final part of the ICES datacall is to associate the species caught and their value as registered in the logbooks with the VMS pings. We've done that before using the splitAmongPings function which we'll use here now again.

#We split eflalo into two sets: one with only those eflalo records that have VMS associated per trip and one with all remaining trips (without VMS, e.g. small vessels)
eflalo$FT_DDATIM 	<- as.POSIXct(paste(eflalo$FT_DDAT,  eflalo$FT_DTIME,   sep=" "), tz="GMT", format="%d/%m/%Y  %H:%M") #turn date and time into one date-time stamp
eflalo$FT_LDATIM 	<- as.POSIXct(paste(eflalo$FT_LDAT,  eflalo$FT_LTIME,   sep=" "), tz="GMT", format="%d/%m/%Y  %H:%M") #turn date and time into one date-time stamp
eflaloM			<- subset(eflalo, FT_REF %in% unique(tacsatp$FT_REF)) 
eflaloNM		<- subset(eflalo,!FT_REF %in% unique(tacsatp$FT_REF)) 
# Split the euro's and kg's over the pings
tacsatEflalo  		  <- splitAmongPings(tacsat=tacsatp,eflalo=eflaloM,variable="all",level=c("day","ICESrectangle","trip"),conserve=T)

tacsatEflalo$LE_MET <- eflalo$LE_MET_6[match(tacsatEflalo$FT_REF,eflaloM$FT_REF)]

Spatial references

Finally, ICES uses the CSquare notation which we add including some other spatial-temporal aspects of the data. We use these settings to aggregate at the right scale later on

# 6. Assign some temporal and spatial references
data(ICESareas)
#- Subset to ICES areas 4a-4c
inICESarea            <- ICESareas$Area_27[ICESarea(tacsatEflalo,ICESareas,st_crs=4326)]
table(inICESarea) #goot to see, all are in ICES area 4b and 4c already
tacsatEflalo$Csquare  <- CSquare(tacsatEflalo$SI_LONG,tacsatEflalo$SI_LATI,degrees=0.05)
tacsatEflalo$Year     <- year(tacsatEflalo$SI_DATIM)
tacsatEflalo$Quarter  <- quarter(tacsatEflalo$SI_DATIM)
tacsatEflalo$Month    <- month(tacsatEflalo$SI_DATIM)
tacsatEflalo$Area     <- inICESarea
tacsatEflalo$kwHour   <- tacsatEflalo$VE_KW * tacsatEflalo$INTV/60

eflaloNM$Year           <- year(eflaloNM$FT_LDATIM)
eflaloNM$Quarter        <- quarter(eflaloNM$FT_LDATIM)
eflaloNM$Month          <- month(eflaloNM$FT_LDATIM)
eflaloNM$Area           <- "ICESarea" #shortcut to limit converting time between ICES rectangles and ICES areas
eflaloNM$LE_MET		      <- eflaloNM$LE_MET_level6

#- Sometimes its relevant to get the effort per eflalo row, that's what happens below
eflaloNM$INTV           <- difftime(eflaloNM$FT_LDATIM,eflaloNM$FT_DDATIM,units="hours")
eflaloNM$dummy          <- 1
res                     <- aggregate(eflaloNM$dummy,by=list(eflaloNM$FT_REF),FUN=sum,na.rm=T)
colnames(res)           <- c("FT_REF","nrRecords")
eflaloNM                <- merge(eflaloNM,res,by="FT_REF")
eflaloNM$INTV           <- eflaloNM$INTV / eflaloNM$nrRecords
eflaloNM$kwHour         <- eflaloNM$VE_KW * eflaloNM$INTV

Calculating Swept Area Ratios

The ICES products showing the Swept Area Ratios can be found here: https://ices-library.figshare.com/articles/dataset/Data_for_OSPAR_request_on_the_production_of_spatial_data_layers_of_fishing_intensity_pressure/18601508

We need to get information on metier conversion from the level 6 definition to the benthis metier definition. All this info is stored at ICES TAF.

library(icesConnect)
#install.packages(c("sfdSAR", "icesVMS"), repos = 'https://ices-tools-prod.r-universe.dev')
library(icesVMS)
library(sfdSAR)
#- Preparations to get gear width information for countries submitting 2009-2020 data or for countries that do not hold gear width information
# metier level 6 -> benthis metiers (updated 2022)
metier_lookup <-
  get_metier_lookup() %>%
  tibble() %>%
  select(leMetLevel6, benthisMetiers)

# benthis metiers -> gear width model and coefficients
gear_widths_models <-
  get_benthis_parameters() %>%
  tibble() %>%
  select(-id)

# metier level 6 -> gear width model and coefficients
metier6_gear_model_lookup <-
  metier_lookup %>%
  filter(!is.na(benthisMetiers)) %>%
  left_join(
    gear_widths_models, by = c(benthisMetiers = "benthisMet")
  ) #%>%
  #rename(
  #  a = "firstFactor", b = "secondFactor"
  #)
metier6_gear_model_lookup$a <- metier6_gear_model_lookup$firstFactor
metier6_gear_model_lookup$b <- metier6_gear_model_lookup$secondFactor
metier6_gear_model_lookup   <- metier6_gear_model_lookup[,-grep("firstFactor",colnames(metier6_gear_model_lookup))]
metier6_gear_model_lookup   <- metier6_gear_model_lookup[,-grep("secondFactor",colnames(metier6_gear_model_lookup))]

#- Get only the relevant lookup table records
metier6_gear_model_lookup   <- subset(metier6_gear_model_lookup,leMetLevel6 %in% unique(tacsatEflalo$LE_MET))

Now that we have the conversion matrix, we have to merge it to our tacsatEflalo dataset

tacsatEflalo$contactModel <- merge(tacsatEflalo,metier6_gear_model_lookup,by.x="LE_MET",by.y="leMetLevel6",all.x=T)$contactModel

And it also lets us predict the swept area. tacsatEflalo$SweptArea <- predict_surface_contact(tacsatEflalo$contactModel,tacsatEflalo$INTV/60,tacsatEflalo$LE_WIDTH/1000,tacsatEflalo$SI_SP)

In this case, we had gear width information already available as these beam trawlers fish with a fixed gear width but if we have to consider e.g. otter trawls, we need to predict gear width.

tacsatEflalo$AverageGearWidth <- merge(tacsatEflalo,metier6_gear_model_lookup,by.x="LE_MET",by.y="leMetLevel6",all.x=T)$gearWidth*1000

To finally go from swept areas to ratios we need to grid our data first, aggregate the swept areas within a grid cell, and then divide the swept areas with the grid cell surface

xrange  <- range(pretty(range(tacsatEflalo$SI_LONG)))
yrange  <- range(pretty(range(tacsatEflalo$SI_LATI)))

grd     <- createGrid(xrange,yrange,resx=1/4,resy=1/8,type="GridDF")
st_crs(grd) <- 4326

grID     <- st_over(st_as_sf(tacsatEflalo,coords=c("SI_LONG","SI_LATI"),crs=st_crs(4326)),grd)

tacsatEflalo$grID <- grID
aggData  <- aggregate(tacsatEflalo$SweptArea,by=list(tacsatEflalo$grID),FUN=sum,na.rm=T)
colnames(aggData) <- c("grID","SweptArea")
grd$SA   <- NA
grd$SA[aggData$grID] <- aggData$SweptArea
grd$Area <- st_area(grd)
library(units)
grd$Area <- set_units(grd$Area,km^2)
grd$SAR  <- grd$SA / grd$Area

Exercise 2

  1. Find the latest ICES advice for OSPAR on spatial layers and make a map of the SAR for the year 2020
Show answer

  1. Navigate to ices.dk -> advice -> latest advice. Then search on the page for 'OSPAR' and look up the technical service with reference to spatial layers. Click the link and then a new window pops up which displays the actual advice. But, we're looking for data files. If you scroll down on the page (not the advice), you can find a link to the supplemented data files which can be found at https://doi.org/10.17895/ices.data.8294. If we go there, we can download the shapefiles. In that list, there is a 'total-2020.shp' file, let's read that one in and make a pretty map.
t2020 <- st_read("C:/Downloads/ICES.2021.OSPAR_production_of_spatial_fishing_pressure_data_layers/shapefiles/total-2020.shp")

#- Set some colour scheme and elimate empty cells
require(RColorBrewer)
colintens   <- brewer.pal(9,"YlOrRd")  #colour for fishing intensity
cols <- c(NA,colintens)[cut(t2020$sar,breaks=c(-1,0,0.05,0.1,0.2,0.5,1,2,5,10,500))]
isNA <- which(is.na(cols))


#- Let's plot the map of europa, with on top the gridcells with info and we style our own legend

data(europa)
ggplot() + geom_sf(data=europa,fill="darkgrey",colour=1) +
           geom_sf(data=st_geometry(t2020[-isNA,]),aes(fill=cols[-isNA]),colour=NA) + theme_bw() + xlim(-4,10) + ylim(50,60) + geom_sf(data=europa,fill="darkgrey",colour=1) + xlab("Longitude") + ylab("Latitude") + 
           guides(fill=guide_legend(title="Swept Area Ratio (SAR)")) + 
           scale_fill_manual(values = c(rev(colintens)),
                    labels=rev(c("0 <= 0.05","0.05 <= 0.1", "0.1 <= 0.2","0.2 <= 0.5","0.5 <= 1","1 <= 2","2 <= 5","5 <= 10","> 10")))
Add a custom footer
⚠️ **GitHub.com Fallback** ⚠️