R/optimizeParams.R

In rLakeHabitat: Interpolate Bathymetry and Quantify Physical Aquatic Habitat

#' Optimize Point Density and DEM Resolution From Existing Depth Data
#'
#' Given a dense depth dataset, determine what point spacing the data can be rarified to, and what DEM cell resolution to interpolate at, without meaningfully losing accuracy. Unlike samplingDensity(), this works directly with measuremed depth data: at every combination of candidate rarefaction spacing and candidate DEM resolution, the data are thinned with rarify() and scored with crossValidate(), building an RMSE surface across the full spacing x resolution grid.
#'
#' @param outline shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read.
#' @param df dataframe of coordinates and depths for a given waterbody. Coordinates are assumed to be in the same CRS as 'outline'.
#' @param x character giving name of longitude column
#' @param y character giving name of latitude column
#' @param z character giving name of depth column
#' @param spacings numeric vector of candidate rarify() target spacings to test, in meters.
#' @param res_values numeric vector of candidate DEM cell resolutions to test, in meters.
#' @param k numeric value describing the number of cross-validation folds to use at each grid point, default = 5
#' @param zeros logical describing if bounding zeros are needed (FALSE) or provided (TRUE), default = FALSE
#' @param separation number describing distance between boundary points, in meters (required if zeros = FALSE)
#' @param method character describing method of interpolation, "IDW", "OK", or "UK. Default = "IDW"
#' @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, resolution) pair: among all grid points whose mean RMSE is within 'tolerance' of the single best (lowest-RMSE) grid point, the coarsest spacing is chosen (ties broken by coarsest resolution). Default = 0.1 (10%).
#' @param plot logical: should a diagnostic heatmap of RMSE across the spacing x resolution grid be drawn? Default = TRUE.
#' @param seed optional numeric value used to seed the random number generator (fold assignment in crossValidate() is stochastic), for reproducible results across runs. Default = NULL (not seeded).
#' @details
#' Runtime scales with length(spacings) x length(res_values) x k, since each grid point runs a full crossValidate() call (which itself fits k separate DEMs). Start with a small grid (e.g. 3x3) to confirm the analysis runs end-to-end before scaling up to a finer search.
#' @return a list with:
#' \describe{
#'   \item{results}{a data frame of spacing, res, mean RMSE, and the number of points remaining after rarefaction at each spacing}
#'   \item{recommended_spacing}{the recommended rarefaction spacing, in meters}
#'   \item{recommended_res}{the recommended DEM cell resolution, in meters}
#' }
#' @author Tristan Blechinger, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @import dplyr
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @examples
#' \donttest{
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' data <- read.csv(system.file("extdata", "example_depths.csv", package = 'rLakeHabitat'))
#' optimizeParams(outline, data, "x", "y", "z", spacings = c(10, 25, 50),
#' res_values = c(5, 10, 20), seed = 123)}

optimizeParams <- function(outline, df, x, y, z, spacings, res_values, k = 5, zeros = FALSE, separation = NULL,
                             method = "IDW", nmax = 20, idp = 2, model = "Sph", psill = NULL, range = NULL,
                             nugget = NULL, kappa = NULL, trend_order = 1, zero_threshold = 0.05, tolerance = 0.1, plot = TRUE, seed = NULL){

  if(!is.null(seed)){
    if(!is.numeric(seed))
      stop("seed must be numeric")
    set.seed(seed)
  }

  #checks
  if(!inherits(df, "data.frame"))
    stop("df must be a dataframe")
  if(!is.numeric(spacings) || any(spacings <= 0))
    stop("spacings must be positive numeric values")
  if(!is.numeric(res_values) || any(res_values <= 0))
    stop("res_values must be positive numeric values")
  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")

  spacings <- sort(unique(spacings))
  res_values <- sort(unique(res_values))

  n_iter <- length(spacings) * length(res_values)
  pb <- utils::txtProgressBar(min = 0, max = n_iter, style = 3)
  iter <- 0

  grid_df <- expand.grid(spacing = spacings, res = res_values, KEEP.OUT.ATTRS = FALSE)
  grid_df$mean_rmse <- NA_real_
  grid_df$n_points <- NA_integer_

  for(s in spacings){

    rarified <- tryCatch(
      rarify(outline, df, x = x, y = y, z = z, res = s),
      error = function(e){
        warning("rarify() failed for spacing = ", s, ": ", conditionMessage(e))
        return(NULL)
      }
    )

    n_pts <- if(is.null(rarified)) NA_integer_ else nrow(rarified)

    for(r in res_values){

      rmse_val <- NA_real_

      if(!is.null(rarified)){
        rmse_val <- tryCatch({
          out <- utils::capture.output(
            cv <- crossValidate(outline, rarified, x = "x", y = "y", z = "z", k = k,
                                zeros = zeros, separation = separation, res = r,
                                method = method, nmax = nmax, idp = idp, model = model,
                                psill = psill, range = range, nugget = nugget, kappa = kappa,
                                zero_threshold = zero_threshold, trend_order = trend_order)
          )
          cv[["RMSE"]]
        }, error = function(e){
          warning("crossValidate() failed for spacing = ", s, ", res = ", r, ": ", conditionMessage(e))
          return(NA_real_)
        })
      }

      idx <- which(grid_df$spacing == s & grid_df$res == r)
      grid_df$mean_rmse[idx] <- rmse_val
      grid_df$n_points[idx] <- n_pts

      iter <- iter + 1
      utils::setTxtProgressBar(pb, iter)
    }
  }
  close(pb)

  best_rmse <- min(grid_df$mean_rmse, na.rm = TRUE)
  within_tol <- grid_df[!is.na(grid_df$mean_rmse) & grid_df$mean_rmse <= best_rmse * (1 + tolerance), ]
  within_tol <- within_tol[order(-within_tol$spacing, -within_tol$res), ]
  recommended_spacing <- within_tol$spacing[1]
  recommended_res <- within_tol$res[1]

  if(plot){
    rmse_matrix <- matrix(grid_df$mean_rmse, nrow = length(spacings), ncol = length(res_values))
    graphics::image(x = spacings, y = res_values, z = rmse_matrix,
                    xlab = "Rarefaction spacing (m)", ylab = "DEM resolution (m)",
                    main = "Cross-validated RMSE across spacing x resolution",
                    col = grDevices::heat.colors(20))
    graphics::points(recommended_spacing, recommended_res, pch = 4, cex = 2, lwd = 2, col = "blue")
    # graphics::legend("topright", legend = "recommended", pch = 4, col = "blue", bty = "n")
  }

  return(list(results = grid_df,
              recommended_spacing = recommended_spacing,
              recommended_res = recommended_res))
}