---
title: "data_overview"
author: "Bram Parmentier"
output: html_document
---

# Data figures for Parmentier et al (2025): "Small fish biomass in the North Sea is far greater than previously estimated"

Objective:
- Overview of the fish data 
- Overview of depth and mud for sampling locations

load libraries
```{r, warning=FALSE}
library(xlsx);library(ggplot2)
library(sf);library(raster)
library(dplyr);library(png)
library(grid);library(scales)
library(latticeExtra)
```

load trawl and fish data
```{r}
#set path
path <- c("")

trawldata <- read.csv(file.path(path,"trawldata.csv"))
fishdata <- read.csv(file.path(path,"fishbiomass.csv"))

trawldata <- trawldata[trawldata$StationID %in%  fishdata$StationID,] #2 proposed locations are not fished
trawldata <- trawldata[trawldata$Lat_mid < 56,] #outside area of interest
trawldata <- trawldata[trawldata$StationID != "COSTRA2008#44",] #sorting fraction too low (1/11 is sorted)

fishdata <- fishdata[fishdata$StationID %in% trawldata$StationID,] #Exclude fish samples for which trawling information is unavailable

#Year, month, and week
Sys.setlocale("LC_TIME", "C") # to make sure month names are in English
trawldata$Date <- as.Date(trawldata$Date)
trawldata$Year <- as.numeric(format(trawldata$Date, '%Y'))
trawldata$Month <- format(trawldata$Date, "%b")
trawldata$Month_1 <- as.numeric(format(trawldata$Date, "%m"))
trawldata$Year_Month <- paste0(trawldata$Year, trawldata$Month)

#depth data
depth <- read.csv(file.path(path,"otherdata/depth_NorthSea.csv"))
depth <- depth[depth$depth < 1,]
depth$depth <- depth$depth *-1 
```

One of the Surveys was used to compare sand mining pits with a reference area. This Survey is split in two (sand pits and reference area).
```{r}
mining <- trawldata[trawldata$Survey == "SANDEXTR22",]
selected_stations <- mining$StationID[grepl("-[^-]*-P$", mining$StationID)]
selected_stations <- mining[mining$StationID %in% selected_stations  | mining$StationID == "ARCA22#M9J-01-P_b",]$StationID
trawldata[trawldata$StationID %in% selected_stations,]$Survey <- paste0(trawldata[trawldata$StationID %in% selected_stations,]$Survey, "P")
```

Create a more informative Survey name which includes year, month (of most samples), and number of samples
```{r}
Year_Month <- as.data.frame(table(trawldata$Survey, trawldata$Year_Month))
Year_Month <- Year_Month[Year_Month$Freq != 0,]
Year_Month_max <- aggregate(Freq ~ Var1, Year_Month, max)

trawldata$Survey2 <- NA
trawldata$numbersamples <- NA

#create name with month based on most samples
for (i in unique(trawldata$Survey)) {
print(i)
trawldata[trawldata$Survey == i,]$Survey2 <- as.character(Year_Month[Year_Month$Var1 == i & Year_Month$Freq == Year_Month_max[Year_Month_max$Var1 == i,]$Freq,]$Var2)

trawldata[trawldata$Survey == i,]$numbersamples <- nrow(trawldata[trawldata$Survey == i,])
}

trawldata$Survey2 <- paste0(trawldata$Survey2, sep= "_", trawldata$numbersamples)
trawldata[trawldata$Survey2 == "2022Feb_34",]$Survey2 <- paste0(trawldata[trawldata$Survey2 == "2022Feb_34",]$Survey2, "P")
```

#trawldata
```{r}
#order trawldata based on date
earliest_dates <- aggregate(Date ~ Survey2, data = trawldata, FUN = min)
Survey2_order <- earliest_dates[order(earliest_dates$Date), "Survey2"]
trawldata$Survey2 <- factor(trawldata$Survey2, levels = Survey2_order)

#Year
table(trawldata$Survey2, trawldata$Year)
plot(table(trawldata$Year), las=2, xlab="Year", ylab="Number of samples")

#Month
table(trawldata$Survey2, trawldata$Month)
plot(table(trawldata$Month_1), xlab="Year", ylab="Number of samples")

#DOY
trawldata$Year <- as.numeric(format(trawldata$Date, '%Y'))
trawldata$Month <- format(trawldata$Date, "%b")
trawldata$Month_1 <- as.numeric(format(trawldata$Date, "%m"))
trawldata$Year_Month <- paste0(trawldata$Year, trawldata$Month)
trawldata$DOY <- as.POSIXlt(trawldata$Date)$yday + 1 

# Adjust DOY to start from October 1st
trawldata$DOY_adjusted <- ifelse(trawldata$DOY >= 274,  
                                 trawldata$DOY - 273,   
                                 trawldata$DOY + 92)

# Fig S2
#pdf(file = file.path(path, "figures/supplement/S2_samples_during_year.pdf"), width = 170 / 25.4, height = 80 / 25.4)  # mm to inches
par(cex = 0.67)
plot(table(trawldata$DOY_adjusted), xlab="Day from 1st October", ylab="Triple-D samples", xaxt="n", type="h", lwd=.3, xlim=c(0,364))
axis_ticks <- seq(0, 360, by=30)  
axis(side=1, at=axis_ticks, labels=axis_ticks)
abline(h=c(seq(5,30,5)), lty=2, col="grey")
abline(v=92, lty=2, col="red")
#dev.off()

#Haul distance
table(trawldata$Dist)
table(trawldata$Survey2, trawldata$Dist)
```

load spatial information
```{r}
land <- st_read(file.path(path,"shapefiles/WGS84/WGS84geoKust500.shp"))
land <- as(land, "Spatial")
projection(land) <- CRS("+proj=longlat +datum=WGS84")
common.projection <- "+proj=utm +zone=31 +ellps=WGS84 
                    +datum=WGS84 +units=km +no_defs +towgs84=0,0,0"
land.utm <- spTransform(land,CRS(common.projection))

#Economic Exclusive Zone 
EEZ <- st_read(file.path(path,"shapefiles/World_EEZ_v12_20231025/eez_boundaries_v12.shp"))
EEZ <- as(EEZ, "Spatial")

#convert SpatialPolygonsDataFrame and to SpatialLinesDataFrames to a sf object
land_sf <- st_as_sf(land)
EEZ_sf <- st_as_sf(EEZ)

# Crop spatial data
bbox <- st_bbox(c(xmin = 0, xmax = 7, ymin = 52, ymax = 56), crs = st_crs(EEZ_sf))
EEZ_cropped <- st_crop(EEZ_sf, bbox)
```


#Fig. 1 
```{r}
#custum midpoint
custom_gradient <- gradient_n_pal(c("aliceblue", "lightblue1","dodgerblue", "blue",  "darkblue", "black"), values = c(0, 0.1, 0.2, 0.25, 0.5, 0.7)) #dark 
custom_gradient <- gradient_n_pal(c("aliceblue", "lightblue1","dodgerblue", "blue",  "darkblue", "black"), values = c(0, 0.1, 0.2, 0.5, 0.6, 0.7)) #light

xlim <- c(0,7.)
ylim <- c(52,56)

p <- ggplot() +
  geom_point(data = depth, aes(x = long, y = lat, color = depth), show.legend =F) +
  theme_bw() +
  scale_color_gradientn(colors = custom_gradient(seq(0, 1, length = 100)), name = "Depth (m)") +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  geom_sf(data = EEZ_sf, linetype = "dashed", size = 0.3, color="black") +
  coord_sf(xlim = xlim, ylim = ylim) + labs(x= "Longitude", y= "Latitude")  +
  geom_point(aes(x = trawldata$Lon_mid, y = trawldata$Lat_mid), color = alpha("brown", 1.0), shape=16, size=0.8) +
  theme(text = element_text(size = 7.34)) 

#inset
inset_map <- ggplot() + theme_void() +
  geom_sf(data = land_sf, fill = "grey90", color = "black") +
  theme(panel.background = element_rect(fill = "white"), panel.border = element_rect(color = "black", fill = NA, linewidth = 0.4)) +
  geom_point(aes(x = 3.6, y = 54), color = "red", size = 19, shape = 0) + 
  coord_sf(xlim = c(-10,11), ylim = c(48,60)) 

#text made in powerpoint
wadden <- readPNG(paste0(path,"/area_names/waddencoast_text.png"))
wadden <- annotation_custom(
  rasterGrob(wadden, interpolate = TRUE),
  xmin = 4.2, xmax = 6.2,  
  ymin = 53.0, ymax = 53.7  
)

holland <- readPNG(paste0(path, "/area_names/hollandcoast_text.png"))
holland <- annotation_custom(
  rasterGrob(holland, interpolate = TRUE),
  xmin = 3.6, xmax = 4.6,  
  ymin = 51.8, ymax = 53  
)

frisian <- readPNG(paste0(path, "/area_names/frisianfront_text.png"))
frisian <- annotation_custom(
  rasterGrob(frisian, interpolate = TRUE),
  xmin = 3.9, xmax = 5.45,  
  ymin = 53.45, ymax = 54.2  
)


oyster <- readPNG(paste0(path, "/area_names/oystergrounds_text.png"))
oyster <- annotation_custom(
  rasterGrob(oyster, interpolate = TRUE),
  xmin = 3.3, xmax = 5.3,  
  ymin = 54.15, ymax = 54.95  
)

dogger <- readPNG(paste0(path, "/area_names/doggerbank_text.png"))
dogger <- annotation_custom(
  rasterGrob(dogger, interpolate = TRUE),
  xmin = 1.8, xmax = 4,  
  ymin = 54.7, ymax = 55.8  
)

p <- p + wadden + holland + frisian + oyster + dogger  +
  annotation_custom(ggplotGrob(inset_map), 
                    xmin = 0, xmax = 2, 
                    ymin = 51, ymax = 54)

#ggsave(file = file.path(path, "figures/main/F1_locations.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)
```


#Fig. S3a
```{r}
#nearest point
#to sf
trawldata_sf <- st_as_sf(trawldata, coords = c("Lon_mid", "Lat_mid"), crs = 4326)
depth_sf <- st_as_sf(depth, coords = c("long", "lat"), crs = 4326)

nearest_indices <- st_nearest_feature(trawldata_sf, depth_sf) 

# Extract nearest information
nearest_depth <- depth_sf[nearest_indices, ]

# Calculate the distance to the nearest point
distances <- st_distance(trawldata_sf, nearest_depth, by_element = TRUE)

# Add nearest 
trawldata <- cbind(
  trawldata, 
  depth = nearest_depth$depth  # Nearest depth value
  #nearest_Lon = st_coordinates(nearest_depth)[, 1],  # Nearest Lon
  #nearest_Lat = st_coordinates(nearest_depth)[, 2],  # Nearest Lat
  #distance_m = round(as.numeric(distances),0)  # Distance in meters
)

ggplot() + 
  geom_point(data = trawldata, aes(x = Lon_mid, y = Lat_mid, color = depth), size=0.01) + 
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +  # Plot land polygons from land_sf
  scale_color_gradient(low = "cyan", high = "darkblue") + 
  theme_bw() + theme(text = element_text(size = 7.34)) +
  labs(
    title = "",
    x = "Longitude",
    y = "Latitude",
    color = "Depth"
  ) + xlim(xlim) + ylim(ylim)
#ggsave(file = file.path(path, "figures/supplement/S3_depth.pdf"), units = "mm", width = 85, height = 85, dpi = 1000)
```

#Fig. S3b
source: Bockelmann 2018: "Mapping mud content and median grain-size of North Sea sediments – A geostatistical approach"
```{r}
mud <- read.csv(file.path(path,"otherdata/mud.csv"))
mud <- mud[, !(colnames(mud)) == "X"]

#nearest point
#to sf
trawldata_sf <- st_as_sf(trawldata, coords = c("Lon_mid", "Lat_mid"), crs = 4326)
mud_sf <- st_as_sf(mud, coords = c("Longitude", "Latitude"), crs = 4326)

nearest_indices <- st_nearest_feature(trawldata_sf, mud_sf) 

# Extract nearest information
nearest_mud <- mud_sf[nearest_indices, ]

# Calculate the distance to the nearest point
distances <- st_distance(trawldata_sf, nearest_mud, by_element = TRUE)

# Add nearest
trawldata <- cbind(
  trawldata, 
  mud_perc = nearest_mud$mud_perc#,  # Nearest mud value
  #nearest_Lon = st_coordinates(nearest_mud)[, 1],  # Nearest Lon
  #nearest_Lat = st_coordinates(nearest_mud)[, 2],  # Nearest Lat
  #distance_m = round(as.numeric(distances),0)  # Distance in meters
)

ggplot() + 
  geom_point(data = trawldata, aes(x = Lon_mid, y = Lat_mid, color = mud_perc), size = 0.01) + 
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +  # Plot land polygons from land_sf
  scale_color_gradient(low = "yellow", high = "black") + 
  theme_bw() + theme(text = element_text(size = 7.34)) +
  labs(
    title = "",
    x = "Longitude",
    y = "Latitude",
    color = "Mud (%)"
  ) + xlim(xlim) + ylim(ylim)
#ggsave(file = file.path(path, "figures/supplement/S3_mud.pdf"), units = "mm", width = 85, height = 85, dpi = 1000)
```

```{r}
#save
#write.csv(trawldata, file = file.path(path, "trawldata_depth_mud.csv"))
```

#Fig. S1
spatial survey
```{r}
xlim <- c(0,7)
ylim <- c(52.2,55.7)

p <- ggplot() +
  geom_point(data = trawldata, aes(x = Lon_mid, y = Lat_mid), colour = "brown", size = 0.1, show.legend = FALSE) +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.2) +
  geom_sf(data = EEZ_cropped, fill = "grey90", color = "black", linewidth=0.3, linetype= "dashed") +
  # Facet by species
  facet_wrap(~ Survey2, ncol = 4) +
  labs(x = "Longitude", y = "Latitude") +
  coord_sf(xlim = xlim, ylim = ylim) + 
  theme_bw() +
  theme(
    strip.text = element_text(size = 6),  
    axis.text = element_text(size = 5),   
    axis.title = element_text(size = 7)   
  )

print(p)
#ggsave(file = file.path(path, "figures/supplement/S1_spatial_survey.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)
```

#Fish data
Individuals larger than 30 cm, whether flatfish or roundfish, are discarded. While the Triple-D catches them infrequently, they can still make a substantial local contribution to the observed biomass.
```{r}
large_individuals <- fishdata[fishdata$length_cm >= 30 & !is.na(fishdata$length_cm),]
stayin <- c("Ammodytidae", "Syngnathiformes", "Myxine glutinosa")
large_individuals <- large_individuals[!large_individuals$Species_reported %in% stayin,]
table(large_individuals$Species_reported)

#large fish in sample with smaller fish
Sample <- large_individuals[large_individuals$Weight_type == "Sample" & !is.na(large_individuals$Weight_type),]
Sample <- Sample[Sample$Species_reported != "Amblyraja radiata",] #these ray are both larger than 30cm

#update WW_g by removing the expected weight of the large individual
fishdata$StationID_fish <- paste0(fishdata$StationID, fishdata$Species_reported) #will be removed 
Sample$StationID_fish <- paste0(Sample$StationID, Sample$Species_reported) 

for ( i in unique(Sample$StationID_fish)) {
  fishdata[fishdata$StationID_fish == i,]$WW_g[1] <- sum(fishdata[fishdata$StationID_fish == i,]$WW_g, na.rm=T) - fishdata[fishdata$StationID_fish == i & fishdata$length_cm >= 30,]$expected_WW_g
}

#update fishdata
fishdata <- fishdata[!fishdata$unique_code %in% large_individuals$unique_code,-32]
```

Since not all entries have a WW_g (either 'Sample' or 'Entry'), a new column is made in which also includes the predicted weight of the fish where only length was measured (previously called 'wel_lengte'), only count data, and the broken fish.
```{r}
fishdata$Weight <- fishdata$WW_g
weight_type <- c("Entry", NA)
fishdata[is.na(fishdata$Weight) & fishdata$Weight_type %in% weight_type,]$Weight <- fishdata[is.na(fishdata$Weight) & fishdata$Weight_type %in% weight_type,]$WW_g_pred_group
```

Sometimes, a sample station contains multiple fractions. For example, this occurred when there were many gobies, but only half of the catch was specifically searched for gobies. While for other species the whole sample haul was sorted. 
```{r}
#which stations do have multiple 
multipleFractions <- as.data.frame(table(fishdata$StationID, fishdata$Fraction))
colnames(multipleFractions) <- c("StationID","Fraction", "Freq")

multipleFractions[multipleFractions$Freq > 0,]$Freq <- 1 #to create only present value
multipleFractions <- aggregate(Freq~StationID, multipleFractions, FUN="sum")
multipleFractions <- multipleFractions[multipleFractions$Freq > 1,]
```

Only (?) solution is to divide weight (count is not analysed) by fraction in the case there are multiple fractions
```{r}
fishdata[fishdata$StationID %in% multipleFractions$StationID,]$Weight <- fishdata[fishdata$StationID %in% multipleFractions$StationID,]$Weight/fishdata[fishdata$StationID %in% multipleFractions$StationID,]$Fraction

fishdata[fishdata$StationID %in% multipleFractions$StationID,]$Fraction <- 1
```

what are most caught fish in weight?
```{r}
#merge dist
fishdata <- merge(fishdata, trawldata[,c("StationID", "Dist")], by="StationID", all.x=T, all.y=F)
fishdata$Weight_cor <- fishdata$Weight/fishdata$Fraction/fishdata$Dist*5

#weight
weightperspecies <- aggregate(Weight_cor~Species_reported, fishdata, FUN="sum")
weightperspecies <- weightperspecies[order(weightperspecies$Weight_cor, decreasing = T),]
weightperspecies$WW_kg <- round(weightperspecies$Weight_cor/1000,1)
weightperspecies$relative_contribution <- round(weightperspecies$Weight_cor/sum(weightperspecies$Weight_cor)*100,1)
head(weightperspecies)

#numbers
fishdata$Count_cor <- fishdata$Count/fishdata$Fraction/fishdata$Dist*50
countperspecies <- aggregate(Count_cor~Species_reported, fishdata, FUN="sum")
countperspecies <- countperspecies[order(countperspecies$Count, decreasing = T),]
countperspecies$Count_cor <- round(countperspecies$Count_cor,0)
head(countperspecies,10)
```

Create English name of the 8 most common fish species
```{r}
fishdata$Species_English <- "Remaining"
fishdata[fishdata$Species_reported == "Gobiidae",]$Species_English <- "Goby spp."
fishdata[fishdata$Species_reported == "Ammodytidae",]$Species_English <- "Sandeel spp."
fishdata[fishdata$Species_reported == "Buglossidium luteum",]$Species_English <- "Solenette"
fishdata[fishdata$Species_reported == "Limanda limanda",]$Species_English <- "Dab"
fishdata[fishdata$Species_reported == "Arnoglossus laterna",]$Species_English <- "Scaldfish"
fishdata[fishdata$Species_reported == "Callionymus sp.",]$Species_English <- "Dragonet spp."
fishdata[fishdata$Species_reported == "Pleuronectes platessa",]$Species_English <- "Plaice"
```

#Fig. S5
Contribution to total observed weight per length. Weight per individual size is needed.
```{r}
#weight
fishdata$Weight_persize <- NA
fishdata[fishdata$Weight_type == "Entry" & !is.na(fishdata$Weight_type),]$Weight_persize <- fishdata[fishdata$Weight_type == "Entry" & !is.na(fishdata$Weight_type),]$WW_g
fishdata[is.na(fishdata$Weight_persize),]$Weight_persize <- fishdata[is.na(fishdata$Weight_persize),]$WW_g_pred_group

#length
fishdata$length_cm_all <- NA
fishdata$length_cm_all <- fishdata$length_cm
fishdata[is.na(fishdata$length_cm_all),]$length_cm_all <- fishdata[is.na(fishdata$length_cm_all),]$length_pred

#fraction and distance correction
fishdata$Weight_persize <-  fishdata$Weight_persize/fishdata$Fraction/fishdata$Dist*50 

length <- seq(0,max(fishdata$length_cm_all,na.rm = T)+1, 1)
cumalative <- data.frame("Length" = NA, Species_English= NA, "Proportion" = NA)
cumalative_species <- data.frame("Length" = length,Species_English = NA, "Proportion" = NA)

weight_type <- c(NA, "Sample")
for (j in unique(fishdata$Species_English)) {
for (i in length) {
cumalative_species$Species_English <- j 
cumalative_species[cumalative_species$Length == i & cumalative_species$Species_English == j,]$Proportion <- round(sum(fishdata[fishdata$length_cm_all <= i & fishdata$Species_English == j,]$Weight_persize,na.rm=T)  / sum(fishdata$Weight_persize, na.rm = T)*100, 1)

}
cumalative <- rbind(cumalative, cumalative_species)  
}

# order based on contribution
cumalative_max <- cumalative[cumalative$Length == 43 & !is.na(cumalative$Length),]
order <- cumalative_max$Species_English[order(cumalative_max$Proportion, decreasing = F)]
cumalative$Species_English <- factor(cumalative$Species_English,levels = order)

order <- c("Remaining", order[order != "Remaining"])
cumalative$Species_English <- factor(cumalative$Species_English,levels = order)

# Plot
ggplot(cumalative, aes(x = Length, y = Proportion / 100, fill = Species_English)) + 
  geom_area() + 
  labs(fill = 'Species') +  
  xlab("Fish length (cm)") + 
  ylab("Proportion") + 
  theme_bw() + 
  xlim(1, max(cumalative$Length, na.rm = TRUE)) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.1), position = "right") +
  theme(
    legend.position = c(0.01, 0.95), 
    legend.justification = c(0, 1) 
  )

#ggsave(file = file.path(path, "figures/supplement/S5_cumulative_biomass.pdf"), units = "mm", width = 170, height = 140, dpi = 1000)
```
Species specific length histogram
```{r}
#including length_pred
fishdata.expanded <- fishdata[rep(row.names(fishdata), fishdata$Count),]
fishdata.expanded$length_cm_all <- fishdata.expanded$length_cm
summary(fishdata.expanded$length_cm_all)
fishdata.expanded[is.na(fishdata.expanded $length_cm_all),]$length_cm_all <- fishdata.expanded[is.na(fishdata.expanded$length_cm_all),]$length_pred
summary(fishdata.expanded $length_cm_all)

ggplot(fishdata.expanded[fishdata.expanded$Species_English != "Remaining",], 
       aes(x = length_cm_all)) +
  # Place the lines before the histogram so that they appear behind the bars
  geom_vline(xintercept = seq(5, 30, by = 5), color = "grey", linetype = "dashed") +  # Manually set xintercepts
  geom_hline(yintercept = seq(0.1, 1, by = 0.1), color = "grey", linetype = "dashed") +  # Set yintercepts
  geom_histogram(aes(y = after_stat(density)), breaks = seq(0, 30, 1), 
                 fill = "white", color = "black") +  # Adjust density
  labs(x = "Length (cm)", y = "Proportion") +
  facet_wrap(~Species_English, scales = "free_y", nrow = 4, ncol = 2) +  # Facet by Species_English
  theme_bw() +
  theme(
    strip.background = element_rect(fill = "#E5E5E5"),  # Custom background for strip labels
    strip.text = element_text(size = 10),  # Font size of strip text
    axis.text.x = element_text(angle = 45, hjust = 1),  # Rotate x-axis labels
    panel.grid = element_blank(),  # Remove grid lines
    axis.ticks.length = unit(0.2, "cm")  # Adjust the length of ticks
  ) +
  scale_x_continuous(breaks = seq(0, 30, 5)) +  # Custom x-axis breaks
  scale_y_continuous(limits = c(0, 0.5), expand = c(0, 0))  # Remove space below the histogram
#ggsave(file = file.path(path, "figures/main/F2_species_hist.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)
```