---
title: "temporal variance fishbiomass"
author: "Bram Parmentier"
date: "2024-06-27"
output: html_document
---

# Temporal patterns of observed biomass for Parmentier et al (2025): "Small fish biomass in the North Sea is far greater than previously estimated"

Objective:
- gain insight in the temporal patterns of the observed Triple-D data. This is done in areas which are sample multiple times.
  -over years
  -seasonal effects

load libraries
```{r, warning=F}
library(xlsx);library(ggplot2)
library(sf)
library(INLA);library(mapdata)
library(stringr); library(fields)
library(gstat); library(lattice);library(latticeExtra)
library(dplyr); library(VGAM);library(raster); library(sp)
```

load data
```{r}
# Set locale to English
Sys.setlocale("LC_TIME", "C")

path <- ""
path <- "C:/Users/bparmentier/OneDrive - NIOZ/Documents/NIOZ/code/visbiomassa/DOI"

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

fishdata <- fishdata[!is.na(fishdata$StationID),]
fishdata <- fishdata[fishdata$StationID %in% trawldata$StationID,]
trawldata <- trawldata[trawldata$StationID %in%  fishdata$StationID,] 

fishdata$Fraction <- as.numeric(fishdata$Fraction)
```
Flat or round fish larger than 30cm are discarded. The used gear (Triple-D) caught them only very occasionally but they contribute largely to the observed biomass. However large (>=30cm) individuals of Ammodytidae, Syngnathiformes, and Myxine glutinosa are not removed. 
```{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,!names(fishdata) %in% ("StationID_fish")]
```

Not all records have a WW_g (either 'Sample' or 'Entry'), a new column is made which also includes the predicted weight of fish. This predicted weight (WW_g_pred_group) consist of weights derived from the length-weight species specific relationship (when a length was available), based on the length distribution of measured fish, or broken estimated weight. 
```{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
```

In some cases the catch is subsampled and also the fished distance differs between samples. In the models this is dealt with by applying an offset.
```{r}
fishdata$Fraction <- as.numeric(fishdata$Fraction)
table(fishdata$Fraction)
table(trawldata$Dist)
```

However, in some cases there are multiple fractions within one sample station. This happened, for example, when there were many gobies, and only in half of the catch is processed for gobies while other species the whole fraction 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,]
```

We solved this by dividing the weight (count is not analysed) by fraction in the case there are multiple fractions
```{r}
fishdata_fraction_cor <- fishdata
fishdata_fraction_cor[fishdata_fraction_cor$StationID %in% multipleFractions$StationID,]$Weight <- fishdata_fraction_cor[fishdata_fraction_cor$StationID %in% multipleFractions$StationID,]$Weight/fishdata_fraction_cor[fishdata_fraction_cor$StationID %in% multipleFractions$StationID,]$Fraction

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

combine trawl and fishdata_fraction_cor
```{r}
combined <- merge(fishdata_fraction_cor, trawldata, by="StationID", all.x = T) 
```

Weight per station
```{r}
#all fish
combined[is.na(combined$Weight),]$Weight <- 0 
totperstation <- aggregate(Weight ~ StationID, combined, sum, drop = FALSE) 

#add fraction and Date
totperstation <- merge(totperstation, combined[!duplicated(combined$StationID),c("StationID", "Fraction"),], by= "StationID", all.x = T )

#create utm
trawldata_sf <- st_as_sf(trawldata, coords = c("Lon_mid", "Lat_mid"), crs = 4326)
trawldata_utm <- st_transform(trawldata_sf, crs = 32631)
utm_coords <- st_coordinates(trawldata_utm)
trawldata$Xkm <- utm_coords[,1]/1000
trawldata$Ykm <- utm_coords[,2]/1000

#include trawl information with utm
totperstation <- merge(totperstation, trawldata, by= "StationID", all.x = T)
totperstation <- totperstation[!is.na(totperstation$Lat_mid),]

#exclude Survey which goes to the north part of north sea
totperstation <- totperstation[totperstation$Lat_mid < 56,]
```

observed weight per station, but corrected for fraction and fished distance
```{r}
totperstation$Weight_cor <- totperstation$Weight/totperstation$Fraction
totperstation$Weight_cor <- totperstation$Weight_cor/totperstation$Dist*5 #width of dregde is 0.2m, so times 5 to get to m2

summary(totperstation$Weight_cor)

# histogram
hist(totperstation$Weight_cor, breaks=100, xlab= "Wet weight (g m⁻²)", main="", ylab="Triple-D samples", 
     xaxt='n')
     axis(1, at=seq(0, max(totperstation$Weight_cor +10 , na.rm = TRUE), by=10), las=1)
     abline(v=c(quantile(totperstation$Weight_cor, probs = c(0.05, 0.95))), col="red", lty=2)
     abline(v=quantile(totperstation$Weight_cor, probs = 0.5), col="blue", lty=2)
```

map
```{r}
land <- st_read(file.path(path,"shapefiles/WGS84/WGS84geoKust500.shp"))
land <- as(land, "Spatial")
projection(land) <- CRS("+proj=longlat +datum=WGS84")
land_sf <- st_as_sf(land)
```
spatial Survey
```{r}
colors <- rainbow(length(unique(trawldata$Survey2)))
names(colors) <- unique(trawldata$Survey2)
trawldata <- trawldata[order(trawldata$Month),]
trawldata$Color <- colors[trawldata$Survey2]

plot(trawldata$Lon_mid, trawldata$Lat_mid, asp=1.5, col=alpha(trawldata$Color,0.5) , pch=16, cex=0.5)
lines(land)
```
#vlieland/texel area
+ polygon
```{r}
x <- c(4.95, 4.71, 4.66, 4.83)
y <- c(53.28, 53.31, 53.26, 53.22)
good_order <- order(atan2(x-mean(x),y-mean(y)))

plot(trawldata$Lon_mid, trawldata$Lat_mid, col=alpha(trawldata$Color,0.5) , pch=16, cex=0.5, ylim=c(53,53.5), xlim=c(4,5), asp=1.5, xlab="Longitute", ylab="Latitude")
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))

# Create a matrix of ordered coordinates & close the polygon by adding the first point again
coords <- rbind(cbind(x[good_order], y[good_order]), c(x[good_order[1]], y[good_order[1]]))

# Create an sf polygon
polygon_geom <- st_polygon(list(coords))  # Create the polygon geometry
polygon_sf <- st_sfc(polygon_geom, crs = 4326)  # Convert to simple feature
Vlieland_poly <- st_sf(geometry = polygon_sf)  # Convert to sf object
```


```{r}
# Convert trawldata points into an sf object
trawldata_sf <- st_as_sf(trawldata, coords = c("Lon_mid", "Lat_mid"), crs = 4326)

# Check which points are inside the polygon
inside <- st_within(trawldata_sf, Vlieland_poly, sparse = FALSE)

# Add a column to indicate whether the point is inside or outside
trawldata$inside_Vlieland <- inside

# Plot with different colors for inside/outside points
plot(trawldata$Lon_mid, trawldata$Lat_mid, col=ifelse(inside, "green", "red"),
     pch=16, cex=0.5, ylim=c(53,53.5), xlim=c(4,5), asp=1.5, xlab="Longitude", ylab="Latitude")
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))
```

Which Surveys?
```{r}
Vlieland <- trawldata[trawldata$inside_Vlieland  == T,]
table(Vlieland$Survey)
Vlieland <- Vlieland[Vlieland$Survey != "BSIK2007",] #only 1 sample

#Colors
Vlieland$Color <- "red"
Vlieland[Vlieland$Survey == "ZZ20",]$Color <- "blue"
Vlieland[Vlieland$Survey == "ZZ21",]$Color <- "green"

plot(Vlieland$Lon_mid, Vlieland$Lat_mid, asp=1.5, pch=16, col=alpha(Vlieland$Color,0.5))
lines(land)
```
```{r}
totperstation_vlieland <- totperstation[totperstation$StationID %in% Vlieland$StationID,]
totperstation_vlieland$Area <- "5 Vlieland"

table(totperstation_vlieland$Fraction)
table(totperstation_vlieland$Survey, totperstation_vlieland$Month)
table(totperstation_vlieland$Survey, totperstation_vlieland$Year) #binnen 1.5 jaar van elkaar -> vooral seizoens effect

ggplot(totperstation_vlieland, aes(Survey, Weight/Dist*5)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

aggregate(totperstation_vlieland$Weight/totperstation_vlieland$Dist*5~totperstation_vlieland$Survey, FUN= "mean") 
aggregate(totperstation_vlieland$Weight/totperstation_vlieland$Dist*5~totperstation_vlieland$Survey, FUN= "sd")
table(totperstation_vlieland$Survey2)
```


#FrieseFront
```{r}
x <- c(5.14, 4.23, 4.23, 5.14)
y <- c(54.18, 53.82, 53.40, 53.80)
good_order <- order(atan2(x-mean(x),y-mean(y)))

plot(trawldata$Lon_mid, trawldata$Lat_mid, col=alpha(trawldata$Color,0.5) , pch=16, cex=0.5, ylim=c(53.3,54.5), xlim=c(4,5), asp=1.5)
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))

# Create a matrix of ordered coordinates & close the polygon by adding the first point again
coords <- rbind(cbind(x[good_order], y[good_order]), c(x[good_order[1]], y[good_order[1]]))

# Create an sf polygon
polygon_geom <- st_polygon(list(coords))  # Create the polygon geometry
polygon_sf <- st_sfc(polygon_geom, crs = 4326)  # Convert to simple feature
FrisianFront_poly <- st_sf(geometry = polygon_sf)  # Convert to sf object
```

```{r}
# Check which points are inside the polygon
inside <- st_within(trawldata_sf, FrisianFront_poly, sparse = FALSE)

# Add a column to indicate whether the point is inside or outside
trawldata$inside_FrisianFront <- inside

# Plot with different colors for inside/outside points
plot(trawldata$Lon_mid, trawldata$Lat_mid, col=ifelse(inside, "green", "red"),
     pch=16, cex=0.5, ylim=c(53,54.5), xlim=c(4,5), asp=1.5, xlab="Longitude", ylab="Latitude")
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))
```

which Surveys?
```{r}
friesefront <- trawldata[trawldata$inside_FrisianFront == T,]
table(friesefront$Survey)

niet <- c("64PE340", "64PE438", "COSTRA2007") #only 4,1,1 hauls
friesefront <- friesefront[!friesefront$Survey %in% niet,]
table(friesefront$Survey)

#duidelijkere kleuren
friesefront$Color <- "red"
friesefront[friesefront$Survey == "ARCA2006",]$Color <- "blue"
friesefront[friesefront$Survey == "BSIK2007",]$Color <- "green"
friesefront[friesefront$Survey == "KRM21",]$Color <- "black"
friesefront[friesefront$Survey == "KRM22",]$Color <- "orange"
friesefront[friesefront$Survey == "RWS31144108",]$Color <- "purple"
friesefront[friesefront$Survey == "ZZ19",]$Color <- "brown"

plot(friesefront$Lon_mid, friesefront$Lat_mid, asp=1.5, pch=16, col=alpha(friesefront$Color,0.5))
lines(land)
```

```{r}
totperstation_friesefront <- totperstation[totperstation$StationID %in% friesefront$StationID,]
totperstation_friesefront$Area <- "3 Frisian Front"
table(totperstation_friesefront$Fraction)

totperstation_friesefront$Year_month <- paste0(totperstation_friesefront$Year, totperstation_friesefront$Month)
table(totperstation_friesefront$Survey, totperstation_friesefront$Year)
table(totperstation_friesefront$Survey, totperstation_friesefront$Month) 
table(totperstation_friesefront$Survey, totperstation_friesefront$Year_month)

totperstation_friesefront$month_Survey <- paste0(totperstation_friesefront$Month, totperstation_friesefront$Survey)
totperstation_friesefront$month_year <- paste0(totperstation_friesefront$Month, totperstation_friesefront$Year)

ggplot(totperstation_friesefront, aes(Survey , Weight/Dist*5, WW_g)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

ggplot(totperstation_friesefront, aes(month_year , Weight/Dist*5, WW_g)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

aggregate(totperstation_friesefront$Weight/totperstation_friesefront$Dist*5~totperstation_friesefront$Survey2, FUN= "mean") 
aggregate(totperstation_friesefront$Weight/totperstation_friesefront$Dist*5~totperstation_friesefront$Survey2, FUN= "sd") 
table(totperstation_friesefront$Survey2)
```

#FrieseFront_zuid
```{r}
x <- c(5.326, 5.43, 4.976, 4.93)
y <- c(53.85, 53.64, 53.50, 53.65)
good_order <- order(atan2(x-mean(x),y-mean(y)))

plot(trawldata$Lon_mid, trawldata$Lat_mid, col=alpha(trawldata$Color,0.5) , pch=16, cex=0.5, ylim=c(53.3,54.5), xlim=c(4,5), asp=1.5)
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))

# Create a matrix of ordered coordinates & close the polygon by adding the first point again
coords <- rbind(cbind(x[good_order], y[good_order]), c(x[good_order[1]], y[good_order[1]]))

# Create an sf polygon
polygon_geom <- st_polygon(list(coords))  # Create the polygon geometry
polygon_sf <- st_sfc(polygon_geom, crs = 4326)  # Convert to simple feature
FrisianFrontsouth_poly <- st_sf(geometry = polygon_sf)  # Convert to sf object
```

```{r}
# Check which points are inside the polygon
inside <- st_within(trawldata_sf, FrisianFrontsouth_poly , sparse = FALSE)

# Add a column to indicate whether the point is inside or outside
trawldata$inside_FrisianFrontsouth <- inside

# Plot with different colors for inside/outside points
plot(trawldata$Lon_mid, trawldata$Lat_mid, col=ifelse(inside, "green", "red"),
     pch=16, cex=0.5, ylim=c(53,54.5), xlim=c(4,5), asp=1.5, xlab="Longitude", ylab="Latitude")
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))
```

```{r}
friesefront_zuid <- trawldata[trawldata$inside_FrisianFrontsouth  == 1,]
table(friesefront_zuid$Survey)
niet <- c("64PE338", "BSIK2007", "Normomap2008", "ZZ21") 
friesefront_zuid <- friesefront_zuid[!friesefront_zuid$Survey %in% niet,]

#colors
friesefront_zuid$Color <- "red"
friesefront_zuid[friesefront_zuid$Survey == "RWS31144108",]$Color <- "purple"

plot(friesefront_zuid$Lon_mid, friesefront_zuid$Lat_mid, asp=1.5, pch=16, col=alpha(friesefront_zuid$Color,0.5))
lines(land)
```
```{r}
totperstation_friesefront_zuid <- totperstation[totperstation$StationID %in% friesefront_zuid$StationID,]
totperstation_friesefront_zuid$Area <- "4 Frisian Front south"
table(totperstation_friesefront_zuid$Fraction)

totperstation_friesefront_zuid$Year_month <- paste0(totperstation_friesefront_zuid$Year, totperstation_friesefront_zuid$Month)
table(totperstation_friesefront_zuid$Survey, totperstation_friesefront_zuid$Year)
table(totperstation_friesefront_zuid$Survey, totperstation_friesefront_zuid$Month)
table(totperstation_friesefront_zuid$Survey, totperstation_friesefront_zuid$Year_month)

totperstation_friesefront_zuid$month_Survey <- paste0(totperstation_friesefront_zuid$Month, totperstation_friesefront_zuid$Survey)
totperstation_friesefront_zuid$month_year <- paste0(totperstation_friesefront_zuid$Month, totperstation_friesefront_zuid$Year)

ggplot(totperstation_friesefront_zuid, aes(Survey , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

ggplot(totperstation_friesefront_zuid, aes(month_Survey , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

ggplot(totperstation_friesefront_zuid, aes(month_year , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

aggregate(totperstation_friesefront_zuid$Weight/totperstation_friesefront_zuid$Dist*5~totperstation_friesefront_zuid$Survey2, FUN= "mean") 
aggregate(totperstation_friesefront_zuid$Weight/totperstation_friesefront_zuid$Dist*5~totperstation_friesefront_zuid$Survey2, FUN= "sd") 
table(totperstation_friesefront_zuid$Survey2)
```


#Oyster grounds
```{r}
x <- c(3.99, 4.55, 4.55, 4.1, 3.62)
y <- c(54.3, 54.4, 54.9, 55.1, 54.76)
good_order <- order(atan2(x-mean(x),y-mean(y)))

plot(trawldata$Lon_mid, trawldata$Lat_mid, col=alpha(trawldata$Color,0.5) , pch=16, cex=0.5, ylim=c(53.8,55.3), xlim=c(4,5), asp=1.5)
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))

# Create a matrix of ordered coordinates & close the polygon by adding the first point again
coords <- rbind(cbind(x[good_order], y[good_order]), c(x[good_order[1]], y[good_order[1]]))

# Create an sf polygon
polygon_geom <- st_polygon(list(coords))  # Create the polygon geometry
polygon_sf <- st_sfc(polygon_geom, crs = 4326)  # Convert to simple feature
oystergrounds_poly <- st_sf(geometry = polygon_sf)  # Convert to sf object
```

```{r}
# Check which points are inside the polygon
inside <- st_within(trawldata_sf, oystergrounds_poly  , sparse = FALSE)

# Add a column to indicate whether the point is inside or outside
trawldata$inside_Oystergrounds <- inside

# Plot with different colors for inside/outside points
plot(trawldata$Lon_mid, trawldata$Lat_mid, col=ifelse(inside, "green", "red"),
     pch=16, cex=0.5, ylim=c(53,55.5), xlim=c(4,5), asp=1.5, xlab="Longitude", ylab="Latitude")
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))
```
```{r}
oestergronden <- trawldata[trawldata$inside_Oystergrounds  == T,]
table(oestergronden$Survey)
niet <- c("64PE340", "64PE438", "Normomap2008", "COSTRA2007", "KRM21")
oestergronden <- oestergronden[!oestergronden$Survey %in% niet,]
table(oestergronden$Survey)

#Colors
oestergronden$Color <- "red"
oestergronden[oestergronden$Survey == "KRM22",]$Color <- "blue"
oestergronden[oestergronden$Survey == "RWS31144108",]$Color <- "green"

plot(oestergronden$Lon_mid, oestergronden$Lat_mid, asp=1.5, pch=16, col=alpha(oestergronden$Color,0.5))
lines(land)
```
```{r}
totperstation_oestergronden <- totperstation[totperstation$StationID %in% oestergronden$StationID,]
totperstation_oestergronden$Area <- "2 Oyster Grounds"

table(totperstation_oestergronden$Fraction)

totperstation_oestergronden$Year_month <- paste0(totperstation_oestergronden$Year, totperstation_oestergronden$Month)
table(totperstation_oestergronden$Survey, totperstation_oestergronden$Year)
table(totperstation_oestergronden$Survey, totperstation_oestergronden$Month)
table(totperstation_oestergronden$Survey, totperstation_oestergronden$Year_month)

totperstation_oestergronden$month_Survey <- paste0(totperstation_oestergronden$Month, totperstation_oestergronden$Survey)
totperstation_oestergronden$month_year <- paste0(totperstation_oestergronden$Month, totperstation_oestergronden$Year)

ggplot(totperstation_oestergronden, aes(month_Survey , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

ggplot(totperstation_oestergronden, aes(month_year , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

aggregate(totperstation_oestergronden$Weight/totperstation_oestergronden$Dist*5~totperstation_oestergronden$Survey2, FUN= "mean")

aggregate(totperstation_oestergronden$Weight/totperstation_oestergronden$Dist*5~totperstation_oestergronden$Survey2, FUN= "sd")
table(totperstation_oestergronden$Survey2)
```

#NL Doggerbank
```{r}
x <- c(2.875, 4.25, 4.25, 2.875)
y <- c(54.451, 55.35, 55.7, 55.7)
good_order <- order(atan2(x-mean(x),y-mean(y)))

plot(trawldata$Lon_mid, trawldata$Lat_mid, col=alpha(trawldata$Color,0.5) , pch=16, cex=0.5, ylim=c(54.5,56.3), xlim=c(4,5), asp=1.5)
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))

# Create a matrix of ordered coordinates & close the polygon by adding the first point again
coords <- rbind(cbind(x[good_order], y[good_order]), c(x[good_order[1]], y[good_order[1]]))

# Create an sf polygon
polygon_geom <- st_polygon(list(coords))  # Create the polygon geometry
polygon_sf <- st_sfc(polygon_geom, crs = 4326)  # Convert to simple feature
doggerbankNL_poly <- st_sf(geometry = polygon_sf)  # Convert to sf object
```
```{r}
# Check which points are inside the polygon
inside <- st_within(trawldata_sf, doggerbankNL_poly  , sparse = FALSE)

# Add a column to indicate whether the point is inside or outside
trawldata$inside_DoggerBank <- inside

# Plot with different colors for inside/outside points
plot(trawldata$Lon_mid, trawldata$Lat_mid, col=ifelse(inside, "green", "red"),
     pch=16, cex=0.5, ylim=c(53,55.5), xlim=c(4,5), asp=1.5, xlab="Longitude", ylab="Latitude")
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))
```

```{r}
doggerbank <- trawldata[trawldata$inside_DoggerBank  == T,]
table(doggerbank$Survey)
niet <- c("64PE340", "64PE438")
doggerbank <- doggerbank[!doggerbank$Survey %in% niet,]
table(doggerbank$Survey)

#duidelijkere kleuren
doggerbank$Color <- "red"
doggerbank[doggerbank$Survey == "KRM21",]$Color <- "blue"
doggerbank[doggerbank$Survey == "KRM22",]$Color <- "green"
doggerbank[doggerbank$Survey == "RWS31144108",]$Color <- "black"

plot(doggerbank$Lon_mid, doggerbank$Lat_mid, asp=1.5, pch=16, col=alpha(doggerbank$Color,0.5))
lines(land)

plot(doggerbank[doggerbank$Survey == "KRM22",]$Lon_mid, doggerbank[doggerbank$Survey == "KRM22",]$Lat_mid, asp=1.5, pch=16, col=alpha(doggerbank[doggerbank$Survey == "KRM22",]$Color,0.5))
lines(land)
```
```{r}
totperstation_doggerbank <- totperstation[totperstation$StationID %in% doggerbank$StationID,]
totperstation_doggerbank$Area <- "1 Dogger Bank NL"

table(totperstation_doggerbank$Fraction)
totperstation_doggerbank$Weight_fraction <- totperstation_doggerbank$Weight/totperstation_doggerbank$Fraction

totperstation_doggerbank$Year_month <- paste0(totperstation_doggerbank$Year, totperstation_doggerbank$Month)
table(totperstation_doggerbank$Survey, totperstation_doggerbank$Year)
table(totperstation_doggerbank$Survey, totperstation_doggerbank$Month)
table(totperstation_doggerbank$Survey, totperstation_doggerbank$Year_month)

totperstation_doggerbank$month_Survey <- paste0(totperstation_doggerbank$Month, totperstation_doggerbank$Survey)
totperstation_doggerbank$month_year <- paste0(totperstation_doggerbank$Month, totperstation_doggerbank$Year)

ggplot(totperstation_doggerbank, aes(month_Survey , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

ggplot(totperstation_doggerbank, aes(month_year , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

aggregate(totperstation_doggerbank$Weight_fraction/totperstation_doggerbank$Dist*5~totperstation_doggerbank$Survey2, FUN= "mean") 
aggregate(totperstation_doggerbank$Weight_fraction/totperstation_doggerbank$Dist*5~totperstation_doggerbank$Survey2, FUN= "sd")
table(totperstation_doggerbank$Survey2)
```


#Holland kust noord
```{r}
x <- c(4.71, 4.478, 4.55, 4.21)
y <- c(52.95, 52.98, 52.47, 52.52)
good_order <- order(atan2(x-mean(x),y-mean(y)))

plot(trawldata$Lon_mid, trawldata$Lat_mid, col=alpha(trawldata$Color,0.5) , pch=16, cex=0.5, ylim=c(52,53), xlim=c(3.5,4.5), asp=1.5)
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))

# Create a matrix of ordered coordinates & close the polygon by adding the first point again
coords <- rbind(cbind(x[good_order], y[good_order]), c(x[good_order[1]], y[good_order[1]]))

# Create an sf polygon
polygon_geom <- st_polygon(list(coords))  # Create the polygon geometry
polygon_sf <- st_sfc(polygon_geom, crs = 4326)  # Convert to simple feature
holland_north_poly <- st_sf(geometry = polygon_sf)  # Convert to sf object
```

```{r}
# Check which points are inside the polygon
inside <- st_within(trawldata_sf, holland_north_poly   , sparse = FALSE)

# Add a column to indicate whether the point is inside or outside
trawldata$inside_Holland_north <- inside

# Plot with different colors for inside/outside points
plot(trawldata$Lon_mid, trawldata$Lat_mid, col=ifelse(inside, "green", "red"),
     pch=16, cex=0.5, ylim=c(52,53.5), xlim=c(4,5), asp=1.5, xlab="Longitude", ylab="Latitude")
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))
```
```{r}
coast_west <- trawldata[trawldata$inside_Holland_north  == T,]
table(coast_west$Survey)
niet <- c("SANDEXTR22") #, "PAWP22")
coast_west <- coast_west[!coast_west$Survey %in% niet,]
table(coast_west$Survey)

#duidelijkere kleuren
coast_west$Color <- "red"
coast_west[coast_west$Survey == "COSTRA2007",]$Color <- "blue"
coast_west[coast_west$Survey == "Noordzeewind2007",]$Color <- "green"
coast_west[coast_west$Survey == "ZZ20",]$Color <- "black"
coast_west[coast_west$Survey == "ZZ21",]$Color <- "purple"
coast_west[coast_west$Survey == "PAWP22",]$Color <- "gold"
#coast_west[coast_west$Survey == "SANDEXTR22",]$Color <- "yellow"

plot(coast_west$Lon_mid, coast_west$Lat_mid, asp=1.5, pch=16, col=alpha(coast_west$Color,0.5))
lines(land)
```
```{r}
totperstation_coast_west <- totperstation[totperstation$StationID %in% coast_west$StationID,]
totperstation_coast_west$Area <- "6 Holland coast north"

table(totperstation_coast_west$Fraction)
totperstation_coast_west$Weight_fraction <- totperstation_coast_west$Weight/totperstation_coast_west$Fraction

totperstation_coast_west$Year_month <- paste0(totperstation_coast_west$Year, totperstation_coast_west$Month)
table(totperstation_coast_west$Survey, totperstation_coast_west$Year)
table(totperstation_coast_west$Survey, totperstation_coast_west$Month)
table(totperstation_coast_west$Survey, totperstation_coast_west$Year_month)

totperstation_coast_west$month_Survey <- paste0(totperstation_coast_west$Month, totperstation_coast_west$Survey)
totperstation_coast_west$month_year <- paste0(totperstation_coast_west$Month, totperstation_coast_west$Year)

ggplot(totperstation_coast_west, aes(month_Survey , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

ggplot(totperstation_coast_west, aes(month_year , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

aggregate(totperstation_coast_west$Weight_fraction/totperstation_coast_west$Dist*5~totperstation_coast_west$Survey2, FUN= "mean") # 
aggregate(totperstation_coast_west$Weight_fraction/totperstation_coast_west$Dist*5~totperstation_coast_west$Survey2, FUN= "sd")
table(totperstation_coast_west$Survey2)
```

#Westkust_zuid
```{r}
x <- c(4.52, 4.28, 3.90, 4.17, 4.45)
y <- c(52.44, 52.5, 52.07, 52.0, 52.3)
good_order <- order(atan2(x-mean(x),y-mean(y)))

plot(trawldata$Lon_mid, trawldata$Lat_mid, col=alpha(trawldata$Color,0.5) , pch=16, cex=0.5, ylim=c(52,53), xlim=c(3.5,4.5), asp=1.5)
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))

# Create a matrix of ordered coordinates & close the polygon by adding the first point again
coords <- rbind(cbind(x[good_order], y[good_order]), c(x[good_order[1]], y[good_order[1]]))

# Create an sf polygon
polygon_geom <- st_polygon(list(coords))  # Create the polygon geometry
polygon_sf <- st_sfc(polygon_geom, crs = 4326)  # Convert to simple feature
holland_south_poly <- st_sf(geometry = polygon_sf)  # Convert to sf object
```

```{r}
# Check which points are inside the polygon
inside <- st_within(trawldata_sf, holland_south_poly   , sparse = FALSE)

# Add a column to indicate whether the point is inside or outside
trawldata$inside_Holland_south <- inside

# Plot with different colors for inside/outside points
plot(trawldata$Lon_mid, trawldata$Lat_mid, col=ifelse(inside, "green", "red"),
     pch=16, cex=0.5, ylim=c(52,53.5), xlim=c(4,5), asp=1.5, xlab="Longitude", ylab="Latitude")
lines(land)
polygon(x[good_order], y[good_order], col=alpha("blue", 0.5))
```
```{r}
coast_west_south <- trawldata[trawldata$inside_Holland_south  == T,]
table(coast_west_south$Survey)
niet <- c("SANDEXTR22", "Noordzeewind2007", "64PE330")
coast_west_south <- coast_west_south[!coast_west_south$Survey %in% niet,]
table(coast_west_south$Survey)

#duidelijkere kleuren
coast_west_south$Color <- "red"
coast_west_south[coast_west_south$Survey == "ZZ22",]$Color <- "blue"
coast_west_south[coast_west_south$Survey == "64PE338",]$Color <- "green"

plot(coast_west_south$Lon_mid, coast_west_south$Lat_mid, asp=1.5, pch=16, col=alpha(coast_west_south$Color,0.5))
lines(land)
```
```{r}
totperstation_coast_west_south <- totperstation[totperstation$StationID %in% coast_west_south$StationID,]
totperstation_coast_west_south$Area <- "7 Holland coast south"
table(totperstation_coast_west_south$Fraction)
totperstation_coast_west_south$Weight_fraction <- totperstation_coast_west_south$Weight/totperstation_coast_west_south$Fraction

totperstation_coast_west_south$Year_month <- paste0(totperstation_coast_west_south$Year, totperstation_coast_west_south$Month)
table(totperstation_coast_west_south$Survey, totperstation_coast_west_south$Year)
table(totperstation_coast_west_south$Survey, totperstation_coast_west_south$Month)
table(totperstation_coast_west_south$Survey, totperstation_coast_west_south$Year_month)

totperstation_coast_west_south$month_Survey <- paste0(totperstation_coast_west_south$Month, totperstation_coast_west_south$Survey)
totperstation_coast_west_south$month_year <- paste0(totperstation_coast_west_south$Month, totperstation_coast_west_south$Year)

ggplot(totperstation_coast_west_south, aes(month_Survey , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

ggplot(totperstation_coast_west_south, aes(month_year , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

ggplot(totperstation_coast_west_south, aes(Survey , Weight/Dist*5, Weight)) + 
  geom_jitter(width=0.13, size = 3) + xlab("Survey") + ylab("WW (g)")

aggregate(totperstation_coast_west_south$Weight_fraction/totperstation_coast_west_south$Dist*5~totperstation_coast_west_south$Survey2, FUN= "mean") 
aggregate(totperstation_coast_west_south$Weight_fraction/totperstation_coast_west_south$Dist*5~totperstation_coast_west_south$Survey2, FUN= "sd") 
table(totperstation_coast_west_south$Survey2)
```
#Area overview
```{r}
selectie_Surveys <- unique(c(totperstation_vlieland$Survey2, totperstation_friesefront$Survey2, totperstation_friesefront_zuid$Survey2, totperstation_oestergronden$Survey2, totperstation_doggerbank$Survey2, totperstation_coast_west$Survey2, totperstation_coast_west_south$Survey2))

xlim <- c(2,6)
ylim <- c(52,56)

#Fig S12
ggplot() + 
  geom_point(data = trawldata[trawldata$Survey2 %in% selectie_Surveys,], aes(x = Lon_mid, y = Lat_mid), color = "black", size = 0.6) +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +  # Plot land polygons from land_sf
  geom_sf(data = doggerbankNL_poly, fill = alpha("#1B9E77", 0.5), color = "#1B9E77") +
  geom_sf(data = oystergrounds_poly, fill = alpha("#D95F02", 0.5), color = "#D95F02") +
  geom_sf(data = FrisianFront_poly, fill = alpha("#7570B3", 0.5), color = "#7570B3") +
  geom_sf(data = FrisianFrontsouth_poly, fill = alpha("#E7298A", 0.5), color = "#E7298A") +
  geom_sf(data = Vlieland_poly, fill = alpha("#66A61E", 0.5), color = "#66A61E") +
  geom_sf(data = holland_south_poly, fill = alpha("#A6761D", 0.5), color = "#A6761D") +
  geom_sf(data = holland_north_poly, fill = alpha("#E6AB02", 0.5), color = "#E6AB02") +
  theme_bw() +
  labs(
    x = "Longitude",
    y = "Latitude"
  ) + xlim(xlim) + ylim(ylim)
ggsave(file.path(path, "figures/supplement/S12_temporal_biomass_areas.pdf"), width = 170, height = 255, dpi = 1000, units = "mm")
```

#Temporal plots
fraction
```{r}
#combine areas
totperstation_areas <- rbind(totperstation_vlieland, totperstation_friesefront[, names(totperstation_friesefront) %in% names(totperstation_vlieland)], totperstation_friesefront_zuid[, names(totperstation_friesefront_zuid) %in% names(totperstation_vlieland)], totperstation_friesefront_zuid[, names(totperstation_friesefront_zuid) %in% names(totperstation_vlieland)], totperstation_friesefront_zuid[, names(totperstation_friesefront_zuid) %in% names(totperstation_vlieland)], totperstation_coast_west[, names(totperstation_coast_west) %in% names(totperstation_vlieland)], totperstation_coast_west_south[, names(totperstation_coast_west_south) %in% names(totperstation_vlieland)], totperstation_doggerbank[,names(totperstation_doggerbank) %in% names(totperstation_vlieland)], totperstation_oestergronden[,names(totperstation_oestergronden) %in% names(totperstation_vlieland)])

#weight per area and Survey
meanperSurvey <- aggregate(Weight/Dist*5 ~ Area + Survey2, data=totperstation_areas, FUN="mean")
names(meanperSurvey )[names(meanperSurvey ) == "Weight/Dist * 5"] <- 'Weight_cor'

#add date
totperstation_areas$Date <- as.Date(totperstation_areas$Date)
mean_date <- aggregate(Date ~ Area + Survey2, data=totperstation_areas, FUN="mean")
meanperSurvey <- merge(meanperSurvey, mean_date, by=c("Survey2", "Area"), all.x=T)

#plot
# Create a sequence for the start of each year
year_starts <- seq(from = as.Date("2006-01-01"), to = as.Date("2024-01-01"), by = "1 year")

#Fig 13
ggplot(meanperSurvey, aes(x = Date, y = Weight_cor, color = Area, group = Area)) +
  geom_point(size = 1.2) +  
  geom_line(size = 0.5, linetype = "dotted") +   
  theme_bw(base_size = 6) +  
  labs(
    x = "Date",
    y = expression(Weight~observed~(g~m^{-2})),
    color = "Area",
    title = ""
  ) +
  scale_color_brewer(palette = "Dark2") +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.text = element_text(size = 5),      
    legend.title = element_text(size = 6),
    legend.key.height = unit(0.5, "lines"),
    axis.text = element_text(size = 5),        
    axis.title = element_text(size = 6),       
    plot.title = element_text(size = 6)       
  ) +
  geom_vline(xintercept = year_starts, linetype = "solid", color = alpha("black", 0.2))

ggsave(file.path(path, "figures/supplement/S13_annual_biomass.pdf"), width = 170, height = 85, dpi = 1000, units = "mm")


#within a year
meanperSurvey$DOY <- as.POSIXlt(meanperSurvey$Date)$yday + 1 
meanperSurvey$DOY_adjusted <- ifelse(meanperSurvey$DOY >= 274,  
                                 meanperSurvey$DOY - 273,   
                                 meanperSurvey$DOY + 92)

ggplot(meanperSurvey, aes(x = DOY_adjusted, y = Weight_cor, color = Area, group = Area)) +
  geom_point(size = 1.2) +  
  geom_line(size = 0.5, linetype = "dotted") +   
  theme_bw(base_size = 6) +  
  labs(
    x = "Day from 1st October",
    y = expression(Weight~observed~(g~m^{-2})),
    color = "Area",
    title = ""
  ) +
  scale_color_brewer(palette = "Dark2") +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.text = element_text(size = 5),      
    legend.title = element_text(size = 6),
    legend.key.height = unit(0.5, "lines"),
    axis.text = element_text(size = 5),        
    axis.title = element_text(size = 6),       
    plot.title = element_text(size = 6)       
  ) 
ggsave(file.path(path, "figures/supplement/S13_seasonal_biomass.pdf"), width = 170, height = 85, dpi = 1000, units = "mm")
```