This document briefly demonstrates the basic functionalities of the
telemac package. The package is an R interface for the modelling suite OpenTELEMAC. The focus of this vignette is the TELEMAC-2D module for 2-dimensional hydrodynamic modelling. The example demonstrates rainfall-induced inundation modelling for an urban area in the northwest of the megacity of Lagos, Nigeria, including mesh generation based on a freely available DEM.
R and the
telemac package (and its package dependencies) the OpenTELEMAC software is needed. The model suite can be obtained and downloaded for free after registration from http://www.opentelemac.org. Install it and make sure it is running on your system.
library(raster) library(sf) library(tidyverse) library(telemac)
For a basic TELEMAC-2D setup three input files are required: a geometry file (
*.slf) with mesh information, boundary conditions (
*.cli), and steering parameters (
*.cas). In the following the basic setup of each input is explained in more detail.
The geometry of a TELEMAC setup is stored in a binary file of type SELAFIN (also referred to as SERAFIN). The mesh is typically stored as a triangulated irregular network (TIN). There are several possible approaches to obtain such a mesh, which shall not be discussed here in greater detail. In this demonstration we will start with a typical DEM represented by a regular grid and convert it into a TIN by constrained Delaunay triangulation using the well-known Triangle algorithm.
First we read an example DEM. This shows an urban area in the northwestern region of the city of Lagos, Nigeria, and is derived from the (after registration) freely available MERIT Hydro DEM. Despite its coarse resolution of about 90 x 90 m the DEM is hydrologically conditioned and is known for a relatively good performance in hydrodynamic modelling in comparison to other freely available DEM products. However, in real-world applications, high resolution LiDAR-based DEMs are certainly required to obtain optimal model results.
dem_rast <- raster(system.file("dem/dem_merit_lagos.tif", package = "telemac")) NAvalue(dem_rast) <- -999 # I defined -999 to be NA when creating the dem file plot(dem_rast, col = colorRampPalette(c("blue", "green", "red"))(255))
In addition we need the catchment boundary. Vector data in R are best handled using package
sp or the successor
telemac can handle both).
bnd <- st_read(system.file("dem/boundary_lagos.gpkg", package = "telemac")) ## Reading layer `aoi_line_small' from data source `/tmp/RtmpEAitar/Rinst181ff82b192/telemac/dem/boundary_lagos.gpkg' using driver `GPKG' ## Simple feature collection with 1 feature and 1 field ## geometry type: LINESTRING ## dimension: XY ## bbox: xmin: 530280 ymin: 729900 xmax: 536310 ymax: 744300 ## projected CRS: WGS 84 / UTM zone 31N
telemac package provides the function
tin() to initialise a TIN the model mesh will be based on. Internally it uses the R package RTriangle to perform triangulation with the Triangle algorithm. In the example the catchment boundary is taken and triangulation performed starting with a resolution of 90 m (
s = 90) of TIN points along the catchment boundary. Furthermore, the resulting triangles shall not be larger than 10,000 m^2 (
a = 100^2) and the angles within the triangles not smaller than 30 degrees (
q = 30).
There is also the possibility to include breaklines to further refine the triangulation (see
?tin). Holes are, however, not yet supported. Besides, there are no checks implemented to assure the success of triangulation and validity of the resulting TIN for model application (e.g. all breaklines must be within the boundary, breaklines should not intersect, etc.). A good summary about successful mesh generation for TELEMAC is presented in the documentation of the pputils python tools.
tin_obj <- tin(list(boundary = bnd), s = 90, a = 100^2, q = 30)
The function creates an object of class
t2d_tin that contains all TIN points and their relations to define the triangles.
str(tin_obj) ## List of 5 ## $ points : num [1:4257, 1:2] 531990 531910 531843 531780 531691 ... ## $ triangles : int [1:7991, 1:3] 526 22 878 23 522 630 22 1270 2085 25 ... ## $ edges : int [1:12247, 1:2] 526 1153 911 22 21 19 878 751 618 23 ... ## $ boundaries: int [1:521] 21 20 19 18 17 16 15 14 13 12 ... ## $ breaklines: NULL ## - attr(*, "class")= chr [1:2] "t2d_tin" "list"
There are also a print and a plot method for
tin_obj ## Object of class t2d_tin: TELEMAC mesh (TIN) ## The mesh is composed of 7991 elements (triangles) and 4257 points. plot(tin_obj, pch = ".")