#' Make a depth stratified survey grid
#'
#' This function sets up a depth stratified survey grid. A simple gradient in depth
#' is simulated using [`stats::spline`] (default) with a shallow portion, shelf and
#' deep portion. Adding covariance to the depth simulation is an option.
#'
#' @param x_range Range (min x, max x) in x dimension in km
#' @param y_range Range (min y, max y) in y dimension in km
#' @param res Resolution, in km, of the grid cells
#' @param shelf_depth Approximate depth of the shelf in m
#' @param shelf_width Approximate width of the shelf in km
#' @param depth_range Range (min depth, max depth) in depth in m
#' @param n_div Number of divisions to include
#' @param strat_breaks Define strata given these depth breaks
#' @param strat_splits Number of times to horizontally split strat (i.e. easy way to increase the number of strata)
#' @param method Use a "spline", "loess" or "bezier" to generate a smooth gradient or simply use "linear" interpolation?
#'
#' @return Returns a stars object with 2 dimensions (x and y) and 4 attributes (depth, cell, division, strat).
#'
#' @seealso [`survey_grid`]
#'
#' @export
#'
#' @examples
#'
#' r <- make_grid(res = c(10, 10))
#' plot(r)
#'
#' p <- sf::st_as_sf(r["strat"], as_points = FALSE, merge = TRUE)
#' plot(p)
#'
#' @import sf
#' @import stars
#'
make_grid <- function(x_range = c(-140, 140), y_range = c(-140, 140),
res = c(3.5, 3.5), shelf_depth = 200,
shelf_width = 100, depth_range = c(0, 1000),
n_div = 1, strat_breaks = seq(0, 1000, by = 40),
strat_splits = 2, method = "spline") {
cell <- NULL
## set-up raster
utm_proj <- "+proj=utm +zone=21 +datum=WGS84 +units=km +no_defs"
bbox <- sf::st_bbox(c(xmin = x_range[1], ymin = y_range[1],
xmax = x_range[2], ymax = y_range[2]),
crs = utm_proj)
r <- stars::st_as_stars(bbox, dy = res[1], dx = res[2])
xy <- as.data.frame(r)
names(xy) <- c("x", "y", "depth")
## simulate depth
sx <- c(x_range[1], -shelf_width, -shelf_width / 2,
0, shelf_width / 2, shelf_width, x_range[2])
sy <- c(depth_range[1], rep(shelf_depth, 5), depth_range[2])
if (method == "loess") {
lo <- stats::loess(sy ~ sx)
px <- seq(min(sx), max(sx), length.out = 100)
py <- stats::predict(lo, data.frame(sx = xy$x))
s <- list(x = xy$x, y = py)
}
if (method == "spline") {
s <- stats::spline(sx, sy, xout = xy$x)
}
if (method == "linear") {
s <- stats::approx(sx, sy, xout = xy$x)
}
if (method == "bezier") {
if (requireNamespace("bezier", quietly = TRUE)) {
t <- seq(0, 1, length = 100)
p <- cbind(sx, sy)
s <- bezier::bezier(t = t, p = p)
s <- stats::approx(s[, 1], s[, 2], xout = xy$x)
} else {
stop("The bezier package is needed for to use the bezier method in make_grid. Please install it.", call. = FALSE)
}
}
depth <- s$y
depth[depth < depth_range[1]] <- depth_range[1] + 1 # impose depth range
depth[depth > depth_range[2]] <- depth_range[2] - 1
depth <- round(depth) # 1 m res ought to be good
xy$depth <- depth
## add cell number
xy$cell <- seq(nrow(xy))
## define divisions and/or splits
xy$division <- xy$split <- 1
if (n_div > 1) {
xy$division <- as.numeric(cut(xy$y, n_div))
}
if (strat_splits > 1) {
xy$split <- as.numeric(cut(xy$y, strat_splits))
}
## define strata and use unique labels along the x-axis
## (i.e. don't duplicate strat labels for areas with the same depth
## that are not adjacent to each other)
xy$strat <- as.numeric(cut(xy$depth, strat_breaks))
xy <- xy[order(xy$x), ]
rl <- rle(xy$strat)$length
xy$strat <- rep(seq_along(rl), rl)
## ensure strata within a division or a split have different numbers
div_id <- xy$division * 10 ^ (max(nchar(xy$split)) + max(nchar(xy$strat)))
split_id <- xy$split * 10 ^ max(nchar(xy$strat))
xy$strat <- as.numeric(factor(div_id + split_id + xy$strat))
## convert xyz data to a raster
xy$split <- NULL # drop - was useful for defining strata
r <- stars::st_as_stars(xy, crs = utm_proj)
r
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.