R/interpBathy.R

In rLakeHabitat: Interpolate Bathymetry and Quantify Physical Aquatic Habitat

#' Interpolate bathymetry
#'
#' Generate a bathymetric digital elevation model (DEM) for a given waterbody using Inverse Distance Weighting (IDW), Ordinary Kriging (OK), or Universal Kriging (UK) interpolation. For high densities of point data, we recommend rarifying prior to interpolation to improve accuracy and reduce computation time (see rarify function).
#'
#' @param outline shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read (e.g., a file path).
#' @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 zeros logical describing if bounding zeros are needed (FALSE) or provided (TRUE), default = FALSE
#' @param separation number describing distance between points, in meters
#' @param res number describing desired cell resolution in meters, default = 10
#' @details
#' The function automatically detects whether 'outline' (and therefore 'df', which is assumed to share its CRS) is in a geographic (decimal degree) or projected (meters) coordinate system. If geographic, the outline and point data are internally reprojected to their best-fit UTM zone so that all distance-based calculations (resolution, nmax neighbor selection, IDW power, kriging variogram parameters, and boundary point separation) operate on meters rather than degrees. The final DEM is reprojected back to the original CRS of 'outline' before being returned. The CRS used for interpolation, and progress through the major steps, are printed/reported as the function runs.
#' 'res' is required and is always in meters, regardless of the CRS 'outline' was originally supplied in.
#' @param method character describing method of interpolation: Inverse Distance Weighted ("IDW"), Ordinary Kriging ("OK"), or Universal Kriging ("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/Universal Kriging, options include 'Sph', 'Exp', 'Gau', 'Mat', default = 'Sph'
#' @param psill numeric value describing the partial sill value for OK/UK interpolation, default = NULL
#' @param range numeric describing distance beyond which there is no spatial correlation in Ordinary/Universal Kriging models, default = NULL
#' @param nugget numeric describing variance at zero distance in Ordinary/Universal Kriging models, default = NULL
#' @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.
#' @details
#' For the model argument there are four different methods included here that are supported by gstat::vgm ("Sph", "Exp", "Gau", "Mat").
#' "Sph" = The default gstat::vgm method. Spherical model characterized by a curve that rises steeply to defined range then flattens, indicates no spatial correlation between points beyond that range.
#' "Exp" = Exponential model characterized by spatial correlation decaying rapidly with distance, results in a rougher surface.
#' "Gau" = Gaussian model similar to spatial model but with slower decay over distance, results in a smoother surface.
#' "Mat" = Matern model that uses kappa to define the variogram relationship. High kappa values approach a Guassian model (smooth surface), and low kappa values approach the Exponential model (kappa = 0.5 is equivalent to Exponential).
#' Three parameters (psill, range, kappa) are incorporated from a fitted variogram (default = NULL). If specified in function input, chosen values will overwrite variogram values - and any parameter that is auto-fit is fit with knowledge of the others you did supply (including nugget), rather than fitting as if the rest were still at their gstat defaults.
#' Universal Kriging ("UK") differs from Ordinary Kriging in that it fits a polynomial trend surface across x/y (see
#' 'trend_order') and models spatial correlation in the residuals from that trend, rather than assuming a constant mean
#' across the whole waterbody. This can help for reservoirs with a strong directional depth gradient (e.g. a river-fed
#' arm sloping steadily toward a dam), where OK's constant-mean assumption doesn't hold well.
#'
#' @return the interpolated DEM. For "IDW", a single-layer SpatRaster. For "OK" and "UK", a two-layer SpatRaster:
#' layer 'depth' (the interpolated values) and layer 'error' (the associated standard error of each estimate).
#' @author Tristan Blechinger & Sean Bertalot, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @import dplyr
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @import gstat
#' @examples
#' \donttest{
#' #load example outline
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' #load example xyz data
#' data <- read.csv(system.file("extdata", "example_depths.csv", package = 'rLakeHabitat'))
#' #run function
#' interpBathy(outline, data, "x", "y", "z", zeros = FALSE, separation = 10,
#' res = 5, method = "IDW", nmax = 4, idp = 2)}

interpBathy <- function(outline, df, x, y, z, zeros = FALSE, separation = NULL, res = 10, method = "IDW", nmax = 20, idp = 2, model = "Sph", psill = NULL, range = NULL, nugget = NULL, kappa = NULL, trend_order = 1, zero_threshold = 0.05){

  # transform outline shapefile into vector
  # accepts a SpatVector, sf object, or anything else terra::vect() can read (e.g. a file path to a shapefile)
  if(!inherits(outline, "SpatVector")){
    if(inherits(outline, "sf") && requireNamespace("sf", quietly = TRUE)){
      # drop Z/M dimensions if present
      outline <- sf::st_zm(outline, drop = TRUE, what = "ZM")
    }
    outline <- terra::vect(outline)
  }
  else{
    outline <- outline
  }

  # store original CRS for final DEM reprojection
  original_crs <- terra::crs(outline)

  # report progress
  pb_steps <- 5
  pb <- utils::txtProgressBar(min = 0, max = pb_steps, style = 3)
  pb_step <- 0
  update_pb <- function(){
    pb_step <<- pb_step + 1
    utils::setTxtProgressBar(pb, pb_step)
  }
  on.exit(close(pb), add = TRUE)

  old_terra_opts <- tryCatch(terra::terraOptions(print = FALSE), error = function(e) NULL)
  old_progress <- if(!is.null(old_terra_opts) && !is.null(old_terra_opts$progress)) old_terra_opts$progress else 3
  terra::terraOptions(progress = 1)
  on.exit(terra::terraOptions(progress = old_progress), add = TRUE)

  # checks
  if(!inherits(df, "data.frame"))
    stop("df must be a dataframe")
  if(!inherits(x, "character"))
    stop("x must be a character giving the longitude column name")
  if(!inherits(y, "character"))
    stop("y must be a character giving the latitude column name")
  if(!inherits(z, "character"))
    stop("z must be a character giving the depth column name")
  if(x %in% names(df) == FALSE)
    stop("The value of x does not appear to be a valid column name")
  if(y %in% names(df) == FALSE)
    stop("The value of y does not appear to be a valid column name")
  if(z %in% names(df) == FALSE)
    stop("The value of z does not appear to be a valid column name")
  if(!inherits(df[, x], "numeric"))
    stop("data in x column is not formatted as numeric")
  if(!inherits(df[, y], "numeric"))
    stop("data in y column is not formatted as numeric")
  if(!inherits(df[, z], "numeric"))
    stop("data in z column is not formatted as numeric")
  if(!inherits(outline, "SpatVector"))
    stop("outline is not a SpatVector or cannot be transformed")
  if(!is.logical(zeros))
    stop("zeros must be either 'T', 'F', TRUE, or FALSE")
  if(zeros == F){
    if(is.null(separation) || is.na(separation))
      stop("separation value must be specified")
  }
  if(zeros == T){
    if(!is.null(separation))
      stop("separation must be null if zeros = T")
  }
  if(is.null(res) || !is.numeric(res)){
    stop("res must be specified as a numeric value")
  }
  if(!method %in% c("IDW", "OK", "UK"))
    stop("method misspecified. Please choose 'IDW', 'OK', or 'UK'")
  if(method == "IDW"){
    if(!is.numeric(nmax))
      stop("nmax must be numeric")
    if(!is.numeric(idp))
      stop("idp must be numeric")
    if(is.null(nmax) || is.na(nmax))
      stop("nmax must be specified")
    if(is.null(idp) || is.na(idp))
      stop("idp must be specified")
    if(nmax > nrow(df))
      stop("nmax cannot exceed number of observations in df")
    model <- NULL
    psill <- NULL
    range <- NULL
    nugget <- NULL
    kappa <- NULL
  }
  if(method %in% c("OK", "UK")){
    if(!is.character(model) || !model %in% c("Sph", "Exp", "Gau", "Mat"))
      stop("model must be character string of either 'Sph', 'Exp', 'Gau', or 'Mat'")
    if(!is.numeric(nmax))
      stop("nmax must be numeric")
    if(nmax > nrow(df))
      stop("nmax cannot exceed number of observations in df")
    if(!is.null(nugget)){
      if(!is.numeric(nugget))
        stop("nugget must be numeric")
    }
    else{
      nugget <- NA
    }
    if(!is.null(range)){
      if(!is.numeric(range))
        stop("range must be numeric")
    }
    else{
      range <- NA
    }
    if(!is.null(psill)){
      if(!is.numeric(psill))
        stop("psill must be numeric")
    }
    else{
      psill <- NA
    }
    if(!is.null(kappa)){
      if(!is.numeric(kappa))
        stop("kappa must be numeric")
    }
    else{
      kappa <- NA
    }
    idp <- NULL
  }
  if(method == "UK"){
    if(!trend_order %in% c(1, 2))
      stop("trend_order must be 1 or 2")
  }
  test_crs <- terra::crs(outline)
  if(is.na(test_crs) || test_crs == ""){
    stop("CRS of 'outline' is unable to be defined.")
  }


  #### Reproject outline to UTM meters for interpolation

  # Function to determine the best UTM CRS for a given vector
  get_best_utm <- function(outline) {
    # Get centroid of the input vector (assuming it's a SpatVector or SpatRaster)
    centroid <- terra::centroids(outline) # Get a sample point

    # Extract longitude and latitude
    lon <- terra::crds(centroid)[1]
    lat <- terra::crds(centroid)[2]

    # Compute UTM zone
    utm_zone <- base::floor((lon + 180) / 6) + 1

    # Determine Northern or Southern Hemisphere
    hemisphere <- base::ifelse(lat >= 0, 32600, 32700)  # 326xx for North, 327xx for South

    # Construct the EPSG code
    epsg_code <- hemisphere + utm_zone

    # Return the CRS in terra format
    return(terra::crs(paste0("EPSG:", epsg_code)))
  }

  # Automatically detect whether outline is in a geographic (decimal degree) or
  # already-projected (meters) CRS. If geographic, project to its best-fit UTM
  # zone so distance-based calculations run in true meters, not degrees.
  if(terra::is.lonlat(outline)){
    best_crs <- get_best_utm(outline)
    outline <- terra::project(outline, best_crs)
  }

  message("Interpolating in CRS: ", terra::crs(outline, describe = TRUE)$name)
  update_pb()

  # Identify rows/columns needed for desired resolution
  get_res <- function(outline, res) {
    ext <- terra::ext(outline) # gets extent, now in meters
    ext_length <- base::abs(ext$xmin - ext$xmax)
    ext_height <- base::abs(ext$ymax - ext$ymin) # Finds height of extent in m

    set_ext_x <- ext_length / res # Divides by res
    set_ext_y <- ext_height / res

    xy <- c(set_ext_x, set_ext_y) # List of length and height
    return(xy)
  }

  if(!is.null(res)){
    xy <- get_res(outline, res)

    empty_raster <- terra::rast(ext(outline), ncol = xy[1], nrow = xy[2], crs = terra::crs(outline))
  }

  # select and order df columns
  df <- df %>%
    dplyr::select(all_of(c(x, y, z))) %>%
    dplyr::rename(x = x, y = y, z = z)

  # if outline was reprojected, reproject point data to match
  if(!identical(terra::crs(outline), original_crs)){
    pts <- terra::vect(df, geom = c("x", "y"), crs = original_crs)
    pts <- terra::project(pts, terra::crs(outline))
    proj_coords <- terra::crds(pts)
    df$x <- proj_coords[, 1]
    df$y <- proj_coords[, 2]
  }

  # add bounding zeros to dataframe if not included
  if(zeros == F){
    #segment line
    line_segmented <- terra::densify(outline, interval = separation)

    # convert segmented line to points
    points <- terra::as.points(line_segmented)

    # convert points to coordinates, forcing every boundary point to z = 0
    # (this must not use terra::geom()'s 'hole' column as a stand-in for depth -
    # hole is a ring-type flag (0 = outer boundary, 1 = an interior ring, i.e.
    # an island), not a depth value)
    coords <- terra::geom(points)
    boundary_zeros <- as.data.frame(coords)
    boundary_zeros <- boundary_zeros %>% dplyr::select(all_of(c("x", "y")))
    boundary_zeros$z <- 0

    df <- base::rbind(df, boundary_zeros)
  }
  else{
    df <- df
  }

  # generate empty raster of waterbody at the requested resolution
  disagg.ras <- empty_raster

  update_pb()

  # creates a raster of the shape outline in the grid dissagg.ras
  ras <- terra::rasterize(outline, disagg.ras)

  # masks the raster for the shapefile (everything outside the reservoir = NA)
  grid <- terra::mask(ras, outline)

  # inverse distance weighted interpolation
  if (method == "IDW") {

    message("Running IDW interpolation in CRS: ", terra::crs(outline, describe = TRUE)$name)

    # Interpolation function for deriving contours
    gs <- gstat::gstat(formula = z ~ 1,
                       locations = ~x + y,
                       data = df,
                       nmax = nmax,
                       set = list(idp = idp))

    # Remove NA cells from the grid
    grid <- terra::na.omit(grid)

    # Create DEM with interpolate function
    DEM <- terra::interpolate(grid, gs)
    update_pb()

    # Mask to lake interior only
    lake_interior <- terra::mask(DEM, outline)
    lake_interior <- lake_interior[["var1.pred"]]

    # Identify zero values, and what fraction of the waterbody's surface area they represent
    r_zero <- terra::ifel(lake_interior == 0, 0, NA)
    total_cells <- terra::global(lake_interior, fun = "notNA")$notNA
    zero_cells_n <- terra::global(r_zero, fun = "notNA")$notNA
    zero_fraction <- if(total_cells > 0) zero_cells_n / total_cells else 0

    message(zero_cells_n, " of ", total_cells, " interior cells (", signif(zero_fraction * 100, 3),
            "%) interpolated to exactly 0; zero_threshold = ", zero_threshold * 100, "%. ",
            if(zero_fraction >= zero_threshold) "Re-interpolating flagged cells." else "Leaving as-is.")

    if (zero_fraction >= zero_threshold) {

      # Convert zeros to NA in lake_interior
      lake_interior[lake_interior == 0] <- NA
      lake_interior_clean <- lake_interior

      # Use valid (non-zero, non-NA) points as input
      valid_points <- as.data.frame(lake_interior_clean, xy = TRUE, na.rm = TRUE)
      names(valid_points) <- c("x", "y", "z")

      sample_size <- min(nrow(df) * 2, nrow(valid_points))
      valid_points <- valid_points[sample(nrow(valid_points), sample_size), ]

      # Second IDW pass with valid points
      gs2 <- gstat::gstat(formula = z ~ 1,
                          locations = ~x + y,
                          data = valid_points,
                          nmax = nmax,
                          set = list(idp = idp))

      # Grid of just the flagged cells to predict to (lake_interior is
      # already NA exactly at those cells)
      DEM_filled <- terra::interpolate(lake_interior, gs2)

      # Merge filled values back into original
      DEM <- terra::cover(lake_interior, DEM_filled)
    } else {
      DEM <- lake_interior
    }
    update_pb()

    # Mask to reintroduce NA's where there's no water
    mask.na <- terra::init(grid, NA)
    # Create a grid with only 0 where the land is
    mask.2 <- terra::init(grid, 0)

    # Replace NA values with the max height
    replace.na <- terra::cover(grid, mask.2, values = NA)
    # Fill the lake with NA values
    replace.water <- terra::cover(replace.na, mask.na, values = 1)

    # Merge the two rasters
    final_DEM <- terra::merge(replace.water, DEM)

    # Isolate layer 1, bathymetry
    final_DEM <- final_DEM[[1]]
    final_DEM[[1]][final_DEM[[1]] < 0] <- 0
    # Remove all external zeros from original lake shape
    final_DEM <- terra::mask(final_DEM, outline)

    # Reproject the DEM back to the CRS the user originally supplied
    final_DEM <- terra::project(final_DEM, original_crs)
    update_pb()

    print(paste("Raster Interpolated in CRS: ", terra::crs(final_DEM), sep = " "))

    return(final_DEM)
  }

  # Ordinary or Universal Kriging interpolation
  if (method %in% c("OK", "UK")) {

    message("Running ", method, " interpolation in CRS: ", terra::crs(outline, describe = TRUE)$name)

    # OK assumes a constant (unknown) mean everywhere; UK instead fits a
    # polynomial trend surface across x/y and krige's the residuals from it
    if(method == "OK"){
      krige_formula <- z ~ 1
    }
    if(method == "UK"){
      if(trend_order == 1){
        krige_formula <- z ~ x + y
      }
      if(trend_order == 2){
        krige_formula <- z ~ x + y + I(x^2) + I(x*y) + I(y^2)
      }
    }

    # Run and save variogram output, informed by any explicitly supplied
    # nugget/range/psill/kappa rather than fitting as if they were still at
    # gstat's defaults
    if (is.na(psill) || is.na(range) || is.na(kappa) || is.na(nugget)) {
      vgram <- gstat::variogram(krige_formula, locations = ~x + y, data = df)
      # use 0 as the starting value for nugget if it's being auto-fit;
      # if the user supplied a fixed nugget, that value is used as-is and
      # is not overwritten below
      init_vgm <- gstat::vgm(psill = psill, model = model, range = range,
                             nugget = if(is.na(nugget)) 0 else nugget, kappa = kappa)
      gramParam <- gstat::fit.variogram(vgram, model = init_vgm)

      # fit.variogram() returns one row per structure: a "Nug" row (nugget)
      # and a row for the requested model (psill/range/kappa) - always pull
      # each parameter from its own row, not the whole column
      model_row <- which(gramParam$model != "Nug")
      if(length(model_row) == 0) model_row <- nrow(gramParam)
      nug_row <- which(gramParam$model == "Nug")
    }
    if (is.na(psill)) {
      psill <- gramParam$psill[model_row]
    }
    if (is.na(range)) {
      range <- gramParam$range[model_row]
    }
    if (is.na(kappa)) {
      kappa <- gramParam$kappa[model_row]
    }
    if (is.na(nugget)) {
      nugget <- if(length(nug_row) > 0) gramParam$psill[nug_row] else 0
    }

    # Define bathymetry model
    gs <- gstat::gstat(formula = krige_formula,
                       locations = ~x + y,
                       data = df,
                       nmax = nmax,
                       model = gstat::vgm(model = model,
                                          range = range,
                                          nugget = nugget,
                                          psill = psill,
                                          kappa = kappa))

    # Remove NA cells from the grid
    grid <- terra::na.omit(grid)

    # Create DEM with interpolate function
    DEM <- terra::interpolate(grid, gs)
    update_pb()

    # Mask to lake interior only
    lake_interior <- terra::mask(DEM[[1]], outline)

    # Identify zero values, and what fraction of the waterbody's surface area they represent
    r_zero <- terra::ifel(lake_interior == 0, 0, NA)
    total_cells <- terra::global(lake_interior, fun = "notNA")$notNA
    zero_cells_n <- terra::global(r_zero, fun = "notNA")$notNA
    zero_fraction <- if(total_cells > 0) zero_cells_n / total_cells else 0

    message(zero_cells_n, " of ", total_cells, " interior cells (", signif(zero_fraction * 100, 3),
            "%) interpolated to exactly 0; zero_threshold = ", zero_threshold * 100, "%. ",
            if(zero_fraction >= zero_threshold) "Re-interpolating flagged cells." else "Leaving as-is.")

    if (zero_fraction >= zero_threshold) {

      # Convert zeros to NA in lake_interior
      lake_interior[lake_interior == 0] <- NA
      lake_interior_clean <- lake_interior

      # Use valid (non-zero, non-NA) points as input
      valid_points <- as.data.frame(lake_interior_clean, xy = TRUE, na.rm = TRUE)
      names(valid_points) <- c("x", "y", "z")

      sample_size <- min(nrow(df) * 2, nrow(valid_points))
      valid_points <- valid_points[sample(nrow(valid_points), sample_size), ]

      # Second kriging pass with valid points
      gs2 <- gstat::gstat(formula = krige_formula,
                          locations = ~x + y,
                          data = valid_points,
                          nmax = nmax,
                          model = gstat::vgm(model = model,
                                             range = range,
                                             nugget = nugget,
                                             psill = psill,
                                             kappa = kappa))

      # Grid of just the flagged cells to predict to (lake_interior is
      # already NA exactly at those cells)
      DEM_filled <- terra::interpolate(lake_interior, gs2)

      # Merge filled values back into original
      lake_interior <- terra::cover(lake_interior_clean, DEM_filled[["var1.pred"]])

      var_clean <- terra::mask(DEM[[2]], lake_interior_clean)
      var_combined <- terra::cover(var_clean, DEM_filled[["var1.var"]])

      DEM <- c(lake_interior, var_combined)
      names(DEM) <- c("var1.pred", "var1.var")
    }
    update_pb()

    # Mask to reintroduce NA's where there's no water
    mask.na <- terra::init(grid, NA)
    # Create a grid with only 0 where the land is
    mask.2 <- terra::init(grid, 0)

    # Replace NA values with the max height
    replace.na <- terra::cover(grid, mask.2, values = NA)
    # Fill the lake with NA values
    replace.water <- terra::cover(replace.na, mask.na, values = 1)

    # Get error raster
    error <- sqrt(DEM[[2]])

    # Merge the two rasters
    final_DEM <- terra::merge(replace.water, DEM)
    final_DEM_error <- terra::merge(replace.water, error)

    # Isolate layer 1, bathymetry
    final_DEM[[1]][final_DEM[[1]] < 0] <- 0

    final_DEM <- final_DEM[[1]]
    final_DEM_error <- final_DEM_error[[1]]

    # Remove all external zeros from original lake outline
    final_DEM <- terra::mask(final_DEM, outline)
    final_DEM_error <- terra::mask(final_DEM_error, outline)

    # Reproject the DEM and error raster back to the CRS the user originally supplied
    final_DEM <- terra::project(final_DEM, original_crs)
    final_DEM_error <- terra::project(final_DEM_error, original_crs)
    update_pb()

    final_stack <- c(final_DEM, final_DEM_error)
    names(final_stack) <- c("depth", "error")

    print(paste("Raster Interpolated in CRS: ", terra::crs(final_DEM), sep = " "))

    message("Kriging parameters used - method: ", method, ", model: ", model,
            ", nugget: ", signif(nugget, 4), ", psill: ", signif(psill, 4),
            ", range: ", signif(range, 4),
            if(method == "UK") paste0(", trend_order: ", trend_order) else "")

    return(final_stack)
  }
}