The code for the bathymetry function I have had knocking around in my head has been brought into existence. It is still pretty rough, but it runs. It could probably be sped up a decent amount as running this on ~600 lines of data takes about 12 minutes… Not great. Besides that, it’s also completely dependent on a coastline shapefile and bathymetry .asc file saved on my computer. Meaning one can only calculate shelf width and and slope if it is within the boundaries of the southern African coast. The ideal situation would be for this to be able to run on any place in the world and would use downloading on demand to retrieve the bathy data. Or perhaps I could just get the bathy data for the whole world and save that as an .RData file… That could work actually. Hmm, something to look into.

Anyway, the code shows how to calculate the width of the continental shelf (calculated here as everything above -200m) at a certain point along the coast. This is done by drawing a shore normal transect from the point. Besides shelf width, this function also calculates how wide the shelf slope is (-200 to -1000m) and uses this to come to the angle of the continental shelf. This last bit I am not 100% sure is calculated correctly as the angles are exceedingly small. With the largest being only 0.2 degrees whereas the global average is 4.0 degrees. Though perhaps I am calculating this differently.

shelfWidthSlope <- function(site){
# Find the site on the coastline and it's nearest neighbour points
coords <- data.frame(lon = site$lon, lat = site$lat)
coords2 <- knnx.index(southern_africa_coast[,1:2], as.matrix(coords), k = 1)
coords3 <- data.frame(site = site$site, southern_africa_coast
<pre></pre>
) # This line controls how wide the spread of points used to determine the shore-normal transect is. A wider spread gives a more even, less sensitive transect
coords3 <- coords3[2:1,1:3]
# Define the shore normal transect bearing
heading <- earth.bear(coords3[1,2], coords3[1,3], coords3[2,2], coords3[2,3]) + 90
if(heading >= 360){
heading <- heading-360
} else {
heading <- heading
}
# Find the depths along this bearing for 300km
distances <- seq(from = 0, to = 300000, by = 1000)
distances2 <- destPoint(p = coords, b = heading, d = distances)
sitesIdx <- knnx.index(bathy[,1:2], as.matrix(distances2), k = 1)
bathy2 <- bathy[sitesIdx,]
bathy2 <- bathy2[complete.cases(bathy2$depth),]
if(length(bathy2$depth) < 2){
info <- data.frame(site = site$site, coords, shelf = NA, slope = NA, depth_diff = NA, angle = NA)
} else {
bathy3 <- data.frame(sapply(bathy2[1:2], deg2rad), depth = bathy2$depth)
# Find the point at which the depth first hits/ exceeds 200m
shelf_end <- bathy3[bathy3$depth <= -200,][1,]
# Calculate how far it takes to reach the 200m isobath
shelf_width <- gcd.hf(bathy3$lon[1], bathy3$lat[1], shelf_end$lon, shelf_end$lat)
# Find the point at which the depth first hits/ exceeds 1000m
slope_end <- bathy3[bathy3$depth <= -1000,][1,]
# Calculate how far it takes to reach the 1000m isobath
slope_width <- gcd.hf(shelf_end$lon, shelf_end$lat, slope_end$lon, slope_end$lat)
# From this calculate the shelf slope angle
slope_depth <- abs(slope_end$depth - shelf_end$depth)
slope_angle <- tan(slope_width/slope_depth)
# Put it all together
info <- data.frame(site = site$site, coords, shelf = shelf_width, slope = slope_width, depth_diff = slope_depth, angle = slope_angle)
}
return(info)
}