#import required libraries - eventually install if not available 
import os
import ilwis
import openeo
import numpy as np
import matplotlib.pyplot as plt
from collections import deque
from scipy.ndimage import binary_dilation
from scipy.ndimage import distance_transform_edt

#suppress warnings
import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter('ignore')

#within the notebook folder a directory is created to store your processing results
work_dir = (os.getcwd())+ r'/flood_modelling'

print("current dir is: %s" % (os.getcwd()))
print("current working directory is:",work_dir) 

if os.path.isdir(work_dir):
    print("Folder exists")
else:
    print("Folder doesn't exists")
    os.mkdir(work_dir)

current dir is: d:\jupyter\notebook_scripts\Special\DEM_FloodMask
current working directory is: d:\jupyter\notebook_scripts\Special\DEM_FloodMask/flood_modelling
Folder exists

#set the working directory for ILWISPy
ilwis.setWorkingCatalog(work_dir)
print(work_dir)

d:\jupyter\notebook_scripts\Special\DEM_FloodMask/flood_modelling

ilwis.version()

'1.0 build 20260708'

#remove previous login credentials
from openeo.rest.auth.config import RefreshTokenStore
RefreshTokenStore().remove()

#ensure openeo is installed and you have registered
connection = openeo.connect("openeo.dataspace.copernicus.eu").authenticate_oidc()

Authenticated using device code flow.

#connection = openeo.connect("openeofed.dataspace.copernicus.eu").authenticate_oidc()
#connection = openeo.connect("openeo.vito.be").authenticate_oidc()
#connection = openeo.connect("openeo.cloud").authenticate_oidc()

print(connection.list_collection_ids())

['SENTINEL3_OLCI_L1B', 'SENTINEL3_SLSTR', 'SENTINEL_5P_L2', 'SENTINEL1_GRD', 'COPERNICUS_30', 'LANDSAT_BIMONTHLY_MOSAIC', 'SENTINEL3_SYN_L2_SYN', 'SENTINEL3_SLSTR_L2_LST', 'SENTINEL5P_L2_CO', 'COPERNICUS_VEGETATION_PHENOLOGY_PRODUCTIVITY_10M_SEASON1', 'COPERNICUS_VEGETATION_PHENOLOGY_PRODUCTIVITY_10M_SEASON2', 'COPERNICUS_PLANT_PHENOLOGY_INDEX', 'ESA_WORLDCOVER_10M_2020_V1', 'ESA_WORLDCOVER_10M_2021_V2', 'SENTINEL1_GLOBAL_MOSAICS', 'SENTINEL2_L2A', 'SENTINEL2_L1C', 'SENTINEL3_OLCI_L2_LAND', 'SENTINEL3_OLCI_L2_WATER', 'SENTINEL3_SYN_L2_AOD', 'CLMS_LCM_GLOBAL_10M_YEARLY_V1', 'CLMS_TCD_PANTROPICAL_10M_YEARLY_V1', 'CLMS_BA_GLOBAL_300M_DAILY_V4', 'CLMS_NDVI_GLOBAL_300M_10DAILY_V3', 'CLMS_ETA_GLOBAL_300M_10DAILY_V1', 'CLMS_HF_GLOBAL_300M_DAILY_V1', 'CLMS_LSP_GLOBAL_300M_YEARLY_V2', 'CLMS_SSM_EUROPE_1KM_DAILY_V1', 'CLMS_NPP_GLOBAL_300M_10DAILY_V2', 'CLMS_VLCC_CROP_TYPES_EUROPE_10M_YEARLY_V1', 'CLMS_LST_TCI_GLOBAL_3KM_10DAILY_V3', 'CLMS_FAPAR_GLOBAL_1KM_10DAILY_V2', 'CLMS_FCOVER_GLOBAL_1KM_10DAILY_V2', 'CLMS_GDMP_GLOBAL_1KM_10DAILY_V2', 'CLMS_LAI_GLOBAL_1KM_10DAILY_V2', 'CLMS_NDVI_LTS_GLOBAL_1KM_10DAILY_V3', 'CLMS_NDVI_GLOBAL_1KM_10DAILY_V3', 'CLMS_SWI_EUROPE_1KM_DAILY_V1', 'CLMS_SWI_EUROPE_1KM_DAILY_V2', 'CLMS_SWI_GLOBAL_12_5KM_DAILY_V3', 'CLMS_SWI_GLOBAL_12_5KM_DAILY_V4', 'CLMS_SWI_GLOBAL_12_5KM_10DAILY_V3', 'CLMS_SWI_GLOBAL_12_5KM_10DAILY_V4', 'CLMS_SWI_STATIC_GLOBAL_12_5KM_DAILY_V3', 'CLMS_SCE_EUROPE_500M_DAILY_V1', 'CLMS_SCE_NORTHERNHEMISPHERE_1KM_DAILY_V1', 'CLMS_SCE_GLOBAL_1KM_DAILY_V1', 'CLMS_SWE_NORTHERNHEMISPHERE_5KM_DAILY_V1', 'CLMS_SWE_NORTHERNHEMISPHERE_5KM_DAILY_V2', 'CLMS_LST_GLOBAL_5KM_HOURLY_V1', 'CLMS_LST_GLOBAL_5KM_HOURLY_V2', 'CLMS_LST_DAILY_CYCLE_GLOBAL_5KM_10DAILY_V1', 'CLMS_LST_DAILY_CYCLE_GLOBAL_5KM_10DAILY_V2', 'CLMS_LST_TCI_GLOBAL_5KM_10DAILY_V1', 'CLMS_LST_TCI_GLOBAL_5KM_10DAILY_V2', 'CLMS_LSWT_NRT_GLOBAL_1KM_10DAILY_V1', 'CLMS_LSWT_OFFLINE_GLOBAL_1KM_10DAILY_V1', 'CLMS_BA_GLOBAL_300M_MONTHLY_V4', 'CLMS_FCOVER_GLOBAL_300M_10DAILY_V2', 'CLMS_FAPAR_GLOBAL_300M_10DAILY_V2', 'CLMS_LAI_GLOBAL_300M_10DAILY_V2', 'CLMS_NDVI_STS_GLOBAL_1KM_10DAILY_V3_NC', 'CLMS_DMP_GLOBAL_300M_10DAILY_V2', 'CLMS_GDMP_GLOBAL_300M_10DAILY_V2', 'CLMS_DMP_GLOBAL_1KM_10DAILY_V2', 'CLMS_GPP_GLOBAL_300M_10DAILY_V2', 'CLMS_WB_GLOBAL_100M_MONTHLY_V1', 'CLMS_WB_GLOBAL_300M_MONTHLY_V2', 'CLMS_WB_GLOBAL_300M_10DAILY_V1', 'CLMS_WB_GLOBAL_1KM_10DAILY_V2', 'CLMS_LIE_EUROPE_250M_DAILY_V2', 'CLMS_LIE_BALTIC_250M_DAILY_V1', 'CLMS_LIE_NORTHERNHEMISPHERE_500M_DAILY_V1', 'CLMS_LIE_GLOBAL_500M_DAILY_V2', 'CLMS_LWQ_NRT_GLOBAL_100M_10DAILY_V2', 'CLMS_LWQ_NRT_GLOBAL_100M_10DAILY_V1', 'CLMS_LWQ_NRT_GLOBAL_300M_10DAILY_V2', 'CLMS_LWQ_NRT_GLOBAL_300M_10DAILY_V1', 'CLMS_LWQ_REPROC_GLOBAL_300M_10DAILY_V1', 'CLMS_LC_GLOBAL_100M_YEARLY_V3', 'CLMS_TOC_GLOBAL_300M_DAILY_V2', 'CLMS_LST_GLOBAL_3KM_HOURLY_V3', 'CLMS_LST_DAILY_CYCLE_GLOBAL_3KM_10DAILY_V3']

connection.describe_collection("COPERNICUS_30")

DEM = connection.load_collection("COPERNICUS_30",
    spatial_extent={'west': 13.00, 'east': 13.37, 'south': 11.63, 'north': 12.00,"crs": "EPSG:4326"},
    bands=['DEM'],
)

DEM.download(work_dir+"/Maiduguri_dem.tif")

#import the DEM in ilwispy
DEM = ilwis.RasterCoverage ('Maiduguri_dem.tif')
print(DEM.size())
print(DEM.envelope())
coordSys = DEM.coordinateSystem()
coordSys.toWKT()

Size(1333, 1333, 1)
12.999861 11.629861 13.370139 12.000139

'PROJCS["maiduguri_dem.tif",GEOCS["maiduguri_dem.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION[""],,UNIT[degrees,1.0]]'

#calculatate DEM statistics
stats = DEM.statistics(ilwis.PropertySets.pHISTOGRAM)
#print(stats.histogram())
print('minimum =',(stats[ilwis.PropertySets.pMIN])) # minimum value of the dem
print('maximum =',(stats[ilwis.PropertySets.pMAX])) # maximum value of the dem

minimum = 297.906494140625
maximum = 346.7610168457031

#Transform to a numpy array for visualization using imshow
row = DEM.size().ysize 
col = DEM.size().xsize
dim = DEM.size().linearSize()
print(row, col, dim)

dem_2np = np.fromiter(iter(DEM), np.float64, dim) 
dem_2np = dem_2np.reshape((row, col))

1333 1333 1776889

#create plot
plt.figure(figsize=(8, 8))
plt.imshow(dem_2np, cmap='terrain')
plt.colorbar(shrink=0.5)
plt.title('DEM retrieved from OpenEO service provider');

#note the date - after the flood
t = ["2024-09-17", "2024-09-18"]
s1_cube = connection.load_collection(
  "SENTINEL1_GRD",
  spatial_extent={'west': 13.00, 'east': 13.37, 'south': 11.63, 'north': 12.00,"crs": "EPSG:4326"},
  temporal_extent=t,
  bands=["VV", "VH"]
)

s1_cube_red = s1_cube.median_time()

#Use of the “sigma0-ellipsoid” coefficient for SAR backscattering computation.
s1 = s1_cube.sar_backscatter(coefficient="sigma0-ellipsoid")

#Radar Vegetation Index (RVI) 
rvi = (4 * s1.band("VH")) / (s1.band("VV") + s1.band("VH"))

rvi.download(work_dir+'/RVI.tif')

#import the RVI in ilwispy
RVI = ilwis.RasterCoverage ('RVI.tif')
print(RVI.size())
print(RVI.envelope())
coordSys = RVI.coordinateSystem()
coordSys.toWKT()

Size(4060, 4121, 1)
281940.000000 1286140.000000 322540.000000 1327350.000000

'PROJCS["rvi.tif",GEOCS["rvi.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Transverse_Mercator"],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],PARAMETER["scale",0.9996],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],UNIT[meter,1.0]]'

#Transform to a numpy array for visualization using imshow
row = RVI.size().ysize 
col = RVI.size().xsize
dim = RVI.size().linearSize()
print(row, col, dim)

rvi_2np = np.fromiter(iter(RVI), np.float64, dim) 
rvi_2np = rvi_2np.reshape((row, col))

4121 4060 16731260

#calculatate statistics
stats = RVI.statistics(ilwis.PropertySets.pHISTOGRAM)

Val_minimum = stats[ilwis.PropertySets.pMIN]
Val_maximum = stats[ilwis.PropertySets.pMAX]

plt.figure(figsize=(7, 7))
plt.imshow(rvi_2np,cmap='jet', vmin=Val_minimum, vmax=Val_maximum-1.5)
plt.axis("off")
plt.title('S1 Radar Vegetation Index');

s1_cube_red.download(work_dir+'/S1_after.tif')

s1 = ilwis.RasterCoverage('S1_after.tif')
print(s1.size())
print(s1.envelope())
coordSys = s1.coordinateSystem()
coordSys.toWKT()

Size(4060, 4121, 2)
281940.000000 1286140.000000 322540.000000 1327350.000000

'PROJCS["s1_after.tif",GEOCS["s1_after.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Transverse_Mercator"],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],PARAMETER["scale",0.9996],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],UNIT[meter,1.0]]'

VV = ilwis.do('selection',s1,'rasterbands(0)')
VH = ilwis.do('selection',s1,'rasterbands(1)')
VV.store('VV.mpr')
VH.store('VH.mpr')

#calculate the cross polarized image (VV/VH)
VC = ilwis.do('mapcalc','(@1/@2)', VV, VH)
VC.store('VC.mpr')

#create histogram using large number of bins, given the fact that the input image is a float64, with both positive and negative numbers
hist = VC.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = hist.calcStretchRange(1)
VCs = ilwis.do('linearstretch',VC, minPerc1, maxPerc1)
VCs = ilwis.do('mapcalc','iff(@1==?,255,@1)', VCs) #modify pixels over water
VCs = ilwis.do('setvaluerange', VCs, 0, 255, 1)
VCs.store('VCs.mpr')

VC_2np = np.fromiter(iter(VCs), np.ubyte, VCs.size().linearSize()) 
VC_2np = VC_2np.reshape((VCs.size().ysize, VCs.size().xsize))

#stretch the VV and VH bands
hist = VV.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = hist.calcStretchRange(1)
VVs = ilwis.do('linearstretch',VV, minPerc1, maxPerc1)
VVs = ilwis.do('mapcalc','iff(@1==?,0,@1)', VVs) #remove bad pixels over water
VVs = ilwis.do('setvaluerange', VVs, 0, 255, 1)
VVs.store('VVs.mpr')

VV_2np = np.fromiter(iter(VVs), np.ubyte, VVs.size().linearSize()) 
VV_2np = VV_2np.reshape((VVs.size().ysize, VVs.size().xsize))

hist = VH.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = hist.calcStretchRange(1)
VHs = ilwis.do('linearstretch',VH, minPerc1, maxPerc1)
VHs = ilwis.do('mapcalc','iff(@1==?,0,@1)', VHs) #remove bad pixels over water
VHs = ilwis.do('setvaluerange', VHs, 0, 255, 1)
VHs.store('VHs.mpr')

VH_2np = np.fromiter(iter(VHs), np.ubyte, VHs.size().linearSize()) 
VH_2np = VH_2np.reshape((VHs.size().ysize, VHs.size().xsize))

# create  color composite (in RGB)
cc = np.dstack((VV_2np, VH_2np, VC_2np))

#band interleaved
S1ALL = np.array([VC_2np, VH_2np, VV_2np])

plt.figure(figsize=(7, 7))
plt.imshow(cc)
plt.axis("off")
plt.title('S1 polarization composite');

#create empty ilwis raster 
S1_ilw = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1))
S1_ilw.setDataDef(defNumr)
S1_ilw.setSize(ilwis.Size(4121, 4060, 3)) 
S1_ilw.setGeoReference(s1.geoReference())

#read the numpy array
S1_ilw.array2raster(S1ALL.flatten())

#comment if you don't want to store the image in an ILWIS maplist format
S1_ilw.store('S1_after_processed.mpl')

#uncomment if you want to store the image in geotif format
#S1_ilw.store("S1_after.tif", "GTiff", "gdal")

#connection.describe_collection("SENTINEL2_L2A")
connection.describe_collection("SENTINEL2_L1C")

#note the date - after the flood
t = ["2024-09-10", "2024-09-15"]
s2_cube = connection.load_collection("SENTINEL2_L1C",
    spatial_extent={'west': 13.00, 'east': 13.37, 'south': 11.63, 'north': 12.00,"crs": "EPSG:4326"},
    temporal_extent= t,
    bands=['B02', 'B03', 'B04', 'B08'],
    max_cloud_cover=100,
)

s2_cube.download(work_dir+"/s2_directafter.tif")

#read the S2 image in ilwispy
s2 = ilwis.RasterCoverage ('s2_directafter.tif')
print(s2.size())

Size(4060, 4121, 4)

#note index starts from 0
Blue_temp = ilwis.do('selection',s2,"rasterbands(0)")
Green_temp = ilwis.do('selection',s2,"rasterbands(1)")
Red_temp = ilwis.do('selection',s2,"rasterbands(2)")
IR_temp = ilwis.do('selection',s2,"rasterbands(3)")

Blue = ilwis.do('mapcalc','iff(@1==?, 0, @1)', Blue_temp)
Green = ilwis.do('mapcalc','iff(@1==?, 0, @1)', Green_temp)
Red = ilwis.do('mapcalc','iff(@1==?, 0, @1)', Red_temp)
IR = ilwis.do('mapcalc','iff(@1==?, 0, @1)', IR_temp)

# calculate statistics, retrieve cut-off at lower and upper tail of 1%, conduct linear stretch and ensure that all data is transformed to a byte range
# note that this can also be done in a loop, but here each band is handled seperately
stat_Blue = Blue.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = stat_Blue.calcStretchRange(1)
Blues = ilwis.do('linearstretch',Blue, minPerc1, maxPerc1)
Blues = ilwis.do('setvaluerange', Blues, 0, 255, 1)

stat_Green = Green.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = stat_Green.calcStretchRange(1)
Greens = ilwis.do('linearstretch',Green, minPerc1, maxPerc1)
Greens = ilwis.do('setvaluerange', Greens, 0, 255, 1)

stat_Red = Red.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = stat_Red.calcStretchRange(1)
Reds = ilwis.do('linearstretch',Red, minPerc1, maxPerc1)
Reds = ilwis.do('setvaluerange', Reds, 0, 255, 1)

stat_IR = IR.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)
minPerc1, maxPerc1 = stat_IR.calcStretchRange(1)
IRs = ilwis.do('linearstretch',IR, minPerc1, maxPerc1)
IRs = ilwis.do('setvaluerange', IRs, 0, 255, 1)

#Transform to a numpy array for visualization using imshow
row = IRs.size().ysize 
col = IRs.size().xsize
dim = IRs.size().linearSize()
print(row, col, dim)

Blues_2np = np.fromiter(iter(Blues), np.ubyte, dim) 
Blues_2np = Blues_2np.reshape((row, col))

Greens_2np = np.fromiter(iter(Greens), np.ubyte, dim) 
Greens_2np = Greens_2np.reshape((row, col))

Reds_2np = np.fromiter(iter(Reds), np.ubyte, dim) 
Reds_2np = Reds_2np.reshape((row, col))

IRs_2np = np.fromiter(iter(IRs), np.ubyte, dim) 
IRs_2np = IRs_2np.reshape((row, col))

4121 4060 16731260

# create  color composite (in RGB)
#pixel interleaved
fc_S2 = np.dstack((IRs_2np, Reds_2np, Greens_2np))
nc_S2 = np.dstack((Reds_2np, Greens_2np, Blues_2np))

#band interleaved
S2ALL = np.array([Blues_2np, Greens_2np, Reds_2np, IRs_2np])

fig1 = plt.figure(figsize=(15, 10))

plt.subplot(1, 2, 1)
plt.imshow(nc_S2)
plt.axis("off")
plt.title('NC - S2 image retrieved from OpenEO service provider')

plt.subplot(1, 2, 2)
plt.imshow(fc_S2)
plt.axis("off")
plt.title('FC - S2 image retrieved from OpenEO service provider');

#create empty ilwis raster 
S2_ilw = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1))
S2_ilw.setDataDef(defNumr)
S2_ilw.setSize(ilwis.Size(4121, 4060, 4)) 
S2_ilw.setGeoReference(s2.geoReference())

#read the numpy array
S2_ilw.array2raster(S2ALL.flatten())

#comment if you don't want to store the image in an ILWIS maplist format
S2_ilw.store('S2_directafter_processed.mpl')

#uncomment if you want to store the image in geotif format
#S2_ilw.store("S2_after_processed.tif", "GTiff", "gdal")

#calculate the ndvi
Red = ilwis.do('selection',s2,"rasterbands(2)")
IR = ilwis.do('selection',s2,"rasterbands(3)")

ndvi = ilwis.do('mapcalc','(@1 - @2)/(@1 + @2)', IR, Red)
ndvi = ilwis.do('setvaluerange', ndvi, -1, 1, 0.001)

#check size and value range
print(ndvi.size())
ndv_min = ndvi.min()
ndv_max =ndvi.max()
print(ndv_min, ndv_max)

Size(4060, 4121, 1)
-0.356 0.8190000000000001

ndvi.store('S2_after_ndvi.mpr')

#convert to numpy array for visualization
ndvi_2np = np.fromiter(iter(ndvi), np.float64, ndvi.size().linearSize()) 
ndvi_2np = ndvi_2np.reshape((ndvi.size().ysize, ndvi.size().xsize))

# Keep only negative NDVI values (representing water), mask the rest
ndvi_water = np.ma.masked_where(ndvi_2np >= 0, ndvi_2np)

fig2 = plt.figure(figsize=(15, 10))

plt.subplot(1, 2, 1)
plt.imshow(ndvi_2np, interpolation ='nearest', vmin=0, vmax=ndv_max, cmap='turbo')
plt.colorbar(shrink=0.5, label = 'NDVI')
plt.title('NDVI derived from Sentinel-2 L1C, Maiduguri town and surrounding area', fontsize = 10)
plt.axis("off")

plt.subplot(1, 2, 2)
plt.imshow(ndvi_water, cmap='Blues_r', vmin=ndv_min, vmax=0)
plt.colorbar(label="NDVI (water range only)", shrink=0.5)
plt.title('Water areas (NDVI < 0) - note clouded / shaded areas!!', fontsize = 10)
plt.axis("off")

plt.show()

#once more import the Copernicus DEM in ilwispy
copdem = ilwis.RasterCoverage ('Maiduguri_dem.tif')
print(copdem.size())
print(copdem.envelope())
coordSys = copdem.coordinateSystem()
coordSys.toWKT()

Size(1333, 1333, 1)
12.999861 11.629861 13.370139 12.000139

'PROJCS["maiduguri_dem.tif",GEOCS["maiduguri_dem.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION[""],,UNIT[degrees,1.0]]'

# target georeference (the grid you want to match)
s2_grf = s2.geoReference()

# choose resampling method: nearestneighbour | bilinear | bicubic
copdem_res = ilwis.do("resample", copdem, s2_grf, "nearestneighbour")

# sub map creation
copdem_ressub = ilwis.do('selection',copdem_res,'boundingbox(50 50, 3999 3999)')
print(copdem_ressub.size().xsize)
print(copdem_ressub.size().ysize)
print(copdem_ressub.size().zsize)

3950
3950
1

copdem_ressub.store("copdem_ressub.mpr")

stats = copdem_ressub.statistics(ilwis.PropertySets.pHISTOGRAM)
#print(stats.histogram())
print('minimum =',(stats[ilwis.PropertySets.pMIN])) # minimum value of the dem
print('maximum =',(stats[ilwis.PropertySets.pMAX])) # maximum value of the dem

minimum = 297.906494140625
maximum = 346.1675109863281

#uncomment line below if running the notebook for the first time
fillsink = ilwis.do('MapFillSinks', copdem_ressub, 'cut')
fillsink = ilwis.do('MapFillSinks', fillsink, 'cut')
fillsink = ilwis.do('MapFillSinks', fillsink, 'cut')
fillsink = ilwis.do('MapFillSinks', fillsink, 'fill')

#uncomment line below if running the notebook for the first time and the fill sink operation has been executed
fillsink.store('copdem_fill.mpr')

#uncomment if you want to store the image in geotif format
#fillsink.store("fillsink.tif", "GTiff", "gdal")

dem_fd = ilwis.do('MapFlowDirection', fillsink, 'slope', 'yes')

dem_fd.store('copdem_fd.mpr')

dem_fa = ilwis.do('MapFlowAccumulation', dem_fd)

dem_fa.store('copdem_fa.mpr')

#drainage raster map creation using single flow accumulation threshold - just downstream of dam breach location 
breach_inlet = ilwis.do('ExtractDrainageUseThresholdValue', dem_fa, 3051741)

breach_inlet.store('br_inlet.mpr')

#uncomment if you want to store the image in geotif format
#breach_inlet.store("br_inlet.tif", "GTiff", "gdal")

#read the fabdem in ilwispy - note the extent from the file name (tile of 1 by 1 degree)
fabdem = ilwis.RasterCoverage ('N11E013_FABDEM_V1-2.tif')
print(fabdem.size())

Size(3600, 3600, 1)

#convert to numpy array for visualization
fabdem_2np = np.fromiter(iter(fabdem), np.float64, fabdem.size().linearSize()) 
fabdem_2np = fabdem_2np.reshape((fabdem.size().ysize, fabdem.size().xsize))

#create plot - note area of interest is in upper left hand corner
plt.figure(figsize=(8, 8))
plt.imshow(fabdem_2np, cmap='terrain', vmin=300, vmax=350)
plt.colorbar(shrink=0.5)
plt.title('FABDEM provided by the University of Bristol');

# target georeference (the grid you want to match)
s2_grf = s2.geoReference()

# choose resampling method: nearestneighbour | bilinear | bicubic
fabdem_res = ilwis.do("resample", fabdem, s2_grf, "nearestneighbour")

# sub map creation
fabdem_ressub = ilwis.do('selection',fabdem_res,'boundingbox(50 50, 3999 3999)')
print(fabdem_ressub.size().xsize)
print(fabdem_ressub.size().ysize)
print(fabdem_ressub.size().zsize)

3950
3950
1

#convert to numpy array for visualization
fabdemressub_2np = np.fromiter(iter(fabdem_ressub), np.float64, fabdem_ressub.size().linearSize()) 
fabdemressub_2np = fabdemressub_2np.reshape((fabdem_ressub.size().ysize, fabdem_ressub.size().xsize))

#create plot
plt.figure(figsize=(8, 8))
plt.imshow(fabdemressub_2np, cmap='terrain')
plt.colorbar(shrink=0.5)
plt.title('Resampled FABDEM');

fabdem_ressub.store("fabdem_ressub.mpr")

stats1 = fabdem_ressub.statistics(ilwis.PropertySets.pHISTOGRAM)
#print(stats.histogram())
print('minimum =',(stats1[ilwis.PropertySets.pMIN])) # minimum value of the dem
print('maximum =',(stats1[ilwis.PropertySets.pMAX])) # maximum value of the dem

minimum = 300.0299987792969
maximum = 344.489990234375

difference = ilwis.do('mapcalc','@1-@2', fabdem_ressub, copdem_ressub)

#convert to numpy array for visualization
difference_2np = np.fromiter(iter(difference), np.float64, difference.size().linearSize()) 
difference_2np = difference_2np.reshape((difference.size().ysize, difference.size().xsize))

#create plot
plt.figure(figsize=(8, 8))
plt.imshow(difference_2np, cmap='jet', vmin=-3, vmax=3)
plt.colorbar(shrink=0.5)
plt.title('difference = FABDEM - COPDEM');

drconditioned = ilwis.do('mapcalc','iff(@1==1,@2-1,@2)', breach_inlet, fabdem_ressub)

#uncomment line below if running the notebook for the first time
#cut operation only done once since urban artefacts are already removed
fillsink = ilwis.do('MapFillSinks', drconditioned, 'cut')
fillsink = ilwis.do('MapFillSinks', fillsink, 'fill')

fillsink.store("fabdem_fill_dr_cond.mpr")

#import dem and selected river reach based on flow accumulation threshold and create visualization

# ---initialize size and dimensions---
row = fillsink.size().ysize 
col = fillsink.size().xsize
dim = fillsink.size().linearSize()

# --- DEM read raster to numpy array ---
dem = np.fromiter(iter(fillsink), np.float64, dim) 
dem = dem.reshape((row, col))
rows, cols = dem.shape

# --- process the river reach and create a raster ---
river_raster = np.fromiter(iter(breach_inlet), np.float64, dim) 
river_raster = river_raster.reshape((row, col))
river_mask = river_raster > 0

print("River pixels:", np.sum(river_mask))


# --- dilate river to make it wider ---
river_mask_wide = binary_dilation(river_mask, iterations=4)  # Approximate river width (pixels) = 1 + 2 * iterations

# --- create river water surface for plotting ---
water_surface_river = np.full_like(dem, np.nan)
water_surface_river[river_mask_wide] = river_raster[river_mask_wide] 

# plot the figure
plt.figure(figsize=(8,8))

# 1. DEM
plt.imshow(dem, cmap='terrain')
plt.colorbar(label='DEM Elevation (m)', shrink=0.5)

# 2. Wider river for visibility
plt.imshow(np.ma.masked_invalid(water_surface_river), cmap='gray')

plt.title("DEM with selected river section")
plt.axis('off')
plt.show();

River pixels: 5371

# Input paramter settings
pixel_size = 10               # DEM resolution (m)

#Adjust accoring to your requirements
water_depth_river = 1.50      # water above river (m) - initial value 1.5 
spread_decay = 0.00085        # friction across floodplain per pixel - lower friction results in further spread -initial value 0.00085
max_spread_distance = 2000    # max flood distance (m) from river - initial value 2000

#Preparation of the other required input

# --- retrieve georeference and coordinate system info ---
georef = fillsink.geoReference()

crs = fillsink.coordinateSystem()
print(crs.toWKT())

# --- raster size ---
ny = fillsink.size().ysize
nx = fillsink.size().xsize

PROJCS["s2_directafter.tif",GEOCS["s2_directafter.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Transverse_Mercator"],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],PARAMETER["scale",0.9996],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],UNIT[meter,1.0]]

#Initial water surface on river, set up of the starting water level inside the river network.
#It creates the initial water surface that will later spread across the floodplain. The code creates a raster of water 
#surface elevation and fills it only at river locations using DEM elevation + assumed river water depth. The first line 
#creates a new raster having the same size and grid as the DEM and is filled with NaN (no value) everywhere. 
#line 2 assignes the water level to river pixels through means of boolean indexing, if the value is true (river pixel)
#then water surface elevation = terrain elevation + water depth. Note the model assumes that the river everywhere contains
#a uniform water depth above the channel bed!

water_surface = np.full_like(dem, np.nan)
water_surface[river_mask] = dem[river_mask] + water_depth_river

#Distance to river (for the max. spread). This code field computes how far every pixel is from the river
#and limits flooding to a maximum distance acting as a very important geomorphologic constraint in the model.
#(~river_mask) = NOT river, edt uses the Euclidean Distance Transform: The EDT computes, for every non-river pixel the 
#distance to the nearest river pixel. It results in a map representing a kind of distance gradient away from the river. 

distance_pixels = distance_transform_edt(~river_mask)
distance_m = distance_pixels * pixel_size # convert pixels to meters
distance_m[distance_m > max_spread_distance] = np.nan #limit the flood spreading distance

#Breadth-First Search flood propagation from river. It spreads a water surface outward from the river across the 
#DEM until the water runs out of energy. This code section spreads the river water surface across the terrain using 
#a breadth-first search algorithm, reducing water level with distance and stopping where the ground becomes 
#too high or too far from the river.

#This raster will store the final water surface elevation after spreading. At the start all values = NaN.
water_surface_spread = np.full_like(dem, np.nan)

#Initialize the spreading queue (flood sources)
queue = deque()

#np.argwhere(river_mask) returns the coordinates of all river pixels. Each river pixel will act as a 
#starting source of floodwater
river_indices = np.argwhere(river_mask)

#Seed the model with river water. Copy the river water surface into the spreading map and add all river 
#cells to the queue. The queue now contains the initial flood front. 

for r, c in river_indices:
    water_surface_spread[r, c] = water_surface[r, c]
    queue.append((r, c))

#Define 8-neighbour connectivity to approximates 2-D overland flow
neighbors = [(-1,0),(1,0),(0,-1),(0,1),(-1,-1),(-1,1),(1,-1),(1,1)]

#Start BFS flood propagation. The queue works like a moving flood front, take one flooded cell and try
#to spread water to its neighbours and continues until no more cells can be flooded. 

while queue:
    r, c = queue.popleft()
    current_ws = water_surface_spread[r, c]

    for dr, dc in neighbors:
        rr, cc = r + dr, c + dc #Check neighbour cells

        if rr < 0 or rr >= rows or cc < 0 or cc >= cols:
            continue

        # skip if too far from river
        if distance_m[rr, cc] > max_spread_distance: #Limit flood distance from river
            continue

        #apply friction decay, representing the hydraulic energy loss on the floodplain
        #each step away from the river the water surface drops slightly (to mimics friction
        #and shallow flow resistance). So the flood surface gradually slopes downward.
       
        new_ws = current_ws - spread_decay * pixel_size #Apply floodplain friction loss

        #Check terrain blocking, stop if water too shallow. Water can only flood a cell if
        #water surface>ground elevation. This lets the DEM control ridges to block water, depressions fill
        #and valleys to guide the flow This is the key terrain constraint.
        
        if new_ws <= dem[rr, cc]:
            continue

        #Update if higher water level. It ensures the model always keeps the highest possible water surface 
        #at each cell.This allows multiple flow paths, gradual correction and equilibrium solution, making
        #the algorithm diffusive and stable.
        
        if np.isnan(water_surface_spread[rr, cc]) or new_ws > water_surface_spread[rr, cc]:
            
            #continue spreading, newly flooded cell becomes a new source. The flood keeps expanding until
            #it runs out of energy.The loop ends /  the queue becomes empty when water cannot climb terrain anymore
            #or friction removed too much energy or max river distance is reached. At that moment the flood reaches 
            #equilibrium obtaining the final water surface elevation everywhere the flood can reach. 
            water_surface_spread[rr, cc] = new_ws
            queue.append((rr, cc))

# Create flood mask and depth
flood_mask = water_surface_spread > dem
flood_depth = np.where(flood_mask, water_surface_spread - dem, 0)

print("Flooded pixels:", np.sum(flood_mask))

Flooded pixels: 3275111

# Visualization

plt.figure(figsize=(12,6))

plt.subplot(1,2,1)
plt.title("Flood Mask")
plt.imshow(flood_mask, cmap="Blues")
plt.colorbar(shrink=0.5)

plt.subplot(1,2,2)
plt.title("Flood Depth")
plt.imshow(flood_depth/2, cmap="Blues", vmin=0, vmax=2.5)
plt.colorbar(shrink=0.5)

plt.tight_layout()
plt.show()

# --- dilate river to make it wider ---
river_mask_wide = binary_dilation(river_mask, iterations=4)  # Approximate river width (pixels) = 1 + 2 * iterations

# --- create river water surface for plotting ---
water_surface_river = np.full_like(dem, np.nan)
water_surface_river[river_mask_wide] = river_raster[river_mask_wide]# + water_depth_river

# --- plot ---
plt.figure(figsize=(8,8))

# 1. DEM
plt.imshow(dem, cmap='terrain')
plt.colorbar(label='DEM Elevation (m)', shrink=0.5)

# 2. Floodplain water surface
plt.imshow(np.ma.masked_invalid(water_surface_spread), cmap='Blues', alpha=0.75)

# 3. Wider river for visibility
plt.imshow(np.ma.masked_invalid(water_surface_river), cmap='gray')

plt.title("DEM with water surface mask and selected river section")
plt.axis('off')
plt.show()

#create empty ilwis raster for flood depth map
fl_depth = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 100, 0.001))
fl_depth.setDataDef(defNumr)
fl_depth.setSize(ilwis.Size(fabdem_ressub.size().ysize, fabdem_ressub.size().xsize, 1)) 
fl_depth.setGeoReference(fabdem_ressub.geoReference())

fl_depth.array2raster(flood_depth.flatten())

#adjust depth - as depth is exaggerated - to be further eleaborated upon as no temporal flood wave dynamics
#(hydrograph) is currently implemented
fl_depth_adj = ilwis.do('mapcalc','iff(@1>0,@1/2,?)', fl_depth)

#comment if you don't want to store the image in an ILWIS maplist format
fl_depth_adj.store('M1_flood_depth.mpr')

#uncomment if you want to store the image in geotif format
#fl_depth_adj.store("flood_depth.tif", "GTiff", "gdal")

#create empty ilwis raster for flood mask map
fl_mask = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 1,1))
fl_mask.setDataDef(defNumr)
fl_mask.setSize(ilwis.Size(fabdem_ressub.size().ysize, fabdem_ressub.size().xsize, 1)) 
fl_mask.setGeoReference(fabdem_ressub.geoReference())

fl_mask.array2raster(flood_mask.flatten())

#comment if you don't want to store the image in an ILWIS maplist format
fl_mask.store('M1_flood_mask.mpr')

#uncomment if you want to store the image in geotif format
#fl_mask.store("flood_mask.tif", "GTiff", "gdal")

# (Hydraulic) input parameter settings

# --- raster pixel size (in meters) of your resampled DEM ---
cellsize = 10

inflows = [(3005, 2976)]  #note index offset as python starts from 0
                         

# --- other hydraulic parameter settings --- Adjust upon your requirement
flood_depth_inlet = 2.25
longitudinal_loss = 0.000285 #m/m
lateral_decay = 0.00085 #m/m
bankfull_depth = 1.85  # bankfull channel depth in m

#Preparation of required input raster data and transformation of sink-filled DEM 
#and river raster (drainage line starting at inlet) to numpy arrays
#load the original datasets once more

# --- DEM check size and dimensions---
row = fillsink.size().ysize 
col = fillsink.size().xsize
dim = fillsink.size().linearSize()
print(row, col, dim)

# --- DEM read raster to numpy array ---
dem = np.fromiter(iter(fillsink), np.float64, dim) 
dem = dem.reshape((row, col))

# --- retrieve georeference and coordinate system info ---
georef = fillsink.geoReference()

crs = fillsink.coordinateSystem()
print(crs.toWKT())

# --- raster size ---
ny = fillsink.size().ysize
nx = fillsink.size().xsize

# --- RIVER MASK read ilwis raster to numpy array ---
river_raster = np.fromiter(iter(breach_inlet), np.float64, dim) 
river_raster = river_raster.reshape((row, col))

river = np.array(river_raster) == 1

3950 3950 15602500
PROJCS["s2_directafter.tif",GEOCS["s2_directafter.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Transverse_Mercator"],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],PARAMETER["scale",0.9996],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],UNIT[meter,1.0]]

# Creating the model state grids — the maps that will store the simulation results while the flood propagates 
# At the start of the model for all pixels: river_stage = NaN, 	flood_depth = 0 and flood_mask = False

river_stage = np.full_like(dem, np.nan, dtype=np.float32) # River routing fills river_stage
flood_depth = np.zeros_like(dem, dtype=np.float32) # Bank spill starts filling flood_depth
flood_mask = np.zeros_like(dem, dtype=bool) #Floodplain spreading fills flood_mask

# River routing, computes the water surface profile along the river network before any flooding happens.In hydraulic terms it represents
# a simplified steady 1-D river routing step. Spread water downstream through river cells while gradually losing energy.
# The algorithm used is Breadth-First Search (BFS) on the river network.

#Create routing memory
visited = np.zeros_like(dem, dtype=bool)
queue = deque()

# Insert inflow boundary conditions
# At each inflow point: Compute starting water surface elevation (stage=ground elevation+water depth)
# Store this in the river map. Put the cell into the routing queue.
# This defines the upstream boundary condition. Now the queue contains the first water cells.

for r,c in inflows:
    stage = dem[r,c] + flood_depth_inlet
    river_stage[r,c] = stage
    queue.append((r,c,stage))
    visited[r,c] = True

# While there are river cells with water to process: take the next cell from the queue. Try to send water to neighbouring river cells, 
# water moves cell by cell through the network.
while queue:
    r,c,stage = queue.popleft()

#Explore neighbouring cells (D8 connectivity), skip the centre cell:
    for dr in [-1,0,1]:
        for dc in [-1,0,1]:
            if dr==0 and dc==0:
                continue

# Check if neighbour is valid. Water can only move if the neighbour is inside the raster, is a river cell and has not already been reached.
# So water stays inside the river network.
            nr,nc = r+dr, c+dc
            if not (0<=nr<ny and 0<=nc<nx):
                continue
            if not river[nr,nc]:
                continue
            if visited[nr,nc]:
                continue

# Compute hydraulic energy loss. This is the physics of the model. Water loses elevation as it moves downstream: 
# new water level=old level−(slope×distance), Where: distance = cell size (or diagonal length), longitudinal_loss = hydraulic gradient (m/m)
# This mimics: friction, bed slope, energy dissipation. So the river water surface tilts downstream.
            
            dist = np.hypot(dr,dc)*cellsize
            new_stage = stage - longitudinal_loss*dist

# Prevent uphill river flow. Water cannot flow if the water surface would be below the ground. This prevents impossible uphill routing.
            
            min_stage = dem[nr,nc]
            if new_stage < min_stage:
                continue

# Accept the new river cell. If all checks pass: store water level in the new river cell, mark it visited and add it to the queue to continue routing
# The flood wave keeps moving.
            
            river_stage[nr,nc] = new_stage
            visited[nr,nc] = True
            queue.append((nr,nc,new_stage))

# Bankfull spill - bankfull exceedance / overbank flow. Computation of the downstream river water surface profile by routing water through river cells 
# while applying a hydraulic slope (energy loss). Hydraulically, this is the channel and floodplain coupling step.
# Up to now the model only knows the water level inside the river. This block checks: Where does the river overflow its banks?
# This code section finds where the river water level exceeds the bank height and converts those river cells into sources of floodplain flooding.
# This step simulates: channel capacity, levee overtopping, bank overflow, start of overbank flooding. 

# Create a new queue for flooding. This queue will store starting points of floodplain flooding. 
# River routing queue: when water moving in channel. Spill queue: when water moving onto land

spill_queue = deque()

# Loop over all river cells. np.where(river) finds every cell where the river mask = True. So we now inspect every river pixel in the map.

for r,c in zip(*np.where(river)):

# Skip river cells without water. If the routing step never reached this river cell, it has no water, it cannot spill, thenignore it.
    if np.isnan(river_stage[r,c]):
        continue
        
# Compute bankfull elevation, his defines the top of the river banks (bankfull elevation=ground level+channel depth)
# Example: ground = 12 m, bankfull_depth = 1.5 m, banks at 13.5 m. This is the maximum water level the channel can hold.
    
    bankfull_elev = dem[r,c] + bankfull_depth

# Check if the river overtops, compare: river water surface and bank height. If water level is below banks nothing happens.
# If water level is above banks then flooding starts.

    spill_depth = river_stage[r,c] - bankfull_elev

    if spill_depth <= 0:
        continue

# Mark the river cell as flooded, Once overtopping occurs, the river cell becomes part of the flooded terrain.
# Flood depth = water surface − ground elevation. So the river becomes the first flooded area.
    
    flood_mask[r,c] = True
    flood_depth[r,c] = river_stage[r,c] - dem[r,c]

# Add spill locations to flood queue. These cells are now the starting points of floodplain spreading.
# These act like “holes in the river banks” where water pours onto land.
    
    spill_queue.append((r,c,river_stage[r,c]))

# Lateral floodplain propagation which is spreading floodwater across the terrain from overtopped river cells, losing energy with distance 
# and stopping where the ground becomes too high.This is the actual floodplain flooding step — where water spreads over land after leaving the river.
# Hydraulically this represents 2-D overland flood spreading with friction losses. If Step 1 = river flow and Step 2 = river overflow
# then in the code field below = water is moving across the terrain. Spread water from spill locations across the DEM while losing energy
# and respecting terrain elevations represents basically a 2-D hydraulic fill model using the Breadth-First Search (BFS) algorithm.
# When the queue becomes empty it has created the full flood extent, flood depth and floodplain water surface.
# The flood stops when water loses too much energy to climb terrain. This step mimics sheet flow, overbank spreading, ponding in depressions
# and floodplain hydraulics, it is similar to a simplified 2-D shallow water model, but much faster.


# Create floodplain visitation map. Just like in river routing, it tracks which floodplain cells have already been processed.
# This prevents: infinite loops, repeated calculations and water bouncing back and forth

visited_fp = np.zeros_like(dem, dtype=bool)

# Start flood spreading from spill cells. The previous step added overtopped river cells to this queue. Now each of those cells becomes
# a source of overland flow and Water spreads outward.

while spill_queue:
    r,c,stage = spill_queue.popleft()

# Explore neighbouring terrain cells (D8). Water can spread in all 8 directions. This is a simplified 2-D hydraulic connectivity.
    for dr in [-1,0,1]:
        for dc in [-1,0,1]:
            if dr==0 and dc==0:
                continue

#Skip invalid neighbours, skip cells that are outside the raster and already received flood water
            
            nr,nc = r+dr, c+dc
            if not (0<=nr<ny and 0<=nc<nx):
                continue
            if visited_fp[nr,nc]:
                continue
                
# Compute floodplain energy loss. Water loses elevation due to vegetation, buildings, rough terrain, shallow flow friction, etc.
# so the flood surface gradually slopes downward away from the river.
            
            dist = np.hypot(dr,dc)*cellsize
            new_stage = stage - lateral_decay*dist

# Check terrain blocking. Water can only flood a cell if the water surface is above the ground. This automatically creates natural levees, 
# terraces, ridges, flow barriers, etc. In this way the DEM controls where water can go.
            
            if new_stage <= dem[nr,nc]:
                continue

# Mark the cell as flooded. Compute flood depth on land: depth = water_surface − ground. This creates the final flood map. 
            
            flood_mask[nr,nc] = True
            flood_depth[nr,nc] = new_stage - dem[nr,nc]

# Continue flood propagation. The flooded cell becomes a new source of flooding. This is how the flood wave expands outward.
            
            visited_fp[nr,nc] = True
            spill_queue.append((nr,nc,new_stage))

# visualization
print("Flooded pixels:", np.sum(flood_mask))

plt.figure(figsize=(12,6))

plt.subplot(1,2,1)
plt.title("Flood Mask")
plt.imshow(flood_mask, cmap="Blues")
plt.colorbar(shrink=0.5)

plt.subplot(1,2,2)
plt.title("Flood Depth")
plt.imshow(flood_depth/2, cmap="Blues", vmin=0, vmax=2.5)
plt.colorbar(shrink=0.5)

plt.tight_layout()
plt.show()

Flooded pixels: 2740069

# Dilate river to make it wider (3 pixels)
river_mask_wide = binary_dilation(river, iterations=4)  # Approximate river width (pixels) = 1 + 2 * iterations

# Create river water surface for plotting
water_surface_river = np.full_like(dem, np.nan)
water_surface_river[river_mask_wide] = river[river_mask_wide]

# Plot
plt.figure(figsize=(10,8))

# 1. DEM
plt.imshow(dem, cmap='terrain')
plt.colorbar(label='DEM Elevation (m)', shrink=0.5)

# 2. Flood mask
mask_masked = np.ma.masked_where(flood_mask <= 0, flood_mask)
plt.imshow(mask_masked, cmap='Blues', alpha=0.5)


# 3. Wider river for visibility
plt.imshow(np.ma.masked_invalid(water_surface_river), cmap='gray')

plt.title("DEM with flood mask and selected river section")
plt.axis('off')
plt.show()

#create empty ilwis raster for flood depth map
fl_depth = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 100, 0.001))
fl_depth.setDataDef(defNumr)
fl_depth.setSize(ilwis.Size(fillsink.size().ysize, fillsink.size().xsize, 1)) 
fl_depth.setGeoReference(fillsink.geoReference())

fl_depth.array2raster(flood_depth.flatten())

#adjust depth - as depth is exaggerated - to be further eleaborated upon as no temporal flood wave dynamics (hydrograph) is implemented
fl_depth_adj = ilwis.do('mapcalc','iff(@1>0,@1/4,?)', fl_depth)

#comment if you don't want to store the image in an ILWIS format
fl_depth_adj.store('M2_flood_depth.mpr')

#uncomment if you want to store the image in geotif format
#fl_depth_adj.store("flood_depth.tif", "GTiff", "gdal")

#create empty ilwis raster for flood mask map
fl_mask = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 1,1))
fl_mask.setDataDef(defNumr)
fl_mask.setSize(ilwis.Size(fillsink.size().ysize, fillsink.size().xsize, 1)) 
fl_mask.setGeoReference(fillsink.geoReference())

fl_mask.array2raster(flood_mask.flatten())

#comment if you don't want to store the image in an ILWIS format
fl_mask.store('M2_flood_mask.mpr')

#uncomment if you want to store the image in geotif format
#fl_mask.store("flood_mask.tif", "GTiff", "gdal")

#Mask dry pixels
depth = flood_depth.copy()/4
depth[depth <= 0.05] = np.nan

# Choose depth color scaling
max_depth = (np.nanmax(depth))  # adjust to your flood
norm_depth = np.clip((depth) / max_depth, 0, 1)
print(max_depth)

2.045902

# Convert to RGBA overlay,  create a blue colormap manually. deep water = darker + more opaque, 
# shallow water = lighter + transparent. Dry pixels fully transparent.

overlay = np.zeros((*depth.shape, 4))

# Blue color
overlay[..., 2] = 1.0

# Transparency based on depth
overlay[..., 3] = np.nan_to_num(norm_depth) * 0.6

# visualization
plt.figure(figsize=(12,6))

plt.subplot(1,2,1)
plt.imshow(fc_S2)

im = plt.imshow(depth, cmap="Blues", alpha=0.45, vmin=0, vmax=max_depth)
plt.title("Flood depth mask over Sentinel-2 FCC", fontsize = 10)
plt.axis("off")

plt.subplot(1,2,2)
plt.imshow(ndvi_water, cmap='Blues_r', vmin=ndv_min, vmax=0)
im = plt.imshow(depth, cmap="Blues", alpha=0.65, vmin=0, vmax=max_depth)
plt.title('Water areas (NDVI < 0) - note clouded / shaded areas and flood depth mask', fontsize = 10)
plt.axis("off")

plt.tight_layout()
plt.show()

Parameter	Units	Physical meaning
Inflow locations	raster index	upstream boundary
Initial flood depth	m	starting water depth
Longitudinal loss	m/m	downstream energy slope
Lateral decay	m/m	floodplain friction slope
Bankfull depth	m	channel capacity

Model 1	Model 2
River = flood source everywhere	River → spill → floodplain
Assumes river already flooded	Simulates bank overtopping
No river routing	Has river routing (longitudinal hydraulics)
Single-stage spreading	Two-stage hydraulics
More GIS / cost-distance style	Physically closer to hydraulics

Dam breach modeling using ILWISPy in conjunction with Python¶

Download and pre-processing data¶

Other OpenEO data providers / data repositories¶

Selection of images¶

Retrieve Sentinel 1 GRD image from after the flood event¶

Create radar RGB composite¶

Retrieve Sentinel 2 level 2A image from after the flood event¶

Retrieval of NDVI¶

Some further Hydro-DEM preprocessing¶

Fill Sinks (routine requires some time)¶

Flow Direction¶

Flow Accumulation¶

Drainage Network extraction - define single threshold¶

FABDEM¶

Some modeling using the data pre-processed¶

Model 1: A simplified flood water spreading model¶

Adjust the following parameters to change the flood extent:¶

Store the results obained¶

Model 2: Drainage controlled flow model - a raster hydraulic routing & floodplain propagation model¶

Model input parameters (what drives the simulation)¶

Initial values given to the hydraulic parameters¶

Store the results obtained¶

Plot flood depth map on top of Sentinel-2 FCC image and NDVI¶

Concluding remarks:¶

Step	Physics
River routing	1-D channel hydraulics
Bankfull spill	river–floodplain coupling
Floodplain propagation	2-D overland flow