diff --git a/moneos_2024/080_hyperbenthos/080_hyperbenthos_data.Rmd b/moneos_2024/080_hyperbenthos/080_hyperbenthos_data.Rmd index 165cad6..575c0fe 100644 --- a/moneos_2024/080_hyperbenthos/080_hyperbenthos_data.Rmd +++ b/moneos_2024/080_hyperbenthos/080_hyperbenthos_data.Rmd @@ -26,6 +26,7 @@ library(lubridate) library(readxl) library(writexl) library(RODBC) +library(ggpubr) ``` @@ -48,7 +49,7 @@ pad_tabellen <- maak_pad(params$hoofdstuk, "tabellen") ```{r data} #eerst recente kopie van op Amazon op G gezet -MDB <- odbcConnectAccess2007("G:/.shortcut-targets-by-id/0B0xcP-eNvJ9dZDBwVVJOVk5Ld2s/PRJ_SCHELDE/Benthos/HyperEpibenthos/DATA/HYPERBENTHOS_SCHELDEaug2024.accdb") +MDB <- odbcConnectAccess2007("G:/.shortcut-targets-by-id/0B0xcP-eNvJ9dZDBwVVJOVk5Ld2s/PRJ_SCHELDE/Benthos/HyperEpibenthos/DATA/HYPERBENTHOS_SCHELDEdec2024.accdb") # querry met mulitply (substaal nr) al gebruikt voor aantal, WW en AFDW, daardoor is AFDW niet meer DW-AW in uiteindelijke tabel! sqlCode <- " @@ -87,27 +88,140 @@ tel %>% ``` + + ```{r} -# paar ontbrekende AFDW aanvullen op basis van geschatte regressie WW-AFDW, vooral relaties gebruiken met hogere aantallen, want vaak afwijkingen bij zr lage gewichten. OOk enkele neg AFDW vervangen op zelfde manier -# Pomatoschistus soms als sp. (kleintjes en soms als microps, maar wrsch bijna altijd microps. Omdat dit voor taxonrijkdom en procentuele bijdrage lastig is, hier samen gevoegd) - -telc <- tel %>% - dplyr::mutate(AFDW = if_else(DW==0 & soort=="Gasterosteus aculeatus", WW/6, AFDW), - AFDW = if_else(DW==0 & soort == "Gammarus tigrinus", WW/5.5, AFDW), - AFDW = if_else(DW==0 & soort == "Pomatoschistus microps", WW/5.8, AFDW), - AFDW = if_else(DW==0 & soort == "Limnomysis benedeni", WW/4.7, AFDW), - AFDW = if_else(DW==0 & soort =="Bathyporeia pilosa", WW/6, AFDW), - AFDW = if_else(AFDW<0 & soort == "Abramis brama", WW/7.7, AFDW), - AFDW = if_else(AFDW<0 & soort=="Synidotea laticauda", WW/7, AFDW), - AFDW = if_else(WW% - dplyr::filter(WW% + group_by(soort) %>% + mutate( + # Calculate WW/AFDW ratio + ratio_WW_AFDW = WW / AFDW, + + # Calculate species-level quantiles and IQR for WW/AFDW ratio + Q1_ratio = as.numeric(quantile(ratio_WW_AFDW, 0.25, na.rm = TRUE)), + Q3_ratio = as.numeric(quantile(ratio_WW_AFDW, 0.75, na.rm = TRUE)), + IQR_ratio = Q3_ratio - Q1_ratio, + lower_bound_ratio = Q1_ratio - 1.5 * IQR_ratio, + upper_bound_ratio = Q3_ratio + 1.5 * IQR_ratio, + + # Define outliers for WW/AFDW ratio + outlier_AFDW_WW = ifelse( + ratio_WW_AFDW < lower_bound_ratio | ratio_WW_AFDW > upper_bound_ratio, + "Outlier", + "Non-Outlier" + ), + + # Identify inliers for slope calculation + inlier_WW_AFDW = !is.na(ratio_WW_AFDW) & + ratio_WW_AFDW >= lower_bound_ratio & + ratio_WW_AFDW <= upper_bound_ratio + ) %>% + mutate( + # Fit linear model for inliers only and extract slope + slopeAFDW_WW = if (any(inlier_WW_AFDW, na.rm = TRUE)) { + as.numeric(coef(lm(AFDW ~ WW, data = cur_data()[inlier_WW_AFDW, , drop = FALSE]))[2]) + } else { + NA_real_ + } + ) %>% + ungroup() %>% + group_by(soort) %>% + mutate( + # Calculate AFDW/n + ratio_AFDW_n = AFDW / n, + + # Calculate species-level quantiles and IQR for AFDW/n + Q1_ratio_n = as.numeric(quantile(ratio_AFDW_n, 0.25, na.rm = TRUE)), + Q3_ratio_n = as.numeric(quantile(ratio_AFDW_n, 0.75, na.rm = TRUE)), + IQR_ratio_n = Q3_ratio_n - Q1_ratio_n, + lower_bound_ratio_n = Q1_ratio_n - 1.5 * IQR_ratio_n, + upper_bound_ratio_n = Q3_ratio_n + 1.5 * IQR_ratio_n, + + # Define outliers for AFDW/n ratio + outlier_AFDW_n = ifelse( + ratio_AFDW_n < lower_bound_ratio_n | ratio_AFDW_n > upper_bound_ratio_n, + 1,0), + inlier_AFDW_n = !is.na(ratio_AFDW_n) & + ratio_AFDW_n >= lower_bound_ratio_n & + ratio_AFDW_n <= upper_bound_ratio_n) %>% + + mutate( + slopeAFDW_n = if (any(inlier_AFDW_n, na.rm = TRUE)) { + as.numeric(coef(lm(AFDW ~ n, data = cur_data()[inlier_AFDW_n, , drop = FALSE]))[2]) + } else { + NA_real_ + } + ) %>% + ungroup() %>% + dplyr::mutate(AFDWcor.1 = dplyr::if_else(outlier_AFDW_WW == "Outlier" & outlier_AFDW_n == 1 & !is.na(WW) & WW !=0 & slopeAFDW_WW >0.1, WW*slopeAFDW_WW, dplyr::if_else(AFDW == 0 | AFDW<0 | is.na(AFDW) & WW !=0, WW*slopeAFDW_WW, AFDW))) %>% + dplyr::mutate(AFDWcor.2 = dplyr::if_else(is.na(AFDWcor.1) & WW !=0 & slopeAFDW_WW >0.1, WW*slopeAFDW_WW, AFDWcor.1)) %>% + dplyr::mutate(AFDWcor.3 = dplyr::if_else(is.na(AFDWcor.2) & WW !=0 & is.na(slopeAFDW_WW), WW*0.17, AFDWcor.2)) %>% + dplyr::mutate(AFDWcor.4 = dplyr::if_else(is.na(AFDWcor.3) & slopeAFDW_n>0 & n>0, n*slopeAFDW_n, AFDWcor.3)) %>% +dplyr::mutate(AFDWcor.5 = dplyr::if_else(outlier_AFDW_WW == "Outlier" & soort == "Dicentrarchus labrax" & !is.na(WW) & !is.na(AFDW), WW*slopeAFDW_WW, AFDWcor.4)) %>% + dplyr::mutate(AFDWcor.7 = dplyr::if_else(WW/AFDWcor.5 < 2 & slopeAFDW_WW > 0.1 & !is.na(slopeAFDW_WW) & !is.na(WW), WW*slopeAFDW_WW, AFDWcor.5)) %>% + dplyr::mutate(AFDWcor.8 = dplyr::if_else(WW/AFDWcor.7 > 15 & slopeAFDW_WW > 0.1 & !is.na(WW) & !is.na(slopeAFDW_WW) & outlier_AFDW_n == 0, WW*slopeAFDW_WW, AFDWcor.7)) %>% + dplyr::mutate(AFDWcor.9 = dplyr::if_else(WW < AFDWcor.8 & !is.na(WW) & !is.na(slopeAFDW_WW) & slopeAFDW_WW > 0.1, WW*slopeAFDW_WW, AFDWcor.8)) + #relocate(c(AFDWcor.4, AFDWcor.5, AFDWcor.6, AFDWcor.7, AFDWcor.8, AFDWcor.9, soort, outlier_AFDW_WW, outlier_AFDW_n), .after = AFDW) + +telc.1 <- tel.1 %>% + dplyr::select(campagne, gebied, event_name, datum, fractie, hoger_taxon, soort, vis_NL, exoot, n, WW, DW, AW, AFDW = AFDWcor.9, materiaal, Invoerder, Aanmaakdatum, `Laatst gewijzigd`, Jaar, Maand) +``` + +```{r check 2: dubbels, vreemde 0-en etc} +#data die moeten w aangepast in hyperdatabase +tel.2 <- tel.1 %>% + dplyr::filter(AFDW != AFDWcor.9 | n == 0 | is.na(n)) %>% + dplyr::select(campagne, gebied, datum, fractie, soort, n, WW, DW, AW, AFDW, AFDWcorrected = AFDWcor.9, materiaal, Invoerder, Aanmaakdatum) + +file_name1 <- + paste0(pad_data, "Hyperbenthos_records_NeedRevision", ".xlsx") + + +write_xlsx(tel.2, + path = file_name1) ``` +```{r om data te checken, figuren maken v ratio ww/afdw per soort} +# Loop through each unique species and save a separate plot +soort_lijst <- unique(tel.1$soort) + +for (soort in soort_lijst) { + # Subset data for the current species + soort_data <- tel.1[tel.1$soort == soort, ] + + # Identify outliers using the 1.5 * IQR rule + soort_data$ratio <- soort_data$WW / soort_data$AFDWcor.9 + Q1 <- quantile(soort_data$ratio, 0.25, na.rm = TRUE) + Q3 <- quantile(soort_data$ratio, 0.75, na.rm = TRUE) + IQR <- Q3 - Q1 + lower_bound <- Q1 - 1.5 * IQR + upper_bound <- Q3 + 1.5 * IQR + soort_data$outlier <- ifelse(soort_data$ratio < lower_bound | soort_data$ratio > upper_bound, "Outlier", "Non-Outlier") + + # Create the plot + plot <- ggplot(soort_data, aes(x = WW, y = AFDWcor.9, color=outlier)) + + geom_point() + # Add points + scale_color_manual(values = c("Outlier" = "red", "Non-Outlier" = "blue")) + # Color mapping + geom_smooth(method = "lm", se = FALSE) + # Add linear smoother + stat_regline_equation(aes(label = ..eq.label..), label.x.npc = "left", label.y.npc = "top") + # Add equation + theme_minimal() + # Use a minimal theme + labs( + title = paste("WW vs AFDW for", soort), + x = "Wet Weight (WW)", + y = "Ash-Free Dry Weight (AFDW)" + ) + + # Save the plot to a JPG file + filename <- paste0("G:/.shortcut-targets-by-id/0B0xcP-eNvJ9dZDBwVVJOVk5Ld2s/PRJ_SCHELDE/VNSC/Rapportage_INBO/2024/080_hyperbenthos/data/ratio_figuren_ww_afdw", "/", soort, "_plotWW_AFDW.jpg") # Create a filename + ggsave(filename, plot = plot, width = 6, height = 4, dpi = 300) # Save the plot +} + +``` + + + data gecreëerd op `r Sys.time()` data weggeschreven naar `r paste0(pad_data, "template_data.csv")` @@ -115,7 +229,7 @@ data weggeschreven naar `r paste0(pad_data, "template_data.csv")` ```{r jaren} jaren <- - telc %>% + telc.1 %>% distinct(Jaar) %>% pull(Jaar) @@ -126,9 +240,10 @@ jaar_range <- ```{r wegschrijven-data} file_name <- - paste0(pad_data, "hyperbenthos_data_", paste(jaar_range, collapse = "_"), ".xlsx") + paste0(pad_data, "hyperbenthos_data_revised", paste(jaar_range, collapse = "_"), ".xlsx") + -write_xlsx(telc, +write_xlsx(telc.1, path = file_name) ``` diff --git a/moneos_2025/080_hyperbenthos/080_hyperbenthos_analyse.Rmd b/moneos_2025/080_hyperbenthos/080_hyperbenthos_analyse.Rmd new file mode 100644 index 0000000..77f479f --- /dev/null +++ b/moneos_2025/080_hyperbenthos/080_hyperbenthos_analyse.Rmd @@ -0,0 +1,2177 @@ +--- +params: + hoofdstuk: "080_hyperbenthos" +knit: (function(inputFile, ...) { + rmarkdown::render(inputFile, + output_dir = paste0(rmarkdown::yaml_front_matter(inputFile)$params$hoofdstuk, "/output"))}) +title: "Hyperbenthos data" +output: word_document +editor_options: + chunk_output_type: console +--- + +```{r 080-setup, include=FALSE} + +knitr::opts_chunk$set(echo = FALSE, error=FALSE, warning=FALSE, message=FALSE, cache=FALSE) + +``` + + +```{r 080-libraries} + +library(tidyverse) +library(readxl) +library(writexl) +library(ggpubr) +library(INBOtheme) +library(rprojroot) +library(lubridate) +library(ggforce) +library(grid) +library(forcats) +library(gghighlight) +library(scales) + +``` + + +```{r 080-pad} + +# inlezen van variabelen +# pad naar data : pad_data +# pad naar tabellen : pad_tabellen +# pad naar figuren : pad_figuren + +source(find_root_file("../pad.R", criterion = is_rstudio_project)) #../ + +pad_data <- maak_pad(params$hoofdstuk, "data") +pad_figuren <- maak_pad(params$hoofdstuk, "figuren") +pad_tabellen <- maak_pad(params$hoofdstuk, "tabellen") + +source("G:/Gedeelde drives/PRJ_SCHELDE/MONEOS_rapportage/2022/Function_CalculateShannonIndex.R") +``` + +# Dataframes voorbereiden +```{r 080-data} + +gebied_order <- + c("Paardenschor","St-Anna","Rupel", "Ballooi","Dendermonde","Brede Schoren") + +zone_order <- c("Sterke Saliniteitsgradiënt","Oligohalien", "Zoet") + +################################ +# VOORBEREIDING +################################ +#ruwe datafile met alles erin, selectie van jaren en maanden en echte hyperbenthische soorten (s.l.) + +data_hyperbenthos <- + read_excel(paste0(pad_data, "hyperbenthos_data_revised2013_2024.xlsx")) %>% + dplyr::filter(soort != "Potamocorbula amurensis") %>% # is geen hyperb, volgens EMSE moeten alle soorten behalve garnalen, aasgarnalen en steurgarnalen eruit, daarom hieronder ook een versie voor deze groepen alleen + dplyr::filter(Maand %in% c(4:10)) %>% # data bevatten nu nog alle maanden, maar monitoring is enkel voor april-oktober. Tem 2018 werden ook andere maanden gedaan. + dplyr::filter(Jaar > 2013) %>% #2013 was een try-out jaar niet-conform de methode en spreiding van stations en periodes + dplyr::mutate(zone = if_else(gebied %in% c("Brede Schoren","Dendermonde"), "Zoet", if_else(gebied %in% c("Ballooi","Rupel"), "Oligohalien", "Sterke Saliniteitsgradiënt"))) %>% + dplyr::mutate(gebied = factor(gebied, + levels = gebied_order)) %>% + dplyr::mutate(zone = factor(zone, + levels = zone_order)) + + +# voor EMSE is men strikter en worden enkel (aas, steur)garnalen meegerekend +data_hyperbenthos.EMSE <- data_hyperbenthos %>% + dplyr::filter(hoger_taxon %in% c("Decapoda", "Mysida")) %>% #(steur)garnalen & aasgarnalen + dplyr::filter(!str_detect(soort, "Carcinus|Eriocheir|Hemigrap")) #zonder de krabben want die zitten in ankerkuildata + +unique(data_hyperbenthos.EMSE$soort) #controle of er geen krabbensoorten meer tussen zitten + +################################ +# SOORTENRIJKDOM +################################ + +#cf EMSE moet (enkel) soortenrijkdom zonder exoten, daarom hiervoor aparte basis dframes +data_hpb.gnexoot <- data_hyperbenthos %>% #algemeen + dplyr::filter(exoot == 0) + +data_hpb.EMSE.gnexoot <- data_hyperbenthos.EMSE %>% #EMSE + dplyr::filter(exoot == 0) + + +################################ +# ABUNDANTIE & BIOMASSA +################################ + +# Som biomassa & abundantie alle soorten algemeen +data_hyperbenthos_totaal <- + data_hyperbenthos %>% + dplyr::group_by(Jaar, Maand, gebied) %>% + dplyr::summarise_at(vars(n, AFDW), ~sum(.,na.rm=TRUE)) %>% + dplyr::mutate(zone = if_else(gebied == "Brede Schoren" | gebied == "Dendermonde", "Zoet", if_else(gebied == "Ballooi" | gebied == "Rupel", "Oligohalien", "Sterke Saliniteitsgradiënt"))) %>% + ungroup() + +# Som biomassa & abundantie soorten EMSE +data_hyperbenthos_totaal.EMSE <- + data_hyperbenthos.EMSE %>% + complete(soort, gebied, Maand, Jaar, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar, Maand, gebied) %>% + dplyr::summarise_at(vars(n, AFDW), ~sum(.,na.rm=TRUE)) %>% + dplyr::mutate(zone = if_else(gebied == "Brede Schoren" | gebied == "Dendermonde", "Zoet", if_else(gebied == "Ballooi" | gebied == "Rupel", "Oligohalien", "Sterke Saliniteitsgradiënt"))) %>% + ungroup() + +################################ +# VARIA +################################ + +vroegste_jaar <- + data_hyperbenthos %>% + pull(Jaar) %>% + min() + +laatste_jaar <- + data_hyperbenthos %>% + pull(Jaar) %>% + max() + + + +``` + + +```{r 080- EXTRA totaal-NIET GEBRUIKT, eval=FALSE, include=FALSE} + +#som per soort, per jaar, ook exoten, en in percent. Stomme INBOtheme kan maar 9 kleuren, dus selecteren van 8 grootste +#eerst jaarsom om % te knn berekenen in volgende stap +hpb.totsom <- data_hyperbenthos %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(totbiom = sum(na.omit(AFDW)), + totdens = sum(na.omit(n))) + +data.h.b <- data_hyperbenthos %>% + left_join(hpb.totsom, by = "Jaar") %>% + dplyr::group_by(soort, Jaar) %>% + dplyr::summarise(biom_perc = sum(AFDW)/totbiom) %>% + ungroup() %>% + distinct() %>% + arrange(desc(biom_perc)) %>% + dplyr::group_by(Jaar) %>% + slice(1:8) %>% + ungroup() + +#om negende fractie ("rest") te krijgen +tussenstap.1 <- data.h.b %>% + group_by(Jaar) %>% + summarise(biom_perct = 1- sum(biom_perc)) %>% + ungroup() %>% + mutate(soort = "rest") + +#de hoogste 8 samenvoegen met de "rest" +data.hpb.biomsoort<- tussenstap.1 %>% + dplyr::select(soort, Jaar, biom_perct) %>% + rename(biom_perc = biom_perct) %>% + rbind(data.h.b) + + +## probleem blijft want je hebt over meerdere jaren 21 soorten, dus soorten krijgen versch kleuren in versch jaren. Daarom selectie van 9 meest abundante taxa maken, en die dan voor elk jaar plotten. + + + + +``` + +```{r 080-EXTRA dframe trends voor 6 algemeenste soorten} +#trends voor een selectie van 6 soorten die de bulk van biomassa en aantallen uitmaken + +data.h.bsp <- data_hyperbenthos %>% + dplyr::mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + dplyr::mutate(soort = recode(soort, "Pomatoschistus minutus" = "Pomatoschistus sp")) %>% + dplyr::filter(soort %in% c("Platichthys flesus", "Crangon crangon", "Neomysis integer", "Pomatoschistus sp", "Palaemon longirostris", "Mesopodopsis slabberi")) %>% + dplyr::group_by(soort, Jaar, zone) %>% + dplyr::summarise(biom_sp = sum(na.omit(AFDW))) %>% + ungroup() %>% + complete(soort, zone, Jaar, fill = list(biom_sp = 0)) #geen Palaemons in 2023!!! + +ggplot(data.h.bsp, aes(x=Jaar, y=biom_sp, color=soort))+ + geom_line(aes(color=soort), linewidth=2)+ + scale_x_continuous(labels = scales::number_format(accuracy = 1))+ + ylab("Biomassa AFDW jaarsom 2 locaties")+ + facet_grid(~zone) + +ggsave(paste0(pad_figuren, "080-figuur-biomassa_soorten_jaren.jpg"), height=6, width=9) + +levels(ordered(data_hyperbenthos$soort)) + +``` + + +```{r 080-EXTRA dataframe taxa - selectie 8 ifv pie-diagr} +# selecteer 8 taxa met hoogste proc jaarbijdrage aan biomassa overheen alle jaren. Dus als ze in 1 jaar extreem talrijk was, dan telt dat, ook al is ze gemiddeld overheen alle jaren niet bij de talrijkste. +# note: Keuze voor 8 omdat INBO_theme maar 9 kleur levels heeft. + + +hpb.sel.8sp<- data_hyperbenthos %>% + left_join(hpb.totsom, by = "Jaar") %>% + dplyr::group_by(soort, Jaar) %>% + dplyr::reframe(biom_perc = sum(na.omit(AFDW))/totbiom*100) %>% + ungroup() %>% + distinct() %>% + dplyr::group_by(soort) %>% + dplyr::summarize(perc_biom = max(na.omit(biom_perc))) %>% + arrange(desc(perc_biom)) %>% + slice(1:8) %>% + ungroup() + +sel <- hpb.sel.8sp$soort + +#dframe met alleen de 8 sp en hun perc_biom +hpb.sel.8sp.biom<- data_hyperbenthos %>% + left_join(hpb.totsom, by = "Jaar") %>% + dplyr::group_by(soort, Jaar) %>% + dplyr::reframe(biom_perc = sum(na.omit(AFDW))/totbiom*100) %>% + ungroup() %>% + complete(soort, Jaar, fill = list(biom_perc = 0)) %>% + dplyr::filter(soort %in% sel) %>% + distinct() + +# Nu berekenen wat biom proc is van rest voor elk jaar +restbiom8sp <-data_hyperbenthos %>% + dplyr::filter(soort %in% sel) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(biom_8sp = sum(na.omit(AFDW))) %>% + left_join(hpb.totsom, by ="Jaar") %>% + dplyr::mutate(perc_biom = (totbiom-biom_8sp)/totbiom*100, + soort = "rest") %>% + dplyr::select(soort, Jaar, perc_biom) %>% + rename(biom_perc = perc_biom) + +hpb.sel.8sp.compl <- hpb.sel.8sp.biom %>% + rbind(restbiom8sp) + +``` + + +```{r 080-per-waterloop-en-tidaal} + +data_hyperbenthos_ZSmaand <- + data_hyperbenthos_totaal %>% + group_by(Jaar, Maand) %>% + summarise_at(vars(n, AFDW), + list(mean = ~max(0, mean(., na.rm = TRUE)), + med = ~max(0, median(., na.rm = TRUE)), + lwr1 = ~max(0, quantile(., 0.25, na.rm = TRUE)), + upr1 = ~max(0, quantile(., 0.75, na.rm = TRUE)), + lwr2 = ~max(0, quantile(., 0.05, na.rm = TRUE)), + upr2 = ~max(0, quantile(., 0.95, na.rm = TRUE)))) %>% + ungroup() + +#jaartotalen voor ZS - som van 6 gebieden, alleen april-okt, +data_hyperbenthos_ZS <- + data_hyperbenthos_totaal %>% + dplyr::group_by(Jaar) %>% + summarise_at(vars(n, AFDW), + list(tot = ~sum(., na.rm = TRUE), mean = ~max(0, mean(., na.rm = TRUE)), + med = ~max(0, median(., na.rm = TRUE)), + lwr1 = ~max(0, quantile(., 0.25, na.rm = TRUE)), + upr1 = ~max(0, quantile(., 0.75, na.rm = TRUE)), + lwr2 = ~max(0, quantile(., 0.05, na.rm = TRUE)), + upr2 = ~max(0, quantile(., 0.95, na.rm = TRUE)))) %>% + ungroup() +``` + +# Beschrijvende figuren trends alle stations: DENSITEIT +```{r 080-figuur-densiteit-totaal-gebied-maandverloop-jaren} + +## Alle hyperbenthos samen +ylb <- expression(paste("densiteit ", "(ind/", '40', m^3, ")")) + +fnt <- 8 + + +bxp_hpbaprokt <- + data_hyperbenthos_totaal %>% + mutate(gebied = factor(gebied, + levels = gebied_order), + Maand = ordered(Maand)) + +bxp_hpbaprokt24 <- + data_hyperbenthos_totaal %>% + mutate(gebied = factor(gebied, + levels = gebied_order), + Maand = ordered(Maand)) %>% + dplyr::filter(Jaar == 2024) + +bxp_hpb <- bxp_hpbaprokt %>% + ggplot(aes(x=Maand, y= n, group=Jaar)) + + geom_line(aes(colour=Jaar), linewidth = 0.8) + + geom_line(data = bxp_hpbaprokt24, aes(x=Maand, y= n), linewidth = 3, colour = "lightblue") + + #gghighlight(Jaar == 2023, + # unhighlighted_params = list(linewidth = 0.5, colour= NULL), keep_scales = TRUE)+ + #geom_line(data = dplyr::filter(bxp_hpbaprokt, Jaar == 2023), linewidth = 2)+ + scale_y_log10(breaks = c(0,10,1000,100000)+1, labels = c(0,10,1000,100000)) + + scale_color_continuous(label = function(x) sprintf("%.0f", x)) + + labs(x = "Maand", + y = ylb) + +bxp_hpb1 <- bxp_hpb+ + facet_grid_paginate(~gebied, ncol=3, nrow=1, page=1) + + +bxp_hpb2 <- bxp_hpb+ + facet_grid_paginate(~gebied, ncol=3, nrow=1, page=2) + + +ggarrange(bxp_hpb1 + rremove("xlab")+ font("xy.text", size = fnt), + bxp_hpb2 + font("xy.text", size = fnt), + nrow = 2, common.legend = TRUE, legend = "right") + +ggsave(paste0(pad_figuren, "080-figuur-densiteit_totaal-gebied-maandverloop_jaren.jpg"), height=6, width=9) + +############################################################### +# EMSE +############################################################### + +## Alle hyperbenthos samen + +bxp_hpbaprokt.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::mutate(gebied = factor(gebied, + levels = gebied_order), + Maand = ordered(Maand)) + +bxp_hpbaprokt24.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::mutate(gebied = factor(gebied, + levels = gebied_order), + Maand = ordered(Maand)) %>% + dplyr::filter(Jaar == 2024) + +bxp_hpb.EMSE <- bxp_hpbaprokt.EMSE %>% + ggplot(aes(x=Maand, y= n, group=Jaar)) + + geom_line(aes(colour=Jaar), linewidth = 0.8) + + geom_line(data = bxp_hpbaprokt24, aes(x=Maand, y= n), linewidth = 3, colour = "lightblue") + + #gghighlight(Jaar == 2023, + # unhighlighted_params = list(linewidth = 0.5, colour= NULL), keep_scales = TRUE)+ + #geom_line(data = dplyr::filter(bxp_hpbaprokt, Jaar == 2023), linewidth = 2)+ + scale_y_log10(breaks = c(0,10,1000,100000)+1, labels = c(0,10,1000,100000)) + + scale_color_continuous(label = function(x) sprintf("%.0f", x)) + + labs(x = "Maand", + y = ylb) + +bxp_hpb1 <- bxp_hpb.EMSE + + facet_grid_paginate(~gebied, ncol=3, nrow=1, page=1) + + +bxp_hpb2 <- bxp_hpb.EMSE + + facet_grid_paginate(~gebied, ncol=3, nrow=1, page=2) + + +ggarrange(bxp_hpb1 + rremove("xlab")+ font("xy.text", size = fnt), + bxp_hpb2 + font("xy.text", size = fnt), + nrow = 2, common.legend = TRUE, legend = "right") + +ggsave(paste0(pad_figuren, "080-figuur-densiteit_totaal-gebied-maandverloop_jaren.EMSEsoorten.jpg"), height=6, width=9) + +### VERHOUDING EMSE versus REST DOORHEEN DE JAREN ##### +####################################################### + +verhouding <- data_hyperbenthos %>% + dplyr::mutate(EMSE = ifelse(hoger_taxon %in% c("Mysida", "Decapoda"), "EMSE", "REST")) %>% # EMSE taxa (aas|steur)garnalen + dplyr::mutate(EMSE = ifelse(str_detect(soort, "Carcinus|Eriocheir|takanoi"), "REST", EMSE)) %>% #geen krabben + dplyr::select(zone, Jaar, soort, hoger_taxon, n, AFDW, EMSE) %>% + dplyr::group_by(zone, Jaar) %>% + dplyr::mutate(totbiom = sum(na.omit(AFDW)), # jaarsom per zone alle soorten + totN = sum(na.omit(n)), .groups = "drop") %>% + dplyr::group_by(zone, Jaar, EMSE) %>% + dplyr::summarise(tot.biom = sum(na.omit(AFDW)), # jaarsom per zone per EMSE/REST + tot.N = sum(na.omit(n)), .groups = "drop") + +EMSE.REST.fig.biom <- ggplot(verhouding, aes(x = Jaar, y = tot.biom, fill = EMSE)) + + geom_area() + + scale_fill_manual(values = c("EMSE" = "darkseagreen3", "REST" = "cadetblue")) + + ylab(" Jaarsom Biomassa per zone") + + scale_x_continuous(breaks=seq(2014,2024,2)) + + facet_grid(rows = vars(zone), scales = "free_y") + +EMSE.REST.fig.biom +ggsave(paste0(pad_figuren, "080-figuur-biomassa_EMSEvsREST.jpg"), height=6, width=9) + +EMSE.REST.fig.N <- ggplot(verhouding, aes(x = Jaar, y = tot.N, fill = EMSE)) + + geom_area() + + scale_fill_manual(values = c("EMSE" = "darkseagreen3", "REST" = "cadetblue")) + + ylab("Jaarsom abundantie per zone") + + scale_x_continuous(breaks=seq(2014,2024,2)) + + facet_grid(rows = vars(zone), scales = "free_y") + +EMSE.REST.fig.N +ggsave(paste0(pad_figuren, "080-figuur-densiteit_EMSEvsREST.jpg"), height=6, width=9) + +``` + +# Beschrijvende figuren trends EMSE soorten: BIOMASSA + +```{r 80-figuur-EMSE-soorten-trends-zone } + +## Avg Biomassa per zone doorheen de jaren +Crangon <- data_hyperbenthos.EMSE %>% + complete(soort, Maand, Jaar, gebied, zone, fill = list(n = 0, AFDW = 0)) %>% + dplyr::filter(soort == "Crangon crangon") %>% + dplyr::group_by(zone, Jaar) %>% + dplyr::summarise(N.avg = mean(n), + Biom.avg = mean(AFDW)) %>% + dplyr::ungroup() %>% + dplyr::mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) + + +crangon.fig <- ggplot(Crangon, aes(x = Jaar, y = Biom.avg, fill = zone)) + + geom_area() + + scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) + + annotate("text", x = 2019, y=1.4, label = expression(italic("Crangon crangon")), size = unit(4.5, "pt")) + + scale_fill_manual(values = c("Sterke Saliniteitsgradiënt" = "darkseagreen3", "Oligohalien" = "cadetblue", "Zoet" = "darkturquoise")) + +labs(x = "", y="") + + theme(plot.margin = unit(c(0,0.2,0,1), 'lines')) +crangon.fig +######### + +Neomysis.integer <- data_hyperbenthos.EMSE %>% + complete(soort, Maand, Jaar, gebied, zone, fill = list(n = 0, AFDW = 0)) %>% + dplyr::filter(soort == "Neomysis integer") %>% + dplyr::group_by(zone, Jaar) %>% + dplyr::summarise(N.avg = mean(n), + Biom.avg = mean(AFDW)) %>% + dplyr::ungroup() + + +Neomysis.integer.fig <- ggplot(Neomysis.integer, aes(x = Jaar, y = Biom.avg, fill = zone)) + + geom_area() + + scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) + + annotate("text", x = 2019, y=0.23, label = expression(italic("Neomysis integer")), size = unit(4.5, "pt")) + + scale_fill_manual(values = c("Sterke Saliniteitsgradiënt" = "darkseagreen3", "Oligohalien" = "cadetblue", "Zoet" = "darkturquoise"))+ +labs(x = "", y="") + + theme(plot.margin = unit(c(0,0.2,0,1), 'lines')) +Neomysis.integer.fig + +####### + +Mesopodopsis <- data_hyperbenthos.EMSE %>% + complete(soort, Maand, Jaar, gebied, zone, fill = list(n = 0, AFDW = 0)) %>% + dplyr::filter(soort == "Mesopodopsis slabberi") %>% + dplyr::group_by(zone, Jaar) %>% + dplyr::summarise(N.avg = mean(n), + Biom.avg = mean(AFDW)) %>% + dplyr::ungroup() + + +Mesopodopsis.fig <- ggplot(Mesopodopsis, aes(x = Jaar, y = Biom.avg, fill = zone)) + + geom_area() + + scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) + + annotate("text", x = 2019, y=1, label = expression(italic("Mesopodopsis slabberi")), size = unit(4.5, "pt")) + + scale_fill_manual(values = c("Sterke Saliniteitsgradiënt" = "darkseagreen3", "Oligohalien" = "cadetblue", "Zoet" = "darkturquoise"))+ +labs(x = "", y="") + + theme(plot.margin = unit(c(0,0.2,0,1), 'lines')) +Mesopodopsis.fig + +####### +Palaemon <- data_hyperbenthos.EMSE %>% + complete(soort, Maand, Jaar, gebied, zone, fill = list(n = 0, AFDW = 0)) %>% + dplyr::filter(soort == "Palaemon longirostris") %>% + dplyr::group_by(zone, Jaar) %>% + dplyr::summarise(N.avg = mean(n), + Biom.avg = mean(AFDW)) %>% + dplyr::ungroup() + + +Palaemon.fig <- ggplot(Palaemon, aes(x = Jaar, y = Biom.avg, fill = zone)) + + geom_area() + + scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) + + annotate("text", x = 2019, y=2.2, label = expression(italic("Palaemon longirostris")), size = unit(4.5, "pt")) + + scale_fill_manual(values = c("Sterke Saliniteitsgradiënt" = "darkseagreen3", "Oligohalien" = "cadetblue", "Zoet" = "darkturquoise")) + +labs(x = "", y="") + + theme(plot.margin = unit(c(0.1,0.2,0,1), 'lines')) +Palaemon.fig + +t <- ggarrange(crangon.fig + font("xy.text", size = fnt), + Palaemon.fig + font("xy.text", size = fnt), + Mesopodopsis.fig + font("xy.text", size = fnt), + Neomysis.integer.fig + font("xy.text", size = fnt), ncol=1, nrow = 4) +annotate_figure(t, left = text_grob("Gemiddelde biomassa per jaar per zone", rot = 90, size = 16), bottom = text_grob("Jaar", size = 16)) + +ggsave(paste0(pad_figuren, "080-EMSE-soorten-biomassa-zone-jaren.jpg"), height = 9, width = 9) + + +############################## +# Zelfde voor biomassa +crangon.fig.n <- ggplot(Crangon, aes(x = Jaar, y = N.avg, fill = zone)) + + geom_area() + + scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) + + annotate("text", x = 2019, y=650, label = expression(italic("Crangon crangon")), size = unit(4.5, "pt")) + + scale_fill_manual(values = c("Sterke Saliniteitsgradiënt" = "darkseagreen3", "Oligohalien" = "cadetblue", "Zoet" = "darkturquoise"))+ +labs(x = "", y="") + + theme(plot.margin = unit(c(0,0.2,0,1), 'lines')) +crangon.fig.n + +Palaemon.fig.n <- ggplot(Palaemon, aes(x = Jaar, y = N.avg, fill = zone)) + + geom_area() + + scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) + + annotate("text", x = 2019, y=68, label = expression(italic("Palaemon longirostris")), size = unit(4.5, "pt")) + + scale_fill_manual(values = c("Sterke Saliniteitsgradiënt" = "darkseagreen3", "Oligohalien" = "cadetblue", "Zoet" = "darkturquoise")) + +labs(x = "", y="") + + theme(plot.margin = unit(c(0,0.2,0,1), 'lines')) +Palaemon.fig.n + +Mesopodopsis.fig.n <- ggplot(Mesopodopsis, aes(x = Jaar, y = N.avg, fill = zone)) + + geom_area() + + scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) + + annotate("text", x = 2019, y=2100, label = expression(italic("Mesopodopsis slabberi")), size = unit(4.5, "pt")) + + scale_fill_manual(values = c("Sterke Saliniteitsgradiënt" = "darkseagreen3", "Oligohalien" = "cadetblue", "Zoet" = "darkturquoise"))+ +labs(x = "", y="") + + theme(plot.margin = unit(c(0,0.2,0,1), 'lines')) +Mesopodopsis.fig.n + +Neomysis.integer.fig.n <- ggplot(Neomysis.integer, aes(x = Jaar, y = N.avg, fill = zone)) + + geom_area() + + scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) + + annotate("text", x = 2019, y=175, label = expression(italic("Neomysis integer")), size = unit(4.5, "pt")) + + scale_fill_manual(values = c("Sterke Saliniteitsgradiënt" = "darkseagreen3", "Oligohalien" = "cadetblue", "Zoet" = "darkturquoise"))+ +labs(x = "", y="") + + theme(plot.margin = unit(c(0,0.2,0,1), 'lines')) +Neomysis.integer.fig.n + +t.n <- ggarrange(crangon.fig.n + font("xy.text", size = fnt), + Palaemon.fig.n + font("xy.text", size = fnt), + Mesopodopsis.fig.n + font("xy.text", size = fnt), + Neomysis.integer.fig.n + font("xy.text", size = fnt), ncol=1, nrow = 4) +annotate_figure(t.n, left = text_grob("Gemiddelde abundantie per jaar per zone", rot = 90, size=16), bottom = text_grob("Jaar", size = 16)) + +ggsave(paste0(pad_figuren, "080-EMSE-soorten-abundantie-zone-jaren.jpg"), height = 9, width = 9) + +``` + + +```{r 080-figuur-densiteit-totaal-ZS-maandverloop-jaren} +ylb <- expression(paste("densiteit ", "(ind/", '40', m^3, ")")) + +fnt <- 8 + +bxp_hpbZSaprokt <- + data_hyperbenthos_totaal %>% + mutate(Maand = ordered(Maand)) %>% + dplyr::group_by(Jaar, Maand) %>% + dplyr::summarise(nZS = sum(na.omit(n)), + AFDWZS = sum(na.omit(AFDW))) + +bxp_hpbZSaprokt24 <- bxp_hpbZSaprokt %>% + dplyr::filter(Jaar == 2024) + +bxp_hpbZS <- bxp_hpbZSaprokt %>% + ggplot(aes(x = Maand, y = nZS, group = Jaar)) + + geom_line(aes(colour=Jaar), linewidth = 0.8) + + geom_line(data= bxp_hpbZSaprokt24, aes(x = Maand, y = nZS), linewidth = 3, colour= "lightblue") + + scale_y_log10(breaks = c(0,10,1000,10000) + 1, labels = c(0,10,1000,10000)) + + labs(x = "Maand", + y = ylb) +bxp_hpbZS + +ggsave(paste0(pad_figuren, "080-figuur-densiteit_totaal-ZS-maandverloop_jaren.jpg"), height = 4, width = 6) + +################################################################# +##deze versie bijgewerkt met ribbon en hline +bxp_hpbZSj_ALL <- + data_hyperbenthos_totaal %>% + group_by(Jaar, Maand) %>% # stations en zones samen, zodat enkel maandelijks verschillen overblijven. + dplyr::summarise(n=mean(n), + AFDW =mean(AFDW), .groups = "drop") %>% + dplyr::group_by(Jaar) %>% + summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)), + med = ~max(0, median(., na.rm = TRUE)), + lwr1 = ~max(0, quantile(., 0.25, na.rm = TRUE)), + upr1 = ~max(0, quantile(., 0.75, na.rm = TRUE)), + lwr2 = ~max(0, quantile(., 0.05, na.rm = TRUE)), + upr2 = ~max(0, quantile(., 0.95, na.rm = TRUE)) + )), .groups = "drop") %>% + ggplot(aes(x=Jaar, y= AFDW_mean)) + + geom_ribbon(aes(ymin = n_lwr2, ymax = n_upr2), alpha = 0.3) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDW_mean, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(0,6000)+ ggtitle("Volledige Zeeschelde")+ + theme(title =element_text(size=10)) + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) +data_summary <- data_hyperbenthos_totaal %>% + dplyr::group_by(Jaar, Maand) %>% + dplyr::summarise(n = mean(n), AFDW = mean(AFDW), .groups = "drop") %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)), + med = ~max(0, median(., na.rm = TRUE)), + lwr1 = ~max(0, quantile(., 0.25, na.rm = TRUE)), + upr1 = ~max(0, quantile(., 0.75, na.rm = TRUE)), + lwr2 = ~max(0, quantile(., 0.05, na.rm = TRUE)), + upr2 = ~max(0, quantile(., 0.95, na.rm = TRUE)) + )), .groups = "drop") + +# Compute 3-year rollmean of AFDW_mean +data_summary$rollmean <- zoo::rollmean(data_summary$AFDW_mean, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016 <- data_summary$rollmean[data_summary$Jaar == 2016] + +bxp_hpbZSj_ALL + + geom_hline(yintercept = rollmean_2016, linetype = "dashed", colour = "red") + +bxp_hpbZSj_SS <- + data_hyperbenthos_totaal %>% + dplyr::filter(zone == "Sterke Saliniteitsgradiënt") %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(nZS = sum(na.omit(n)/2), + AFDWZS = sum(na.omit(AFDW))/2) %>% + dplyr::mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + ggplot(aes(x=Jaar, y= nZS)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(nZS, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(100,80000)+ ggtitle("Sterke Saliniteitsgradiënt")+ + theme(title =element_text(size=10)) + +bxp_hpbZSj_OL <- + data_hyperbenthos_totaal %>% + dplyr::filter(zone == "Oligohalien") %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(nZS = sum(na.omit(n)/2), + AFDWZS = sum(na.omit(AFDW))/2) %>% + mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + ggplot(aes(x=Jaar, y= nZS)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(nZS, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(100,3000)+ ggtitle("Oligohalien") + + theme(title =element_text(size=10)) + +bxp_hpbZSj_ZO <- + data_hyperbenthos_totaal %>% + dplyr::filter(zone == "Zoet") %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(nZS = sum(na.omit(n))/2, + AFDWZS = sum(na.omit(AFDW))/2) %>% + mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + ggplot(aes(x = Jaar, y = nZS)) + + geom_line(size = 1.2) + + geom_line(aes(y = zoo::rollmean(nZS, 3, na.pad = TRUE)), size = 1.2, colour="red") + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45)) + + ylim(100, 3000) + + ggtitle("Zoet") + + theme(title = element_text(size=10)) + +t <-ggarrange(bxp_hpbZSj_SS + font("xy.text", size = fnt), + bxp_hpbZSj_OL + font("xy.text", size = fnt), + bxp_hpbZSj_ZO + font("xy.text", size = fnt), + bxp_hpbZSj_ALL + font("xy.text", size = fnt), ncol=2, nrow = 2) +annotate_figure(t, left = text_grob("Gemiddelde densiteit per jaar per zone", rot=90), bottom = text_grob("Jaar")) + + +ggsave(paste0(pad_figuren, "080-figuur-densiteit_ZS_jaarverloop_zones1.jpg"), height = 4, width = 6) + + + +bxp_hpbZSb <- + data_hyperbenthos_totaal %>% + mutate(Maand = ordered(Maand)) %>% + dplyr::group_by(Jaar, Maand) %>% + dplyr::summarise(nZS = sum(na.omit(n)), + AFDWZS = sum(na.omit(AFDW))) %>% + ggplot(aes(x = Jaar, y = nZS, group=Maand)) + + geom_line(aes(colour = Maand), size = 0.8) + + scale_y_log10(breaks = c(0,10,1000,10000) + 1, labels = c(0,10,1000,10000)) + + labs(x = "Maand", + y = ylb) +bxp_hpbZSb + + +``` + +#TOETSPARAMETER abundantie EMSE +```{r 080 - TOETSPARAMETER ABUNDANTIE} +fnt <- 8 + +trend_ABUN_ZS.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::group_by(Jaar, Maand) %>% # stations en zones samen, zodat enkel maandelijks verschillen overblijven. + dplyr::summarise(n=mean(n), + AFDW =mean(AFDW), .groups = "drop") %>% + complete(Jaar, Maand, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar) %>% + summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") %>% + ggplot(aes(x=Jaar, y= n_mean)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(n_mean, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(0,4000)+ ggtitle("Volledige Zeeschelde")+ + theme(title = element_text(size=10)) +trend_ABUN_ZS.EMSE + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) +data_summary.EMSE <- data_hyperbenthos_totaal.EMSE %>% + dplyr::group_by(Jaar, Maand) %>% + dplyr::summarise(n = mean(n), AFDW = mean(AFDW), .groups = "drop") %>% + complete(Jaar, Maand, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") + +## For illustration add a ribbon of variation based on 1.96se of simple smoothed gam +library(mgcv) + +gam_model <- mgcv::gam(n_mean ~ s(Jaar, k = 5), data = data_summary.EMSE) + +# Predict & add predicted values & CI +newdata <- data.frame(Jaar = data_summary.EMSE$Jaar) +pred <- predict(gam_model, newdata = newdata, se.fit = TRUE) +newdata$fit <- pred$fit +newdata$se <- pred$se.fit +newdata$upper <- newdata$fit + 1.96 * newdata$se +newdata$lower <- newdata$fit - 1.96 * newdata$se + +# Filter to only positive fitted values for plotting purposes +newdata_pos <- newdata %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(n_mean = fit) + +# Compute 3-year rollmean of AFDW_mean +data_summary.EMSE$rollmean <- zoo::rollmean(data_summary.EMSE$n_mean, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.EMSE <- data_summary.EMSE$rollmean[data_summary.EMSE$Jaar == 2016] + +ABUND.ZS <- trend_ABUN_ZS.EMSE + + geom_hline(yintercept = rollmean_2016.EMSE, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_pos, aes(x=Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +ABUND.ZS + + + +######## Saliniteitsgradiënt ########## + +trend_ABUN_SS.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Sterke Saliniteitsgradiënt") %>% + dplyr::mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% #zowel variatie ts stations als ts maanden wordt meegenomen + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") %>% + ggplot(aes(x=Jaar, y= n_mean)) + + geom_point(size = 1.2) + + geom_line(aes(y=zoo::rollmean(n_mean, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(0,10000)+ ggtitle("Saliniteitsgradiënt")+ + theme(title = element_text(size=10)) + +### Add threshold based on rollmean, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) +data_summary.SS.EMSE <- data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Sterke Saliniteitsgradiënt") %>% + mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") + +##### Add some indication of variation of data, purely illustrative. Based on simple GAM +gam_modelSS <- mgcv::gam(n_mean ~ s(Jaar, k = 5), data = data_summary.SS.EMSE) + +# Predict & add predicted values & CI +newdataSS <- data.frame(Jaar = data_summary.SS.EMSE$Jaar) +predSS <- predict(gam_modelSS, newdata = newdataSS, se.fit = TRUE) +newdataSS$fit <- predSS$fit +newdataSS$se <- predSS$se.fit +newdataSS$upper <- newdataSS$fit + 1.96 * newdataSS$se +newdataSS$lower <- newdataSS$fit - 1.96 * newdataSS$se + +# Filter to only positive fitted values for plotting purposes +newdata_posSS <- newdataSS %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(n_mean = fit) + +##### Compute 3-year rollmean of AFDW_mean +data_summary.SS.EMSE$rollmean <- zoo::rollmean(data_summary.SS.EMSE$n_mean, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.SS.EMSE <- data_summary.SS.EMSE$rollmean[data_summary.SS.EMSE$Jaar == 2016] + +ABUND.SS <- trend_ABUN_SS.EMSE + + geom_hline(yintercept = rollmean_2016.SS.EMSE, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_posSS, aes(x=Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +ABUND.SS + + +######## OLIGOHALIEN ########## + +trend_ABUN_OLIG.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Oligohalien") %>% + dplyr::mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% #zowel variatie ts stations als ts maanden wordt meegenomen + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") %>% + ggplot(aes(x=Jaar, y= n_mean)) + + geom_point(size = 1.2) + + geom_line(aes(y=zoo::rollmean(n_mean, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(0,300)+ ggtitle("Oligohalien")+ + theme(title = element_text(size=10)) +trend_ABUN_OLIG.EMSE + +### Add threshold based on rollmean, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) +data_summary.OLIG.EMSE <- data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Oligohalien") %>% + mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") + +# Compute 3-year rollmean of AFDW_mean +data_summary.OLIG.EMSE$rollmean <- zoo::rollmean(data_summary.OLIG.EMSE$n_mean, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.OLIG.EMSE <- data_summary.OLIG.EMSE$rollmean[data_summary.OLIG.EMSE$Jaar == 2016] + +##### Add some indication of variation of data, purely illustrative. Based on simple GAM +gam_modelOLIG <- mgcv::gam(n_mean ~ s(Jaar, k = 5), data = data_summary.OLIG.EMSE) + +# Predict & add predicted values & CI +newdataOLIG <- data.frame(Jaar = data_summary.OLIG.EMSE$Jaar) +predOLIG <- predict(gam_modelOLIG, newdata = newdataOLIG, se.fit = TRUE) +newdataOLIG$fit <- predOLIG$fit +newdataOLIG$se <- predOLIG$se.fit +newdataOLIG$upper <- newdataOLIG$fit + 1.96 * newdataOLIG$se +newdataOLIG$lower <- newdataOLIG$fit - 1.96 * newdataOLIG$se + +# Filter to only positive fitted values for plotting purposes +newdata_posOLIG <- newdataOLIG %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(n_mean = fit) + +# FIGUUR +ABUND.OLIG <- trend_ABUN_OLIG.EMSE + + geom_hline(yintercept = rollmean_2016.OLIG.EMSE, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_posOLIG, aes(x=Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +ABUND.OLIG + + +######## ZOET ########## + +trend_ABUN_ZOET.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Zoet") %>% + dplyr::mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% #zowel variatie ts stations als ts maanden wordt meegenomen + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") %>% + ggplot(aes(x=Jaar, y= n_mean)) + + geom_point(size = 1.2) + + geom_line(aes(y=zoo::rollmean(n_mean, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(0,300)+ ggtitle("Zoet")+ + theme(title = element_text(size=10)) +trend_ABUN_ZOET.EMSE + +### Add threshold based on rollmean, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) +data_summary.ZOET.EMSE <- data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Zoet") %>% + mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") + +# Compute 3-year rollmean of AFDW_mean +data_summary.ZOET.EMSE$rollmean <- zoo::rollmean(data_summary.ZOET.EMSE$n_mean, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.ZOET.EMSE <- data_summary.ZOET.EMSE$rollmean[data_summary.ZOET.EMSE$Jaar == 2016] + +##### Add some indication of variation of data, purely illustrative. Based on simple GAM +gam_modelZOET <- mgcv::gam(n_mean ~ s(Jaar, k = 5), data = data_summary.ZOET.EMSE) + +# Predict & add predicted values & CI +newdataZOET <- data.frame(Jaar = data_summary.ZOET.EMSE$Jaar) +predZOET <- predict(gam_modelZOET, newdata = newdataZOET, se.fit = TRUE) +newdataZOET$fit <- predZOET$fit +newdataZOET$se <- predZOET$se.fit +newdataZOET$upper <- newdataZOET$fit + 1.96 * newdataZOET$se +newdataZOET$lower <- newdataZOET$fit - 1.96 * newdataZOET$se + +# Filter to only positive fitted values for plotting purposes +newdata_posZOET <- newdataZOET %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(n_mean = fit) + + +# FIGUUR + +ABUND.ZOET <- trend_ABUN_ZOET.EMSE + + geom_hline(yintercept = rollmean_2016.ZOET.EMSE, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_posZOET, aes(x=Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +ABUND.ZOET + + +####FIGUREN SAMEN + +t.EMSE <- ggarrange(ABUND.ZS + font("xy.text", size = fnt), + ABUND.SS + font("xy.text", size = fnt), + ABUND.OLIG + font("xy.text", size = fnt), + ABUND.ZOET + font("xy.text", size = fnt), ncol=2, nrow = 2, + font.label = list(face = "plain")) +annotate_figure(t.EMSE, left = text_grob(" Gemiddelde densiteit per bemonstering per zone", rot=90), bottom = text_grob("Jaar")) + +ggsave(paste0(pad_figuren, "080-figuur-densiteit_zones.EMSE2.jpg"), height=4, width=6) +``` + + +```{r 080-figuur-biomassa-totaal-ZS-maandverloop-jaren} +ylbb <- expression(paste("biomassa ", "(g droge stof/", '40', m^3, ")")) + +fnt <- 8 + +biom_hpbZSOA <- + data_hyperbenthos_totaal %>% + mutate(Maand = ordered(Maand)) %>% + dplyr::group_by(Jaar, Maand) %>% + dplyr::summarise(nZS = sum(na.omit(n)), + AFDWZS = sum(na.omit(AFDW))) %>% + dplyr::ungroup() + +biom_hpbZSOA24 <- biom_hpbZSOA %>% + dplyr::filter(Jaar == 2024) + + +biom_hpbZS <- biom_hpbZSOA %>% + ggplot(aes(x=Maand, y= AFDWZS, group=Jaar)) + + geom_line(aes(colour=Jaar), size = 0.8) + + geom_line(data= biom_hpbZSOA24, aes(x=Maand, y= AFDWZS), linewidth = 3, colour= "lightblue") + + scale_y_log10(breaks = c(0,10,100, 1000,10000)+1, labels = c(0,10,100,1000,10000)) + + labs(x = "Maand", + y = ylbb) + + scale_color_continuous(label = function(x) sprintf("%.0f", x)) + + theme(axis.text.x = element_text(angle = 45)) +biom_hpbZS + +ggsave(paste0(pad_figuren, "080-figuur-biomassa_totaal-ZS-maandverloop_jaren.jpg"), height=4, width=6) + +biom_hpbZSj_ALL <- + data_hyperbenthos_totaal %>% + dplyr::filter(Maand %in% c(4:10), Jaar !=2013) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(nZS = sum(na.omit(n)/6), + AFDWZS = sum(na.omit(AFDW))/6) %>% + ggplot(aes(x=Jaar, y= AFDWZS)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDWZS, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="")+ + ylim(0,60)+ ggtitle("Gehele Zeeschelde")+ + theme(title =element_text(size=10)) +biom_hpbZSj_ALL + +biom_hpbZSj_SS <- + data_hyperbenthos_totaal %>% + dplyr::filter(Maand %in% c(4:10), Jaar !=2013, zone == "Sterke Saliniteitsgradiënt") %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(nZS = sum(na.omit(n)/2), + AFDWZS = sum(na.omit(AFDW))/2) %>% + mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + ggplot(aes(x=Jaar, y= AFDWZS)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDWZS, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="")+ + ylim(0,60)+ ggtitle("Sterke Saliniteitsgradiënt")+ + theme(title =element_text(size=10)) + +biom_hpbZSj_OL <- + data_hyperbenthos_totaal %>% + dplyr::filter(Maand %in% c(4:10), Jaar !=2013, zone == "Oligohalien") %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(nZS = sum(na.omit(n)/2), + AFDWZS = sum(na.omit(AFDW))/2) %>% + mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + ggplot(aes(x=Jaar, y= AFDWZS)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDWZS, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="")+ + ylim(0,60)+ ggtitle("Oligohalien")+ + theme(title =element_text(size=10)) + +biom_hpbZSj_ZO <- + data_hyperbenthos_totaal %>% + dplyr::filter(Maand %in% c(4:10), Jaar !=2013, zone == "Zoet") %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(nZS = sum(na.omit(n)/2), + AFDWZS = sum(na.omit(AFDW))/2) %>% + mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + ggplot(aes(x=Jaar, y= AFDWZS)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDWZS, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="")+ + ylim(0,60)+ ggtitle("Zoet")+ + theme(title =element_text(size=10)) + +q <-ggarrange(biom_hpbZSj_SS + font("xy.text", size = fnt), + biom_hpbZSj_OL + font("xy.text", size = fnt), + biom_hpbZSj_ZO + font("xy.text", size = fnt), + biom_hpbZSj_ALL + font("xy.text", size = fnt), ncol=2, nrow = 2) +annotate_figure(q, left=text_grob("Gemiddelde biomassa (2 locaties)", rot=90), bottom=text_grob("Jaar")) + + + +ggsave(paste0(pad_figuren, "080-figuur-biomassa_ZS_jaarverloop_zones.jpg"), height=4, width=6) + +``` + +#TOETSPARAMETER Biomassa EMSE + +```{r 080 - TOETSPARAMETER BIOMASSA} +trend_BIOM_ZS.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::group_by(Jaar, Maand) %>% # stations en zones samen, zodat enkel maandelijks verschillen overblijven. + dplyr::summarise(n = mean(n), + AFDW = mean(AFDW), .groups = "drop") %>% + complete(Jaar, Maand, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar) %>% + summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") %>% + ggplot(aes(x=Jaar, y= AFDW_mean)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDW_mean, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(0,5)+ + ggtitle("Volledige Zeeschelde") + + theme(title = element_text(size=10)) + +## For illustration add a ribbon of variation based on 1.96se of simple smoothed gam +library(mgcv) + +gam_model <- mgcv::gam(AFDW_mean ~ s(Jaar, k = 5), data = data_summary.EMSE) + +# Predict & add predicted values & CI +newdata <- data.frame(Jaar = data_summary.EMSE$Jaar) +pred <- predict(gam_model, newdata = newdata, se.fit = TRUE) +newdata$fit <- pred$fit +newdata$se <- pred$se.fit +newdata$upper <- newdata$fit + 1.96 * newdata$se +newdata$lower <- newdata$fit - 1.96 * newdata$se + +# Filter to only positive fitted values for plotting purposes +newdata_pos <- newdata %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(AFDW_mean = fit) + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) +data_summary.EMSE <- data_hyperbenthos_totaal.EMSE %>% + dplyr::group_by(Jaar, Maand) %>% + dplyr::summarise(n = mean(n), AFDW = mean(AFDW), .groups = "drop") %>% + complete(Jaar, Maand, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") + +# Compute 3-year rollmean of AFDW_mean +data_summary.EMSE$rollmean <- zoo::rollmean(data_summary.EMSE$AFDW_mean, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.EMSE <- data_summary.EMSE$rollmean[data_summary.EMSE$Jaar == 2016] + +BIOM.ZS <-trend_BIOM_ZS.EMSE + + geom_hline(yintercept = rollmean_2016.EMSE, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_pos, aes(x=Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +BIOM.ZS + + + +######## Saliniteitsgradiënt ########## + +trend_BIOM_SS.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Sterke Saliniteitsgradiënt") %>% + dplyr::mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% #zowel variatie ts stations als ts maanden wordt meegenomen + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") %>% + ggplot(aes(x=Jaar, y= AFDW_mean)) + + geom_point(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDW_mean, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(0,8)+ ggtitle("Saliniteitsgradiënt")+ + theme(title = element_text(size=10)) + +### Add threshold based on rollmean, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) +data_summary.SS.EMSE <- data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Sterke Saliniteitsgradiënt") %>% + mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") + +##### Add some indication of variation of data, purely illustrative. Based on simple GAM +gam_modelSS <- mgcv::gam(AFDW_mean ~ s(Jaar, k = 5), data = data_summary.SS.EMSE) + +# Predict & add predicted values & CI +newdataSS <- data.frame(Jaar = data_summary.SS.EMSE$Jaar) +predSS <- predict(gam_modelSS, newdata = newdataSS, se.fit = TRUE) +newdataSS$fit <- predSS$fit +newdataSS$se <- predSS$se.fit +newdataSS$upper <- newdataSS$fit + 1.96 * newdataSS$se +newdataSS$lower <- newdataSS$fit - 1.96 * newdataSS$se + +# Filter to only positive fitted values for plotting purposes +newdata_posSS <- newdataSS %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(AFDW_mean = fit) + +##### Compute 3-year rollmean of AFDW_mean +data_summary.SS.EMSE$rollmean <- zoo::rollmean(data_summary.SS.EMSE$AFDW_mean, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.SS.EMSE <- data_summary.SS.EMSE$rollmean[data_summary.SS.EMSE$Jaar == 2016] + +BIOM.SS <- trend_BIOM_SS.EMSE + + geom_hline(yintercept = rollmean_2016.SS.EMSE, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_posSS, aes(x=Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) + +BIOM.SS + +######## OLIGOHALIEN ########## + +trend_BIOM_OLIG.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Oligohalien") %>% + dplyr::mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% #zowel variatie ts stations als ts maanden wordt meegenomen + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") %>% + ggplot(aes(x=Jaar, y= AFDW_mean)) + + geom_point(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDW_mean, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(0,5)+ ggtitle("Oligohalien")+ + theme(title = element_text(size=10)) +trend_BIOM_OLIG.EMSE + +### Add threshold based on rollmean, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) +data_summary.OLIG.EMSE <- data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Oligohalien") %>% + mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") + +# Compute 3-year rollmean of AFDW_mean +data_summary.OLIG.EMSE$rollmean <- zoo::rollmean(data_summary.OLIG.EMSE$AFDW_mean, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.OLIG.EMSE <- data_summary.OLIG.EMSE$rollmean[data_summary.OLIG.EMSE$Jaar == 2016] + +##### Add some indication of variation of data, purely illustrative. Based on simple GAM +gam_modelOLIG <- mgcv::gam(AFDW_mean ~ s(Jaar, k = 5), data = data_summary.OLIG.EMSE) + +# Predict & add predicted values & CI +newdataOLIG <- data.frame(Jaar = data_summary.OLIG.EMSE$Jaar) +predOLIG <- predict(gam_modelOLIG, newdata = newdataOLIG, se.fit = TRUE) +newdataOLIG$fit <- predOLIG$fit +newdataOLIG$se <- predOLIG$se.fit +newdataOLIG$upper <- newdataOLIG$fit + 1.96 * newdataOLIG$se +newdataOLIG$lower <- newdataOLIG$fit - 1.96 * newdataOLIG$se + +# Filter to only positive fitted values for plotting purposes +newdata_posOLIG <- newdataOLIG %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(AFDW_mean = fit) + + +# FIGUUR + +BIOM.OLIG <- trend_BIOM_OLIG.EMSE + + geom_hline(yintercept = rollmean_2016.OLIG.EMSE, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_posOLIG, aes(x=Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +BIOM.OLIG + + +######## ZOET ########## + +trend_BIOM_ZOET.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Zoet") %>% + dplyr::mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% #zowel variatie ts stations als ts maanden wordt meegenomen + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") %>% + ggplot(aes(x=Jaar, y= AFDW_mean)) + + geom_point(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDW_mean, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(0,5)+ ggtitle("Zoet")+ + theme(title = element_text(size=10)) +trend_BIOM_ZOET.EMSE + +### Add threshold based on rollmean, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) +data_summary.ZOET.EMSE <- data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Zoet") %>% + mutate(gebied = droplevels(factor(gebied))) %>% + complete(Jaar, Maand, gebied, fill = list(n = 0, AFDW = 0)) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(across(c(n, AFDW), list( + tot = ~sum(., na.rm = TRUE), + mean = ~max(0, mean(., na.rm = TRUE)))), .groups = "drop") + +# Compute 3-year rollmean of AFDW_mean +data_summary.ZOET.EMSE$rollmean <- zoo::rollmean(data_summary.ZOET.EMSE$AFDW_mean, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.ZOET.EMSE <- data_summary.ZOET.EMSE$rollmean[data_summary.ZOET.EMSE$Jaar == 2016] + +##### Add some indication of variation of data, purely illustrative. Based on simple GAM +gam_modelZOET <- mgcv::gam(AFDW_mean ~ s(Jaar, k = 5), data = data_summary.ZOET.EMSE) + +# Predict & add predicted values & CI +newdataZOET <- data.frame(Jaar = data_summary.ZOET.EMSE$Jaar) +predZOET <- predict(gam_modelZOET, newdata = newdataZOET, se.fit = TRUE) +newdataZOET$fit <- predZOET$fit +newdataZOET$se <- predZOET$se.fit +newdataZOET$upper <- newdataZOET$fit + 1.96 * newdataZOET$se +newdataZOET$lower <- newdataZOET$fit - 1.96 * newdataZOET$se + +# Filter to only positive fitted values for plotting purposes +newdata_posZOET <- newdataZOET %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(AFDW_mean = fit) + + +# FIGUUR + +BIOM.ZOET <- trend_BIOM_ZOET.EMSE + + geom_hline(yintercept = rollmean_2016.ZOET.EMSE, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_posZOET, aes(x=Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +BIOM.ZOET + + + +####FIGUREN SAMEN + +t.EMSE.BIOM <- ggarrange(BIOM.ZS + font("xy.text", size = fnt), + BIOM.SS + font("xy.text", size = fnt), + BIOM.OLIG + font("xy.text", size = fnt), + BIOM.ZOET + font("xy.text", size = fnt), ncol = 2, nrow = 2) +annotate_figure(t.EMSE.BIOM, left=text_grob("Gemiddelde biomassa (g AFDW) per 2bemonstering per zone", rot = 90, size = 9), bottom = text_grob("Jaar", size = 12)) + + +ggsave(paste0(pad_figuren, "080-figuur-biomassa_zones.EMSE2.jpg"), height=4, width=6) +``` + + + + + + + +```{r TOETSPARAMETER EMSE Biomassa OUDE VERSIE, eval=FALSE, include=FALSE} +biom_hpbZSj_ALL.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(nZS = sum(na.omit(n)/6), + AFDWZS = sum(na.omit(AFDW))/6) %>% + ggplot(aes(x=Jaar, y= AFDWZS)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDWZS, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + ylim(0,50) + + ggtitle("Gehele Zeeschelde") + + theme(title =element_text(size=10)) +biom_hpbZSj_ALL.EMSE + +biom_hpbZSj_SS.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Sterke Saliniteitsgradiënt") %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(nZS = sum(na.omit(n)/2), + AFDWZS = sum(na.omit(AFDW))/2) %>% + mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + ggplot(aes(x=Jaar, y= AFDWZS)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDWZS, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="")+ + ylim(0,50)+ ggtitle("Sterke Saliniteitsgradiënt")+ + theme(title =element_text(size=10)) + +biom_hpbZSj_SS.EMSE + +biom_hpbZSj_OL.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Oligohalien") %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(nZS = sum(na.omit(n)/2), + AFDWZS = sum(na.omit(AFDW))/2) %>% + mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + ggplot(aes(x=Jaar, y= AFDWZS)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDWZS, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="")+ + ylim(0,50)+ ggtitle("Oligohalien")+ + theme(title =element_text(size=10)) + +biom_hpbZSj_OL.EMSE + +biom_hpbZSj_ZO.EMSE <- + data_hyperbenthos_totaal.EMSE %>% + dplyr::filter(zone == "Zoet") %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(nZS = sum(na.omit(n)/2), + AFDWZS = sum(na.omit(AFDW))/2) %>% + mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + ggplot(aes(x=Jaar, y= AFDWZS)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(AFDWZS, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="")+ + ylim(0,50)+ ggtitle("Zoet")+ + theme(title =element_text(size=10)) + +biom_hpbZSj_ZO.EMSE + +q <-ggarrange(biom_hpbZSj_SS.EMSE + font("xy.text", size = fnt), + biom_hpbZSj_OL.EMSE + font("xy.text", size = fnt), + biom_hpbZSj_ZO.EMSE + font("xy.text", size = fnt), + biom_hpbZSj_ALL.EMSE + font("xy.text", size = fnt), ncol=2, nrow = 2) +annotate_figure(q, left=text_grob("Gemiddelde biomassa (2 locaties)", rot=90), bottom=text_grob("Jaar")) + + + +ggsave(paste0(pad_figuren, "080-figuur-biomassa_ZS_jaarverloop_zones.EMSE.jpg"), height=4, width=6) + +``` + + + +```{r 080-figuur-biomassa-gebied-maandverloop-jaar, eval=FALSE, include=FALSE} + + +ylbb <- expression(paste(" biomassa ", "(mg droge stof/", '40', m^3, ")")) + +fnt <- 8 + + +bxp_hpbbOA <- + data_hyperbenthos_totaal %>% + dplyr::filter(Maand %in% c(4:10)) %>% + mutate(gebied = factor(gebied, + levels = gebied_order), + Maand = ordered(Maand)) +bxp_hpbbOA23 <- bxp_hpbbOA %>% + dplyr::filter(Jaar == 2023) + +bxp_hpbb <- bxp_hpbbOA %>% + ggplot(aes(x=Maand, y= AFDW+1, group=Jaar)) + + geom_line(aes(colour=Jaar)) + + geom_line(data= bxp_hpbbOA23, aes(x=Maand, y= AFDW+1), linewidth = 2, colour= "lightblue") + + scale_y_log10(breaks = c(0,10,1000,100000)+1, labels = c(0,10,1000,100000)) + + labs(x = "Maand", + y = ylbb) + + scale_color_continuous(label = function(x) sprintf("%.0f", x)) + + theme(axis.text.x = element_text(angle = 45)) + +bxp_hpbb1 <- bxp_hpbb+ + labs(x = "Maand", + y = "") + + facet_grid_paginate(~gebied, ncol=3, nrow=1, page=1) + +bxp_hpbb2 <- bxp_hpbb+ + facet_grid_paginate(~gebied, ncol=3, nrow=1, page=2) + + +ggarrange(bxp_hpbb1 + rremove("xlab")+ font("xy.text", size = fnt), + bxp_hpbb2 + font("xy.text", size = fnt), + nrow = 2, common.legend = TRUE, legend = "right") + +ggsave(paste0(pad_figuren, "080_figuur_biomassa_gebied_maandverloop_jaar.jpg"), height=6, width=9) + + +``` + + +```{r 080-figuur-biomassa-biomassa-perc-soorten-perjaar, eval=FALSE, include=FALSE} + +pie_up <- hpb.sel.8sp.compl %>% + ggplot(aes(x="", y=biom_perc, fill=soort)) + + geom_bar(stat="identity", width=1, color="white") + + coord_polar("y", start=0) + + theme_void() + + facet_grid_paginate(~Jaar, ncol=3, nrow=1, page=1) +pie_middle<- hpb.sel.8sp.compl %>% + ggplot(aes(x="", y=biom_perc, fill=soort)) + + geom_bar(stat="identity", width=1, color="white") + + coord_polar("y", start=0) + + theme_void() + + facet_grid_paginate(~Jaar, ncol=3, nrow=1, page=2) +pie_bottom<- hpb.sel.8sp.compl %>% + ggplot(aes(x="", y=biom_perc, fill=soort)) + + geom_bar(stat="identity", width=1, color="white") + + coord_polar("y", start=0) + + theme_void() + + facet_grid_paginate(~Jaar, ncol=3, nrow=1, page=3) +pie_last<- hpb.sel.8sp.compl %>% + ggplot(aes(x="", y=biom_perc, fill=soort)) + + geom_bar(stat="identity", width=1, color="white") + + coord_polar("y", start=0) + + theme_void() + + facet_grid_paginate(~Jaar, ncol=3, nrow=1, page=4) + +ggarrange(pie_up, pie_middle, pie_bottom, pie_last, + nrow = 4, common.legend = TRUE, legend ="right") + +ggsave(paste0(pad_figuren, "080_figuur_biomassa_biomassa_perc_soorten_perjaar.jpg"), height=4, width=6) + +``` + + +#TOETSPARAMETER SOORTENRIJKDOM +```{r 080-figuur-soortenrijkdomZS-perjaar} +#met en zonder exoten ter vgl + + +hpb_specrich.ZS <- data_hyperbenthos %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(specrich = length(unique(soort))) %>% + ungroup() %>% + dplyr::mutate(zone = "Zeeschelde") %>% + dplyr::select(Jaar, zone, specrich) + +hpb_specrich <- data_hyperbenthos %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(specrich = length(unique(soort))) %>% + ungroup() %>% + bind_rows(hpb_specrich.ZS) %>% + dplyr::mutate(zone = + factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet", "Zeeschelde"))) + + +welex <-ggplot(hpb_specrich, aes(x = Jaar, y = specrich, color = zone)) + + geom_line(size=1) + + ylab("Taxa rijkdom") + + theme_bw() + + theme(axis.title=element_text(size=15)) + + ylim(13,60) + + theme(legend.position = "none") +welex + +###### ZONDER EXOTEN ##### +########################## + +hpb_specrichgnex.ZS <- data_hyperbenthos %>% + dplyr::filter(exoot == 0) %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(specrich = length(unique(soort))) %>% + ungroup() %>% + dplyr::mutate(zone = "Zeeschelde") %>% + dplyr::select(Jaar, zone, specrich) + +hpb_specrichgnex <- data_hyperbenthos %>% + dplyr::filter(exoot == 0) %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(specrich = length(unique(soort))) %>% + ungroup() %>% + bind_rows(hpb_specrichgnex.ZS) %>% + dplyr::mutate(zone = + factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet", "Zeeschelde"))) + +gnex <- ggplot(hpb_specrichgnex, aes(x = Jaar, y = specrich, color = zone)) + + geom_line(size=1) + + ylab("Taxa rijkdom zonder exoten") + + ylim(13,60) + + theme_bw() + + theme(axis.title=element_text(size=15)) + + theme(legend.text = element_text(size=12)) + + theme(legend.key.size = unit(0.44, "cm")) + + theme(legend.position = "inside", + legend.position.inside = c(0.7,0.85)) + + +gnex +ggarrange(welex, gnex, + nrow = 1) + +ggsave(paste0(pad_figuren, "080_figuur_soortenrijkdomZS_perjaar.jpg"), height = 7, width = 10) + +################## met avg soortenrijkdom ########### +####################################################### + +# dframe maken met spec richness +hpb_specrichgnex.avgzs <- data_hyperbenthos %>% + dplyr::filter(exoot == 0) %>% + dplyr::group_by(Jaar, Maand) %>% + dplyr::summarise(specrich = length(unique(soort)), .groups = "drop") %>% + dplyr::mutate(zone = "Zeeschelde") %>% + dplyr::select(Jaar, Maand, zone, specrich) + +hpb_specrichgnex.avg.month <- data_hyperbenthos %>% + dplyr::filter(exoot == 0) %>% + dplyr::group_by(Jaar, Maand, zone) %>% + dplyr::summarise(specrich = length(unique(soort)), .groups = "drop") %>% + bind_rows(hpb_specrichgnex.avgzs) +hpb_specrichgnex.avg <- hpb_specrichgnex.avg.month %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(Nsp.avg = mean(specrich), .groups = "drop") + +# dframes voor elke zone +Sp.rich.ZS <- hpb_specrichgnex.avg %>% + dplyr::filter(zone == "Zeeschelde") +Sp.rich.SS <- hpb_specrichgnex.avg %>% + dplyr::filter(zone == "Sterke Saliniteitsgradiënt") +Sp.rich.ZOET <- hpb_specrichgnex.avg %>% + dplyr::filter(zone == "Zoet") +Sp.rich.OLIG <- hpb_specrichgnex.avg %>% + dplyr::filter(zone == "Oligohalien") + + +Sp.rich.ZS.month <- hpb_specrichgnex.avg.month %>% + dplyr::filter(zone == "Zeeschelde") +Sp.rich.SS.month <- hpb_specrichgnex.avg.month %>% + dplyr::filter(zone == "Sterke Saliniteitsgradiënt") +Sp.rich.ZOET.month <- hpb_specrichgnex.avg.month %>% + dplyr::filter(zone == "Zoet") +Sp.rich.OLIG.month <- hpb_specrichgnex.avg.month %>% + dplyr::filter(zone == "Oligohalien") + + + +######### figuren voor elke zone ########### +library(mgcv) + + +## Sterke Saliniteitsgradiënt +brk <- function(x) seq(ceiling(x[1]), floor(x[2]), by = 1) + +Navg.SS <- ggplot(data = Sp.rich.SS, aes(x = Jaar, y = Nsp.avg)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(Nsp.avg, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ggtitle("Sterke Saliniteitsgradiënt") + + theme(title = element_text(size=10)) + + scale_y_continuous(breaks=seq(4,15,2), limits = c(4, 15)) +Navg.SS +## For illustration add a ribbon of variation based on 1.96se of simple smoothed gam + +gam_model <- mgcv::gam(specrich ~ s(Jaar, k = 5), data = Sp.rich.SS.month) + +# Predict & add predicted values & CI +newdata <- data.frame(Jaar = Sp.rich.SS$Jaar) +pred <- predict(gam_model, newdata = newdata, se.fit = TRUE) +newdata$fit <- pred$fit +newdata$se <- pred$se.fit +newdata$upper <- newdata$fit + 1.96 * newdata$se +newdata$lower <- newdata$fit - 1.96 * newdata$se + +# Filter to only positive fitted values for plotting purposes +newdata_posSSrich <- newdata %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(Nsp.avg = fit) + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) + +# Compute 3-year rollmean of AFDW_mean +Sp.rich.SS$rollmean <- zoo::rollmean(Sp.rich.SS$Nsp.avg, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.EMSE.SS <- Sp.rich.SS$rollmean[Sp.rich.SS$Jaar == 2016] + +Navg.SS.x <- Navg.SS + + geom_hline(yintercept = rollmean_2016.EMSE.SS, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_posSSrich, aes(x = Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +Navg.SS.x + +## ZEESCHELDE + +Navg.ZS <- ggplot(data = Sp.rich.ZS, aes(x = Jaar, y = Nsp.avg)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(Nsp.avg, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(4,20)+ + ggtitle("Zeeschelde") + + theme(title = element_text(size=10)) +Navg.ZS +## For illustration add a ribbon of variation based on 1.96se of simple smoothed gam + +gam_modelZS <- mgcv::gam(Nsp.avg ~ s(Jaar, k = 5), data = Sp.rich.ZS) + +# Predict & add predicted values & CI +newdataZS <- data.frame(Jaar = Sp.rich.ZS$Jaar) +predZS <- predict(gam_modelZS, newdata = newdataZS, se.fit = TRUE) +newdataZS$fit <- predZS$fit +newdataZS$se <- predZS$se.fit +newdataZS$upper <- newdataZS$fit + 1.96 * newdataZS$se +newdataZS$lower <- newdataZS$fit - 1.96 * newdataZS$se + +# Filter to only positive fitted values for plotting purposes +newdata_posZSrich <- newdataZS %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(Nsp.avg = fit) + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) + +# Compute 3-year rollmean of AFDW_mean +Sp.rich.ZS$rollmean <- zoo::rollmean(Sp.rich.ZS$Nsp.avg, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.EMSE.ZS <- Sp.rich.ZS$rollmean[Sp.rich.ZS$Jaar == 2016] + +Navg.ZS.x <- Navg.ZS + + geom_hline(yintercept = rollmean_2016.EMSE.ZS, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_posZSrich, aes(x = Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +Navg.ZS.x + +## OLIGOHALIEN + +Navg.OLIG <- ggplot(data = Sp.rich.OLIG, aes(x = Jaar, y = Nsp.avg)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(Nsp.avg, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(4,10)+ + ggtitle("Oligohalien") + + theme(title = element_text(size=10)) +Navg.OLIG +## For illustration add a ribbon of variation based on 1.96se of simple smoothed gam + +gam_modelOLIG <- mgcv::gam(Nsp.avg ~ s(Jaar, k = 5), data = Sp.rich.OLIG) + +# Predict & add predicted values & CI +newdataOLIG <- data.frame(Jaar = Sp.rich.OLIG$Jaar) +predOLIG <- predict(gam_modelOLIG, newdata = newdataOLIG, se.fit = TRUE) +newdataOLIG$fit <- predOLIG$fit +newdataOLIG$se <- predOLIG$se.fit +newdataOLIG$upper <- newdataOLIG$fit + 1.96 * newdataOLIG$se +newdataOLIG$lower <- newdataOLIG$fit - 1.96 * newdataOLIG$se + +# Filter to only positive fitted values for plotting purposes +newdata_posOLIGrich <- newdataOLIG %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(Nsp.avg = fit) + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) + +# Compute 3-year rollmean of AFDW_mean +Sp.rich.OLIG$rollmean <- zoo::rollmean(Sp.rich.OLIG$Nsp.avg, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.EMSE.OLIG <- Sp.rich.OLIG$rollmean[Sp.rich.OLIG$Jaar == 2016] + +Navg.OLIG.x <- Navg.OLIG + + geom_hline(yintercept = rollmean_2016.EMSE.OLIG, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_posOLIGrich, aes(x = Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +Navg.OLIG.x + +## ZOET + +Navg.ZOET <- ggplot(data = Sp.rich.ZOET, aes(x = Jaar, y = Nsp.avg)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(Nsp.avg, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ylim(4,10)+ + ggtitle("Zoet") + + theme(title = element_text(size=10)) +Navg.ZOET +## For illustration add a ribbon of variation based on 1.96se of simple smoothed gam + +gam_modelZOET <- mgcv::gam(Nsp.avg ~ s(Jaar, k = 5), data = Sp.rich.ZOET) + +# Predict & add predicted values & CI +newdataZOET <- data.frame(Jaar = Sp.rich.ZOET$Jaar) +predZOET <- predict(gam_modelZOET, newdata = newdataZOET, se.fit = TRUE) +newdataZOET$fit <- predZOET$fit +newdataZOET$se <- predZOET$se.fit +newdataZOET$upper <- newdataZOET$fit + 1.96 * newdataZOET$se +newdataZOET$lower <- newdataZOET$fit - 1.96 * newdataZOET$se + +# Filter to only positive fitted values for plotting purposes +newdata_posZOETrich <- newdataZOET %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(Nsp.avg = fit) + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) + +# Compute 3-year rollmean of AFDW_mean +Sp.rich.ZOET$rollmean <- zoo::rollmean(Sp.rich.ZOET$Nsp.avg, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.EMSE.ZOET <- Sp.rich.ZOET$rollmean[Sp.rich.ZOET$Jaar == 2016] + +Navg.ZOET.x <- Navg.ZOET + + geom_hline(yintercept = rollmean_2016.EMSE.ZOET, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_posZOETrich, aes(x = Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +Navg.ZOET.x + + +t.EMSE.SPECRICH <-ggarrange(Navg.ZS.x + font("xy.text", size = fnt), + Navg.SS.x + font("xy.text", size = fnt), + Navg.OLIG.x + font("xy.text", size = fnt), + Navg.ZOET.x + font("xy.text", size = fnt), ncol=2, nrow = 2) +annotate_figure(t.EMSE.SPECRICH, left = text_grob("Gemiddelde soortenrijkdom per maand per zone", rot=90), bottom=text_grob("Jaar")) + +ggsave(paste0(pad_figuren, "080-figuur-Soortenrijkdom.per.maand_ZONES2.jpg"), height = 4, width = 6) + + +########################################################## +#### Statistiek #### +# rijkdom verschilt weinig per maand, daarom als herhaalde meting gebruiken voor bepalen richness per jaar (niet echt juist, kweet het...) +ggplot(Sp.rich.ZS.month, aes(x=Maand, y = specrich))+ + geom_bar(stat = "summary", fun.y = "mean") + +## testen van de flating mean tov de referentie (ook een floating mean, maar dan rond 2016) +## test met lmer met als vb Oligohalien: +# Required packages +library(zoo) +library(lme4) +library(emmeans) + + +# Step 1: Sort and calculate 3-year floating mean per plot (Maand) +df_ma <- Sp.rich.OLIG.month %>% + arrange(Maand, Jaar) %>% + dplyr::group_by(Maand) %>% + dplyr::mutate( + specrich_ma3 = rollapply(specrich, width = 3, FUN = mean, align = "center", fill = NA), + center_year = rollapply(Jaar, width = 3, FUN = function(x) x[2], align = "center", fill = NA) + ) %>% + dplyr::ungroup() %>% + dplyr::filter(!is.na(specrich_ma3)) + +# Step 2: Create year_window factor and set reference (e.g., center year = 2001) +df_ma <- df_ma %>% + dplyr::mutate( + year_window = factor(center_year), + year_window = relevel(year_window, ref = "2016") # set 2016 as reference + ) + +# Step 3: Fit linear mixed model +model <- lmer(specrich_ma3 ~ year_window + (1 | Maand), data = df_ma) + +# Step 4: Estimated marginal means and comparison to reference window +em <- emmeans(model, ~ year_window) +results <- contrast(em, method = "trt.vs.ctrl", ref = which(levels(df_ma$year_window) == "2016")) + +# View results +summary(results) # voor floating mean van 2023 is p-waarde 0.12 voor contrast met floating mean van 2016 (referentie) + + +``` + +```{r 080-figuur-soortenrijkdomZS-perjaargebied} +#In jaar 2013 niet volledig gesampled wat taxa rijkdom beinvloedt, dus die niet +#met exoten +hpb_specrichzon <- data_hyperbenthos %>% + dplyr::mutate(zone = if_else(gebied == "Brede Schoren" | gebied == "Dendermonde", "Zoet", if_else(gebied == "Ballooi" | gebied == "Rupel", "Oligohalien", "Sterke Saliniteitsgradiënt"))) %>% + mutate(zone = factor(zone, + levels = c("Sterke Saliniteitsgradiënt", "Oligohalien", "Zoet"))) %>% + dplyr::filter(Maand %in% c(4:10)) %>% + dplyr::filter(Jaar != 2013) %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(specrich = length(unique(soort))) %>% + ungroup() + +welexgeb <-ggplot(hpb_specrichzon, aes(Jaar, specrich)) + + geom_line(size=1)+ + ylab("Taxa rijkdom")+ + ylim(0,50)+ + facet_grid(~zone)+ + theme(axis.text.x = element_text(angle = 45))+ + theme(strip.text = element_text(size=14)) +welexgeb + +ggsave(paste0(pad_figuren, "080_figuur_soortenrijkdom_perjaarzone.OUD.jpg"), height=5, width=10) +``` + +#TOETSPARAMETER Shannon diversity +```{r 080-figuur-Shannon-diversiteit} +shannon_hpbZS <- data_hyperbenthos %>% + dplyr::group_by(soort, Jaar) %>% + dplyr::summarise(nZS = sum(na.omit(n)), + afdwZS = sum(na.omit(AFDW)), .groups = "drop") %>% + dplyr::group_by(Jaar) %>% + dplyr::summarise(shannon_afdw = calc_shannon_index(afdwZS), + shannon_n = calc_shannon_index(nZS), .groups = "drop") %>% + dplyr::mutate(zone = "Zeeschelde") + +## Datafile Shannon: Zeeschelde + de 3 evaluatie zones +shannon_hpb <- data_hyperbenthos %>% + dplyr::group_by(soort, zone, Jaar) %>% + dplyr::summarise(nzone = sum(na.omit(n)), + afdwzone = sum(na.omit(AFDW)), .groups = "drop") %>% + dplyr::group_by(zone, Jaar) %>% + dplyr::summarise(shannon_afdw = calc_shannon_index(afdwzone), + shannon_n = calc_shannon_index(nzone), , .groups = "drop") %>% + rbind(shannon_hpbZS) + + + +### Figuren +shan_n <- ggplot(shannon_hpb, aes(x=Jaar, colour=zone))+ + geom_line(aes(y=shannon_n), size=1)+ + ylab("Shannon diversiteit densiteit")+ + xlab("Jaar") + + theme_bw() + + theme(legend.text = element_text(size=14))+ + theme(legend.title = element_text(size=14))+ + scale_x_continuous(labels = scales::number_format(accuracy = 1)) + + +shan_b <- ggplot(shannon_hpb, aes(x=Jaar, colour=zone))+ + geom_line(aes(y=shannon_afdw), size=1)+ + ylab("Shannon diversiteit biomassa")+ + xlab("Jaar") + + theme_bw() + + theme(legend.text = element_text(size=14))+ + theme(legend.title = element_text(size=14))+ + scale_x_continuous(labels = scales::number_format(accuracy = 1)) + +ggarrange(shan_n, shan_b, + nrow = 1, common.legend=T) + +ggsave(paste0(pad_figuren, "080_figuur_shannon.jpg"), height=5, width=10) + + + +#################################################### +# met gemiddeldes/maand*zone + +### Eerst dframes makes + +# voor Zeeschelde +shannon_hpbZS.avg <- data_hyperbenthos %>% + dplyr::group_by(soort, Jaar, Maand) %>% + dplyr::summarise(nZS = sum(na.omit(n)), + afdwZS = sum(na.omit(AFDW)), .groups = "drop") %>% + dplyr::group_by(Jaar, Maand) %>% + dplyr::summarise(shannon_afdw = calc_shannon_index(afdwZS), + shannon_n = calc_shannon_index(nZS), .groups = "drop") %>% + dplyr::mutate(zone = "Zeeschelde") + +## voor Zeeschelde + 3 zones +shannon_hpb.avg <- data_hyperbenthos %>% + dplyr::group_by(soort, zone, Jaar, Maand) %>% + dplyr::summarise(nzone = sum(na.omit(n)), + afdwzone = sum(na.omit(AFDW)), .groups = "drop") %>% + dplyr::group_by(zone, Jaar, Maand) %>% + dplyr::summarise(shannon_afdw = calc_shannon_index(afdwzone), + shannon_n = calc_shannon_index(nzone), , .groups = "drop") %>% + rbind(shannon_hpbZS.avg) %>% + dplyr::group_by(Jaar, zone) %>% + dplyr::summarise(Shannon.avg = mean(shannon_n), .groups = "drop") + +Shannon.ZS <- shannon_hpb.avg %>% + dplyr::filter(zone == "Zeeschelde") +Shannon.SS <- shannon_hpb.avg %>% + dplyr::filter(zone == "Sterke Saliniteitsgradiënt") +Shannon.OLIG <- shannon_hpb.avg %>% + dplyr::filter(zone == "Oligohalien") +Shannon.ZOET <- shannon_hpb.avg %>% + dplyr::filter(zone == "Zoet") + +########################################### +### Figuur voor Zeeschelde #### + +Shan.avg.ZS <- ggplot(data = Shannon.ZS, aes(x = Jaar, y = Shannon.avg)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(Shannon.avg, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ggtitle("Zeeschelde") + + theme(title = element_text(size=10)) + + scale_y_continuous(breaks=seq(0,1.8,0.2), limits = c(0, 1.8)) +Shan.avg.ZS +## For illustration add a ribbon of variation based on 1.96se of simple smoothed gam + +gam_model <- mgcv::gam(Shannon.avg ~ s(Jaar, k = 5), data = Shannon.ZS) + +# Predict & add predicted values & CI +newdata.sh.ZS <- data.frame(Jaar = Shannon.ZS$Jaar) +pred.sh.ZS <- predict(gam_model, newdata.sh.ZS = newdata.sh.ZS, se.fit = TRUE) +newdata.sh.ZS$fit <- pred.sh.ZS$fit +newdata.sh.ZS$se <- pred.sh.ZS$se.fit +newdata.sh.ZS$upper <- newdata.sh.ZS$fit + 1.96 * newdata.sh.ZS$se +newdata.sh.ZS$lower <- newdata.sh.ZS$fit - 1.96 * newdata.sh.ZS$se + +# Filter to only positive fitted values for plotting purposes +newdata_ribbon.sh.ZS <- newdata.sh.ZS %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(Shannon.avg = fit) + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) + +# Compute 3-year rollmean of AFDW_mean +Shannon.ZS$rollmean <- zoo::rollmean(Shannon.ZS$Shannon.avg, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.sh.ZS <- Shannon.ZS$rollmean[Shannon.ZS$Jaar == 2016] + +Shan.avg.ZS.x <- Shan.avg.ZS + + geom_hline(yintercept = rollmean_2016.sh.ZS, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_ribbon.sh.ZS, aes(x = Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +Shan.avg.ZS.x + +### Figuur voor Sterke Saliniteitsgradiënt #### + +Shan.avg.SS <- ggplot(data = Shannon.SS, aes(x = Jaar, y = Shannon.avg)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(Shannon.avg, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ggtitle("Sterke Saliniteitsgradiënt") + + theme(title = element_text(size=10)) + + scale_y_continuous(breaks=seq(0,1.8,0.2), limits = c(0, 1.8)) +Shan.avg.SS +## For illustration add a ribbon of variation based on 1.96se of simple smoothed gam + +gam_model <- mgcv::gam(Shannon.avg ~ s(Jaar, k = 5), data = Shannon.SS) + +# Predict & add predicted values & CI +newdata.sh.SS <- data.frame(Jaar = Shannon.SS$Jaar) +pred.sh.SS <- predict(gam_model, newdata.sh.SS = newdata.sh.SS, se.fit = TRUE) +newdata.sh.SS$fit <- pred.sh.SS$fit +newdata.sh.SS$se <- pred.sh.SS$se.fit +newdata.sh.SS$upper <- newdata.sh.SS$fit + 1.96 * newdata.sh.SS$se +newdata.sh.SS$lower <- newdata.sh.SS$fit - 1.96 * newdata.sh.SS$se + +# Filter to only positive fitted values for plotting purposes +newdata_ribbon.sh.SS <- newdata.sh.SS %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(Shannon.avg = fit) + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) + +# Compute 3-year rollmean of AFDW_mean +Shannon.SS$rollmean <- zoo::rollmean(Shannon.SS$Shannon.avg, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.sh.SS <- Shannon.SS$rollmean[Shannon.SS$Jaar == 2016] + +Shan.avg.SS.x <- Shan.avg.SS + + geom_hline(yintercept = rollmean_2016.sh.SS, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_ribbon.sh.SS, aes(x = Jaar, ymin = lower, ymax = upper), alpha = 0.2)+ + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +Shan.avg.SS.x + +### Figuur voor Oligohalien #### + +Shan.avg.OLIG <- ggplot(data = Shannon.OLIG, aes(x = Jaar, y = Shannon.avg)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(Shannon.avg, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ggtitle("Oligohalien") + + theme(title = element_text(size=10)) + + scale_y_continuous(breaks=seq(0.4,1.8,0.2), limits = c(0.4, 1.8)) +Shan.avg.OLIG +## For illustration add a ribbon of variation based on 1.96se of simple smoothed gam + +gam_model <- mgcv::gam(Shannon.avg ~ s(Jaar, k = 5), data = Shannon.OLIG) + +# Predict & add predicted values & CI +newdata.sh.OLIG <- data.frame(Jaar = Shannon.OLIG$Jaar) +pred.sh.OLIG <- predict(gam_model, newdata.sh.OLIG = newdata.sh.OLIG, se.fit = TRUE) +newdata.sh.OLIG$fit <- pred.sh.OLIG$fit +newdata.sh.OLIG$se <- pred.sh.OLIG$se.fit +newdata.sh.OLIG$upper <- newdata.sh.OLIG$fit + 1.96 * newdata.sh.OLIG$se +newdata.sh.OLIG$lower <- newdata.sh.OLIG$fit - 1.96 * newdata.sh.OLIG$se + +# Filter to only positive fitted values for plotting purposes +newdata_ribbon.sh.OLIG <- newdata.sh.OLIG %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(Shannon.avg = fit) + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) + +# Compute 3-year rollmean of AFDW_mean +Shannon.OLIG$rollmean <- zoo::rollmean(Shannon.OLIG$Shannon.avg, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.sh.OLIG <- Shannon.OLIG$rollmean[Shannon.OLIG$Jaar == 2016] + +Shan.avg.OLIG.x <- Shan.avg.OLIG + + geom_hline(yintercept = rollmean_2016.sh.OLIG, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_ribbon.sh.OLIG, aes(x = Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +Shan.avg.OLIG.x + + +### Figuur voor Zoet #### + +Shan.avg.ZOET <- ggplot(data = Shannon.ZOET, aes(x = Jaar, y = Shannon.avg)) + + geom_line(size = 1.2) + + geom_line(aes(y=zoo::rollmean(Shannon.avg, 3, na.pad=TRUE)), size = 1.2, colour="red") + + geom_line(size = 1.2) + + labs(x = "", y="") + + #theme(axis.text.x = element_text(angle = 45))+ + ggtitle("Zoet") + + theme(title = element_text(size=10)) + + scale_y_continuous(breaks=seq(0.4,1.8,0.2), limits = c(0.4, 1.8)) +Shan.avg.ZOET +## For illustration add a ribbon of variation based on 1.96se of simple smoothed gam + +gam_model <- mgcv::gam(Shannon.avg ~ s(Jaar, k = 5), data = Shannon.ZOET) + +# Predict & add predicted values & CI +newdata.sh.ZOET <- data.frame(Jaar = Shannon.ZOET$Jaar) +pred.sh.ZOET <- predict(gam_model, newdata.sh.ZOET = newdata.sh.ZOET, se.fit = TRUE) +newdata.sh.ZOET$fit <- pred.sh.ZOET$fit +newdata.sh.ZOET$se <- pred.sh.ZOET$se.fit +newdata.sh.ZOET$upper <- newdata.sh.ZOET$fit + 1.96 * newdata.sh.ZOET$se +newdata.sh.ZOET$lower <- newdata.sh.ZOET$fit - 1.96 * newdata.sh.ZOET$se + +# Filter to only positive fitted values for plotting purposes +newdata_ribbon.sh.ZOET <- newdata.sh.ZOET %>% dplyr::mutate(lower = ifelse(lower<0, 0, lower)) %>% rename(Shannon.avg = fit) + +### Add threshold, strong bias by exceptional 2014, and single year is always prone to variation so extract rollmean of 2016 (first year out of influence of 2014 with rollmean=3) + +# Compute 3-year rollmean of AFDW_mean +Shannon.ZOET$rollmean <- zoo::rollmean(Shannon.ZOET$Shannon.avg, 3, na.pad = TRUE) + +# Extract rollmean value for year 2016 +rollmean_2016.sh.ZOET <- Shannon.ZOET$rollmean[Shannon.ZOET$Jaar == 2016] + +Shan.avg.ZOET.x <- Shan.avg.ZOET + + geom_hline(yintercept = rollmean_2016.sh.ZOET, linetype = "dashed", colour = "red") + + geom_ribbon(data = newdata_ribbon.sh.ZOET, aes(x = Jaar, ymin = lower, ymax = upper), alpha = 0.2) + + theme_bw(base_family = "sans") + + theme(panel.border = element_blank()) + + theme(plot.title = element_text(hjust = 0.5)) + + theme(plot.margin = unit(c(0,0,0,0), 'lines')) +Shan.avg.ZOET.x + +####Gezamenlijk figuur + +t.SHANNON <- ggarrange(Shan.avg.ZS.x + font("xy.text", size = fnt), + Shan.avg.SS.x + font("xy.text", size = fnt), + Shan.avg.OLIG.x + font("xy.text", size = fnt), + Shan.avg.ZOET.x + font("xy.text", size = fnt), ncol=2, nrow = 2) +annotate_figure(t.SHANNON, left = text_grob("Gemiddelde Shannon-diversity per maand per zone", rot=90), bottom=text_grob("Jaar")) + +ggsave(paste0(pad_figuren, "080-figuur-Shannon.per.maand_ZONES2.jpg"), height=4, width=6) +``` + + + +```{r 080-meta-data} + +meta_data <- + enframe(c(laatstejaar = laatste_jaar, + vroegstejaar = vroegste_jaar, + aantal_stalen = n_staal), + name = "naam", value = "waarde") + +meta_data %>% + write_delim(paste0(pad_data, "meta_data.csv"), + delim = ";") + +``` + diff --git a/moneos_2025/080_hyperbenthos/080_hyperbenthos_data.Rmd b/moneos_2025/080_hyperbenthos/080_hyperbenthos_data.Rmd new file mode 100644 index 0000000..3fa494d --- /dev/null +++ b/moneos_2025/080_hyperbenthos/080_hyperbenthos_data.Rmd @@ -0,0 +1,255 @@ +--- +params: + hoofdstuk: "080_hyperbenthos" +knit: (function(inputFile, ...) { + rmarkdown::render(inputFile, + output_dir = paste0(rmarkdown::yaml_front_matter(inputFile)$params$hoofdstuk, "/output"))}) +title: "Hyperbenthos data" +output: word_document +editor_options: + chunk_output_type: console +--- + + +```{r setup, include=FALSE} + +knitr::opts_chunk$set(echo = FALSE, error=FALSE, warning=FALSE, message=FALSE, cache=FALSE) + +``` + + +```{r libraries} + +library(rprojroot) +library(tidyverse) +library(lubridate) +library(readxl) +library(writexl) +library(RODBC) +library(ggpubr) + +``` + + +```{r pad} + +# inlezen van variabelen +# pad naar data : pad_data +# pad naar tabellen : pad_tabellen +# pad naar figuren : pad_figuren + +source(find_root_file("pad.R", criterion = is_rstudio_project)) #../ weg gedaan want gaf foutpad na Gedeelde Drives? + +pad_data <- maak_pad(params$hoofdstuk, "data") +pad_figuren <- maak_pad(params$hoofdstuk, "figuren") +pad_tabellen <- maak_pad(params$hoofdstuk, "tabellen") + +``` + + +```{r data} +#eerst recente kopie van op Amazon op G gezet +MDB <- odbcConnectAccess2007("G:/Gedeelde drives/PRJ_SCHELDE/Benthos/HyperEpibenthos/DATA/HYPERBENTHOS_SCHELDEjuli2025.accdb") + +# querry met mulitply (substaal nr) al gebruikt voor aantal, WW en AFDW, daardoor is AFDW niet meer DW-AW in uiteindelijke tabel! +sqlCode <- " +SELECT tblSample.campagne, tblStaalnameGebied.Naam AS gebied, tblSample.event_name, tblSample.datum, tblHyperBiot.fractie, tblTaxaList.hoger_taxon, tblHyperBiot.soort, tblTaxaList.vis_NL, tblTaxaList.exoot, Sum([aantal]*[multiply]) AS n, Sum([WW]*[multiply]) AS WWs, tblHyperBiot.DW, tblHyperBiot.AW, Sum([tblHyperBiot].[multiply]*([tblHyperBiot].[DW]-[tblHyperBiot].[AW])) AS AFDW, tblSample.materiaal, tblHyperBiot.Invoerder, tblHyperBiot.Aanmaakdatum, tblHyperBiot.[Laatst gewijzigd] +FROM tblStaalnameGebied INNER JOIN (tblCampagne INNER JOIN (tblTaxaList INNER JOIN (tblSample INNER JOIN tblHyperBiot ON (tblSample.Id = tblHyperBiot.IdSample) AND (tblSample.repeat = tblHyperBiot.repeat)) ON tblTaxaList.soortnaam = tblHyperBiot.soort) ON tblCampagne.campagnecode = tblSample.campagne) ON tblStaalnameGebied.Id = tblSample.IdStaalnameGebied +GROUP BY tblSample.campagne, tblStaalnameGebied.Naam, tblSample.event_name, tblSample.datum, tblHyperBiot.fractie, tblTaxaList.hoger_taxon, tblHyperBiot.soort, tblTaxaList.vis_NL, tblTaxaList.exoot, tblHyperBiot.WW, tblHyperBiot.DW, tblHyperBiot.AW, tblSample.materiaal, tblHyperBiot.Invoerder, tblHyperBiot.Aanmaakdatum, tblHyperBiot.[Laatst gewijzigd], tblTaxaList.taxonredux, tblSample.repeat, tblHyperBiot.contaminatie +HAVING (((tblTaxaList.taxonredux)='in') AND ((tblSample.repeat)=1) AND ((tblHyperBiot.contaminatie) Is Null)); +" + + +#, tblTaxaList.exoot tblTaxaList.exoot + +tel.0 <- sqlQuery(channel = MDB, sqlCode, stringsAsFactors=FALSE) +moniloc <- c("Paardenschor","St-Anna","Ballooi","Dendermonde","Brede Schoren","Rupel") + +#data van BoVO campagnes, van 6 stations, dubbelcampagne van 2014 wegdoen +tel <- tel.0 %>% + dplyr::mutate(Jaar = year(datum), + Maand = month(datum)) %>% + dplyr::mutate(gebied = if_else(gebied == "Konkelschoor - KS" & Jaar == 2013, "Brede Schoren", gebied), + gebied = if_else(gebied == "Vlassenbroek" & Jaar == 2013, "Dendermonde", gebied)) %>% + dplyr::filter(gebied %in% moniloc, str_detect(campagne, "BoVo")) %>% + dplyr::filter(!str_detect(campagne,"_1")) %>% + dplyr::filter(Jaar!=2025) %>% + rename(WW = WWs) + + +jaarorder <- c("2013", "2014", "2015", "2016","2017","2018","2019", "2020", "2021", "2022", "2023", "2024") + + +overzicht <- tel %>% + distinct(datum, gebied, Maand, Jaar) %>% + count(gebied, Maand, Jaar) %>% + dplyr::mutate(Jaar = ordered(Jaar, levels=jaarorder)) %>% + pivot_wider(names_from = Maand, values_from= n) %>% + knitr::kable() + +overzicht + +``` + + + +```{r} +# Correcties via groepsregressies WW-AFDW en n-AFDW voor niet-outliers, etc +tel.1 <- tel %>% + dplyr::group_by(soort) %>% + dplyr::mutate( + # Calculate WW/AFDW ratio + ratio_WW_AFDW = WW / AFDW, + + # Calculate species-level quantiles and IQR for WW/AFDW ratio + Q1_ratio = as.numeric(quantile(ratio_WW_AFDW, 0.25, na.rm = TRUE)), + Q3_ratio = as.numeric(quantile(ratio_WW_AFDW, 0.75, na.rm = TRUE)), + IQR_ratio = Q3_ratio - Q1_ratio, + lower_bound_ratio = Q1_ratio - 1.5 * IQR_ratio, + upper_bound_ratio = Q3_ratio + 1.5 * IQR_ratio, + + # Define outliers for WW/AFDW ratio + outlier_AFDW_WW = ifelse( + ratio_WW_AFDW < lower_bound_ratio | ratio_WW_AFDW > upper_bound_ratio, + "Outlier", + "Non-Outlier" + ), + + # Identify inliers for slope calculation + inlier_WW_AFDW = !is.na(ratio_WW_AFDW) & + ratio_WW_AFDW >= lower_bound_ratio & + ratio_WW_AFDW <= upper_bound_ratio + ) %>% + mutate( + # Fit linear model for inliers only and extract slope + slopeAFDW_WW = if (any(inlier_WW_AFDW, na.rm = TRUE)) { + as.numeric(coef(lm(AFDW ~ WW, data = pick(everything())[inlier_WW_AFDW, , drop = FALSE]))[2]) + } else { + NA_real_ + } + ) %>% + ungroup() %>% + group_by(soort) %>% + mutate( + # Calculate AFDW/n + ratio_AFDW_n = AFDW / n, + + # Calculate species-level quantiles and IQR for AFDW/n + Q1_ratio_n = as.numeric(quantile(ratio_AFDW_n, 0.25, na.rm = TRUE)), + Q3_ratio_n = as.numeric(quantile(ratio_AFDW_n, 0.75, na.rm = TRUE)), + IQR_ratio_n = Q3_ratio_n - Q1_ratio_n, + lower_bound_ratio_n = Q1_ratio_n - 1.5 * IQR_ratio_n, + upper_bound_ratio_n = Q3_ratio_n + 1.5 * IQR_ratio_n, + + # Define outliers for AFDW/n ratio + outlier_AFDW_n = ifelse( + ratio_AFDW_n < lower_bound_ratio_n | ratio_AFDW_n > upper_bound_ratio_n, + 1,0), + inlier_AFDW_n = !is.na(ratio_AFDW_n) & + ratio_AFDW_n >= lower_bound_ratio_n & + ratio_AFDW_n <= upper_bound_ratio_n) %>% + + mutate( + slopeAFDW_n = if (any(inlier_AFDW_n, na.rm = TRUE)) { + as.numeric(coef(lm(AFDW ~ n, data = cur_data()[inlier_AFDW_n, , drop = FALSE]))[2]) + } else { + NA_real_ + } + ) %>% + ungroup() %>% + dplyr::mutate(AFDWcor.1 = dplyr::if_else(outlier_AFDW_WW == "Outlier" & outlier_AFDW_n == 1 & !is.na(WW) & WW !=0 & slopeAFDW_WW >0.1, WW*slopeAFDW_WW, dplyr::if_else(AFDW == 0 | AFDW<0 | is.na(AFDW) & WW !=0, WW*slopeAFDW_WW, AFDW))) %>% + dplyr::mutate(AFDWcor.2 = dplyr::if_else(is.na(AFDWcor.1) & WW !=0 & slopeAFDW_WW >0.1, WW*slopeAFDW_WW, AFDWcor.1)) %>% + dplyr::mutate(AFDWcor.3 = dplyr::if_else(is.na(AFDWcor.2) & WW !=0 & is.na(slopeAFDW_WW), WW*0.17, AFDWcor.2)) %>% + dplyr::mutate(AFDWcor.4 = dplyr::if_else(is.na(AFDWcor.3) & slopeAFDW_n>0 & n>0, n*slopeAFDW_n, AFDWcor.3)) %>% +dplyr::mutate(AFDWcor.5 = dplyr::if_else(outlier_AFDW_WW == "Outlier" & soort == "Dicentrarchus labrax" & !is.na(WW) & !is.na(AFDW), WW*slopeAFDW_WW, AFDWcor.4)) %>% + dplyr::mutate(AFDWcor.7 = dplyr::if_else(WW/AFDWcor.5 < 2 & slopeAFDW_WW > 0.1 & !is.na(slopeAFDW_WW) & !is.na(WW), WW*slopeAFDW_WW, AFDWcor.5)) %>% + dplyr::mutate(AFDWcor.8 = dplyr::if_else(WW/AFDWcor.7 > 15 & slopeAFDW_WW > 0.1 & !is.na(WW) & !is.na(slopeAFDW_WW) & outlier_AFDW_n == 0, WW*slopeAFDW_WW, AFDWcor.7)) %>% + dplyr::mutate(AFDWcor.9 = dplyr::if_else(WW < AFDWcor.8 & !is.na(WW) & !is.na(slopeAFDW_WW) & slopeAFDW_WW > 0.1, WW*slopeAFDW_WW, AFDWcor.8)) + #relocate(c(AFDWcor.4, AFDWcor.5, AFDWcor.6, AFDWcor.7, AFDWcor.8, AFDWcor.9, soort, outlier_AFDW_WW, outlier_AFDW_n), .after = AFDW) + +telc.1 <- tel.1 %>% + dplyr::select(campagne, gebied, event_name, datum, fractie, hoger_taxon, soort, vis_NL, exoot, n, WW, DW, AW, AFDW = AFDWcor.9, materiaal, Invoerder, Aanmaakdatum, `Laatst gewijzigd`, Jaar, Maand) +``` + +```{r check 2: dubbels, vreemde 0-en etc} +#data die moeten w aangepast in hyperdatabase +tel.2 <- tel.1 %>% + dplyr::filter(AFDW != AFDWcor.9 | n == 0 | is.na(n)) %>% + dplyr::select(campagne, gebied, datum, fractie, soort, n, WW, DW, AW, AFDW, AFDWcorrected = AFDWcor.9, materiaal, Invoerder, Aanmaakdatum) + +file_name1 <- + paste0(pad_data, "Hyperbenthos_records_NeedRevision", ".xlsx") + + +write_xlsx(tel.2, + path = file_name1) + +``` + +```{r om data te checken, figuren maken v ratio ww/afdw per soort} +# Loop through each unique species and save a separate plot +soort_lijst <- unique(tel.1$soort) + +for (soort in soort_lijst) { + # Subset data for the current species + soort_data <- tel.1[tel.1$soort == soort, ] + + # Identify outliers using the 1.5 * IQR rule + soort_data$ratio <- soort_data$WW / soort_data$AFDWcor.9 + Q1 <- quantile(soort_data$ratio, 0.25, na.rm = TRUE) + Q3 <- quantile(soort_data$ratio, 0.75, na.rm = TRUE) + IQR <- Q3 - Q1 + lower_bound <- Q1 - 1.5 * IQR + upper_bound <- Q3 + 1.5 * IQR + soort_data$outlier <- ifelse(soort_data$ratio < lower_bound | soort_data$ratio > upper_bound, "Outlier", "Non-Outlier") + + # Create the plot + plot <- ggplot(soort_data, aes(x = WW, y = AFDWcor.9, color=outlier)) + + geom_point() + # Add points + scale_color_manual(values = c("Outlier" = "red", "Non-Outlier" = "blue")) + # Color mapping + geom_smooth(method = "lm", se = FALSE) + # Add linear smoother + stat_regline_equation(aes(label = ..eq.label..), label.x.npc = "left", label.y.npc = "top") + # Add equation + theme_minimal() + # Use a minimal theme + labs( + title = paste("WW vs AFDW for", soort), + x = "Wet Weight (WW)", + y = "Ash-Free Dry Weight (AFDW)" + ) + + # Save the plot to a JPG file + filename <- paste0("G:/Gedeelde drives/PRJ_SCHELDE/MONEOS_rapportage/2025/080_hyperbenthos/data/ratio_figuren_ww_afdw", "/", soort, "_plotWW_AFDW.jpg") # Create a filename + ggsave(filename, plot = plot, width = 6, height = 4, dpi = 300) # Save the plot +} + +``` + + + +data gecreëerd op `r Sys.time()` + +data weggeschreven naar `r paste0(pad_data, "template_data.csv")` + +```{r jaren} + +jaren <- + telc.1 %>% + distinct(Jaar) %>% + pull(Jaar) + +jaar_range <- + range(jaren) + +``` + +```{r wegschrijven-data} +file_name <- + paste0(pad_data, "hyperbenthos_data_revised", paste(jaar_range, collapse = "_"), ".xlsx") + + +write_xlsx(telc.1, + path = file_name) +``` + + + + diff --git a/moneos_2025/150_geintegreerd_rapport/080_hyperbenthos.Rmd b/moneos_2025/150_geintegreerd_rapport/080_hyperbenthos.Rmd new file mode 100644 index 0000000..fa757cf --- /dev/null +++ b/moneos_2025/150_geintegreerd_rapport/080_hyperbenthos.Rmd @@ -0,0 +1,354 @@ +--- +editor_options: + markdown: + wrap: sentence +--- + +```{r 080-hoofdstuk, include=FALSE} + +hoofdstuk <- "080_hyperbenthos" + +``` + +```{r 080-setup, include=FALSE} + +knitr::opts_chunk$set(echo = FALSE, error=FALSE, warning=FALSE, message=FALSE, cache=FALSE, fig.pos = "H") +knitr::opts_knit$set(eval.after = "fig.cap") + +``` + +```{r 080-libraries} + +library(tidyverse) +library(readxl) +library(kableExtra) +library(INBOtheme) +library(rprojroot) ## workaround pad + +``` + +```{r 080-pad} + +# pad naar data : pad_data +# pad naar tabellen : pad_tabellen +# pad naar figuren : pad_figuren + +source(find_root_file("../pad.R", criterion = is_rstudio_project)) + +pad_data <- maak_pad(hoofdstuk, "data") +pad_figuren <- maak_pad(hoofdstuk, "figuren") +pad_tabellen <- maak_pad(hoofdstuk, "tabellen") +``` + +```{r 080-meta_data, eval=FALSE, include=FALSE} +##metadata nog niet aangepast in 2022, wat moet dat zijn? +#meta_data <- + # read_delim(paste0(pad_data, "meta_data.csv"), + # delim = ";") + +#for(i in 1:nrow(meta_data)){ + ##first extract the object value + # tempobj=meta_data$waarde[i] + ##now create a new variable with the original name of the list item + # eval(parse(text=paste(meta_data$naam[i],"= tempobj"))) +#} +``` + +# Hyperbenthos + +Fichenummer: S-DS-V-003 - Hyperbenthos + +**Frank Van de Meutter**, Dimitri Buerms, Ada Coudenys, Charles Lefranc, Bram Loos, Anouk Organe, Vincent Smeekens, Jan Soors + +## Inleiding + +Onder hyperbenthos verstaan we alle kleine fauna (1 mm tot enkele cm) die op en net boven de bodem leeft. +In de Zeeschelde betreft het vooral garnalen en krabben (Decapoda), aasgarnalen (Mysida) met daarnaast ook een groot aandeel juveniele vis. +De monitoring van het hyperbenthos in de Zeeschelde op zes vaste locaties startte in 2013. +Vóór 2013 periode gebeurden op (sommige) van deze zes stations al vangsten met een andere frequentie (zie verder) maar dezelfde methode. +Bij de rapportage gebruiken we doorgaans 2014 als aanvangsjaar, omdat toen voor het eerst een volledig seizoen bemonsterd werd, wat de vergelijkingen en trendbepalingen vergemakkelijkt. + +Een belangrijk verschil met eerdere rapportages is dat we — met enige vertraging — de rapportage in lijn brengen met de nieuwste EMSE evaluatiecriteria (Consortium Schelde in Beeld, 2022). +Daarbij wordt bijvoorbeeld de definitie van hyperbenthos voor abundantie en biomassa-evaluatie verengd tot aasgarnalen, garnalen en steurgarnalen. +We toetsen en illustreren nu ook expliciet de meest recente evaluatiecriteria, die gelden per saliniteitszone. +Hierna vermelden we per toetsparameter steeds de specifieke evaluatiecriteria, en eventuele wijzigingen van de rapportage tegenover eerdere rapportages. + +De gegevens van 2013 tot en met 2024 voor alle hyperbenthos soorten worden geleverd in een Excel-bestand (S_DS_V_003_hyperbenthos_data2013-2024_rapportage2025.xlsx). + +## Materiaal en methode + +### Strategie + +Vijf vaste locaties langsheen de Zeeschelde en één langs de Rupel worden vanaf 2014 maandelijks bemonsterd van april tot oktober. +Volgens de nieuwste criteria (Consortium Schelde in Beeld, 2022) worden de evaluatiecriteria voor de Zeeschelde opgesteld voor vier deelzones. +Van onder- naar bovenstrooms zijn dit: de zone Sterke Saliniteitsgradiënt, de zone Oligohalien, de zone Zoet met lange verblijftijd en de zone Zoet met korte verblijftijd. +Onze vaste monitoringlocaties Paardenschor en Sint-Anna liggen in de zone Sterke Salineitsgradiënt en de stations Ballooi en ook Rupel rekenen we hier tot het Oligohalien. +Voor beide Zoete zones van de Zeeschelde hebben we maar 1 station, waarbij Brede Schoren (bij Berlare) in Zoet met korte verblijftijd ligt, en station Dendermonde in de zone Zoet met lange verblijftijd, maar wel tegen de grens met de zone Zoet korte verblijftijd aan (zie kaart Figuur \@ref(fig:080-figuur10-kaart)). +Vanuit pragmatisch oogpunt, en omdat de ecologische verschillen er beperkt zijn, evalueren we daarom beide zoete zones samen op basis van data afkomstig van deze twee stations. +In de rapportage behandelen we hierna dus steeds 3 zones: de zone Sterke Saliniteitsgradiënt, de zone Oligohalien en de zone Zoet. + +```{r 080-figuur10-kaart, fig.cap=caption_fig1kaart, fig.height=3, fig.width=4.5, out.width="80%"} +caption_fig1kaart <- "Situering staalnamelocaties hyperbenthos. Sampling stations worden aangeduid door een driehoek, het cijfer in de driehoek is de afstand tot de monding (km). Naamgeving: M1=Paardenschor, M2=St. Anna, O1=Ballooi, O2=Rupel, F1=Dendermonde, F2=Brede Schoren." + + +knitr::include_graphics(paste0(pad_figuren, "080-Kaart_hypersamplings.png")) +``` + +
+ +### Staalname + +De bemonstering gebeurt telkens rond het laagwatertijdstip in de dagen rond springtij. +Twee personen slepen een net met cirkelvormige opening (diameter: 50 cm) over een vast traject van 2 x 100 m (heen en terug). +Het net heeft een maaswijdte van 1 mm. +Een stroomsnelheidsmeter wordt in het net opgehangen om het watervolume dat door het net gaat (en dat bemonsterd werd) te kwantificeren. +Na de sleep wordt de vangst gefixeerd met F-Solv (glutaaraldehyde). +Bijkomende metingen van omgevingsvariabelen worden verricht met een multimeter ter bepaling van de saliniteit, het zuurstofgehalte en de watertemperatuur en de gemeten waarden worden genoteerd. +Per bemonstering wordt een waterstaal verzameld om het gehalte aan zwevende stof en de organische fractie ervan achteraf te bepalen. +Dit staal wordt bij laag water rond de waterkering genomen waarbij de persoon op heupdiepte in het water staat en water verzamelt op ca. +20 cm onder het wateroppervlak. + +### Verwerking + +De stalen worden in het labo gespoeld over een 1mm-zeef en alle organismen worden uitgeselecteerd, tot op soort gedetermineerd (tenzij dat niet mogelijk is, in dat geval tot op maximale taxonomische resolutie) en per soort geteld. +Als finale variabele voor analyse werden vroeger de getelde aantallen gestandaardiseerd naar aantal per m³, door de vansgtaantallen te delen door het gemeten watervolume dat door het net is gegaan, indien gegevens over dit volume beschikbaar zijn. +Deze correctie is echter niet aangewezen voor organismen die op de bodem leven (epibenthische soorten, bijvoorbeeld veel garnalen), omdat hun aantallen en biomassa in relatie tot de lengte van het transect staan, en niet in relatie tot het bemonsterd watervolume. +De vangstmethode zelf is bovendien zo opgesteld dat het watervolume bij elk vangbeurt zeer vergelijkbaar is: er wordt gevangen bij de tijkering met minimale stroming, en er wordt een gelijke lengte stroomop- en stroomaf gewandeld met het bongonet (zodat eventuele verschillen als gevolg van stroming elkaar opheffen). +De stroomsnelheidsmeters geven bovendien een minder accuraat beeld wanneer het net zeer traag getrokken wordt of bij frequente stops (bij moeilijk bewandelbare bodems) en wanneer het net stroomafwaarts getrokken wordt (bij lage effectieve stroming door het net stopt de propeller soms). +In deze gevallen werden onderschattingen tot 30% van het bemonsterd watervolume opgemerkt (INBO, niet gepubliceerde gegevens). +Al deze argumenten samen leidden ons tot de conclusie dat het met de gebruikte vangstmethode en de grote vertegenwoordiging van epibenthische taxa wellicht correcter is om uit te gaan van een vast vangvolume van 40m³. +In deze en volgende rapportages gebruiken we daarom de niet-gecorrigeerde vangstaantallen en biomassa (per 40 m³). + +Om de biomassa te bepalen worden de dieren vervolgens per soort verzameld in een kroes, gedroogd, gewogen (ter bepaling van droog gewicht (DW in g)), verast en opnieuw gewogen (ter bepaling van het asgewicht (AW in g)) waarna de uiteindelijke biomassa als asvrij drooggewicht (AFDW in g) berekend wordt door DW-AW (zie ook procedure biomassabepaling macrobenthos). + +## Resultaten: data-analyse hyperbenthos + +Zoals eerder aangehaald volgt deze rapportage de nieuwste evaluatiestrategie, die opgesteld is voor de deelzones, en die het hyperbenthos verengt tot aasgarnalen, garnalen en steurgarnalen (verder: de EMSE-soorten). +Uit de Figuur \@ref(fig:080-figuur-densiteit-EMSEvsREST) blijkt dat voor wat betreft de aantallen het aandeel EMSE-soorten dominant is in de zone Sterke Saliniteitsgradiënt, maar geleidelijk afneemt naar de zone Oligohalien en de Zoete Zeeschelde. +In die laatste twee zones maken de EMSE-soorten vaak minder dan de helft uit van het hyperbenthos. +Voor biomassa (zie \@ref(fig:080-figuur-biomassa-EMSEvsREST)) is het beeld veel genuanceerder. +Het biomassa-aandeel van de EMSE soorten is minder sterk verschillend tussen de zones, en de REST groep van soorten heeft overal een belangrijk aandeel in de totale biomassa van het hyperbenthos (in ruime zin). + +```{r 080-figuur-densiteit-EMSEvsREST, fig.cap=caption_figdensEMSE, out.width="100%"} +caption_figdensEMSE <- "Jaarsom per zone (overheen 2 stations en 7 maanden) van hyperbenthos aantallen voor de periode 2014-2024, voor EMSE-soorten (EMSE) en de overige soorten (REST)." + +knitr::include_graphics(paste0(pad_figuren, "080-figuur-densiteit_EMSEvsREST.jpg")) +``` + +```{r 080-figuur-biomassa-EMSEvsREST, fig.cap=caption_figbiomEMSE, out.width="100%"} +caption_figbiomEMSE <- "Jaarsom per zone (overheen 2 stations en 7 maanden) van hyperbenthos biomassa (g AFDW) voor de periode 2014—2024 voor EMSE-soorten (EMSE) en de overige soorten (REST)." + +knitr::include_graphics(paste0(pad_figuren, "080-figuur-biomassa_EMSEvsREST.jpg")) +``` + +### Toetsparameter: Abundantie + +De toetsparameter abundantie vereist volgens de evaluatiestrategie een abundantie van EMSE-soorten van hyperbenthos (zie eerder) per m², waarop een gemiddelde (EDIT: vermoedelijk per vangst) berekend wordt. +Dit gemiddelde mag niet dalen sinds de T2015, noch toenemen met meer dan 25%. +Voor de Zeeschelde wordt hyperbenthos niet per oppervlakte (m²) bepaald (zie Materiaal en methode), maar per samplevolume (40 m³). +Omdat steeds een vaste transectlengte bemonsterd wordt, is onze methode ook gestandaardiseerd naar oppervlakte, zodat ze ook aan de vereiste van een gestandaardiseerde oppervlakte voldoet. + +Omdat aantallen van hyperbenthos in de Zeeschelde sterk kunnen verschillen tussen jaren, met name door weersextremen (De Neve et al. 2020), zijn een ijkpunt (referentiejaar) vastgesteld op 1 jaar en evaluaties op basis van de vangst van slechts 1 jaar – zoals voorgesteld binnen Consortium Schelde in Beeld (2022) – minder aangewezen. +Naast de effectieve jaarlijkse abundanties (per zone) geven we daarom ook een rollend gemiddelde overheen 3 jaren, en nemen we als ijkpunt het rollend gemiddelde van 2016. +De keuze voor 2016 (en niet voor 2015, wat het oudste beschikbare rollend gemiddelde is), komt doordat 2014 een bijzonder afwijkend jaar was, met enorme aantallen hyperbenthos in grote delen van de Zeeschelde. +Dat jaar staat bijvoorbeeld ook gekend als een anomalie voor wat vissen betreft, door een enorme densiteit aan spiering in dat jaar (zie hoofdstuk Vissen). +Hoewel de keuze voor een referentietoestand altijd vatbaar is voor discussie + +De voor de Zeeschelde aangepaste evaluatie van de abundantie van hyperbenthos gebeurt dus door per zone het actuele rollend gemiddelde te vergelijken met het rollend gemiddelde van 2016 (weergegeven in de figuur als een rode horizontale streepjeslijn). +Een illustratie van de variatie op de gegevens wordt weergegeven als 1.96\*se van een GAM model (in R mgcv:gam(gem_abundantie \~ s(Jaar, k = 5)). +De trend in de abundantie van hyperbenthos, het rollend gemiddelde overheen 3 jaar en de evaluatiegrens worden weergegeven in Figuur \@ref(fig:080-figuur-dens). + +We zien dat voor alle zones het meest recente rollend gemiddelde overheen 3 jaar onder het abundantiecriterium (rollend gemiddelde 2016) ligt. +Uitgezonderd voor de zone Oligohalien is ook de gemiddelde abundantie in 2024 lager dan de referentie. +Meer diepgaande analyse is nodig naar de oorzaken, maar de afgelopen jaren dienden zich aan als jaren met extremen tijdens en voorafgaand aan het bemonsteringsseizoen, met zowel langdurige droogtes als zeer natte jaren/piekperiodes (zoals in 2024, of de waterbom in 2021). +Zeker voor natte jaren (hoge bovenafvoer, De Neve. et al. 2020) is geweten dat dit de aanwezigheid van hyperbenthos in de Zeeschelde negatief beïnvloedt. + +De EMSE evaluatiemethodiek schrijft voor dat wanneer er zich belangrijke ontwikkelingen voordoen in het hyperbenthos, er dan meer in detail kan gekeken worden naar afzonderlijke soorten en groepen, om meer inzicht te krijgen in deze veranderingen. +In Figuur \@ref(fig:080-figuur-abundantie-soorten) is het langjarig abundantieverloop voor de vier belangrijkste EMSE hyperbenthische soorten weergegeven: de grijze garnaal (*Crangon crangon*), de langsneussteurgarnaal (*Palaemon longirostris*), de steeloog-aasgarnaal (*Mesopodopsis slabberi*) en de brakwateraasgarnaal (*Neomysis integer*). +De figuur maakt ook onderscheid per zone. +Uit de figuur blijkt duidelijk dat de grijze garnaal en de steeloog-aasgarnaal soorten zijn die enkel in de zone Sterke Saliniteitsgradiënt voorkomen, maar dat de langneussteurgarnaal en de brakwateraasgarnaal ook of zelfs bij voorkeur in de Oligohaliene en Zoete zone voorkomen. +De verspreiding van deze soorten langsheen de Zeeschelde verschilt tussen jaren, waarschijnlijk onder invloed van het weer. +Zo was de brakwateraasgarnaal vrij talrijk in 2024, maar werd ze nauwelijks in de zoete Zeeschelde gezien, terwijl ze in 2017, 2020 en 2022 hier het talrijkst was. +De meest opmerkelijke evolutie zien we echter bij de grijze garnaal en de langneussteurgarnaal. +De eerste soort kende in 2024 een historisch dieptepunt met minder dan 300 exs. +in de Zeeschelde op een heel jaar (doorgaans \>2000 exs.). +Van de langneussteurgarnaal zijn in 2023 en 2024 respectievelijk 1 en 13 exemplaren gezien. +Het is onduidelijk waar dit aan ligt. +Naast bovenafvoer (debiet) wat in de Zeeschelde voor een "flush"-effect zorgt, waarbij soorten uitspoelen naar meer zeewaarts gelegen zones (De Neve et al. 2020), zijn ook hoge temperatuursextremen een mogelijke oorzaak, doordat dit de overleving van jonge garnalen vermindert (Consortium Schelde in Beeld 2022). + +```{r 080-figuur-dens, fig.cap=caption_figdens, out.width="100%"} +caption_figdens <- "Gemiddelde densiteit van hyperbenthos per 200 m sleep (blauwe lijn) per jaar, per zone en voor de gehele Zeeschelde. Bijkomend worden een rollend gemiddelde overheen 3 jaar (rode lijn) en een evaluatiegrens (=rollend gemiddelde in 2016, als rode streepjeslijn) getoond." + +knitr::include_graphics(paste0(pad_figuren, "080-figuur-densiteit_zones.EMSE2.jpg")) +``` + +
+ +```{r 080-figuur-abundantie-soorten, fig.cap=caption_fig.biom, out.width="100%"} +caption_fig.biom <- "Gemiddelde densiteit (per sleepvangst) voor vier hyperbenthos-soorten per zone voor de verschillende monitoringsjaren." + +knitr::include_graphics(paste0(pad_figuren, "080-EMSE-soorten-abundantie-zone-jaren.jpg")) +``` + +
+ +Hyperbenthos densiteiten kunnen jaarlijks sterk wisselen in de Zeeschelde, vermoedelijk deels natuurlijk en deels door omgevingsvariabelen die (mee) door de mens bepaald worden (bv. zwevende stof gehaltes). +Hyperbenthos is dus inherent een volatiele groep in de Zeeschelde. +Daarbij komt nog dat de monitoring voor de drie deelzones elk slechts gebaseerd is op 2 stations. +Dit draagt verder bij tot de vrij grote variatie tussen opeenvolgende meetjaren en meetmaanden. +Daar tegenover staat dat herhaalde, en systeemwijde negatieve trends, een duidelijke indicatie zijn voor kwaliteitsverlies. + +
+ +### Toetsparameter: Biomassa + +De toetsparameter biomassa vereist volgens de evaluatiestrategie een biomassa van hyperbenthos EMSE-soorten per m², waarop een gemiddelde (EDIT: vermoedelijk per vangst) berekend wordt. +Dit gemiddelde mag niet dalen sinds de T2015, noch toenemen met meer dan 25%. +Voor de Zeeschelde wordt hyperbenthos niet per oppervlakte (m²) bepaald (zie Materiaal en methode), maar per samplevolume (40 m³). +Omdat steeds een vaste transectlengte bemonsterd wordt, is onze methode ook gestandaardiseerd naar oppervlakte, zodat ze ook aan de vereiste van een gestandaardiseerde oppervlakte voldoet. + +Net als bij abundantie passen we de evaluatiecriteria voor de Zeeschelde lichtjes aan ten opzichte van de richtlijnen in Consortium Schelde in Beeld (2022). +Omdat aantallen van hyperbenthos in de Zeeschelde sterk kunnen verschillen tussen jaren, met name door weersextremen (De Neve et al. 2020), zijn een ijkpunt (referentiejaar) vastgesteld op 1 jaar en evaluaties op basis van de vangst van slechst 1 jaar, ook voor de biomassa, minder aangewezen. +Naast de gemiddelde jaarlijkse biomassa (per zone) geven we daarom ook een rollend gemiddelde overheen 3 jaren, en nemen we als ijkpunt het rollend gemiddelde van 2016. +De keuze voor 2016 (en niet voor 2015, wat het oudste beschikbare rollend gemiddelde is), komt doordat 2014 een bijzonder afwijkend jaar was, met enorme aantallen hyperbenthos in grote delen van de Zeeschelde. +Dat jaar staat bijvoorbeeld ook gekend als een anomalie voor wat vissen betreft, door een enorme densiteit aan spiering in dat jaar (zie hoofdstuk Vissen). + +De voor de Zeeschelde aangepaste evaluatie van de biomassa van het hyperbenthos gebeurt dus analoog aan deze voor densiteit (zie daar voor meer uitleg). +De trend in de biomassa van het hyperbenthos, het rollend gemiddelde ervan overheen 3 jaar en de evaluatiegrens, worden weergegeven in Figuur \@ref(fig:080-figuur-biom). + +We zien dat voor alle zones het meest recente rollend gemiddelde overheen 3 jaar én de gemiddelde biomassa voor 2024 onder het biomassacriterium (rollend gemiddelde van 2016) liggen. +De negatieve trend voor biomassa tekent zich nog duidelijker af dan voor densiteiten. +De reden is dat vooral de grotere (en zwaardere) grijze garnalen en langneussteurgarnalen een aantal opeenvolgende slechte jaren kenden, net als de vrij kleine maar soms erg talrijke steeloog-aasgarnaal (zie Figuur \@ref(080-figuur-biom-soorten)). +De evaluatie van de biomassa is dus uitgesproken negatief. +Voor een korte bespreking verwijzen we naar de toetsparameter abundantie. +In ieder geval is meer onderzoek nodig om inzicht te krijgen in de invloed van klimaat en van eventuele lokale antropogene factoren. + +```{r 080-figuur-biom, fig.cap=caption_figbiom, out.width="100%"} +caption_figbiom <- "Gemiddelde biomassa van hyperbenthos per 200 m sleep (blauwe lijn) per jaar, per zone en voor de gehele Zeeschelde. Bijkomend worden een rollend gemiddelde overheen 3 jaar (rode lijn) en een evaluatiegrens (=rollend gemiddelde in 2016, als rode streepjeslijn) getoond." + +knitr::include_graphics(paste0(pad_figuren, "080-figuur-biomassa_zones.EMSE2.jpg")) +``` + +
+ +```{r 080-figuur-biom-soorten, fig.cap=caption_figbiomsoorten, out.width="100%"} +caption_figbiomsoorten <- "Gemiddelde biomassa van hyperbenthos per 200 m sleep (blauwe lijn) per jaar, per zone en voor de gehele Zeeschelde. Bijkomend worden een rollend gemiddelde overheen 3 jaar (rode lijn) en een evaluatiegrens (=rollend gemiddelde in 2016, als rode streepjeslijn) getoond." + +knitr::include_graphics(paste0(pad_figuren, "080-EMSE-soorten-biomassa-zone-jaren.jpg")) +``` + +
+ +### Toetsparameter: Soortenrijkdom + +De toetsparameter soortenrijkdom kijkt naar het totaal aantal soorten (dus niet enkel de EMSE-soorten) exclusief de exoten, per deelzone. +Dit soortenaantal mag niet significant dalen volgens het EMSE criterium (Consortium Schelde in Beeld 2022). +Totale soortenrijdom op jaarbasis, per deelzone en voor de gehele Zeeschelde, staat weergegeven in Figuur \@ref(fig:080-figuur8). +De figuur toont soortenrijkdom met en zonder exoten, waaruit blijkt dat vooral in de zones Sterke Saliniteitsgradiënt en Oligohalien er een belangrijk aandeel exoten aanwezig is (15—20%). + +```{r 080-figuur8, fig.cap=caption_fig8, out.width="100%"} +caption_fig8 <- "Taxa rijkdom per jaar per deelzone, en voor de gehele Zeeschelde, mét en zonder exoten." + +knitr::include_graphics(paste0(pad_figuren, "080_figuur_soortenrijkdomZS_perjaar.jpg")) +``` + +
+ +Soortenrijkdom in de Zeeschelde wordt sterk beïnvloed door het occasioneel opduiken ("zwerven") van soorten die normaal meer zeewaarts voorkomen. +Een toetsparameter kent bij voorkeur een meer stabiel en minder door toeval gedreven gedrag, en we kozen daarom voor de gemiddelde maandelijkse soortenrijkdom per zone. +De manier waarop we deze rapporteren en evalueren is analoog aan de toetsparameters abundantie en biomassa, en we verwijzen naar die onderdelen voor meer uitleg over de methode. +Dit wil zeggen dat we een rollend gemiddelde overheen 3 jaren als evaluatieparameter gebruiken. +Op dit manier zijn de data gebufferd tegen de inherente jaar-tot-jaar variatie van de hyperbenthos gemeenschap in de Zeeschelde, en ligt de focus meer op herhaalde, repetitieve gebeurtenissen (bv. weersanomalieën) of trends. + +De gemiddelde maandelijkse soortenrijkdom of taxa rijkdom (in het geval niet alle exemplaren tot op soort herkend werden), kent een vrij stabiel verloop doorheen de tijd in de Zeeschelde. +In alle zones was 2024 een jaar met gemiddeld minder soorten dan de referentie (het rollend gemiddelde van 2016), en ook het meest recente rollend gemiddelde scoort onder de evaluatiegrens. +De meest duidelijke daling doorheen de tijd noteren we voor de zone Oligohalien. +Een van de redenen voor een afname van de gemiddelde maandelijkse soortenrijkdom is de eerder gemelde crash van garnalen en steurgarnalen, die daardoor in de meeste maandelijkse staalnames afwezig waren. +In de zone Sterke Saliniteitsgradiënt is er geen afname. +Wel valt het haaientandprofiel op, met afwisselend tussen jaren verschillen in de gemiddelde soortenrijkdom van ca. +3 soorten, ofwel bijna 25% van het soortaantal. +Die jaarlijkse verschillen zijn opmerkelijk, aangezien het een gemiddeld verschil betreft, op basis van monsters van april tot en met oktober op 2 stations. + +De evaluatie van de toetsparameter soortenrijkdom is gunstig voor de Zeeschelde als geheel en de zone Sterke Saliniteitsgradiënt, maar evolueert negatief voor de Oligohaliene en (in mindere mate) de Zoete zone van de Zeeschelde. + +```{r 080-figuur9, fig.cap=caption_fig9, out.width="100%"} +caption_fig9 <- "Gemiddelde maandelijkse taxa rijkdom van hyperbenthos voor de onderzoeksjaren, per deelzone en voor de gehele Zeeschelde." + +knitr::include_graphics(paste0(pad_figuren, "080-figuur-Soortenrijkdom.per.maand_ZONES2.jpg")) +``` + +
+ +### Toetsparameter: Shannon-index abundantie + +Om de structuur en samenstelling van de hyperbenthos gemeenschap op te volgen, is door Consortium Schelde in Beeld (2022) de Shannon-diversiteits-index voorgesteld. +Deze index geeft een waarde tussen 0 en het natuurlijk logaritme (ln) van soortenrijkdom, en geeft een hoge waarden indien er veel soorten in het monster zitten én indien deze soorten gelijk verdeeld zijn. +Bij afwijkingen, bijvoorbeeld de dominantie door 1 soort, verlaagt de index. +Deze index mag niet significant dalen volgens het EMSE criterium ten opzichte van de T2015 (Consortium Schelde in Beeld 2022), en moet bepaald worden per zone. +Deze toetsparameter wordt bepaald op basis van abundantie. + +We volgen dezelfde evaluatiestrategie als bij de toetsparameter Soortenrijkdom, op basis van een gemiddelde maandelijkse Shannon-index per jaar, en per zone (zie hoger). +Als toetsparameter nemen we het rollend gemiddelde overheen 3 jaar en als referentie nemen we het rollend gemiddelde van het jaar 2016 (voor argumentatie: zie de andere toetsparameters). + +De Shannon diversiteit is een zone-specifieke eigenschap van het hyperbenthos, zo blijkt uit Figuur \@ref(fig:080-figuur10). +In de Zeeschelde schommelt het rollend gemiddelde rond 1.4, in de zone Sterke Saliniteitsgradiënt rond de 1, en in de overige twee zones rond de 1.2. +De Shannon-index is behoorlijk stabiel overheen de tijd, in alle zones. + +```{r 080-figuur10, fig.cap=caption_fig10, out.width="100%"} +caption_fig10 <- "Gemiddelde maandelijkse Shannon diversiteit, per deelzone en voor de volledige Zeeschelde, voor de verschillende monitoringsjaren. De shannon diversiteit werd zowel berekend op densiteiten als voor biomassa." + +knitr::include_graphics(paste0(pad_figuren, "080-figuur-Shannon.per.maand_ZONES2.jpg")) +``` + +
+ +## Algemene conclusie + +**Toetsparameters Abundantie en Biomassa** + +De evaluatie van densiteiten en biomassa van hyperbenthos wordt vanaf het evaluatiejaar 2024 beperkt tot de EMSE-soorten: de aasgarnalen, garnalen en steurgarnalen. +Deze aanpak verschilt van de voorgaande rapportages waarbij ook juveniele vis, vlokreeftjes, etc... +meegerekend werd. +De evolutie van abundantie en biomassa overheen de monitoringsjaren is onderhevig aan grote schommelingen, die samen hangen met goede en slechte jaren van specifieke soorten, wellicht voor een deel aangedreven door weersextremen. +We ontwikkelden daarom een aangepaste evaluatiemethodiek, op basis van een rollend gemiddelde, waardoor jaar-tot-jaar variatie minder tot uiting komt en de nadruk meer ligt op herhaalde of continue trends. +De evaluatie toont aan dat voor abundantie, maar nog meer uitgesproken voor biomassa, er recent een negatieve trend is van het hyperbenthos (EMSE-soorten) in de Zeeschelde. +Deze trend is algemeen, maar meest uitgesproken in de zone Oligohalien. +Een analyse per soort toont aan dat in 2024 zowel grijze garnalen als langneussteurgarnalen, twee soorten die doorgaans een belangrijk aandeel hebben in de biomassa van het hyperbenthos in de Zeeschelde, vrijwel afwezig waren. +Voor de langneussteurgarnaal was dit al het twee jaar op rij waarbij de soort grotendeels ontbrak. +Meer onderzoek is nodig naar de mogelijke oorzaken, waarbij zeker aandacht moet gaan naar de invloed van weersextremen, ook met het oog op toekomstige effecten van klimaatverstoring op de ecologie van de Zeeschelde (De Neve et al. 2020, Consortium Schelde in Beeld 2022). + +**Toetsparameter soortenrijkdom** + +Het rollend gemiddelde van de soortenrijkdom (of taxa rijkdom) van het hyperbenthos in de deelzones van de Zeeschelde kende een vrij stabiel verloop doorheen de tijd. +Enkel in de zone Oligohalien lijkt het huidig rollend gemiddeld lager dan de referentie rond 2016. +In alle zones was de effectieve gemiddelde maandelijkse soortenrijkdom in 2024 lager dan de referentie, mogelijk mede door het vaak ontbreken van garnalen en steurgarnalen in de monsters. +Voor de zone Oligohalien evolueert de toestand voor deze parameter naar ongunstig, voor de andere zones valt deze voorlopig binnen de evaluatiemarge. + +**Toetsparameter Shannon-index diversiteit** + +De gemiddelde maandelijkse Shannon diversiteit is een vrij stabiele, zone-specifieke eigenschap van de hyperbenthos gemeenschappen in de Zeeschelde. +Het is ook een vrij stabiele parameter, die geen duidelijke veranderingen overheen de monitoringsperiode laat zien. +We evalueren deze daarom als gunstig. + +## Referenties + +Consortium Schelde in Beeld (2022). +Evaluatiemethodiek Schelde-estuarium. +Update 2021. +HKV/universiteit Gent/Bureau Waardenburg/ Antea Group: Nederland, Bergen-op-Zoom. +p. 396. + +De Neve L., Van Ryckegem G., Vanoverbeke J., Van de Meutter F., Van Braeckel A., Van den Bergh E., & Speybroeck, J. +(2020). +Hyperbenthos in the upper reaches of the Scheldt estuary (Belgium): Spatiotemporal patterns and ecological drivers of a recovered community. +Estuarine, Coastal and Shelf Science 245: 106967. +DOI: 10.1016/j.ecss.2020.106967. + +Van Ryckegem G., Vanoverbeke J., Van Braeckel A., Van de Meutter F., Mertens W. Mertens A. +& Breine J. +(2021). +MONEOS-Datarapport INBO: toestand Zeeschelde 2020. +Monitoringsoverzicht en 1ste lijnsrapportage Geomorfologie, diversiteit Habitats en diversiteit Soorten. +Rapporten van het Instituut voor Natuur- en Bosonderzoek 2021 (47). +Instituut voor Natuur- en Bosonderzoek, Brussel. +DOI: doi.org/10.21436/inbor.52484672.