Nothing
R/samplingDensity.R
In rLakeHabitat: Interpolate Bathymetry and Quantify Physical Aquatic Habitat
Defines functions
samplingDensity
samplingDensity
samplingDensity
Documented in
samplingDensity
#' Pre-Survey Power Analysis for Bathymetric Sampling Design
#'
#' Before any depth data are collected, estimate how far apart survey transects can be spaced and still produce an accurate bathymetric DEM. Because no real depth data exist yet, this works as a simulation study: a plausible "true" bathymetry is generated for the waterbody, a survey at each candidate transect spacing is simulated by sampling that truth, a DEM is reconstructed from those simulated samples (via interpBathy), and the reconstruction is compared back against the generated DEM. This is repeated many times per spacing to average out randomness, and across all candidate spacings, to build a curve of expected DEM accuracy vs. transect spacing - along with a recommended spacing and a map of what that survey design looks like on the actual waterbody outline.
#'
#' @param outline shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read.
#' @param max_depth expected maximum depth of the waterbody, in meters. This drives the synthetic bathymetry and has a large effect on results.
#' @param spacings numeric vector of candidate transect spacings to test, in meters. Default = NULL, in which case a sequence scaled to the waterbody's size is generated automatically (2%-40% of the lake's characteristic length).
#' @param n_sim number of independent "true" bathymetries to simulate and test each spacing against, default = 10. Higher values give a more stable accuracy estimate at the cost of runtime.
#' @param shape numeric value controlling the synthetic depth profile's shape: depth = max_depth * (relative distance to shore) ^ shape. shape = 1 (default) gives a linear, cone-like profile; shape < 1 gives a broader, flatter deep basin; shape > 1 gives a narrower, steeper-sided basin.
#' @param truth_model character variogram model used to generate the synthetic bathymetry's random roughness (see gstat::vgm options; "Sph", "Exp", "Gau", or "Mat"), default = "Exp".
#' @param truth_range numeric range parameter (in meters) for the synthetic roughness's spatial autocorrelation, default = NULL, in which case it is set to one quarter of the waterbody's characteristic length (sqrt(area)).
#' @param truth_sill numeric sill (variance, in squared meters) for the synthetic roughness, default = NULL, in which case it is set to (max_depth/5)^2.
#' @param truth_nugget numeric nugget for the synthetic roughness variogram, default = 0.
#' @param along_track_interval numeric spacing (in meters) at which simulated samples are drawn along each transect line, representing how frequently a depth is recorded. Default = NULL, in which case it is set to 1/10th of the smallest tested spacing.
#' @param orientation numeric angle, in degrees (0 = along the x-axis/east, increasing counterclockwise), giving the orientation of transect lines. Default = NULL, in which case the orientation is chosen automatically to align with the waterbody's long axis.
#' @param res numeric DEM cell resolution (in meters) used both for generating the synthetic truth surface and for reconstructing the DEM from simulated samples. Default = NULL, in which case it is set to 1/5th of the smallest tested spacing. This should generally be left fine relative to 'spacings'.
#' @param method character describing method of interpolation used to reconstruct DEMs, "IDW" or "OK". Default = "IDW" (recommended for this use, since it is run many times and OK's variogram fitting adds considerable runtime).
#' @param nmax numeric value describing number of neighbors used in interpolation, default = 20
#' @param idp numeric value describing inverse distance power value for IDW interpolation
#' @param model character describing type of model used in Ordinary Kriging, options include 'Sph', 'Exp', 'Gau', 'Mat', default = 'Sph'
#' @param psill numeric value describing the partial sill value for OK interpolation, default = NULL
#' @param range numeric describing distance beyond which there is no spatial correlation in Ordinary Kriging models, default = NULL
#' @param nugget numeric describing variance at zero distance in Ordinary Kriging models, default = 0
#' @param kappa numeric value describing model smoothness, default = NULL
#' @param trend_order numeric value (1 or 2) giving the order of the polynomial trend surface fit for Universal Kriging. 1 = linear trend (z ~ x + y), 2 = quadratic trend. Default = 1.
#' @param zero_threshold numeric proportion (0-1) of the waterbody's surface area that must interpolate to exactly 0 before the automatic zero re-interpolation pass runs, default = 0.05 (5%). A handful of scattered zero cells won't trigger it; a large contiguous block collapsing to zero (typically an interpolation artifact, often from the shoreline zero ring dominating a narrow bay or inlet) will.
#' @param tolerance numeric value (proportion) used to pick the recommended spacing: the coarsest tested spacing whose mean RMSE is still within 'tolerance' of the best (finest-spacing) RMSE. Default = 0.1 (10%).
#' @param n_truth_examples numeric value describing how many of the n_sim synthetic "true" bathymetries to save and return/plot, so you can visually check whether the assumed depth profile and roughness look like your lake before trusting the rest of the results. Default = 3 (capped at n_sim).
#' @param plot logical: should diagnostic plots (accuracy-vs-spacing curve, and outline + recommended transects map) be drawn? Default = TRUE.
#' @param seed optional numeric value used to seed the random number generator, for reproducible results across runs. Default = NULL (not seeded).
#' @details
#' This is a power analysis under assumed, not observed, bathymetry - it can only be as realistic as the 'max_depth', 'shape', 'truth_range', and 'truth_sill' inputs. Once real depth data have been collected, use crossValidate(), interpBathy(), and optimizeParams() on the actual data to check whether observed accuracy matches what was predicted here, and adjust future survey effort accordingly.
#' Runtime scales with n_sim x length(spacings) x (cost of one synthetic-surface simulation + one interpBathy call). For a first look, consider a smaller n_sim (e.g. 5) and a coarser 'res' before committing to a longer run.
#' @return a list with:
#' \describe{
#' \item{results}{a data frame of spacing, mean RMSE, and SD of RMSE across simulations}
#' \item{recommended_spacing}{the recommended transect spacing, in meters}
#' \item{transects}{a SpatVector of the recommended transect lines, in the original CRS of 'outline'}
#' \item{total_transect_length}{total length of the recommended transects, in meters}
#' \item{truth_examples}{a multi-layer SpatRaster (one layer per saved example) of synthetic "true" bathymetries used
#' during the simulation, in the original CRS of 'outline', for visually sanity-checking the assumed depth profile}
#' }
#' @author Tristan Blechinger, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @import dplyr
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @import gstat
#' @examples
#' \donttest{
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' samplingDensity(outline, max_depth = 40, n_sim = 5, seed = 123)}
samplingDensity <- function(outline, max_depth, spacings = NULL, n_sim = 10, shape = 1,
truth_model = "Exp", truth_range = NULL, truth_sill = NULL, truth_nugget = 0,
along_track_interval = NULL, orientation = NULL, res = NULL,
method = "IDW", nmax = 20, idp = 2, model = "Sph", psill = NULL, range = NULL,
nugget = 0, kappa = NULL, trend_order = 1, zero_threshold = 0.5, tolerance = 0.1, n_truth_examples = 3,
plot = TRUE, seed = NULL){
if(!is.null(seed)){
if(!is.numeric(seed))
stop("seed must be numeric")
set.seed(seed)
}
# transform outline shapefile into vector
if(!inherits(outline, "SpatVector")){
if(inherits(outline, "sf") && requireNamespace("sf", quietly = TRUE)){
outline <- sf::st_zm(outline, drop = TRUE, what = "ZM")
}
outline <- terra::vect(outline)
}
else{
outline <- outline
}
# checks
if(!inherits(outline, "SpatVector"))
stop("outline is not a SpatVector or cannot be transformed")
if(missing(max_depth) || !is.numeric(max_depth) || max_depth <= 0)
stop("max_depth must be specified as a positive numeric value")
if(!method %in% c("IDW", "OK"))
stop("method misspecified. Please choose either 'IDW' or 'OK'")
if(!is.numeric(tolerance) || tolerance < 0)
stop("tolerance must be a non-negative numeric value")
#store the original CRS so the recommended transects can be returned in it
original_crs <- terra::crs(outline)
test_crs <- terra::crs(outline)
if(is.na(test_crs) || test_crs == ""){
stop("CRS of 'outline' is unable to be defined.")
}
# Function to determine the best UTM CRS for a given vector
get_best_utm <- function(outline) {
centroid <- terra::centroids(outline)
lon <- terra::crds(centroid)[1]
lat <- terra::crds(centroid)[2]
utm_zone <- base::floor((lon + 180) / 6) + 1
hemisphere <- base::ifelse(lat >= 0, 32600, 32700)
epsg_code <- hemisphere + utm_zone
return(terra::crs(paste0("EPSG:", epsg_code)))
}
# run analysis in meters (UTM)
if(terra::is.lonlat(outline)){
best_crs <- get_best_utm(outline)
outline <- terra::project(outline, best_crs)
}
message("Simulating in CRS: ", terra::crs(outline, describe = TRUE)$name)
# characteristic length scale of the waterbody, used to scale several defaults
area <- as.numeric(terra::expanse(outline))
L <- sqrt(area)
if(is.null(spacings)){
spacings <- round(L * c(.02, .05, .1, .2, .4))
spacings <- unique(spacings[spacings > 0])
}
if(!is.numeric(spacings) || any(spacings <= 0))
stop("spacings must be positive numeric values")
spacings <- sort(spacings)
if(is.null(res)){
res <- max(1, min(spacings) / 5)
}
if(is.null(along_track_interval)){
along_track_interval <- max(5, min(spacings) / 10)
}
if(is.null(truth_range)){
truth_range <- L / 4
}
if(is.null(truth_sill)){
truth_sill <- (max_depth / 5)^2
}
#### helper: build a masked raster template over the waterbody at a given resolution
build_grid <- function(outline, res){
ext_o <- terra::ext(outline)
ncol_r <- ceiling((ext_o$xmax - ext_o$xmin) / res)
nrow_r <- ceiling((ext_o$ymax - ext_o$ymin) / res)
empty <- terra::rast(ext_o, ncol = ncol_r, nrow = nrow_r, crs = terra::crs(outline))
ras <- terra::rasterize(outline, empty)
grid <- terra::mask(ras, outline)
return(grid)
}
#### helper: rotate xy points by theta (radians, counterclockwise) about a center point
rotate_pts <- function(xy, theta, center){
xy_c <- sweep(xy, 2, center, FUN = "-")
rot_mat <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), nrow = 2)
rotated <- xy_c %*% rot_mat
sweep(rotated, 2, center, FUN = "+")
}
#### helper: build parallel transect lines across outline at a given spacing/orientation,
#### clipped to the outline polygon
make_transects <- function(outline, spacing, orientation_deg = NULL){
verts <- terra::crds(outline)
center <- colMeans(verts)
if(is.null(orientation_deg)){
pc <- stats::prcomp(scale(verts, scale = FALSE))
theta <- atan2(pc$rotation[2, 1], pc$rotation[1, 1])
} else {
theta <- orientation_deg * pi / 180
}
rot_verts <- rotate_pts(verts, -theta, center)
xr <- range(rot_verts[, 1])
yr <- range(rot_verts[, 2])
pad <- spacing
y_lines <- seq(yr[1] - pad, yr[2] + pad, by = spacing)
line_list <- lapply(y_lines, function(yy){
pts_rot <- rbind(c(xr[1] - pad, yy), c(xr[2] + pad, yy))
pts_real <- rotate_pts(pts_rot, theta, center)
terra::vect(pts_real, type = "lines", crs = terra::crs(outline))
})
lines_all <- do.call(rbind, line_list)
lines_clip <- terra::intersect(lines_all, outline)
return(lines_clip)
}
#### helper: generate one synthetic "true" bathymetry raster
simulate_truth <- function(){
grid <- build_grid(outline, res)
dist_to_shore <- terra::distance(grid, terra::as.lines(outline))
dist_to_shore <- terra::mask(dist_to_shore, outline)
max_dist <- terra::global(dist_to_shore, "max", na.rm = TRUE)[[1]]
bowl <- max_depth * (dist_to_shore / max_dist)^shape
# simulate the random roughness component on a coarse grid sized off its own spatial autocorrelation range, then resample to match the fine grid
coarse_res <- max(res, truth_range / 4)
coarse_grid <- build_grid(outline, coarse_res)
coarse_df <- as.data.frame(coarse_grid, xy = TRUE, na.rm = TRUE)
names(coarse_df)[1:2] <- c("x", "y")
g_dummy <- gstat::gstat(formula = z ~ 1, locations = ~x + y, dummy = TRUE, beta = 0,
model = gstat::vgm(psill = truth_sill, model = truth_model,
range = truth_range, nugget = truth_nugget),
nmax = 30)
sim <- stats::predict(g_dummy, newdata = coarse_df, nsim = 1)
roughness_coarse <- terra::rast(data.frame(x = coarse_df$x, y = coarse_df$y, roughness = sim$sim1),
type = "xyz", crs = terra::crs(outline))
roughness_fine <- terra::resample(roughness_coarse, bowl, method = "bilinear")
true_depth <- bowl + roughness_fine
true_depth <- terra::ifel(true_depth < 0, 0, true_depth)
true_depth <- terra::mask(true_depth, outline)
names(true_depth) <- "true_depth"
return(true_depth)
}
#### helper: simulate a survey at a given spacing against a given truth, reconstruct, and score
score_spacing <- function(truth_raster, spacing){
transects <- make_transects(outline, spacing, orientation)
densified <- terra::densify(transects, interval = along_track_interval)
sample_pts <- terra::as.points(densified)
vals <- terra::extract(truth_raster, sample_pts, ID = FALSE)
xy <- terra::crds(sample_pts)
sample_df <- data.frame(x = xy[, 1], y = xy[, 2], z = vals[[1]])
sample_df <- sample_df[!is.na(sample_df$z), ]
result <- tryCatch({
out <- utils::capture.output(
dem <- suppressMessages(interpBathy(outline, sample_df, x = "x", y = "y", z = "z",
zeros = FALSE, separation = along_track_interval,
res = res, method = method, nmax = nmax, idp = idp,
model = model, psill = psill, range = range,
nugget = nugget, kappa = kappa, trend_order = trend_order,
zero_threshold = zero_threshold))
)
if(method == "OK") dem <- dem[["depth"]]
dem_resampled <- terra::resample(dem, truth_raster)
sq_err <- (dem_resampled - truth_raster)^2
sqrt(terra::global(sq_err, "mean", na.rm = TRUE)[[1]])
}, error = function(e){
warning("Reconstruction failed for spacing = ", spacing, ": ", conditionMessage(e))
return(NA_real_)
})
return(result)
}
# Run simulation
n_iter <- n_sim * length(spacings)
pb <- utils::txtProgressBar(min = 0, max = n_iter, style = 3)
iter <- 0
rmse_matrix <- matrix(NA_real_, nrow = n_sim, ncol = length(spacings))
colnames(rmse_matrix) <- as.character(spacings)
n_truth_examples <- min(n_truth_examples, n_sim)
truth_list <- list()
for(i in 1:n_sim){
truth_raster <- simulate_truth()
if(i <= n_truth_examples){
truth_list[[i]] <- truth_raster
}
for(j in seq_along(spacings)){
rmse_matrix[i, j] <- score_spacing(truth_raster, spacings[j])
iter <- iter + 1
utils::setTxtProgressBar(pb, iter)
}
}
close(pb)
results <- data.frame(spacing = spacings,
mean_rmse = apply(rmse_matrix, 2, mean, na.rm = TRUE),
sd_rmse = apply(rmse_matrix, 2, stats::sd, na.rm = TRUE))
best_rmse <- min(results$mean_rmse, na.rm = TRUE)
recommended_spacing <- max(results$spacing[results$mean_rmse <= best_rmse * (1 + tolerance)])
#build the recommended transects once more, cleanly, and return them in the original CRS
transects_final <- make_transects(outline, recommended_spacing, orientation)
total_transect_length <- sum(terra::perim(transects_final))
transects_final <- terra::project(transects_final, original_crs)
outline_orig <- terra::project(outline, original_crs)
truth_examples <- do.call(c, truth_list)
names(truth_examples) <- paste0("sim_", seq_len(n_truth_examples))
truth_examples <- terra::project(truth_examples, original_crs)
if(plot){
graphics::plot(results$spacing, results$mean_rmse, type = "b", pch = 19,
xlab = "Transect spacing (m)", ylab = "RMSE (m)",
main = "Simulated DEM accuracy vs. transect spacing")
graphics::abline(v = recommended_spacing, lty = 2, col = "blue")
graphics::legend("topleft", legend = paste0("Recommended: ", recommended_spacing, " m"),
lty = 2, col = "blue", bty = "n")
terra::plot(outline_orig, main = paste0("Recommended survey design (", recommended_spacing, " m spacing)"))
terra::plot(transects_final, add = TRUE, col = "blue", lwd = 1.5)
terra::plot(truth_examples, col = grDevices::hcl.colors(100, "Blues"),
mar = c(2, 2, 2, 4))
}
return(list(results = results,
recommended_spacing = recommended_spacing,
transects = transects_final,
total_transect_length = total_transect_length,
truth_examples = truth_examples))
}
#' Pre-Survey Power Analysis for Bathymetric Sampling Design
#'
#' Before any depth data are collected, estimate how far apart survey transects can be spaced and still produce an accurate bathymetric DEM. Because no real depth data exist yet, this works as a simulation study: a plausible "true" bathymetry is generated for the waterbody, a survey at each candidate transect spacing is simulated by sampling that truth, a DEM is reconstructed from those simulated samples (via interpBathy), and the reconstruction is compared back against the generated DEM. This is repeated many times per spacing to average out randomness, and across all candidate spacings, to build a curve of expected DEM accuracy vs. transect spacing - along with a recommended spacing and a map of what that survey design looks like on the actual waterbody outline.
#'
#' @param outline shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read.
#' @param max_depth expected maximum depth of the waterbody, in meters. This drives the synthetic bathymetry and has a large effect on the results.
#' @param spacings numeric vector of candidate transect spacings to test, in meters. Default = NULL, in which case a sequence scaled to the waterbody's size is generated automatically (roughly 2%-40% of the lake's characteristic length).
#' @param n_sim number of independent synthetic "true" bathymetries to simulate and test each spacing against, default = 10.
#' @param shape numeric value controlling the synthetic depth profile's shape: depth = max_depth * (relative distance to shore) ^ shape. shape = 1 (default) gives a linear, cone-like profile; shape > 1 gives a broader, flatter deep basin; shape < 1 gives a narrower, steeper-sided basin.
#' @param truth_model character variogram model used to generate the synthetic bathymetry's random roughness (see gstat::vgm options), default = "Exp".
#' @param truth_range numeric range parameter (in meters) for the synthetic roughness's spatial autocorrelation, default = NULL, in which case it is set to one quarter of the waterbody's characteristic length (sqrt(area)).
#' @param truth_sill numeric sill (variance, in squared meters) for the synthetic roughness, default = NULL, in which case it is set to (max_depth/5)^2.
#' @param truth_nugget numeric nugget for the synthetic roughness variogram, default = 0.
#' @param along_track_interval numeric spacing (in meters) at which simulated samples are drawn along each transect line, representing how frequently a sonar unit records depth. Default = NULL, in which case it is set to 1/10th of the smallest tested spacing.
#' @param orientation numeric angle, in degrees (0 = along the x-axis/east, increasing counterclockwise), giving the orientation of transect lines. Default = NULL, in which case the orientation is chosen automatically to align with the waterbody's long axis.
#' @param res numeric DEM cell resolution (in meters) used both for generating the synthetic truth surface and for reconstructing the DEM from simulated samples. Default = NULL, in which case it is set to 1/5th of the smallest tested spacing.
#' @param method character describing method of interpolation used to reconstruct DEMs, "IDW", "OK", or "UK". Default = "IDW" (recommended for this use, since it is run many times and OK's variogram fitting adds considerable runtime).
#' @param nmax numeric value describing number of neighbors used in interpolation, default = 20
#' @param idp numeric value describing inverse distance power value for IDW interpolation
#' @param model character describing type of model used in Ordinary Kriging, options include 'Sph', 'Exp', 'Gau', 'Sta', default = 'Sph'
#' @param psill numeric value describing the partial sill value for OK interpolation, default = NULL
#' @param range numeric describing distance beyond which there is no spatial correlation in Ordinary Kriging models, default = NULL
#' @param nugget numeric describing variance at zero distance in Ordinary Kriging models, default = 0
#' @param kappa numeric value describing model smoothness, default = NULL
#' @param trend_order numeric value (1 or 2) giving the polynomial trend order for Universal Kriging ("UK" only), default = 1
#' @param zero_threshold numeric proportion (0-1) of surface area that must interpolate to exactly 0 before the automatic zero re-interpolation pass runs - passed through to interpBathy(). Default = 0.05.
#' @param tolerance numeric value (proportion) used to pick the recommended spacing: the coarsest tested spacing whose mean RMSE is still within 'tolerance' of the best (finest-spacing) RMSE. Default = 0.1 (10%).
#' @param n_truth_examples numeric value giving how many of the n_sim synthetic "true" bathymetries to save and return/plot, so you can visually check whether the assumed depth profile and roughness look like your lake before trusting the rest of the results. Default = 3 (capped at n_sim).
#' @param plot logical: should diagnostic plots (accuracy-vs-spacing curve, and outline + recommended transects map) be drawn? Default = TRUE.
#' @param seed optional numeric value used to seed the random number generator, for reproducible results across runs. Default = NULL (not seeded).
#' @details
#' This is a power analysis under assumed, not observed, bathymetry - it can only be as realistic as the 'max_depth',
#' 'shape', 'truth_range', and 'truth_sill' inputs. Once real depth data have been collected, use crossValidate() and
#' interpBathy() on the actual data to check whether observed accuracy matches what was predicted here, and adjust
#' future survey effort accordingly.
#' @return a list with:
#' \describe{
#' \item{results}{a data frame of spacing, mean RMSE, and SD of RMSE across simulations}
#' \item{recommended_spacing}{the recommended transect spacing, in meters}
#' \item{transects}{a SpatVector of the recommended transect lines, in the original CRS of 'outline'}
#' \item{total_transect_length}{total length of the recommended transects, in meters}
#' \item{truth_examples}{a multi-layer SpatRaster (one layer per saved example) of synthetic "true" bathymetries used
#' during the simulation, in the original CRS of 'outline', for visually sanity-checking the assumed depth profile}
#' }
#' @author Tristan Blechinger, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @import dplyr
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @import gstat
#' @examples
#' \donttest{
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' samplingDensity(outline, max_depth = 40, n_sim = 5, seed = 123)}
samplingDensity <- function(outline, max_depth, spacings = NULL, n_sim = 10, shape = 1,
truth_model = "Exp", truth_range = NULL, truth_sill = NULL, truth_nugget = 0,
along_track_interval = NULL, orientation = NULL, res = NULL,
method = "IDW", nmax = 20, idp = 2, model = "Sph", psill = NULL, range = NULL,
nugget = 0, kappa = NULL, trend_order = 1, zero_threshold = 0.05,
tolerance = 0.1, n_truth_examples = 3,
plot = TRUE, seed = NULL){
if(!is.null(seed)){
if(!is.numeric(seed))
stop("seed must be numeric")
set.seed(seed)
}
#transform outline shapefile into vector (accepts SpatVector, sf, or anything terra::vect() reads)
if(!inherits(outline, "SpatVector")){
if(inherits(outline, "sf") && requireNamespace("sf", quietly = TRUE)){
outline <- sf::st_zm(outline, drop = TRUE, what = "ZM")
}
outline <- terra::vect(outline)
}
else{
outline <- outline
}
#checks
if(!inherits(outline, "SpatVector"))
stop("outline is not a SpatVector or cannot be transformed")
if(missing(max_depth) || !is.numeric(max_depth) || max_depth <= 0)
stop("max_depth must be specified as a positive numeric value")
if(!method %in% c("IDW", "OK", "UK"))
stop("method misspecified. Please choose 'IDW', 'OK', or 'UK'")
if(!is.numeric(tolerance) || tolerance < 0)
stop("tolerance must be a non-negative numeric value")
#store the original CRS so the recommended transects can be returned in it
original_crs <- terra::crs(outline)
test_crs <- terra::crs(outline)
if(is.na(test_crs) || test_crs == ""){
stop("CRS of 'outline' is unable to be defined.")
}
# Function to determine the best UTM CRS for a given vector
get_best_utm <- function(outline) {
centroid <- terra::centroids(outline)
lon <- terra::crds(centroid)[1]
lat <- terra::crds(centroid)[2]
utm_zone <- base::floor((lon + 180) / 6) + 1
hemisphere <- base::ifelse(lat >= 0, 32600, 32700)
epsg_code <- hemisphere + utm_zone
return(terra::crs(paste0("EPSG:", epsg_code)))
}
#run the whole analysis in meters (UTM), same reasoning as interpBathy/crossValidate
if(terra::is.lonlat(outline)){
best_crs <- get_best_utm(outline)
outline <- terra::project(outline, best_crs)
}
message("Simulating in CRS: ", terra::crs(outline, describe = TRUE)$name)
#characteristic length scale of the waterbody, used to scale several defaults
area <- as.numeric(terra::expanse(outline))
L <- sqrt(area)
if(is.null(spacings)){
spacings <- round(L * c(.02, .05, .1, .2, .4))
spacings <- unique(spacings[spacings > 0])
}
if(!is.numeric(spacings) || any(spacings <= 0))
stop("spacings must be positive numeric values")
spacings <- sort(spacings)
if(is.null(res)){
res <- max(1, min(spacings) / 5)
}
if(is.null(along_track_interval)){
along_track_interval <- max(5, min(spacings) / 10)
}
if(is.null(truth_range)){
truth_range <- L / 4
}
if(is.null(truth_sill)){
truth_sill <- (max_depth / 5)^2
}
#### helper: build a masked raster template over the waterbody at a given resolution
build_grid <- function(outline, res){
ext_o <- terra::ext(outline)
ncol_r <- ceiling((ext_o$xmax - ext_o$xmin) / res)
nrow_r <- ceiling((ext_o$ymax - ext_o$ymin) / res)
empty <- terra::rast(ext_o, ncol = ncol_r, nrow = nrow_r, crs = terra::crs(outline))
ras <- terra::rasterize(outline, empty)
grid <- terra::mask(ras, outline)
return(grid)
}
#### helper: rotate xy points by theta (radians, counterclockwise) about a center point
rotate_pts <- function(xy, theta, center){
xy_c <- sweep(xy, 2, center, FUN = "-")
rot_mat <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), nrow = 2)
rotated <- xy_c %*% rot_mat
sweep(rotated, 2, center, FUN = "+")
}
#### helper: build parallel transect lines across outline at a given spacing/orientation,
#### clipped to the outline polygon
make_transects <- function(outline, spacing, orientation_deg = NULL){
verts <- terra::crds(outline)
center <- colMeans(verts)
if(is.null(orientation_deg)){
pc <- stats::prcomp(scale(verts, scale = FALSE))
theta <- atan2(pc$rotation[2, 1], pc$rotation[1, 1])
} else {
theta <- orientation_deg * pi / 180
}
rot_verts <- rotate_pts(verts, -theta, center)
xr <- range(rot_verts[, 1])
yr <- range(rot_verts[, 2])
pad <- spacing
y_lines <- seq(yr[1] - pad, yr[2] + pad, by = spacing)
line_list <- lapply(y_lines, function(yy){
pts_rot <- rbind(c(xr[1] - pad, yy), c(xr[2] + pad, yy))
pts_real <- rotate_pts(pts_rot, theta, center)
terra::vect(pts_real, type = "lines", crs = terra::crs(outline))
})
lines_all <- do.call(rbind, line_list)
lines_clip <- terra::intersect(lines_all, outline)
return(lines_clip)
}
#### helper: generate one synthetic "true" bathymetry raster
simulate_truth <- function(){
grid <- build_grid(outline, res)
dist_to_shore <- terra::distance(grid, terra::as.lines(outline))
dist_to_shore <- terra::mask(dist_to_shore, outline)
max_dist <- terra::global(dist_to_shore, "max", na.rm = TRUE)[[1]]
bowl <- max_depth * (dist_to_shore / max_dist)^shape
#simulate the random roughness component on a coarse grid sized off its own
#spatial autocorrelation range (truth_range), then resample up to match the
#fine bowl-shape grid. The roughness has no real structure finer than
#truth_range anyway, so this avoids running gstat's sequential Gaussian
#simulation over potentially hundreds of thousands to millions of cells on
#large waterbodies - the simulation cost grows much faster than linearly
#with cell count, so this is the main lever for keeping runtime practical.
coarse_res <- max(res, truth_range / 4)
coarse_grid <- build_grid(outline, coarse_res)
coarse_df <- as.data.frame(coarse_grid, xy = TRUE, na.rm = TRUE)
names(coarse_df)[1:2] <- c("x", "y")
g_dummy <- gstat::gstat(formula = z ~ 1, locations = ~x + y, dummy = TRUE, beta = 0,
model = gstat::vgm(psill = truth_sill, model = truth_model,
range = truth_range, nugget = truth_nugget),
nmax = 30)
sim <- stats::predict(g_dummy, newdata = coarse_df, nsim = 1)
roughness_coarse <- terra::rast(data.frame(x = coarse_df$x, y = coarse_df$y, roughness = sim$sim1),
type = "xyz", crs = terra::crs(outline))
roughness_fine <- terra::resample(roughness_coarse, bowl, method = "bilinear")
true_depth <- bowl + roughness_fine
true_depth <- terra::ifel(true_depth < 0, 0, true_depth)
true_depth <- terra::mask(true_depth, outline)
names(true_depth) <- "true_depth"
return(true_depth)
}
#### helper: simulate a survey at a given spacing against a given truth, reconstruct, and score
score_spacing <- function(truth_raster, spacing){
transects <- make_transects(outline, spacing, orientation)
densified <- terra::densify(transects, interval = along_track_interval)
sample_pts <- terra::as.points(densified)
vals <- terra::extract(truth_raster, sample_pts, ID = FALSE)
xy <- terra::crds(sample_pts)
sample_df <- data.frame(x = xy[, 1], y = xy[, 2], z = vals[[1]])
sample_df <- sample_df[!is.na(sample_df$z), ]
result <- tryCatch({
out <- utils::capture.output(
dem <- suppressMessages(interpBathy(outline, sample_df, x = "x", y = "y", z = "z",
zeros = FALSE, separation = along_track_interval,
res = res, method = method, nmax = nmax, idp = idp,
model = model, psill = psill, range = range,
nugget = nugget, kappa = kappa, trend_order = trend_order,
zero_threshold = zero_threshold))
)
if(method %in% c("OK", "UK")) dem <- dem[["depth"]]
dem_resampled <- terra::resample(dem, truth_raster)
sq_err <- (dem_resampled - truth_raster)^2
sqrt(terra::global(sq_err, "mean", na.rm = TRUE)[[1]])
}, error = function(e){
warning("Reconstruction failed for spacing = ", spacing, ": ", conditionMessage(e))
return(NA_real_)
})
return(result)
}
#### run the simulation: one truth surface per rep, tested against every spacing
n_iter <- n_sim * length(spacings)
pb <- utils::txtProgressBar(min = 0, max = n_iter, style = 3)
iter <- 0
rmse_matrix <- matrix(NA_real_, nrow = n_sim, ncol = length(spacings))
colnames(rmse_matrix) <- as.character(spacings)
n_truth_examples <- min(n_truth_examples, n_sim)
truth_list <- list()
for(i in 1:n_sim){
truth_raster <- simulate_truth()
if(i <= n_truth_examples){
truth_list[[i]] <- truth_raster
}
for(j in seq_along(spacings)){
rmse_matrix[i, j] <- score_spacing(truth_raster, spacings[j])
iter <- iter + 1
utils::setTxtProgressBar(pb, iter)
}
}
close(pb)
results <- data.frame(spacing = spacings,
mean_rmse = apply(rmse_matrix, 2, mean, na.rm = TRUE),
sd_rmse = apply(rmse_matrix, 2, stats::sd, na.rm = TRUE))
best_rmse <- min(results$mean_rmse, na.rm = TRUE)
recommended_spacing <- max(results$spacing[results$mean_rmse <= best_rmse * (1 + tolerance)])
#build the recommended transects once more, cleanly, and return them in the original CRS
transects_final <- make_transects(outline, recommended_spacing, orientation)
total_transect_length <- sum(terra::perim(transects_final))
transects_final <- terra::project(transects_final, original_crs)
outline_orig <- terra::project(outline, original_crs)
truth_examples <- do.call(c, truth_list)
names(truth_examples) <- paste0("sim_", seq_len(n_truth_examples))
truth_examples <- terra::project(truth_examples, original_crs)
if(plot){
graphics::plot(results$spacing, results$mean_rmse, type = "b", pch = 19,
xlab = "Transect spacing (m)", ylab = "RMSE (m)",
main = "Simulated DEM accuracy vs. transect spacing")
graphics::abline(v = recommended_spacing, lty = 2, col = "blue")
graphics::legend("topleft", legend = paste0("Recommended: ", recommended_spacing, " m"),
lty = 2, col = "blue", bty = "n")
terra::plot(outline_orig, main = paste0("Recommended survey design (", recommended_spacing, " m spacing)"))
terra::plot(transects_final, add = TRUE, col = "blue", lwd = 1.5)
terra::plot(truth_examples, col = grDevices::hcl.colors(100, "Blues"),
mar = c(2, 2, 2, 4))
}
return(list(results = results,
recommended_spacing = recommended_spacing,
transects = transects_final,
total_transect_length = total_transect_length,
truth_examples = truth_examples))
}
#' Pre-Survey Power Analysis for Bathymetric Sampling Design
#'
#' Before any depth data are collected, estimate how far apart survey transects can be spaced
#' and still produce an accurate bathymetric DEM. Because no real depth data exist yet, this
#' works as a simulation study: a plausible synthetic "true" bathymetry is generated for the
#' waterbody (a bowl-shaped depth profile plus spatially-correlated random roughness), a survey
#' at each candidate transect spacing is simulated by sampling that synthetic truth, a DEM is
#' reconstructed from those simulated samples (via interpBathy), and the reconstruction is
#' compared back against the known synthetic truth. This is repeated many times per spacing to
#' average out randomness, and across all candidate spacings, to build a curve of expected DEM
#' accuracy vs. transect spacing - along with a recommended spacing and a map of what that
#' survey design looks like on the actual lake outline.
#'
#' @param outline shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read.
#' @param max_depth expected maximum depth of the waterbody, in meters. This drives the synthetic bathymetry and has a large
#' effect on the results - use the best estimate available (a known max depth, an old chart, or a knowledgeable guess).
#' @param spacings numeric vector of candidate transect spacings to test, in meters. Default = NULL, in which case a
#' sequence scaled to the waterbody's size is generated automatically (roughly 2%%-40%% of the lake's characteristic length).
#' @param n_sim number of independent synthetic "true" bathymetries to simulate and test each spacing against, default = 10.
#' Higher values give a more stable accuracy estimate at the cost of runtime.
#' @param shape numeric value controlling the synthetic depth profile's shape: depth = max_depth * (relative distance to
#' shore) ^ shape. shape = 1 (default) gives a linear, cone-like profile; shape < 1 gives a broader, flatter deep basin;
#' shape > 1 gives a narrower, steeper-sided basin.
#' @param truth_model character variogram model used to generate the synthetic bathymetry's random roughness (see
#' gstat::vgm options), default = "Exp".
#' @param truth_range numeric range parameter (in meters) for the synthetic roughness's spatial autocorrelation, default =
#' NULL, in which case it is set to one quarter of the waterbody's characteristic length (sqrt(area)).
#' @param truth_sill numeric sill (variance, in squared meters) for the synthetic roughness, default = NULL, in which case
#' it is set to (max_depth/5)^2.
#' @param truth_nugget numeric nugget for the synthetic roughness variogram, default = 0.
#' @param along_track_interval numeric spacing (in meters) at which simulated samples are drawn along each transect line,
#' representing how frequently a sonar unit records depth while underway. Default = NULL, in which case it is set to
#' 1/10th of the smallest tested spacing.
#' @param orientation numeric angle, in degrees (0 = along the x-axis/east, increasing counterclockwise), giving the
#' orientation of transect lines. Default = NULL, in which case the orientation is chosen automatically to align with the
#' waterbody's long axis.
#' @param res numeric DEM cell resolution (in meters) used both for generating the synthetic truth surface and for
#' reconstructing the DEM from simulated samples. Default = NULL, in which case it is set to 1/5th of the smallest tested
#' spacing. This should generally be left fine relative to 'spacings'.
#' @param method character describing method of interpolation used to reconstruct DEMs, "IDW", "OK", or "UK". Default = "IDW"
#' (recommended for this use, since it is run many times and OK's variogram fitting adds considerable runtime).
#' @param nmax numeric value describing number of neighbors used in interpolation, default = 20
#' @param idp numeric value describing inverse distance power value for IDW interpolation
#' @param model character describing type of model used in Ordinary Kriging, options include 'Sph', 'Exp', 'Gau', 'Sta', default = 'Sph'
#' @param psill numeric value describing the partial sill value for OK interpolation, default = NULL
#' @param range numeric describing distance beyond which there is no spatial correlation in Ordinary Kriging models, default = NULL
#' @param nugget numeric describing variance at zero distance in Ordinary Kriging models, default = 0
#' @param kappa numeric value describing model smoothness, default = NULL
#' @param trend_order numeric value (1 or 2) giving the polynomial trend order for Universal Kriging ("UK" only), default = 1
#' @param zero_threshold numeric proportion (0-1) of surface area that must interpolate to exactly 0 before the
#' automatic zero re-interpolation pass runs - passed through to interpBathy(). Default = 0.05.
#' @param tolerance numeric value (proportion) used to pick the recommended spacing: the coarsest tested spacing whose mean
#' RMSE is still within 'tolerance' of the best (finest-spacing) RMSE. Default = 0.1 (10%%).
#' @param n_truth_examples numeric value giving how many of the n_sim synthetic "true" bathymetries to save and return/plot,
#' so you can visually check whether the assumed depth profile and roughness look like your lake before trusting the
#' rest of the results. Default = 3 (capped at n_sim).
#' @param plot logical: should diagnostic plots (accuracy-vs-spacing curve, and outline + recommended transects map) be
#' drawn? Default = TRUE.
#' @param seed optional numeric value used to seed the random number generator, for reproducible results across runs.
#' Default = NULL (not seeded).
#' @details
#' This is a power analysis under assumed, not observed, bathymetry - it can only be as realistic as the 'max_depth',
#' 'shape', 'truth_range', and 'truth_sill' inputs. Once real depth data have been collected, use crossValidate() and
#' interpBathy() on the actual data to check whether observed accuracy matches what was predicted here, and adjust
#' future survey effort accordingly.
#' Runtime scales with n_sim x length(spacings) x (cost of one synthetic-surface simulation + one interpBathy call).
#' For a first look, consider a smaller n_sim (e.g. 5) and a coarser 'res' before committing to a long run.
#' @return a list with:
#' \describe{
#' \item{results}{a data frame of spacing, mean RMSE, and SD of RMSE across simulations}
#' \item{recommended_spacing}{the recommended transect spacing, in meters}
#' \item{transects}{a SpatVector of the recommended transect lines, in the original CRS of 'outline'}
#' \item{total_transect_length}{total length of the recommended transects, in meters}
#' \item{truth_examples}{a multi-layer SpatRaster (one layer per saved example) of synthetic "true" bathymetries used
#' during the simulation, in the original CRS of 'outline', for visually sanity-checking the assumed depth profile}
#' }
#' @author Tristan Blechinger, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @import dplyr
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @import gstat
#' @examples
#' \donttest{
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' samplingDensity(outline, max_depth = 40, n_sim = 5, seed = 123)}
samplingDensity <- function(outline, max_depth, spacings = NULL, n_sim = 10, shape = 1,
truth_model = "Exp", truth_range = NULL, truth_sill = NULL, truth_nugget = 0,
along_track_interval = NULL, orientation = NULL, res = NULL,
method = "IDW", nmax = 20, idp = 2, model = "Sph", psill = NULL, range = NULL,
nugget = 0, kappa = NULL, trend_order = 1, zero_threshold = 0.05,
tolerance = 0.1, n_truth_examples = 3,
plot = TRUE, seed = NULL){
if(!is.null(seed)){
if(!is.numeric(seed))
stop("seed must be numeric")
set.seed(seed)
}
#transform outline shapefile into vector (accepts SpatVector, sf, or anything terra::vect() reads)
if(!inherits(outline, "SpatVector")){
if(inherits(outline, "sf") && requireNamespace("sf", quietly = TRUE)){
outline <- sf::st_zm(outline, drop = TRUE, what = "ZM")
}
outline <- terra::vect(outline)
}
else{
outline <- outline
}
#checks
if(!inherits(outline, "SpatVector"))
stop("outline is not a SpatVector or cannot be transformed")
if(missing(max_depth) || !is.numeric(max_depth) || max_depth <= 0)
stop("max_depth must be specified as a positive numeric value")
if(!method %in% c("IDW", "OK", "UK"))
stop("method misspecified. Please choose 'IDW', 'OK', or 'UK'")
if(!is.numeric(tolerance) || tolerance < 0)
stop("tolerance must be a non-negative numeric value")
#store the original CRS so the recommended transects can be returned in it
original_crs <- terra::crs(outline)
test_crs <- terra::crs(outline)
if(is.na(test_crs) || test_crs == ""){
stop("CRS of 'outline' is unable to be defined.")
}
# Function to determine the best UTM CRS for a given vector
get_best_utm <- function(outline) {
centroid <- terra::centroids(outline)
lon <- terra::crds(centroid)[1]
lat <- terra::crds(centroid)[2]
utm_zone <- base::floor((lon + 180) / 6) + 1
hemisphere <- base::ifelse(lat >= 0, 32600, 32700)
epsg_code <- hemisphere + utm_zone
return(terra::crs(paste0("EPSG:", epsg_code)))
}
#run the whole analysis in meters (UTM), same reasoning as interpBathy/crossValidate
if(terra::is.lonlat(outline)){
best_crs <- get_best_utm(outline)
outline <- terra::project(outline, best_crs)
}
message("Simulating in CRS: ", terra::crs(outline, describe = TRUE)$name)
#characteristic length scale of the waterbody, used to scale several defaults
area <- as.numeric(terra::expanse(outline))
L <- sqrt(area)
if(is.null(spacings)){
spacings <- round(L * c(.02, .05, .1, .2, .4))
spacings <- unique(spacings[spacings > 0])
}
if(!is.numeric(spacings) || any(spacings <= 0))
stop("spacings must be positive numeric values")
spacings <- sort(spacings)
if(is.null(res)){
res <- max(1, min(spacings) / 5)
}
if(is.null(along_track_interval)){
along_track_interval <- max(5, min(spacings) / 10)
}
if(is.null(truth_range)){
truth_range <- L / 4
}
if(is.null(truth_sill)){
truth_sill <- (max_depth / 5)^2
}
#### helper: build a masked raster template over the waterbody at a given resolution
build_grid <- function(outline, res){
ext_o <- terra::ext(outline)
ncol_r <- ceiling((ext_o$xmax - ext_o$xmin) / res)
nrow_r <- ceiling((ext_o$ymax - ext_o$ymin) / res)
empty <- terra::rast(ext_o, ncol = ncol_r, nrow = nrow_r, crs = terra::crs(outline))
ras <- terra::rasterize(outline, empty)
grid <- terra::mask(ras, outline)
return(grid)
}
#### helper: rotate xy points by theta (radians, counterclockwise) about a center point
rotate_pts <- function(xy, theta, center){
xy_c <- sweep(xy, 2, center, FUN = "-")
rot_mat <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), nrow = 2)
rotated <- xy_c %*% rot_mat
sweep(rotated, 2, center, FUN = "+")
}
#### helper: build parallel transect lines across outline at a given spacing/orientation,
#### clipped to the outline polygon
make_transects <- function(outline, spacing, orientation_deg = NULL){
verts <- terra::crds(outline)
center <- colMeans(verts)
if(is.null(orientation_deg)){
pc <- stats::prcomp(scale(verts, scale = FALSE))
theta <- atan2(pc$rotation[2, 1], pc$rotation[1, 1])
} else {
theta <- orientation_deg * pi / 180
}
rot_verts <- rotate_pts(verts, -theta, center)
xr <- range(rot_verts[, 1])
yr <- range(rot_verts[, 2])
pad <- spacing
y_lines <- seq(yr[1] - pad, yr[2] + pad, by = spacing)
line_list <- lapply(y_lines, function(yy){
pts_rot <- rbind(c(xr[1] - pad, yy), c(xr[2] + pad, yy))
pts_real <- rotate_pts(pts_rot, theta, center)
terra::vect(pts_real, type = "lines", crs = terra::crs(outline))
})
lines_all <- do.call(rbind, line_list)
lines_clip <- terra::intersect(lines_all, outline)
return(lines_clip)
}
#### helper: generate one synthetic "true" bathymetry raster
simulate_truth <- function(){
grid <- build_grid(outline, res)
dist_to_shore <- terra::distance(grid, terra::as.lines(outline))
dist_to_shore <- terra::mask(dist_to_shore, outline)
max_dist <- terra::global(dist_to_shore, "max", na.rm = TRUE)[[1]]
bowl <- max_depth * (dist_to_shore / max_dist)^shape
#simulate the random roughness component on a coarse grid sized off its own
#spatial autocorrelation range (truth_range), then resample up to match the
#fine bowl-shape grid. The roughness has no real structure finer than
#truth_range anyway, so this avoids running gstat's sequential Gaussian
#simulation over potentially hundreds of thousands to millions of cells on
#large waterbodies - the simulation cost grows much faster than linearly
#with cell count, so this is the main lever for keeping runtime practical.
coarse_res <- max(res, truth_range / 4)
coarse_grid <- build_grid(outline, coarse_res)
coarse_df <- as.data.frame(coarse_grid, xy = TRUE, na.rm = TRUE)
names(coarse_df)[1:2] <- c("x", "y")
g_dummy <- gstat::gstat(formula = z ~ 1, locations = ~x + y, dummy = TRUE, beta = 0,
model = gstat::vgm(psill = truth_sill, model = truth_model,
range = truth_range, nugget = truth_nugget),
nmax = 30)
sim <- stats::predict(g_dummy, newdata = coarse_df, nsim = 1)
roughness_coarse <- terra::rast(data.frame(x = coarse_df$x, y = coarse_df$y, roughness = sim$sim1),
type = "xyz", crs = terra::crs(outline))
roughness_fine <- terra::resample(roughness_coarse, bowl, method = "bilinear")
true_depth <- bowl + roughness_fine
true_depth <- terra::ifel(true_depth < 0, 0, true_depth)
true_depth <- terra::mask(true_depth, outline)
names(true_depth) <- "true_depth"
return(true_depth)
}
#### helper: simulate a survey at a given spacing against a given truth, reconstruct, and score
score_spacing <- function(truth_raster, spacing){
transects <- make_transects(outline, spacing, orientation)
densified <- terra::densify(transects, interval = along_track_interval)
sample_pts <- terra::as.points(densified)
vals <- terra::extract(truth_raster, sample_pts, ID = FALSE)
xy <- terra::crds(sample_pts)
sample_df <- data.frame(x = xy[, 1], y = xy[, 2], z = vals[[1]])
sample_df <- sample_df[!is.na(sample_df$z), ]
result <- tryCatch({
out <- utils::capture.output(
dem <- suppressMessages(interpBathy(outline, sample_df, x = "x", y = "y", z = "z",
zeros = FALSE, separation = along_track_interval,
res = res, method = method, nmax = nmax, idp = idp,
model = model, psill = psill, range = range,
nugget = nugget, kappa = kappa, trend_order = trend_order,
zero_threshold = zero_threshold))
)
if(method %in% c("OK", "UK")) dem <- dem[["depth"]]
dem_resampled <- terra::resample(dem, truth_raster)
sq_err <- (dem_resampled - truth_raster)^2
sqrt(terra::global(sq_err, "mean", na.rm = TRUE)[[1]])
}, error = function(e){
warning("Reconstruction failed for spacing = ", spacing, ": ", conditionMessage(e))
return(NA_real_)
})
return(result)
}
#### run the simulation: one truth surface per rep, tested against every spacing
n_iter <- n_sim * length(spacings)
pb <- utils::txtProgressBar(min = 0, max = n_iter, style = 3)
iter <- 0
rmse_matrix <- matrix(NA_real_, nrow = n_sim, ncol = length(spacings))
colnames(rmse_matrix) <- as.character(spacings)
n_truth_examples <- min(n_truth_examples, n_sim)
truth_list <- list()
for(i in 1:n_sim){
truth_raster <- simulate_truth()
if(i <= n_truth_examples){
truth_list[[i]] <- truth_raster
}
for(j in seq_along(spacings)){
rmse_matrix[i, j] <- score_spacing(truth_raster, spacings[j])
iter <- iter + 1
utils::setTxtProgressBar(pb, iter)
}
}
close(pb)
results <- data.frame(spacing = spacings,
mean_rmse = apply(rmse_matrix, 2, mean, na.rm = TRUE),
sd_rmse = apply(rmse_matrix, 2, stats::sd, na.rm = TRUE))
best_rmse <- min(results$mean_rmse, na.rm = TRUE)
recommended_spacing <- max(results$spacing[results$mean_rmse <= best_rmse * (1 + tolerance)])
#build the recommended transects once more, cleanly, and return them in the original CRS
transects_final <- make_transects(outline, recommended_spacing, orientation)
total_transect_length <- sum(terra::perim(transects_final))
transects_final <- terra::project(transects_final, original_crs)
outline_orig <- terra::project(outline, original_crs)
truth_examples <- do.call(c, truth_list)
names(truth_examples) <- paste0("sim_", seq_len(n_truth_examples))
truth_examples <- terra::project(truth_examples, original_crs)
if(plot){
graphics::plot(results$spacing, results$mean_rmse, type = "b", pch = 19,
xlab = "Transect spacing (m)", ylab = "RMSE (m)",
main = "Simulated DEM accuracy vs. transect spacing")
graphics::abline(v = recommended_spacing, lty = 2, col = "blue")
graphics::legend("topleft", legend = paste0("Recommended: ", recommended_spacing, " m"),
lty = 2, col = "blue", bty = "n")
terra::plot(outline_orig, main = paste0("Recommended survey design (", recommended_spacing, " m spacing)"))
terra::plot(transects_final, add = TRUE, col = "blue", lwd = 1.5)
terra::plot(truth_examples, col = grDevices::hcl.colors(100, "Blues"),
mar = c(2, 2, 2, 4))
}
return(list(results = results,
recommended_spacing = recommended_spacing,
transects = transects_final,
total_transect_length = total_transect_length,
truth_examples = truth_examples))
}
Try the rLakeHabitat package in your browser
Any scripts or data that you put into this service are public.
rLakeHabitat documentation built on July 30, 2026, 5:11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions samplingDensity samplingDensity samplingDensity
Documented in samplingDensity
#' Pre-Survey Power Analysis for Bathymetric Sampling Design
#'
#' Before any depth data are collected, estimate how far apart survey transects can be spaced and still produce an accurate bathymetric DEM. Because no real depth data exist yet, this works as a simulation study: a plausible "true" bathymetry is generated for the waterbody, a survey at each candidate transect spacing is simulated by sampling that truth, a DEM is reconstructed from those simulated samples (via interpBathy), and the reconstruction is compared back against the generated DEM. This is repeated many times per spacing to average out randomness, and across all candidate spacings, to build a curve of expected DEM accuracy vs. transect spacing - along with a recommended spacing and a map of what that survey design looks like on the actual waterbody outline.
#'
#' @param outline shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read.
#' @param max_depth expected maximum depth of the waterbody, in meters. This drives the synthetic bathymetry and has a large effect on results.
#' @param spacings numeric vector of candidate transect spacings to test, in meters. Default = NULL, in which case a sequence scaled to the waterbody's size is generated automatically (2%-40% of the lake's characteristic length).
#' @param n_sim number of independent "true" bathymetries to simulate and test each spacing against, default = 10. Higher values give a more stable accuracy estimate at the cost of runtime.
#' @param shape numeric value controlling the synthetic depth profile's shape: depth = max_depth * (relative distance to shore) ^ shape. shape = 1 (default) gives a linear, cone-like profile; shape < 1 gives a broader, flatter deep basin; shape > 1 gives a narrower, steeper-sided basin.
#' @param truth_model character variogram model used to generate the synthetic bathymetry's random roughness (see gstat::vgm options; "Sph", "Exp", "Gau", or "Mat"), default = "Exp".
#' @param truth_range numeric range parameter (in meters) for the synthetic roughness's spatial autocorrelation, default = NULL, in which case it is set to one quarter of the waterbody's characteristic length (sqrt(area)).
#' @param truth_sill numeric sill (variance, in squared meters) for the synthetic roughness, default = NULL, in which case it is set to (max_depth/5)^2.
#' @param truth_nugget numeric nugget for the synthetic roughness variogram, default = 0.
#' @param along_track_interval numeric spacing (in meters) at which simulated samples are drawn along each transect line, representing how frequently a depth is recorded. Default = NULL, in which case it is set to 1/10th of the smallest tested spacing.
#' @param orientation numeric angle, in degrees (0 = along the x-axis/east, increasing counterclockwise), giving the orientation of transect lines. Default = NULL, in which case the orientation is chosen automatically to align with the waterbody's long axis.
#' @param res numeric DEM cell resolution (in meters) used both for generating the synthetic truth surface and for reconstructing the DEM from simulated samples. Default = NULL, in which case it is set to 1/5th of the smallest tested spacing. This should generally be left fine relative to 'spacings'.
#' @param method character describing method of interpolation used to reconstruct DEMs, "IDW" or "OK". Default = "IDW" (recommended for this use, since it is run many times and OK's variogram fitting adds considerable runtime).
#' @param nmax numeric value describing number of neighbors used in interpolation, default = 20
#' @param idp numeric value describing inverse distance power value for IDW interpolation
#' @param model character describing type of model used in Ordinary Kriging, options include 'Sph', 'Exp', 'Gau', 'Mat', default = 'Sph'
#' @param psill numeric value describing the partial sill value for OK interpolation, default = NULL
#' @param range numeric describing distance beyond which there is no spatial correlation in Ordinary Kriging models, default = NULL
#' @param nugget numeric describing variance at zero distance in Ordinary Kriging models, default = 0
#' @param kappa numeric value describing model smoothness, default = NULL
#' @param trend_order numeric value (1 or 2) giving the order of the polynomial trend surface fit for Universal Kriging. 1 = linear trend (z ~ x + y), 2 = quadratic trend. Default = 1.
#' @param zero_threshold numeric proportion (0-1) of the waterbody's surface area that must interpolate to exactly 0 before the automatic zero re-interpolation pass runs, default = 0.05 (5%). A handful of scattered zero cells won't trigger it; a large contiguous block collapsing to zero (typically an interpolation artifact, often from the shoreline zero ring dominating a narrow bay or inlet) will.
#' @param tolerance numeric value (proportion) used to pick the recommended spacing: the coarsest tested spacing whose mean RMSE is still within 'tolerance' of the best (finest-spacing) RMSE. Default = 0.1 (10%).
#' @param n_truth_examples numeric value describing how many of the n_sim synthetic "true" bathymetries to save and return/plot, so you can visually check whether the assumed depth profile and roughness look like your lake before trusting the rest of the results. Default = 3 (capped at n_sim).
#' @param plot logical: should diagnostic plots (accuracy-vs-spacing curve, and outline + recommended transects map) be drawn? Default = TRUE.
#' @param seed optional numeric value used to seed the random number generator, for reproducible results across runs. Default = NULL (not seeded).
#' @details
#' This is a power analysis under assumed, not observed, bathymetry - it can only be as realistic as the 'max_depth', 'shape', 'truth_range', and 'truth_sill' inputs. Once real depth data have been collected, use crossValidate(), interpBathy(), and optimizeParams() on the actual data to check whether observed accuracy matches what was predicted here, and adjust future survey effort accordingly.
#' Runtime scales with n_sim x length(spacings) x (cost of one synthetic-surface simulation + one interpBathy call). For a first look, consider a smaller n_sim (e.g. 5) and a coarser 'res' before committing to a longer run.
#' @return a list with:
#' \describe{
#' \item{results}{a data frame of spacing, mean RMSE, and SD of RMSE across simulations}
#' \item{recommended_spacing}{the recommended transect spacing, in meters}
#' \item{transects}{a SpatVector of the recommended transect lines, in the original CRS of 'outline'}
#' \item{total_transect_length}{total length of the recommended transects, in meters}
#' \item{truth_examples}{a multi-layer SpatRaster (one layer per saved example) of synthetic "true" bathymetries used
#' during the simulation, in the original CRS of 'outline', for visually sanity-checking the assumed depth profile}
#' }
#' @author Tristan Blechinger, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @import dplyr
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @import gstat
#' @examples
#' \donttest{
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' samplingDensity(outline, max_depth = 40, n_sim = 5, seed = 123)}
samplingDensity <- function(outline, max_depth, spacings = NULL, n_sim = 10, shape = 1,
truth_model = "Exp", truth_range = NULL, truth_sill = NULL, truth_nugget = 0,
along_track_interval = NULL, orientation = NULL, res = NULL,
method = "IDW", nmax = 20, idp = 2, model = "Sph", psill = NULL, range = NULL,
nugget = 0, kappa = NULL, trend_order = 1, zero_threshold = 0.5, tolerance = 0.1, n_truth_examples = 3,
plot = TRUE, seed = NULL){
if(!is.null(seed)){
if(!is.numeric(seed))
stop("seed must be numeric")
set.seed(seed)
}
# transform outline shapefile into vector
if(!inherits(outline, "SpatVector")){
if(inherits(outline, "sf") && requireNamespace("sf", quietly = TRUE)){
outline <- sf::st_zm(outline, drop = TRUE, what = "ZM")
}
outline <- terra::vect(outline)
}
else{
outline <- outline
}
# checks
if(!inherits(outline, "SpatVector"))
stop("outline is not a SpatVector or cannot be transformed")
if(missing(max_depth) || !is.numeric(max_depth) || max_depth <= 0)
stop("max_depth must be specified as a positive numeric value")
if(!method %in% c("IDW", "OK"))
stop("method misspecified. Please choose either 'IDW' or 'OK'")
if(!is.numeric(tolerance) || tolerance < 0)
stop("tolerance must be a non-negative numeric value")
#store the original CRS so the recommended transects can be returned in it
original_crs <- terra::crs(outline)
test_crs <- terra::crs(outline)
if(is.na(test_crs) || test_crs == ""){
stop("CRS of 'outline' is unable to be defined.")
}
# Function to determine the best UTM CRS for a given vector
get_best_utm <- function(outline) {
centroid <- terra::centroids(outline)
lon <- terra::crds(centroid)[1]
lat <- terra::crds(centroid)[2]
utm_zone <- base::floor((lon + 180) / 6) + 1
hemisphere <- base::ifelse(lat >= 0, 32600, 32700)
epsg_code <- hemisphere + utm_zone
return(terra::crs(paste0("EPSG:", epsg_code)))
}
# run analysis in meters (UTM)
if(terra::is.lonlat(outline)){
best_crs <- get_best_utm(outline)
outline <- terra::project(outline, best_crs)
}
message("Simulating in CRS: ", terra::crs(outline, describe = TRUE)$name)
# characteristic length scale of the waterbody, used to scale several defaults
area <- as.numeric(terra::expanse(outline))
L <- sqrt(area)
if(is.null(spacings)){
spacings <- round(L * c(.02, .05, .1, .2, .4))
spacings <- unique(spacings[spacings > 0])
}
if(!is.numeric(spacings) || any(spacings <= 0))
stop("spacings must be positive numeric values")
spacings <- sort(spacings)
if(is.null(res)){
res <- max(1, min(spacings) / 5)
}
if(is.null(along_track_interval)){
along_track_interval <- max(5, min(spacings) / 10)
}
if(is.null(truth_range)){
truth_range <- L / 4
}
if(is.null(truth_sill)){
truth_sill <- (max_depth / 5)^2
}
#### helper: build a masked raster template over the waterbody at a given resolution
build_grid <- function(outline, res){
ext_o <- terra::ext(outline)
ncol_r <- ceiling((ext_o$xmax - ext_o$xmin) / res)
nrow_r <- ceiling((ext_o$ymax - ext_o$ymin) / res)
empty <- terra::rast(ext_o, ncol = ncol_r, nrow = nrow_r, crs = terra::crs(outline))
ras <- terra::rasterize(outline, empty)
grid <- terra::mask(ras, outline)
return(grid)
}
#### helper: rotate xy points by theta (radians, counterclockwise) about a center point
rotate_pts <- function(xy, theta, center){
xy_c <- sweep(xy, 2, center, FUN = "-")
rot_mat <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), nrow = 2)
rotated <- xy_c %*% rot_mat
sweep(rotated, 2, center, FUN = "+")
}
#### helper: build parallel transect lines across outline at a given spacing/orientation,
#### clipped to the outline polygon
make_transects <- function(outline, spacing, orientation_deg = NULL){
verts <- terra::crds(outline)
center <- colMeans(verts)
if(is.null(orientation_deg)){
pc <- stats::prcomp(scale(verts, scale = FALSE))
theta <- atan2(pc$rotation[2, 1], pc$rotation[1, 1])
} else {
theta <- orientation_deg * pi / 180
}
rot_verts <- rotate_pts(verts, -theta, center)
xr <- range(rot_verts[, 1])
yr <- range(rot_verts[, 2])
pad <- spacing
y_lines <- seq(yr[1] - pad, yr[2] + pad, by = spacing)
line_list <- lapply(y_lines, function(yy){
pts_rot <- rbind(c(xr[1] - pad, yy), c(xr[2] + pad, yy))
pts_real <- rotate_pts(pts_rot, theta, center)
terra::vect(pts_real, type = "lines", crs = terra::crs(outline))
})
lines_all <- do.call(rbind, line_list)
lines_clip <- terra::intersect(lines_all, outline)
return(lines_clip)
}
#### helper: generate one synthetic "true" bathymetry raster
simulate_truth <- function(){
grid <- build_grid(outline, res)
dist_to_shore <- terra::distance(grid, terra::as.lines(outline))
dist_to_shore <- terra::mask(dist_to_shore, outline)
max_dist <- terra::global(dist_to_shore, "max", na.rm = TRUE)[[1]]
bowl <- max_depth * (dist_to_shore / max_dist)^shape
# simulate the random roughness component on a coarse grid sized off its own spatial autocorrelation range, then resample to match the fine grid
coarse_res <- max(res, truth_range / 4)
coarse_grid <- build_grid(outline, coarse_res)
coarse_df <- as.data.frame(coarse_grid, xy = TRUE, na.rm = TRUE)
names(coarse_df)[1:2] <- c("x", "y")
g_dummy <- gstat::gstat(formula = z ~ 1, locations = ~x + y, dummy = TRUE, beta = 0,
model = gstat::vgm(psill = truth_sill, model = truth_model,
range = truth_range, nugget = truth_nugget),
nmax = 30)
sim <- stats::predict(g_dummy, newdata = coarse_df, nsim = 1)
roughness_coarse <- terra::rast(data.frame(x = coarse_df$x, y = coarse_df$y, roughness = sim$sim1),
type = "xyz", crs = terra::crs(outline))
roughness_fine <- terra::resample(roughness_coarse, bowl, method = "bilinear")
true_depth <- bowl + roughness_fine
true_depth <- terra::ifel(true_depth < 0, 0, true_depth)
true_depth <- terra::mask(true_depth, outline)
names(true_depth) <- "true_depth"
return(true_depth)
}
#### helper: simulate a survey at a given spacing against a given truth, reconstruct, and score
score_spacing <- function(truth_raster, spacing){
transects <- make_transects(outline, spacing, orientation)
densified <- terra::densify(transects, interval = along_track_interval)
sample_pts <- terra::as.points(densified)
vals <- terra::extract(truth_raster, sample_pts, ID = FALSE)
xy <- terra::crds(sample_pts)
sample_df <- data.frame(x = xy[, 1], y = xy[, 2], z = vals[[1]])
sample_df <- sample_df[!is.na(sample_df$z), ]
result <- tryCatch({
out <- utils::capture.output(
dem <- suppressMessages(interpBathy(outline, sample_df, x = "x", y = "y", z = "z",
zeros = FALSE, separation = along_track_interval,
res = res, method = method, nmax = nmax, idp = idp,
model = model, psill = psill, range = range,
nugget = nugget, kappa = kappa, trend_order = trend_order,
zero_threshold = zero_threshold))
)
if(method == "OK") dem <- dem[["depth"]]
dem_resampled <- terra::resample(dem, truth_raster)
sq_err <- (dem_resampled - truth_raster)^2
sqrt(terra::global(sq_err, "mean", na.rm = TRUE)[[1]])
}, error = function(e){
warning("Reconstruction failed for spacing = ", spacing, ": ", conditionMessage(e))
return(NA_real_)
})
return(result)
}
# Run simulation
n_iter <- n_sim * length(spacings)
pb <- utils::txtProgressBar(min = 0, max = n_iter, style = 3)
iter <- 0
rmse_matrix <- matrix(NA_real_, nrow = n_sim, ncol = length(spacings))
colnames(rmse_matrix) <- as.character(spacings)
n_truth_examples <- min(n_truth_examples, n_sim)
truth_list <- list()
for(i in 1:n_sim){
truth_raster <- simulate_truth()
if(i <= n_truth_examples){
truth_list[[i]] <- truth_raster
}
for(j in seq_along(spacings)){
rmse_matrix[i, j] <- score_spacing(truth_raster, spacings[j])
iter <- iter + 1
utils::setTxtProgressBar(pb, iter)
}
}
close(pb)
results <- data.frame(spacing = spacings,
mean_rmse = apply(rmse_matrix, 2, mean, na.rm = TRUE),
sd_rmse = apply(rmse_matrix, 2, stats::sd, na.rm = TRUE))
best_rmse <- min(results$mean_rmse, na.rm = TRUE)
recommended_spacing <- max(results$spacing[results$mean_rmse <= best_rmse * (1 + tolerance)])
#build the recommended transects once more, cleanly, and return them in the original CRS
transects_final <- make_transects(outline, recommended_spacing, orientation)
total_transect_length <- sum(terra::perim(transects_final))
transects_final <- terra::project(transects_final, original_crs)
outline_orig <- terra::project(outline, original_crs)
truth_examples <- do.call(c, truth_list)
names(truth_examples) <- paste0("sim_", seq_len(n_truth_examples))
truth_examples <- terra::project(truth_examples, original_crs)
if(plot){
graphics::plot(results$spacing, results$mean_rmse, type = "b", pch = 19,
xlab = "Transect spacing (m)", ylab = "RMSE (m)",
main = "Simulated DEM accuracy vs. transect spacing")
graphics::abline(v = recommended_spacing, lty = 2, col = "blue")
graphics::legend("topleft", legend = paste0("Recommended: ", recommended_spacing, " m"),
lty = 2, col = "blue", bty = "n")
terra::plot(outline_orig, main = paste0("Recommended survey design (", recommended_spacing, " m spacing)"))
terra::plot(transects_final, add = TRUE, col = "blue", lwd = 1.5)
terra::plot(truth_examples, col = grDevices::hcl.colors(100, "Blues"),
mar = c(2, 2, 2, 4))
}
return(list(results = results,
recommended_spacing = recommended_spacing,
transects = transects_final,
total_transect_length = total_transect_length,
truth_examples = truth_examples))
}
#' Pre-Survey Power Analysis for Bathymetric Sampling Design
#'
#' Before any depth data are collected, estimate how far apart survey transects can be spaced and still produce an accurate bathymetric DEM. Because no real depth data exist yet, this works as a simulation study: a plausible "true" bathymetry is generated for the waterbody, a survey at each candidate transect spacing is simulated by sampling that truth, a DEM is reconstructed from those simulated samples (via interpBathy), and the reconstruction is compared back against the generated DEM. This is repeated many times per spacing to average out randomness, and across all candidate spacings, to build a curve of expected DEM accuracy vs. transect spacing - along with a recommended spacing and a map of what that survey design looks like on the actual waterbody outline.
#'
#' @param outline shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read.
#' @param max_depth expected maximum depth of the waterbody, in meters. This drives the synthetic bathymetry and has a large effect on the results.
#' @param spacings numeric vector of candidate transect spacings to test, in meters. Default = NULL, in which case a sequence scaled to the waterbody's size is generated automatically (roughly 2%-40% of the lake's characteristic length).
#' @param n_sim number of independent synthetic "true" bathymetries to simulate and test each spacing against, default = 10.
#' @param shape numeric value controlling the synthetic depth profile's shape: depth = max_depth * (relative distance to shore) ^ shape. shape = 1 (default) gives a linear, cone-like profile; shape > 1 gives a broader, flatter deep basin; shape < 1 gives a narrower, steeper-sided basin.
#' @param truth_model character variogram model used to generate the synthetic bathymetry's random roughness (see gstat::vgm options), default = "Exp".
#' @param truth_range numeric range parameter (in meters) for the synthetic roughness's spatial autocorrelation, default = NULL, in which case it is set to one quarter of the waterbody's characteristic length (sqrt(area)).
#' @param truth_sill numeric sill (variance, in squared meters) for the synthetic roughness, default = NULL, in which case it is set to (max_depth/5)^2.
#' @param truth_nugget numeric nugget for the synthetic roughness variogram, default = 0.
#' @param along_track_interval numeric spacing (in meters) at which simulated samples are drawn along each transect line, representing how frequently a sonar unit records depth. Default = NULL, in which case it is set to 1/10th of the smallest tested spacing.
#' @param orientation numeric angle, in degrees (0 = along the x-axis/east, increasing counterclockwise), giving the orientation of transect lines. Default = NULL, in which case the orientation is chosen automatically to align with the waterbody's long axis.
#' @param res numeric DEM cell resolution (in meters) used both for generating the synthetic truth surface and for reconstructing the DEM from simulated samples. Default = NULL, in which case it is set to 1/5th of the smallest tested spacing.
#' @param method character describing method of interpolation used to reconstruct DEMs, "IDW", "OK", or "UK". Default = "IDW" (recommended for this use, since it is run many times and OK's variogram fitting adds considerable runtime).
#' @param nmax numeric value describing number of neighbors used in interpolation, default = 20
#' @param idp numeric value describing inverse distance power value for IDW interpolation
#' @param model character describing type of model used in Ordinary Kriging, options include 'Sph', 'Exp', 'Gau', 'Sta', default = 'Sph'
#' @param psill numeric value describing the partial sill value for OK interpolation, default = NULL
#' @param range numeric describing distance beyond which there is no spatial correlation in Ordinary Kriging models, default = NULL
#' @param nugget numeric describing variance at zero distance in Ordinary Kriging models, default = 0
#' @param kappa numeric value describing model smoothness, default = NULL
#' @param trend_order numeric value (1 or 2) giving the polynomial trend order for Universal Kriging ("UK" only), default = 1
#' @param zero_threshold numeric proportion (0-1) of surface area that must interpolate to exactly 0 before the automatic zero re-interpolation pass runs - passed through to interpBathy(). Default = 0.05.
#' @param tolerance numeric value (proportion) used to pick the recommended spacing: the coarsest tested spacing whose mean RMSE is still within 'tolerance' of the best (finest-spacing) RMSE. Default = 0.1 (10%).
#' @param n_truth_examples numeric value giving how many of the n_sim synthetic "true" bathymetries to save and return/plot, so you can visually check whether the assumed depth profile and roughness look like your lake before trusting the rest of the results. Default = 3 (capped at n_sim).
#' @param plot logical: should diagnostic plots (accuracy-vs-spacing curve, and outline + recommended transects map) be drawn? Default = TRUE.
#' @param seed optional numeric value used to seed the random number generator, for reproducible results across runs. Default = NULL (not seeded).
#' @details
#' This is a power analysis under assumed, not observed, bathymetry - it can only be as realistic as the 'max_depth',
#' 'shape', 'truth_range', and 'truth_sill' inputs. Once real depth data have been collected, use crossValidate() and
#' interpBathy() on the actual data to check whether observed accuracy matches what was predicted here, and adjust
#' future survey effort accordingly.
#' @return a list with:
#' \describe{
#' \item{results}{a data frame of spacing, mean RMSE, and SD of RMSE across simulations}
#' \item{recommended_spacing}{the recommended transect spacing, in meters}
#' \item{transects}{a SpatVector of the recommended transect lines, in the original CRS of 'outline'}
#' \item{total_transect_length}{total length of the recommended transects, in meters}
#' \item{truth_examples}{a multi-layer SpatRaster (one layer per saved example) of synthetic "true" bathymetries used
#' during the simulation, in the original CRS of 'outline', for visually sanity-checking the assumed depth profile}
#' }
#' @author Tristan Blechinger, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @import dplyr
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @import gstat
#' @examples
#' \donttest{
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' samplingDensity(outline, max_depth = 40, n_sim = 5, seed = 123)}
samplingDensity <- function(outline, max_depth, spacings = NULL, n_sim = 10, shape = 1,
truth_model = "Exp", truth_range = NULL, truth_sill = NULL, truth_nugget = 0,
along_track_interval = NULL, orientation = NULL, res = NULL,
method = "IDW", nmax = 20, idp = 2, model = "Sph", psill = NULL, range = NULL,
nugget = 0, kappa = NULL, trend_order = 1, zero_threshold = 0.05,
tolerance = 0.1, n_truth_examples = 3,
plot = TRUE, seed = NULL){
if(!is.null(seed)){
if(!is.numeric(seed))
stop("seed must be numeric")
set.seed(seed)
}
#transform outline shapefile into vector (accepts SpatVector, sf, or anything terra::vect() reads)
if(!inherits(outline, "SpatVector")){
if(inherits(outline, "sf") && requireNamespace("sf", quietly = TRUE)){
outline <- sf::st_zm(outline, drop = TRUE, what = "ZM")
}
outline <- terra::vect(outline)
}
else{
outline <- outline
}
#checks
if(!inherits(outline, "SpatVector"))
stop("outline is not a SpatVector or cannot be transformed")
if(missing(max_depth) || !is.numeric(max_depth) || max_depth <= 0)
stop("max_depth must be specified as a positive numeric value")
if(!method %in% c("IDW", "OK", "UK"))
stop("method misspecified. Please choose 'IDW', 'OK', or 'UK'")
if(!is.numeric(tolerance) || tolerance < 0)
stop("tolerance must be a non-negative numeric value")
#store the original CRS so the recommended transects can be returned in it
original_crs <- terra::crs(outline)
test_crs <- terra::crs(outline)
if(is.na(test_crs) || test_crs == ""){
stop("CRS of 'outline' is unable to be defined.")
}
# Function to determine the best UTM CRS for a given vector
get_best_utm <- function(outline) {
centroid <- terra::centroids(outline)
lon <- terra::crds(centroid)[1]
lat <- terra::crds(centroid)[2]
utm_zone <- base::floor((lon + 180) / 6) + 1
hemisphere <- base::ifelse(lat >= 0, 32600, 32700)
epsg_code <- hemisphere + utm_zone
return(terra::crs(paste0("EPSG:", epsg_code)))
}
#run the whole analysis in meters (UTM), same reasoning as interpBathy/crossValidate
if(terra::is.lonlat(outline)){
best_crs <- get_best_utm(outline)
outline <- terra::project(outline, best_crs)
}
message("Simulating in CRS: ", terra::crs(outline, describe = TRUE)$name)
#characteristic length scale of the waterbody, used to scale several defaults
area <- as.numeric(terra::expanse(outline))
L <- sqrt(area)
if(is.null(spacings)){
spacings <- round(L * c(.02, .05, .1, .2, .4))
spacings <- unique(spacings[spacings > 0])
}
if(!is.numeric(spacings) || any(spacings <= 0))
stop("spacings must be positive numeric values")
spacings <- sort(spacings)
if(is.null(res)){
res <- max(1, min(spacings) / 5)
}
if(is.null(along_track_interval)){
along_track_interval <- max(5, min(spacings) / 10)
}
if(is.null(truth_range)){
truth_range <- L / 4
}
if(is.null(truth_sill)){
truth_sill <- (max_depth / 5)^2
}
#### helper: build a masked raster template over the waterbody at a given resolution
build_grid <- function(outline, res){
ext_o <- terra::ext(outline)
ncol_r <- ceiling((ext_o$xmax - ext_o$xmin) / res)
nrow_r <- ceiling((ext_o$ymax - ext_o$ymin) / res)
empty <- terra::rast(ext_o, ncol = ncol_r, nrow = nrow_r, crs = terra::crs(outline))
ras <- terra::rasterize(outline, empty)
grid <- terra::mask(ras, outline)
return(grid)
}
#### helper: rotate xy points by theta (radians, counterclockwise) about a center point
rotate_pts <- function(xy, theta, center){
xy_c <- sweep(xy, 2, center, FUN = "-")
rot_mat <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), nrow = 2)
rotated <- xy_c %*% rot_mat
sweep(rotated, 2, center, FUN = "+")
}
#### helper: build parallel transect lines across outline at a given spacing/orientation,
#### clipped to the outline polygon
make_transects <- function(outline, spacing, orientation_deg = NULL){
verts <- terra::crds(outline)
center <- colMeans(verts)
if(is.null(orientation_deg)){
pc <- stats::prcomp(scale(verts, scale = FALSE))
theta <- atan2(pc$rotation[2, 1], pc$rotation[1, 1])
} else {
theta <- orientation_deg * pi / 180
}
rot_verts <- rotate_pts(verts, -theta, center)
xr <- range(rot_verts[, 1])
yr <- range(rot_verts[, 2])
pad <- spacing
y_lines <- seq(yr[1] - pad, yr[2] + pad, by = spacing)
line_list <- lapply(y_lines, function(yy){
pts_rot <- rbind(c(xr[1] - pad, yy), c(xr[2] + pad, yy))
pts_real <- rotate_pts(pts_rot, theta, center)
terra::vect(pts_real, type = "lines", crs = terra::crs(outline))
})
lines_all <- do.call(rbind, line_list)
lines_clip <- terra::intersect(lines_all, outline)
return(lines_clip)
}
#### helper: generate one synthetic "true" bathymetry raster
simulate_truth <- function(){
grid <- build_grid(outline, res)
dist_to_shore <- terra::distance(grid, terra::as.lines(outline))
dist_to_shore <- terra::mask(dist_to_shore, outline)
max_dist <- terra::global(dist_to_shore, "max", na.rm = TRUE)[[1]]
bowl <- max_depth * (dist_to_shore / max_dist)^shape
#simulate the random roughness component on a coarse grid sized off its own
#spatial autocorrelation range (truth_range), then resample up to match the
#fine bowl-shape grid. The roughness has no real structure finer than
#truth_range anyway, so this avoids running gstat's sequential Gaussian
#simulation over potentially hundreds of thousands to millions of cells on
#large waterbodies - the simulation cost grows much faster than linearly
#with cell count, so this is the main lever for keeping runtime practical.
coarse_res <- max(res, truth_range / 4)
coarse_grid <- build_grid(outline, coarse_res)
coarse_df <- as.data.frame(coarse_grid, xy = TRUE, na.rm = TRUE)
names(coarse_df)[1:2] <- c("x", "y")
g_dummy <- gstat::gstat(formula = z ~ 1, locations = ~x + y, dummy = TRUE, beta = 0,
model = gstat::vgm(psill = truth_sill, model = truth_model,
range = truth_range, nugget = truth_nugget),
nmax = 30)
sim <- stats::predict(g_dummy, newdata = coarse_df, nsim = 1)
roughness_coarse <- terra::rast(data.frame(x = coarse_df$x, y = coarse_df$y, roughness = sim$sim1),
type = "xyz", crs = terra::crs(outline))
roughness_fine <- terra::resample(roughness_coarse, bowl, method = "bilinear")
true_depth <- bowl + roughness_fine
true_depth <- terra::ifel(true_depth < 0, 0, true_depth)
true_depth <- terra::mask(true_depth, outline)
names(true_depth) <- "true_depth"
return(true_depth)
}
#### helper: simulate a survey at a given spacing against a given truth, reconstruct, and score
score_spacing <- function(truth_raster, spacing){
transects <- make_transects(outline, spacing, orientation)
densified <- terra::densify(transects, interval = along_track_interval)
sample_pts <- terra::as.points(densified)
vals <- terra::extract(truth_raster, sample_pts, ID = FALSE)
xy <- terra::crds(sample_pts)
sample_df <- data.frame(x = xy[, 1], y = xy[, 2], z = vals[[1]])
sample_df <- sample_df[!is.na(sample_df$z), ]
result <- tryCatch({
out <- utils::capture.output(
dem <- suppressMessages(interpBathy(outline, sample_df, x = "x", y = "y", z = "z",
zeros = FALSE, separation = along_track_interval,
res = res, method = method, nmax = nmax, idp = idp,
model = model, psill = psill, range = range,
nugget = nugget, kappa = kappa, trend_order = trend_order,
zero_threshold = zero_threshold))
)
if(method %in% c("OK", "UK")) dem <- dem[["depth"]]
dem_resampled <- terra::resample(dem, truth_raster)
sq_err <- (dem_resampled - truth_raster)^2
sqrt(terra::global(sq_err, "mean", na.rm = TRUE)[[1]])
}, error = function(e){
warning("Reconstruction failed for spacing = ", spacing, ": ", conditionMessage(e))
return(NA_real_)
})
return(result)
}
#### run the simulation: one truth surface per rep, tested against every spacing
n_iter <- n_sim * length(spacings)
pb <- utils::txtProgressBar(min = 0, max = n_iter, style = 3)
iter <- 0
rmse_matrix <- matrix(NA_real_, nrow = n_sim, ncol = length(spacings))
colnames(rmse_matrix) <- as.character(spacings)
n_truth_examples <- min(n_truth_examples, n_sim)
truth_list <- list()
for(i in 1:n_sim){
truth_raster <- simulate_truth()
if(i <= n_truth_examples){
truth_list[[i]] <- truth_raster
}
for(j in seq_along(spacings)){
rmse_matrix[i, j] <- score_spacing(truth_raster, spacings[j])
iter <- iter + 1
utils::setTxtProgressBar(pb, iter)
}
}
close(pb)
results <- data.frame(spacing = spacings,
mean_rmse = apply(rmse_matrix, 2, mean, na.rm = TRUE),
sd_rmse = apply(rmse_matrix, 2, stats::sd, na.rm = TRUE))
best_rmse <- min(results$mean_rmse, na.rm = TRUE)
recommended_spacing <- max(results$spacing[results$mean_rmse <= best_rmse * (1 + tolerance)])
#build the recommended transects once more, cleanly, and return them in the original CRS
transects_final <- make_transects(outline, recommended_spacing, orientation)
total_transect_length <- sum(terra::perim(transects_final))
transects_final <- terra::project(transects_final, original_crs)
outline_orig <- terra::project(outline, original_crs)
truth_examples <- do.call(c, truth_list)
names(truth_examples) <- paste0("sim_", seq_len(n_truth_examples))
truth_examples <- terra::project(truth_examples, original_crs)
if(plot){
graphics::plot(results$spacing, results$mean_rmse, type = "b", pch = 19,
xlab = "Transect spacing (m)", ylab = "RMSE (m)",
main = "Simulated DEM accuracy vs. transect spacing")
graphics::abline(v = recommended_spacing, lty = 2, col = "blue")
graphics::legend("topleft", legend = paste0("Recommended: ", recommended_spacing, " m"),
lty = 2, col = "blue", bty = "n")
terra::plot(outline_orig, main = paste0("Recommended survey design (", recommended_spacing, " m spacing)"))
terra::plot(transects_final, add = TRUE, col = "blue", lwd = 1.5)
terra::plot(truth_examples, col = grDevices::hcl.colors(100, "Blues"),
mar = c(2, 2, 2, 4))
}
return(list(results = results,
recommended_spacing = recommended_spacing,
transects = transects_final,
total_transect_length = total_transect_length,
truth_examples = truth_examples))
}
#' Pre-Survey Power Analysis for Bathymetric Sampling Design
#'
#' Before any depth data are collected, estimate how far apart survey transects can be spaced
#' and still produce an accurate bathymetric DEM. Because no real depth data exist yet, this
#' works as a simulation study: a plausible synthetic "true" bathymetry is generated for the
#' waterbody (a bowl-shaped depth profile plus spatially-correlated random roughness), a survey
#' at each candidate transect spacing is simulated by sampling that synthetic truth, a DEM is
#' reconstructed from those simulated samples (via interpBathy), and the reconstruction is
#' compared back against the known synthetic truth. This is repeated many times per spacing to
#' average out randomness, and across all candidate spacings, to build a curve of expected DEM
#' accuracy vs. transect spacing - along with a recommended spacing and a map of what that
#' survey design looks like on the actual lake outline.
#'
#' @param outline shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read.
#' @param max_depth expected maximum depth of the waterbody, in meters. This drives the synthetic bathymetry and has a large
#' effect on the results - use the best estimate available (a known max depth, an old chart, or a knowledgeable guess).
#' @param spacings numeric vector of candidate transect spacings to test, in meters. Default = NULL, in which case a
#' sequence scaled to the waterbody's size is generated automatically (roughly 2%%-40%% of the lake's characteristic length).
#' @param n_sim number of independent synthetic "true" bathymetries to simulate and test each spacing against, default = 10.
#' Higher values give a more stable accuracy estimate at the cost of runtime.
#' @param shape numeric value controlling the synthetic depth profile's shape: depth = max_depth * (relative distance to
#' shore) ^ shape. shape = 1 (default) gives a linear, cone-like profile; shape < 1 gives a broader, flatter deep basin;
#' shape > 1 gives a narrower, steeper-sided basin.
#' @param truth_model character variogram model used to generate the synthetic bathymetry's random roughness (see
#' gstat::vgm options), default = "Exp".
#' @param truth_range numeric range parameter (in meters) for the synthetic roughness's spatial autocorrelation, default =
#' NULL, in which case it is set to one quarter of the waterbody's characteristic length (sqrt(area)).
#' @param truth_sill numeric sill (variance, in squared meters) for the synthetic roughness, default = NULL, in which case
#' it is set to (max_depth/5)^2.
#' @param truth_nugget numeric nugget for the synthetic roughness variogram, default = 0.
#' @param along_track_interval numeric spacing (in meters) at which simulated samples are drawn along each transect line,
#' representing how frequently a sonar unit records depth while underway. Default = NULL, in which case it is set to
#' 1/10th of the smallest tested spacing.
#' @param orientation numeric angle, in degrees (0 = along the x-axis/east, increasing counterclockwise), giving the
#' orientation of transect lines. Default = NULL, in which case the orientation is chosen automatically to align with the
#' waterbody's long axis.
#' @param res numeric DEM cell resolution (in meters) used both for generating the synthetic truth surface and for
#' reconstructing the DEM from simulated samples. Default = NULL, in which case it is set to 1/5th of the smallest tested
#' spacing. This should generally be left fine relative to 'spacings'.
#' @param method character describing method of interpolation used to reconstruct DEMs, "IDW", "OK", or "UK". Default = "IDW"
#' (recommended for this use, since it is run many times and OK's variogram fitting adds considerable runtime).
#' @param nmax numeric value describing number of neighbors used in interpolation, default = 20
#' @param idp numeric value describing inverse distance power value for IDW interpolation
#' @param model character describing type of model used in Ordinary Kriging, options include 'Sph', 'Exp', 'Gau', 'Sta', default = 'Sph'
#' @param psill numeric value describing the partial sill value for OK interpolation, default = NULL
#' @param range numeric describing distance beyond which there is no spatial correlation in Ordinary Kriging models, default = NULL
#' @param nugget numeric describing variance at zero distance in Ordinary Kriging models, default = 0
#' @param kappa numeric value describing model smoothness, default = NULL
#' @param trend_order numeric value (1 or 2) giving the polynomial trend order for Universal Kriging ("UK" only), default = 1
#' @param zero_threshold numeric proportion (0-1) of surface area that must interpolate to exactly 0 before the
#' automatic zero re-interpolation pass runs - passed through to interpBathy(). Default = 0.05.
#' @param tolerance numeric value (proportion) used to pick the recommended spacing: the coarsest tested spacing whose mean
#' RMSE is still within 'tolerance' of the best (finest-spacing) RMSE. Default = 0.1 (10%%).
#' @param n_truth_examples numeric value giving how many of the n_sim synthetic "true" bathymetries to save and return/plot,
#' so you can visually check whether the assumed depth profile and roughness look like your lake before trusting the
#' rest of the results. Default = 3 (capped at n_sim).
#' @param plot logical: should diagnostic plots (accuracy-vs-spacing curve, and outline + recommended transects map) be
#' drawn? Default = TRUE.
#' @param seed optional numeric value used to seed the random number generator, for reproducible results across runs.
#' Default = NULL (not seeded).
#' @details
#' This is a power analysis under assumed, not observed, bathymetry - it can only be as realistic as the 'max_depth',
#' 'shape', 'truth_range', and 'truth_sill' inputs. Once real depth data have been collected, use crossValidate() and
#' interpBathy() on the actual data to check whether observed accuracy matches what was predicted here, and adjust
#' future survey effort accordingly.
#' Runtime scales with n_sim x length(spacings) x (cost of one synthetic-surface simulation + one interpBathy call).
#' For a first look, consider a smaller n_sim (e.g. 5) and a coarser 'res' before committing to a long run.
#' @return a list with:
#' \describe{
#' \item{results}{a data frame of spacing, mean RMSE, and SD of RMSE across simulations}
#' \item{recommended_spacing}{the recommended transect spacing, in meters}
#' \item{transects}{a SpatVector of the recommended transect lines, in the original CRS of 'outline'}
#' \item{total_transect_length}{total length of the recommended transects, in meters}
#' \item{truth_examples}{a multi-layer SpatRaster (one layer per saved example) of synthetic "true" bathymetries used
#' during the simulation, in the original CRS of 'outline', for visually sanity-checking the assumed depth profile}
#' }
#' @author Tristan Blechinger, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @import dplyr
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @import gstat
#' @examples
#' \donttest{
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' samplingDensity(outline, max_depth = 40, n_sim = 5, seed = 123)}
samplingDensity <- function(outline, max_depth, spacings = NULL, n_sim = 10, shape = 1,
truth_model = "Exp", truth_range = NULL, truth_sill = NULL, truth_nugget = 0,
along_track_interval = NULL, orientation = NULL, res = NULL,
method = "IDW", nmax = 20, idp = 2, model = "Sph", psill = NULL, range = NULL,
nugget = 0, kappa = NULL, trend_order = 1, zero_threshold = 0.05,
tolerance = 0.1, n_truth_examples = 3,
plot = TRUE, seed = NULL){
if(!is.null(seed)){
if(!is.numeric(seed))
stop("seed must be numeric")
set.seed(seed)
}
#transform outline shapefile into vector (accepts SpatVector, sf, or anything terra::vect() reads)
if(!inherits(outline, "SpatVector")){
if(inherits(outline, "sf") && requireNamespace("sf", quietly = TRUE)){
outline <- sf::st_zm(outline, drop = TRUE, what = "ZM")
}
outline <- terra::vect(outline)
}
else{
outline <- outline
}
#checks
if(!inherits(outline, "SpatVector"))
stop("outline is not a SpatVector or cannot be transformed")
if(missing(max_depth) || !is.numeric(max_depth) || max_depth <= 0)
stop("max_depth must be specified as a positive numeric value")
if(!method %in% c("IDW", "OK", "UK"))
stop("method misspecified. Please choose 'IDW', 'OK', or 'UK'")
if(!is.numeric(tolerance) || tolerance < 0)
stop("tolerance must be a non-negative numeric value")
#store the original CRS so the recommended transects can be returned in it
original_crs <- terra::crs(outline)
test_crs <- terra::crs(outline)
if(is.na(test_crs) || test_crs == ""){
stop("CRS of 'outline' is unable to be defined.")
}
# Function to determine the best UTM CRS for a given vector
get_best_utm <- function(outline) {
centroid <- terra::centroids(outline)
lon <- terra::crds(centroid)[1]
lat <- terra::crds(centroid)[2]
utm_zone <- base::floor((lon + 180) / 6) + 1
hemisphere <- base::ifelse(lat >= 0, 32600, 32700)
epsg_code <- hemisphere + utm_zone
return(terra::crs(paste0("EPSG:", epsg_code)))
}
#run the whole analysis in meters (UTM), same reasoning as interpBathy/crossValidate
if(terra::is.lonlat(outline)){
best_crs <- get_best_utm(outline)
outline <- terra::project(outline, best_crs)
}
message("Simulating in CRS: ", terra::crs(outline, describe = TRUE)$name)
#characteristic length scale of the waterbody, used to scale several defaults
area <- as.numeric(terra::expanse(outline))
L <- sqrt(area)
if(is.null(spacings)){
spacings <- round(L * c(.02, .05, .1, .2, .4))
spacings <- unique(spacings[spacings > 0])
}
if(!is.numeric(spacings) || any(spacings <= 0))
stop("spacings must be positive numeric values")
spacings <- sort(spacings)
if(is.null(res)){
res <- max(1, min(spacings) / 5)
}
if(is.null(along_track_interval)){
along_track_interval <- max(5, min(spacings) / 10)
}
if(is.null(truth_range)){
truth_range <- L / 4
}
if(is.null(truth_sill)){
truth_sill <- (max_depth / 5)^2
}
#### helper: build a masked raster template over the waterbody at a given resolution
build_grid <- function(outline, res){
ext_o <- terra::ext(outline)
ncol_r <- ceiling((ext_o$xmax - ext_o$xmin) / res)
nrow_r <- ceiling((ext_o$ymax - ext_o$ymin) / res)
empty <- terra::rast(ext_o, ncol = ncol_r, nrow = nrow_r, crs = terra::crs(outline))
ras <- terra::rasterize(outline, empty)
grid <- terra::mask(ras, outline)
return(grid)
}
#### helper: rotate xy points by theta (radians, counterclockwise) about a center point
rotate_pts <- function(xy, theta, center){
xy_c <- sweep(xy, 2, center, FUN = "-")
rot_mat <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), nrow = 2)
rotated <- xy_c %*% rot_mat
sweep(rotated, 2, center, FUN = "+")
}
#### helper: build parallel transect lines across outline at a given spacing/orientation,
#### clipped to the outline polygon
make_transects <- function(outline, spacing, orientation_deg = NULL){
verts <- terra::crds(outline)
center <- colMeans(verts)
if(is.null(orientation_deg)){
pc <- stats::prcomp(scale(verts, scale = FALSE))
theta <- atan2(pc$rotation[2, 1], pc$rotation[1, 1])
} else {
theta <- orientation_deg * pi / 180
}
rot_verts <- rotate_pts(verts, -theta, center)
xr <- range(rot_verts[, 1])
yr <- range(rot_verts[, 2])
pad <- spacing
y_lines <- seq(yr[1] - pad, yr[2] + pad, by = spacing)
line_list <- lapply(y_lines, function(yy){
pts_rot <- rbind(c(xr[1] - pad, yy), c(xr[2] + pad, yy))
pts_real <- rotate_pts(pts_rot, theta, center)
terra::vect(pts_real, type = "lines", crs = terra::crs(outline))
})
lines_all <- do.call(rbind, line_list)
lines_clip <- terra::intersect(lines_all, outline)
return(lines_clip)
}
#### helper: generate one synthetic "true" bathymetry raster
simulate_truth <- function(){
grid <- build_grid(outline, res)
dist_to_shore <- terra::distance(grid, terra::as.lines(outline))
dist_to_shore <- terra::mask(dist_to_shore, outline)
max_dist <- terra::global(dist_to_shore, "max", na.rm = TRUE)[[1]]
bowl <- max_depth * (dist_to_shore / max_dist)^shape
#simulate the random roughness component on a coarse grid sized off its own
#spatial autocorrelation range (truth_range), then resample up to match the
#fine bowl-shape grid. The roughness has no real structure finer than
#truth_range anyway, so this avoids running gstat's sequential Gaussian
#simulation over potentially hundreds of thousands to millions of cells on
#large waterbodies - the simulation cost grows much faster than linearly
#with cell count, so this is the main lever for keeping runtime practical.
coarse_res <- max(res, truth_range / 4)
coarse_grid <- build_grid(outline, coarse_res)
coarse_df <- as.data.frame(coarse_grid, xy = TRUE, na.rm = TRUE)
names(coarse_df)[1:2] <- c("x", "y")
g_dummy <- gstat::gstat(formula = z ~ 1, locations = ~x + y, dummy = TRUE, beta = 0,
model = gstat::vgm(psill = truth_sill, model = truth_model,
range = truth_range, nugget = truth_nugget),
nmax = 30)
sim <- stats::predict(g_dummy, newdata = coarse_df, nsim = 1)
roughness_coarse <- terra::rast(data.frame(x = coarse_df$x, y = coarse_df$y, roughness = sim$sim1),
type = "xyz", crs = terra::crs(outline))
roughness_fine <- terra::resample(roughness_coarse, bowl, method = "bilinear")
true_depth <- bowl + roughness_fine
true_depth <- terra::ifel(true_depth < 0, 0, true_depth)
true_depth <- terra::mask(true_depth, outline)
names(true_depth) <- "true_depth"
return(true_depth)
}
#### helper: simulate a survey at a given spacing against a given truth, reconstruct, and score
score_spacing <- function(truth_raster, spacing){
transects <- make_transects(outline, spacing, orientation)
densified <- terra::densify(transects, interval = along_track_interval)
sample_pts <- terra::as.points(densified)
vals <- terra::extract(truth_raster, sample_pts, ID = FALSE)
xy <- terra::crds(sample_pts)
sample_df <- data.frame(x = xy[, 1], y = xy[, 2], z = vals[[1]])
sample_df <- sample_df[!is.na(sample_df$z), ]
result <- tryCatch({
out <- utils::capture.output(
dem <- suppressMessages(interpBathy(outline, sample_df, x = "x", y = "y", z = "z",
zeros = FALSE, separation = along_track_interval,
res = res, method = method, nmax = nmax, idp = idp,
model = model, psill = psill, range = range,
nugget = nugget, kappa = kappa, trend_order = trend_order,
zero_threshold = zero_threshold))
)
if(method %in% c("OK", "UK")) dem <- dem[["depth"]]
dem_resampled <- terra::resample(dem, truth_raster)
sq_err <- (dem_resampled - truth_raster)^2
sqrt(terra::global(sq_err, "mean", na.rm = TRUE)[[1]])
}, error = function(e){
warning("Reconstruction failed for spacing = ", spacing, ": ", conditionMessage(e))
return(NA_real_)
})
return(result)
}
#### run the simulation: one truth surface per rep, tested against every spacing
n_iter <- n_sim * length(spacings)
pb <- utils::txtProgressBar(min = 0, max = n_iter, style = 3)
iter <- 0
rmse_matrix <- matrix(NA_real_, nrow = n_sim, ncol = length(spacings))
colnames(rmse_matrix) <- as.character(spacings)
n_truth_examples <- min(n_truth_examples, n_sim)
truth_list <- list()
for(i in 1:n_sim){
truth_raster <- simulate_truth()
if(i <= n_truth_examples){
truth_list[[i]] <- truth_raster
}
for(j in seq_along(spacings)){
rmse_matrix[i, j] <- score_spacing(truth_raster, spacings[j])
iter <- iter + 1
utils::setTxtProgressBar(pb, iter)
}
}
close(pb)
results <- data.frame(spacing = spacings,
mean_rmse = apply(rmse_matrix, 2, mean, na.rm = TRUE),
sd_rmse = apply(rmse_matrix, 2, stats::sd, na.rm = TRUE))
best_rmse <- min(results$mean_rmse, na.rm = TRUE)
recommended_spacing <- max(results$spacing[results$mean_rmse <= best_rmse * (1 + tolerance)])
#build the recommended transects once more, cleanly, and return them in the original CRS
transects_final <- make_transects(outline, recommended_spacing, orientation)
total_transect_length <- sum(terra::perim(transects_final))
transects_final <- terra::project(transects_final, original_crs)
outline_orig <- terra::project(outline, original_crs)
truth_examples <- do.call(c, truth_list)
names(truth_examples) <- paste0("sim_", seq_len(n_truth_examples))
truth_examples <- terra::project(truth_examples, original_crs)
if(plot){
graphics::plot(results$spacing, results$mean_rmse, type = "b", pch = 19,
xlab = "Transect spacing (m)", ylab = "RMSE (m)",
main = "Simulated DEM accuracy vs. transect spacing")
graphics::abline(v = recommended_spacing, lty = 2, col = "blue")
graphics::legend("topleft", legend = paste0("Recommended: ", recommended_spacing, " m"),
lty = 2, col = "blue", bty = "n")
terra::plot(outline_orig, main = paste0("Recommended survey design (", recommended_spacing, " m spacing)"))
terra::plot(transects_final, add = TRUE, col = "blue", lwd = 1.5)
terra::plot(truth_examples, col = grDevices::hcl.colors(100, "Blues"),
mar = c(2, 2, 2, 4))
}
return(list(results = results,
recommended_spacing = recommended_spacing,
transects = transects_final,
total_transect_length = total_transect_length,
truth_examples = truth_examples))
}
Try the rLakeHabitat package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.