--- title: "Introduction to rLakeHabitat" author: "Tristan Blechinger & Sean Bertalot" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 2 vignette: > %\VignetteIndexEntry{Introduction to rLakeHabitat} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 5, message = FALSE) ``` ## 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: ```{r} library(rLakeHabitat) library(terra) library(tidyterra) library(ggspatial) library(sf) library(dplyr) library(httr) library(spbal) library(ggplot2) library(mapview) ``` ```{r} 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. ```{r} # Load Alcova outline shapefile outline <- read_sf(extdata_path("example_outline.shp")) crs(outline) # verify input CRS 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) 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 design$total_transect_length ``` Once a spacing is chosen, the resulting transects can be saved to a `.gpx` file for use in most common GPS units: ```{r, eval=F} # 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: 1. Load data 2. Convert contours to points (if necessary) 3. Clean data: remove outlying depths (e.g., 0 or 1000) 4. 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. ```{r} # load contours and sample points contours <- read_sf(extdata_path("example_contour.shp")) head(contours) 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) crs(outline) #EPSG:4326 crs(depths) #No CRS, not projected yet -- recorded in WGS84 (EPSG:4326) -- good # # 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. ```{r} # Basic IDW interpolation IDW <- interpBathy(outline, rare, "x", "y", "z", separation = 10, res = 50, nmax = 8) plot(IDW, main = "IDW-interpolated DEM") # Basic OK interpolation OK <- interpBathy(outline, rare, "x", "y", "z", separation = 20, res = 50, method = "OK", nmax = 8) 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. ```{r} 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 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") } ``` ### 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. ```{r} # 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) ``` ### 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. ```{r} 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) tuning$results # 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 tuning$recommended_res ``` --- ## 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`): ```{r} # 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") ``` 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. ```{r} 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. ```{r} # 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 # shoreline development index (how convoluted the shoreline is relative to a circle) sdi <- calcSDI(IDW, units = "m", by = 1) sdi # littoral vs. pelagic volume littoralVol(IDW, photic = photic, depthUnits = "m", by = 1) ``` Base `terra` functions can also be used directly for metrics like average slope or true (terrain-corrected) surface area. ```{r} # 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 # 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 ``` ### 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: ```{r} vol <- calcVolume(IDW, thermo_depth = 3, depthUnits = "m", by = 1) vol # 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) ``` --- ## 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. ```{r} # 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. ```{r, eval = FALSE} 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. ```{r} 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') ``` --- ## 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"`). ```{r} # generate raster stack with layers at 2 m intervals stack <- genStack(IDW, by = 2, save = FALSE) stack # 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: ```{r, eval = FALSE} # 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. ```{r} 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) sites$locations # 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 # 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. ```{r, eval = FALSE} # 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.