slope <- terrain(dem, opt = 'slope')
aspect <- terrain(dem, opt = 'aspect')
hillshade <- hillShade(slope, aspect, angle = 45, direction = 315)
# plot the hillshade, overlay raster with transparency (alpha)
plot(hillshade, col = grey(0:100/100), legend = FALSE,
ylim=c(aoi[3], aoi[4]), xlim=c(aoi[1], aoi[2]),
)
plot(rough_raster,
col = col_palette,
ylim=c(aoi[3], aoi[4]), xlim=c(aoi[1], aoi[2]), legend=F,
zlim = zlim,
add = TRUE,
alpha = 0.8)  # Adjust alpha for transparency
# Add roughness value text in the top-right corner of each plot
text(x = 0.05, y = 0.17, labels = bquote(italic(k[{"s, bed"}]) == .(signif(rough_value * 100, 2)) * " cm"),
col = "black", pos = 4, cex = 3)
# Add y-axis labels only to the plots in the left column
if (i %% 2 == 1) {
mtext("y (m)", side = 2, line = 4, cex=2.2)
}
# Add x-axis labels only to the plots in the bottom row
if (i > 4) {  # Bottom row plots (5, 6)
mtext("x (m)", side = 1, line = 4 , cex=2.2)
}
if (i == 4) {
image.plot( legend.only = TRUE, col = col_palette, zlim = zlim,
legend.args = list(text = paste('Roughness (m)'), side = 3,  #text = paste('Roughness', '(m)', sep = "\n")
line = 3, adj=0.1, cex = 2 ), #cex = 1.5), #
legend.width = 1.5,
legend.mar = 10, # push legend further right
smallplot = c(0.97, 1, 0.3, 0.7))  # adjust as needed
}
}
setwd("~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/")
files<- list.files(path = "~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/data/DataProcessed/Beds/DSMs/")
#files<- list.files(path = "data/DataProcessed/Beds/DSMs/")
select <- c("data/Beds/DEMs/P2_C_0_00_R1_t120_DSM.tif", "data/Beds/DEMs/P2_C_0_00_R2_t120_DSM.tif",
"data/Beds/DEMs/P2_C_E_20_R1_t120_DSM.tif", "data/Beds/DEMs/P2_C_E_20_R2_t120_DSM.tif",
"data/Beds/DEMs/P2_C_S_20_R1_t120_DSM.tif", "data/Beds/DEMs/P2_C_S_20_R2_t120_DSM.tif") #files[c(8, 9, 12, 13, 18, 19)]
plotnames <- c("Control 0.35 m/s", "Control 0.45 m/s",
"Ensis 20% 0.35 m/s","Ensis 20% 0.45 m/s",
"Spisula 20% 0.35 m/s","Ensis 20% 0.45 m/s")
# setwd('../')
# png(filename = paste0("figs/Fig5.png"), width=1400, height=1000)
par(cex.axis = 2.5, cex.lab = 2, cex.main = 3.2)
par(mfrow = c(3, 2), mar = c(3, 4, 2, 3), oma = c(4, 8, 6, 16))
windowsFonts(A=windowsFont("Arial Unicode MS"))
par(family = "A")
col_palette <- colorRampPalette(c("lightyellow3", "#ffa500", "#8B2D2D"))(50)
zlim <- c(0, 0.03)  # Define zlim (limits for roughness values)
breaks <- seq(zlim[1], zlim[2], length.out = length(col_palette) + 1)
i=1
for (i in 1:length(select)){
file <- select[i]
expname <- sub("_t120_DSM.tif", "", file)
DEM <-raster(select[i])
dem <- crop(DEM, aoi)
rough_raster <- focal(dem, w = w, fun = rough_fun, pad = TRUE, padValue = NA)
rough_value <- cellStats(rough_raster, stat = "mean")   # mean roughness value for entire bed
# print(summary(values(dem)) )
# print(summary(values(rough_raster)))
# print(rough_value)
## add hillshading
projection(dem) <- CRS("+proj=utm +zone=33 +datum=WGS84")
slope <- terrain(dem, opt = 'slope')
aspect <- terrain(dem, opt = 'aspect')
hillshade <- hillShade(slope, aspect, angle = 45, direction = 315)
# plot the hillshade, overlay raster with transparency (alpha)
plot(hillshade, col = grey(0:100/100), legend = FALSE,
ylim=c(aoi[3], aoi[4]), xlim=c(aoi[1], aoi[2]),
)
plot(rough_raster,
col = col_palette,
ylim=c(aoi[3], aoi[4]), xlim=c(aoi[1], aoi[2]), legend=F,
zlim = zlim,
add = TRUE,
alpha = 0.8)  # Adjust alpha for transparency
# Add roughness value text in the top-right corner of each plot
text(x = 0.05, y = 0.17, labels = bquote(italic(k[{"s, bed"}]) == .(signif(rough_value * 100, 2)) * " cm"),
col = "black", pos = 4, cex = 3)
# Add y-axis labels only to the plots in the left column
if (i %% 2 == 1) {
mtext("y (m)", side = 2, line = 4, cex=2.2)
}
# Add x-axis labels only to the plots in the bottom row
if (i > 4) {  # Bottom row plots (5, 6)
mtext("x (m)", side = 1, line = 4 , cex=2.2)
}
if (i == 4) {
image.plot( legend.only = TRUE, col = col_palette, zlim = zlim,
legend.args = list(text = paste('Roughness (m)'), side = 3,  #text = paste('Roughness', '(m)', sep = "\n")
line = 3, adj=0.1, cex = 2 ), #cex = 1.5), #
legend.width = 1.5,
legend.mar = 10, # push legend further right
smallplot = c(0.97, 1, 0.3, 0.7))  # adjust as needed
}
}
files
setwd("~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/")
files<- list.files(path = "~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/data/DataProcessed/Beds/DSMs/")
#files<- list.files(path = "data/DataProcessed/Beds/DSMs/")
select <- c("P2_C_0_00_R1_t120_DSM.tif", "P2_C_0_00_R2_t120_DSM.tif",
"P2_C_E_20_R1_t120_DSM.tif", "P2_C_E_20_R2_t120_DSM.tif",
"P2_C_S_20_R1_t120_DSM.tif", "P2_C_S_20_R2_t120_DSM.tif") #files[c(8, 9, 12, 13, 18, 19)]
plotnames <- c("Control 0.35 m/s", "Control 0.45 m/s",
"Ensis 20% 0.35 m/s","Ensis 20% 0.45 m/s",
"Spisula 20% 0.35 m/s","Ensis 20% 0.45 m/s")
# setwd('../')
# png(filename = paste0("figs/Fig5.png"), width=1400, height=1000)
par(cex.axis = 2.5, cex.lab = 2, cex.main = 3.2)
par(mfrow = c(3, 2), mar = c(3, 4, 2, 3), oma = c(4, 8, 6, 16))
windowsFonts(A=windowsFont("Arial Unicode MS"))
par(family = "A")
col_palette <- colorRampPalette(c("lightyellow3", "#ffa500", "#8B2D2D"))(50)
zlim <- c(0, 0.03)  # Define zlim (limits for roughness values)
breaks <- seq(zlim[1], zlim[2], length.out = length(col_palette) + 1)
i=1
for (i in 1:length(select)){
file <- select[i]
expname <- sub("_t120_DSM.tif", "", file)
DEM <-raster(select[i])
dem <- crop(DEM, aoi)
rough_raster <- focal(dem, w = w, fun = rough_fun, pad = TRUE, padValue = NA)
rough_value <- cellStats(rough_raster, stat = "mean")   # mean roughness value for entire bed
# print(summary(values(dem)) )
# print(summary(values(rough_raster)))
# print(rough_value)
## add hillshading
projection(dem) <- CRS("+proj=utm +zone=33 +datum=WGS84")
slope <- terrain(dem, opt = 'slope')
aspect <- terrain(dem, opt = 'aspect')
hillshade <- hillShade(slope, aspect, angle = 45, direction = 315)
# plot the hillshade, overlay raster with transparency (alpha)
plot(hillshade, col = grey(0:100/100), legend = FALSE,
ylim=c(aoi[3], aoi[4]), xlim=c(aoi[1], aoi[2]),
)
plot(rough_raster,
col = col_palette,
ylim=c(aoi[3], aoi[4]), xlim=c(aoi[1], aoi[2]), legend=F,
zlim = zlim,
add = TRUE,
alpha = 0.8)  # Adjust alpha for transparency
# Add roughness value text in the top-right corner of each plot
text(x = 0.05, y = 0.17, labels = bquote(italic(k[{"s, bed"}]) == .(signif(rough_value * 100, 2)) * " cm"),
col = "black", pos = 4, cex = 3)
# Add y-axis labels only to the plots in the left column
if (i %% 2 == 1) {
mtext("y (m)", side = 2, line = 4, cex=2.2)
}
# Add x-axis labels only to the plots in the bottom row
if (i > 4) {  # Bottom row plots (5, 6)
mtext("x (m)", side = 1, line = 4 , cex=2.2)
}
if (i == 4) {
image.plot( legend.only = TRUE, col = col_palette, zlim = zlim,
legend.args = list(text = paste('Roughness (m)'), side = 3,  #text = paste('Roughness', '(m)', sep = "\n")
line = 3, adj=0.1, cex = 2 ), #cex = 1.5), #
legend.width = 1.5,
legend.mar = 10, # push legend further right
smallplot = c(0.97, 1, 0.3, 0.7))  # adjust as needed
}
}
setwd("~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/data/DataProcessed/Beds/DSMs/")
files<- list.files(path = "~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/data/DataProcessed/Beds/DSMs/")
#files<- list.files(path = "data/DataProcessed/Beds/DSMs/")
select <- c("P2_C_0_00_R1_t120_DSM.tif", "P2_C_0_00_R2_t120_DSM.tif",
"P2_C_E_20_R1_t120_DSM.tif", "P2_C_E_20_R2_t120_DSM.tif",
"P2_C_S_20_R1_t120_DSM.tif", "P2_C_S_20_R2_t120_DSM.tif") #files[c(8, 9, 12, 13, 18, 19)]
plotnames <- c("Control 0.35 m/s", "Control 0.45 m/s",
"Ensis 20% 0.35 m/s","Ensis 20% 0.45 m/s",
"Spisula 20% 0.35 m/s","Ensis 20% 0.45 m/s")
# setwd('../')
# png(filename = paste0("figs/Fig5.png"), width=1400, height=1000)
par(cex.axis = 2.5, cex.lab = 2, cex.main = 3.2)
par(mfrow = c(3, 2), mar = c(3, 4, 2, 3), oma = c(4, 8, 6, 16))
windowsFonts(A=windowsFont("Arial Unicode MS"))
par(family = "A")
col_palette <- colorRampPalette(c("lightyellow3", "#ffa500", "#8B2D2D"))(50)
zlim <- c(0, 0.03)  # Define zlim (limits for roughness values)
breaks <- seq(zlim[1], zlim[2], length.out = length(col_palette) + 1)
i=1
for (i in 1:length(select)){
file <- select[i]
expname <- sub("_t120_DSM.tif", "", file)
DEM <-raster(select[i])
dem <- crop(DEM, aoi)
rough_raster <- focal(dem, w = w, fun = rough_fun, pad = TRUE, padValue = NA)
rough_value <- cellStats(rough_raster, stat = "mean")   # mean roughness value for entire bed
# print(summary(values(dem)) )
# print(summary(values(rough_raster)))
# print(rough_value)
## add hillshading
projection(dem) <- CRS("+proj=utm +zone=33 +datum=WGS84")
slope <- terrain(dem, opt = 'slope')
aspect <- terrain(dem, opt = 'aspect')
hillshade <- hillShade(slope, aspect, angle = 45, direction = 315)
# plot the hillshade, overlay raster with transparency (alpha)
plot(hillshade, col = grey(0:100/100), legend = FALSE,
ylim=c(aoi[3], aoi[4]), xlim=c(aoi[1], aoi[2]),
)
plot(rough_raster,
col = col_palette,
ylim=c(aoi[3], aoi[4]), xlim=c(aoi[1], aoi[2]), legend=F,
zlim = zlim,
add = TRUE,
alpha = 0.8)  # Adjust alpha for transparency
# Add roughness value text in the top-right corner of each plot
text(x = 0.05, y = 0.17, labels = bquote(italic(k[{"s, bed"}]) == .(signif(rough_value * 100, 2)) * " cm"),
col = "black", pos = 4, cex = 3)
# Add y-axis labels only to the plots in the left column
if (i %% 2 == 1) {
mtext("y (m)", side = 2, line = 4, cex=2.2)
}
# Add x-axis labels only to the plots in the bottom row
if (i > 4) {  # Bottom row plots (5, 6)
mtext("x (m)", side = 1, line = 4 , cex=2.2)
}
if (i == 4) {
image.plot( legend.only = TRUE, col = col_palette, zlim = zlim,
legend.args = list(text = paste('Roughness (m)'), side = 3,  #text = paste('Roughness', '(m)', sep = "\n")
line = 3, adj=0.1, cex = 2 ), #cex = 1.5), #
legend.width = 1.5,
legend.mar = 10, # push legend further right
smallplot = c(0.97, 1, 0.3, 0.7))  # adjust as needed
}
}
# Add column headers
mtext(expression(paste(0.35," ", m, " ",s^-1)), side = 3, line = 1, outer = TRUE, at = 0.25,  cex = 3)
mtext(expression(paste(0.45," ", m, " ",s^-1)), side = 3, line = 1, outer = TRUE, at = 0.75,  cex = 3)
# Add row headers
mtext("Control (0%)", side = 2, line = 3, outer = TRUE, at = 0.85, cex = 3)
mtext("Ensis 20%", side = 2, line = 3, outer = TRUE, at = 0.5, cex = 3)
mtext("Spisula 20%", side = 2, line = 3 , outer = TRUE, at = 0.15, cex = 3)
# dev.off()
## working directory ###############################################################################
setwd("data/DataRaw/Field")  # "~/1. PhD/Ch3 Shells and sand/Data"
## working directory ###############################################################################
setwd("~/data/DataRaw/Field")  # "~/1. PhD/Ch3 Shells and sand/Data" "C:\Users\tkooistra\OneDrive - NIOZ\Documenten\1. PhD\Ch4 TraSSh - shells and sand\Paper\Data and code\data\DataRaw\Field\202108_shell content_BD.csv"
## working directory ###############################################################################
setwd("~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/data/DataRaw/Field")  # "~/1. PhD/Ch3 Shells and sand/Data" "C:\Users\tkooistra\OneDrive - NIOZ\Documenten\1. PhD\Ch4 TraSSh - shells and sand\Paper\Data and code\data\DataRaw\Field\202108_shell content_BD.csv"
## load data
shells <- read.csv("202108_shellcontent_BD.csv")
## working directory ###############################################################################
setwd("~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/data/DataRaw/Field/")  # "~/1. PhD/Ch3 Shells and sand/Data" "C:\Users\tkooistra\OneDrive - NIOZ\Documenten\1. PhD\Ch4 TraSSh - shells and sand\Paper\Data and code\data\DataRaw\Field\202108_shell content_BD.csv"
## load data
shells <- read.csv("202108_shellcontent_BD.csv")
## load data
shells <- read.csv("202108_shell content_BD.csv")
fieldform <- read.csv("MOSAIC_2021_TRAILS_02_BD_2-8aug.csv")
# prepare data #####################################################################################
#shells <- shells[order(shells$Poskey), ]
fieldform$Shells[is.na(fieldform$Shells)] <- 0
colnames(fieldform)[10] <- "Poskey"
data <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
data <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
dplyr::select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
library(tidyverse)
#transform coordinates to XY
data$Y = as.numeric(data$Nbg)+(as.numeric(data$Nbmin)/60)
data <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
dplyr::select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
# calculate shell fraction (volume %)
data$Vcore <- (56-data$cm_to_top_box)*20*30 # 0.06 m2
data$shell_fraction <- data$Vshell_total_ml/data$Vcore
#transform coordinates to XY
data$Y = as.numeric(data$Nbg)+(as.numeric(data$Nbmin)/60)
data$X = as.numeric(data$olg)+(as.numeric(data$olmin)/60)
data_sp <- data[!is.na(data$X),]
data_sp <- st_as_sf(data_sp, coords = c("X", "Y"), crs = 4326)# convert to spatial object
?st_as_sf
library(ggplot2)
library(sf)
library(terra)
library(viridis)
library(tidyverse)
data_sp <- st_as_sf(data_sp, coords = c("X", "Y"), crs = 4326)# convert to spatial object
my_dpi = 500
w_full = 3740   # Combinations bitmapped line/halftones (color or grayscale): full page
w_single = 1772  # single column
# page size in inches for pdf
w_full = w_full/my_dpi
w_single = w_single/my_dpi
# Map all points over year per month ####
data_sp$shell_perc <- data_sp$shell_fraction*100
p1 <- ggplot(data_sp %>% filter(!is.na(shell_perc))) + #env_select
geom_sf(data=it, linewidth=0.1)+ geom_sf(data=ws) +
geom_sf(aes(colour=shell_perc), alpha = 0.75, size=2)+ #aes(colour=seasons),
lims(x=c(5.45, 5.95), y=c(53.31, 53.51)) +
scale_colour_viridis(name=expression(paste(italic(P[shell]), " (%)"))) +
#scale_colour_manual(breaks=c('winter', 'spring', 'summer', 'autumn'),values = cols) + #breaks=c("23884","22606", "23885", "22621")
#facet_wrap( ~Year) +
theme(legend.position = "bottom", legend.text = element_text((size=35))) +
theme_bw(base_size = 12)
p1 <- ggplot(data_sp %>% filter(!is.na(shell_perc))) + #env_select
#geom_sf(data=it, linewidth=0.1)+ geom_sf(data=ws) +
geom_sf(aes(colour=shell_perc), alpha = 0.75, size=2)+ #aes(colour=seasons),
lims(x=c(5.45, 5.95), y=c(53.31, 53.51)) +
scale_colour_viridis(name=expression(paste(italic(P[shell]), " (%)"))) +
#scale_colour_manual(breaks=c('winter', 'spring', 'summer', 'autumn'),values = cols) + #breaks=c("23884","22606", "23885", "22621")
#facet_wrap( ~Year) +
theme(legend.position = "bottom", legend.text = element_text((size=35))) +
theme_bw(base_size = 12)
p1
p1
## working directory and data #####################################################################
setwd("~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/data/DataRaw/Field/")
#setwd("/data/DataRaw/Field/")
## load data
shells <- read.csv("202108_shell content_BD.csv")
fieldform <- read.csv("MOSAIC_2021_TRAILS_02_BD_2-8aug.csv")
# prepare data #####################################################################################]
fieldform$Shells[is.na(fieldform$Shells)] <- 0
colnames(fieldform)[10] <- "Poskey"
data <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
dplyr::select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
# mapping#####################################################################################
#transform coordinates to XY
data$Y = as.numeric(data$Nbg)+(as.numeric(data$Nbmin)/60)
data$X = as.numeric(data$olg)+(as.numeric(data$olmin)/60)
data_sp <- data[!is.na(data$X),]
data_sp <- st_as_sf(data_sp, coords = c("X", "Y"), crs = 4326)# convert to spatial object
my_dpi = 500
w_full = 3740
# page size in inches for png
w_full = w_full/my_dpi
# Calculate fraction in %
data_sp$shell_perc <- data_sp$shell_fraction*100
setwd("~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/data/DataRaw/Field/")
#setwd("/data/DataRaw/Field/")
## load data
shells <- read.csv("202108_shell content_BD.csv")
fieldform <- read.csv("MOSAIC_2021_TRAILS_02_BD_2-8aug.csv")
# prepare data #####################################################################################]
fieldform$Shells[is.na(fieldform$Shells)] <- 0
colnames(fieldform)[10] <- "Poskey"
data <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
dplyr::select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
# mapping#####################################################################################
#transform coordinates to XY
data$Y = as.numeric(data$Nbg)+(as.numeric(data$Nbmin)/60)
data$X = as.numeric(data$olg)+(as.numeric(data$olmin)/60)
data_sp <- data[!is.na(data$X),]
data_sp <- st_as_sf(data_sp, coords = c("X", "Y"), crs = 4326)# convert to spatial object
my_dpi = 500
w_full = 3740
# page size in inches for png
w_full = w_full/my_dpi
# Calculate fraction in %
data_sp$shell_perc <- data_sp$shell_fraction*100
# calculate shell fraction (volume %)
data$Vcore <- (56-data$cm_to_top_box)*20*30 # 0.06 m2
data$shell_fraction <- data$Vshell_total_ml/data$Vcore
# mapping#####################################################################################
#transform coordinates to XY
data$Y = as.numeric(data$Nbg)+(as.numeric(data$Nbmin)/60)
data$X = as.numeric(data$olg)+(as.numeric(data$olmin)/60)
data_sp <- data[!is.na(data$X),]
data_sp <- st_as_sf(data_sp, coords = c("X", "Y"), crs = 4326)# convert to spatial object
my_dpi = 500
w_full = 3740
# page size in inches for png
w_full = w_full/my_dpi
# Calculate fraction in %
data_sp$shell_perc <- data_sp$shell_fraction*100
# prepare data #####################################################################################]
fieldform$Shells[is.na(fieldform$Shells)] <- 0
colnames(fieldform)[10] <- "Poskey"
data <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
dplyr::select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
# calculate shell fraction (volume %)
data$Vcore <- (56-data$cm_to_top_box)*20*30 # 0.06 m2
data$shell_perc <- data$Vshell_total_ml/data$Vcore * 100
# # Calculate fraction in %
# data_sp$shell_perc <- data_sp$shell_fraction*100
# mapping#####################################################################################
#transform coordinates to XY
data$Y = as.numeric(data$Nbg)+(as.numeric(data$Nbmin)/60)
data$X = as.numeric(data$olg)+(as.numeric(data$olmin)/60)
data_sp <- data[!is.na(data$X),]
data_sp <- st_as_sf(data_sp, coords = c("X", "Y"), crs = 4326)# convert to spatial object
my_dpi = 500
w_full = 3740
# page size in inches for png
w_full = w_full/my_dpi
p1 <- ggplot(data_sp %>% filter(!is.na(shell_perc))) + #env_select
#geom_sf(data=it, linewidth=0.1)+ geom_sf(data=ws) +
geom_sf(aes(colour=shell_perc), alpha = 0.75, size=2)+ #aes(colour=seasons),
lims(x=c(5.45, 5.95), y=c(53.31, 53.51)) +
scale_colour_viridis(name=expression(paste(italic(P[shell]), " (%)"))) +
#scale_colour_manual(breaks=c('winter', 'spring', 'summer', 'autumn'),values = cols) + #breaks=c("23884","22606", "23885", "22621")
#facet_wrap( ~Year) +
theme(legend.position = "bottom", legend.text = element_text((size=35))) +
theme_bw(base_size = 12)
p1
## working directory and data #####################################################################
setwd("~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/data/DataRaw/Field/")
#setwd("/data/DataRaw/Field/")
## load data
shells <- read.csv("202108_shell content_BD.csv")
fieldform <- read.csv("MOSAIC_2021_TRAILS_02_BD_2-8aug.csv")
# prepare data #####################################################################################]
fieldform$Shells[is.na(fieldform$Shells)] <- 0
colnames(fieldform)[10] <- "Poskey"
data21 <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
dplyr::select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
# prepare data 2022 #####################################################################################
shells <- read.csv("TRAILS2202_2205_Vshells.csv")
fieldform <- read.csv("TRAILS2202_2205_2210_boxcores.csv")
#shells <- shells[order(shells$Poskey), ]
colnames(fieldform)[10] <- "Poskey"
fieldform$Shells[is.na(fieldform$Shells)] <- 0
#data <- cbind(fieldform[which(is.element(fieldform$Worms, shells$Poskey)),10:44], shells)
data22 <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
data <- rbind(data21, data22)
# calculate shell fraction (volume %)
data$Vcore <- (56-data$cm_to_top_box)*20*30 # 0.06 m2
data$shell_perc <- data$Vshell_total_ml/data$Vcore * 100
# mapping#####################################################################################
#transform coordinates to XY
data$Y = as.numeric(data$Nbg)+(as.numeric(data$Nbmin)/60)
data$X = as.numeric(data$olg)+(as.numeric(data$olmin)/60)
data_sp <- data[!is.na(data$X),]
data_sp <- st_as_sf(data_sp, coords = c("X", "Y"), crs = 4326)# convert to spatial object
my_dpi = 500
w_full = 3740
# page size in inches for png
w_full = w_full/my_dpi
p1 <- ggplot(data_sp %>% filter(!is.na(shell_perc))) + #env_select
#geom_sf(data=it, linewidth=0.1)+ geom_sf(data=ws) +
geom_sf(aes(colour=shell_perc), alpha = 0.75, size=2)+ #aes(colour=seasons),
lims(x=c(5.45, 5.95), y=c(53.31, 53.51)) +
scale_colour_viridis(name=expression(paste(italic(P[shell]), " (%)"))) +
#scale_colour_manual(breaks=c('winter', 'spring', 'summer', 'autumn'),values = cols) + #breaks=c("23884","22606", "23885", "22621")
#facet_wrap( ~Year) +
theme(legend.position = "bottom", legend.text = element_text((size=35))) +
theme_bw(base_size = 12)
p1
View(data21)
data <- data21
# calculate shell fraction (volume %)
data$Vcore <- (56-data$cm_to_top_box)*20*30 # 0.06 m2
data$shell_perc <- data$Vshell_total_ml/data$Vcore * 100
View(data)
## packages ########################################################################################
library(ggplot2)
## Maps and spatial
library("sf")
#library(sp) # replaced by sf/terra
library(terra)
library(viridis)
library(tidyverse)
## working directory ###############################################################################
setwd("~/1. PhD/Ch4 TraSSh - shells and sand/Data/Field data")  # "~/1. PhD/Ch3 Shells and sand/Data"
## load data
shells <- read.csv("202108_shellcontent_BD.csv")
fieldform <- read.csv("MOSAIC_2021_TRAILS_02_BD_2-8aug.csv")
# prepare data #####################################################################################
#shells <- shells[order(shells$Poskey), ]
fieldform$Shells[is.na(fieldform$Shells)] <- 0
colnames(fieldform)[10] <- "Poskey"
#data <- cbind(fieldform[which(is.element(fieldform$Worms, shells$Poskey)),10:44], shells)
data21 <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
# prepare data 2022 #####################################################################################
shells <- read.csv("TRAILS2202_2205_Vshells.csv")
fieldform <- read.csv("TRAILS2202_2205_2210_boxcores.csv")
#shells <- shells[order(shells$Poskey), ]
colnames(fieldform)[10] <- "Poskey"
fieldform$Shells[is.na(fieldform$Shells)] <- 0
#data <- cbind(fieldform[which(is.element(fieldform$Worms, shells$Poskey)),10:44], shells)
data22 <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
data <- rbind(data21, data22)
data <- data21
View(data)
## working directory and data #####################################################################
setwd("~/1. PhD/Ch4 TraSSh - shells and sand/Paper/Data and code/data/DataRaw/Field/")
#setwd("/data/DataRaw/Field/")
## load data
shells <- read.csv("202108_shellcontent_BD.csv")
fieldform <- read.csv("MOSAIC_2021_TRAILS_02_BD_2-8aug.csv")
# prepare data #####################################################################################]
fieldform$Shells[is.na(fieldform$Shells)] <- 0
colnames(fieldform)[10] <- "Poskey"
data21 <- left_join(fieldform, shells %>% dplyr::select(Poskey, Vshell_total_ml), by = "Poskey") %>%
dplyr::select(Poskey, Month, Year, Nbg, Nbmin, olg, olmin, Shells, cm_to_top_box, Vshell_total_ml)
data <- data21
# calculate shell fraction (volume %)
data$Vcore <- (56-data$cm_to_top_box)*20*30 # 0.06 m2
data$shell_perc <- data$Vshell_total_ml/data$Vcore * 100
# mapping#####################################################################################
#transform coordinates to XY
data$Y = as.numeric(data$Nbg)+(as.numeric(data$Nbmin)/60)
data$X = as.numeric(data$olg)+(as.numeric(data$olmin)/60)
data_sp <- data[!is.na(data$X),]
data_sp <- st_as_sf(data_sp, coords = c("X", "Y"), crs = 4326)# convert to spatial object
my_dpi = 500
w_full = 3740
# page size in inches for png
w_full = w_full/my_dpi
p1 <- ggplot(data_sp %>% filter(!is.na(shell_perc))) + #env_select
#geom_sf(data=it, linewidth=0.1)+ geom_sf(data=ws) +
geom_sf(aes(colour=shell_perc), alpha = 0.75, size=2)+ #aes(colour=seasons),
lims(x=c(5.45, 5.95), y=c(53.31, 53.51)) +
scale_colour_viridis(name=expression(paste(italic(P[shell]), " (%)"))) +
#scale_colour_manual(breaks=c('winter', 'spring', 'summer', 'autumn'),values = cols) + #breaks=c("23884","22606", "23885", "22621")
#facet_wrap( ~Year) +
theme(legend.position = "bottom", legend.text = element_text((size=35))) +
theme_bw(base_size = 12)
p1
