inst/doc/rLakeHabitat_Introduction.R

In rLakeHabitat: Interpolate Bathymetry and Quantify Physical Aquatic Habitat

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 5, message = FALSE)

## -----------------------------------------------------------------------------
library(rLakeHabitat)
library(terra)
library(tidyterra)
library(ggspatial)
library(sf)
library(dplyr)
library(httr)
library(spbal)
library(ggplot2)
library(mapview)

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

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


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

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


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

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

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

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


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

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

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

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

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

## ----eval = FALSE-------------------------------------------------------------
# animBathy(IDW, units = "m", littoral = TRUE, photic = photic, by = 1)

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

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

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

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

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