{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "144c0f32-91c6-447e-8171-e7ee6d15be1b",
   "metadata": {},
   "source": [
    "# Dam breach modeling using ILWISPy in conjunction with Python"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13045c13-4d10-4230-9dd2-3020ea2a7f13",
   "metadata": {},
   "source": [
    "Notebook prepared by Ben Maathuis. ITC-University of Twente, Enschede. The Netherlands"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb6851d7-5a50-4802-a773-aeee0f07c807",
   "metadata": {},
   "source": [
    "Within this Notebook your are going to use ILWISPy in conjunction with a number of common used Python libraries (like Numpy and Matplotlib) as well as OpenEO for initial data retrieval. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6185be4-0b5f-48ad-9df1-d90f93917729",
   "metadata": {},
   "source": [
    "The case introduced here is a flood resulting from a dam breach. In late August 2024, northeastern Nigeria was struck by intense rainfall, leading to severe flooding. The city of Maiduguri in Borno State (North-East of Nigeria) was particularly affected when the Alau Dam, situated 20 kilometres south on the Ngadda River, collapsed on 10 September. This disaster has impacted a large number of people, especially those living in Maiduguri city."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed281c91-bb95-4978-b983-e16bfa02ddb8",
   "metadata": {},
   "source": [
    "First we are going to download an elevation model (the Copernicus 30 meter DEM) as well as a post disaster Sentinel-2 image. The DEM is used as flood model input and the Sentinel-2 image for model validation. Two flood spreading models are implemented. Both models create a flood that spreads pixel by pixel from river pixels using a queue (Breadth-First Search) together with friction decay coeficients. The flood stops when the resulting water_surface <= DEM.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c47c931-3288-4b56-b525-77204bcc01e7",
   "metadata": {},
   "source": [
    "The models starts at the dam breach location. An initial set of parameters is provided, these can be modified and the model can be executed using these new values."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e01b103f-d0b1-40b9-bfd3-b7934fc54f21",
   "metadata": {},
   "source": [
    "### Download and pre-processing data\n",
    "\n",
    "Here use is made of the Copernicus DataSpace Ecosystem (https://dataspace.copernicus.eu/). Before you continue register and create an account (for free). Also download the FABDEM data used, available at: https://filetransfer.itc.nl/pub/52n/ilwis_py/sample_data_V2/flood_modelling.zip. Unzip the file, note that the folder 'flood_modelling' created under this notebook is the output folder for all data that is subsequently downloaded / processed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "030661a2-9791-4b84-9551-12d22077bff8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import required libraries - eventually install if not available \n",
    "import os\n",
    "import ilwis\n",
    "import openeo\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import deque\n",
    "from scipy.ndimage import binary_dilation\n",
    "from scipy.ndimage import distance_transform_edt\n",
    "\n",
    "#suppress warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f91300e3-977c-4541-ba85-adb8fea0334b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#within the notebook folder a directory is created to store your processing results\n",
    "work_dir = (os.getcwd())+ r'/flood_modelling'\n",
    "\n",
    "print(\"current dir is: %s\" % (os.getcwd()))\n",
    "print(\"current working directory is:\",work_dir) \n",
    "\n",
    "if os.path.isdir(work_dir):\n",
    "    print(\"Folder exists\")\n",
    "else:\n",
    "    print(\"Folder doesn't exists\")\n",
    "    os.mkdir(work_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3512453-d4b1-40a6-a329-0d7e84773f75",
   "metadata": {},
   "outputs": [],
   "source": [
    "#set the working directory for ILWISPy\n",
    "ilwis.setWorkingCatalog(work_dir)\n",
    "print(work_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ec744dd-7a7d-4837-8b19-6b1d330eb074",
   "metadata": {},
   "outputs": [],
   "source": [
    "ilwis.version()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecc691d2-31ad-45c1-834b-32836e60e3ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "#remove previous login credentials\n",
    "from openeo.rest.auth.config import RefreshTokenStore\n",
    "RefreshTokenStore().remove()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f08fa67-eebe-4833-841b-fad346ea79a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ensure openeo is installed and you have registered\n",
    "connection = openeo.connect(\"openeo.dataspace.copernicus.eu\").authenticate_oidc()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98710d35-c4fa-447f-af45-57d6b2d9006e",
   "metadata": {},
   "source": [
    "#### Other OpenEO data providers / data repositories\n",
    "Uncomment the lines below if you want to review the data provided by other openEO service providers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2f6ca5f-ef05-433e-ad3d-f4297a56ed78",
   "metadata": {},
   "outputs": [],
   "source": [
    "#connection = openeo.connect(\"openeofed.dataspace.copernicus.eu\").authenticate_oidc()\n",
    "#connection = openeo.connect(\"openeo.vito.be\").authenticate_oidc()\n",
    "#connection = openeo.connect(\"openeo.cloud\").authenticate_oidc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05c2e405-5445-47fc-b1c7-227ceeb8d024",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(connection.list_collection_ids())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "514f813f-2ba9-4109-b2b4-ea240be51749",
   "metadata": {},
   "outputs": [],
   "source": [
    "connection.describe_collection(\"COPERNICUS_30\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d4f7fe8-a079-4987-a64b-7406788dfe8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "DEM = connection.load_collection(\"COPERNICUS_30\",\n",
    "    spatial_extent={'west': 13.00, 'east': 13.37, 'south': 11.63, 'north': 12.00,\"crs\": \"EPSG:4326\"},\n",
    "    bands=['DEM'],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20fa6ae5-f34c-4d99-9c91-392e1fa224fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "DEM.download(work_dir+\"/Maiduguri_dem.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82f96b6e-3cf1-4ee5-81f9-f0b5c8ec102f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import the DEM in ilwispy\n",
    "DEM = ilwis.RasterCoverage ('Maiduguri_dem.tif')\n",
    "print(DEM.size())\n",
    "print(DEM.envelope())\n",
    "coordSys = DEM.coordinateSystem()\n",
    "coordSys.toWKT()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "741a739f-ae0b-4ee3-b351-c5c7e56f73aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculatate DEM statistics\n",
    "stats = DEM.statistics(ilwis.PropertySets.pHISTOGRAM)\n",
    "#print(stats.histogram())\n",
    "print('minimum =',(stats[ilwis.PropertySets.pMIN])) # minimum value of the dem\n",
    "print('maximum =',(stats[ilwis.PropertySets.pMAX])) # maximum value of the dem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3a3ee3b-10e6-41f2-97c3-1f46c88e5981",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transform to a numpy array for visualization using imshow\n",
    "row = DEM.size().ysize \n",
    "col = DEM.size().xsize\n",
    "dim = DEM.size().linearSize()\n",
    "print(row, col, dim)\n",
    "\n",
    "dem_2np = np.fromiter(iter(DEM), np.float64, dim) \n",
    "dem_2np = dem_2np.reshape((row, col))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "973f6a74-1de6-4533-a219-a050ed5d8ea3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create plot\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.imshow(dem_2np, cmap='terrain')\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.title('DEM retrieved from OpenEO service provider');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9476e2b9-2bf6-40f4-b144-5b62c5ff984c",
   "metadata": {},
   "source": [
    "### Selection of images\n",
    "\n",
    "see: https://browser.dataspace.copernicus.eu/ for manual selection of a Sentinel image"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12b8a22c-2b33-4a82-befd-f3f78dbcac0a",
   "metadata": {},
   "source": [
    "### Retrieve Sentinel 1 GRD image from after the flood event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19d70cac-a477-4b2e-a882-09216a88ceea",
   "metadata": {},
   "outputs": [],
   "source": [
    "#note the date - after the flood\n",
    "t = [\"2024-09-17\", \"2024-09-18\"]\n",
    "s1_cube = connection.load_collection(\n",
    "  \"SENTINEL1_GRD\",\n",
    "  spatial_extent={'west': 13.00, 'east': 13.37, 'south': 11.63, 'north': 12.00,\"crs\": \"EPSG:4326\"},\n",
    "  temporal_extent=t,\n",
    "  bands=[\"VV\", \"VH\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e070f06-0520-4c72-ac80-95ff50846a37",
   "metadata": {},
   "outputs": [],
   "source": [
    "s1_cube_red = s1_cube.median_time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6d1bab6-00b2-41e7-81e3-6e389e72fb5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Use of the “sigma0-ellipsoid” coefficient for SAR backscattering computation.\n",
    "s1 = s1_cube.sar_backscatter(coefficient=\"sigma0-ellipsoid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9660bdf2-e059-439f-b3f7-94d7df63e0c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Radar Vegetation Index (RVI) \n",
    "rvi = (4 * s1.band(\"VH\")) / (s1.band(\"VV\") + s1.band(\"VH\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e9a26a7-96ad-487b-9c08-5b0693e898ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "rvi.download(work_dir+'/RVI.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "208c9d61-a373-4379-a088-fa1ab44bc65e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import the RVI in ilwispy\n",
    "RVI = ilwis.RasterCoverage ('RVI.tif')\n",
    "print(RVI.size())\n",
    "print(RVI.envelope())\n",
    "coordSys = RVI.coordinateSystem()\n",
    "coordSys.toWKT()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e34e406-702f-402a-be8c-b6d718fe03a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transform to a numpy array for visualization using imshow\n",
    "row = RVI.size().ysize \n",
    "col = RVI.size().xsize\n",
    "dim = RVI.size().linearSize()\n",
    "print(row, col, dim)\n",
    "\n",
    "rvi_2np = np.fromiter(iter(RVI), np.float64, dim) \n",
    "rvi_2np = rvi_2np.reshape((row, col))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76b48edb-4a98-4e55-b946-b494455b5155",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculatate statistics\n",
    "stats = RVI.statistics(ilwis.PropertySets.pHISTOGRAM)\n",
    "\n",
    "Val_minimum = stats[ilwis.PropertySets.pMIN]\n",
    "Val_maximum = stats[ilwis.PropertySets.pMAX]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe8b6486-9ef4-4c4c-8a69-497c948240ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7, 7))\n",
    "plt.imshow(rvi_2np,cmap='jet', vmin=Val_minimum, vmax=Val_maximum-1.5)\n",
    "plt.axis(\"off\")\n",
    "plt.title('S1 Radar Vegetation Index');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68d19b86-6254-46cd-93a6-dfa399321370",
   "metadata": {},
   "source": [
    "#### Create radar RGB composite \n",
    "+ VV - red\n",
    "+ VH - green\n",
    "+ VC - blue (cross polarized)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c35a19e-b4c8-42b4-98fe-3374681f870e",
   "metadata": {},
   "outputs": [],
   "source": [
    "s1_cube_red.download(work_dir+'/S1_after.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "597488f8-d09a-4aaa-9ec2-64d4543be524",
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = ilwis.RasterCoverage('S1_after.tif')\n",
    "print(s1.size())\n",
    "print(s1.envelope())\n",
    "coordSys = s1.coordinateSystem()\n",
    "coordSys.toWKT()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58c5fb51-c600-418c-a49e-db1928b87cff",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "VV = ilwis.do('selection',s1,'rasterbands(0)')\n",
    "VH = ilwis.do('selection',s1,'rasterbands(1)')\n",
    "VV.store('VV.mpr')\n",
    "VH.store('VH.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b18a1af0-08a1-41f2-bcda-7836099152ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate the cross polarized image (VV/VH)\n",
    "VC = ilwis.do('mapcalc','(@1/@2)', VV, VH)\n",
    "VC.store('VC.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fa6dbe6-96a7-4e97-81b7-b34a25e52c52",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create histogram using large number of bins, given the fact that the input image is a float64, with both positive and negative numbers\n",
    "hist = VC.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n",
    "minPerc1, maxPerc1 = hist.calcStretchRange(1)\n",
    "VCs = ilwis.do('linearstretch',VC, minPerc1, maxPerc1)\n",
    "VCs = ilwis.do('mapcalc','iff(@1==?,255,@1)', VCs) #modify pixels over water\n",
    "VCs = ilwis.do('setvaluerange', VCs, 0, 255, 1)\n",
    "VCs.store('VCs.mpr')\n",
    "\n",
    "VC_2np = np.fromiter(iter(VCs), np.ubyte, VCs.size().linearSize()) \n",
    "VC_2np = VC_2np.reshape((VCs.size().ysize, VCs.size().xsize))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f74ff83b-6679-4426-97f8-1529ad0f4d6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#stretch the VV and VH bands\n",
    "hist = VV.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n",
    "minPerc1, maxPerc1 = hist.calcStretchRange(1)\n",
    "VVs = ilwis.do('linearstretch',VV, minPerc1, maxPerc1)\n",
    "VVs = ilwis.do('mapcalc','iff(@1==?,0,@1)', VVs) #remove bad pixels over water\n",
    "VVs = ilwis.do('setvaluerange', VVs, 0, 255, 1)\n",
    "VVs.store('VVs.mpr')\n",
    "\n",
    "VV_2np = np.fromiter(iter(VVs), np.ubyte, VVs.size().linearSize()) \n",
    "VV_2np = VV_2np.reshape((VVs.size().ysize, VVs.size().xsize))\n",
    "\n",
    "hist = VH.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n",
    "minPerc1, maxPerc1 = hist.calcStretchRange(1)\n",
    "VHs = ilwis.do('linearstretch',VH, minPerc1, maxPerc1)\n",
    "VHs = ilwis.do('mapcalc','iff(@1==?,0,@1)', VHs) #remove bad pixels over water\n",
    "VHs = ilwis.do('setvaluerange', VHs, 0, 255, 1)\n",
    "VHs.store('VHs.mpr')\n",
    "\n",
    "VH_2np = np.fromiter(iter(VHs), np.ubyte, VHs.size().linearSize()) \n",
    "VH_2np = VH_2np.reshape((VHs.size().ysize, VHs.size().xsize))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f140ce80-8883-4d5a-8a5a-6bab1055ff3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create  color composite (in RGB)\n",
    "cc = np.dstack((VV_2np, VH_2np, VC_2np))\n",
    "\n",
    "#band interleaved\n",
    "S1ALL = np.array([VC_2np, VH_2np, VV_2np])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92d91713-67cd-4af7-ad49-c9732f5d53d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7, 7))\n",
    "plt.imshow(cc)\n",
    "plt.axis(\"off\")\n",
    "plt.title('S1 polarization composite');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44281f2c-b9d5-4954-8490-18013b521f19",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create empty ilwis raster \n",
    "S1_ilw = ilwis.RasterCoverage()\n",
    "defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1))\n",
    "S1_ilw.setDataDef(defNumr)\n",
    "S1_ilw.setSize(ilwis.Size(4121, 4060, 3)) \n",
    "S1_ilw.setGeoReference(s1.geoReference())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e0e583f-a619-4aa3-a797-00f02b6c9f72",
   "metadata": {},
   "outputs": [],
   "source": [
    "#read the numpy array\n",
    "S1_ilw.array2raster(S1ALL.flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd2e0a89-5f3e-4345-9652-40008f2a6d30",
   "metadata": {},
   "outputs": [],
   "source": [
    "#comment if you don't want to store the image in an ILWIS maplist format\n",
    "S1_ilw.store('S1_after_processed.mpl')\n",
    "\n",
    "#uncomment if you want to store the image in geotif format\n",
    "#S1_ilw.store(\"S1_after.tif\", \"GTiff\", \"gdal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66661470-e9ad-46aa-ab23-b7d838806123",
   "metadata": {},
   "source": [
    "### Retrieve Sentinel 2 level 2A image from after the flood event\n",
    "(see also: https://www.copernicus.eu/en/media/image-day-gallery/collapse-alau-dam-nigeria)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01ee7982-c0b8-4d57-835e-7489aeb224f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#connection.describe_collection(\"SENTINEL2_L2A\")\n",
    "connection.describe_collection(\"SENTINEL2_L1C\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12805550-e400-4e83-9c9d-a6c5655c8edd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#note the date - after the flood\n",
    "t = [\"2024-09-10\", \"2024-09-15\"]\n",
    "s2_cube = connection.load_collection(\"SENTINEL2_L1C\",\n",
    "    spatial_extent={'west': 13.00, 'east': 13.37, 'south': 11.63, 'north': 12.00,\"crs\": \"EPSG:4326\"},\n",
    "    temporal_extent= t,\n",
    "    bands=['B02', 'B03', 'B04', 'B08'],\n",
    "    max_cloud_cover=100,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae9e1a1e-d9bd-4682-ad57-d197bf4e4ec9",
   "metadata": {},
   "outputs": [],
   "source": [
    "s2_cube.download(work_dir+\"/s2_directafter.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bf75db0-4394-42db-a687-726610d2d8c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#read the S2 image in ilwispy\n",
    "s2 = ilwis.RasterCoverage ('s2_directafter.tif')\n",
    "print(s2.size())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27071974-aac6-42bd-813d-2364a820ef2b",
   "metadata": {},
   "source": [
    "From the image select and assign the spectral bands as seperate variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dfdc0fd-1070-43e7-99de-5bdb7101035c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#note index starts from 0\n",
    "Blue_temp = ilwis.do('selection',s2,\"rasterbands(0)\")\n",
    "Green_temp = ilwis.do('selection',s2,\"rasterbands(1)\")\n",
    "Red_temp = ilwis.do('selection',s2,\"rasterbands(2)\")\n",
    "IR_temp = ilwis.do('selection',s2,\"rasterbands(3)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c32985fa-bdfd-443c-ab98-588a94a10bb3",
   "metadata": {},
   "source": [
    "Check for corrupt data and re-assign to a value of '0'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de687bf4-9a6a-4cc6-8284-8049745702ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "Blue = ilwis.do('mapcalc','iff(@1==?, 0, @1)', Blue_temp)\n",
    "Green = ilwis.do('mapcalc','iff(@1==?, 0, @1)', Green_temp)\n",
    "Red = ilwis.do('mapcalc','iff(@1==?, 0, @1)', Red_temp)\n",
    "IR = ilwis.do('mapcalc','iff(@1==?, 0, @1)', IR_temp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f14f762d-6834-4d43-be3f-1406b1112110",
   "metadata": {},
   "source": [
    "Conduct some further image processing for each of the spectral channels, calculation the histogram (now using a large number of bins), conducting a linear stretch omitting the lower and upper tail of the histogram (1 % cut-off) and assigning the stretched data to a data range from 0 to 255, whole numbers only (precision is set to 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91765392-84be-48a5-92e4-6116c59652aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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\n",
    "# note that this can also be done in a loop, but here each band is handled seperately\n",
    "stat_Blue = Blue.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n",
    "minPerc1, maxPerc1 = stat_Blue.calcStretchRange(1)\n",
    "Blues = ilwis.do('linearstretch',Blue, minPerc1, maxPerc1)\n",
    "Blues = ilwis.do('setvaluerange', Blues, 0, 255, 1)\n",
    "\n",
    "stat_Green = Green.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n",
    "minPerc1, maxPerc1 = stat_Green.calcStretchRange(1)\n",
    "Greens = ilwis.do('linearstretch',Green, minPerc1, maxPerc1)\n",
    "Greens = ilwis.do('setvaluerange', Greens, 0, 255, 1)\n",
    "\n",
    "stat_Red = Red.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n",
    "minPerc1, maxPerc1 = stat_Red.calcStretchRange(1)\n",
    "Reds = ilwis.do('linearstretch',Red, minPerc1, maxPerc1)\n",
    "Reds = ilwis.do('setvaluerange', Reds, 0, 255, 1)\n",
    "\n",
    "stat_IR = IR.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n",
    "minPerc1, maxPerc1 = stat_IR.calcStretchRange(1)\n",
    "IRs = ilwis.do('linearstretch',IR, minPerc1, maxPerc1)\n",
    "IRs = ilwis.do('setvaluerange', IRs, 0, 255, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1144a8dd-a676-45c0-b910-ea8e8eb910a4",
   "metadata": {},
   "source": [
    "For visualization in Matplotlib (using Imshow) the data is transformed into a numpy array, using the ILWISPY operation 'iter', each pixel is assigned into a 1 D array, which by numpy is reshaped back into a 2-D array - for each of the spectral channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "578827d8-e562-4d7a-9262-5bbd9074642c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transform to a numpy array for visualization using imshow\n",
    "row = IRs.size().ysize \n",
    "col = IRs.size().xsize\n",
    "dim = IRs.size().linearSize()\n",
    "print(row, col, dim)\n",
    "\n",
    "Blues_2np = np.fromiter(iter(Blues), np.ubyte, dim) \n",
    "Blues_2np = Blues_2np.reshape((row, col))\n",
    "\n",
    "Greens_2np = np.fromiter(iter(Greens), np.ubyte, dim) \n",
    "Greens_2np = Greens_2np.reshape((row, col))\n",
    "\n",
    "Reds_2np = np.fromiter(iter(Reds), np.ubyte, dim) \n",
    "Reds_2np = Reds_2np.reshape((row, col))\n",
    "\n",
    "IRs_2np = np.fromiter(iter(IRs), np.ubyte, dim) \n",
    "IRs_2np = IRs_2np.reshape((row, col))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2656779d-4ff0-4b54-a2fc-584354fec42a",
   "metadata": {},
   "source": [
    "A numpy multi dimensional data stack is created, for visualization in Matplotlib the data is organized pixel interleaved, to use the results later as an ILWIS Maplist, the data stack is transformed into a multi-dimensional array which is band interleaved"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb2a7d61-4c64-4920-907b-3e3078aa319a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create  color composite (in RGB)\n",
    "#pixel interleaved\n",
    "fc_S2 = np.dstack((IRs_2np, Reds_2np, Greens_2np))\n",
    "nc_S2 = np.dstack((Reds_2np, Greens_2np, Blues_2np))\n",
    "\n",
    "#band interleaved\n",
    "S2ALL = np.array([Blues_2np, Greens_2np, Reds_2np, IRs_2np])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfd0eb9f-1099-4a10-a406-7f85aa7b05e8",
   "metadata": {},
   "source": [
    "Create a RGB plot showing a true color composite as well as a false colour composite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fe89489-8498-4b5b-a624-65a6ab88b0df",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1 = plt.figure(figsize=(15, 10))\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(nc_S2)\n",
    "plt.axis(\"off\")\n",
    "plt.title('NC - S2 image retrieved from OpenEO service provider')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(fc_S2)\n",
    "plt.axis(\"off\")\n",
    "plt.title('FC - S2 image retrieved from OpenEO service provider');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb81ccf0-97d1-4d3b-bf2a-016d20975a22",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create empty ilwis raster \n",
    "S2_ilw = ilwis.RasterCoverage()\n",
    "defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1))\n",
    "S2_ilw.setDataDef(defNumr)\n",
    "S2_ilw.setSize(ilwis.Size(4121, 4060, 4)) \n",
    "S2_ilw.setGeoReference(s2.geoReference())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4e7c4d2-5af1-48c5-bd3a-5db7e45cf6ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "#read the numpy array\n",
    "S2_ilw.array2raster(S2ALL.flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1333891a-4dc1-4422-8ef1-e16e62eb0698",
   "metadata": {},
   "outputs": [],
   "source": [
    "#comment if you don't want to store the image in an ILWIS maplist format\n",
    "S2_ilw.store('S2_directafter_processed.mpl')\n",
    "\n",
    "#uncomment if you want to store the image in geotif format\n",
    "#S2_ilw.store(\"S2_after_processed.tif\", \"GTiff\", \"gdal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6952bf42-a7d9-4b30-bca1-9091bf5a5ed5",
   "metadata": {},
   "source": [
    "### Retrieval of NDVI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ceb9cc0-cc94-4717-a84c-7f25e459a722",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate the ndvi\n",
    "Red = ilwis.do('selection',s2,\"rasterbands(2)\")\n",
    "IR = ilwis.do('selection',s2,\"rasterbands(3)\")\n",
    "\n",
    "ndvi = ilwis.do('mapcalc','(@1 - @2)/(@1 + @2)', IR, Red)\n",
    "ndvi = ilwis.do('setvaluerange', ndvi, -1, 1, 0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b793162-ee93-49fb-ba91-80f305855d89",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check size and value range\n",
    "print(ndvi.size())\n",
    "ndv_min = ndvi.min()\n",
    "ndv_max =ndvi.max()\n",
    "print(ndv_min, ndv_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d1b7aea-daa4-4321-86b2-f2fdfebf8f99",
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi.store('S2_after_ndvi.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c2632f8-4ec7-4fea-9fb7-9032507e6371",
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert to numpy array for visualization\n",
    "ndvi_2np = np.fromiter(iter(ndvi), np.float64, ndvi.size().linearSize()) \n",
    "ndvi_2np = ndvi_2np.reshape((ndvi.size().ysize, ndvi.size().xsize))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97f32902-7a67-416c-8864-b7918715d905",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Keep only negative NDVI values (representing water), mask the rest\n",
    "ndvi_water = np.ma.masked_where(ndvi_2np >= 0, ndvi_2np)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8b66c69-5e62-4d12-84e9-c68e15e2a94f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig2 = plt.figure(figsize=(15, 10))\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(ndvi_2np, interpolation ='nearest', vmin=0, vmax=ndv_max, cmap='turbo')\n",
    "plt.colorbar(shrink=0.5, label = 'NDVI')\n",
    "plt.title('NDVI derived from Sentinel-2 L1C, Maiduguri town and surrounding area', fontsize = 10)\n",
    "plt.axis(\"off\")\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(ndvi_water, cmap='Blues_r', vmin=ndv_min, vmax=0)\n",
    "plt.colorbar(label=\"NDVI (water range only)\", shrink=0.5)\n",
    "plt.title('Water areas (NDVI < 0) - note clouded / shaded areas!!', fontsize = 10)\n",
    "plt.axis(\"off\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9041c220-4664-488b-8373-ca8de6cab586",
   "metadata": {},
   "source": [
    "### Some further Hydro-DEM preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "864e5ee8-0eb8-4a3d-9ca9-ca6d3a82aa4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#once more import the Copernicus DEM in ilwispy\n",
    "copdem = ilwis.RasterCoverage ('Maiduguri_dem.tif')\n",
    "print(copdem.size())\n",
    "print(copdem.envelope())\n",
    "coordSys = copdem.coordinateSystem()\n",
    "coordSys.toWKT()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85d29c5e-c397-4079-afd6-c386114ccf08",
   "metadata": {},
   "outputs": [],
   "source": [
    "# target georeference (the grid you want to match)\n",
    "s2_grf = s2.geoReference()\n",
    "\n",
    "# choose resampling method: nearestneighbour | bilinear | bicubic\n",
    "copdem_res = ilwis.do(\"resample\", copdem, s2_grf, \"nearestneighbour\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e7596d6-cbac-4ca6-a08b-98e0467adc43",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sub map creation\n",
    "copdem_ressub = ilwis.do('selection',copdem_res,'boundingbox(50 50, 3999 3999)')\n",
    "print(copdem_ressub.size().xsize)\n",
    "print(copdem_ressub.size().ysize)\n",
    "print(copdem_ressub.size().zsize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f1d4354-3f35-4820-94b8-86883e4bf760",
   "metadata": {},
   "outputs": [],
   "source": [
    "copdem_ressub.store(\"copdem_ressub.mpr\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6a23f75-23d2-4c53-99d4-9ada5eb1cdec",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats = copdem_ressub.statistics(ilwis.PropertySets.pHISTOGRAM)\n",
    "#print(stats.histogram())\n",
    "print('minimum =',(stats[ilwis.PropertySets.pMIN])) # minimum value of the dem\n",
    "print('maximum =',(stats[ilwis.PropertySets.pMAX])) # maximum value of the dem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b522da4c-2dbf-4778-bd4e-5fa51463c40c",
   "metadata": {},
   "source": [
    "#### Fill Sinks (routine requires some time)\n",
    "\n",
    "MapFillSinks(inputraster, fill|cut)\n",
    "+ Inputraster is the input DEM map\n",
    "+ fill|cut is a method indicate to remove local depressions or cut terrain from a DEM raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36984928-2257-4696-9f28-350b79f4a758",
   "metadata": {},
   "outputs": [],
   "source": [
    "#uncomment line below if running the notebook for the first time\n",
    "fillsink = ilwis.do('MapFillSinks', copdem_ressub, 'cut')\n",
    "fillsink = ilwis.do('MapFillSinks', fillsink, 'cut')\n",
    "fillsink = ilwis.do('MapFillSinks', fillsink, 'cut')\n",
    "fillsink = ilwis.do('MapFillSinks', fillsink, 'fill')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ba9b2ad-39ce-42b4-bb8b-cea037534094",
   "metadata": {},
   "outputs": [],
   "source": [
    "#uncomment line below if running the notebook for the first time and the fill sink operation has been executed\n",
    "fillsink.store('copdem_fill.mpr')\n",
    "\n",
    "#uncomment if you want to store the image in geotif format\n",
    "#fillsink.store(\"fillsink.tif\", \"GTiff\", \"gdal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61574463-8962-4441-af7f-4a64f3f7cd9c",
   "metadata": {},
   "source": [
    "#### Flow Direction\n",
    "\n",
    "MapFlowDirection(inputraster, slope|height, yes|no)\n",
    "+ Inputraster is the input raster map (sink-free DEM)\n",
    "+ slope|height is a method to indicate how flow direction should be calculated (the steepst slope or the smallest height)\n",
    "+ yes|no is an option for the Parallel drainage correction algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37be5354-7094-41ea-814e-ead49a06775e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_fd = ilwis.do('MapFlowDirection', fillsink, 'slope', 'yes')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9d6fe1e-c722-446a-96cc-22a0e48b0176",
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_fd.store('copdem_fd.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99b9d17e-a8f3-4a50-9fa6-a296cb71ead0",
   "metadata": {},
   "source": [
    "#### Flow Accumulation\n",
    "\n",
    "MapFlowAccumulation(inputraster)\n",
    "+ Inputraster is the input flow direction map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bfda692-1b8d-4567-a376-4c393b85ffce",
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_fa = ilwis.do('MapFlowAccumulation', dem_fd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9e66b0b-7caf-48b6-ae86-63ecf65b210a",
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_fa.store('copdem_fa.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b58e521-982f-4c06-aaac-fcdd8e5d5a6f",
   "metadata": {},
   "source": [
    "#### Drainage Network extraction - define single threshold\n",
    "\n",
    "ExtractDrainageUseThresholdValue(FlowAccumulationRaster,thresholdval) \n",
    "+ FlowAccumulationRaster is the flow accumulation map\n",
    "+ Thresholdval is a threshold value (integer > 0) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7075f11-e12a-4547-99c9-04deb26dcb01",
   "metadata": {},
   "outputs": [],
   "source": [
    "#drainage raster map creation using single flow accumulation threshold - just downstream of dam breach location \n",
    "breach_inlet = ilwis.do('ExtractDrainageUseThresholdValue', dem_fa, 3051741) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce68e321-329e-47ea-9f63-730a0df663d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "breach_inlet.store('br_inlet.mpr')\n",
    "\n",
    "#uncomment if you want to store the image in geotif format\n",
    "#breach_inlet.store(\"br_inlet.tif\", \"GTiff\", \"gdal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2f2caa1-c526-4e5a-b433-f1f37c75d944",
   "metadata": {},
   "source": [
    "### FABDEM\n",
    "\n",
    "FABDEM (Forest And Buildings removed Copernicus DEM) is a global 30-meter resolution digital elevation model (DEM). By using machine learning to strip away the heights of trees and urban structures, it provides an accurate \"bare-earth\" terrain model (DTM) that is widely used for global flood and natural hazard modeling. Source: University of Bristol (https://research-information.bris.ac.uk/en/datasets/fabdem-2/). The downloaded DEM is available in your data folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bb846c2-577b-4655-9040-30c5b21d2ad1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#read the fabdem in ilwispy - note the extent from the file name (tile of 1 by 1 degree)\n",
    "fabdem = ilwis.RasterCoverage ('N11E013_FABDEM_V1-2.tif')\n",
    "print(fabdem.size())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6794d3f0-9d2f-4bd5-91c0-92a1c53d6eac",
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert to numpy array for visualization\n",
    "fabdem_2np = np.fromiter(iter(fabdem), np.float64, fabdem.size().linearSize()) \n",
    "fabdem_2np = fabdem_2np.reshape((fabdem.size().ysize, fabdem.size().xsize))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4152c65f-fa4d-497d-8925-4d9145113e89",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create plot - note area of interest is in upper left hand corner\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.imshow(fabdem_2np, cmap='terrain', vmin=300, vmax=350)\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.title('FABDEM provided by the University of Bristol');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c89240f1-37df-4570-82ed-b3fabcaa1e7d",
   "metadata": {},
   "source": [
    "Ensure that the elevation model is having the same projection, extent and resolution (as the Copernicus DEM and the satellite images)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f30a45a-cd8b-4502-b61d-0fb08b7c1536",
   "metadata": {},
   "outputs": [],
   "source": [
    "# target georeference (the grid you want to match)\n",
    "s2_grf = s2.geoReference()\n",
    "\n",
    "# choose resampling method: nearestneighbour | bilinear | bicubic\n",
    "fabdem_res = ilwis.do(\"resample\", fabdem, s2_grf, \"nearestneighbour\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c5bc836-fb49-4e6b-b582-f30d05f53e2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sub map creation\n",
    "fabdem_ressub = ilwis.do('selection',fabdem_res,'boundingbox(50 50, 3999 3999)')\n",
    "print(fabdem_ressub.size().xsize)\n",
    "print(fabdem_ressub.size().ysize)\n",
    "print(fabdem_ressub.size().zsize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07aab987-d885-46b0-a5eb-5a2529ab1425",
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert to numpy array for visualization\n",
    "fabdemressub_2np = np.fromiter(iter(fabdem_ressub), np.float64, fabdem_ressub.size().linearSize()) \n",
    "fabdemressub_2np = fabdemressub_2np.reshape((fabdem_ressub.size().ysize, fabdem_ressub.size().xsize))\n",
    "\n",
    "#create plot\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.imshow(fabdemressub_2np, cmap='terrain')\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.title('Resampled FABDEM');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "761baa79-41c6-49c1-a551-16be8d931e32",
   "metadata": {},
   "outputs": [],
   "source": [
    "fabdem_ressub.store(\"fabdem_ressub.mpr\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf8f8f91-8b3d-42b5-a86d-05554cf28cf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats1 = fabdem_ressub.statistics(ilwis.PropertySets.pHISTOGRAM)\n",
    "#print(stats.histogram())\n",
    "print('minimum =',(stats1[ilwis.PropertySets.pMIN])) # minimum value of the dem\n",
    "print('maximum =',(stats1[ilwis.PropertySets.pMAX])) # maximum value of the dem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8b299ba-0724-44ac-9ee7-4152d70e2bde",
   "metadata": {},
   "source": [
    "To see the impact of especially the urban artifact removal the 2 elevation model differences are derived"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afa4516f-a70a-47ea-bd77-6fc76ec69484",
   "metadata": {},
   "outputs": [],
   "source": [
    "difference = ilwis.do('mapcalc','@1-@2', fabdem_ressub, copdem_ressub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6aa88491-0247-4ea1-bfc0-057897888f31",
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert to numpy array for visualization\n",
    "difference_2np = np.fromiter(iter(difference), np.float64, difference.size().linearSize()) \n",
    "difference_2np = difference_2np.reshape((difference.size().ysize, difference.size().xsize))\n",
    "\n",
    "#create plot\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.imshow(difference_2np, cmap='jet', vmin=-3, vmax=3)\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.title('difference = FABDEM - COPDEM');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8d9f598-12a7-40ce-926e-bccbc98ee4f7",
   "metadata": {},
   "source": [
    "Given the impact of the urban area adjustments we are going to use the fabdem. We need to ensure hydrological consistency, so we do the cut / fill operation once more. First we 'adjust' the main channel, by enforcing the previous retrieved channel (breach_inlet) and lower the FABDEM with 1 meter for this drainage pathway. Note that bathymetry is not included in the DTM. Furthermore, as can be seen on the difference map above, along the main channel some pixels have obtained higher values!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f21f0798-2627-48f4-b353-3e1978e4936d",
   "metadata": {},
   "outputs": [],
   "source": [
    "drconditioned = ilwis.do('mapcalc','iff(@1==1,@2-1,@2)', breach_inlet, fabdem_ressub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ed71ac1-fb23-42dc-8f8d-8daaedc0473a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#uncomment line below if running the notebook for the first time\n",
    "#cut operation only done once since urban artefacts are already removed\n",
    "fillsink = ilwis.do('MapFillSinks', drconditioned, 'cut')\n",
    "fillsink = ilwis.do('MapFillSinks', fillsink, 'fill')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d0b44b0-c523-48b5-90a1-b3800a159878",
   "metadata": {},
   "outputs": [],
   "source": [
    "fillsink.store(\"fabdem_fill_dr_cond.mpr\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "313ef3c1-06ef-4e28-b7c9-f5a369cf59e0",
   "metadata": {},
   "source": [
    "### Some modeling using the data pre-processed"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e9c72f1-b546-4252-bbc1-4f40fd5c185f",
   "metadata": {},
   "source": [
    "### Model 1: A simplified flood water spreading model\n",
    "\n",
    "The model is a distance-limited flood spreading model starting from the river network. It assumes the river is already at a high water level everywhere, then uses a breadth-first search (BFS) to spread a gradually lowering water surface over the DEM while accounting for terrain elevation and floodplain friction. Flooding is restricted to areas within a maximum distance from the river, and water can only propagate where the simulated water surface remains above the ground. The result is an equilibrium flood extent and depth map that represents potential floodplain inundation rather than a time-dependent river flood event."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6bc27bb-0f7f-4544-bb97-bfe5fb51f140",
   "metadata": {},
   "source": [
    "Within this model the start location of the flood is upstream of where one expects the flood wave to spread into the floodplain. Once more have a look at the visualization of the elevation model and the drainage line created, based upon the flow accumulation threshold selected. Note that the upstream section of the river - below the dam - is confined to an incised narrow river reach and a flood wave of a few meters would not overtop the banks. This is different once river is passing throught the town, see also the Sentinel-2 image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a144592-17b7-4a79-b641-fb9d6ef9e5fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import dem and selected river reach based on flow accumulation threshold and create visualization\n",
    "\n",
    "# ---initialize size and dimensions---\n",
    "row = fillsink.size().ysize \n",
    "col = fillsink.size().xsize\n",
    "dim = fillsink.size().linearSize()\n",
    "\n",
    "# --- DEM read raster to numpy array ---\n",
    "dem = np.fromiter(iter(fillsink), np.float64, dim) \n",
    "dem = dem.reshape((row, col))\n",
    "rows, cols = dem.shape\n",
    "\n",
    "# --- process the river reach and create a raster ---\n",
    "river_raster = np.fromiter(iter(breach_inlet), np.float64, dim) \n",
    "river_raster = river_raster.reshape((row, col))\n",
    "river_mask = river_raster > 0\n",
    "\n",
    "print(\"River pixels:\", np.sum(river_mask))\n",
    "\n",
    "\n",
    "# --- dilate river to make it wider ---\n",
    "river_mask_wide = binary_dilation(river_mask, iterations=4)  # Approximate river width (pixels) = 1 + 2 * iterations\n",
    "\n",
    "# --- create river water surface for plotting ---\n",
    "water_surface_river = np.full_like(dem, np.nan)\n",
    "water_surface_river[river_mask_wide] = river_raster[river_mask_wide] \n",
    "\n",
    "# plot the figure\n",
    "plt.figure(figsize=(8,8))\n",
    "\n",
    "# 1. DEM\n",
    "plt.imshow(dem, cmap='terrain')\n",
    "plt.colorbar(label='DEM Elevation (m)', shrink=0.5)\n",
    "\n",
    "# 2. Wider river for visibility\n",
    "plt.imshow(np.ma.masked_invalid(water_surface_river), cmap='gray')\n",
    "\n",
    "plt.title(\"DEM with selected river section\")\n",
    "plt.axis('off')\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09f34ea0-3d53-4370-b6d1-649b4514752d",
   "metadata": {},
   "source": [
    "#### Adjust the following parameters to change the flood extent:\n",
    "+ water_depth_river to raise river water\n",
    "+ spread_decay to lower friction to spreads farther\n",
    "+ max_spread_distance to allow more floodplain area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46429ef2-8f05-4542-95f1-77c842258330",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input paramter settings\n",
    "pixel_size = 10               # DEM resolution (m)\n",
    "\n",
    "#Adjust accoring to your requirements\n",
    "water_depth_river = 1.50      # water above river (m) - initial value 1.5 \n",
    "spread_decay = 0.00085        # friction across floodplain per pixel - lower friction results in further spread -initial value 0.00085\n",
    "max_spread_distance = 2000    # max flood distance (m) from river - initial value 2000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc1a2f46-6138-4f2f-aa05-f74c96bdb0ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Preparation of the other required input\n",
    "\n",
    "# --- retrieve georeference and coordinate system info ---\n",
    "georef = fillsink.geoReference()\n",
    "\n",
    "crs = fillsink.coordinateSystem()\n",
    "print(crs.toWKT())\n",
    "\n",
    "# --- raster size ---\n",
    "ny = fillsink.size().ysize\n",
    "nx = fillsink.size().xsize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54857bc8-a5df-4c57-96e3-544677277b2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Initial water surface on river, set up of the starting water level inside the river network.\n",
    "#It creates the initial water surface that will later spread across the floodplain. The code creates a raster of water \n",
    "#surface elevation and fills it only at river locations using DEM elevation + assumed river water depth. The first line \n",
    "#creates a new raster having the same size and grid as the DEM and is filled with NaN (no value) everywhere. \n",
    "#line 2 assignes the water level to river pixels through means of boolean indexing, if the value is true (river pixel)\n",
    "#then water surface elevation = terrain elevation + water depth. Note the model assumes that the river everywhere contains\n",
    "#a uniform water depth above the channel bed!\n",
    "\n",
    "water_surface = np.full_like(dem, np.nan)\n",
    "water_surface[river_mask] = dem[river_mask] + water_depth_river"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a35dc79d-6e51-4803-a712-d6ba1ff9b0fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Distance to river (for the max. spread). This code field computes how far every pixel is from the river\n",
    "#and limits flooding to a maximum distance acting as a very important geomorphologic constraint in the model.\n",
    "#(~river_mask) = NOT river, edt uses the Euclidean Distance Transform: The EDT computes, for every non-river pixel the \n",
    "#distance to the nearest river pixel. It results in a map representing a kind of distance gradient away from the river. \n",
    "\n",
    "distance_pixels = distance_transform_edt(~river_mask)\n",
    "distance_m = distance_pixels * pixel_size # convert pixels to meters\n",
    "distance_m[distance_m > max_spread_distance] = np.nan #limit the flood spreading distance "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7da680b-009c-46f7-b23e-3b1e7c83bc36",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Breadth-First Search flood propagation from river. It spreads a water surface outward from the river across the \n",
    "#DEM until the water runs out of energy. This code section spreads the river water surface across the terrain using \n",
    "#a breadth-first search algorithm, reducing water level with distance and stopping where the ground becomes \n",
    "#too high or too far from the river.\n",
    "\n",
    "#This raster will store the final water surface elevation after spreading. At the start all values = NaN.\n",
    "water_surface_spread = np.full_like(dem, np.nan)\n",
    "\n",
    "#Initialize the spreading queue (flood sources)\n",
    "queue = deque()\n",
    "\n",
    "#np.argwhere(river_mask) returns the coordinates of all river pixels. Each river pixel will act as a \n",
    "#starting source of floodwater\n",
    "river_indices = np.argwhere(river_mask)\n",
    "\n",
    "#Seed the model with river water. Copy the river water surface into the spreading map and add all river \n",
    "#cells to the queue. The queue now contains the initial flood front. \n",
    "\n",
    "for r, c in river_indices:\n",
    "    water_surface_spread[r, c] = water_surface[r, c]\n",
    "    queue.append((r, c))\n",
    "\n",
    "#Define 8-neighbour connectivity to approximates 2-D overland flow\n",
    "neighbors = [(-1,0),(1,0),(0,-1),(0,1),(-1,-1),(-1,1),(1,-1),(1,1)]\n",
    "\n",
    "#Start BFS flood propagation. The queue works like a moving flood front, take one flooded cell and try\n",
    "#to spread water to its neighbours and continues until no more cells can be flooded. \n",
    "\n",
    "while queue:\n",
    "    r, c = queue.popleft()\n",
    "    current_ws = water_surface_spread[r, c]\n",
    "\n",
    "    for dr, dc in neighbors:\n",
    "        rr, cc = r + dr, c + dc #Check neighbour cells\n",
    "\n",
    "        if rr < 0 or rr >= rows or cc < 0 or cc >= cols:\n",
    "            continue\n",
    "\n",
    "        # skip if too far from river\n",
    "        if distance_m[rr, cc] > max_spread_distance: #Limit flood distance from river\n",
    "            continue\n",
    "\n",
    "        #apply friction decay, representing the hydraulic energy loss on the floodplain\n",
    "        #each step away from the river the water surface drops slightly (to mimics friction\n",
    "        #and shallow flow resistance). So the flood surface gradually slopes downward.\n",
    "       \n",
    "        new_ws = current_ws - spread_decay * pixel_size #Apply floodplain friction loss\n",
    "\n",
    "        #Check terrain blocking, stop if water too shallow. Water can only flood a cell if\n",
    "        #water surface>ground elevation. This lets the DEM control ridges to block water, depressions fill\n",
    "        #and valleys to guide the flow This is the key terrain constraint.\n",
    "        \n",
    "        if new_ws <= dem[rr, cc]:\n",
    "            continue\n",
    "\n",
    "        #Update if higher water level. It ensures the model always keeps the highest possible water surface \n",
    "        #at each cell.This allows multiple flow paths, gradual correction and equilibrium solution, making\n",
    "        #the algorithm diffusive and stable.\n",
    "        \n",
    "        if np.isnan(water_surface_spread[rr, cc]) or new_ws > water_surface_spread[rr, cc]:\n",
    "            \n",
    "            #continue spreading, newly flooded cell becomes a new source. The flood keeps expanding until\n",
    "            #it runs out of energy.The loop ends /  the queue becomes empty when water cannot climb terrain anymore\n",
    "            #or friction removed too much energy or max river distance is reached. At that moment the flood reaches \n",
    "            #equilibrium obtaining the final water surface elevation everywhere the flood can reach. \n",
    "            water_surface_spread[rr, cc] = new_ws\n",
    "            queue.append((rr, cc))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "129af076-f729-49d1-b3a5-96f05e171a36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create flood mask and depth\n",
    "flood_mask = water_surface_spread > dem\n",
    "flood_depth = np.where(flood_mask, water_surface_spread - dem, 0)\n",
    "\n",
    "print(\"Flooded pixels:\", np.sum(flood_mask))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b4c2c3b-e53c-4b3b-9e6e-5641cdce197c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualization\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.title(\"Flood Mask\")\n",
    "plt.imshow(flood_mask, cmap=\"Blues\")\n",
    "plt.colorbar(shrink=0.5)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.title(\"Flood Depth\")\n",
    "plt.imshow(flood_depth/2, cmap=\"Blues\", vmin=0, vmax=2.5)\n",
    "plt.colorbar(shrink=0.5)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c66740dc-8406-46a9-80d6-1134dea6e3c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# --- dilate river to make it wider ---\n",
    "river_mask_wide = binary_dilation(river_mask, iterations=4)  # Approximate river width (pixels) = 1 + 2 * iterations\n",
    "\n",
    "# --- create river water surface for plotting ---\n",
    "water_surface_river = np.full_like(dem, np.nan)\n",
    "water_surface_river[river_mask_wide] = river_raster[river_mask_wide]# + water_depth_river\n",
    "\n",
    "# --- plot ---\n",
    "plt.figure(figsize=(8,8))\n",
    "\n",
    "# 1. DEM\n",
    "plt.imshow(dem, cmap='terrain')\n",
    "plt.colorbar(label='DEM Elevation (m)', shrink=0.5)\n",
    "\n",
    "# 2. Floodplain water surface\n",
    "plt.imshow(np.ma.masked_invalid(water_surface_spread), cmap='Blues', alpha=0.75)\n",
    "\n",
    "# 3. Wider river for visibility\n",
    "plt.imshow(np.ma.masked_invalid(water_surface_river), cmap='gray')\n",
    "\n",
    "plt.title(\"DEM with water surface mask and selected river section\")\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66fd10bd-6c53-402a-8977-5bb04270d77b",
   "metadata": {},
   "source": [
    "#### Store the results obained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9df859e6-00be-4b04-b43a-ccc0ebf9a2a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create empty ilwis raster for flood depth map\n",
    "fl_depth = ilwis.RasterCoverage()\n",
    "defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 100, 0.001))\n",
    "fl_depth.setDataDef(defNumr)\n",
    "fl_depth.setSize(ilwis.Size(fabdem_ressub.size().ysize, fabdem_ressub.size().xsize, 1)) \n",
    "fl_depth.setGeoReference(fabdem_ressub.geoReference())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14748963-d175-4542-8f5e-a4187e3459ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "fl_depth.array2raster(flood_depth.flatten())\n",
    "\n",
    "#adjust depth - as depth is exaggerated - to be further eleaborated upon as no temporal flood wave dynamics\n",
    "#(hydrograph) is currently implemented\n",
    "fl_depth_adj = ilwis.do('mapcalc','iff(@1>0,@1/2,?)', fl_depth)\n",
    "\n",
    "#comment if you don't want to store the image in an ILWIS maplist format\n",
    "fl_depth_adj.store('M1_flood_depth.mpr')\n",
    "\n",
    "#uncomment if you want to store the image in geotif format\n",
    "#fl_depth_adj.store(\"flood_depth.tif\", \"GTiff\", \"gdal\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3de175e-5835-4bcd-acac-ce3cfd793eaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create empty ilwis raster for flood mask map\n",
    "fl_mask = ilwis.RasterCoverage()\n",
    "defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 1,1))\n",
    "fl_mask.setDataDef(defNumr)\n",
    "fl_mask.setSize(ilwis.Size(fabdem_ressub.size().ysize, fabdem_ressub.size().xsize, 1)) \n",
    "fl_mask.setGeoReference(fabdem_ressub.geoReference())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d31213c1-3639-4192-a1bd-6fa73e544423",
   "metadata": {},
   "outputs": [],
   "source": [
    "fl_mask.array2raster(flood_mask.flatten())\n",
    "\n",
    "#comment if you don't want to store the image in an ILWIS maplist format\n",
    "fl_mask.store('M1_flood_mask.mpr')\n",
    "\n",
    "#uncomment if you want to store the image in geotif format\n",
    "#fl_mask.store(\"flood_mask.tif\", \"GTiff\", \"gdal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca38e3eb-c457-4d9b-8a5f-58a7d788eb49",
   "metadata": {},
   "source": [
    "### Model 2: Drainage controlled flow model - a raster hydraulic routing & floodplain propagation model\n",
    "\n",
    "Below a hybrid 1D river routing + 2D floodplain spread DEM flood model using frictional water surface decay instead of full hydrodynamics is prepared. This approach features a fast DEM-based hydraulic flood model with river routing, bank overtopping and floodplain spread. It also uses energy loss instead of solving shallow-water equations, which makes it fast and stable.\n",
    "\n",
    "By running the code fields below a 3-stage raster flood model is executed:\n",
    "+ Route water down the river network\n",
    "+ Check where rivers overtop their banks\n",
    "+ Spread water across the floodplain"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd261002-c6ee-4b06-bae5-2513ce832a2f",
   "metadata": {},
   "source": [
    "#### Model input parameters (what drives the simulation)\n",
    "The model uses:\n",
    "+ Topography (represented by the DEM filled raster and it's spatial resolution)\n",
    "+ River network (a binary river mask raster, extracted from the flow accumulation threshold applied)\n",
    "\n",
    "Hydraulic parameters are: \n",
    "+ (Multiple) inflow point(s) (upstream water level source, row / column location of the source(s), representing upstream boundary water level)\n",
    "+ Initial flood depth at inflow (water depth above ground at inflow cell, typical ranges for river flood: 0.5–3 m, for dam breach: 3–20 m) \n",
    "+ Longitudinal energy loss (downstream friction, unit = m/m, representing the longitudinal slope / hydraulic gradient, typical values for flat floodplain:\t0.00002 – 0.0001, for large rivers:\t0.00005 – 0.0003)\n",
    "+ Lateral decay (floodplain friction, same concept as above, unit is m/m, it must be larger than longitudinal loss because floodplains are rougher. Typical ratio: lateral decay is about 3 to 10 times the longitudinal loss\n",
    "+ Bankfull depth (when rivers spills, unit = meter. Typical values for bankfull depth are for small stream:\t0.5–1 m, medium river: 1–3 m and large river: 3–8 m)\n",
    "\n",
    "\n",
    "| Parameter           | Units        | Physical meaning          |\n",
    "| ------------------- | ------------ | ------------------------- |\n",
    "| Inflow locations    | raster index | upstream boundary         |\n",
    "| Initial flood depth | **m**        | starting water depth      |\n",
    "| Longitudinal loss   | **m/m**      | downstream energy slope   |\n",
    "| Lateral decay       | **m/m**      | floodplain friction slope |\n",
    "| Bankfull depth      | **m**        | channel capacity          |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ea60311-deb4-4e08-b41e-97371f6f8037",
   "metadata": {},
   "source": [
    "##### Initial values given to the hydraulic parameters\n",
    "\n",
    "In the code field below modify the parameters according to your requirements:\n",
    "+ flood_depth_inlet initial values: low = 0.45, medium = 1.25 and high = 2.25\n",
    "+ longitudinal_loss: initial value = 0.000285\n",
    "+ lateral_decay: initial value = 0.00085\n",
    "+ bankfull_depth: initial value = 1.85"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4efbae5c-582b-49c3-a14c-c52293bdc5f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# (Hydraulic) input parameter settings\n",
    "\n",
    "# --- raster pixel size (in meters) of your resampled DEM ---\n",
    "cellsize = 10\n",
    "\n",
    "inflows = [(3005, 2976)]  #note index offset as python starts from 0\n",
    "                         \n",
    "\n",
    "# --- other hydraulic parameter settings --- Adjust upon your requirement\n",
    "flood_depth_inlet = 2.25\n",
    "longitudinal_loss = 0.000285 #m/m\n",
    "lateral_decay = 0.00085 #m/m\n",
    "bankfull_depth = 1.85  # bankfull channel depth in m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5070c4a-cd96-4f69-8213-4c7286816fc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Preparation of required input raster data and transformation of sink-filled DEM \n",
    "#and river raster (drainage line starting at inlet) to numpy arrays\n",
    "#load the original datasets once more\n",
    "\n",
    "# --- DEM check size and dimensions---\n",
    "row = fillsink.size().ysize \n",
    "col = fillsink.size().xsize\n",
    "dim = fillsink.size().linearSize()\n",
    "print(row, col, dim)\n",
    "\n",
    "# --- DEM read raster to numpy array ---\n",
    "dem = np.fromiter(iter(fillsink), np.float64, dim) \n",
    "dem = dem.reshape((row, col))\n",
    "\n",
    "# --- retrieve georeference and coordinate system info ---\n",
    "georef = fillsink.geoReference()\n",
    "\n",
    "crs = fillsink.coordinateSystem()\n",
    "print(crs.toWKT())\n",
    "\n",
    "# --- raster size ---\n",
    "ny = fillsink.size().ysize\n",
    "nx = fillsink.size().xsize\n",
    "\n",
    "# --- RIVER MASK read ilwis raster to numpy array ---\n",
    "river_raster = np.fromiter(iter(breach_inlet), np.float64, dim) \n",
    "river_raster = river_raster.reshape((row, col))\n",
    "\n",
    "river = np.array(river_raster) == 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b65e6290-f19a-4e61-8306-6458b6914be5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating the model state grids — the maps that will store the simulation results while the flood propagates \n",
    "# At the start of the model for all pixels: river_stage = NaN, \tflood_depth = 0 and flood_mask = False\n",
    "\n",
    "river_stage = np.full_like(dem, np.nan, dtype=np.float32) # River routing fills river_stage\n",
    "flood_depth = np.zeros_like(dem, dtype=np.float32) # Bank spill starts filling flood_depth\n",
    "flood_mask = np.zeros_like(dem, dtype=bool) #Floodplain spreading fills flood_mask"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "509a9db2-61f2-4b5d-8a56-7698ffb00e2f",
   "metadata": {},
   "source": [
    "A 3-step approach to create a flood depth / flood mask map\n",
    "\n",
    "| Step                   | Physics                   |\n",
    "| ---------------------- | ------------------------- |\n",
    "| River routing          | 1-D channel hydraulics    |\n",
    "| Bankfull spill         | river–floodplain coupling |\n",
    "| Floodplain propagation | 2-D overland flow         |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3817fb3-f5b9-490f-80f7-2bcf7bbb7b61",
   "metadata": {},
   "outputs": [],
   "source": [
    "# River routing, computes the water surface profile along the river network before any flooding happens.In hydraulic terms it represents\n",
    "# a simplified steady 1-D river routing step. Spread water downstream through river cells while gradually losing energy.\n",
    "# The algorithm used is Breadth-First Search (BFS) on the river network.\n",
    "\n",
    "#Create routing memory\n",
    "visited = np.zeros_like(dem, dtype=bool)\n",
    "queue = deque()\n",
    "\n",
    "# Insert inflow boundary conditions\n",
    "# At each inflow point: Compute starting water surface elevation (stage=ground elevation+water depth)\n",
    "# Store this in the river map. Put the cell into the routing queue.\n",
    "# This defines the upstream boundary condition. Now the queue contains the first water cells.\n",
    "\n",
    "for r,c in inflows:\n",
    "    stage = dem[r,c] + flood_depth_inlet\n",
    "    river_stage[r,c] = stage\n",
    "    queue.append((r,c,stage))\n",
    "    visited[r,c] = True\n",
    "\n",
    "# While there are river cells with water to process: take the next cell from the queue. Try to send water to neighbouring river cells, \n",
    "# water moves cell by cell through the network.\n",
    "while queue:\n",
    "    r,c,stage = queue.popleft()\n",
    "\n",
    "#Explore neighbouring cells (D8 connectivity), skip the centre cell:\n",
    "    for dr in [-1,0,1]:\n",
    "        for dc in [-1,0,1]:\n",
    "            if dr==0 and dc==0:\n",
    "                continue\n",
    "\n",
    "# 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.\n",
    "# So water stays inside the river network.\n",
    "            nr,nc = r+dr, c+dc\n",
    "            if not (0<=nr<ny and 0<=nc<nx):\n",
    "                continue\n",
    "            if not river[nr,nc]:\n",
    "                continue\n",
    "            if visited[nr,nc]:\n",
    "                continue\n",
    "\n",
    "# Compute hydraulic energy loss. This is the physics of the model. Water loses elevation as it moves downstream: \n",
    "# new water level=old level−(slope×distance), Where: distance = cell size (or diagonal length), longitudinal_loss = hydraulic gradient (m/m)\n",
    "# This mimics: friction, bed slope, energy dissipation. So the river water surface tilts downstream.\n",
    "            \n",
    "            dist = np.hypot(dr,dc)*cellsize\n",
    "            new_stage = stage - longitudinal_loss*dist\n",
    "\n",
    "# Prevent uphill river flow. Water cannot flow if the water surface would be below the ground. This prevents impossible uphill routing.\n",
    "            \n",
    "            min_stage = dem[nr,nc]\n",
    "            if new_stage < min_stage:\n",
    "                continue\n",
    "\n",
    "# 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\n",
    "# The flood wave keeps moving.\n",
    "            \n",
    "            river_stage[nr,nc] = new_stage\n",
    "            visited[nr,nc] = True\n",
    "            queue.append((nr,nc,new_stage))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bae9ce34-8550-4e4b-bbf8-4f46b2b6bd82",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bankfull spill - bankfull exceedance / overbank flow. Computation of the downstream river water surface profile by routing water through river cells \n",
    "# while applying a hydraulic slope (energy loss). Hydraulically, this is the channel and floodplain coupling step.\n",
    "# Up to now the model only knows the water level inside the river. This block checks: Where does the river overflow its banks?\n",
    "# This code section finds where the river water level exceeds the bank height and converts those river cells into sources of floodplain flooding.\n",
    "# This step simulates: channel capacity, levee overtopping, bank overflow, start of overbank flooding. \n",
    "\n",
    "# Create a new queue for flooding. This queue will store starting points of floodplain flooding. \n",
    "# River routing queue: when water moving in channel. Spill queue: when water moving onto land\n",
    "\n",
    "spill_queue = deque()\n",
    "\n",
    "# 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.\n",
    "\n",
    "for r,c in zip(*np.where(river)):\n",
    "\n",
    "# Skip river cells without water. If the routing step never reached this river cell, it has no water, it cannot spill, thenignore it.\n",
    "    if np.isnan(river_stage[r,c]):\n",
    "        continue\n",
    "        \n",
    "# Compute bankfull elevation, his defines the top of the river banks (bankfull elevation=ground level+channel depth)\n",
    "# Example: ground = 12 m, bankfull_depth = 1.5 m, banks at 13.5 m. This is the maximum water level the channel can hold.\n",
    "    \n",
    "    bankfull_elev = dem[r,c] + bankfull_depth\n",
    "\n",
    "# Check if the river overtops, compare: river water surface and bank height. If water level is below banks nothing happens.\n",
    "# If water level is above banks then flooding starts.\n",
    "\n",
    "    spill_depth = river_stage[r,c] - bankfull_elev\n",
    "\n",
    "    if spill_depth <= 0:\n",
    "        continue\n",
    "\n",
    "# Mark the river cell as flooded, Once overtopping occurs, the river cell becomes part of the flooded terrain.\n",
    "# Flood depth = water surface − ground elevation. So the river becomes the first flooded area.\n",
    "    \n",
    "    flood_mask[r,c] = True\n",
    "    flood_depth[r,c] = river_stage[r,c] - dem[r,c]\n",
    "\n",
    "# Add spill locations to flood queue. These cells are now the starting points of floodplain spreading.\n",
    "# These act like “holes in the river banks” where water pours onto land.\n",
    "    \n",
    "    spill_queue.append((r,c,river_stage[r,c]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b185282f-00c7-479e-94a2-f58c48aa47ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lateral floodplain propagation which is spreading floodwater across the terrain from overtopped river cells, losing energy with distance \n",
    "# and stopping where the ground becomes too high.This is the actual floodplain flooding step — where water spreads over land after leaving the river.\n",
    "# Hydraulically this represents 2-D overland flood spreading with friction losses. If Step 1 = river flow and Step 2 = river overflow\n",
    "# then in the code field below = water is moving across the terrain. Spread water from spill locations across the DEM while losing energy\n",
    "# and respecting terrain elevations represents basically a 2-D hydraulic fill model using the Breadth-First Search (BFS) algorithm.\n",
    "# When the queue becomes empty it has created the full flood extent, flood depth and floodplain water surface.\n",
    "# The flood stops when water loses too much energy to climb terrain. This step mimics sheet flow, overbank spreading, ponding in depressions\n",
    "# and floodplain hydraulics, it is similar to a simplified 2-D shallow water model, but much faster.\n",
    "\n",
    "\n",
    "# Create floodplain visitation map. Just like in river routing, it tracks which floodplain cells have already been processed.\n",
    "# This prevents: infinite loops, repeated calculations and water bouncing back and forth\n",
    "\n",
    "visited_fp = np.zeros_like(dem, dtype=bool)\n",
    "\n",
    "# Start flood spreading from spill cells. The previous step added overtopped river cells to this queue. Now each of those cells becomes\n",
    "# a source of overland flow and Water spreads outward.\n",
    "\n",
    "while spill_queue:\n",
    "    r,c,stage = spill_queue.popleft()\n",
    "\n",
    "# Explore neighbouring terrain cells (D8). Water can spread in all 8 directions. This is a simplified 2-D hydraulic connectivity.\n",
    "    for dr in [-1,0,1]:\n",
    "        for dc in [-1,0,1]:\n",
    "            if dr==0 and dc==0:\n",
    "                continue\n",
    "\n",
    "#Skip invalid neighbours, skip cells that are outside the raster and already received flood water\n",
    "            \n",
    "            nr,nc = r+dr, c+dc\n",
    "            if not (0<=nr<ny and 0<=nc<nx):\n",
    "                continue\n",
    "            if visited_fp[nr,nc]:\n",
    "                continue\n",
    "                \n",
    "# Compute floodplain energy loss. Water loses elevation due to vegetation, buildings, rough terrain, shallow flow friction, etc.\n",
    "# so the flood surface gradually slopes downward away from the river.\n",
    "            \n",
    "            dist = np.hypot(dr,dc)*cellsize\n",
    "            new_stage = stage - lateral_decay*dist\n",
    "\n",
    "# Check terrain blocking. Water can only flood a cell if the water surface is above the ground. This automatically creates natural levees, \n",
    "# terraces, ridges, flow barriers, etc. In this way the DEM controls where water can go.\n",
    "            \n",
    "            if new_stage <= dem[nr,nc]:\n",
    "                continue\n",
    "\n",
    "# Mark the cell as flooded. Compute flood depth on land: depth = water_surface − ground. This creates the final flood map. \n",
    "            \n",
    "            flood_mask[nr,nc] = True\n",
    "            flood_depth[nr,nc] = new_stage - dem[nr,nc]\n",
    "\n",
    "# Continue flood propagation. The flooded cell becomes a new source of flooding. This is how the flood wave expands outward.\n",
    "            \n",
    "            visited_fp[nr,nc] = True\n",
    "            spill_queue.append((nr,nc,new_stage))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be5a62dd-10ea-4477-8e4a-f8032c035520",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualization\n",
    "print(\"Flooded pixels:\", np.sum(flood_mask))\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.title(\"Flood Mask\")\n",
    "plt.imshow(flood_mask, cmap=\"Blues\")\n",
    "plt.colorbar(shrink=0.5)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.title(\"Flood Depth\")\n",
    "plt.imshow(flood_depth/2, cmap=\"Blues\", vmin=0, vmax=2.5)\n",
    "plt.colorbar(shrink=0.5)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d57ef3e5-f84f-4d81-a1b8-f81aa80d8d66",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dilate river to make it wider (3 pixels)\n",
    "river_mask_wide = binary_dilation(river, iterations=4)  # Approximate river width (pixels) = 1 + 2 * iterations\n",
    "\n",
    "# Create river water surface for plotting\n",
    "water_surface_river = np.full_like(dem, np.nan)\n",
    "water_surface_river[river_mask_wide] = river[river_mask_wide]\n",
    "\n",
    "# Plot\n",
    "plt.figure(figsize=(10,8))\n",
    "\n",
    "# 1. DEM\n",
    "plt.imshow(dem, cmap='terrain')\n",
    "plt.colorbar(label='DEM Elevation (m)', shrink=0.5)\n",
    "\n",
    "# 2. Flood mask\n",
    "mask_masked = np.ma.masked_where(flood_mask <= 0, flood_mask)\n",
    "plt.imshow(mask_masked, cmap='Blues', alpha=0.5)\n",
    "\n",
    "\n",
    "# 3. Wider river for visibility\n",
    "plt.imshow(np.ma.masked_invalid(water_surface_river), cmap='gray')\n",
    "\n",
    "plt.title(\"DEM with flood mask and selected river section\")\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6af0d674-031c-474f-a10e-1f977dbe0907",
   "metadata": {},
   "source": [
    "#### Store the results obtained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c32cc2f-cfea-4024-9a7b-b5c0909626c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create empty ilwis raster for flood depth map\n",
    "fl_depth = ilwis.RasterCoverage()\n",
    "defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 100, 0.001))\n",
    "fl_depth.setDataDef(defNumr)\n",
    "fl_depth.setSize(ilwis.Size(fillsink.size().ysize, fillsink.size().xsize, 1)) \n",
    "fl_depth.setGeoReference(fillsink.geoReference())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2eba0c8e-f32f-4a2f-b0c8-e4eea24c86f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fl_depth.array2raster(flood_depth.flatten())\n",
    "\n",
    "#adjust depth - as depth is exaggerated - to be further eleaborated upon as no temporal flood wave dynamics (hydrograph) is implemented\n",
    "fl_depth_adj = ilwis.do('mapcalc','iff(@1>0,@1/4,?)', fl_depth)\n",
    "\n",
    "#comment if you don't want to store the image in an ILWIS format\n",
    "fl_depth_adj.store('M2_flood_depth.mpr')\n",
    "\n",
    "#uncomment if you want to store the image in geotif format\n",
    "#fl_depth_adj.store(\"flood_depth.tif\", \"GTiff\", \"gdal\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fe9f436-c40d-416f-aecf-b4db5c3738ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create empty ilwis raster for flood mask map\n",
    "fl_mask = ilwis.RasterCoverage()\n",
    "defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 1,1))\n",
    "fl_mask.setDataDef(defNumr)\n",
    "fl_mask.setSize(ilwis.Size(fillsink.size().ysize, fillsink.size().xsize, 1)) \n",
    "fl_mask.setGeoReference(fillsink.geoReference())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea9125d7-b090-4ddb-8585-10abc9238a97",
   "metadata": {},
   "outputs": [],
   "source": [
    "fl_mask.array2raster(flood_mask.flatten())\n",
    "\n",
    "#comment if you don't want to store the image in an ILWIS format\n",
    "fl_mask.store('M2_flood_mask.mpr')\n",
    "\n",
    "#uncomment if you want to store the image in geotif format\n",
    "#fl_mask.store(\"flood_mask.tif\", \"GTiff\", \"gdal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d176ac09-55de-4a11-be95-7a8850f7e97c",
   "metadata": {},
   "source": [
    "#### Plot flood depth map on top of Sentinel-2 FCC image and NDVI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5269f783-53a7-463d-9c21-0282e7b39221",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Mask dry pixels\n",
    "depth = flood_depth.copy()/4\n",
    "depth[depth <= 0.05] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0367d1e5-732c-41ae-b8c1-d567a530a7fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose depth color scaling\n",
    "max_depth = (np.nanmax(depth))  # adjust to your flood\n",
    "norm_depth = np.clip((depth) / max_depth, 0, 1)\n",
    "print(max_depth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db210c99-e516-43ff-b041-1ff70df8bf7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to RGBA overlay,  create a blue colormap manually. deep water = darker + more opaque, \n",
    "# shallow water = lighter + transparent. Dry pixels fully transparent.\n",
    "\n",
    "overlay = np.zeros((*depth.shape, 4))\n",
    "\n",
    "# Blue color\n",
    "overlay[..., 2] = 1.0\n",
    "\n",
    "# Transparency based on depth\n",
    "overlay[..., 3] = np.nan_to_num(norm_depth) * 0.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a31fcbd-61e7-430c-8f87-5bed5fc767f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualization\n",
    "plt.figure(figsize=(12,6))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(fc_S2)\n",
    "\n",
    "im = plt.imshow(depth, cmap=\"Blues\", alpha=0.45, vmin=0, vmax=max_depth)\n",
    "plt.title(\"Flood depth mask over Sentinel-2 FCC\", fontsize = 10)\n",
    "plt.axis(\"off\")\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(ndvi_water, cmap='Blues_r', vmin=ndv_min, vmax=0)\n",
    "im = plt.imshow(depth, cmap=\"Blues\", alpha=0.65, vmin=0, vmax=max_depth)\n",
    "plt.title('Water areas (NDVI < 0) - note clouded / shaded areas and flood depth mask', fontsize = 10)\n",
    "plt.axis(\"off\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fca4226-3ca5-454e-bbaa-5dee681cb116",
   "metadata": {},
   "source": [
    "Note the main model differences:\n",
    "| Model 1                                     | Model 2                                    |\n",
    "| ------------------------------------------- | ------------------------------------------ |\n",
    "| River = flood source everywhere             |  River → spill → floodplain                |\n",
    "| Assumes river already flooded               |  Simulates bank overtopping                |\n",
    "| No river routing                            | Has river routing (longitudinal hydraulics)|\n",
    "| Single-stage spreading                      |Two-stage hydraulics                        |\n",
    "| More GIS / cost-distance style              |   Physically closer to hydraulics          |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1226929f-a3a7-4beb-a6c2-b8c5dd7a0335",
   "metadata": {},
   "source": [
    "### Concluding remarks:\n",
    "To conduct a detailed assessment of a dam breach using a hydraulic model a lot of data is required which is often not available. The approach presented here could be used as an alternative to conduct a first assessment of the area that might be affected by a flood resulting from a dam breach. The data required is limited to an elevation model and a number of assumptions on bankfull depth, initial flood depth and decay parameters. The parameters can be tuned and results can be compared with a satellite image for validation. Model 1 is a slightly more simplified approach compared to model 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "662bc6cf-f8cd-4121-bafb-ca80139aabc9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d7c2558-5941-4083-8850-93dd329dae4a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}