R interface for TELEMAC-2D
Tobias Pilz
February 2021
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.
Prerequisites
Aside from 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.
In addition we need some more R packages to run this demonstration, which include tidyverse, raster, and sf. The packages can be installed directly from CRAN.
library(raster)
library(sf)
library(tidyverse)
library(telemac)Input data
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 (mesh) file (*.slf)
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 sf (telemac can handle both).
bnd <- st_read(system.file("dem/boundary_lagos.gpkg", package = "telemac"))
## Reading layer `aoi_line_small' from data source
## `/tmp/Rtmp09aQS5/Rinst228e8cfa7f2/telemac/dem/boundary_lagos.gpkg'
## using driver `GPKG'
## Simple feature collection with 1 feature and 1 field
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 530280 ymin: 729900 xmax: 536310 ymax: 744300
## Projected CRS: WGS 84 / UTM zone 31NThe 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 t2d_tin objects.
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 = ".")
What is still missing to obtain a complete mesh for application within TELEMAC is elevation for each TIN point. This can be inferred by interpolation from the DEM.
For this step the telemac package provides function geo(). The function accepts several input types, e.g. the name of an existing SELAFIN file from which the geometry is read. However, in this demonstration we want to generate a new geometry file based on our triangulated DEM. The interpolation will be conducted using an approach of inverse distance weighting (IDW) of nearest neighbours as implemented in function idw() of package gstat.
In the example, IDW is conducted using the five (n = 5) nearest neighbours around each TIN point and a weighting power of 2 (idp = 2) of inverse distances.
geo_obj <- geo(tin_obj, dem = dem_rast, title = "title", fname = "geometry.slf",
n = 5, idp = 2)The resulting t2d_geo object is a list and individual elements may be accessed (or modified) with list functionalities.
str(geo_obj)
## List of 4
## $ header :List of 5
## ..$ title : chr "title"
## ..$ precision: int 4
## ..$ nelem : int 7991
## ..$ npoin : int 4257
## ..$ ndp : int 3
## $ tin :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"
## $ elevation: num [1:4257] 45.6 45.5 44.8 43.6 43.1 ...
## $ privars : NULL
## - attr(*, "file")= chr "geometry.slf"
## - attr(*, "class")= chr [1:2] "t2d_geo" "list"
geo_obj$header$title
## [1] "title"
geo_obj$header$title <- "Merit, RTriangle, 90m"
geo_obj
## Object of class t2d_geo: TELEMAC geometry (mesh)
## Mesh title: Merit, RTriangle, 90m
## The mesh is composed of 7991 elements (triangles) and 4257 points
## Associated geometry file: geometry.slfThere is also a plot method implemented. This creates a 2d surface plot of the mesh by linearly interpolating the three values of a triangle to grid points of given resolution (s = 30 meter) within the triangle. The higher the resolution (the smaller argument s), the smoother the plot will appear.
plot(geo_obj, s = 30, col = colorRampPalette(c("blue", "green", "red"))(255))
The boundary conditions file (*.cli)
The boundary conditions are stored as a simple text file that can be investigated and modified with any text editor. Each line in the file represents a point along the mesh boundary. Function cli() generates an object of class t2d_cli that is associated with a boundary conditions file. The function either imports an existing file or creates a new configuration when providing a data.frame or a t2d_geo object from which a template configuration will be inferred. The template configuration represents a closed boundary with zero prescribed depths or velocities.
cli_obj <- cli(geo_obj, fname = "boundary.cli")
cli_obj
## Object of class t2d_cli: TELEMAC boundary parameters
## Associated boundary file: boundary.cli
## First 10 of 521 boundary points:
## lihbor liubor livbor hbor ubor vbor aubor litbor tbor atbor btbor n k
## 1 2 2 2 0 0 0 0 2 0 0 0 21 1
## 2 2 2 2 0 0 0 0 2 0 0 0 20 2
## 3 2 2 2 0 0 0 0 2 0 0 0 19 3
## 4 2 2 2 0 0 0 0 2 0 0 0 18 4
## 5 2 2 2 0 0 0 0 2 0 0 0 17 5
## 6 2 2 2 0 0 0 0 2 0 0 0 16 6
## 7 2 2 2 0 0 0 0 2 0 0 0 15 7
## 8 2 2 2 0 0 0 0 2 0 0 0 14 8
## 9 2 2 2 0 0 0 0 2 0 0 0 13 9
## 10 2 2 2 0 0 0 0 2 0 0 0 12 10Internally a t2d_cli object is a data.frame and individual element may be accessed (or modified) with data.frame (in this example tidyverse) functionalities.
str(cli_obj)
## Classes 't2d_cli' and 'data.frame': 521 obs. of 13 variables:
## $ lihbor: num 2 2 2 2 2 2 2 2 2 2 ...
## $ liubor: num 2 2 2 2 2 2 2 2 2 2 ...
## $ livbor: num 2 2 2 2 2 2 2 2 2 2 ...
## $ hbor : num 0 0 0 0 0 0 0 0 0 0 ...
## $ ubor : num 0 0 0 0 0 0 0 0 0 0 ...
## $ vbor : num 0 0 0 0 0 0 0 0 0 0 ...
## $ aubor : num 0 0 0 0 0 0 0 0 0 0 ...
## $ litbor: num 2 2 2 2 2 2 2 2 2 2 ...
## $ tbor : num 0 0 0 0 0 0 0 0 0 0 ...
## $ atbor : num 0 0 0 0 0 0 0 0 0 0 ...
## $ btbor : num 0 0 0 0 0 0 0 0 0 0 ...
## $ n : int 21 20 19 18 17 16 15 14 13 12 ...
## $ k : int 1 2 3 4 5 6 7 8 9 10 ...
## - attr(*, "file")= chr "boundary.cli"
cli_obj <- cli_obj %>%
# open boundary with free depth and velocities
mutate(lihbor = 4, liubor = 4, livbor = 4)The steering file (*.cas)
The steering file is a simple text file that contains a number of key–value pairs (delimited by characters : or =), the steering parameters for a TELEMAC simulation. The telemac package comes with function cas() to initialise a t2d_cas object with steering parameters associated with a steering file. The function provides different options, e.g. providing the name of an existing file or creating a template configuration if no argument is given. The template configuration contains a number of arbitrary parameters, in this case defining a strong rainfall event of 66 mm over 4 hours, dry initial conditions, a total simulation period of one day with adaptive temporal resolution using a finite-volume scheme, and output printout every 3600 timesteps. Note that there is a huge number of possible steering parameters that are not included in the template file (for these the respective default settings are used by TELEMAC).
cas_obj <- cas(fname = "steering.cas")
cas_obj
## Object of class t2d_cas: TELEMAC steering parameters
## Associated steering file: steering.cas
## First 10 of 20 parameters:
## key value
## 1 BOUNDARY CONDITIONS FILE boundary.cli
## 2 GEOMETRY FILE geometry.slf
## 3 RESULTS FILE results.slf
## 4 TITLE 'template'
## 5 VARIABLES FOR GRAPHIC PRINTOUTS 'U,V,S,B,H,M'
## 6 DURATION 86400
## 7 GRAPHIC PRINTOUT PERIOD 3600
## 8 LISTING PRINTOUT PERIOD 3600
## 9 MASS-BALANCE YES
## 10 COMPUTATION CONTINUED NOInternally a t2d_cas object is a list and individual element may be accessed (or modified) with list functionalities.
str(cas_obj)
## List of 20
## $ BOUNDARY CONDITIONS FILE : chr "boundary.cli"
## $ GEOMETRY FILE : chr "geometry.slf"
## $ RESULTS FILE : chr "results.slf"
## $ TITLE : chr "'template'"
## $ VARIABLES FOR GRAPHIC PRINTOUTS : chr "'U,V,S,B,H,M'"
## $ DURATION : chr "86400"
## $ GRAPHIC PRINTOUT PERIOD : chr "3600"
## $ LISTING PRINTOUT PERIOD : chr "3600"
## $ MASS-BALANCE : chr "YES"
## $ COMPUTATION CONTINUED : chr "NO"
## $ INITIAL CONDITIONS : chr "'ZERO DEPTH'"
## $ EQUATIONS : chr "'SAINT-VENANT FV'"
## $ FINITE VOLUME SCHEME : chr "5"
## $ VARIABLE TIME-STEP : chr "YES"
## $ DESIRED COURANT NUMBER : chr "0.9"
## $ LAW OF BOTTOM FRICTION : chr "3"
## $ FRICTION COEFFICIENT : chr "50"
## $ RAIN OR EVAPORATION : chr "YES"
## $ RAIN OR EVAPORATION IN MM PER DAY : chr "400.0"
## $ DURATION OF RAIN OR EVAPORATION IN HOURS: chr "4."
## - attr(*, "file")= chr "steering.cas"
## - attr(*, "class")= chr [1:2] "t2d_cas" "list"
cas_obj[["VARIABLES FOR GRAPHIC PRINTOUTS"]] <- "H"TELEMAC-2D project setup
Finally we can complete our TELEMAC-2D setup. To do so we create a t2d object using function t2d(). Adjust argument wdir (path to working, i.e. TELEMAC-2D project directory) according to your needs.
t2d_obj <- t2d(title = "Test setup", wdir = "t2d_test",
cas = cas_obj, geo = geo_obj, cli = cli_obj)
t2d_obj
## Object of class t2d: a TELEMAC-2D project
##
## Title: Test setup
## Project directory: t2d_test
## TELEMAC-2D executable: <unspecified>
## Steering parameters: A t2d_cas object pointing to steering.cas
## Geometry / mesh: A t2d_geo object pointing to geometry.slf
## Boundary conditions: A t2d_cli object pointing to boundary.cli
## Optional input: None given
## Simulation results: A t2d_res object pointing to results.slf (file does not yet exist)Note that we have not yet written any files. As a last step to prepare a TELEMAC-2D run we need to write the model’s input files. All prepared input data will be written into their associated files relative to wdir given in t2d_obj.
write_t2d(t2d_obj)Model simulation
The model can be run either directly or via R interface simulate_t2d(). However, TELEMAC-2D runs, depending on domain size, spatial and temporal resolution, and numerical settings, can take quite a while. Using the R interface to conduct simulations makes only sense if you are sure the model run will not take too long, e.g. for test runs of minimal examples.
The simulate_t2d() function accepts a selection of output variables to be imported into R via argument vars. Depending on the number of mesh points, output intervals, and variables specified in the steering file the required storage capacity in R can be high. Therefore it makes sense to only import variables for immediate analyses.
t2d_obj <- simulate_t2d(t2d_obj, log = "test_run.log", res = "test_results.slf",
vars = "water depth", exec = "telemac2d.py")Note that the results of a previous simulation run can be imported as well. Besides, we import only some selected timesteps for illustration.
t2d_obj$res <- results(system.file("telemac/basics/results.slf", package = "telemac"),
log = system.file("telemac/basics/test_run.log", package = "telemac"),
vars = "water depth", times = 3600 * c(0,1,4,12,24))The simulation results are represented as a t2d_res object with header information, the underlying mesh as t2d_tin object, a log of TELEMAC’s runtime messages, and the actual data as a tidy data.frame. For a first investigation it is always good to have a look at the log file. The runtime log is stored in the t2d_res object and is accessible via t2d_obj$res$log (not shown here). Otherwise there is a print method to show some general information.
t2d_obj$res
## Object of class t2d_res: TELEMAC simulation results
## Associated results file: /tmp/Rtmp09aQS5/Rinst228e8cfa7f2/telemac/telemac/basics/results.slf
## Associated log file: /tmp/Rtmp09aQS5/Rinst228e8cfa7f2/telemac/telemac/basics/test_run.log
## Simulation title: Test setup
## The results comprise the following 1 variables (units):
## WATER DEPTH (M)
## over 25 timesteps.
## Imported 5 timesteps ( 0, 3600, 14400, 43200, 86400 ) and variable water depth as a tidy data.frame.
## First 10 of 21285 values:
## variable timestep x y value
## 1 water depth 0 531990.0 729900.0 0
## 2 water depth 0 531909.5 729940.2 0
## 3 water depth 0 531843.0 730000.9 0
## 4 water depth 0 531780.4 730065.6 0
## 5 water depth 0 531690.5 730070.1 0
## 6 water depth 0 531600.6 730074.6 0
## 7 water depth 0 531510.8 730079.2 0
## 8 water depth 0 531438.9 730133.4 0
## 9 water depth 0 531451.1 730222.6 0
## 10 water depth 0 531381.2 730279.4 0
str(t2d_obj$res)
## List of 4
## $ header:List of 9
## ..$ title : chr "Test setup"
## ..$ precision: int 4
## ..$ nelem : int 7991
## ..$ npoin : int 4257
## ..$ ndp : int 3
## ..$ varnames : chr "WATER DEPTH"
## ..$ varunits : chr "M"
## ..$ date : POSIXct[1:1], format: "1900-02-01"
## ..$ ntimes : num 25
## $ tin :List of 5
## ..$ points : num [1:4257, 1:2] 531990 531910 531843 531780 531690 ...
## ..$ triangles : int [1:7991, 1:3] 526 22 878 23 522 630 22 1270 2085 25 ...
## ..$ edges : int [1:12247, 1:2] 526 911 526 21 19 19 751 618 618 23 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "V1" "V2"
## ..$ boundaries: int [1:521] 21 20 19 18 17 16 15 14 13 12 ...
## ..$ breaklines: NULL
## ..- attr(*, "class")= chr [1:2] "t2d_tin" "list"
## $ log : chr [1:248] "" "" "Loading Options and Configurations" "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" ...
## $ values:'data.frame': 21285 obs. of 5 variables:
## ..$ variable: chr [1:21285] "water depth" "water depth" "water depth" "water depth" ...
## ..$ timestep: num [1:21285] 0 0 0 0 0 0 0 0 0 0 ...
## ..$ x : num [1:21285] 531990 531910 531843 531780 531690 ...
## ..$ y : num [1:21285] 729900 729940 730001 730066 730070 ...
## ..$ value : num [1:21285] 0 0 0 0 0 0 0 0 0 0 ...
## - attr(*, "file")= chr "/tmp/Rtmp09aQS5/Rinst228e8cfa7f2/telemac/telemac/basics/results.slf"
## - attr(*, "log")= chr "/tmp/Rtmp09aQS5/Rinst228e8cfa7f2/telemac/telemac/basics/test_run.log"
## - attr(*, "class")= chr [1:2] "t2d_res" "list"There is also a function write_results() to write data into a SELAFIN file. This can be used, for instance, to reduce the amount of data or write retrospectively corrected data. Note that in this example we have only imported a few selected timesteps that will now be written into a SELAFIN file.
# change associated file (otherwise the original file would be overwritten)
t2d_obj$res <- results(t2d_obj$res, fname = "t2d_test/results_updated.slf")
t2d_obj$res <- write_results(t2d_obj)
t2d_obj$resFor visual inspection a plot method is available. In the following variable v = "water_depth" at timestep t = 4*3600 is shown. As for the mesh explained above, the values are interpolated to a grid with s = 30 meter resolution to generate a 2d surface plot.
plot(t2d_obj$res, s = 30, v = "water depth", t = 4*3600,
col = c('#eff3ff','#c6dbef','#9ecae1','#6baed6','#3182bd','#08519c'))