Raster-based example¶

This example demonstrates how to run the raster-based flood damage assessment using a land-use map for the Kampen region in the Netherlands.

Import Dependencies and setup data paths¶

In [ ]:

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import rioxarray as rxr

from damagescanner.raster import RasterScanner

from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd import rioxarray as rxr from damagescanner.raster import RasterScanner

In [33]:

# Setup data paths relative to this notebook
data_path = Path("../../data")
hazard_file = data_path / "kampen" / "hazard" / "1in1000_inundation_map.tif"
landuse_file = data_path / "landuse" / "landuse_map.tif"
curve_file = data_path / "kampen" / "vulnerability" / "curves.csv"
maxdam_file = data_path / "kampen" / "vulnerability" / "maxdam.csv"

# Setup data paths relative to this notebook data_path = Path("../../data") hazard_file = data_path / "kampen" / "hazard" / "1in1000_inundation_map.tif" landuse_file = data_path / "landuse" / "landuse_map.tif" curve_file = data_path / "kampen" / "vulnerability" / "curves.csv" maxdam_file = data_path / "kampen" / "vulnerability" / "maxdam.csv"

Load data and visualize input¶

In [34]:

curves: pd.DataFrame = pd.read_csv(curve_file)
display(curves.head())

maxdam: pd.DataFrame = pd.read_csv(maxdam_file)
display(maxdam.head())

landuse_data = rxr.open_rasterio(landuse_file).sel(band=1)
plt.figure(figsize=(10, 8))

plt.imshow(landuse_data, cmap="tab20")

plt.title("Land Use map")
plt.axis("off")
plt.show()

hazard_data = rxr.open_rasterio(hazard_file).sel(band=1)
plt.figure(figsize=(10, 8))
plt.imshow(hazard_data, cmap="Blues")
plt.colorbar(label="Water depth (m)", shrink=0.6)
plt.title("Hazard map")
plt.axis("off")
plt.show()

curves: pd.DataFrame = pd.read_csv(curve_file) display(curves.head()) maxdam: pd.DataFrame = pd.read_csv(maxdam_file) display(maxdam.head()) landuse_data = rxr.open_rasterio(landuse_file).sel(band=1) plt.figure(figsize=(10, 8)) plt.imshow(landuse_data, cmap="tab20") plt.title("Land Use map") plt.axis("off") plt.show() hazard_data = rxr.open_rasterio(hazard_file).sel(band=1) plt.figure(figsize=(10, 8)) plt.imshow(hazard_data, cmap="Blues") plt.colorbar(label="Water depth (m)", shrink=0.6) plt.title("Hazard map") plt.axis("off") plt.show()

	Unnamed: 0	110	111	112	120	130	131	132	133	134	...	520	530	540	550	560	610	620	630	640	650
0	0.0	0.0	0.00	0.00	0.0	0.0	0.00	0.00	0.00	0.00	...	0.0	0.0	0.0	0.0	0.0	0.0	0	0	0	0
1	0.1	0.1	0.05	0.05	0.1	0.1	0.05	0.05	0.05	0.05	...	0.2	0.2	0.2	0.2	0.2	0.2	0	0	0	0
2	0.2	0.2	0.09	0.09	0.2	0.2	0.09	0.09	0.09	0.09	...	0.4	0.4	0.4	0.4	0.4	0.4	0	0	0	0
3	0.3	0.3	0.12	0.12	0.3	0.3	0.12	0.12	0.12	0.12	...	0.6	0.6	0.6	0.6	0.6	0.6	0	0	0	0
4	0.4	0.4	0.15	0.15	0.4	0.4	0.15	0.15	0.15	0.15	...	0.8	0.8	0.8	0.8	0.8	0.8	0	0	0	0

5 rows × 53 columns

	object_type	damage
0	110	50.0
1	111	1600.0
2	112	1000.0
3	120	20.0
4	130	50.0

Run Assessment¶

In [35]:

damage_df, damagemap, landuse_data, hazard_data = RasterScanner(
    landuse_file,
    hazard_file,
    curves,
    maxdam,
)

print(f"Total Damage: €{damage_df['damage'].sum():,.0f}")
damage_df.head()

damage_df, damagemap, landuse_data, hazard_data = RasterScanner( landuse_file, hazard_file, curves, maxdam, ) print(f"Total Damage: €{damage_df['damage'].sum():,.0f}") damage_df.head()

Total Damage: €2,235,699,598

Out[35]:

	damage
landuse
110	131816800
111	545526000
112	32535500
120	41192520
130	5526800

Visualize Output¶

In [36]:

# Mask out zero values
damagemap_masked = np.where(damagemap == 0, np.nan, damagemap)

plt.figure(figsize=(10, 8))
plt.imshow(damagemap_masked / 1_000, cmap="YlOrRd", vmin=0)
plt.colorbar(label="Damage (thousands €)", shrink=0.6)
plt.title("Damage map")
plt.axis("off")
plt.show()

# Mask out zero values damagemap_masked = np.where(damagemap == 0, np.nan, damagemap) plt.figure(figsize=(10, 8)) plt.imshow(damagemap_masked / 1_000, cmap="YlOrRd", vmin=0) plt.colorbar(label="Damage (thousands €)", shrink=0.6) plt.title("Damage map") plt.axis("off") plt.show()