Postprocessing TELEMAC-2D results with Python

A typical postprocessing task is the visualization of time dependent water levels e.g. along the river axis.
Such visualizations are not beautiful only but they help to analyze the results or to check the simulation settings.

This post is intented to be a short tutorial how to interpolate the time dependent water levels along a polyline from a Telemac-2D result file and to produce a nice video. We will not go into details about the basics of Python, Matplotlib and Numpy scripting but on how to produce the video inside the Python framework. Such scripts are very cool since, once the setup has been done, it can be quickly used for other cases, too.

For example the following animation shows the evolution of the time dependent water level and erodible bottom in a longitudinal section of a simulated dam failure with Telemac-2D and Sisyphe (for better video quality watch in full screen mode!)


What you need:
Software:
– Python 2.7.x with site-packages Matplotlib and Numpy
– FFmpeg, open source video and audio converter, https://www.ffmpeg.org
– Python scripts, which are usually equipped within an openTelemac installation: ..\vXX\scripts\python27

Data:
– Telemac-2D result file
– Polyline with x and y coordinates, e.g. generated with Blue Kenue

Let’s say we have a Telemac-2D result file of a simulated river which by chance has the shape of flussbüro. The flood propagation on an initially dry bed caused by a sudden release of water from a reservoir was simulated. Now we want to analyze the time dependent water levels along the red polyline, as shown in the following figure:

With the Python script explained below we want to produce such a video:


Some preliminary remarks:
– some Python skills are assumed and hence, the script will not be explained line by line
– the script seems really long but it is equipped with a lot of comments
– most of the lines deal with formatting the layout of the figure / video
– it is advisable to write the defined functions explained below in a separate file and import them
– copy the code below to an empty text file and open it with a Python editor

From line 5 to line 13
Python modules and site-packages are imported. In line 7 we add the Python scripts directory of the Telemac Python scripts in order to import the selafin parser in line 13.

From line 17 to line 23
Telemac result file and the polyline are imported. E.g. here the polyline consisting of x an y coordinates was generated with Blue Kenue and saved as i2s file.

From line 41 to line 49
The function “Get_polyline_distances” calculates the stationing (distances) and the transposed input polyline (poly_new)

From line 52 to line 61
The function “Get_Values” gets the time dependent results of a specific variable and puts them in an array

From line 64 to line 81
The function “Interpolate_1D_profiles” interpolates the time dependent results of a specific variable on the polyline

Line 86: now we calculate with the above defined function “Get_polyline_distances” the distances and the transposed polyline of our input polyline

From line 89 to line 102
First, here we can choose which time steps (frames) we want to extract and visualize. Second, in this example we interpolate the Bottom and the Free Surface to the polyline

From line 105 to line 148
Here we specify the layout options of the figure / video

Lines 150 to 181 finally deal with the production of the video:
In line 156 we specify the path to the installed ffmpeg executable
In line 160 we specify the amazingly good libx264 codec and a constant rate factor of 18 for the video quality. One can set this video quality factor to the default value of 23.
From line 168 to 171 we specify e.g. the line colors, line widths etc.
In line 176 we put the interpolated free surface values in a loop.

My work is done!

#-------------------------------------------------------------------------------
# Author:      Clemens Dorfmann - www.flussbuero.at
# Created:     15/11/2017
#-------------------------------------------------------------------------------
import sys
from os import path
sys.path.append(path.join( path.dirname(sys.argv[0]), r'..\path_to\scripts\python27'))
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import tri
import matplotlib.animation as animation

from parsers.parserSELAFIN import SELAFIN
# from myFunctions import Get_polyline_distances, Get_Values, Interpolate_1D_profiles

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
INPUT / OUTPUT Folder and TELEMAC RESULT FILE
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
InputFolder = r"..\path to Telemac-2D result file and polyline file/"
OutputFolder = r"..\path to your output directory" # output of videos, pictures, etc.

Input_slf = "Telemac-2D_results_file.slf"
Input_polyline = np.genfromtxt(InputFolder + 'polyline.i2s', delimiter=' ', usecols = (0, 1), skip_header=16)

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
LOAD TELEMAC FILE
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
filename_in = (InputFolder + Input_slf)
filename_out = (OutputFolder + Input_slf)
slf = SELAFIN(filename_in)

times = slf.tags["times"]
#print len(times)

var_names=[item.rstrip() for item in slf.VARNAMES]
print var_names

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
FUNCTIONS
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
# Function to calculate distances between points of a polyline
def Get_polyline_distances(polyline):
    poly_new = np.transpose(polyline)
    delta_x_y = np.diff(poly_new)
    segdistances = np.sqrt((delta_x_y**2).sum(axis=0))
    segdistances = np.insert(segdistances,0,0.0) # insert at i=0 -> 0.0
    distances = np.cumsum(segdistances,axis=0)

    return poly_new, distances

# Function to get values in time and space of a specific variable
def Get_Values(slf, var_name, frames):
    var_names=[item.rstrip() for item in slf.VARNAMES]
    z_var = []

    for i in frames:
        z_var.append(slf.getVariablesAt(i, [var_names.index(var_name)]))

    z_var = np.squeeze(z_var, axis=None) # remove one-dimensional entries

    return z_var

# Function for interpolating values to polyline vertices
def Interpolate_1D_profiles(slf, x, y, mesh, z_var, frames, poly_new):

    # Interpolate using tri.LinearTriInterpolator
    triang = tri.Triangulation(x, y, mesh) # triangulation

    if len(frames) > 1: # extract values for more than one frame
        zp = []
        for i in range(len(frames)):
            interpolator = tri.LinearTriInterpolator(triang, z_var[i])
            zp.append(interpolator(poly_new[0], poly_new[1]))  # n.b. zp points in outblocked areas are = nan
    else:
        interpolator = tri.LinearTriInterpolator(triang, z_var)
        zp = interpolator(poly_new[0], poly_new[1])

    # nan_to_num converts list into numpy array, splits the i arrays, converts nan to 0.0 for plotting
    zp = np.nan_to_num(zp)

    return zp

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
POLYLINE: calculate distances / stationing for the x-axis in the plot
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
poly_new, distances = Get_polyline_distances(Input_polyline)

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Create interpolated polylines from bottom and free surface
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
# Which frames to extract, attention frame list starts from 0!
frames = range(0,len(times),1)
frame_bottom = [0] # in this case for the non movable bottom we need only one time step

# Which variables
# BOTTOM
z_var = Get_Values(slf, 'BOTTOM', frame_bottom)
zp_BOTTOM = Interpolate_1D_profiles(slf, slf.MESHX, slf.MESHY, slf.IKLE3, z_var, frame_bottom, poly_new)

# FREE SURFACE
z_var = Get_Values(slf, 'FREE SURFACE', frames)
zp = Interpolate_1D_profiles(slf, slf.MESHX, slf.MESHY, slf.IKLE3, z_var, frames, poly_new)

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
# DEFINE THE PLOT for the video
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

fig_width, fig_height = 10.2, 5.7375 # inch # = 16:9
title = 'flussbuero - simulation'
xlabel = 'Station [m]'
ylabel = 'Bottom and Free Surface [m]'
legend_title = 'Time [s]'

# STYLE FORMATTING, e.g. font size and line width
title_fs = 16
label_fs = 14
ax_fs = 14
legend_title = 'Time [s]'
legend_title_fs = 14
legend_fs = 14
line_width = 1.5

# LAYOUT OPTIONS
fig = plt.figure(figsize=(fig_width,fig_height),facecolor=None,linewidth=1.5,edgecolor='k')

ax1 = plt.subplot(111)
plt.title(title, fontsize = title_fs)

ax1.set_xlabel(xlabel, fontsize = label_fs)
ax1.set_ylabel(ylabel, fontsize = label_fs)

for label in ax1.xaxis.get_ticklabels():
    label.set_fontsize(ax_fs)
for label in ax1.yaxis.get_ticklabels():
    label.set_fontsize(ax_fs)

# Min and Max values of your plot
X_min = 0.0
X_max = distances[-1]
ax1.set_xticks(np.arange(X_min,X_max,5.0))
plt.xlim([X_min-0.1,X_max+0.1])

Y_min = np.amin(zp_BOTTOM)
Y_min = 0.0
Y_max = np.amax(zp)
Y_max = 1.501
ax1.set_yticks(np.arange(Y_min,Y_max,0.1))
plt.ylim([Y_min,Y_max+0.1])

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
# MAKE THE MOVIE
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
print ('Making the movie...')

# Set the path to ffmpeg.exe
plt.rcParams['animation.ffmpeg_path'] = r'..\path_to\bin\ffmpeg.exe'

# Encoding options
fps = 10 # frames per second
extra_args=['-vcodec', 'libx264', '-crf', '18'] # codec and quality

# Writer
FFMpegWriter = animation.writers['ffmpeg']
# metadata = dict(title='Nice movie', artist='cle')
writer = FFMpegWriter(fps=fps, extra_args=extra_args)

# Plot
z_FreeSurface, = plt.plot([], [], linestyle='-', color ='b', linewidth=line_width)

z_Bottom, = plt.plot([], [], linestyle='-', color ='#804000', linewidth=1.2)
z_Bottom.set_data(distances, zp_BOTTOM) # plot bottom

plt.grid(True)
fig.tight_layout()
with writer.saving(fig, "nice_video.mp4", 200): #dpi > 200: maybe performance problems when playing
    for i in range(0,len(zp)):
        z_FreeSurface.set_data(distances, zp[i])
        plt.legend(['Time = ' + str(round(times[frames[i]]/60.0,1)) + ' min'], loc = 'upper right', fontsize = legend_fs)
        writer.grab_frame()

print ('My job is done')
Posted in General and tagged , , .

5 Comments

  1. from myFunctions import Get_polyline_distances, Get_Values, Interpolate_1D_profiles
    ImportError: No module named myFunctions

    • Hi AMANj,

      for tutorial reasons I copied these three functions from the myfunctions file into the main file, see lines 41 to 81. Hence, simply you have to comment the import statement as in the above script (line 14).
      With best regards,
      Clemens

  2. Dear Flussbuero,

    Thank you for sharing this code.

    I have run this code successfully.

    Now I want to add two things.

    1) I want to show the water surface along multiple polylines in the same video (For example, along the river centreline, along the left bankline, along the right bankline. How should I achieve this? Is it possible? Which part of the code do I need to modify?

    2) If I have the river bed changing with time. (i.e. channel bottom is also movable), How should I get this changing bed elevation along with the water level change?

    Thank you.

    Regards,
    Saroj

    • Hello Saroj!

      Thank you for the interest and sorry for the late reply, busy days…

      1) You can read in, e.g. with the “np.genfromtxt” function, multiple polylines and apply exactly the same procedure as in the shown python script!
      2) You do repeat the same procedure as with the water surface, lines 101 and 102 but instead of the variable water surface you have to use the bottom. Take care that in the animation section, lines 175 – 179 you have to put the time dependent bottom variable also into the loop.

      Hope this helps!

      Clemens

  3. Dear Clemens,

    Thank you for your reply. I will try your advice and get back to you if any problem arises.

    Regards,
    Saroj

Leave a Reply

Your email address will not be published. Required fields are marked *