Merge Bathymetry with DEM data
merge-bathy-dem.RmdLake bathymetry data is useful for a range of applications, including habitat mapping, water quality monitoring, and hydrodynamic modelling. However, lakes sit within a landscape, and it is often useful to know how the bathymetry data relates to the surrounding topography. This is particularly important for hydrodynamic modelling of lakes with large fluctuations in water level, or for flood risk assessments in the surrounding area.
This vignette demonstrates how to merge bathymetry data with a
Digital Elevation Model (DEM) raster using the
merge_bathy_dem() function.
Load the data
We will use the bathytools package to merge bathymetry data with a DEM raster. The package includes example data for Lake Rotoma, New Zealand. The bathymetry data is stored as a XYZ data.frame in an .rds file within the package. Whereas the DEM data is stored as a raster in a TIF file. The DEM raster was prepared from LiDAR data provided by Land Information New Zealand and can be downloaded here.
We will also use a shapefile of the lake shoreline and a shapefile of the lake catchment. These were both sourced from the Freshwater Ecosytems of New Zealand database.
library(bathytools)
library(tmap) # Used for plotting spatial data
#> Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
#> remotes::install_github('r-tmap/tmap')
shoreline <- readRDS(system.file("extdata/rotoma_shoreline.rds",
package = "bathytools"))
catchment <- readRDS(system.file("extdata/rotoma_catchment.rds",
package = "bathytools"))The example lake we will be looking at is Lake Rotoma in the Bay of Plenty Region on the North Island of Aotearoa New Zealand. The lake is a popular recreational spot for fishing and boating.
The tmap package is used
to plot spatial data. The tmap_mode("view") function is
used to display the map in the viewer pane in RStudio. The
tmap_options(basemaps = "Esri.WorldImagery") function is
used to set the basemap to a satellite image. It is similar to the
ggplot2 package, but is specifically designed for spatial
data. Here we plot the lake shoreline in light blue and the lake
catchment in pink.
tmap_mode("view")
#> tmap mode set to interactive viewing
tmap_options(basemaps = "Esri.WorldImagery")
tm_shape(shoreline) +
tm_borders(col = "#8DA0CB", lwd = 2) +
tm_shape(catchment) +
tm_borders(col = "#E78AC3", lwd = 2) The depth data is stored in XYZ format, with the x and y coordinates representing the location of the depth point, and the z coordinate representing the depth in meters. The depth data is stored as a data.frame in an .rds file within the package. Here we show the first few rows of the depth data.
point_data <- readRDS(system.file("extdata/depth_points.rds",
package = "bathytools"))
head(point_data)
#> lon lat depth
#> 1 176.5952 -38.03974 -3.88
#> 2 176.5853 -38.03590 -79.90
#> 3 176.5904 -38.03109 -63.30
#> 4 176.5960 -38.06580 -5.06
#> 5 176.5790 -38.03500 -38.51
#> 6 176.5691 -38.02957 -20.06
dem_raster <- terra::rast(system.file("extdata/dem_32m.tif",
package = "bathytools"))
tm_shape(dem_raster) +
tm_raster(alpha = 0.5, style = "cont", palette = "-YlGnBu") +
tm_shape(shoreline) +
tm_borders(col = "#FC8D62", lwd = 2) +
tm_shape(catchment) +
tm_borders(col = "#A6D854", lwd = 2) Generate the bathymetry raster
The first step is to generate a bathymetry raster from the shoreline
and depth data. The rasterise_bathy() function is used to
generate the bathymetry raster. The function takes the shoreline, depth
data, and the coordinate reference system (CRS) as inputs. The function
returns a SpatRaster object representing the bathymetry
raster.
bathy_raster <- rasterise_bathy(shoreline = shoreline,
point_data = point_data, crs = 2193,
res = 8)
#> Generating depth points... [2024-11-13 23:09:49]
#> Finished! [2024-11-13 23:09:49]
#> Interpolating to raster... [2024-11-13 23:09:49]
#> Adjusting depths >= 0 to -0.82 m
#> Finished! [2024-11-13 23:09:50]
tm_shape(dem_raster) +
tm_raster(alpha = 0.5, style = "cont", palette = "-YlGnBu") +
tm_shape(bathy_raster) +
tm_raster(style = "cont", palette = "-viridis", breaks = seq(-90, 0, by = 10)) +
tm_shape(shoreline) +
tm_borders(col = "#FC8D62", lwd = 2) +
tm_shape(catchment) +
tm_borders(col = "#A6D854", lwd = 2) Merge the bathymetry with the DEM
The next step is to merge the bathymetry raster with the DEM raster.
The merge_bathy_dem() function is used to merge the
bathymetry raster with the DEM raster. The function takes the shoreline,
bathymetry raster, DEM raster, and catchment shapefile as inputs. The
function returns a SpatRaster object representing the
merged bathymetry and DEM data.
If the resolution of the bathymetry raster is different from the DEM raster, the bathymetry raster will be resampled to match the resolution of the DEM raster.
dem_bath <- merge_bathy_dem(shoreline = shoreline, bathy_raster = bathy_raster,
dem_raster = dem_raster, catchment = catchment)
#> Warning: Resolutions of bathy_raster [8, 8] and dem_raster [32, 32] do notmatch. Resampling bathy_raster to match dem_raster
#> Lake surface elevation from DEM: 313.3 m
dem_bath
#> class : SpatRaster
#> dimensions : 217, 194, 1 (nrow, ncol, nlyr)
#> resolution : 32, 32 (x, y)
#> extent : 1911744, 1917952, 5779472, 5786416 (xmin, xmax, ymin, ymax)
#> coord. ref. : NZGD2000 / New Zealand Transverse Mercator 2000 (EPSG:2193)
#> source(s) : memory
#> varname : dem_32m
#> name : elevation
#> min value : 232.0779
#> max value : 599.1376Spatial plot of the merged raster
Here we plot the merged raster with the lake and catchment boundaries.
tm_shape(dem_bath) +
tm_raster(alpha = 0.5, style = "cont", palette = "-YlGnBu") +
tm_shape(shoreline) +
tm_borders(col = "#FC8D62", lwd = 2) +
tm_shape(catchment) +
tm_borders(col = "#A6D854", lwd = 2) The colours in this plot do not clearly distinguish between the
bathymetry and DEM data. We will add a break at the surface elevation of
the lake to better distinguish between the two datasets. We can extract
the water surface elevation from the DEM data using the
get_lake_surface_elevation() function.
lake_elev <- get_lake_surface_elevation(dem_raster = dem_raster,
shoreline = shoreline)
#> Lake surface elevation from DEM: 313.3 mWe can now plot the merged raster with the lake surface elevation as the break in the colour palette.
tm_dem_bath(dem_bath = dem_bath, lake_elev = lake_elev)Hypsograph of the merged raster
A hypsograph is a plot of the area of a lake at different depths. The
bathy_to_hypso() function can be used to generate a
hypsograph from the bathymetry raster. The function takes the bathymetry
raster as input and returns a data frame with the depth and area at each
depth.
hyps <- bathy_to_hypso(bathy_raster = bathy_raster)
head(hyps)
#> depth area
#> 1 0 11240576
#> 2 -1 10494272
#> 3 -2 10215744
#> 4 -3 9997504
#> 5 -4 9797120
#> 6 -5 9604800The hypsograph can be plotted to show the area of the lake at different depths.
library(ggplot2)
ggplot(hyps, aes(x = area, y = depth)) +
geom_line() +
geom_point() +
labs(x = "Area (m^2)", y = "Depth (m)") +
theme_bw()
3-D plot of the merged raster
The plot_raster_3d() function can be used to create a
3-D plot of the merged raster. The function takes the merged raster and
the shoreline as inputs. The fact argument controls the
aggregation factor for the raster. A higher factor will result in a
smoother plot, but will take longer to render.
p1 <- plot_raster_3d(x = dem_bath, shoreline = shoreline, split_lake = TRUE)
p1Saving the merged raster
The merged raster can be saved to a file using the
terra::writeRaster() function. The function takes the
raster object and the file path as inputs.
It is important to note that SpatRaster can not be saved
as “.rds” files. They can also be quite large, so it is recommended to
save the raster in a compressed format, such as GeoTIFF.
terra::writeRaster(dem_bath, "dem_bath.tif", overwrite = TRUE)