Tutorial 3: Large rasters#
This tutorial demonstrates floodsr with a larger low-res flood raster and an end-to-end CLI run that fetches HRDEM and writes the super-resolved output in one shot.
What you’ll do:
verify the extended install from Extended Install
download and visualize a larger test raster
run
floodsr tohrwith--fetch-hrdeminspect both the fetched HRDEM and final output
Install floodsr#
This tutorial requires the extended install described in Extended Install. If you followed the instructions for the Basic Install, the simplest way to upgrade is to install a new environment per Extended Install. And you should probably uninstall the old/previous basic environment to avoid confusion and free up disk space.
Proceed with Extended Install based on your execution context:
command line (CLI): copy and paste the commands from Command line (CLI) - Extended into your terminal.
local notebook (Jupyter): follow the same extended install as for the CLI, then run the additional steps from Local notebook (Jupyter) - Extended so this notebook uses the correct kernel.
hosted notebook (Colab): Probably a bad idea. but see Hosted notebook (Colab) - Extended, or uncomment the experimental setup cell below if you’re feeling brave.
# 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
Install additional packages#
For this tutorial, we need the additional Python packages matplotlib and rasterio as shown in Tutorial 2
After you’ve done this, check versions:
import matplotlib
import rasterio
print(f"matplotlib=={matplotlib.__version__}")
print(f"rasterio=={rasterio.__version__}")
matplotlib==3.10.8
rasterio==1.5.0
Imports#
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
Set Cache#
When rasters become large, caching downloads to disk can be a huge time saver.
By default, this notebook follows the same cache behavior as the source code: if no cache override is supplied, floodsr uses its normal platform cache location.
If you want to override that behavior, edit the below cell and set base_cache_dir to a path you want floodsr to reuse.
# Edit this cell if you want to override the default cache behavior.
base_cache_dir = os.getenv("FLOODSR_SHARED_CACHE_DIR")
if base_cache_dir is not None:
base_cache_dir = Path(base_cache_dir).expanduser().resolve()
base_cache_dir.mkdir(parents=True, exist_ok=True)
if base_cache_dir is None:
print("cache_dir=None -> floodsr will use its default platform cache")
else:
print(f"cache_dir={base_cache_dir}")
cache_dir=/home/cefect/LS/09_REPOS/04_TOOLS/floodsr/_cache
Set Output Paths#
# 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()
result_fp = Path("result.tif").resolve()
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-QmtjEu/run/inunriver_historical_000000000WATCH_1980_rp01000_fetched_hrdem.vrt
result_fp=/home/cefect/LS/10_IO/2407_FHIMP/tmp/floodsr-tutorial_3-QmtjEu/run/result.tif
Download Test Data#
We use a clipped tile from the old Aqueduct Floods Hazard Maps dataset.
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}"
Plot the Low-Res Input#
with rasterio.open(lowres_fp) as src:
lowres = src.read(1).astype(float)
if src.nodata is not None:
lowres[lowres == src.nodata] = np.nan
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("Low-res input depth")
plt.axis("off")
plt.tight_layout()
print(f"LR depth shape: {lowres.shape}")
LR depth shape: (30, 60)
yuck
Fetch HRDEM and Run tohr#
For this larger example, we use one CLI command to fetch HRDEM, resolve the model, and run super-resolution in one shot.
The fetched HRDEM is projected (EPSG:3979) while the low-res depth raster is geographic (EPSG:4326), so we pass --crs-policy=use-dem to reproject the depth raster onto the DEM grid.
The --fetch-out argument keeps a copy of the fetched HRDEM so we can inspect it in the next section.
If you left base_cache_dir as None, this command uses the default floodsr cache location.
# Run the one-shot fetch + super-resolution command through the notebook CLI.
cache_dir_arg = "" if base_cache_dir is None else f" --cache-dir {base_cache_dir}"
!floodsr tohr --in {lowres_fp} -f --fetch-out inunriver_historical_000000000WATCH_1980_rp01000_fetched_hrdem.vrt --fetch-force-tiling{cache_dir_arg} --crs-policy use-dem --model-version ResUNet_16x_DEM --out result.tif --min-depth-threshold=0.1
INFO:floodsr.cli:starting DEM fetch
source_id=hrdem
stac_url=https://datacube.services.geo.ca/api
collection=hrdem-mosaic-1m
asset_key=dtm
gdal_available=True
force_tiling=True
fetch_window_size=8192
memory_limit_gib=16.00
stac_query_limit=200
use_project_extent_filter=True
use_cache=True
cache_dir=/home/cefect/LS/09_REPOS/04_TOOLS/floodsr/_cache
show_progress=True
project_extent_url=https://maps-cartes.services.geo.ca/server_serveur/rest/services/NRCan/coverage_HRDEM_en/MapServer/4
depth_lr_fp=
/home/cefect/LS/10_IO/2407_FHIMP/tmp/floodsr-tutorial_3-QmtjEu/run/inunriver_historical_000000000WATCH_1980_rp01000.tif
INFO:floodsr.cli:found 1 HRDEM item(s) intersecting low-res tile bounds after exact intersection filter
WARNING:floodsr.cli:forcing tiled fetch via configuration
INFO:floodsr.cli:raw fetch request grid: width=41,378, height=36,698, pixels=1,518,489,844, non_windowed_peak_estimate=19.09 GiB
INFO:floodsr.cli:HRDEM tile cache state: enabled=1, dir=/home/cefect/LS/09_REPOS/04_TOOLS/floodsr/_cache/floodsr_hrdem_tile_cache, existing_files=30, request_token=efe01d715bf4e2c5fa58b4c7
INFO:floodsr.cli:window tiling plan: source_window_shape=8,192 x 8,192, source_grid=5x6, tiles_total=30, full_tiles=20, edge_tiles=10
INFO:floodsr.cli:downloaded 6 HRDEM project extent feature(s) for clip_bounds_3979=(-1337991.5909947976, 418972.1146557964, -1296614.5561817389, 455669.36264035257), raster_shape=(36698, 41378), cache=none
HRDEM window fetch: 0%| | 0/30 [00:00<?, ?tile/s]
HRDEM window fetch: 3%|▋ | 1/30 [00:00<00:24, 1.17tile/s]
HRDEM window fetch: 7%|█▍ | 2/30 [00:01<00:20, 1.36tile/s]
HRDEM window fetch: 10%|██▏ | 3/30 [00:02<00:21, 1.28tile/s]
HRDEM window fetch: 13%|██▉ | 4/30 [00:03<00:20, 1.28tile/s]
HRDEM window fetch: 17%|███▋ | 5/30 [00:03<00:19, 1.27tile/s]
HRDEM window fetch: 23%|█████▏ | 7/30 [00:04<00:12, 1.85tile/s]
HRDEM window fetch: 27%|█████▊ | 8/30 [00:05<00:13, 1.65tile/s]
HRDEM window fetch: 30%|██████▌ | 9/30 [00:06<00:13, 1.53tile/s]
HRDEM window fetch: 33%|███████ | 10/30 [00:06<00:13, 1.45tile/s]
HRDEM window fetch: 37%|███████▋ | 11/30 [00:07<00:14, 1.35tile/s]
HRDEM window fetch: 43%|█████████ | 13/30 [00:08<00:09, 1.78tile/s]
HRDEM window fetch: 47%|█████████▊ | 14/30 [00:09<00:09, 1.63tile/s]
HRDEM window fetch: 50%|██████████▌ | 15/30 [00:09<00:09, 1.52tile/s]
HRDEM window fetch: 53%|███████████▏ | 16/30 [00:10<00:09, 1.41tile/s]
HRDEM window fetch: 57%|███████████▉ | 17/30 [00:11<00:09, 1.38tile/s]
HRDEM window fetch: 63%|█████████████▎ | 19/30 [00:12<00:06, 1.72tile/s]
HRDEM window fetch: 67%|██████████████ | 20/30 [00:13<00:06, 1.55tile/s]
HRDEM window fetch: 70%|██████████████▋ | 21/30 [00:14<00:06, 1.44tile/s]
HRDEM window fetch: 73%|███████████████▍ | 22/30 [00:14<00:05, 1.35tile/s]
HRDEM window fetch: 77%|████████████████ | 23/30 [00:15<00:05, 1.32tile/s]
HRDEM window fetch: 83%|█████████████████▌ | 25/30 [00:16<00:02, 1.96tile/s]
HRDEM window fetch: 87%|██████████████████▏ | 26/30 [00:16<00:01, 2.14tile/s]
HRDEM window fetch: 90%|██████████████████▉ | 27/30 [00:16<00:01, 2.23tile/s]
HRDEM window fetch: 93%|███████████████████▌ | 28/30 [00:17<00:00, 2.22tile/s]
HRDEM window fetch: 97%|████████████████████▎| 29/30 [00:17<00:00, 2.25tile/s]
HRDEM window fetch: 100%|█████████████████████| 30/30 [00:17<00:00, 1.69tile/s]
INFO:floodsr.cli:wrote fetched HRDEM tile to
/home/cefect/LS/10_IO/2407_FHIMP/tmp/floodsr-tutorial_3-QmtjEu/run/inunriver_historical_000000000WATCH_1980_rp01000_fetched_hrdem.vrt
INFO:floodsr.cli:loaded ORT model 'model_infer.onnx' with providers=['CPUExecutionProvider'] and scale=16
INFO:floodsr.cli:tohr path selection
requested_window_method=feather
dem_float32_bytes=6,073,959,376
windowed_io_threshold_bytes=33,554,432
selected_platform_materialization=simple
WARNING:floodsr.cli:CRS mismatch; applying --crs-policy=use-dem to resolve a canonical target grid
depth=EPSG:4326
dem=EPSG:3979
target=EPSG:3979
INFO:floodsr.cli:starting tohr inference with model_version=ResUNet_16x_DEM
model
/home/cefect/LS/09_REPOS/04_TOOLS/floodsr/_cache/ResUNet_16x_DEM/model_infer.onnx
platform_depth_lr
/home/cefect/LS/10_IO/2407_FHIMP/tmp/floodsr-tutorial_3-QmtjEu/tmp/floodsr-platform-prep-zdcm51v_/inunriver_historical_000000000WATCH_1980_rp01000_platform_depth.tif
platform_dem_hr
/home/cefect/LS/10_IO/2407_FHIMP/tmp/floodsr-tutorial_3-QmtjEu/tmp/floodsr-platform-prep-zdcm51v_/inunriver_historical_000000000WATCH_1980_rp01000_fetched_hrdem_platform_dem.tif
output
/home/cefect/LS/10_IO/2407_FHIMP/tmp/floodsr-tutorial_3-QmtjEu/run/result.tif
INFO:floodsr.cli:platform-preprocessed inputs
depth_lr shape=(30, 60) res=(689.6172468843132, 1223.2415994852063) m/pix
dem_hr shape=(36698, 41378) res=(0.9999766739102613, 0.9999795079992422) m/pix
INFO:floodsr.cli:model preprocessing complete
INFO:floodsr.cli:tohr execution path
window_method=feather
execution_path=simple
model_space_shape=(480, 960)
raw_output_shape=(36698, 41378)
INFO:floodsr.cli:prepared inputs summary:
aligned depth_lr shape=(30, 60) res=(689.6172468843132, 1223.2415994852063) m/pix
aligned dem_hr shape=(480, 960) res=(43.101077930269575, 76.45259996782539) m/pix
max_depth=5.0
dem_pct_clip=95.0
INFO:floodsr.cli:window config
method=feather
overlap_lr=8
overlap_hr=128
tile_size_lr=32
tile_size_hr=512
INFO:floodsr.cli:running feather tiling over 1x3 grid
stride_hr=384
overlap_hr=128
feather pass: 0%| | 0/3 [00:00<?, ?window/s]
feather pass: 33%|█████████ | 1/3 [00:00<00:00, 2.96window/s]
feather pass: 100%|███████████████████████████| 3/3 [00:00<00:00, 7.96window/s]
feather pass: 100%|███████████████████████████| 3/3 [00:00<00:00, 6.77window/s]
INFO:floodsr.cli:post-resampling model output from (480, 960) to (36698, 41378) on raw DEM grid with bilinear interpolation.
post-resample pass: 0%| | 0/1 [00:00<?, ?step/s]
post-resample pass: 100%|███████████████████████| 1/1 [00:18<00:00, 18.26s/step]
post-resample pass: 100%|███████████████████████| 1/1 [00:18<00:00, 18.26s/step]
final write pass: 0%| | 0/23328 [00:00<?, ?block/s]
final write pass: 2%|▎ | 456/23328 [00:00<00:05, 4552.20block/s]
final write pass: 4%|▋ | 981/23328 [00:00<00:04, 4959.40block/s]
final write pass: 6%|▉ | 1497/23328 [00:00<00:04, 5049.20block/s]
final write pass: 9%|█▎ | 2004/23328 [00:00<00:04, 5053.26block/s]
final write pass: 11%|█▌ | 2510/23328 [00:00<00:04, 5034.16block/s]
final write pass: 13%|█▉ | 3017/23328 [00:00<00:04, 5041.90block/s]
final write pass: 15%|██▎ | 3538/23328 [00:00<00:03, 5094.89block/s]
final write pass: 17%|██▌ | 4049/23328 [00:00<00:03, 5096.95block/s]
final write pass: 20%|██▉ | 4563/23328 [00:00<00:03, 5108.10block/s]
final write pass: 22%|███▎ | 5074/23328 [00:01<00:03, 5103.90block/s]
final write pass: 24%|███▌ | 5585/23328 [00:01<00:03, 5080.65block/s]
final write pass: 26%|███▉ | 6105/23328 [00:01<00:03, 5115.72block/s]
final write pass: 28%|████▎ | 6633/23328 [00:01<00:03, 5163.43block/s]
final write pass: 31%|████▌ | 7150/23328 [00:01<00:03, 5044.19block/s]
final write pass: 33%|████▉ | 7662/23328 [00:01<00:03, 5063.60block/s]
final write pass: 35%|█████▎ | 8169/23328 [00:01<00:02, 5060.32block/s]
final write pass: 37%|█████▌ | 8676/23328 [00:01<00:02, 5051.53block/s]
final write pass: 39%|█████▉ | 9182/23328 [00:01<00:02, 5051.90block/s]
final write pass: 42%|██████▏ | 9688/23328 [00:01<00:02, 4945.42block/s]
final write pass: 44%|██████ | 10184/23328 [00:02<00:02, 4898.93block/s]
final write pass: 46%|██████▍ | 10675/23328 [00:02<00:02, 4781.17block/s]
final write pass: 48%|██████▋ | 11191/23328 [00:02<00:02, 4890.68block/s]
final write pass: 50%|███████ | 11681/23328 [00:02<00:02, 4891.10block/s]
final write pass: 52%|███████▎ | 12171/23328 [00:02<00:02, 4862.51block/s]
final write pass: 54%|███████▌ | 12689/23328 [00:02<00:02, 4954.18block/s]
final write pass: 57%|███████▉ | 13193/23328 [00:02<00:02, 4978.81block/s]
final write pass: 59%|████████▏ | 13692/23328 [00:02<00:01, 4920.35block/s]
final write pass: 61%|████████▌ | 14185/23328 [00:02<00:01, 4840.14block/s]
final write pass: 63%|████████▊ | 14685/23328 [00:02<00:01, 4886.07block/s]
final write pass: 65%|█████████ | 15175/23328 [00:03<00:01, 4860.53block/s]
final write pass: 67%|█████████▍ | 15673/23328 [00:03<00:01, 4894.92block/s]
final write pass: 69%|█████████▋ | 16165/23328 [00:03<00:01, 4901.52block/s]
final write pass: 71%|██████████ | 16664/23328 [00:03<00:01, 4925.79block/s]
final write pass: 74%|██████████▎ | 17157/23328 [00:03<00:01, 4895.27block/s]
final write pass: 76%|██████████▌ | 17679/23328 [00:03<00:01, 4991.48block/s]
final write pass: 78%|██████████▉ | 18179/23328 [00:03<00:01, 4932.31block/s]
final write pass: 80%|███████████▏ | 18673/23328 [00:03<00:00, 4816.75block/s]
final write pass: 82%|███████████▌ | 19178/23328 [00:03<00:00, 4884.80block/s]
final write pass: 84%|███████████▊ | 19668/23328 [00:03<00:00, 4822.81block/s]
final write pass: 86%|████████████ | 20151/23328 [00:04<00:00, 4821.28block/s]
final write pass: 88%|████████████▍ | 20634/23328 [00:04<00:00, 3922.83block/s]
final write pass: 90%|████████████▋ | 21069/23328 [00:04<00:00, 4029.34block/s]
final write pass: 92%|████████████▉ | 21493/23328 [00:04<00:00, 4013.72block/s]
final write pass: 94%|█████████████▏| 21974/23328 [00:04<00:00, 4229.29block/s]
final write pass: 96%|█████████████▍| 22492/23328 [00:04<00:00, 4493.81block/s]
final write pass: 98%|█████████████▊| 22952/23328 [00:04<00:00, 4503.69block/s]
final write pass: 100%|██████████████| 23328/23328 [00:04<00:00, 4835.82block/s]
INFO:floodsr.cli:finished tohr inference in 45.739s; wrote 6,115,575,724 bytes to
/home/cefect/LS/10_IO/2407_FHIMP/tmp/floodsr-tutorial_3-QmtjEu/run/result.tif
/home/cefect/LS/10_IO/2407_FHIMP/tmp/floodsr-tutorial_3-QmtjEu/run/result.tif
Visualize the Fetched HRDEM#
The previous command wrote the fetched HRDEM to --fetch-out, so we can inspect the DEM tile directly.
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.bilinear).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)
Visualize the Output#
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)
