The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Introduction to rLakeHabitat

Tristan Blechinger & Sean Bertalot

March 14, 2025

rLakeHabitat Introduction and Guide

This user manual introduces the functions of the rLakeHabitat package and provides an example workflow for generating digital elevation models (DEMs) of waterbody bathymetry from multiple types of depth data. The goal of this package is to provide managers and researchers with a tool to quantify area and volume of physical aquatic habitats based on estimates of photic depth and thermocline temperature (e.g., littoral zone, epilimnion, and hypolimnion). Functions included in this package facilitate obtaining and formatting depth data (rarify, contourPoints), estimating average thermocline depths (estThermo), interpolating depths (interpBathy, crossValidate), quanitfying physical habitats (calcHyps, calcLittoral, calcSDI, calcVolume, littoralVol), and visualizing outputs (animBathy, bathyMap). Creation of this package relies on functions from ggplot, terra, sf, and gstat.

Step 1: Loading Data

Depth data can either be in raw ‘xyz’ format or generated from contour shapefiles with known depths. Coordinates should be in decimal degrees. Depths can be in either meters or feet. If measured depth data has a high density of points, the rarify function can be used to reduce point density by averaging depths within raster cells of a specified resolution. A shapefile of the waterbody outline is required to bound rarification. Examples are included below.

library(rLakeHabitat)

#load example depth data
depths <- read.csv("data/example_depths.csv")
#load example outline
outline <- vect("data/example_outline.shp")

#rarify input data -- example data won't change
clean_depth <- rarify(outline, depths, "x", "y", "z", res = 5)

#obtain xyz point data from depth contours
contours <- read_sf("~/UW_Research/Project5_Habitat_Mapping/Repositories/rLakeHabitat/data/example_contour.shp")

contour_depths <- contourPoints(contours, depths = "Z", geometry = "geometry", density = 50)

Step 2: Interpolation

Inverse Distance Weighted (“IDW”) and Ordinary Kriging (“OK”) interpolation methods are supported in the interpBathy function. Depth data and a shapefile outline to bound interpolation are required. To interpolate at a desired resolution, use the crsUnits and res parameters to define a cell size in meters. If a specified resolution is not desired, the fact parameter will increase output resolution by the specified factor. Interpolated DEMs will be in the same units as input depths and the same CRS as the input outline. DEMs generated using the OK method will have two layers, the interpolated values (1) and variance (2). The crossValidate function uses the interpBathy function on k number of folds to compare subsampled test and train datasets to quantify RMSE of the interpolated DEM. Prior to subsampling, points are stratified by depth and evenly sampled across five depth groups to ensure accurate representation of data. Examples are included below.

#load example depth data
depths <- read.csv("data/example_depths.csv")
#load example outline
outline <- vect("data/example_outline.shp")

#interpolate using Inverse Distance Weighted method
DEM <- interpBathy(outline, depths, "x", "y", "z", zeros = F, separation = 10, 
                   crsUnits = "dd", res = 10, method = "IDW", nmax = 6, idp = 2)
plot(DEM)

#obtain RMSE for IDW interpolation
crossValidate(outline, depths, "x", "y", "z", zeros = F, separation = 10, k = 5, 
              crsUnits = "dd", res = 10, method = "IDW", nmax = 6, idp = 2)

#Interpolate using Ordinary Kriging method
DEM <- interpBathy(outline, depths, "x", "y", "z", zeros = F, separation = 10, 
                   crsUnits = "dd", res = 50, nmax = 6, method = "OK", model = "Sph")
plot(DEM)

#obtain RMSE for OK interpolation
crossValidate(outline, depths, "x", "y", "z", zeros = F, separation = 10, k = 5, 
              crsUnits = "dd", res = 50, nmax = 6, method = "OK", model = "Sph")

Step 3: Habitat Quantification

Physical aquatic habitats from interpolated DEMs can be quantified using multiple functions. Habitats include littoral area, shoreline development index, and epilimnetic, hypolimnetic, metalimnetic, and littoral volume as well as hypsography data. Most functions require an average Secchi or photic depth (photic depth = Secchi x 2.6 when not specified). Average thermocline depth can be averaged across sites and/or dates using the estThermo function, which uses the thermo.depth function from rLakeAnalyzer. Examples are included below.

#load example temperature profile data
thermo_data <- read.csv("data/example_profile_data.csv") %>% 
  mutate(date = as.Date(date))

#estimate average thermocline depth across sites and dates
estThermo(thermo_data, site = "site", date = "date", depth = "depth", temp = "temp", combine = "all")

#generate hypsography data, output = 'values' or 'plot'
calcHyps(DEM, DEMunits = 'm', depthUnits = 'ft', by = 1, output = 'values')

#calculate littoral area
calcLittoral(DEM, secchi = 2, DEMunits = 'm', depthUnits = 'ft', by = 1)

#calculate shoreline development index
calcSDI(DEM, units = 'm', by = 1)

#calculate volume of pelagic habitats
calcVolume(DEM, thermo_depth = 4, DEMunits = 'm', depthUnits = 'm', by = 1)

#calculate volume of pelagic vs. littoral habitat
littoralVol(DEM, secchi = 2, DEMunits = 'm', depthUnits = 'm', by = 1)

Step 4: Visualization

There are two functions to visualize interpolated DEMs. The bathyMap function generates a ggplot object with the option to include contours and contour labels as well as the spacing between plotted contours. The animBathy function is an animated ggplot object meant to visualize changes in littoral or whole surface area of a waterbody across water levels. Examples are included below.

#generate bathymetry map
bathyMap(DEM, contours = T, units = "m", by = 1)

#generate animated bathymetry map of littoral area or whole waterbody area across water levels
animBathy(DEM, units = 'm', littoral = T, secchi = 2, by = 1)

Save DEM

The genStack function is a convenient way to generate raster stack objects with the option to save them as a file. When saving, the default filetype is cloud optimized geotiff (.COG), though this can be specified to other gdal accepted filetypes. Examples are included below.

#generate raster stack from interpolated DEM
raster_stack <- genStack(DEM, by = 1, save = F) #don't save
plot(raster_stack)

genStack(DEM, by = 1, save = T, file_name = "Example_stack") #save

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.