Remplissage de trous dans une matrice

Une matrice de données peut contenir des trous (valeurs manquantes, valeurs masquées, erreurs de calcul…) alors qu’on souhaite pouvoir disposer d’une matrice complète.

Ce cas de figure peut être observé dans les cas suivants (non exhaustif):

• Interpolation de données topographiques sur base d’une tringulation de points (ex: TIN)

• Données altimétriques issues de mesures LIDAR (présence de trous dans les zones d’ombre)

• Données bathymétriques issues de mesures acoustiques (présence de trous dans les zones d’ombre)

• Résultats de modélisation hydrodynamique (présence de trous dans les zones inondées - maisons, routes, …)

On souhaite alors pouvoir remplir ces trous de manière cohérente avec les données environnantes et l’objectif de calcul final.

# import _add_path # for debugging purposes only - must be removed in production

import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt

from wolfhece import is_enough
if not is_enough('2.2.43'):
    raise ImportError("This code requires wolfhece version 2.2.43 or higher. -- try pip install wolfhece --upgrade")

from wolfhece.wolf_array import WolfArray, header_wolf as hw
from wolfhece.pydownloader import toys_dataset
from wolfhece.eikonal import inpaint_array, inpaint_waterlevel

Création d’une matrice et ajout manuel de trous

Pour l’exemple, on crée une matrice de 10x10 avec des valeurs continues entre 1 et 101 puis on ajoute manuellement un trou en masquant une portion de la matrice.

La valeur nulle/NoData est imposée dans le trou. Par défaut, la valeur nulle est de 0.

h = hw()
h.set_origin(10., 20.)
h.set_resolution(0.5, 0.5)
h.shape = (10, 10)

a = WolfArray(srcheader=h)
print(a)

a.array[:,:] = np.arange(1,101).reshape((10, 10))
a.array.mask[2:4, 3:5] = True
a.set_nullvalue_in_mask()

a.plot_matplotlib()

Shape  : 10 x 10
Resolution  : 0.5 x 0.5
Spatial extent :
   - Origin : (10.0 ; 20.0)
   - End : (15.0 ; 25.0)
   - Width x Height : 5.0 x 5.0
   - Translation : (0.0 ; 0.0)
Null value : 0.0

(<Figure size 640x480 with 1 Axes>, <Axes: >)

Comptage du nombre de trous

Il est possible de compter le nombre de trous dans une matrice avec la méthode count_holes.

Si aucun argument n’est fourni, la méthode utilise le masque de la matrice.

En interne, count_holes utilise la fonction scipy.ndimage.label pour identifier les trous.

Voir https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.label.html pour plus de détails.

a.count_holes()

1

Si on fournit une matrice de masque sous forme de booléens (True et False), la méthode compte les trous (True) dans cette matrice.

La matrice de masque doit avoir la même taille que la matrice d’origine.

b = WolfArray(mold = a)
b.mask_greater(50.)
b.set_nullvalue_in_mask()
b.plot_matplotlib()

print(a.count_holes(mask = b.array.mask))

# same as above
print(a.count_holes(mask = (a.array > 50.) | a.array.mask)) # mask can be any boolean array of the same shape as a.array

2
2

Il est donc possible de compter le nombre de trous dans une matrice en utilisant un masque personnalisé.

Remplissage des trous

Le remplissage des trous se fait avec la méthode inpaint.

Elle prend en arguments optionnels:

• un masque (matrice binaire)

• une matrice de test (pour vérifier la pertinence du remplissage)

• un entier ignore_last (nombre de trous à ignorer en partant du plus grand)

• un booléen multiprocess (pour activer le calcul multiprocessus)

time, water_level, water_depth = a.inpaint(mask_array= None,    # Si on ne fournit pas de masque, on utilise le masque de a
                                           test_array= None,    # Si on ne fournit pas de matrice de test, on utilise une matrice temporaire avec comme valeur unique la valeur minimale de "a.array.data" --> sans tenir compte du masque !
                                           ignore_last= 0,      # 0 indique qu'il ne faut pasd ignorer de trous
                                           multiprocess= False) # False indique qu'on n'utilise pas le calcul multiprocessus -- clairement inutile dans cet exemple

assert isinstance(time, np.ndarray)         # time is a numpy array
assert isinstance(water_level, WolfArray)   # water_level is a WolfArray
assert isinstance(water_depth, WolfArray)   # water_depth is a WolfArray

100%|██████████| 1/1 [00:01<00:00,  1.14s/it]

Signification des arguments optionnels

• La matrice de masque a le même effet que dans le cas du compte. Elle peut se substituer au masque de la matrice d’origine.

• La matrice de test est un peu plus subtile. En effet, le remplissage ne sera considéré comme valable que si celui-ci est supérieur à la valeur de test.

• L’argument ignore_last permet d’ignorer un certain nombre de trous en partant du plus grand. Ceci est pariculièrement utile dans le cas du traitement de résultats de modélisation hydrodynamique où le plus grand trou est souvent la zone qui est autour du domaine de calcul (zone non calculée).

• Le calcul multiprocessus peut être utile pour les grandes matrices avec de très nombreux trous (> 1000). En effet, chaque trou est traité indépendamment des autres et le calcul peut donc être parallélisé.

Optimisation du calcul

En interne, la méthode de remplissage est compilée avec Numba pour accélérer le calcul.

Méthode de calcul

Le remplissage est effectué par propagation de l’information pondérée au cours de la résolution de l’équation eikonale (https://en.wikipedia.org/wiki/Eikonal_equation) par méthode Fast Marching (https://en.wikipedia.org/wiki/Fast_marching_method).

\[|\nabla T(x)| = \frac{1}{F(x)}\]

Retour de la méthode

La méthode retourne trois objets:

• une matrice numpy avec les “temps” de propagation (remplissage)

• une matrice WolfArray avec les valeurs interpolées (sur base de la matrice d’origine)

• une matrice WolfArray avec le différentiel entre la matrice interpolée et la matrice de test

time_wa = WolfArray(mold=a)
time_wa.array.data[:,:] = time[:,:]
time_wa.plot_matplotlib(with_legend = True)

(<Figure size 640x480 with 2 Axes>, <Axes: >)

water_level.plot_matplotlib()
water_depth.plot_matplotlib() # Dans ce cas, identique à water_level puisque la valeur minimale de a est 0 et que la matrice de test n'a pas été fournie (soit équivalent à une matrice de test remplie de zéros).

(<Figure size 640x480 with 1 Axes>, <Axes: >)

Application à l’interpolation de résultats hydrodynamiques

Lorsqu’un résultat de modélisation hydrodynamique, via WOLF (GPU et /ou CPU) est généré, étant donné que les bâtiments sont intégrés dans la donnée topographique, il est normal que des trous apparaissent dans les zones inondées (maisons, routes, …).

Ces éléments ne sont en effet que rarement submergés pour des crues de faible à moyenne ampleur.

Si on souhaite utiliser ces résultats pour des analyses ultérieures (ex: analyse de risques, cartographie des zones inondables, …), il est souvent nécessaire de disposer d’une matrice complète.

Une technique souvent utilisée est le remplissage par interpolation/extrapolation sur base de contours vectoriels (PICC, Cadastre…) auxquels on applique un buffer. Une autre technique est le remplissage par interpolation/extrapolation de la matrice de hauteur d’eau dans tous les trous.

Cependant, ces approches cartographiques ne sont pas toujours pertinentes au regard du modèle hydrodynamique.

En effet, si une zone n’est pas inondée, cela peut provenir de plusieurs raisons:

• la zone est en dehors du domaine de calcul

• la zone est protégée par une digue ou un ouvrage hydraulique

• la zone est constituée de bâtiments

• la zone est constituée d’éléments surélevés (ex: autoroutes, voies ferrées, colline…)

Appliquer l’une ou l’autre technique de remplissage amène souvent à des incohérences avec le modèle hydrodynamique.

La méthode d’inpainting proposée ici permet de remplir les trous en respectant les données environnantes et la logique du modèle hydrodynamique.

Pour ce faire, on va utiliser :

• la matrice d’altitudes de surface libre comme matrice de données

• la matrice de MNT (sans les bâtiments) comme matrice de test

• le masque de la matrice de base

Dès lors, les résultats fournis par la méthode inpaint seront :

• une matrice d’altitudes de surface libre complétée

• une matrice de hauteurs d’eau complétée (différence entre la matrice d’altitudes de surface libre et la matrice de MNT)

Sur cette base, on pourra utiliser les outils d’analyse spatiale de WOLF ou d’autres logiciels SIG pour poursuivre les analyses.

Remarques :

• Il faut fournir un MNT sans bâtiment MAIS avec les ouvrages de protection. Pour la Vesdre, cette information est habituellement disponible en complément du MNS ou du modèle de surface de modélisation hydrodynamique (MNT + bâtiments). Cependant, dans les bâtiments, ce MNT a été souvent construit par interpolation du voisinage. On reporte ainsi une certaine responsabilité à la qualité de ce MNT et aux technique employées pour le générer.

• Les zones émergées pour des raisons naturelles (ex: colline, talus, …) ne seront pas remplies puisque la surface libre est inférieure au MNT.

• Proche du cours d’eau, la hauteur d’eau dans le bâtiment n’est pas influencée par la hauteur d’eau importante dans le lit mineur puisqu’elle est contruite sur base de l’altitude de surface libre à laquelle est retranchée la topographie.

• Si le bâtiment est situé près du lit mineur et qu’un de ses murs empêche le débordmeent, il pourra être considéré comme inondé. Ceci n’est pas le cas si un mur anti-crue est présent entre le bâtiment et la rivière et que le mur anti-crue est bien présent dans le MNT de test.

• Si ce post-traitement de résultat est pensé dès le début de la mise en place de la modélisation hydrodynamique, la donnée de MNT peut être améliorée, notamment en incorporant des modifications/aménagements topographiques ou encore des levés de terrain.

• La connaissance précise des niveaux utiles des bâtiments est une donnée manquante connue et déjà mise en évidence à de très nombreuses reprise (TFE, thèses de doctorat, articles scientifiques, …). Des recherches sont en cours pour tenter d’améliorer cette connaissance ce qui est essentiel pour l’application de modèles de dommage où la hauteur d’eau est souvent la variable dominante.

• La méthode d’inpainting sera toujours d’actualité lorsque la donnée du MNT sera améliorée.

Méthode inpaint_waterlevel

Cette méthode permet un calcul plus direct au sens de résultats hydrodynamiques. Ce processus est utilisée dans le calcul de l’acceptabilité des risques d’inondation, employé par le SPW.

Les arguments sont :

• un matrice d’altitudes de surface libre (matrice numpy ou nump.maskedarray)

• une matrice DEM (matrice numpy)

• une matrice DTM (matrice numpy)

• un entier ignore_last_patches (nombre de trous à ignorer en partant du plus grand - valeur par défaut 1)

• un booléen inplace (pour modifier la matrice d’altitudes de surface libre en place - valeur par défaut True)

• des flottants dx et dy (résolution spatiale - valeur par défaut 1)

• un flottant NoDataValue (valeur nulle/NoData - valeur par défaut 0)

• un booléen multiprocess (pour activer le calcul multiprocessus - valeur par défaut True)

• un flottant epsilon (valeur minimale pour considérer qu’une hauteur d’eau est présente - valeur par défaut 1e-3)

Remarques :

• Les matrices doivent avoir les mêmes dimensions et une erreur sera levée si ce n’est pas le cas.

• Etant donné que les martrices sont des matrices numpy, elles ne sont pas géoréférencées. Il est du ressort de l’utilisateur de s’assurer qu’elles sont correctement alignées (même origine, même résolution).

• On peut bien entendu utiliser l’attribut .array d’un WolfArray pour fournir les matrices. Ce sont des matrices Numpy.maskedarray.

• Si la matrice d’altitudes de surface libre est un Numpy.maskedarray, les matrices retournées seront aussi des Numpy.maskedarray.

• Il est possible d’utiliser l’argument optionnel NoDataValue pour définir une valeur à ignorer dans la matrice de base (par défaut égale à 0).

• Dans la routine, les “bâtiments”/”zones à remplir” sont identifiés par différentiel entre la matrice DEM et la matrice DTM.

• Si le DTM est au-dessus du DEM (ce qui ne devrait pas se produire), on ne tiendra pas compte de ces mailles lors du remplissage.

• Un “epsilon” est un argument optionnel de la routine. Toute différence entre le DEM et le DTM inférieure à epsilon sera considérée comme nulle et ne sera donc pas remplie.

Pour info, le code de préparation avant le remplissage (cf eikonal.py) est le suivant :

if inplace:
        if isinstance(water_level, np.ma.MaskedArray):
            base_data = water_level.data
        else:
            base_data = water_level
    else:
        if isinstance(water_level, np.ma.MaskedArray):
            base_data = water_level.data.copy()
        else:
            base_data = water_level.copy()

    assert dem.shape == base_data.shape
    assert dtm.shape == base_data.shape

    # Create the mask where we can fill the water level

    # first we identify the buildings by the difference between the DEM and the DTM
    buildings = dem - dtm
    # If DTM is above DEM, we set the value to 0
    if np.any(buildings < 0.) > 0:
        logging.warning("Some cells in the DTM are above the DEM.")
        logging.info("Setting these values to 0.")
    buildings[buildings < 0.] = 0.

    # If DTM is below DEM, we set the value to 1
    buildings[buildings <= epsilon] = 0.
    buildings[buildings > 0.] = 1.

    # We interpolate only if building cells are not already in the water_level
    comp = np.logical_and(buildings == 1., base_data != NoDataValue)
    if np.any(comp):
        logging.warning("Some building cells are already flooded.")
        logging.info("Ignoring these cells in the interpolation.")

    buildings = np.logical_and(buildings, base_data == NoDataValue)

Documentation de la routine

print(inpaint_waterlevel.__doc__)

 Inpaint the water level using the Fast Marching Method (FMM). Similar to the HOLES.EXE Fortran program.

    Assumptions:
        - The simulations are in a steady state.
        - The flood extent is limited by:
            - natural topography (where DEM == DTM)
            - buildings or blocks of buildings
            - protective structures

    The goal is to propagate the free surface elevations into the buildings.

    We calculate the difference between the DEM (including buildings and walls) and the DTM to identify where the buildings are.

    Specifically:
        - if it is natural topography, the differences will be zero or almost zero
        - if it is a building or anthropogenic structure, the differences are significant

    We set the elevations of the cells where the difference is zero to a value unreachable by the flood.
    Thus, only buildings in contact with the flood will be affected and filled.

    HOLES.EXE vs Python code:
        - In "holes.exe", we must provide "in", "mask", and "dem" files:
            - "in" is the water level
            - "dem" is the digital terrain model (DTM) associated with the simulation's topo-bathymetry (not the topo-bathymetry itself)
            - "mask" is the DTM where the buildings are identified and a large value (above the maximum water level) if the cell is not inside a building.
            - The patches are identified by the water level "0.0" in the "in" file.
            - FMM is used and new water levels are calculated and retained if greater than the "mask" value.
            - The DTM is only used in the final part when evaluating the new water depths (Z-DTM).

        - In Python, we must provide the water level, the DEM, and the DTM:
            - The patches will be identified by the buildings array (filtered DEM - DTM).
            - FMM is used and new water levels are calculated and retained if greater than the DTM value.
            - FMM is only propagated in the cells where "buildings = True".

        - Final results must be the same even if the algorithm is a little bit different.
        - Fortran code demands the "mask" file to be provided (pre-computed/modified by the user), but in Python, we calculate it automatically from the DEM and DTM.
        - We can use "inpaint_array" to be more flexible... by manually providing a "where_compute" and a "test" arrays

    :param water_level: (2D numpy array): The water level to inpaint.
    :param dem: (2D numpy array): The simulation's Digital Elevation Model (DEM).
    :param dtm: (2D numpy array, optional): The digital terrain model (DTM).
    :param ignore_last_patches: (int, optional): The number of last patches to ignore.
    :param inplace: (bool, optional): Whether to update the water_level in place.
    :param dx: (float, optional): The mesh size in x direction.
    :param dy: (float, optional): The mesh size in y direction.
    :param NoDataValue: (float, optional): The NoDataValue, used if mask is not explicitly provided (mask attribute or water_level as a Numpy MaskedArray). Default is 0.
    :param multiprocess: (bool, optional): Whether to use multiprocessing.

Exemple d’utilisation

Cet exemple est une copie d’un test unitaire permettant de vérifier le bon fonctionnement de la méthode.

# Water surface elevation matrix (water level)
water_level = np.ma.zeros((100, 100))

# Digital terrain model matrix (dtm)
dtm = np.zeros((100, 100))
# Digital elevation model matrix (dem)
dem = np.zeros((100, 100))
# Matrix indicating zones to compute (1 = to compute, 0 = do not compute)
to_compute = np.zeros((100, 100))

# Define some patches with water level = 0
patches = [slice(50,52), slice(40,42), slice(60,62)]

# Leave an external border of 2 pixels with water level 0.
# It will be the largest patch --> ignore it.

# Initialize matrices
water_level[2:-2, 2:-2] = 2. # water level = 2 everywhere except on the border
dtm[:,:] = 1. # flat terrain at elevation 1
dem[:,:] = 3. # flat surface at elevation 3

# Add patches with water level = 0
for patch in patches:
    to_compute[patch, patch] = 1.
    water_level[patch, patch] = 0.

# The DTM in the first patch is set to 3, so that the water level inpainting will not be able to fill it.
dtm[patches[0], patches[0]] = 3.

# Create a masked array for water level
np.ma.masked_equal(water_level, 0., copy=False)

assert isinstance(water_level, np.ma.MaskedArray)

# Inpaint water level
time, wl, wh = inpaint_waterlevel(water_level, dem, dtm, ignore_last_patches= 1, inplace= True, dx=1., dy=1.)

assert isinstance(time, np.ndarray)
assert isinstance(wl, np.ma.MaskedArray)
assert isinstance(wh, np.ma.MaskedArray)

# Compare data
assert np.all(water_level.data[patches[0], patches[0]] == 0.)
assert np.all(water_level.data[patches[1], patches[1]] == 2.)
assert np.all(water_level.data[patches[2], patches[2]] == 2.)

# Compare masked array
water_level.mask[:,:] = water_level.data == 0.
assert np.all(water_level.mask[patches[0], patches[0]] == True)
assert np.all(water_level[patches[1], patches[1]] == 2.)
assert np.all(water_level[patches[2], patches[2]] == 2.)

WARNING:root:Some building cells are already flooded.

Multiprocessing : (OUI/NON), à quoi être attentif ?

L’activation du multiprocessing peut être utile pour les grandes matrices avec de très nombreux trous (> 1000). En deça, le gain de temps n’est pas garanti car le surcoût lié à la création des processus peut annuler les bénéfices.

Cependant, un script Python utilisant le multiprocessing doit être protégé par if __name__ == "__main__":.

Pour plus d’informations, voir cet article sur Real Python.

Si ce n’est pas le cas, sous Windows, une erreur sera très probablement levée en provenance de la tentative de création de processus enfants du paquet multiprocessing de Python. Cette erreur n’est pas liée à WOLF mais au fonctionnement même de Python.

Il faut donc prendre soin de toujours protéger le code principal par if __name__ == "__main__": lorsque le multiprocessing est activé.

if __name__ == "__main__":
    print("Démarrage du script principal")

    # Code principal
    #...

    print("Fin du script principal")