The General Lake Model (GLM) requires information on the depth-area relationship for a given lake. In GLM this is referred to as hypsography (rather than bathymetry) and has the following description in the glmtools
documentation for get_hypsography
:
library(glmtools)
library(printr)
?get_hypsography
Let's start by investigating (reverse-engineering) the hypsography provided for Sparkling Lake (Wisconsin) in the run_example_sim
function.
sim_folder <- run_example_sim("inst/extdata", verbose = FALSE) nml_file <- file.path("inst/extdata", "glm2.nml")
nml_file <- file.path("~/R/scripts/glmtools/inst/extdata", "glm3.nml")
nml <- read_nml(nml_file) hp <- get_hypsography(file = nml_file)[,2:1]
# plot(hp, ylim = rev(range(hp[,2])))
hp_diff <- cbind(c(diff(hp$depths), 0), hp$areas) calc_layer_vol <- function(hp_diff, i){ top_area <- hp_diff[i, 1] bot_area <- hp_diff[i + 1, 2] # volume of a layer approximated as a truncated cone (hp_diff[i, 1] / 3) * (top_area + bot_area + sqrt(top_area * bot_area)) } vol_layers <- sapply(seq_len(nrow(hp_diff) - 1), function(i) calc_layer_vol(hp_diff, i)) (vol_estimated <- sum(vol_layers))
(nlayers <- nrow(hp)) (max_depth <- tail(hp$depths, 1)) (surface_area <- head(hp$areas, 1)) res <- suppressWarnings( glmutils::generate_hypsography(nlayers, max_depth, surface_area)) plot(res[,c("area", "depth_linear")], ylim = c(max(res$depth_linear), min(res$depth_linear)), ylab = "depth") points(res[,c("area", "depth_ellipsoid")], col = "red", pch = 19) legend("topleft", legend = c("linear", "ellipsoid"), col = c("black", "red"), pch = c(21, 19))
# Get max depth tail(hp, 1)
This value matches the max depth reported by the Wisconsin DNR.
get_nml_value(nml, "bsn_len") get_nml_value(nml, "bsn_wid")
library(nhdR) # Get NHD Waterbody Polygons # nhd_plus_get(vpu = 7, "NHDSnapshot") dt <- nhd_plus_load(vpu = 7, "NHDSnapshot", "NHDWaterbody") spark_outline <- dt[grep("Sparkling", dt$GNIS_NAME),] plot(spark_outline$geometry, main = "Sparkling Lake") suppressPackageStartupMessages(library(lakemorpho)) library(raster) # Get length and width spark_outline_sp <- as(spark_outline, "Spatial") r <- rasterize(spark_outline_sp, raster(spark_outline_sp)) spark_outline_sp <- spTransform(spark_outline_sp, CRS("+proj=utm +zone=10 +datum=WGS84")) r <- projectRaster(r, crs = "+proj=utm +zone=10 +datum=WGS84") spark_lm <- lakeSurroundTopo(spark_outline_sp, r) calcLakeMetrics(spark_lm, bearing = 45, pointDens = 250)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.