Methodology

Species

Three seabird species were included in this analysis: * Yelkouan shearwater (Puffinus yelkouan) * Scopoli’s shearwater (Calonectris diomedea) * Audouin’s gull (Larus audouinii)

Geographical scope

id	iso3	label
5673	HRV	Croatia
50167	HRV	Joint regime area Croatia / Slovenia

The Croatian exclusive economic zone (EEZ code 5673) was included in this analysis. The joint regime area with Slovenia (EEZ code 50167) was not.

Input data

Fishing effort

This analysis used both VMS data provided by the national competent authority and the AIS “Apparent Fishing Effort” dataset provided by Global Fishing Watch https://globalfishingwatch.org/.

VMS (Vessel Monitoring System)

The VMS data provided for this analysis was limited in several important regards, specifically: - includes only vessels using drifting longlines - the spatial extent of the data is limited (see map) - provided VMS data included all vessel movement, with no explicit way to filter only locations where fishing activity was performed.

vms <- read_delim("data-input/VMS brodice.csv", delim = ";") %>% 
  filter(!grepl("^[0-9\\-]*$", GPS_DTIME)) %>% 
  mutate(dttm = parse_date_time(GPS_DTIME, orders = "%Y-%m-%d %H:%M", exact = TRUE, truncated = 1)) 

vms_hull <- st_read("data-input/vms_brodice_hull.gpkg")

vms %>%
  st_as_sf(coords = c("WGS84_LONG", "WGS84_LAT"), crs = wgs84) %>% 
  st_geometry() %>% 
  terra::plot(pch = ".", col = "red")
lines(vms_hull, col = "blue")
lines(map_base_wgs)

In order to filter out location points that likely don’t correspond to fishing activity (i.e. vessel movement, anchoring, mooring), a coarse speed filter was applied, removing points with a recorded speed below 1 knot and above 10 knots. Following this, the points were rasterised by counting location points per grid cell. Finally, the raster values were clamped to the 99.5 percentile value to remove extreme outliers that tend to appear in ports.

vms %>%
  mutate(year = year(dttm)) %>% 
  group_by(year) %>% 
  summarise(min_date = min(dttm),
            max_date = max(dttm))

vms <- vms %>%
  mutate(month = month(dttm)) %>% 
  # filter(month %in% c(5,6,7,8) %>% 
  filter(KNOTS > 1 & KNOTS < 10) %>% 
  transmute(lat = round(WGS84_LAT, digits = 2),
            lon = round(WGS84_LONG, digits = 2)) %>%
  group_by(lon, lat) %>% 
  summarise(effort = n(), .groups = "keep")

cutoff <- vms$effort %>% quantile(.995)
vms$effort_cl <- clamp(vms$effort, upper = cutoff, values = TRUE)

vms_rast <- rast(vms, crs = wgs84_str)

AIS (Automatic Identification System)

The AIS-based dataset was accessed through the GFW API by using the gfwr package.

feff <- readRDS("data-input/gfw_vessels_2015-2021.rds") %>% 
  dplyr::select(Lon, Lat, Time = `Time Range`, flag = Flag, geartype = `Gear Type`, fishing = `Apparent Fishing Hours`)

This AIS dataset includes fishing gear type and the apparent fishing effort in hours per grid cell. For more information on the methodology used by GFW to transform raw AIS data into the fishing effort dataset, please see: https://globalfishingwatch.org/dataset-and-code-fishing-effort/

AIS data has a known bias for larger vessels. Cumulative effort by fishing gear shows trawlers as the most represented in the dataset, followed by ‘other purse-seines’ but including other gear types with a potential seabird bycatch risk, including drifting and set longlines as well as set gillnets:

feff %>% group_by(geartype) %>% summarise(effort_hrs = sum(fishing)) %>% arrange(desc(effort_hrs)) %>% 
  knitr::kable()

geartype	effort_hrs
trawlers	1.855673e+06
other_purse_seines	5.403368e+05
fishing	3.634301e+04
drifting_longlines	3.541763e+04
set_longlines	2.638757e+04
tuna_purse_seines	2.340727e+04
set_gillnets	2.157518e+04
fixed_gear	1.048694e+04
pole_and_line	4.358924e+03
pots_and_traps	7.202127e+02
NULL	4.939928e+02
purse_seines	3.043678e+02
dredge_fishing	9.129000e+01
trollers	2.340278e+00

Filtering for set and drifting longlines and grouping by vessel flag reveals that for Croatian vessels the GFW dataset has very little data for set longlines and no data for drifting longlines:

feff_ll <- feff %>% filter(grepl("longline", geartype))
feff_ll %>% group_by(flag, geartype) %>% summarise(effort_hrs = sum(fishing)) %>% 
  arrange(geartype, desc(effort_hrs)) %>% 
  knitr::kable()

## `summarise()` has grouped output by 'flag'. You can override using the
## `.groups` argument.

flag	geartype	effort_hrs
ITA	drifting_longlines	34574.4847
CYP	drifting_longlines	843.1492
ITA	set_longlines	25171.2769
HRV	set_longlines	1216.2886

By year and month, the distribution of fishing effort in the dataset was the following:

feff %>% 
  mutate(month = gsub(".*-(\\d\\d)$", "\\1", Time),
         year = gsub("^(\\d\\d\\d\\d)-.*$", "\\1", Time)) %>% 
  ggplot(aes(y = reorder(month, desc(month)), x = fishing, fill = year)) + 
  geom_bar(stat = "identity") +
  ylab("month") +
  xlab("fishing effort in hours")

Cumulative fishing effort in hours by month/year

The years 2016-2021 were included in the analysis, and the estimated cumulative fishing effort was similar for all years:

feff %>% 
  mutate(year = gsub("^(\\d\\d\\d\\d)-.*$", "\\1", Time)) %>% 
  group_by(year) %>% 
  summarise(effort_hrs = sum(fishing)) %>% 
  knitr::kable()

year	effort_hrs
2015	396046.0
2016	387977.4
2017	388305.9
2018	338966.0
2019	350949.4
2020	325870.2
2021	367483.1

The cumulative fishing effort was calculated for each grid cell. For plotting purposes, the top 0.5% of values were restricted (or ‘clamped’) to eliminate outlier cells:

x <- feff %>% 
  group_by(Lon, Lat) %>% 
  summarise(sum = sum(fishing)) %>% 
  rast(crs = wgs84_str)
xtab <- feff %>% 
  group_by(Lon, Lat) %>% 
  summarise(sum = sum(fishing))
ll <- feff_ll %>% 
  group_by(Lon, Lat) %>% 
  summarise(sum = sum(fishing)) %>% 
  rast(crs = wgs84_str)
lltab <- feff_ll %>% 
  group_by(Lon, Lat) %>% 
  summarise(sum = sum(fishing))

map_extent <- ext(x)

cutoff <- xtab$sum %>% quantile(.995)
x_cl <- clamp(x, upper = cutoff, values = TRUE)
terra::plot(x_cl, col = hcl.colors(50, "Inferno"))
lines(vect(map_base_wgs))

Fishing effort across the study area

The same was done for longline fishing effort, but the results are almost exclusively outside Croatian territorial waters, which seems to further confirm that Croatian vessels are excluded from the sample.

cutoff <- lltab$sum %>% quantile(.995)
ll_cl <- clamp(ll, upper = cutoff, values = TRUE)
terra::plot(ll_cl, col = hcl.colors(50, "Inferno"))
lines(vect(map_base_wgs))

Longline fishing effort across the study area

Seabird distribution

Movement data was collected with GPS tracking devices mounted on seabirds from colonies around the island of Lastovo during the breeding seasons in 2019, 2020 and 2021. This data was analysed according to the BirdLife track2KBA methodology to identify core areas for each species based on individual birds’ overlapping utilization distributions (UD). Only representative core areas, as defined by the track2KBA methodology, were considered for this analysis. These areas represent areas used regularly by a significant number of birds from the tagged colonies. For a detailed description of the methodology please see https://github.com/BirdLifeInternational/track2kba.

core_yelk <- st_read("data-input/new-iba-areas/CoreAreas_Puffinus yelkouan.shp")
core_diom <- st_read("data-input/new-iba-areas/CoreAreas_Calonectris diomedea.shp")
core_audo <- st_read("data-input/new-iba-areas/CoreAreas_Larus audouinii.shp")

densityCalc <- function(tbl) {
  tbl <- tbl %>%
    filter(ptntlSt == 1) %>%
    dplyr::select(mature_indiv = Bst_MtI) %>%
    mutate(area_km2 = units::set_units(st_area(.), km^2),
           density = mature_indiv / area_km2) %>%
    arrange(desc(density)) %>%
    mutate(rel_density = density / density[2],
           rel_indiv = mature_indiv / max(mature_indiv))
  tbl$rel_density[1] <- tbl$rel_density[2]
  return(tbl)
}

cy <- rasterize(densityCalc(core_yelk), rast(x), field = "rel_indiv")
cd <- rasterize(densityCalc(core_diom), rast(x), field = "rel_indiv")
ca <- rasterize(densityCalc(core_audo), rast(x), field = "rel_indiv")

An output of this analysis is the estimated number of mature individuals regularly occupying a certain area. The inputs were resampled to the same raster as the GFW fishing data for further analysis.

plot(cy, ext = s_zoom_extent)
lines(vect(map_base_wgs))

Relative bird density for Puffinus yelkouan (closer view)

plot(cd, ext = s_zoom_extent)
lines(vect(map_base_wgs))

Relative bird density for Calonectris diomedea (closer view)

plot(ca, ext = s_zoom_extent)
lines(vect(map_base_wgs))

Relative bird density for Larus audouinii (closer view)

Output data

VMS vs AIS comparison

Because the GFW AIS-based fishing effort dataset doesn’t seem to include drifting longline effort for any Croatian-flag vessels, and the VMS dataset covers only a small part of the Croatian Adriatic, their direct comparison is challenging. We attempted to overlap this data and check to what extent they are similar, primarily to verify whether the AIS dataset, which covers the whole Croatian Adriatic, can reasonably be assumed to represent the fishing effort by drifting longlines. For this purpose, a comparison was performed in the following way:

The AIS data was masked to the area with rasterised VMS data
The values representing effort in both rasters (estimated hours in AIS, number of overlapping vessel location points in VMS) were normalised to the range 0-1
For pairwise combinations of different AIS dataset gear types and the VMS raster, RMSE (root mean square error) was calculated as a measure of difference between the rasters

raster.rmse <- function(rast_a, rast_b) sqrt(mean(values((rast_a - rast_b)^2), na.rm = TRUE))

get_geartype <- function(feff, i) {
  x <- feff %>%
    filter(geartype == i) %>% 
    group_by(Lon, Lat) %>% 
    summarise(sum = sum(fishing), .groups = "keep") %>% 
    rast(crs = wgs84_str)
  
  xtab <- feff %>% 
    filter(geartype == i) %>% 
    group_by(Lon, Lat) %>% 
    summarise(sum = sum(fishing), .groups = "keep")
  
  cutoff <- xtab$sum %>% quantile(.995)
  x_cl <- clamp(x, upper = cutoff, values = TRUE)

  return(x_cl)
}

compare_with_vms <- function(x_cl, vms_rast) {
  x_cl <- crop(x_cl, ext(vms_rast)) %>%
    subst(NA, 0)
  x_cl_masked <- mask(x_cl, vms_rast)
  x_cl_masked$norm <- values(x_cl_masked$lyr1) / max(values(x_cl_masked$lyr1), na.rm = TRUE)
  return(x_cl_masked$norm)
}

x_fishing <- get_geartype(feff, "fishing")
x_trawlers <- get_geartype(feff, "trawlers")
x_other_purse <- get_geartype(feff, "other_purse_seines")
x_tuna_purse <- get_geartype(feff, "tuna_purse_seines")
x_set_longlines <- get_geartype(feff, "set_longlines")
x_all <- get_geartype(feff, feff$geartype)

vms_rast$norm <- values(vms_rast$effort_cl) / max(values(vms_rast$effort_cl), na.rm = TRUE)

par(mfrow = c(1,2))
plot(compare_with_vms(x_fishing, vms_rast), main = paste0("AIS data: fishing"))
lines(vect(map_base_wgs))

plot(vms_rast$norm, main = "VMS data: drifting longlines")
lines(vect(map_base_wgs))

mtext(text = paste0("difference (RMSE): ", raster.rmse(compare_with_vms(x_fishing, vms_rast), vms_rast$norm)), side = 1, outer = FALSE)

Comparison between AIS ‘fishing’ category and VMS drifting longlines raster

par(mfrow = c(1,2))
plot(compare_with_vms(x_trawlers, vms_rast), main = paste0("AIS data: trawlers"))
lines(vect(map_base_wgs))

plot(vms_rast$norm, main = "VMS data: drifting longlines")
lines(vect(map_base_wgs))

mtext(text = paste0("difference (RMSE): ", raster.rmse(compare_with_vms(x_trawlers, vms_rast), vms_rast$norm)), side = 1, outer = FALSE)

Comparison between AIS ‘trawlers’ category and VMS drifting longlines raster

par(mfrow = c(1,2))
plot(compare_with_vms(x_other_purse, vms_rast), main = paste0("AIS data: other purse seines"))
lines(vect(map_base_wgs))

plot(vms_rast$norm, main = "VMS data: drifting longlines")
lines(vect(map_base_wgs))

mtext(text = paste0("difference (RMSE): ", raster.rmse(compare_with_vms(x_other_purse, vms_rast), vms_rast$norm)), side = 1, outer = FALSE)

Comparison between AIS ‘other purse seines’ category and VMS drifting longlines raster

par(mfrow = c(1,2))
plot(compare_with_vms(x_all, vms_rast), main = paste0("AIS data: all fishing gear types"))
lines(vect(map_base_wgs))

plot(vms_rast$norm, main = "VMS data: drifting longlines")
lines(vect(map_base_wgs))

mtext(text = paste0("difference (RMSE): ", raster.rmse(compare_with_vms(x_all, vms_rast), vms_rast$norm)), side = 1, outer = FALSE)

Comparison between AIS total fishing effort category and VMS drifting longlines raster

Based on this analysis, the most similar fishing gear type in the AIS dataset seems to be “other purse seines” (RMSE: 0.212), followed by “trawlers” (RMSE: 0.225). The normalised sum of all fishing gear types seems to have the lowest RMSE (0.2115) of the fishing effort rasters. The overall fishing patterns in this area seem to be broadly similar across the entire range of different fishing gear types.

Risk map

In order to combine the data into a map of elevated seabird bycatch risk, a product of the fishing effort and bird density rasters was calculated for each bird species. This effectively means that only the representative core areas will be assigned a risk value, with the relative bird density acting as a weighting factor – e.g. for a grid cell with a bird density of 0.5, the fishing effort would need to be twice as large to have the same risk value as a grid cell with a bird density of 1 (the highest value). The resulting rasters were downsampled with bilinear interpolation to create more manageable areas and to smooth out the outlier grid cells.

miba <- st_read("data-output/miba_novi_2.gpkg")

## Reading layer `miba_novi_2' from data source 
##   `C:\Users\mzec\Documents\Code\artina-telemetry\data-output\miba_novi_2.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 8 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 12.89215 ymin: 42.37587 xmax: 17.81552 ymax: 45.09187
## Geodetic CRS:  WGS 84

risk_yelk <- (cy * x)
risk_diom <- (cd * x)
risk_audo <- (ca * x)

resampleRaster30 <- function(spatRast){resample(spatRast, 
                                                rast(resolution = c(0.03, 0.03), 
                                                     extent = ext(x), 
                                                     crs = wgs84_str), 
                                                method = "bilinear")}

effort_30 <- resampleRaster30(x)
risk_yelk_30 <- resampleRaster30(risk_yelk)
risk_diom_30 <- resampleRaster30(risk_diom)
risk_audo_30 <- resampleRaster30(risk_audo)

par(mfrow = c(1,2))
plot(risk_yelk, col = hcl.colors(30, "viridis"), ext = s_zoom_extent)
lines(vect(map_base_wgs))
plot(risk_yelk_30, col = hcl.colors(30, "viridis"), ext = s_zoom_extent)
lines(vect(map_base_wgs))

Puffinus yelkouan estimated bycatch risk in the South Adriatic (1/100 vs 1/30 degree resolution)

Results

The North Adriatic seems to be an important foraging site for Yelkouan shearwaters, with some possible pressure from fishing. As this part of the Adriatic is very shallow (up to ~40 m depth), trawling is the main fishing method here, which poses unclear (and likely low) seabird bycatch risk.

plot(effort_30, ext = n_zoom_extent)
lines(vect(map_base_wgs))
lines(vect(miba), col = "purple")

Estimated total fishing effort in the North Adriatic (1/30 degree resolution)

plot(risk_yelk_30, col = hcl.colors(30, "viridis"), ext = n_zoom_extent)
lines(vect(map_base_wgs))
lines(vect(miba), col = "purple")

Puffinus yelkouan estimated bycatch risk in the North Adriatic (1/30 degree resolution)

In the South Adriatic, the bycatch risk seems to be higher in the west part of the Lastovo channel, along the western tip of the island of Korčula, the western part of the Hvar channel and along the southern coast of the island of Šolta. As the Lastovo channel in close proximity to the Yelkouan shearwater colonies, some of the bird density might not be foraging birds, but rather rafting, or simply birds arriving to and departing from the colony. Therefore the bycatch risk in this place might actually be lower than it appears from the map.

plot(effort_30, ext = s_zoom_extent)
lines(vect(map_base_wgs))
lines(vect(miba), col = "purple")

Estimated total fishing effort in the South Adriatic (1/30 degree resolution)

plot(risk_yelk_30, col = hcl.colors(30, "viridis"), ext = s_zoom_extent)
lines(vect(map_base_wgs))
lines(vect(miba), col = "purple")

Puffinus yelkouan estimated bycatch risk in the South Adriatic (1/30 degree resolution)

Scopoli’s shearwaters show the bycatch risk concentrated in a more narrow area, with the western part of the Lastovo channel still showing up as an area of elevated risk. The caveat described above for Yelkouan shearwaters still applies here.

plot(risk_diom_30, col = hcl.colors(30, "viridis"), ext = s_zoom_extent)
lines(vect(map_base_wgs))
lines(vect(miba), col = "purple")

Calonectris diomedea estimated bycatch risk in the South Adriatic (1/30 degree resolution)

Finally, because their movement pattern is much more dispersed and their core areas more restricted, the bycatch risk results for Audouin’s gulls show similarly restricted areas around the Lastovnjaci and Vrhovnjaci arhipelagos and the western tips of the Pelješac peninsula and the island of Mljet.

plot(risk_audo_30, col = hcl.colors(30, "viridis"), ext = s_zoom_extent)
lines(vect(map_base_wgs))
lines(vect(miba), col = "purple")

Larus audouinii estimated bycatch risk in the South Adriatic (1/30 degree resolution)

In order to compare the overall fishing effort inside vs. outside of mIBAs, a number of measures were used, including mean and median values for grid cells, as well as the sum of fishing effort divided by the total area inside vs outside of mIBAs. The mean, median and average by area values are presented in the table below. Overall, all measures seem to suggest that fishing effort is indeed higher inside mIBAs than outside them.

miba_union <- summarise(miba)
miba_union$area <- st_area(miba_union)

marineregion_cro <- st_read("data-input/marineregion_CRO.shp")

## Reading layer `marineregion_CRO' from data source 
##   `C:\Users\mzec\Documents\Code\artina-telemetry\data-input\marineregion_CRO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 34 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 13.00833 ymin: 41.6297 xmax: 18.54925 ymax: 45.56495
## Geodetic CRS:  WGS 84

marineregion_cro_outside_iba <- st_difference(marineregion_cro, miba_union) %>% 
  dplyr::select(MRGID_EEZ)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

marineregion_cro_outside_iba$area <- st_area(marineregion_cro_outside_iba)

geartypes_minus_trawlers <- c("fishing", "set_gillnets", "other_purse_seines", 
                              "tuna_purse_seines", "pole_and_line", "set_longlines", "dredge_fishing", 
                              "fixed_gear")

x <- feff %>%
  filter(geartype %in% geartypes_minus_trawlers) %>% 
  group_by(Lon, Lat) %>% 
  summarise(sum = sum(fishing), .groups = "keep") %>% 
  rast(crs = wgs84_str)

xtab <- feff %>% 
  filter(geartype %in% geartypes_minus_trawlers) %>% 
  group_by(Lon, Lat) %>% 
  summarise(sum = sum(fishing), .groups = "keep")

cutoff <- xtab$sum %>% quantile(.995)
x_cl <- clamp(x, upper = cutoff, values = TRUE)

x_cl_outside_miba <- mask(x_cl, miba, inverse = TRUE)
x_cl_inside_miba <- mask(x_cl, miba, inverse = FALSE)

effort_in_miba <- values(x_cl_inside_miba) %>% as_tibble() %>% transmute(effort = sum, miba = "inside miba")
effort_out_miba <- values(x_cl_outside_miba) %>% as_tibble() %>% transmute(effort = sum, miba = "outside miba")

effort_comp <- tibble(measure = c("mean fishing effort by grid cell", 
                                  "median fishing effort by grid cell", 
                                  "average fishing effort per km^2"),
                      `inside mIBA` = c(mean(values(x_cl_inside_miba), na.rm = TRUE),
                                        median(values(x_cl_inside_miba), na.rm = TRUE),
                                        sum(values(x_cl_inside_miba), na.rm = TRUE) / (as.numeric(miba_union$area) / 1000000)),
                      `outside mIBA` = c(mean(values(x_cl_outside_miba), na.rm = TRUE),
                                         median(values(x_cl_outside_miba), na.rm = TRUE),
                                         sum(values(x_cl_outside_miba), na.rm = TRUE) / (as.numeric(marineregion_cro_outside_iba$area) / 1000000)))

knitr::kable(effort_comp)

measure	inside mIBA	outside mIBA
mean fishing effort by grid cell	19.913044	17.257806
median fishing effort by grid cell	7.758472	5.735138
average fishing effort per km^2	13.407151	11.521646

Assumptions and limitations

This approach includes some limitations which are important to keep in mind when interpreting the results. Inside the area where data was available for all fishing gear types (VMS for longline data and AIS for everything else), the patterns in fishing effort seem to be broadly similar, which suggests that AIS fishing effort data can generally be thought of as a reasonable measure of overall fishing effort, even given its inherent limitations (most prominently the fact that it is biased towards larger vessels).

Because the exact connection between bycatch risk and different types of fishing gear is still not fully understood, we made no distinction between the different kinds of fishing gear as detected/registered in the AIS dataset for the overall bycatch risk calculation.

Finally, no distinction was made based on the conservation status of the species (e.g. the global IUCN red list status of “VU-Vulnerable” for Puffinus yelkouan and “LC-Least Concern” for Calonectris diomedea). An arbitrary choice of weighing factor could be made to reflect this difference in the final risk map.

Conclusions

Fishing effort from the GFW AIS-based dataset follows the same broad activity patterns as VMS data for drifting longlines. In the absence of more comprehensive VMS data for an array of fishing activities, the AIS-based fishing effort is a reasonable approximation of total fishing effort.
For the channel IBAs (Lastovo, Korčula, Hvar), the western parts of the IBAs seem to be under more pressure from fishing activity than the eastern parts.
Overall, the fishing effort within marine IBAs seems to be on average higher than outside them (i.e. in non-mIBA areas of the Croatian Adriatic) when taking total area into account.

A considerable spatial overlap exists between high fishing effort areas as defined by AIS data and areas actively used by seabirds during breeding season, and a risk of seabird bycatch in these areas cannot be ruled out. Further refinements should be based on a more detailed understanding of the risk that different fishing gear types pose to seabirds, and if possible include both AIS and VMS data.

Seabird bycatch risk analysis v4

Mate Zec, Association Biom

2023-05-05

Background