Introduction to rLakeHabitat

Overview

rLakeHabitat provides a complete workflow for turning raw depth measurements into a bathymetric digital elevation model (DEM), quantifying physical habitat metrics from it, and visualizing the results. This vignette walks through that workflow, using an example waterbody to illustrate the package’s functions and the reasoning behind them.

Setup

Beyond rLakeHabitat itself, this vignette uses a few companion packages for spatial data handling, plotting, and site selection:

library(rLakeHabitat)
library(terra)
#> Warning: package 'terra' was built under R version 4.4.3
library(tidyterra)
#> Warning: package 'tidyterra' was built under R version 4.4.3
library(ggspatial)
#> Warning: package 'ggspatial' was built under R version 4.4.3
library(sf)
#> Warning: package 'sf' was built under R version 4.4.3
library(dplyr)
library(httr)
library(spbal)
#> Warning: package 'spbal' was built under R version 4.4.3
library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.4.3
library(mapview)
#> Warning: package 'mapview' was built under R version 4.4.3

extdata_path <- function(file) {
  system.file("extdata", file, package = "rLakeHabitat")
}

Part 1: Study Design

Every bathymetric survey starts with the same question: how far apart can transects be spaced and still produce an accurate map? Most waterbodies are too large to efficiently and accurately map the entire surface, so depth data is generally collected along transects spaced at some fixed distance. That spacing matters: too narrow and you’re spending time and effort collecting more data than you need; too large and you risk missing real features in the shape of the waterbody.

The samplingDensity() function addresses this before data collection by simulating a plausible bathymetry for a lake (a smooth depth profile plus realistic random roughness), testing several candidate transect spacings against it, and reporting how reconstruction accuracy degrades as spacing increases. It does this using sequential Gaussian simulation (SGS) to generate a “true” bathymetry. SGS visits each cell in the outline and uses variogram parameters to estimate a distribution of potential values based on a weighted average of each cell’s neighbors, then draws a random depth from that distribution. The variogram range controls how many surrounding cells influence that estimate: a small range produces a choppier surface (less influence from neighbors), while a larger range produces a smoother one. The sill controls the variance of each cell’s distribution: a larger sill allows a wider range of possible depths at each point, producing a bumpier surface, while a smaller sill produces a flatter one.

samplingDensity() also identifies a recommended spacing: once RMSE is calculated for each candidate spacing, the largest spacing within a specified tolerance (default 10%) of the minimum RMSE is selected. This reflects a real tradeoff between sampling cost and accuracy - for example, increasing transect spacing from 50 to 100 m might cut the total number of transects in half while increasing RMSE by less than 10%.

The function requires a best estimate of maximum depth, and optionally a shape describing whether the basin is a steep bowl (shape < 1) or broad and flat (shape > 1). Below, we simulate a plausible bathymetry for Alcova Reservoir (max depth 55 m, a moderately steep basin) and evaluate three candidate transect spacings against it.

# Load Alcova outline shapefile
outline <- read_sf(extdata_path("example_outline.shp"))

crs(outline) # verify input CRS
#> [1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"

design <- samplingDensity(outline, # shapefile outline
                          max_depth = 6, 
                          shape = 1, # value describing basin shape - steep (shape < 1) or shallow (shape > 1)
                          truth_range = 100, # distance where correlation ends - large range = smooth ridges/valleys, small range = jagged terrain
                          truth_sill = 5, # total variance of field, how big do bumps get - low = flat, high = variable
                          along_track_interval = 5, # distance between points sampled along transect
                          spacings = c(50, 100), # hypothetical distances between transect lines (meters)
                          res = 50, # value describing cell size in meters, lower = more cells and vice versa
                          n_sim = 3, # number of simulations to run - more is better, at the cost of runtime
                          plot = T, 
                          seed = 123)
#>   |                                                                              |                                                                      |   0%[using unconditional Gaussian simulation]
#>   |                                                                              |============                                                          |  17%  |                                                                              |=======================                                               |  33%[using unconditional Gaussian simulation]
#>   |                                                                              |===================================                                   |  50%  |                                                                              |===============================================                       |  67%[using unconditional Gaussian simulation]
#>   |                                                                              |==========================================================            |  83%  |                                                                              |======================================================================| 100%


plot(design$truth_examples) # simulated "true" lake bathymetry


# accuracy vs. transect spacing
plot(design$results$spacing, design$results$mean_rmse, type = "b", pch = 19,
     xlab = "Transect spacing (m)", ylab = "RMSE (m)",
     main = "Simulated DEM accuracy vs. transect spacing")


# recommended spacing, and the total transect length it implies
design$recommended_spacing
#> [1] 50
design$total_transect_length
#> [1] 215294.3

Once a spacing is chosen, the resulting transects can be saved to a .gpx file for use in most common GPS units:

# reproject to WGS84 for .gpx format
transects_wgs84 <- terra::project(design$transects, "EPSG:4326")

transects_sf <- sf::st_as_sf(transects_wgs84)
transects_sf <- sf::st_cast(transects_sf, "MULTILINESTRING")
transects_sf <- sf::st_cast(transects_sf, "LINESTRING")

transects_sf$name <- paste0("Transect_", seq_len(nrow(transects_sf)))
transects_sf <- transects_sf["name"]

sf::st_write(transects_sf, "transects.gpx", driver = "GPX",
             layer_options = "FORCE_GPX_ROUTE=YES", delete_dsn = TRUE)

# Or as a .kml:
# sf::st_write(transects_sf, "transects.kml", driver = "KML")

Part 2: Data Processing

With a transect spacing chosen and sampling complete, the next step is turning raw depth data into a clean interpolation input. Depth data typically arrives one of two ways: raw XYZ points from a sonar unit, or depths encoded as attributes on contour polygons digitized from an existing chart. The general workflow is:

Load data
Convert contours to points (if necessary)
Clean data: remove outlying depths (e.g., 0 or 1000)
Rarify points

Depth data collected using boat-mounted sonar is often very dense - measurements are taken continuously, resulting in multiple points over very small spatial scales. Counterintuitively, this can actually reduce interpolation accuracy by biasing estimates toward data-dense regions. rarify() addresses this by calculating the mean of all depth measurements within each defined cell area. Boat-mounted sonar can also record erroneous depths when the transducer isn’t oriented properly (e.g., while turning), so it’s worth reviewing the data to catch false “islands” before interpolating.

One important consideration with spatial data is the coordinate reference system (CRS), which describes how coordinates are measured and projected. rLakeHabitat’s interpolation and analysis functions handle CRS reprojection internally for their own working calculations, but they assume the depth data and outline are already in the same CRS as one another and do not verify this. If they aren’t, results can be wrong (points may be spatially shifted but still fall within the outline, producing a plausible-looking but incorrect DEM). Always explicitly check and reconcile CRSs before interpolating.

# load contours and sample points
contours <- read_sf(extdata_path("example_contour.shp"))

head(contours)
#> Simple feature collection with 6 features and 5 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -104.7057 ymin: 42.16475 xmax: -104.6948 ymax: 42.16905
#> Geodetic CRS:  WGS 84
#> # A tibble: 6 × 6
#>       Z Field Field1 Shape_Leng   Shape_Area                            geometry
#>   <dbl> <dbl>  <dbl>      <dbl>        <dbl>                  <MULTIPOLYGON [°]>
#> 1    40     0      0   0.00490  0.00000136   (((-104.7003 42.16636, -104.7002 4…
#> 2    40     0      0   0.00552  0.00000203   (((-104.7034 42.16905, -104.7033 4…
#> 3    35     0      0   0.00585  0.00000115   (((-104.6959 42.16705, -104.6958 4…
#> 4    35     0      0   0.00341  0.000000670  (((-104.7044 42.16679, -104.7043 4…
#> 5    35     0      0   0.00126  0.0000000976 (((-104.7056 42.16596, -104.7056 4…
#> 6    35     0      0   0.000943 0.0000000534 (((-104.7045 42.16496, -104.7045 4…

ggplot()+
  geom_sf(data = contours, aes(color = Z))


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


# were going to continue from here with our raw point data, but could replace that with the 'depths' dataframe

depths <- read.csv(extdata_path("example_depths.csv"))

head(depths)
#>           x        y z
#> 1 -89.26792 42.97233 1
#> 2 -89.26081 42.97508 1
#> 3 -89.25180 42.97392 2
#> 4 -89.24628 42.96999 1
#> 5 -89.26208 42.96755 5
#> 6 -89.25275 42.96543 7

crs(outline) #EPSG:4326
#> [1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"
crs(depths) #No CRS, not projected yet -- recorded in WGS84 (EPSG:4326) -- good
#> [1] NA

# # if needed, we could reproject the outline to match the depth data
# outline <- st_transform(AL_outline, crs = 4326)
# crs(outline)

# Before we rarify, we could verify our depths are accurate
outliers <- depths %>%
  filter(z < 1 | z > 10) %>%
  st_as_sf(coords = c("x", "y"), crs = 4326)

mapview::mapview(outline) +
 outliers


# # optional code to remove outliers
# clean <- depths %>%
#   filter(!between(x, -108.17588, -108.17548)) %>%
#   filter(!between(y, 43.38636, 43.38667))

# can manually add maximum and minimum depth points if necessary by getting coordinates from mapview() and entering them in the dataframe

# Now that we're comfortable with the depth data, we can rarify
rare <- rarify(outline, depths, "x", "y", "z", res = 10) # the example data is already sparse, so it won't change

Part 3: Interpolation

Interpolation is the process that fills in all the unknown gaps between our measured depths. The interpBathy() function can interpolate using three methods, Inverse Distance Weighting (IDW), Ordinary Kriging (OK), or Universal Kriging (UK). While both accomplish the same thing, they function in very different ways. IDW is a deterministic weighting rule based purely on distance. Each unknown point is calculated using the weighted average of neighboring points (nmax), where the weight is determined by the inverse distance power (idp). By increasing or decreasing the inverse distance power, you are tuning the model to weight neighboring points stronger or weaker. The relationship is weight = 1/distance^idp. Because IDW is deterministic, there is no uncertainty in the estimates. This method is computationally efficient and works well on both spatially well distributed datasets and sparse datasets where OK may be unable to fit a variogram.

Conversely, kriging is a geostatistical model that learns the spatial correlation structure of the data before predicting unknown points. It does this by fitting a variogram - a model of how much measurements differ as a function of the distance between them. The variogram is calculated by finding the distance and squared difference in depth value for each pair of points. The variogram tells the interpolation how much influence neighbors have at each distance in the actual dataset, so the model can weight points individually based on their distance from the unknown value. The variogram can be tuned to use different models (curve shapes) and specific range, sill, and nugget values. Range is the distance at which a point becomes noninfluential, sill (or partial sill) is the variance at the range, and the nugget is the variance at a distance of zero (the y-intercept). Here, the nugget represents measurement error. By default, interpBathy() automatically fits a nugget from the data’s empirical variogram, the same way it fits sill and range - a fixed value can be supplied instead (e.g., nugget = 0, to force the interpolation to pass exactly through each known depth) if preferred.

Each variogram model suggests a different relationship between variance and distance. Spherical models are often the default and tend to best represent moderate surface smoothness. Exponential models have the steepest behavior near the origin meaning influence of distant points decreases faster, but it asymptotes to the sill so even distant points can have a very small influence. Gaussian models represent the smoothest surface of the three with shallow behavior around the origin of the variogram that increases the influence of distant points. The interpBathy() function defaults to a spherical variogram model and automatically calculates the variogram parameters before using them in the interpolation, but they can be overridden if desired. Kriging is a more computationally intensive method than IDW but can be more accurate for many datasets and provides a variance layer alongside the estimates.

OK and UK differ in how they handle the mean of the unknown depths. OK uses a constant mean, assuming that each unknown depth is the same and spatial patterns in depth are caused only by the spatial correlation between nearby points. UK treats the mean as a trend, so the expected depth differs based on the location of the unknown point. OK is computationally easier than UK because it lacks the extra calclation for the floating mean. UK can be parameterized with two types of trends for the mean - linear or quadratic. A linear trend suggests that mean depth is a flat, tilted plane across the lake while a quadratic trend suggests mean depth follows a parabolic bowl shape across the lake.

Interpolations can be simple or complex, depending on the structure of the data and the desired result. Below, we run some example interpolations using IDW and OK.

# Basic IDW interpolation
IDW <- interpBathy(outline, rare, "x", "y", "z", separation = 10, res = 50, nmax = 8)
#>   |                                                                              |                                                                      |   0%
#>   |                                                                              |==============                                                        |  20%  |                                                                              |============================                                          |  40%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
#> [inverse distance weighted interpolation]
#> |---------|---------|---------|---------|[inverse distance weighted interpolation]
#> =========================================                                            |                                                                              |==========================================                            |  60%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
#> |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          [inverse distance weighted interpolation]
#> |---------|---------|---------|---------|[inverse distance weighted interpolation]
#> =========================================                                          |---------|---------|---------|---------|=========================================                                            |                                                                              |========================================================              |  80%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                            |                                                                              |======================================================================| 100%[1] "Raster Interpolated in CRS:  GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"
plot(IDW, main = "IDW-interpolated DEM")


# Basic OK interpolation
OK <- interpBathy(outline, rare, "x", "y", "z", separation = 20, res = 50, method = "OK", nmax = 8)
#>   |                                                                              |                                                                      |   0%
#>   |                                                                              |==============                                                        |  20%  |                                                                              |============================                                          |  40%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
#> Warning in gstat::fit.variogram(vgram, model = init_vgm): No convergence after
#> 200 iterations: try different initial values?
#> [using ordinary kriging]
#> |---------|---------|---------|---------|[using ordinary kriging]
#> =========================================                                            |                                                                              |==========================================                            |  60%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
#> |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          [using ordinary kriging]
#> |---------|---------|---------|---------|[using ordinary kriging]
#> =========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                            |                                                                              |========================================================              |  80%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                            |                                                                              |======================================================================| 100%[1] "Raster Interpolated in CRS:  GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"
plot(OK, main = "OK-interpolated DEM (depth + error layers)")

You may notice that the OK interpolation above raised a warning: "No convergence after 200 iterations". This means the NLS minimizer couldn’t reach stable values for the nugget, psill, and range parameters while fitting the Gaussian model to the empirical variogram after 200 attempts - a common outcome with low sample density, or when the true variogram shape doesn’t match the model being fit. In this situation, running a variogram manually, comparing model shapes, and identifying the best-fitting parameters directly is a useful diagnostic step. The code below mirrors interpBathy()’s own internal CRS handling by reprojecting to a projected CRS before fitting, allowing model shapes to be compared directly.

if (!inherits(outline, "SpatVector")) {
  outline <- terra::vect(outline)
}

if (terra::is.lonlat(outline)) {
  centroid <- terra::centroids(outline)
  lon <- terra::crds(centroid)[1]; lat <- terra::crds(centroid)[2]
  utm_zone <- floor((lon + 180) / 6) + 1

  epsg <- ifelse(lat >= 0, 32600, 32700) + utm_zone
  best_crs <- terra::crs(paste0("EPSG:", epsg))

  pts <- terra::vect(rare, geom = c("x", "y"), crs = terra::crs(outline))
  pts_proj <- terra::project(pts, best_crs)
  coords <- terra::crds(pts_proj)
  rare$x_m <- coords[, 1]
  rare$y_m <- coords[, 2]
} else {
  rare$x_m <- rare$x
  rare$y_m <- rare$y
}

# calculate empirical points on the variogram and plot
# This will look bad given the sparse data
emp_vgm <- gstat::variogram(z ~ 1, locations = ~x_m + y_m, data = rare)
plot(emp_vgm)


# fit a curve to the variogram
fit_vgm <- gstat::fit.variogram(emp_vgm, model = gstat::vgm(model = "Gau"))
fit_vgm
#>   model    psill    range
#> 1   Gau 11.69316 978.2081
plot(emp_vgm, fit_vgm)


# compare different curves and see which model type fits best
for (m in c("Sph", "Exp", "Gau")) {
  f <- gstat::fit.variogram(emp_vgm, model = gstat::vgm(model = m))
  cat(m, "- SSErr:", attr(f, "SSErr"), "\n")
}
#> Warning in gstat::fit.variogram(emp_vgm, model = gstat::vgm(model = m)):
#> singular model in variogram fit
#> Sph - SSErr: 0.0006850851
#> Warning in gstat::fit.variogram(emp_vgm, model = gstat::vgm(model = m)): No
#> convergence after 200 iterations: try different initial values?
#> Exp - SSErr: 0.0003537983 
#> Gau - SSErr: 0.0003282291

Check Interpolation Accuracy

Now that DEMs have been generated, the next step is checking how accurate they are. While OK and UK interpolations return a raster of variances for each cell, IDW does not. Interpolations can be cross-validated with the crossValidate() function to calculate residual mean square error (RMSE), which is useful both to gauge accuracy and to compare parameters across interpolations to find the best fit. crossValidate() runs k-fold cross-validation with folds blocked spatially and by depth, to ensure the training data still has good coverage of the system. This is more computationally intensive than a single interpolation, since the function splits depth data into k folds and runs k interpolations, leaving one fold out each time to test on. The example below uses IDW interpolation, a low resolution, a low number of folds, and a low nmax to keep runtime short; in practice, at least 5 folds is recommended. This is an excellent method to compare different interpolation parameters to identify the most accurate model.

# calculate RMSE
crossValidate(outline, rare, "x", "y", "z", zeros = FALSE, separation = 20,
              k = 3, res = 10, method = "IDW", nmax = 8, zero_threshold = .5, seed = 123)
#>   |                                                                              |                                                                      |   0%
#>   |                                                                              |==============                                                        |  20%  |                                                                              |============================                                          |  40%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
#> [inverse distance weighted interpolation]
#> |---------|---------|---------|---------|[inverse distance weighted interpolation]
#> =========================================                                            |                                                                              |==========================================                            |  60%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
#>   |                                                                              |========================================================              |  80%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                            |                                                                              |======================================================================| 100%[1] "Raster Interpolated in CRS:  GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"
#> 
#>   |                                                                              |                                                                      |   0%
#>   |                                                                              |==============                                                        |  20%  |                                                                              |============================                                          |  40%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
#> [inverse distance weighted interpolation]
#> |---------|---------|---------|---------|[inverse distance weighted interpolation]
#> =========================================                                            |                                                                              |==========================================                            |  60%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
#>   |                                                                              |========================================================              |  80%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                            |                                                                              |======================================================================| 100%[1] "Raster Interpolated in CRS:  GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"
#> 
#>   |                                                                              |                                                                      |   0%
#>   |                                                                              |==============                                                        |  20%  |                                                                              |============================                                          |  40%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
#> [inverse distance weighted interpolation]
#> |---------|---------|---------|---------|[inverse distance weighted interpolation]
#> =========================================                                            |                                                                              |==========================================                            |  60%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
#>   |                                                                              |========================================================              |  80%|---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                            |                                                                              |======================================================================| 100%[1] "Raster Interpolated in CRS:  GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"
#>     RMSE 
#> 1.917906

Tuning Point Density and Resolution

Two major decisions in the interpolation process are how much to rarify your input data and what resolution to interpolate your DEM at. While you could cross validate different combinations of interpolations, the optimizeParams function does that automatically. Enter a list of rarification densities (spacings), resolutions (res_values), and the parameters for your interpolation and compare RMSEs across those scenarios. This function runs both the rarify() and crossValidate() functions, so this is the same result as cross validating different combinations manually.

tuning <- optimizeParams(outline, depths, "x", "y", "z",
                          spacings = c(50, 100), res_values = c(10, 50),
                          k = 3, zeros = FALSE, separation = 20,
                          nmax = 4, plot = T, seed = 123, zero_threshold = .6)
#>   |                                                                              |                                                                      |   0%
#>   |                                                                              |==================                                                    |  25%
#>   |                                                                              |===================================                                   |  50%
#>   |                                                                              |====================================================                  |  75%
#>   |                                                                              |======================================================================| 100%


tuning$results 
#>   spacing res mean_rmse n_points
#> 1      50  10  1.738981       14
#> 2     100  10  1.436033       14
#> 3      50  50  1.569734       14
#> 4     100  50  1.315579       14

# n_points is the number of points after rarification (one value for each 'spacings' value)
# this allows you to see how many data points are being used in the interpolation

tuning$recommended_spacing
#> [1] 100
tuning$recommended_res
#> [1] 50

Part 4: Quantifying Physical Habitat

With DEMs generated, the next step is quantifying habitat within each waterbody.

Habitat Boundaries

Two important physical aquatic habitat boundaries are considered here: thermocline depth and photic depth.

Thermocline depth regulates the volume of the epilimnion, metalimnion, and hypolimnion. These distinct habitats influence internal nutrient dynamics and are important for aquatic organisms. Given temperature profile data, estThermo() estimates average thermocline depth across sites and/or dates (built on rLakeAnalyzer):

# read in profile data
profiles <- read.csv(extdata_path("example_profile_data.csv")) %>%
  mutate(date = as.Date(date))


# plot profiles
ggplot(profiles, aes(x = temp, y = depth, color = site)) +
  geom_path() +
  geom_point(size = 0.8) +
  scale_y_reverse() +
  facet_wrap(~date) +
  labs(x = "Temperature (°C)", y = "Depth (m)", color = "Site") +
  theme_minimal()



# calculate average thermocline depth across all sites and dates
# change 'combine =' to get values for specific sites or dates
estThermo(profiles, site = "site", date = "date", depth = "depth", temp = "temp", combine = "all")
#>      mean       sd  n
#> 1 10.5195 2.770556 24

Several habitat functions want a photic depth. Photic depth functions as the boundary of the littoral habitat zone in many waterbodies, an important area for nutrient cycling and aquatic organisms where benthic algae and macrophyte growth is possible. Photic depth is most accurate if measured directly with a PAR sensor, but it can be estimated from other measures of water clarity like a Secchi disk. calcPhotic() converts Secchi disk measurements to photic depth using the light-attenuation relationship from Kirk (1994) and Koenings & Edmundson (1991) where Z is the Secchi disk depth and F is a numerical constant that represents optical properties of the water. Standard values for F are 2.76 for stained lakes, 1.99 for clear lakes, and 1.05 for turbid lakes.

photic <- calcPhotic(Z = 1, F = 1.99)

Littoral area, volume, and shoreline complexity

With thermocline and photic depths approximated, the area of littoral habitats and volumes of stratified habitats can be calculated. A hypsographic curve, showing total surface area plotted against depth, is a useful starting point for understanding the general shape of each system.

calcLittoral() estimates the area of littoral and non-littoral habitats across water levels, given a photic depth (the littoral boundary) and a depth interval for the estimates. These functions were designed for reservoir ecosystems, which often exhibit substantial water level fluctuation within and across years - hence quantifying area across a range of depths rather than just at the surface. If only the littoral extent at the current surface level is of interest, only the first row of the output table is needed. littoralVol() integrates across depth slices to estimate littoral and non-littoral (pelagic) volume at each water level.

calcSDI() calculates the Shoreline Development Index across water levels, a useful metric of nearshore habitat complexity.

# hypsography: surface area at each depth increment
# change `output` to "values" for a table of areas by depths
calcHyps(IDW, depthUnits = "m", by = 1, output = "plot")



# littoral (photic-zone) surface area across water levels
lit <- calcLittoral(IDW, photic = photic, depthUnits = "m", by = 1)
lit
#>   depth   tot_area  lit_area perc_lit_tot
#> 1     0 11103145.3 8637865.7     77.79657
#> 2     1  4148466.5 2476081.8     59.68668
#> 3     2  2792721.7 1796866.9     64.34107
#> 4     3  1894285.4 1602022.2     84.57132
#> 5     4  1193399.6 1128453.6     94.55791
#> 6     5   462751.4  462751.4    100.00000
#> 7     6   102831.0  102831.0    100.00000


# shoreline development index (how convoluted the shoreline is relative to a circle)
sdi <- calcSDI(IDW, units = "m", by = 1)
sdi
#>   depth   tot_area perimeter      SDI
#> 1     0 11016553.5 15130.400 1.285946
#> 2     1  4148466.5  9906.101 1.372000
#> 3     2  2792721.7  8221.590 1.387833
#> 4     3  1894285.4  6747.621 1.383003
#> 5     4  1193399.6  5402.620 1.395105
#> 6     5   462751.4  4203.765 1.743251
#> 7     6   102831.0  1441.700 1.268259


# littoral vs. pelagic volume
littoralVol(IDW, photic = photic, depthUnits = "m", by = 1)
#>   depth tot_vol_m3 lit_vol_m3 pel_vol_m3 perc_lit_tot perc_pel_tot tot_vol_acft
#> 1     0 21697600.8 16206872.7 5490728.13        74.69        25.31  17590.53141
#> 2     1 10594455.5  7569007.0 3025448.55        71.44        28.56   8589.06494
#> 3     2  6445989.1  5092925.1 1353063.93        79.01        20.99   5225.84843
#> 4     3  3653267.4  3296058.3  357209.13        90.22         9.78   2961.75210
#> 5     4  1758982.0  1694036.1   64945.92        96.31         3.69   1426.02994
#> 6     5   565582.5   565582.5       0.00       100.00         0.00    458.52517
#> 7     6   102831.0   102831.0       0.00       100.00         0.00     83.36647
#>   lit_vol_acft pel_vol_acft
#> 1  13139.12564   4451.40577
#> 2   6136.29387   2452.77108
#> 3   4128.90165   1096.94679
#> 4   2672.15795    289.59416
#> 5   1373.37743     52.65251
#> 6    458.52517      0.00000
#> 7     83.36647      0.00000

Base terra functions can also be used directly for metrics like average slope or true (terrain-corrected) surface area.

# Before calculating these metrics, the raster needs to be reprojected to a
# projected CRS (here, the appropriate UTM zone), since slope and surface area
# calculations require real linear distance units.

# quick way to get the right UTM zone from the raster's centroid
get_best_utm <- function(r) {
  centroid <- terra::centroids(terra::as.polygons(terra::ext(r), crs = terra::crs(r)))
  centroid <- terra::project(centroid, "EPSG:4326")
  lon <- terra::crds(centroid)[1]
  lat <- terra::crds(centroid)[2]
  utm_zone <- floor((lon + 180) / 6) + 1
  hemisphere <- ifelse(lat >= 0, 32600, 32700)
  paste0("EPSG:", hemisphere + utm_zone)
}

utm_crs <- get_best_utm(IDW)
IDW_utm <- terra::project(IDW[[1]], utm_crs)

# Average slope
slope_deg <- terra::terrain(IDW_utm, v = "slope", unit = "degrees")
plot(slope_deg)


avg_slope <- terra::global(slope_deg, "mean", na.rm = TRUE)
avg_slope
#>            mean
#> slope 0.1820658

# Total true (terrain-corrected) surface area
surf_area_cells <- terra::surfArea(IDW_utm)
total_true_area <- terra::global(surf_area_cells, "sum", na.rm = TRUE)
total_true_area   # in sq meters
#>           sum
#> lyr1 11098516

Thermally defined pelagic habitat volumes

The calcVolume() function estimates the volume of the epilimnion, metalimnion, and hypolimnion by integrating area across depth. It takes either a single thermocline midpoint (to estimate epilimnion and hypolimnion volumes) or an upper and lower thermocline boundary (to additionally estimate metalimnion volume). Using the thermocline depths calculated earlier - 12.45 m for Alcova and 7.18 m for Bull Lake:

vol <- calcVolume(IDW, thermo_depth = 3, depthUnits = "m", by = 1)
vol
#>   depth tot_vol_m3 epi_vol_m3 hyp_vol_m3 perc_epi_tot perc_hyp_tot tot_vol_acft
#> 1     0 21697600.8   19938619  1758982.0        91.89         8.11  17590.53140
#> 2     1 10594455.5   10028873   565582.5        94.66         5.34   8589.06493
#> 3     2  6445989.1    6343158   102831.0        98.40         1.60   5225.84843
#> 4     3  3653267.4    3653267        0.0       100.00         0.00   2961.75209
#> 5     4  1758982.0         NA         NA           NA           NA   1426.02993
#> 6     5   565582.5         NA         NA           NA           NA    458.52517
#> 7     6   102831.0         NA         NA           NA           NA     83.36647
#>   epi_vol_acft hyp_vol_acft
#> 1    16164.501   1426.02993
#> 2     8130.540    458.52517
#> 3     5142.482     83.36647
#> 4     2961.752      0.00000
#> 5           NA           NA
#> 6           NA           NA
#> 7           NA           NA


# Supplying thermo_high/thermo_low instead of a single thermo_depth additionally
# estimates metalimnion volume between the two boundaries
calcVolume(IDW, thermo_high = 2, thermo_low = 4, depthUnits = "m", by = 1)
#>   depth tot_vol_m3 epi_vol_m3 met_vol_m3 hyp_vol_m3 perc_epi_tot perc_met_tot
#> 1     0 21697600.8   18044333    3087685   565582.5        83.16        14.23
#> 2     1 10594455.5    8835474    1656151   102831.0        83.40        15.63
#> 3     2  6445989.1    5880407     565582        0.0        91.23         8.77
#> 4     3  3653267.4    3550436         NA         NA        97.19           NA
#> 5     4  1758982.0    1758982         NA         NA       100.00           NA
#> 6     5   565582.5         NA         NA         NA           NA           NA
#> 7     6   102831.0         NA         NA         NA           NA           NA
#>   perc_hyp_tot tot_vol_acft epi_vol_acft met_vol_acft hyp_vol_acft
#> 1         2.61  17590.53140    14628.779    2503.2270    458.52517
#> 2         0.97   8589.06493     7163.035    1342.6635     83.36647
#> 3         0.00   5225.84843     4767.323     458.5248      0.00000
#> 4           NA   2961.75209     2878.386           NA           NA
#> 5           NA   1426.02993     1426.030           NA           NA
#> 6           NA    458.52517           NA           NA           NA
#> 7           NA     83.36647           NA           NA           NA

Part 5: Visualization

A static bathymetry map

The rLakeHabitat package includes built-in functions for plotting interpolated DEMs, starting with bathyMap(), followed by standalone code for more customized raster plots.

# A basic bathymetry map - the density of contour lines can be adjusted via the 'by' argument
bathyMap(IDW, contours = TRUE, units = "m", by = 1)


# These maps are fairly basic on their own. Helpful references like a scale bar,
# compass, or a background providing geographic context can be added with a few
# additional packages.

# First, reproject the DEM to Web Mercator (EPSG:3857) to match the map tiles used below
IDW_3857 <- terra::project(IDW, "EPSG:3857")

# Then turn the DEM into a data frame to remove NA cells around the lake 
IDW_df <- as.data.frame(IDW_3857, xy = T, na.rm = T)


# Next, get the extent of the raster as a bounding box
bbox_sf <- sf::st_as_sfc(sf::st_bbox(IDW_3857))

# Optional: retrieve the necessary satellite tiles
# sat_tiles <- maptiles::get_tiles(bbox_sf, provider = "OpenTopoMap", zoom = 14, crop = TRUE) #zoom 13 (broader) or zoom 15 (finer)


# Now we can put it all together
ggplot() +
  # tidyterra::geom_spatraster_rgb(data = sat_tiles) + #include for sat_tiles
  geom_raster(data = IDW_df, aes(x = x, y = y, fill = lyr1)) +
  scale_fill_continuous(name = "Depth (m)", low = "blue", high = "lightblue", na.value="transparent", trans = "reverse") +
  tidyterra::geom_spatraster_contour(data = IDW_3857, breaks = seq(0, as.numeric(max(values(IDW_3857, na.rm = T))), 5), color = "black") +
  ggspatial::annotation_scale(location = "br", 
                              width_hint = 0.3,
                              text_col = "black",bar_cols = c("black", "white"),
                              pad_x = unit(5, "cm"), pad_y = unit(2, "cm")) +
  ggspatial::annotation_north_arrow(location = "tl", 
                                    which_north = "true",
                                    style = ggspatial::north_arrow_fancy_orienteering(
                                      text_col = "black", 
                                      fill = c("white", "black")),
                                    pad_x = unit(2.5, "cm"), pad_y = unit(2, "cm"),
                                    height = unit(3, "cm"), width = unit(3, "cm")) +
  coord_sf(crs = 3857) +
  labs(title = "Bathymetric DEM", x = NULL, y = NULL) +
  theme_minimal()



# and we can zoom out if we want more context

# # add buffer around boundary box
# bbox_buffered <- sf::st_buffer(bbox_sf, dist = 4000)  # buffer distance in meters
# 
# # get new map tiles
# sat_tiles <- maptiles::get_tiles(bbox_buffered, provider = "OpenTopoMap", zoom = 13, crop = TRUE)
# 
# # save the extend of our new boundary box to add to our plot
# buff_ext <- sf::st_bbox(bbox_buffered)
# 
# ggplot() +
#   tidyterra::geom_spatraster_rgb(data = sat_tiles) +
#   geom_raster(data = IDW_df, aes(x = x, y = y, fill = lyr1)) +
#   scale_fill_continuous(name = "Depth (m)", low = "blue", high = "lightblue", na.value = "transparent", trans = "reverse") +
#   tidyterra::geom_spatraster_contour(data = IDW_3857, breaks = seq(0, as.numeric(max(values(IDW_3857, na.rm = T))), 5), color = "black") +
#   ggspatial::annotation_scale(location = "br", 
#                               width_hint = 0.3,
#                               text_col = "black", bar_cols = c("black", "white"),
#                               pad_x = unit(5, "cm"), pad_y = unit(2, "cm")) +
#   ggspatial::annotation_north_arrow(location = "tl", 
#                                     which_north = "true",
#                                     style = ggspatial::north_arrow_fancy_orienteering(
#                                       text_col = "black", 
#                                       fill = c("white", "black")),
#                                     pad_x = unit(2.5, "cm"), pad_y = unit(2, "cm"),
#                                     height = unit(3, "cm"), width = unit(3, "cm")) +
#   coord_sf(crs = 3857,
#            xlim = c(buff_ext["xmin"], buff_ext["xmax"]),
#            ylim = c(buff_ext["ymin"], buff_ext["ymax"]),
#            expand = FALSE) +
#   labs(title = "Bathymetric DEM", x = NULL, y = NULL) +
#   theme_minimal()

An animated view across water levels

animBathy() produces a gganimate object showing how littoral (or total) surface area changes as water level drops, useful for reservoirs with substantial drawdown.

animBathy(IDW, units = "m", littoral = TRUE, photic = photic, by = 1)

Plot habitat changes

We can also look at the results from the littoral and pelagic habitat analyses above, combining datasets to see percent change in habitat area and volume across water levels.

habitat <- lit %>%
  merge(c(vol, sdi), by = "depth") %>%
  select(depth, perc_lit_tot, perc_epi_tot, perc_hyp_tot, SDI)

# rescale SDI onto the same 0-100 visual range as the percent variables, computed
# from the actual combined SDI range rather than a fixed assumed scale
sdi_scale_factor <- max(habitat$SDI, na.rm = TRUE) / 100

habitat %>%
  mutate(SDI_scaled = SDI / sdi_scale_factor) %>%
  ggplot() +
  geom_line(aes(x = depth, y = perc_epi_tot, color = 'Epilimnion'), lwd = .8) +
  geom_line(aes(x = depth, y = perc_hyp_tot, color = 'Hypolimnion'), lwd = .8) +
  geom_line(aes(x = depth, y = perc_lit_tot, color = 'Littoral'), lwd = .8) +
  geom_line(aes(x = depth, y = SDI_scaled, color = 'SDI'), lwd = .8) +
  scale_color_manual(values = c('Epilimnion' = 'red', 'Hypolimnion' = 'blue', 'Littoral' = 'green', 'SDI' = 'purple')) +
  scale_y_continuous(name = "Percent of Total",
                     sec.axis = sec_axis(~ . * sdi_scale_factor, name = "SDI")) +
  labs(x = "Water Level Fluctuation (m)", color = "Habitat") +
  theme_bw() +
  theme(strip.text.x = element_text(size = 14),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.position = 'right')
#> Warning: Removed 3 rows containing missing values or values outside the scale range
#> (`geom_line()`).
#> Removed 3 rows containing missing values or values outside the scale range
#> (`geom_line()`).

Part 7: Saving Your Work

genStack() splits a DEM into a stack of binary (wet/dry) layers at each depth increment. A raster stack is a useful format for further GIS work, or for archiving a raster product. Cloud Optimized GeoTIFF ("COG") is the default file type - COGs are formatted so that remote or cloud-hosted tools like web viewers can pull only the portion, or stack layers, needed instead of downloading the whole file. COGs do take up slightly more storage space than standard GeoTIFF files (.tif, "GTiff").

# generate raster stack with layers at 2 m intervals
stack <- genStack(IDW, by = 2, save = FALSE)
stack
#> class       : SpatRaster 
#> size        : 62, 89, 5  (nrow, ncol, nlyr)
#> resolution  : 0.0005464108, 0.0005464108  (x, y)
#> extent      : -89.27918, -89.23055, 42.94888, 42.98276  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> names       :     lyr1,     lyr1,     lyr1,     lyr1,     lyr1 
#> min values  : 0.000000, 0.000000, 2.001398, 4.002190, 6.031895 
#> max values  : 6.961307, 6.961307, 6.961307, 6.961307, 6.961307

# to actually write it to disk as a cloud-optimized GeoTIFF:
# genStack(IDW, by = 2, save = TRUE, file_name = "Lake", file_type = "COG")

Bathymetric contour lines can similarly be generated and exported directly with saveContours(), in GPX, KML, shapefile, or GeoPackage format:

# generate and save 5 m contour lines for our lake as a shapefile
saveContours(IDW, by = 5, units = "m", file_name = "LakeContours", file_type = "shp")

# or as a GPX, for use in a handheld GPS unit
saveContours(IDW, by = 5, units = "m", file_name = "LakeContours", file_type = "gpx")

Part 8: Other Features

Assigning Sampling Locations

Given a DEM, partitionSites() can be used to design a depth-stratified, spatially balanced sampling program - for example, gillnet sets spread evenly across depth zones. It does this using Halton Iterative Partitioning (HIP), a spatially balanced sampling method: within each depth bin specified, it selects locations that are well spread out rather than randomly clumped. HIP accomplishes this by dividing the DEM into n boxes, where n = boxes per bin x number of bins, then numbering the boxes using a Halton sequence, which assigns numbers evenly across space. One sample location is then randomly assigned within each box, ensuring samples are evenly distributed over the whole area. No outline shapefile is needed - a boundary for the plot is derived automatically from the DEM’s own footprint.

If the waterbody’s water level is currently below full pool (the level the DEM represents), water_level_drop rebuilds the DEM before sampling: depths are reduced by the specified amount, and any cell that would now be above the water line is excluded from the candidate pool.

sites <- partitionSites(dem = IDW, 
                        depth_bins = c(0, 2, 4, Inf),
                        n_per_bin = 2,     # or c(1,2,1)
                        min_spacing = 10,
                        plot = T, seed = 123)
#>   |                                                                              |                                                                      |   0%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================================| 100%


sites$locations
#>           x        y depth_bin bin_index
#> 1 -89.26579 42.97374     [0,2)         1
#> 2 -89.24066 42.95571     [0,2)         1
#> 3 -89.26743 42.96227     [2,4)         2
#> 4 -89.24831 42.96773     [2,4)         2
#> 5 -89.26197 42.96664   [4,Inf]         3
#> 6 -89.25486 42.96172   [4,Inf]         3


# These coordinates and their depth bin can also be saved to a file (e.g. .gpx) for external use

sites_df <- sites$locations %>%
  mutate(name = paste0("Site_", row_number()))  # partitionSites() doesn't label sites by name, so add one before export

# Convert to an sf object.
# IMPORTANT: crs must match whatever CRS your x/y coordinates are actually in

# Check the DEM CRS:
crs(IDW) #currently in EPSG:4326 - WGS84
#> [1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"

# Convert to coordinates in the CRS of the DEM
sites_sf <- sf::st_as_sf(sites_df, coords = c("x", "y"), crs = 4326)

# Reproject to WGS84 for .gpx format - don't need to do this for our data
sites_wgs84 <- sf::st_transform(sites_sf, crs = 4326)

# GDAL's GPX driver only recognizes a fixed set of standard waypoint fields (name, desc, cmt, sym, ele, time, etc.)
# "depth" isn't a standard field, so we change it to "ele" (elevation)

sites_gpx <- sites_wgs84 %>%
  tidyr::extract(
    depth_bin,
    into = c("depth_min", "depth_max"),
    regex = "\\[?\\(?(-?(?:[0-9.]+|Inf)),\\s*(-?(?:[0-9.]+|Inf))[\\)\\]]?",
    remove = FALSE,
    convert = TRUE) %>%
  rename(ele = depth_min) %>%
  select(name, ele, geometry)

# save
# sf::st_write(sites_gpx, "sites.gpx", driver = "GPX", delete_dsn = TRUE)

Appendix: Retrieving outlines via API without downloading full datasets

Waterbody outlines can also be retrieved directly from the National Hydrography Dataset via ArcGIS REST API, without needing to download the full dataset. This is included for reference, but is not run as part of this vignette, since it requires an internet connection and can have an extended runtime.

# function to get state lakes
get_state_lakes <- function(state_name) {
  
  read_arcgis <- function(base_query_url) {
  offset <- 0
  page_size <- 2000
  all_features <- list()

  repeat {
    query_url <- paste0(base_query_url, "&resultRecordCount=", page_size,
                         "&resultOffset=", offset)
    response <- GET(query_url, add_headers(`User-Agent` = "Mozilla/5.0"))
    if (status_code(response) != 200) stop("Failed to fetch: ", query_url)

    content_text <- content(response, as = "text", encoding = "UTF-8")
    page <- st_read(content_text, quiet = TRUE)

    if (nrow(page) == 0) break
    all_features[[length(all_features) + 1]] <- page
    offset <- offset + page_size

    if (nrow(page) < page_size) break
  }

  if (length(all_features) == 0) return(NULL)
  do.call(rbind, all_features)
}
  
  # State boundary
  state_query <- paste0(
    "https://services.arcgis.com/P3ePLMYs2RVChkJx/ArcGIS/rest/services/",
    "US_States_boundaries/FeatureServer/3/query",
    "?where=STATE_NAME%3D'", state_name, "'",
    "&outFields=*&returnGeometry=true&f=geojson"
  )
  
  state <- read_arcgis(state_query)
  
  bb <- st_bbox(state)
  bbox_string <- paste(bb["xmin"], bb["ymin"], bb["xmax"], bb["ymax"], sep = ",")
  
  # Waters
  waters_query <- paste0(
    "https://services.arcgis.com/P3ePLMYs2RVChkJx/ArcGIS/rest/services/",
    "NHDPlusV21/FeatureServer/1/query",
    "?where=1%3D1",
    "&geometry=", bbox_string,
    "&geometryType=esriGeometryEnvelope",
    "&inSR=4326",
    "&spatialRel=esriSpatialRelIntersects",
    "&outFields=*&returnGeometry=true&f=geojson"
  )
  
  waters <- read_arcgis(waters_query)
  
  lakes <- waters %>%
    st_make_valid() %>%
    st_filter(state)
  
  return(list(state = state, lakes = lakes))
}


#subset output and verify
wy_data <- get_state_lakes("Wyoming")

# select specific waterbodies of interest
AL_outline <- wy_data$lakes %>%
  filter(GNIS_NAME %in% c("Alcova Reservoir")) %>%
  select(GNIS_NAME, geometry)

For full documentation on any function, use ?functionName from the console.