Tutoriel 3 : rasters de grande taille

Ce tutoriel présente floodsr avec un raster plus grand et une exécution CLI de bout en bout qui récupère le HRDEM et écrit la sortie super-résolue en une seule commande.

Ce que vous allez faire :

• vérifier l’installation étendue décrite dans Installation étendue

• télécharger puis visualiser un raster de test plus grand

• exécuter floodsr tohr avec --fetch-hrdem

• inspecter à la fois le HRDEM récupéré et la sortie finale

Installer floodsr

Ce tutoriel nécessite l’installation étendue décrite dans Installation étendue. floodsr utilise cette approche d’installation à deux niveaux afin de convenir au plus grand nombre d’utilisateurs possible :

• l’installation de base pour celles et ceux qui veulent quelque chose de simple à installer et n’ont pas besoin de l’outillage supplémentaire pour gérer de grands rasters, et

• l’installation étendue pour les personnes plus à l’aise avec Python et qui veulent exploiter toute la puissance de floodsr sur de grands rasters.

Si vous avez suivi les instructions de Installation de base, le moyen le plus simple de passer à niveau consiste à créer un nouvel environnement selon Installation étendue.

Il est aussi conseillé de désinstaller l’ancien environnement de base afin d’éviter toute confusion et de libérer de l’espace disque. Consultez Installation pour les commandes d’installation et de désinstallation correspondant à chaque mode et contexte d’exécution.

Poursuivez avec Installation étendue selon votre contexte d’exécution :

• ligne de commande (CLI) : copiez-collez dans votre terminal les commandes de Ligne de commande (CLI) - étendue.

• carnet local (Jupyter) : suivez les mêmes étapes Ligne de commande (CLI) - étendue que pour la CLI, puis exécutez les étapes supplémentaires de Carnet local (Jupyter) - étendu pour que ce carnet utilise le bon noyau.

• carnet hébergé (Colab) : ce n’est probablement pas une très bonne idée, mais consultez Carnet hébergé (Colab) - étendu, ou décommentez la cellule de configuration expérimentale ci-dessous si vous êtes d’humeur aventureuse.

# Colab experimental setup
# !apt-get update -qq
# !apt-get install -y -qq gdal-bin libgdal-dev
# !pip install -q --upgrade pip
# !pip install -q "gdal[numpy]==$(gdal-config --version).*" rasterio geopandas pyproj shapely fiona
# !pip install -q floodsr

Vérifier l’installation

Vérifions maintenant les versions :

import matplotlib
import rasterio
import sys

print(f"matplotlib=={matplotlib.__version__}")
print(f"rasterio=={rasterio.__version__}")
print(f"python={sys.executable}")

matplotlib==3.10.8
rasterio==1.5.0
python=/opt/conda/envs/dev/bin/python

!floodsr --version

Si cela ne fonctionne pas, vérifiez bien que vous avez sélectionné le bon noyau pour ce carnet (par exemple Python (floodsr-gdal)).

Importations

import os
from pathlib import Path
from urllib.request import urlretrieve

import matplotlib.pyplot as plt
import numpy as np
import rasterio
from rasterio.enums import Resampling

Définir le cache

Quand les rasters (cartes en grille) deviennent volumineux, mettre les téléchargements en cache sur disque peut faire gagner beaucoup de temps.

Ce tutoriel montre comment utiliser un répertoire de cache réutilisable afin que les récupérations répétées de HRDEM et l’exécution du modèle n’aient pas à repartir de zéro à chaque fois. L’espace requis par le cache, et d’ailleurs par votre répertoire de travail, dépendra bien sûr de la taille des rasters que vous manipulez. Pour de grands rasters, cela peut dépasser 10 Go.

Modifiez la cellule ci-dessous si vous voulez que floodsr réutilise un autre répertoire de cache.

# Edit this cell if you want `floodsr` to reuse a different cache directory.
base_cache_dir = "./_cache"

# windows example. un-comment and edit as needed.
# base_cache_dir = r"C:\Users\<you>\floodsr_cache"

# Resolve the cache path (and create if needed) once so later CLI cells can pass it directly.
base_cache_dir = Path(base_cache_dir).expanduser().resolve()
base_cache_dir.mkdir(parents=True, exist_ok=True)
print(f"cache_dir={base_cache_dir}")

cache_dir=/workspace/_cache

Définir les chemins

Commençons par affecter quelques chemins à des variables Python pour plus de commodité.

# Define the low-res input name and the derived output paths.
lowres_name = "inunriver_historical_000000000WATCH_1980_rp01000.tif"
hrdem_fp = Path(f"{Path(lowres_name).stem}_fetched_hrdem.vrt").resolve()
requested_result_fp = Path("result.tif").resolve()
result_fp = requested_result_fp

print(f"hrdem_fp={hrdem_fp}")
print(f"result_fp={result_fp}")

hrdem_fp=/home/cefect/LS/10_IO/2407_FHIMP/tmp/floodsr-tutorial_3-AurJ8Z/run/inunriver_historical_000000000WATCH_1980_rp01000_fetched_hrdem.vrt
result_fp=/home/cefect/LS/10_IO/2407_FHIMP/tmp/floodsr-tutorial_3-AurJ8Z/run/result.tif

Télécharger les données de test

Pour ce tutoriel, nous utilisons une tuile issue de l’ancien jeu de données Aqueduct Floods Hazard Maps.

Un extrait en a été téléversé dans le dépôt du projet pour en faciliter la récupération. Récupérez-le dans votre répertoire de travail en utilisant à nouveau la fonction utilitaire urlretrieve.

lowres_fp = Path(lowres_name).resolve()

urlretrieve(
    "https://github.com/cefect/floodsr/releases/download/v0.0.9/inunriver_historical_000000000WATCH_1980_rp01000.tif",
    lowres_fp,
)
assert lowres_fp.is_file(), f"could not find {lowres_fp}"

Cela récupère un raster comme celui-ci :

On peut aussi le tracer pour se faire une idée des données d’entrée.

with rasterio.open(lowres_fp) as src: # no helper function here because we need fancier loading below
    lowres = src.read(1).astype(float)
    if src.nodata is not None:
        lowres[lowres == src.nodata] = np.nan

    # print out some info about the input raster
    print(f"lowres raster: size={Path(lowres_fp).stat().st_size / 1024**2:.1f} MB, shape={src.shape}, res={src.res}, crs={src.crs}")


plt.figure(figsize=(6, 4))
im = plt.imshow(np.where(lowres > 0, lowres, np.nan), cmap="Blues")
plt.colorbar(im, fraction=0.046, pad=0.04)
_ = plt.title(f"Low-res input depth")
 

lowres raster: size=0.0 MB, shape=(30, 60), res=(0.008333333333333333, 0.008333333333333333), crs=EPSG:4326

beurk

Vous remarquerez que ce raster utilise EPSG:4326, c’est-à-dire des coordonnées latitude/longitude. D’où cette résolution un peu étrange. Il faudra en tenir compte au moment de définir notre --crs-policy.

Récupérer le HRDEM et exécuter tohr

Pour cet exemple plus grand, nous utilisons une seule commande CLI pour récupérer le HRDEM, exécuter le modèle et lancer la super-résolution d’un seul coup, avec les arguments suivants :

• --in {lowres_fp} indique à floodsr le raster de profondeur d’inondation basse résolution que nous avons téléchargé ci-dessus.

• -f / --fetch-hrdem indique à floodsr de récupérer directement le HRDEM au lieu d’exiger un fichier DEM local au préalable.

• --fetch-out {hrdem_fp.name} conserve sur disque une copie du HRDEM récupéré afin que nous puissions l’inspecter dans la section suivante.

• --fetch-force-tiling force l’étape de récupération du HRDEM à s’exécuter par fenêtres tuilées. C’est utile sur de grandes emprises, car cela évite d’essayer d’assembler tout le DEM en une seule grosse lecture.

• --cache-dir {base_cache_dir} indique à la fois à l’étape de récupération et à l’étape d’exécution du modèle de réutiliser le répertoire de cache configuré ci-dessus.

• --window-method hard indique au runtime ToHR d’utiliser des fenêtres rigides non superposées plutôt qu’un lissage par fusion progressive. Pour cet exemple de grand raster, nous le définissons explicitement parce que le chemin mémoire-sûr pour grands rasters est actuellement lié au fenêtrage hard.

• --crs-policy use-dem indique à floodsr de traiter la grille DEM récupérée comme grille cible lorsque le raster d’inondation d’entrée et le DEM utilisent des systèmes de coordonnées différents. Dans cet exemple, le raster d’inondation basse résolution est géographique (EPSG:4326) tandis que le HRDEM récupéré est projeté (EPSG:3979).

• --out result.tif demande un nom de base de sortie prévisible dans le répertoire de travail. Selon le moteur actif, l’artefact final peut être produit sous forme de GeoTIFF ou de VRT, et floodsr affiche le chemin de sortie résolu sur la dernière ligne de sortie.

• --min-depth-threshold=0.1 masque les très faibles profondeurs prédites dans la sortie finale.

Consultez Référence CLI pour la référence CLI complète.

# Run the one-shot fetch + super-resolution command through the notebook CLI.
!floodsr tohr --in {lowres_fp} -f --fetch-out {hrdem_fp.name} --fetch-force-tiling --cache-dir {base_cache_dir} --window-method hard --crs-policy use-dem --out {requested_result_fp.name} --min-depth-threshold=0.1

Tracer les résultats

La commande précédente a écrit le HRDEM récupéré dans --fetch-out, ce qui nous permet d’inspecter directement la tuile DEM.

Pour accélérer le tracé, nous utilisons une version rééchantillonnée du résultat. Le raster est si grand que le rendu visuel reste le même avec ce rééchantillonnage ; prenez toutefois garde à ne pas confondre le travail de visualisation avec votre travail d’analyse.

with rasterio.open(hrdem_fp) as src:
    preview_scale = max(src.height / 1024, src.width / 1024, 1)
    preview_shape = (
        max(1, int(src.height / preview_scale)),
        max(1, int(src.width / preview_scale)),
    )
    dem = src.read(1, out_shape=preview_shape, resampling=Resampling.nearest).astype(float)
    if src.nodata is not None:
        dem[dem == src.nodata] = np.nan

plt.figure(figsize=(6, 4))
im = plt.imshow(dem, cmap="terrain")
plt.colorbar(im, fraction=0.046, pad=0.04)
plt.title("Fetched HRDEM")
plt.axis("off")
plt.tight_layout()

print(f"DEM full shape: {src.shape}")
print(f"DEM preview shape: {dem.shape}")

DEM full shape: (36698, 41378)
DEM preview shape: (908, 1024)

Et maintenant, traçons le résultat haute résolution

with rasterio.open(result_fp) as src:
    preview_scale = max(src.height / 1024, src.width / 1024, 1)
    preview_shape = (
        max(1, int(src.height / preview_scale)),
        max(1, int(src.width / preview_scale)),
    )
    result = src.read(1, out_shape=preview_shape, resampling=Resampling.bilinear).astype(float)
    if src.nodata is not None:
        result[result == src.nodata] = np.nan

plt.figure(figsize=(6, 4))
im = plt.imshow(np.where(result > 0, result, np.nan), cmap="Blues")
plt.colorbar(im, fraction=0.046, pad=0.04)
plt.title("Output depth")
plt.axis("off")
plt.tight_layout()

print(f"Output full shape: {src.shape}")
print(f"Output preview shape: {result.shape}")

Output full shape: (36698, 41378)
Output preview shape: (908, 1024)