interpBathy: Interpolate bathymetry

In rLakeHabitat: Interpolate Bathymetry and Quantify Physical Aquatic Habitat

Interpolate bathymetry

Description

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).

Usage

interpBathy(
  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
)

Arguments

outline	shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read (e.g., a file path).
df	dataframe of coordinates and depths for a given waterbody. Coordinates are assumed to be in the same CRS as 'outline'.
x	character giving name of longitude column
y	character giving name of latitude column
z	character giving name of depth column
zeros	logical describing if bounding zeros are needed (FALSE) or provided (TRUE), default = FALSE
separation	number describing distance between points, in meters
res	number describing desired cell resolution in meters, default = 10
method	character describing method of interpolation: Inverse Distance Weighted ("IDW"), Ordinary Kriging ("OK"), or Universal Kriging ("UK"). Default = "IDW"
nmax	numeric value describing number of neighbors used in interpolation, default = 20
idp	numeric value describing inverse distance power value for IDW interpolation
model	character describing type of model used in Ordinary/Universal Kriging, options include 'Sph', 'Exp', 'Gau', 'Mat', default = 'Sph'
psill	numeric value describing the partial sill value for OK/UK interpolation, default = NULL
range	numeric describing distance beyond which there is no spatial correlation in Ordinary/Universal Kriging models, default = NULL
nugget	numeric describing variance at zero distance in Ordinary/Universal Kriging models, default = NULL
kappa	numeric value describing model smoothness, default = NULL
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.
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

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.

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.

Value

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(s)

Tristan Blechinger & Sean Bertalot, Department of Zoology & Physiology, University of Wyoming

Examples

#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)

rLakeHabitat documentation built on July 30, 2026, 5:11 p.m.