require(raster)
require(RColorBrewer)
require(magick)
require(maps)

#Coastline
setwd("/SET PATH TO DIRECTORY/Europe - Coastlines")
cstl <- shapefile("ne_10m_admin_0_countries")
setwd("/SET PATH TO DIRECTORY/Europe - Coastlines")
lake <- shapefile("ne_10m_lakes_europe")
lake <- lake[lake$WPEU10MV2_ %in% c(182, 188),]
lake <- spTransform(lake, "+proj=longlat +datum=WGS84")

setwd("/SET PATH TO DIRECTORY/Europe - Coastlines")
bat <- raster("bat_fine")
bat <- crop(bat, extent(c(2, 7.5, 51.3, 55.0)))
w <- values(bat)[!is.na(values(bat))]
values(bat)[!is.na(values(bat))] <- -w

setwd("/SET PATH TO DIRECTORY/Europe - Coastlines")
sw <- raster("Sandwave.grd")
values(sw)[!is.na(values(sw))] <- -values(sw)[!is.na(values(sw))]
sta <- read.csv("sta.csv", h = T)
spl <- read.csv("Samples.csv", h = T)
trans <- read.csv("Transect_INFO.csv", h = T, row.names = 1)

pdf("raster1.pdf", width = 10, height = 11)
par(mar = c(3, 5, 1, 6), las = 1, cex.axis = 1.4)
plot(c(2.5, 6.5), c(51.5, 54.5), typ = "n", bty = "n", xaxt = "n", yaxt = "n",
     xlab = "", ylab = "")
col <- colorRampPalette(brewer.pal(9,"Blues"))(100)
plot(bat, col = col, bty = "n", legend.width = 1.5, legend.shrink = 0.5, add = T)
lev <- c(10, 26, 30, 40, 45, 50, 55, 60, 80, 100, 120)
contour(bat, levels = lev, lwd = 1, labcex = 1, col = "gray50", drawlabels = F, add = T)
plot(cstl, add = T)
plot(lake, add = T)
points(4.485076, 53.19163, pch = "+", cex = 4)
text(4.70, 54.30, "Oyster Ground", font = 1, cex = 3, col = 1, srt = 0)
text(4.50, 53.70, "Frisian Front", font = 1, cex = 2, col = 1, srt = 13)
text(3.27, 52.80, "Southern", font = 1, cex = 3, col = 1, srt = 0)
text(3.27, 52.60, "Bight", font = 1, cex = 3, col = 1, srt = 0)
text(5.40, 53.30, "Wadden Sea", font = 1, cex = 2, col = 1, srt = 30)
text(5.45, 52.10, "The Netherlands", font = 1, cex = 2, col = 1, srt = 0)
map.scale(4.30, 51.75, relwidth = 0.2, ratio = F, cex = 1.5, pos = 1)
mtext(side = 1, at = seq(from = 3, to = 6, by = 1), line = 1, cex = 1.6,
      paste(seq(from = 3, to = 6, by = 1), rep("°E", 4), sep = ""), las = 1)
mtext(side = 2, at = seq(from = 52, to = 54, by = 1), line = 0.6,
      paste(seq(from = 52, to = 54, by = 1), rep("°N", 3), sep = ""), cex = 1.6)
abline(v = c(seq(from = 3, to = 6, by = 1)), col = "grey")
abline(h = seq(from = 52, to = 54, by = 1), col = "grey")
box()
title(main = "a", adj = 0.95, font.main = 1, cex.main = 4, line = -2.5)
dev.off()

pdf("raster2.pdf", width = 10, height = 12)
par(mar = c(3, 5, 1, 6), las = 1, cex.axis = 1.6)
plot(c(4.470, 4.500), c(53.180, 53.205), typ = "n", bty = "n",
     xlab = "", ylab = "")
col = colorRampPalette(c("gold", "gold1", "gold3", "gold4", 1))(100)
plot(sw, col = col, legend.width = 1.5, legend.shrink = 0.5, add = T)
points(sta[,c("Longitude", "Latitude")], cex = 1.5)
map.scale(4.49, 53.185, 0.2, relwidth = 0.2, ratio = F, cex = 1.5, pos = 1)
title(main = "b", adj = 0.95, font.main = 1, cex.main = 4, line = -2.5)
dev.off()

r1 <- image_read_pdf("/SET PATH TO DIRECTORY/Europe - Coastlines/raster1.pdf")
r2 <- image_read_pdf("/SET PATH TO DIRECTORY/Europe - Coastlines/raster2.pdf")

pdf("Fig1.pdf", width = 20, height = 15)
layout(matrix(c(1, 2, 3, 3), nc = 2, byrow = T), width = c(1, 1), height = c(1, 0.4))
par(mar = rep(0, 4))
plot(r1)
plot(r2)
par(mar = c(6, 8, 1, 4), las = 1, cex.axis = 2, mgp = c(3, 1.7, 0))
plot(trans$dis, trans$dep, type = "l", bty = "n", xlim = c(0, 1200), ylim = c(-32.5, -27.5),
     xlab = "", ylab = "", yaxt = "n")
w <- read.csv("disdep_box.csv", h = T)
#points(spl[c("dis", "dep")], pch = 21, col = "white", bg = 1, cex = 3)
col = as.character(w$hab)
col[col == "G"] <- "gold3"
col[col == "C"] <- "gold"
col[col == "S"] <- "gold4"
col[col == "T"] <- 1
points(w[c("dis", "dep")], pch = 21, col = "white", bg = col, cex = 3)
mtext(side = 1, "Distance (m)", cex = 2.2, line = 4.2)
axis(side = 2, seq(-28, -32, -2), labels = seq(28, 32, 2))
mtext(side = 2, "Depth (m)", cex = 2.2, line = 4.5, las = 0)
title(main = "c", adj = 0.93, font.main = 1, cex.main = 4, line = -2)
legend(1150, -29, pch = 21, pt.cex = 3, col = "white",
       pt.bg = c("gold", "gold3", "gold4", 1), legend = c("C", "G", "S", "T"),
       bty = "n", cex = 1.5)
dev.off()
