This R tutorial makes it possible to create bathymetry maps for the desired location. Functions are used to extract bathymetry data directly from the EMODNET WCS. By specifying the maximum and minimum longitude and latitude, any area can be. Here we will focus on the Belgian Part of the North Sea. This script was slightly adapted from the EMODNET tutorial for the LifeWatch Datacloud.

Install packages. When working on a Windows, Rtools might have to be installed in order to download packages.

1. EMODnet Bathymetry for the BPNS

1.1. Access raster data through WCS

Start creating bathymetry map for the BPNS (Analogous to the general EMODnet R tutorial). Bathymetry data can be downloaded from the EMODnet Bathymetry WCS for the BPNS

# Fill in the maximum - and minimum
# longitude & latitude coordinates for
# xmax, xmin, ymax & ymin.  For the
# BPNS these can be found on:
# https://www.marineregions.org/gazetteer.php?p=details&id=3293

# Define a function to read in raster
# data from the EMODnet bathymetry WCS
getbathymetry <- function(name = "emodnet:mean", 
    resolution = "0.2km", xmin = 2, xmax = 3.5, 
    ymin = 51, ymax = 52) {
    bbox <- paste(xmin, ymin, xmax, ymax, 
        sep = ",")
    
    con <- paste("https://ows.emodnet-bathymetry.eu/wcs?service=
    wcs&version=1.0.0&request=getcoverage&coverage=", 
        name, "&crs=EPSG:4326&BBOX=", bbox, 
        "&format=image/tiff&interpolation=nearest&resx=0.00208333&resy=
      0.00208333", 
        sep = "")
    
    print(con)
    
    stop
    nomfich <- paste(name, "img.tiff", 
        sep = "_")
    nomfich <- tempfile(nomfich)
    download(con, nomfich, quiet = TRUE, 
        mode = "wb")
    img <- raster(nomfich)
    img[img == 0] <- NA
    img[img < 0] <- 0
    names(img) <- paste(name)
    return(img)
}

# Get the bathymetry data for the BPNS
bathy_img <- getbathymetry(name = "emodnet:mean", 
    resolution = "0.2km", xmin = 2, xmax = 3.5, 
    ymin = 51, ymax = 52)

bathy <- as.data.frame(as(bathy_img, "SpatialPixelsDataFrame"))

1.2. Creating maps

map <- ggplot(aes(x = x, y = y, z = emodnet.mean), 
    data = bathy) + geom_raster(data = bathy, 
    aes(fill = emodnet.mean)) + scale_fill_gradient(low = "white", 
    high = "darkblue", name = "Depth (m)") + 
    coord_quickmap(xlim = range(xmin, xmax), 
        ylim = range(ymin, ymax)) + ggtitle("Bathymetry Belgian Part of the North Sea") + 
    xlab("Longitude") + ylab("Latitude") + 
    theme_bw() + theme(plot.title = element_text(hjust = 0.5))

map

Overlay the bathymetry map with the EEZ shapefile of the BPNS

map_BPNS <- ggplot() + theme_bw() + theme(panel.grid.minor.y = element_blank(), 
    panel.grid.minor.x = element_blank()) + 
    geom_raster(data = bathy, aes(x = x, 
        y = y, fill = emodnet.mean), alpha = 0.75) + 
    scale_fill_gradient(low = "white", 
        high = "darkblue", name = "Depth (m)") + 
    geom_polygon(data = BPNS_EEZ_fort, 
        aes(x = long, y = lat, group = group, 
            fill = "BPNS_EEZ"), colour = "black", 
        fill = "blue", alpha = 0.1) + coord_quickmap(xlim = range(xmin, 
    xmax), ylim = range(ymin, ymax)) + 
    ggtitle("Bathymetry Belgian 
    part of the North Sea") + 
    xlab("Longitude") + ylab("Latitude") + 
    theme(plot.title = element_text(hjust = 0.5))

map_BPNS

Calculate and plot isobaths for the BPNS

map_iso <- map + geom_contour(aes(colour = ..level..), 
    color = "gray30", na.rm = T, show.legend = T)

map_isobath <- direct.label(map_iso, list("far.from.others.borders", 
    "calc.boxes", "enlarge.box", hjust = 1, 
    vjust = 1, box.color = NA, fill = "transparent", 
    "draw.rects"))

map_iso

map_isobath

1.3. Creating 3D visualizations

# Use the persp() function to create 3D
# bathymetry maps
bathymetry.col <- function(n, alpha = 1) {
    if ((n <- as.integer(n[1L])) > 0) {
        j <- n%/%3
        k <- n%/%3
        i <- n - j - k
        c(if (i > 0) hsv(h = seq.int(from = 40/60, 
            to = 40/60, length.out = i), 
            alpha = alpha), if (j > 0) hsv(h = seq.int(from = 40/60, 
            to = 40/60, length.out = j), 
            alpha = alpha), if (k > 0) hsv(h = seq.int(from = 40/60, 
            to = 40/60, length.out = k), 
            alpha = alpha, s = seq.int(from = 1, 
                to = 0.1, length.out = k), 
            v = 1))
    } else {
        character()
    }
}
bathymetry3d <- -1 * bathy_img
bathymetry3d[is.na(bathymetry3d)] <- 0
bathymetry3d[bathymetry3d > 0] <- 0

summary.elev <- summary(bathymetry3d)
setBaseElev <- function(x) {
    x[is.na(x)] <- as.numeric(summary.elev[1])
    return(x)
}

elev3D <- calc(bathymetry3d, setBaseElev)
zBathy <- round(as.matrix(elev3D), 1)
y <- dim(zBathy)[1]
x <- dim(zBathy)[2]

nlat <- y
nlon <- x

nrz <- nrow(zBathy)
ncz <- ncol(zBathy)
zfacet <- (zBathy[-1, -1] + zBathy[-1, 
    -ncz] + zBathy[-nrz, -1] + zBathy[-nrz, 
    -ncz])/4
nbcol <- 254
color <- c("grey", bathymetry.col(nbcol))
facetcol <- cut(zfacet, nbcol + 1)

steplat <- (ymax - ymin)/nlat
steplon <- (xmax - xmin)/nlon

latitude <- numeric(nlat)
for (i in 1:nlat) {
    latitude[i] <- ymin + (i + 0.5) * steplat
}

longitude <- numeric(nlon)
for (i in 1:nlon) {
    longitude[i] <- xmin + (i + 0.5) * 
        steplon
}

img3d <- persp(latitude, longitude, z = zBathy, 
    theta = 90, phi = 45, col = color[facetcol], 
    scale = T, expand = 0.75, ltheta = 75, 
    shade = 0.75, border = NA, box = T, 
    ticktype = "detailed", xlab = "", ylab = "Longitude", 
    zlab = "Bathymetry (m)")

By using the rgl package, tha map can be made interactive

z <- -zBathy
x <- (1:nrow(z))
y <- (1:ncol(z))
zlim <- range(z)
zlen <- zlim[2] - zlim[1] + 1

# colorlut <-
# rev(terrain.colors(zlen,alpha=0))
colorlut <- rev(topo.colors(zlen, alpha = 0))

# colorlut<-rev(topo.colors(220))
col <- colorlut[z - zlim[1] + 1]
open3d()
rgl.surface(x, y, z, color = col)

R Tutorial: Working with EMODnet data: Belgian Part of the North Sea

1. EMODnet Bathymetry for the BPNS

1.1. Access raster data through WCS

1.2. Creating maps

1.3. Creating 3D visualizations