Sunt disponibile două metode pentru a estima volumul total al unui rezervor folosind date batimetrice. Metoda analogică folosește linii de contur pentru a desena secțiuni transversale pe hârtie milimetrică la intervale regulate. Apoi utilizați un planimetru sau formula șiretului pentru a măsura suprafața fiecărei secțiuni și calculați volumul total înmulțind suprafața cu distanța dintre secțiuni. Am folosit această metodă în urmă cu mai bine de treizeci de ani când lucram la proiecte de dragare. Dar tehnologia a mers mai departe.
A doua metodă este mai numerică și există două metode: interpolați o grilă sau construiți o rețea neregulată triunghiulară.
Batimetria furnizează date numai pentru liniile de pe hartă, așa cum este vizualizat mai sus. Pentru a estima analitic aprovizionarea completă, putem suprapune rezervorul cu o grilă și interpola adâncimea pentru fiecare celulă. Putem apoi estima volumul complet de aprovizionare prin însumarea produsului fiecarei suprafețe de celule și adâncime pentru fiecare element de grilă. Când avem o grilă cu n elemente, și x şi y indicați dimensiunea elementului și di adâncimea unui element, volumul se calculează cu:
$$V = sum_{i=1}^{n} {x times y times d_i}$$
Pachetul akima oferă funcționalitate pentru interpolarea datelor cu distanță neregulată și regulată, făcându-l ideal pentru problema noastră. Mai întâi, setăm distanța la 10 metri. Cu cât intervalul este mai mic, cu atât estimarea este mai precisă. Dar spațierea mică are un cost, deoarece interpolarea va dura ceva timp pentru a finaliza.
Odată ce avem o grilă cu o adâncime pentru fiecare locație, putem estima volumul complet de aprovizionare prin înmulțirea adâncimii fiecărei locații cu aria sa grilei și însumând toate grilele.
Această analiză oferă un volum de 71,50 × 106 m3care este aproape de 19 miliarde de galoane americane (71,92 × 106 m3) raportate de Maryland Geological Survey.
# Extrapolate gridded version
library(akima)
# Grid size (m)
dx <- 10
dy <- 10
# Regular grid coordinates
x_seq <- seq(min(prettyboy$easting), max(prettyboy$easting), by = dx)
y_seq <- seq(min(prettyboy$northing), max(prettyboy$northing), by = dy)
# Interpolate scattered (x, y, depth) -> regular grid
# no extrapolation outside convex hull
interp_res <- with(prettyboy,
interp(x = easting,
y = northing,
z = depth_m,
xo = x_seq,
yo = y_seq,
duplicate = "strip",
linear = TRUE,
extrap = FALSE))
# Volume
prettyboy_grid <- interp2xyz(interp_res, data.frame = TRUE) |>
as_tibble() |>
rename(easting = x, northing = y, depth_m = z)
grid_volume_m3 <- sum(prettyboy_grid$depth_m * dx * dy, na.rm = TRUE)
grid_volume_m3 / 100^3
Putem folosi valorile interpolate pentru a crea o hartă frumoasă a rezervorului cu linii de contur și secțiuni transversale. Secțiunile transversale sunt luate în patru puncte aleatorii ale rețelei.
# Visualise Cross sections
set.seed(1234)
northing_sections <- sample(prettyboy_grid$northing, 4)
prettyboy_sections <- filter(prettyboy_grid, northing %in% northing_sections)
max_depth <- round(max(prettyboy_sections$depth_m, na.rm = TRUE))
map <- ggplot(filter(prettyboy_grid)) +
aes(easting, northing) +
geom_raster(aes(fill = depth_m), interpolate = TRUE) +
geom_contour(aes(z = depth_m), colour = "white", alpha = 0.4) +
geom_hline(yintercept = northing_sections, col = "orange") +
coord_equal() +
scale_x_continuous(labels = scales::label_comma()) +
scale_y_continuous(labels = scales::label_comma()) +
scale_fill_gradientn(
colours = c(NA, bathymetry_colours(-1)),
na.value = "transparent",
guide = guide_colourbar(reverse = TRUE),
name = "Depth (m)") +
theme_minimal(base_size = 28) +
labs(title = "Prettyboy Reservoir",
caption = "Source: Maryland Geological Survey",
x = "Easting (m)", y = "Northing (m)") +
theme(legend.position = "bottom",
legend.key.width = unit(2, 'cm'))
sections <- ggplot(prettyboy_sections) +
aes(easting, (fsl - depth_m)) +
geom_area(fill = bathymetry_colours(1), alpha = 0.5) +
scale_x_continuous(labels = scales::label_comma()) +
facet_wrap(~rev(northing), ncol = 1) +
coord_cartesian(ylim = c(fsl - 1.1 * max(prettyboy_sections$depth_m, na.rm = TRUE), fsl)) +
theme_minimal(base_size = 28) +
labs(x = "Easting (m)", y = "Elevation (m)")
library(gridExtra)
png(filename = "prettyboy-sections.png", width = 1920, height = 940, units = "px")
grid.arrange(map, sections, ncol = 2)
dev.off()
Secțiuni transversale ale rezervorului Prettyboy.
A doua metodă pentru calcularea volumului total este utilizarea unei rețele neregulate triangulate (TIN). Un TIN este o reprezentare a unei suprafețe formate din fațete triunghiulare. Principala diferență cu valorile interpolate este că acest model este neregulat, astfel încât fiecare triunghi are o zonă proiectată diferită.
Pachetul de geometrie implementează triangulația Delaunay, unde fiecare rând conține indici ai trei puncte care formează un triunghi. În pasul următor, folosim map_dfr() funcția de a crea un cadru de date cu numărul de index al fiecărui triunghi, cele trei coordonate și adâncimile asociate.
Pentru a măsura volumul, putem folosi formula șiretului pentru a calcula aria fiecărui triunghi:
$$A_i = frac{1}{2}| x_{i_1} (y_{1_2} – y_{i_3}) + x_{i_2} (y_{1_3} – y_{i_1}) + x_{i_3} (y_{1_1} – y_{i_2})|$$
Volumul elementului este aria înmulțită cu adâncimea medie:
$$V = sum{A_i frac{sum d_i}{3}}$$
Această metodă oferă același volum, dar este puțin mai rapid de procesat, deoarece nu creăm date noi.
# Triangulated Irregular Network
library(geometry)
library(purrr)
# 2D coordinates
coords2d <- as.matrix(prettyboy(, c("easting", "northing")))
# Delaunay triangulation
tri <- delaunayn(coords2d)
# Convert the triangles into a long data frame for ggplot
tin_df <- map_dfr(seq_len(nrow(tri)), function(i) {
idx <- tri(i, )
data.frame(
tri_id = i,
easting = prettyboy$easting(idx),
northing = prettyboy$northing(idx),
depth_m = prettyboy$depth_m(idx))})
triangle_area_shoelace <- function(x, y) {
abs(x(1) * (y(2) - y(3)) +
x(2) * (y(3) - y(1)) +
x(3) * (y(1) - y(2))) / 2}
tin_volume <- tin_df |>
dplyr::group_by(tri_id) |>
dplyr::summarise(
area = triangle_area_shoelace(easting, northing),
mean_depth = mean(depth_m),
volume = area * mean_depth,
.groups = "drop")
tin_volume_m3 <- sum(tin_volume$volume)
tin_volume_m3 / 1e6
