Lake Rangsdorf is located at the southern margin of the Teltow Plateau and a former proglacial drainage system, the Rangsdorf-Thyrow-Abflussbahn. Its basin has originally emerged as part of a subglacial drainage system, still visible at the panhandle Krumme Lanke and its northern continuation–the Glasow-Bach valley, whereas the main basin has been designed as a kettle hole. The northern part of Krumme Lanke still exhibits its basin shape as a glacial tunnel lake (max. lake depth 9 m with and average depth of roughly 3.5 m), where the southern main lake shows the typical morphology of a kettle lake silting up to an average depth of approx. 0.9 m (max. 2.7 m). As most of the other shallow lakes in Brandenburg, lake Rangsdorf is an eutrophic to hypertrophic lake suffering from heavy algae blooms in summer season.
Lake bathymetry and chemical sediment data of Lake Rangsdorf was obtained during a field campaign in 2017 led by Kai Hartmann, Department of Earth Sciences at Freie Universität Berlin.
Spatial data
The spatial data for lake Rangsdorf is provided as a zipped file,
LakeRangsdorf.zip. The file contains several ESRI
shapefiles that we download from the server.
# build a temporary folder on disk
temp <- tempfile()
download_url <- "https://userpage.fu-berlin.de/soga/data/raw-data/spatial/"
zipfile <- "LakeRangsdorf.zip"
## download the file
download.file(paste0(download_url, zipfile), temp, mode = "wb")
## unzip the file(s)
unzip(temp)
## close file connection
unlink(temp)The file shoreline_poly.shp is a polygon shapefile of
lake Rangsdorf. We load the file into memory using the
st_read() function from the sf package and plot its geometry to get
a first visual impression of the data.
library(sf)
## read in the data
lake <- st_read("shoreline_poly.shp", quiet = TRUE)
plot(lake$geometry)
In addition to the polygon shapefile, the zip-file contains the file
shoreline.shp which is a shapefile of the shoreline of lake
Rangsdorf. We load the file into memory using again the
st_read() function.
## shoreline multiline geometry
shoreline <- st_read("shoreline.shp", quiet = TRUE)Then we transform it into a POINT geometry type using
the st_cast() function from the sf package.
This allows us to assign each point on the shoreline a depth of zero,
which becomes useful later, when we compute a bathymetric model of the
lake.
## shoreline points
shoreline_point_geometry <- st_cast(x = shoreline, to = "POINT")We add a feature vector denoted as depth_m to the
POINT object and assign the vector values of zero.
## add vector of zero depth
shoreline_pts <- st_sf(
depth_m = rep(0, length(shoreline_point_geometry$geometry)),
geometry = shoreline_point_geometry$geometry
)There is one more spatial data set we load into memory. The file
sonar.shp contains point measurements of the lake depth
obtained during a sonar survey. We load the file and plot it. Note that
we make sure that the data does not contain any invalid values, such as
zero depth.
sonar <- st_read("sonar.shp", quiet = TRUE)
sonar <- sonar[sonar$depth_m > 0, ]
plot(sonar)
Chemical sediment data
In addition to spatial data we provide a chemical sediment data set
(LakeRangsdorf.csv).
download_url <- "https://userpage.fu-berlin.de/soga/data/raw-data/"
lake_Rangsdorf <- read.csv(paste0(
download_url,
"LakeRangsdorf.csv"
))
head(lake_Rangsdorf, 8)## ID x y depth_m visual_depth_cm elec_cond_muS pH LOI550 LOI990
## 1 10a 391829 5794052 0.68 34 305 7.19 1.04 0.20
## 2 11a 391824 5794187 0.66 28 274 7.33 1.27 0.19
## 3 12a 391184 5795337 1.93 15 680 7.25 27.32 17.46
## 4 13b 391312 5795462 3.98 15 825 7.22 36.93 18.25
## 5 14b 391729 5795736 4.95 18 903 7.05 38.33 15.50
## 6 15b 391892 5795771 0.93 20 416 7.31 2.33 0.08
## 7 16a 391430 5795473 3.92 18 1060 6.82 37.79 17.67
## 8 17a 391846 5795878 8.40 18 1278 7.02 36.74 16.42
## LOI C N Ca_mg_g Fe_mg_g K_mg_g Pb_mug_g Mg_mg_g Mn_mug_g Na_mg_g
## 1 98.51 0.85 NA 9.83 0.52 0.22 1.37 0.15 15.39 0.26
## 2 98.29 0.93 NA 3.78 0.42 0.13 1.81 0.13 12.02 0.22
## 3 32.99 19.00 1.18 172.12 4.73 0.34 18.49 1.41 481.61 0.54
## 4 21.60 24.40 1.85 182.74 8.58 0.43 34.19 1.47 454.44 0.70
## 5 26.44 23.70 1.89 193.42 10.12 0.50 31.00 1.42 546.69 0.88
## 6 97.49 1.79 NA 4.22 0.81 0.14 4.14 0.20 61.32 0.32
## 7 22.06 24.70 1.94 175.28 8.07 0.55 30.92 1.46 479.68 0.74
## 8 25.94 23.10 1.82 187.53 9.26 0.64 28.96 1.45 653.76 0.72
## PO4_mg_g S_mg_g Sr_mug_g Zi_mug_g
## 1 0.30 0.77 7.87 11.46
## 2 0.37 0.43 6.77 17.41
## 3 4.83 8.24 238.28 84.19
## 4 4.93 17.18 261.05 94.92
## 5 5.31 20.03 245.53 98.52
## 6 1.18 1.26 8.78 30.30
## 7 5.34 17.21 261.54 87.90
## 8 5.97 17.07 257.63 90.17The data set contains 32 measurements on sediment samples taken at
different locations within the lake. We use the st_as_sf()
function to convert the data.frame object into a
sf object. Note that add the EPSG identifier \(32633\) to the function call,
accounting for the fact that the coordinates are given as Universal Transverse Mercator, zone 33
coordinates.
lake_data_sf <- st_as_sf(lake_Rangsdorf,
coords = c("x", "y"),
crs = 32633
)
head(lake_data_sf)## Simple feature collection with 6 features and 21 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 391184 ymin: 5794052 xmax: 391892 ymax: 5795771
## Projected CRS: WGS 84 / UTM zone 33N
## ID depth_m visual_depth_cm elec_cond_muS pH LOI550 LOI990 LOI C N
## 1 10a 0.68 34 305 7.19 1.04 0.20 98.51 0.85 NA
## 2 11a 0.66 28 274 7.33 1.27 0.19 98.29 0.93 NA
## 3 12a 1.93 15 680 7.25 27.32 17.46 32.99 19.00 1.18
## 4 13b 3.98 15 825 7.22 36.93 18.25 21.60 24.40 1.85
## 5 14b 4.95 18 903 7.05 38.33 15.50 26.44 23.70 1.89
## 6 15b 0.93 20 416 7.31 2.33 0.08 97.49 1.79 NA
## Ca_mg_g Fe_mg_g K_mg_g Pb_mug_g Mg_mg_g Mn_mug_g Na_mg_g PO4_mg_g S_mg_g
## 1 9.83 0.52 0.22 1.37 0.15 15.39 0.26 0.30 0.77
## 2 3.78 0.42 0.13 1.81 0.13 12.02 0.22 0.37 0.43
## 3 172.12 4.73 0.34 18.49 1.41 481.61 0.54 4.83 8.24
## 4 182.74 8.58 0.43 34.19 1.47 454.44 0.70 4.93 17.18
## 5 193.42 10.12 0.50 31.00 1.42 546.69 0.88 5.31 20.03
## 6 4.22 0.81 0.14 4.14 0.20 61.32 0.32 1.18 1.26
## Sr_mug_g Zi_mug_g geometry
## 1 7.87 11.46 POINT (391829 5794052)
## 2 6.77 17.41 POINT (391824 5794187)
## 3 238.28 84.19 POINT (391184 5795337)
## 4 261.05 94.92 POINT (391312 5795462)
## 5 245.53 98.52 POINT (391729 5795736)
## 6 8.78 30.30 POINT (391892 5795771)Finally, in order to contextualize the spatial data for lake
Rangsdorf we make use of the mapview library.
library(mapview)
mapview(shoreline, color = "red", highlight = FALSE, label = "Shoreline") +
mapview(sonar, color = "black", highlight = FALSE, label = "Sonar transect", cex = 0.1) +
mapview(lake_data_sf, cex = 3)Lake bathymetry
So now that we loaded the data we can do some interesting
computations. Let us compute a bathymetric model by interpolation of the
measured depths. First, we need to prepare the data, but this is easy
due to the structure of sf objects. We can simply combine
the data sets row-wise by calling the rbind() function.
lake_depth <- rbind(
sonar,
lake_data_sf[, "depth_m"],
shoreline_pts
)
head(lake_depth, 3)## Simple feature collection with 3 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 391578.1 ymin: 5794653 xmax: 391585.9 ymax: 5794663
## Projected CRS: WGS 84 / UTM zone 33N
## depth_m geometry
## 1 1.62 POINT (391585.9 5794653)
## 2 1.57 POINT (391582.3 5794657)
## 3 1.57 POINT (391578.1 5794663)The depth values are stored as positive numbers, which is fine, but for the computation of a bathymetric model we want them as negative numbers, indicating the depth below the lake surface.
lake_depth$depth_m <- lake_depth$depth_m * -1Moreover, we want to make sure that we do not have any missing values
in our data set. Hence, we subset the data using the
complete.cases() function.
lake_depth <- lake_depth[complete.cases(lake_depth$depth_m), ]Recall that we want to compute a bathymetric model for lake
Rangsdorf. Therefore, we need to provide a bunch of coordinates at which
we actually want to interpolate the depth. Therefore, we create a
regular grid of points using the expand.grid() function. We
specify the grid based on the bounding box of the lake polygon
(lake). Furthermore, we specify the grid step size (cell
size) to be 10 meter.
extent_lake <- st_bbox(lake)
grid_lake <- expand.grid(
x = seq(
from = extent_lake[1],
to = extent_lake[3],
by = 10
),
y = seq(
from = extent_lake[2],
to = extent_lake[4],
by = 10
)
)
head(grid_lake)## x y
## 1 390329.9 5792973
## 2 390339.9 5792973
## 3 390349.9 5792973
## 4 390359.9 5792973
## 5 390369.9 5792973
## 6 390379.9 5792973tail(grid_lake)## x y
## 55396 392039.9 5796093
## 55397 392049.9 5796093
## 55398 392059.9 5796093
## 55399 392069.9 5796093
## 55400 392079.9 5796093
## 55401 392089.9 5796093We now add a vector of zeros to the data frame, representing the not
yet known depth. Then we convert the data.frame object into
an sp object, to be more precisely into a
SpatialPointsDataFrame object.
library(sp)
grid_lake["depth_m"] <- 0
coordinates(grid_lake) <- ~ x + y
proj4string(grid_lake) <- st_crs(lake)$proj4string
grid_lake## class : SpatialPointsDataFrame
## features : 55401
## extent : 390329.9, 392089.9, 5792973, 5796093 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs
## variables : 1
## names : depth_m
## min values : 0
## max values : 0Now, we are ready for interpolation. We use the
interpp() function from the akima package. The
function implements a Delaunay triangulation on a set of points for
irregularly spaced input data. The applied algorithm is also known as
Algorithm 751, and was published by R.Renka in 1996.
library(akima)
lake_interpolation <- interpp(as(lake_depth, "Spatial"),
z = "depth_m",
xo = grid_lake,
linear = TRUE,
duplicate = "mean"
)Done!
In the way we specified the interpp() function it
returns a SpatialPointsDataFrame object. For the purpose of
visualization we transform this object into a
SpatialPixelsDataFrame object by applying the
gridded() function.
gridded(lake_interpolation) <- TRUENow, we may mask the interpolated SpatialPixelsDataFrame
object by the lake polygon (lake) and plot it using the
levelplot() function.
library(raster)
library(rasterVis)
masked_lake_interpolation <- mask(
raster(lake_interpolation),
as(lake, "Spatial")
)
masked_lake_interpolation[masked_lake_interpolation > 0] <- NA
levelplot(masked_lake_interpolation, margin = FALSE) +
latticeExtra::layer(sp.polygons(as(lake$geometry, "Spatial")))
Nice plot! But we may even plot the interpolated data in 3D, using
the wireframe() function. Please note that
wireframe() is a very versatile function, providing a lot
of optional arguments. Look up the documentation for further details
(type ?wireframe into your console). Here we feed the
function with a matrix object, denoted as
mat_lake_interpolation.
mat_lake_interpolation <- as.matrix(masked_lake_interpolation)
wireframe(mat_lake_interpolation,
shade = TRUE,
pretty = TRUE,
xlab = "",
ylab = "",
zlab = "",
main = "Bathymetric model of Lake Rangsdorf",
screen = list(z = -90, x = -55),
scales = list(z = list(arrows = FALSE)),
aspect = c(1, 0.25),
zoom = 1.3,
par.settings = list(
axis.line = list(col = "transparent"),
clip = list(panel = "off")
)
)
Imputation of missing depth values
Now that we have a bathymetric model, we us it to impute missing
depth values. If we print the depth values in our data set we see that
exactly one of them is missing, indicated by NA.
lake_data_sf$depth_m## [1] 0.68 0.66 1.93 3.98 4.95 0.93 3.92 8.40 5.15 4.00 2.97 3.64 2.50 1.90 2.38
## [16] 2.30 2.12 1.38 1.70 1.85 1.68 1.57 1.38 1.08 NA 2.03 1.53 0.50 0.68 0.55
## [31] 0.92 0.93First, we localize the row of the missing value.
idx <- which(is.na(lake_data_sf$depth_m))
idx## [1] 25lake_data_sf[idx, ]## Simple feature collection with 1 feature and 21 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 390807 ymin: 5794369 xmax: 390807 ymax: 5794369
## Projected CRS: WGS 84 / UTM zone 33N
## ID depth_m visual_depth_cm elec_cond_muS pH LOI550 LOI990 LOI C N
## 25 38a NA 38 709 7.04 22.33 26.53 17.37 18.6 1.19
## Ca_mg_g Fe_mg_g K_mg_g Pb_mug_g Mg_mg_g Mn_mug_g Na_mg_g PO4_mg_g S_mg_g
## 25 308.04 6.17 0.27 30.78 1.36 498.38 0.5 3.75 14.98
## Sr_mug_g Zi_mug_g geometry
## 25 370.78 90.46 POINT (390807 5794369)Then, we apply the extract() function from the
raster package to extract the depth of the missing point
value from the bathymetric model and assign the result to the missing
entry in the data set. Please note that the extract()
function expects an sp object, thus we cast the
sf object into an sp object by applying the
as() function. Further note that in the bathymetric model
the depth is given as negative number, whereas in the
lake_data_sf the depth is given as positive numbers. Hence,
we multiply the depth from the bathymetric model with \(-1\).
# define point of interest
poi <- lake_data_sf[idx, "depth_m"]$geometry
# extract depth value
lake_data_sf[idx, "depth_m"] <- extract(
masked_lake_interpolation,
as(poi, "Spatial")
) * -1Now the missing value is imputed!
lake_data_sf$depth_m## [1] 0.68000 0.66000 1.93000 3.98000 4.95000 0.93000 3.92000 8.40000 5.15000
## [10] 4.00000 2.97000 3.64000 2.50000 1.90000 2.38000 2.30000 2.12000 1.38000
## [19] 1.70000 1.85000 1.68000 1.57000 1.38000 1.08000 1.54654 2.03000 1.53000
## [28] 0.50000 0.68000 0.55000 0.92000 0.93000Lake morphometry
Using the lake polygon we can compute the area of lake
Rangsdorf by applying the st_area() function from the
sf package.
st_area(lake)## 2448172 [m^2]Based on the bathymetric model we computed in the previous section,
we may as well calculate the mean and maximum depth using the
cellStats() function from the raster
package.
mean_depth <- cellStats(masked_lake_interpolation, "mean") * -1
mean_depth## [1] 1.393101max_depth <- cellStats(masked_lake_interpolation, "min") * -1
max_depth## [1] 9.604918In addition, we may look at the histogram of depth values by applying
the values() function and plot them with
ggplot.
library(tidyverse)
depth <- values(masked_lake_interpolation * -1)
# exclude invalid data
depth <- depth[complete.cases(depth)]
depth <- depth[depth != 0]
# make a data.frame object for easy plotting
df_depth <- data.frame(depth = depth)
ggplot(data = df_depth, aes(x = depth)) +
geom_histogram(stat = "bin", binwidth = 0.1) +
## vertical line for mean depth
geom_vline(
xintercept = mean_depth,
color = "blue", lwd = 1.25
) +
annotate("text",
x = 0.3, y = 1700,
label = "mean depth", color = "blue", cex = 3
) +
annotate("text",
x = 1.15, y = 1710,
label = sprintf("\u2192"), color = "blue"
) +
## vertical line for max depth
geom_vline(xintercept = max_depth, color = "red", lwd = 1.25) +
annotate("text",
x = 8.45, y = 1700,
label = "max depth", color = "red", cex = 3
) +
annotate("text",
x = 9.25, y = 1710,
label = sprintf("\u2192"), color = "red"
) +
ggtitle("Depth Histogram Lake Rangsdorf")
Beside the mean and maximum depth we may also calculate the volume of lake Rangsdorf.
Exercise: Calculate the volume of Lake Rangsdorf. This can be done by calculating the sum of the depth of all grid cells of the
masked_lake_interpolationobject and then multiplying the sum with the area of the grid cell. Applying theres()function will give the resolution of therasterobject.
## Your code here...
depth_sum <- NULL
volume <- NULLShow code
# resolution of the interpolation grid
res(masked_lake_interpolation)## [1] 10 10depth_sum <- cellStats(masked_lake_interpolation, "sum") * -1
volume <- res(masked_lake_interpolation)[1] * res(masked_lake_interpolation)[2] * depth_sum
volume## [1] 3410730The volume of lake Rangsdorf is approximately \(3,411,000\, m^3\).
Minimum distance to shore
An interesting property of each sample point location is its distance
to the shoreline. We apply the st_distance() function to
compute the distance of each sample point to shoreline and append the
result to the data set.
lake_data_sf$dist <- as.numeric(st_distance(x = lake_data_sf, y = shoreline))
lake_data_sf$dist## [1] 27.536302 6.854090 10.094554 68.471164 93.630593 5.510838
## [7] 35.082980 81.272256 85.280820 68.260142 70.608362 57.225353
## [13] 365.238433 185.789339 169.247017 296.955821 528.424971 5.721908
## [19] 105.451041 162.692419 260.329215 78.925367 149.579911 261.430696
## [25] 254.903516 415.961256 16.745054 101.917240 126.494442 47.202636
## [31] 37.807563 22.524395For later usage we write the several objects to disk.
save(
file = "Lake_Rangsdorf.RData",
lake_data_sf, lake_interpolation, lake, shoreline
)Citation
The E-Learning project SOGA-R was developed at the Department of Earth Sciences by Kai Hartmann, Joachim Krois and Annette Rudolph. You can reach us via mail by soga[at]zedat.fu-berlin.de.
You may use this project freely under the Creative Commons Attribution-ShareAlike 4.0 International License.
Please cite as follow: Hartmann, K., Krois, J., Rudolph, A. (2023): Statistics and Geodata Analysis using R (SOGA-R). Department of Earth Sciences, Freie Universitaet Berlin.