Cet article a pour but de présenter une utilisation possible de Lekan, sous sa version 2.3, avec un cas concret de construction d’une modèle hydraulique 2D d’une cours d’eau, depuis la préparation des données topographiques à la visualisation des résultats.
A la fin de cet article, il est également présenté les nouvelles possibilités qui s’ouvre avec le travail en cours (disponible prochainement en version expérimental) sur la création d’une API en language Python. Dans notre cas, la résolution du modèle et l’exploitation de ces résultats sont d’automatisés pour produire des données raster de champs d’inondation en temps réel et éventuellement envoyer une alerte en cas de dépassement d’un seuil.
La zone d’étude, le Fango fleuve côtier de Corse
Le choix de la zone d’étude a été fait de façon à pouvoir rassembler des données aujourd’hui exploitables sans nécessiter d’acquisition supplémentaire. Il est bien évident que tous secteurs ne peuvent disposer de données exploitables telles que dans ce cas d’usage. Ce choix orienté est délibéré dans un but de démonstration. L’application de ce cas d’usage sur d’autres secteurs peut ainsi demander d’acquérir des données complémentaires.
La zone choisie correspond à un tronçon de cours d’eau sur un secteur :
– Inclus dans la couverture actuelle des données LIDAR HD de l’IGN
– Avec couvert végétal relativement faible, permettant ainsi d’obtenir un modèle numérique de terrain convenable à parti des données LIDAR
– Disposant d’une station hydrométrique fournissant les débits “temps réel” à partir du service Hub’Eau
Préparation d’un MNT
L’IGN a commencé à fournir des levés LIDAR haute définition sous forme de nuage de points, permettant une restitution extrêmement précise de la topographie terrestre. Les données aujourd’hui fournies sont cependant brutes et non classifiées. On se retrouve ainsi avec un ensemble de points contenant aussi bien le terrain naturel que la végétation, les bâtiments, et toute structures existantes. Si ces données peuvent constituer un modèle numérique de surface (MNS), la donnée capitale pour tout modèle hydraulique est le modèle numérique de terrain (MNT), c’est à dire uniquement les surfaces du terrain sans la végétation et les infrastructures plus ou moins transparentes aux écoulements.
Avec un levé LIDAR classifié, c’est à dire que chaque point est associé à une nature (terrain, végétation, bâtiment, …), il est aisé d’obtenir un MNT après filtrage de ces catégories.
Malheureusement, le levé LIDAR de l’IGN n’est pour l’instant pas classifié. Il demande donc un traitement plus complexe. Pour ce traitement, le logiciel libre et open source Cloud Compare propose certains algorithmes permettant de faire le tri. Cela tombe bien, car Cloud Compare permet également de créer un MNT raster à partir de nuages de points.
Une fois le MNT raster exporté au format GeoTIFF, un petit nettoyage sous QGIS s’impose pour supprimer les triangulations sur les bords et imposer le système de coordonnées.
L’utilisation de Cloud Compare et de QGIS n’étant pas le sujet de cet article, nous ne nous attarderons pas sur ce traitement, et nous contentons de monter l’avant / l’après.
Limite de la méthode
Le traitement effectué ici est global et peut entraîner des imprécisions et erreurs dans le MNT. En effet, le tri effectué supprime un nombre important de points suivant certains critères, certains points du terrain naturel pouvant être supprimés et des “trous” pouvant apparaître dans le nuage de point. La conversion en raster interpole ensuite les valeurs d’élévation dans ces trous en fonction de ce qui l’entoure, généralement par triangulation. Ces interpolations sont ensuite bien visibles sur le MNT raster, et peuvent en fonction de leur importance, influencer la modélisation.
Construction du modèle sous Lekan
Maintenant que nous avons un MNT avec une précision plus que convenable, malgré certaines interpolations hasardeuses, nous pouvons l’importer dans Lekan.
Géométrie
La première étape pour la construction d’un modèle 2D est la saisie du domaine du modèle.
Lorsque l’utilisateur ajoute une structure 2D (nom des modèles 2D sous Lekan), la saisie du domaine se fait directement sur la carte.
Il est également possible d’utiliser un polygone d’une couche vectorielle existante, la sélection se faisant par clique droit sur le polygone lors de la saisie du domaine.
Une fois la saisie du domaine terminée, l’entrée dans la fenêtre d’édition, via la fenêtre des propriétés de ce nouvel élément du réseau hydraulique, va permettre d’éditer le modèle.
La fenêtre d’éditions de la structure 2D permet de réaliser toutes les opérations de construction du modèle à partir de plusieurs pages, chacune et dédiée à :
- Edition de la structure générale
- Réglage de la résolution du maillage
- Topographie à appliquer sur le maillage
- Edition des faces et sommets individuellement ou en groupe et vérification de la qualité du maillage
- Rugosité du modèle par zone
- Réglage de la configuration de la simulation
- Définition de la fenêtre temporelle, c’est à dire les dates et heures de début et de fin de simulation
A partir de la page “Structure “, l’utilisateur peut éditer les lignes de structure du modèle, c’est à dire le domaine mais également des lignes de structure interne, ajouter des trous dans le domaine et insérer des conditions aux limites.
C’est également ici que se fait le choix de l’algorithme pour la génération du maillage. Lekan utilise en interne le générateur Gmsh. Ce générateur propose différents algorithmes pour générer des maillages. Nous utiliserons ici algorithme BAMG qui permet d’obtenir une triangulation très régulière.
Un clique sur “Générer maillage” va débuter la génération du maillage avec la résolution par défaut.
Even though this mesh could benefit from further refinement through structure lines, it is now possible to apply a topography for an initial visualization. This is done on the “Topography” page. Here, the user adds one or more DTM to a list that will be used to apply an elevation value to each mesh vertex.
La superposition de MNT dans une liste permet de mêler plusieurs MNT : si le premier MNT n’a pas de valeur valide sous le sommet, le deuxième MNT est utilisé ; si celui-là n’a pas de valeur valide, le troisième est utilisé, et ainsi de suite …
Dans notre cas, nous utilisons le MNT issu du LIDAR HD, et un MNT beaucoup plus grossier (résolution 25m) pour les quelques sommets en périphérie qui ne sont pas couverts par le MNT fin.
Le bouton “Appliquer Topographie” démarre le processus et permet de donner du relief à notre maillage. Le bouton “Style Terrain” peut être utilisé pour personnalisé le rendu 2D de l’altimétrie du terrain.
La vue 3D est alors très utile pour contrôler la qualité du maillage.
Comme attendu, la résolution du maillage ne paraît pas suffisante en particulier dans lit mineur de notre cours d’eau. Il est donc temps d’affiner tout ça.
L’ajout de lignes de structure permet de dessiner le contour du lit mineur et de ces bifurcations. La résolution du maillage se définit par des classes de résolutions associées à des polygones dessinés sur la carte, un clique droit permet de dessiner directement un polygone défini par des lignes de structure.
La page “Edition / Qualité du maillage” identifie les éléments (faces et sommets) qui dépassent certains seuils de qualité. L’utilisateur peut alors éditer les éléments individuellement pour corriger ce que bon lui semble en fonction de ses critères. Dans notre cas nous ne changerons rien et passons à la rugosité.
Rugosité
Cette rugosité se définit de la même façon que la résolution du maillage, c’est à dire par des polygones associés à des classes de valeur.
Nous créons ici arbitrairement deux classes de valeurs à 0,035 et 0,05 (coefficient de Manning) pour le lit mineur et le lit moyen et donnant une valeur de 0,1 pour la valeur défaut qui sera appliquée partout ailleurs.
Un clique droit à l’intérieur d’un polygone formé par des lignes de structure permet d’associer directement ce polygone à une classe de rugosité.
Notez que les valeurs de rugosité sont totalement arbitraires et, dans un cas réel d’étude, devraient être fixées avec plus de sérieux. Dans notre exemple, elles devraient tout de même nous permettre d’avoir des résultats.
Conditions aux limites
La saisie de la condition aux limites se fait dans la page “Structure” de l’édition du modèle en choisissant tout simplement les deux sommets du contour du domaine délimitant la conditions aux limites.
La structure 2D étant un élément d’un réseau hydraulique qui est connecté aux autres éléments par ses conditions aux limites, le type des conditions aux limites dépend de ce qui y est connecté. Trois cas sont possibles :
- rien n’y ait connecté : le type de la conditions aux limites doit être défini dans ses propriétés par l’utilisateur.
- un lien du réseau y arrive : la condition aux limites est une entrée de débit.
- un lien du réseau la quitte : la condition aux limites est une condition de niveau d’eau.
Dans notre cas, nous allons utiliser le débit provenant du service Hub’Eau. Dans Lekan, ce débit sera produit dans un nœud que nous connecterons à la condition aux limites amont.
Pour la limite aval, ici nous fixerons arbitrairement un niveau d’eau constant raisonnable.
Station du service Hub’Eau
La zone d’étude a été en partie choisie pour être à l’aval immédiat d’une station hydrométrique en fonctionnement et présente dans la base de données du service distant Hub-Eau. Il serait donc dommage de pas l’utiliser dans notre modèle.
Pour créer directement un nœud correspondant à cette station, il faut passer par “Ajouter un nœud depuis…”, sélectionner la source “Hub’Eau”, choisir la bonne station sur la carte et l’ajouter.
Pour envoyer son débit dans la structure, il suffit alors de créer un lien depuis le nœud de la station jusqu’à la condition aval.
Astuces :
Le solveur que nous allons utiliser est TELEMAC 2D avec les équations des “volumes finis”. Avec ce schéma, lorsque l’écoulement est torrentiel, la condition d’entrée de débit dans le modèle doit être accompagnée de renseignements sur le profile de vitesse. Or, le tronçon de cours d’eau possède ici une pente de l’ordre de 1%, ce qui engendrera certainement un écoulement torrentiel. De plus, nous allons faire démarrer le modèle à “sec”, la section de la limite amont y sera alors certainement en régime super critique dès le démarrage.
Ce surcreusement peut être créé en éditant les éléments par la modifications des valeurs d’élévations des sommets. Dans ce cas, il serait nécessaire de récréer ce surcreusement à chaque fois que le maillage sera régénéré et/ou la topographie appliquée. Il est préférable de créer un petit MNT raster qui sera ajouté en haut de la liste de MNT du modèle, et donc le surcreusement sera construit à chaque application de la topographie.
Paramétrage de la simulation
Une simulation d’un modèle correspond à un solveur et ces paramètres intrinsèques. Ici, il s’agit d’une simulation TELEMAC 2D pour laquelle il faut préciser :
- le pas de temps de la simulation
- les conditions initiales
- le schéma numérique (éléments finis ou volumes finis) et d’éventuels paramètres associés
Pour la définition des conditions initiales, il y a le choix entre :
- fixer un niveau constant à vitesse nulle sur l’ensemble du domaine
- utiliser le résultat d’une autre simulation du même modèle
- définir le niveau d’eau à vitesse nulle en fonction d’une ligne d’interpolation sur l’ensemble du domaine.
C’est cette dernière solution que nous utiliserons ici. En effet, même si quasiment tout le modèle est à “sec” au démarrage, nous avons vu que nous souhaitons que la section amont soit “mouillée” au démarrage au niveau du surcreusement artificiel. de l’autre coté à l’aval, il est conseillé de partir avec niveau égal au niveau constant que nous avons fixé. Compte tenu de la pente, le niveau amont est environ 10 m au dessus du niveau aval. Nous ne pouvons donc utiliser un niveau constant sur l’ensemble du modèle.
Le principe de la ligne d’interpolation est exposé ici. En plaçant judicieusement cette ligne d’interpolation, et en y mettant la valeur de niveau aval d’un côté et un niveau permettant de “remplir sans débordement” le surcreusement, nous aurons nos conditions initiales qu’il nous faut en amont et en aval.
Fenêtre temporelle
La fenêtre temporelle correspond aux dates/heures de début et de fin. Elle peut être définie soit automatiquement en fonction des conditions aux limites, soit directement par l’utilisateur.
Dans notre cas, l’hydrogramme du serveur Hub’Eau indique l’apparition d’un pic intéressant le 15 avril 2023. Nous allons simuler ce pic en définissant la fenêtre temporelle manuelle de façon appropriée.
Arrivée à ce point, la construction du modèle est terminée. Il ne reste plus que la simulation qui ne prendra que quelques minutes.
Résultat de la simulation
Once the simulation is completed, the user can visualize the flood field with:
- vue en plan avec les valeurs de niveau d’eau, hauteur d’eau, vitesse d’écoulement ou élévation du terrain suivant des palettes de couleurs, visualisation du champ de vitesse vectoriel par des flèches, lignes de courant ou traces (statiques ou dynamiques)
- profils en long ou en travers avec niveau du terrain, niveau d’eau et vitesse d’écoulement
- vue 3D avec mappage des valeurs de niveau d’eau, hauteur d’eau ou vitesse d’écoulement
Ci-dessous un échantillon de visualisation de résultats issus de captures d’écran de Lekan (cliquez dessus pour agrandir).
Lekan et Python : contrôle et personnalisation
En cours d’expérimentation, une API pour le language Python est en cours d’élaboration. Il sera alors possible de créer, paramétrer tous éléments de Lekan avec du code Python, mais également de les contrôler et d’exploiter les résultats.
Par exemple, si nous reprenons le modèle objet de cet article, nous pouvons créer un scripts en Python permettant de lancer une simulation chaque fois qu’une nouvelle valeurs de débits est disponible sur le service Hub’Eau de la station amont. En fin de chaque simulation, les résultats sont enregistrés au format raster et peuvent être envoyés ,par exemple, vers une base de données (pouvant être affiché sur le web). La valeur du niveau d’eau est surveillée en un point spécifique et une alerte peut être envoyée quelque part en cas dépassement.
Ce script de 150 lignes présent ci-dessous constitue ainsi un système de surveillance de crue minimaliste et peut être exécuter en dehors de Lekan.
from reos.core import *
from PyQt5.QtCore import Qt, QTimer, QDateTime
import sys
import os
def to_std_out(string):
print(string, file=sys.stdout)
# A Class that will run a model when a specified hydrograph has new value
class RunController:
def __init__(self, nw, struct, hydr, run_ini, run1, run2, crs):
self.network = nw
self.structure = struct
self.hydrograph = hydr
self.run_ini = run_ini
self.next_run_scheme_id = run1
self.last_run_scheme_id = run2
self.crs = crs
self.timer = QTimer()
self.last_time = self.hydrograph.timeWindow().end()
self.runCount = 0
self.water_level_position = ReosSpatialPosition(6049297.8,7886801.5, self.crs)
self.water_level_threshold = 28.5
def start(self):
self.timer.timeout.connect(self.new_run)
self.export_to_tif(self.run_ini)
self.network.setCurrentScheme(self.next_run_scheme_id)
to_std_out('Current scheme is now {}'.format(self.network.currentSchemeName()))
self.structure.currentSimulation().setHotStartSchemeId(self.run_ini)
self.swap_scheme()
self.timer.start(60000)
def swap_scheme(self):
self.next_run_scheme_id, self.last_run_scheme_id = self.last_run_scheme_id, self.next_run_scheme_id
def export_to_tif(self, run_id):
print('Export to GeoTIFF for time: ', str(self.last_time.toString(Qt.ISODate)))
file_name_wd = '/home/where/to/store/waterdepth_{}.tif'.format(str(self.runCount))
self.structure.rasterizeResult(self.last_time, ReosHydraulicSimulationResults.DatasetType.WaterDepth, run_id,
file_name_wd, self.crs, 1)
file_name_vel = '/home/where/to/store/velocity_{}.tif'.format(str(self.runCount))
self.structure.rasterizeResult(self.last_time, ReosHydraulicSimulationResults.DatasetType.Velocity, run_id,
file_name_vel, self.crs, 1)
self.runCount = self.runCount + 1
send_raster_somewhere(file_name_wd, file_name_vel)
def check_water_level(self, run_id):
water_level=self.structure.resultsValueAt(self.last_time,
self.water_level_position,
ReosHydraulicSimulationResults.DatasetType.WaterLevel,
run_id)
if water_level>=self.water_level_threshold:
send_alert_somewhere()
def new_run(self):
self.timer.stop() # we stop the timer until we finished the new run
self.hydrograph.reloadBlocking(60000, True)
new_hydrograph_time = self.hydrograph.timeWindow().end()
to_std_out('*******************************************************************')
if new_hydrograph_time <= self.last_time:
to_std_out(
'At {} , the last hydrograph time is not updated ({})'.format(QDateTime.currentDateTime().toString(),
self.last_time.toString()))
self.timer.start(60000)
return
to_std_out('*******************************************************************')
to_std_out('New run will start')
run_duration = ReosDuration(-self.last_time.msecsTo(new_hydrograph_time), ReosDuration.millisecond)
self.last_time = new_hydrograph_time
self.structure.timeWindowSettings().setStartOffset(run_duration)
self.structure.currentSimulation().setHotStartUseLastTimeStep(True)
self.structure.runSimulation(network.calculationContext())
self.export_to_tif(self.last_run_scheme_id)
self.network.setCurrentScheme(self.next_run_scheme_id)
to_std_out('Current scheme is now {}'.format(self.network.currentSchemeName()))
self.structure.currentSimulation().setHotStartSchemeId(self.last_run_scheme_id)
self.swap_scheme()
self.timer.start(60000)
# Creation of the ReosApplication that contains the modules
app = ReosApplication(list(map(os.fsencode, sys.argv)))
reos_core = app.coreModule()
# We open the project and get the network
reos_core.openProject('/home/path/to/project/Fango-use-case/fango.lkn')
network = reos_core.hydraulicNetwork()
crs = reos_core.gisEngine().crs()
# Get the junction upstream of the structure 2D and the associated hydrograph
junctions = network.hydraulicNetworkElements(ReosHydrographJunction.staticType())
for junct in junctions:
if (junct.elementName() == 'Le Fango à Galéria'):
upstream_junction = junct
if upstream_junction is None:
sys.exit('Upstream junction not found')
hydrograph = upstream_junction.internalHydrograph()
if hydrograph is None:
sys.exit('No hydrograph associated with the upstream junction')
# make sure we have all the loaded
hydrograph.reloadBlocking(60000, True)
# get the 2D structure
structures = network.hydraulicNetworkElements(ReosHydraulicStructure2D.staticType())
if len(structures) == 0:
sys.exit('No 2D structure found')
structure2d = structures[0]
to_std_out('Structure name: {}'.format(structure2d.elementName()))
# The aim of the first run is to initialize the model so we take a time window with 120 mn before end to end
# to be sure that all the model will be well "wet".
# we use the hydraulic scheme with index 1 that is configured to start with a interpolation line (done under Lekan)
network.setCurrentScheme(1)
init_run_scheme_id = network.currentSchemeId() # id of the scheme is stored to be used with the controller
time_window_settings = structure2d.timeWindowSettings()
time_window_settings.setAutomaticallyDefined(True)
time_window_settings.setOriginStart(ReosTimeWindowSettings.End)
offset_from_end = ReosDuration(-120, ReosDuration.minute)
time_window_settings.setStartOffset(offset_from_end)
last_time = hydrograph.timeWindow().end()
structure2d.runSimulation(network.calculationContext())
# We change to schemes that have a hot start configuration and set them to use the last time step of the used scheme
network.setCurrentScheme(2)
run_id_1 = network.currentSchemeId() # id of the scheme is stored to be used with the controller
structure2d.currentSimulation().setHotStartUseLastTimeStep(True)
network.setCurrentScheme(3)
run_id_2 = network.currentSchemeId() # id of the scheme is stored to be used with the controller
structure2d.currentSimulation().setHotStartUseLastTimeStep(True)
network.setCurrentScheme(init_run_scheme_id)
controller = RunController(network, structure2d, hydrograph, init_run_scheme_id, run_id_1, run_id_2, crs)
controller.start()
app.exec()
Laisser un commentaire