Type: Package
Title: Interpolate Bathymetry and Quantify Physical Aquatic Habitat
Version: 1.0.1
Date: 2025-11-06
URL: https://gitlab.com/tristanblechinger/rlakehabitat
BugReports: https://gitlab.com/tristanblechinger/rlakehabitat/-/issues
Depends: R (≥ 4.3.0)
Imports: dplyr, terra, gstat, sf, ggplot2, gganimate, tidyterra, rLakeAnalyzer, isoband
Maintainer: Tristan Blechinger <tblechin@uwyo.edu>
Description: Offers bathymetric interpolation using Inverse Distance Weighted and Ordinary Kriging via the 'gstat' and 'terra' packages. Other functions focus on quantifying physical aquatic habitats (e.g., littoral, epliminion, metalimnion, hypolimnion) from interpolated digital elevation models (DEMs). Functions were designed to calculate these metrics across water levels for use in reservoirs but can be applied to any DEM and will provide values for fixed conditions. Parameters like Secchi disk depth or estimated photic zone, thermocline depth, and water level fluctuation depth are included in most functions.
License: GPL (≥ 3)
VignetteBuilder: knitr
Encoding: UTF-8
RoxygenNote: 7.3.3
Suggests: testthat (≥ 3.0.0), knitr, rmarkdown
Config/testthat/edition: 3
NeedsCompilation: no
Packaged: 2025-11-06 22:09:29 UTC; tblechin
Author: Tristan Blechinger ORCID iD [aut, cre], Sean Bertalot ORCID iD [aut]
Repository: CRAN
Date/Publication: 2025-11-07 10:00:13 UTC

Generate Animated Plot

Description

Generate an animated plot of littoral area at different water level increments from a raster digital elevation model (DEM).

Usage

animBathy(
  DEM,
  units = "ft",
  littoral = TRUE,
  secchi = NULL,
  photic = NULL,
  stop = NULL,
  by = 1
)

Arguments

DEM

SpatRaster object of a given waterbody, rasters can be transformed to SpatRaster via the rast() function in 'terra'

units

character describing depth units of DEM. Can be meters ("m") or feet ("ft"). Default = "ft"

littoral

logical indicating if littoral zone should be plotted (T) or entire waterbody (F), default = TRUE

secchi

number giving the average Secchi depth of the waterbody, photic zone estimated as 2.6m * secchi

photic

number giving the average photic depth of the waterbody, overwrites Secchi

stop

optional numeric value specifying depth at which to stop animation, default = NULL (all depths)

by

numeric value specifying depth increments between plots. Higher values will result in lower resolution. Default = 1

Value

an animated ggplot object

Author(s)

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

Examples


#load raster
DEM <- terra::rast(system.file("extdata", "example_raster.tif", package = 'rLakeHabitat'))

#run function
animBathy(DEM, units = 'm', littoral = TRUE, secchi = 1, by = 5)

Plot Bathymetry Map

Description

Generate a bathymetry map from a provided DEM raster with optional contours and depth labels.

Usage

bathyMap(
  DEM,
  contours = TRUE,
  start = NULL,
  end = NULL,
  by = 5,
  breaks = NULL,
  units = "ft",
  labels = TRUE,
  textSize = 1.5,
  plotTitle = NULL
)

Arguments

DEM

SpatRaster object of a given waterbody, rasters can be transformed to SpatRaster via the rast() function in 'terra'

contours

logical indicating whether contours should included (TRUE) or not (FALSE), default = TRUE

start

numeric value describing what value contours should start at, default = 0

end

numeric value describing what value contours should end at, default = max depth

by

numeric value describing contour intervals, default = 5

breaks

optional numeric vector describing specific contours to include if contours = T, default = NULL

units

character describing units of depth measurement, default = "ft"

labels

logical indicating whether labels should be included (TRUE) or not (FALSE), default = TRUE

textSize

number describing text size of contour labels if included, default = 1.5

plotTitle

optional character string adding title to output plot

Value

ggplot object

Author(s)

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

Examples

#load raster
DEM <- terra::rast(system.file("extdata", "example_raster.tif", package = 'rLakeHabitat'))
#run function
bathyMap(DEM, contours = TRUE, units = 'm', labels = TRUE)

Calculate Hypsography

Description

Calculates area at each depth for a given waterbody.

Usage

calcHyps(DEM, DEMunits = "m", depthUnits = "ft", by = 1, output = "values")

Arguments

DEM

SpatRaster object of a given waterbody, rasters can be transformed to SpatRaster via the rast() function in 'terra'

DEMunits

character describing units of raster coordinate system. Can be meters, kilometers, or hectares ("m", "km", "ha"), default = "m"

depthUnits

character describing units of depth measurement. Can be either feet or meters ("ft", "m"), default = "ft"

by

numeric increment per unit by which volumes are calculated. Higher values will result in lower resolution. Default = 1

output

character describing desired output, can either be a data frame of values ("values") or a hypsography plot ("plot"). Default = "values"

Value

data frame of areas at each depth unit ("values") or a hypsography plot ("plot")

Author(s)

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

Examples

#load raster
DEM <- terra::rast(system.file("extdata", "example_raster.tif", package = 'rLakeHabitat'))
#run function
calcHyps(DEM, DEMunits = 'm', depthUnits = 'm', by = 1, output = 'values')

Calculate Littoral Area

Description

Calculates littoral surface area (2D) of a given waterbody across water levels based on an average photic depth value.

Usage

calcLittoral(
  DEM,
  photic = NULL,
  secchi = NULL,
  DEMunits = "m",
  depthUnits = "ft",
  by = 1,
  stop = NULL
)

Arguments

DEM

SpatRaster object of a given waterbody, rasters can be transformed to SpatRaster via the rast() function in 'terra'

photic

number giving the average photic depth, overwrites Secchi depth

secchi

number giving the average secchi depth, photic zone estimated as 2.6m * secchi.

DEMunits

character describing units of raster coordinate system. Can be meters, kilometers, or hectares ("m", "km", "ha"), default = "m"

depthUnits

character describing units of depth measurement (secchi and DEM). Can be either feet or meters ("ft", "m"), default = "ft"

by

numeric increment per unit depth by which areas are calculated. Higher values will result in lower resolution. Default = 1

stop

optional numeric value specifying depth at which to stop calculations, default = NULL

Value

data frame of areas in specified units for each depth, as well as the littoral percentage of total surface area

Author(s)

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

Examples

#load raster
DEM <- terra::rast(system.file("extdata", "example_raster.tif", package = 'rLakeHabitat'))
#run function
calcLittoral(DEM, secchi = 1, depthUnits = "m", DEMunits = "m")

Calculate Shoreline Development Index

Description

Calculates Shoreline Development Index value across water levels for a given waterbody.

Usage

calcSDI(DEM, units = "m", by = 1, stop = NULL)

Arguments

DEM

SpatRaster object of a given waterbody, rasters can be transformed to SpatRaster via the rast() function in 'terra'

units

character describing units of raster coordinate system. Can be meters, kilometers, or hectares ("m", "km", "ha"), default = "m"

by

numeric increment per unit depth by which areas are calculated. Higher values will result in lower resolution. Default = 1

stop

optional numeric value specifying depth at which to stop calculations, default = NULL

Value

data frame of perimeter lengths and SDI values for given depths

Author(s)

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

Examples

#load raster
DEM <- terra::rast(system.file("extdata", "example_raster.tif", package = 'rLakeHabitat'))
#run function
calcSDI(DEM, units = 'm')

Calculate Pelagic Habitat Volumes

Description

Calculates epilimnion, metalimnion, and hypolimnion volumes based on defined thermocline depths across water levels.

Usage

calcVolume(
  DEM,
  thermo_depth = NULL,
  thermo_high,
  thermo_low,
  DEMunits = "m",
  depthUnits = "ft",
  by = 1,
  stop = NULL
)

Arguments

DEM

SpatRaster object of a given waterbody, rasters can be transformed to SpatRaster via the rast() function in 'terra'

thermo_depth

number giving the estimated middle of thermocline, results in calculation of only epilimnion and hypolimnion volumes. Default = NULL, cannot use in conjunction with thermo_low and thermo_high

thermo_high

number giving the upper bound of thermocline depth, results in calculation of epilimnion, metalimnion, and hypolimnion values

thermo_low

number giving the lower bound of thermocline depth, results in calculation of epilimnion, metalimnion, and hypolimnion values

DEMunits

character describing units of raster coordinate system. Can be meters, kilometers, or hectares ("m", "km", "ha"), default = "m"

depthUnits

character describing units of depth measurement. Can be either feet or meters ("ft", "m"), default = "ft"

by

numeric increment per unit by which volumes are calculated. Higher values will result in lower resolution. Default = 1

stop

optional numeric value specifying depth at which to stop habitat volume calculations, default = NULL

Value

a data frame of volumes in cubic meters calculated for each habitat (epilimnion, metalimnion, hypolimnion)

Author(s)

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

Examples

#load raster
DEM <- terra::rast(system.file("extdata", "example_raster.tif", package = 'rLakeHabitat'))
#run function
calcVolume(DEM, thermo_depth = 3, DEMunits = 'm', depthUnits = 'm')

Contour Lines to Points

Description

Get point coordinates and depth values along predetermined contours at a specified density.

Usage

contourPoints(object, depths = NULL, geometry = "geometry", density = 10)

Arguments

object

polygon or multipolygon shapefile (.shp) with depths included as an attribute column. Can be an sf or spatVector object.

depths

character string describing column name of depth attribute

geometry

character string describing column name of geometries. Default = "geometry"

density

numeric value describing distance between points in meters, default = 10m

Value

dataframe of coordinates and associated depths

Author(s)

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

Examples

# load test data
data <- sf::read_sf(system.file("extdata", "example_contour.shp", package = 'rLakeHabitat'))
#run function
contourPoints(data, depths = "Z", geometry = "geometry", density = 50)

Cross Validate Interpolated Bathymetry

Description

Obtain residual mean square error (RMSE) from K-fold cross validation of bathymetry interpolation.

Usage

crossValidate(
  outline,
  df,
  x,
  y,
  z,
  zeros = FALSE,
  separation = NULL,
  k = 5,
  crsUnits = "dd",
  res = 5,
  method = "IDW",
  fact = NULL,
  nmax = 20,
  idp = 2,
  model = "Sph",
  psill = NULL,
  range = NULL,
  nugget = 0,
  kappa = NULL
)

Arguments

outline

shapefile outline of a waterbody

df

dataframe of coordinates and depths for a given waterbody

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, units from CRS

k

numeric value describing the number of folds to test, default = 5

crsUnits

character describing CRS units of input outline, either "dd" (decimal degrees) or "m" (meters), default = "dd"

res

number describing desired cell resolution in meters, default = 5

method

character describing method of interpolation, options include Inverse Distance Weighted ("IDW") or Ordinary Kriging ("OK"). Default = "IDW"

fact

numeric value describing the factor by which raster resolution should be increased, default = NULL. If 'crsUnits' and 'res' are defined, fact = NULL

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 Kriging, options include 'Sph', 'Exp', 'Gau', 'Sta', default = 'Sph'

psill

numeric value describing the partial sill value for OK interpolation, default = NULL

range

numeric describing distance beyond which there is no spatial correlation in Ordinary Kriging models, default = NULL

nugget

numeric describing variance at zero distance in Ordinary Kriging models, default = 0

kappa

numeric value describing model smoothness, default = NULL

Details

If both 'crsUnit' and 'res' = NULL, the output raster will be in the same CRS and units as the input 'outline' and the resolution will be increased by 'fact' (default = 10). If both 'crsUnit' and 'res' are defined, fact = NULL and the output raster will be projected to the most appropriate UTM zone at the specified resolution.

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 with distance. "Gau" = Gaussian model similar to spatial model but with stronger decay at shorter distances. "Mat" = Matern model Three parameters (psill, range, kappa) are incorporated from a fitted variogram (default = NULL). If specified in function input, chosen values will overwrite variogram values.

Value

mean RMSE value across k number of folds

Author(s)

Tristan Blechinger, 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
crossValidate(outline, data, "x", "y", "z", zeros = FALSE, separation = 10, k = 5, crsUnit = "dd",
res = 50, method = "IDW", nmax = 4, idp = 1.5)

Estimate Average Thermocline Depth

Description

Estimate average thermocline depth across multiple sites and dates.

Usage

estThermo(data, site, date, depth, temp, combine = "all")

Arguments

data

data frame of water column temperature profiles

site

character giving the name of the site column

date

character giving the name of the date column

depth

character giving the name of the depth column

temp

character giving the name of the temp column

combine

logical indicating whether or not to average across sites ("sites"), dates ("dates"), or sites and dates ("all"), default = "all"

Value

either numeric value of average thermocline depth, standard deviation, and n or data frame of thermocline depths, standard deviations, and n across sites or dates

Author(s)

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

Examples

# load test profile data
data <-  read.csv(system.file("extdata", "example_profile_data.csv", package = 'rLakeHabitat'))
data$date <- base::as.Date(data$date)
#run function
estThermo(data = data, site = "site", date = "date",
depth = "depth", temp = "temp", combine = "all")

Create Raster Stack

Description

Create a raster stack from a single raster, option to save as file.

Usage

genStack(
  DEM,
  by = 1,
  stop = NULL,
  save = TRUE,
  file_name = NULL,
  file_type = "COG"
)

Arguments

DEM

raster object

by

numeric increment per unit depth by which layers are split. Default = 1

stop

optional numeric value specifying depth at which to stop stacking rasters, default = NULL

save

logical, save raster stack (TRUE) or not (FALSE), default = TRUE

file_name

character string used to name saved raster stack

file_type

character string defining file type to save, default = "COG"

Value

a raster stack of specified depth increments for a given waterbody. Raster stack is either stored as an object (save = FALSE) or written to a file in the directory (save = TRUE).

Author(s)

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

Examples

#load raster
DEM <- terra::rast(system.file("extdata", "example_raster.tif", package = 'rLakeHabitat'))
#run function
genStack(DEM, by = 1, save = FALSE)

Interpolate bathymetry

Description

Generate a bathymetric digital elevation model (DEM) for a given waterbody using either Inverse Distance Weighting or Ordinary Kriging 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,
  crsUnits = "dd",
  res = 10,
  method = "IDW",
  fact = NULL,
  nmax = 20,
  idp = 2,
  model = "Sph",
  psill = NULL,
  range = NULL,
  nugget = 0,
  kappa = NULL
)

Arguments

outline

shapefile outline of a waterbody

df

dataframe of coordinates and depths for a given waterbody

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, units from CRS

crsUnits

character describing CRS units of input outline, either "dd" (decimal degrees) or "m" (meters), default = "dd"

res

number describing desired cell resolution in meters, default = 10

method

character describing method of interpolation, options include Inverse Distance Weighted ("IDW") or Ordinary Kriging ("OK"). Default = "IDW"

fact

numeric value describing the factor by which raster resolution should be increased, default = NULL If 'crsUnits' and 'res' are defined, fact = NULL

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 Kriging, options include 'Sph', 'Exp', 'Gau', 'Sta', default = 'Sph'

psill

numeric value describing the partial sill value for OK interpolation, default = NULL

range

numeric describing distance beyond which there is no spatial correlation in Ordinary Kriging models, default = NULL

nugget

numeric describing variance at zero distance in Ordinary Kriging models, default = 0

kappa

numeric value describing model smoothness, default = NULL

Details

If 'res' and 'crsUnits' are specified (recommended), the output raster is returned in the original projection at the specified resolution. If 'res' and 'crsUnits' are not specified, 'fact' must be defined as a numeric value by which the resolution of the output DEM will be increased (1 = no change). Output raster will be returned in the original projection of the input.

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 with distance. "Gau" = Gaussian model similar to spatial model but with stronger decay at shorter distances. "Mat" = Matern model Three parameters (psill, range, kappa) are incorporated from a fitted variogram (default = NULL). If specified in function input, chosen values will overwrite variogram values.

DEMs generated using OK method will have two layers; the first are the interpolated values and the second are the variances associated with each measurement

Value

DEM of waterbody bathymetry

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,
crsUnit = "dd", res = 5, method = "IDW", nmax = 4, idp = 2)

Calculate Littoral Volume

Description

Calculate littoral and pelagic volume across water levels from a DEM based on estimated photic depth.

Usage

littoralVol(
  DEM,
  photic,
  secchi = NULL,
  DEMunits = "m",
  depthUnits = "ft",
  by = 1
)

Arguments

DEM

SpatRaster object of a given waterbody, rasters can be transformed to SpatRaster via the rast() function in 'terra'

photic

number giving the average photic depth, overwrites Secchi depth

secchi

number giving the average secchi depth, photic zone estimated as 2.6m * secchi

DEMunits

character describing units of raster coordinate system. Can be meters, kilometers, or hectares ("m", "km", "ha"), default = "m"

depthUnits

character describing units of depth measurement. Can be either feet or meters ("ft", "m"), default = "ft"

by

numeric increment per unit by which volumes are calculated. Higher values will result in lower resolution. Default = 1

Value

data frame of littoral and pelagic volume estimates

Author(s)

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

Examples

#load raster
DEM <- terra::rast(system.file("extdata", "example_raster.tif", package = 'rLakeHabitat'))
#run function
littoralVol(DEM, photic = 2, DEMunits = "m", depthUnits = "m", by = 1)

Rarify Depth Data

Description

Reduce density of mapped depth data to improve accuracy and computation time.

Usage

rarify(outline, df, x, y, z, res = 10, crsUnits = "dd")

Arguments

outline

shapefile outline of a waterbody

df

dataframe of coordinates and depths for a given waterbody

x

character giving name of longitude column

y

character giving name of latitude column

z

character giving name of depth column

res

number describing by how much to increase point resolution, default = 10

crsUnits

character describing CRS units of input outline, either "dd" (decimal degrees) or "m" (meters), default = "dd"

Value

dataframe of rarified xyz coordinates

Author(s)

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

Examples

#load test data
outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
depths <- read.csv(system.file("extdata", "example_depths.csv", package = 'rLakeHabitat'))
#run function
rarify(outline = outline, df = depths, x = "x", y = "y", z = "z", res = 100, crsUnits = "dd")