{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5f05c7a0",
   "metadata": {},
   "source": [
    "## Elevation data processing for Hydrology\n",
    "\n",
    "Notebook prepared by Ben Maathuis, Lichun Wang and Bas Retsios. ITC-University of Twente, Enschede. The Netherlands\n",
    "\n",
    "The routines provided below are an implementation in ilwispy of the 'DEM Hydro-processing' routines available in the desktop version of ILWIS. For details on the routines also consult the article: 'Digital Elevation Model based Hydro-processing', available at: https://www.researchgate.net/publication/228661426_Digital_Elevation_Model_Based_Hydro-processing."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b051d7fe-0489-4a96-b21b-9669ef43dbad",
   "metadata": {},
   "source": [
    "Sample data is available at: https://filetransfer.itc.nl/pub/52n/ilwis_py/sample_data_V2/DEM_processing.zip. Unzip the file, note the content, a tif file containing the elevation data and two outlet location point maps in ilwis format. Note that a UTM coordinate system is used.  Here it is assumed that the data folder '/DEM_processing' is situated within the notebook folder! It is furthermore assumed that you have locally installed ILWIS386 for data visualization when data is written to disk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3827fd5e-9726-4f84-a5c0-3abf715c2311",
   "metadata": {},
   "outputs": [],
   "source": [
    "# if relevant update ilwispy by uncommenting the line below, use at least ilwispy version of 20260708\n",
    "#!pip install ilwis --upgrade"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1aee7b0",
   "metadata": {},
   "source": [
    "#### Load libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fa11ae3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import ilwis\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import LineString"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7467544d-6f0e-423a-84fc-c3a7faa8c230",
   "metadata": {},
   "outputs": [],
   "source": [
    "#suppress warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51497714-49e3-495b-ba3c-41e2a181fb66",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check your ilwispy version used\n",
    "ilwis.version()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b3a8264",
   "metadata": {},
   "outputs": [],
   "source": [
    "work_dir = os.getcwd()+'/DEM_processing'\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": "57900538",
   "metadata": {},
   "outputs": [],
   "source": [
    "cd_folder = os.getcwd() #check current directory folder\n",
    "print(cd_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f024ee35",
   "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": "17bee2a3-cde2-4186-8e3c-c38d4a59bc65",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import elevation model\n",
    "dem = ilwis.RasterCoverage('Thiba_DEM_UTM.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4faf2f43-dd1e-4ae3-9c84-7b51403ffdc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#retrieve meta data from imported map\n",
    "print(dem.size())\n",
    "coordSys = dem.coordinateSystem()\n",
    "print(coordSys.toWKT())\n",
    "print()\n",
    "print(dem.envelope())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63ff33ad-842f-4350-98bc-8d91ed23d22e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#retrieve some basic statistics from imported map\n",
    "stats_dem = dem.statistics(ilwis.PropertySets.pHISTOGRAM)\n",
    "#print(stats_dem.histogram())\n",
    "dem_min = stats_dem[ilwis.PropertySets.pMIN] # minimum value on the map\n",
    "dem_max = stats_dem[ilwis.PropertySets.pMAX] # maximum value on the map\n",
    "print(dem_min)\n",
    "print(dem_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbeefced-c441-40a6-ad18-4c8b20d0f89e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#transform imported map from ilwis format into a numpy array using the (ilwispy) pixel itereator\n",
    "dem_2np = np.fromiter(iter(dem), np.float64, dem.size().linearSize()) \n",
    "#now we overwrite the initial variable created\n",
    "dem_2np = dem_2np.reshape((dem.size().ysize, dem.size().xsize))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c9c0bad-7361-47f9-bf16-9bbd58982729",
   "metadata": {},
   "outputs": [],
   "source": [
    "#display the numpy array created in the previous step using matplotlib \n",
    "fig = plt.figure(figsize =(10, 7))\n",
    "plt.imshow(dem_2np, interpolation='none', vmin=dem_min, vmax=dem_max, cmap='jet')\n",
    "\n",
    "plt.axis('on')\n",
    "plt.colorbar(shrink=0.45)\n",
    "plt.title('DEM - south of Mount Kenya');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fbcc32d-f9ae-4edd-9084-34630fdaf122",
   "metadata": {},
   "source": [
    "#### Fill Sinks\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": "aac00855-6bee-4212-864a-34e7b12e0c5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#first we apply the cut routine followed by the fill\n",
    "cut = ilwis.do('MapFillSinks', dem, 'cut')\n",
    "fillsink = ilwis.do('MapFillSinks', cut, 'fill')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23f2c127-34e0-4d3a-a2ce-d8385f98f2e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#store your results obtained, check the maps produced using the ilwis desktop version\n",
    "fillsink.store('dem_fill.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2e5ca77-52ba-436c-a66c-75e88d5d13b5",
   "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": "3615220b-65cf-4171-a1e1-47f6b94be2cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_fd = ilwis.do('MapFlowDirection', fillsink, 'slope', 'yes')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc9de143-4fe2-4b3e-8e23-309f8c81f764",
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_fd.store('dem_fd.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcf97b09-13b3-4d7f-8b08-11fc2347b497",
   "metadata": {},
   "source": [
    "### Flow Accumulation\n",
    "\n",
    "MapFlowAccumulation(inputraster)\n",
    "+ Inputraster is the input flow direction map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1248adc2-ee59-494b-a9b2-19e6ab529bff",
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_fa = ilwis.do('MapFlowAccumulation', dem_fd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e55791f-2245-41f0-ad91-114f1c7dc44f",
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_fa.store('dem_fa.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55342089-91fa-40e6-98c4-5b7cc7bd5afb",
   "metadata": {},
   "source": [
    "### Internal relief map calculation\n",
    "\n",
    "This command is used to generate an internal relief raster based on the specified filter size, can be used to evaluate the threshold boundaries and number of classes when computing the 'variable threshold map' (see below)\n",
    "+ internalrelief(inputraster,filtersize)\n",
    "+ Inputraster is the input DEM map\n",
    "+ Filtersize is an odd value to define the filter size\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84fe4822-f667-4fbc-bfb6-808dafa12731",
   "metadata": {},
   "outputs": [],
   "source": [
    "int_relief = ilwis.do('InternalRelief', fillsink, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3e4ddd7-e943-489d-99cb-81f5850a55e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "int_relief.store ('int_relief.mpr') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3411aaca-0bdb-404f-a03b-d36f8c70429c",
   "metadata": {},
   "source": [
    "### Variable threshold computation \n",
    "\n",
    "+ variablethreshold(inputraster,filtersize,NrOfClasses,UpperBounds_and_ThresholdVals)\n",
    "+ inputraster is the input DEM map\n",
    "+ filtersize is a value for the size of the window\n",
    "+ NrOfClasses is a value for the number of classes\n",
    "+ UpperBounds_and_ThresholdVals is a string in which, for each class at a time, the upper boundary values and threshold values are specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fea0857-a83a-40eb-882c-c413ab3bf1a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create threshold raster, 5 is the filtersize, 4 are the number of threshold classes used\n",
    "thres4 = ilwis.do('VariableThreshold', dem, 5, 4, '10, 15000, 25, 10000, 60, 5000, 1000, 1500')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddcb4dc3-7b94-4351-8a80-36ca675be2f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "thres4.store('thres4.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27f85a81-a88d-4009-9dad-599744935f86",
   "metadata": {},
   "source": [
    "### Drainage Network extraction(1) - multiple thresholds\n",
    "\n",
    "ExtractDrainageUseThresholdRaster(inputraster1,inputraster2,inputraster3) \n",
    "+ This command is used to extract the drainage according to the input threshold map\n",
    "+ Inputraster1is the flow accumulation map\n",
    "+ Inputraster2 is the threshold map\n",
    "+ Inputraster3 is the flow direction map "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87ad4eeb-da3a-4325-ad08-6167b670fb84",
   "metadata": {},
   "outputs": [],
   "source": [
    "#drainage raster extraction using raster threshold map (using multiple flow accumulation thresholds)\n",
    "dr_thres4 = ilwis.do('ExtractDrainageUseThresholdRaster', dem_fa, thres4, dem_fd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ab7b6e1-77d3-448e-9de2-50ddad8f3da8",
   "metadata": {},
   "outputs": [],
   "source": [
    "dr_thres4.store('drain_thres4.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a1750d7-bdfa-410d-becf-0678a65947fe",
   "metadata": {},
   "source": [
    "### Drainage Network extraction(2) - 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": "83eecc92-0358-4d10-a23e-2ff7c131ec72",
   "metadata": {},
   "outputs": [],
   "source": [
    "#drainage raster map creation using single flow accumulation threshold\n",
    "dr_single7500 = ilwis.do('ExtractDrainageUseThresholdValue', dem_fa, 7500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1804a85c-44ef-48cb-9974-aa1190351d0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "dr_single7500.store('dr_7500.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce270615-a0c1-4d2f-8550-e4936d289d9d",
   "metadata": {},
   "source": [
    "### Drainage Network Ordering\n",
    "\n",
    "DrainageNetworkOrdering(DEMmap,FlowDiractionmap,DrainageNetworkmap,MinimumDrainageLength)\n",
    "+ DEMmap is the input DEM map\n",
    "+ FlowDirectionmap is the input Flow Direction map\n",
    "+ DrainageNetworkmap is the input drainage network map\n",
    "+ MinimumDrainageLength is a value for the minimum length (in m)\n",
    "\n",
    "This operation will generate the output drainage network raster and segment map. The output attribute table will be stored automatically. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d521edf-eadc-4028-a9e6-f44084c6b1b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "outTable, outSegment, outDNO = ilwis.do('DrainageNetworkOrdering',fillsink, dem_fd, dr_thres4, 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33e2d660-9b9a-49fc-ba9f-1387dfb2aff6",
   "metadata": {},
   "outputs": [],
   "source": [
    "outDNO.store('dr_net.mpr')\n",
    "outSegment.store('dr_net.mps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75b4fc63-3dd4-4b6a-a0e6-5659f97d9943",
   "metadata": {},
   "outputs": [],
   "source": [
    "#store segments as shape file\n",
    "outSegment.store('dr_net1.shp', 'ESRI Shapefile', 'gdal') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a959fe7-020e-4b75-96ce-4b9b70c8c6fa",
   "metadata": {},
   "source": [
    "### Catchment extraction\n",
    "\n",
    "MapCatchmentExtraction(DrainageNetworkOrderingMap,FlowDirectionMap)\n",
    "+ DrainageNetworkOrderingMap is the input drainage network ordering map\n",
    "+ FlowDirectionMap is the input flow direction map\n",
    "\n",
    "It produces the output catchment raster and polygon map. The attribute table will be stored automaticcaly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33990c58-647f-4bbf-a853-9614e773ccee",
   "metadata": {},
   "outputs": [],
   "source": [
    "catchm = ilwis.do('MapCatchmentExtraction', outDNO, dem_fd) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "122615d6-bf11-4fe3-b12e-971e4033b749",
   "metadata": {},
   "outputs": [],
   "source": [
    "#note the layers for polygons and raster\n",
    "catchm[0].store('catchment.mpa') \n",
    "catchm[1].store('catchment.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8800b96-3ac5-4e0e-8e49-491cb23c18c5",
   "metadata": {},
   "source": [
    "### Catchment merge - using outlet location\n",
    "This operation is used to merge the catchment according the outlet location. The operation produces an output raster map for the new merged catchments and related attribute table; an output segment map that contains only segments that fall within the new catchments and related attribute; an output segment map and related attribute table, that will contain the longest possible flow path in each new catchment.\n",
    "\n",
    "MapCatchmentMergeWithOutlet (DrainageNetworkOrderMap,FlowDiractionMap,FlowaccumulationMap,DEM,OutletPointMap,IncludeUndefinedPixels)\n",
    "\n",
    "+ DrainageNetworkOrderMap is the input drainage network oedering raster map \n",
    "+ FlowDirectionMap is the input flow direction map\n",
    "+ DEM is the input DEM map\n",
    "+ OutletPointMap is the input point map that contains the outlet location \n",
    "+ IncludeUndefinedPixels – yes|no, specify whether include the undefined pixels in the DEM into a catchment\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c3d64af-19d7-430f-bf02-dc33419e94b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First import the outlet point map\n",
    "outlet = ilwis.FeatureCoverage('outlet_location.mpp') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96527412-528a-4a46-be9a-2bd042531c71",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check imported point map\n",
    "for i in range(outlet.attributeCount()):\n",
    "    print(str(outlet[i]))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2415dd4-f2a6-463b-98a0-48d9586c0b30",
   "metadata": {},
   "outputs": [],
   "source": [
    "CatchMerge_outlet = ilwis.do('MapCatchmentMergeWithOutlet','dr_net.mpr', dem_fd,dem_fa,fillsink,outlet,'no')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f715457-e671-4b4b-9631-74da3e13e977",
   "metadata": {},
   "outputs": [],
   "source": [
    "#note different polygon and vector maps created\n",
    "CatchMerge_outlet[0].store('catchmerge_polygon.mpa') \n",
    "CatchMerge_outlet[1].store('catchmerge_raster.mpr') \n",
    "CatchMerge_outlet[2].store('catchmerge_long_flowpath.mps')\n",
    "CatchMerge_outlet[3].store('catchmerge_segments.mps')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58374e1d-aa81-442d-98e6-82179d47cb44",
   "metadata": {},
   "source": [
    "### Catchment merge - using stream order method\n",
    "\n",
    "This operation is used to merge the catchment according to the stream order method. The operation produces an output raster map for the new merged catchments and related attribute table. \n",
    "\n",
    "MapCatchmentMergeWithStreamOrder(DrainageNetworkOrderMap,FlowDirectionMap,FlowaccumulationMap,DEM,StreamOrderSystem,StreamOrderValue,ExtractOriginalOrder)\n",
    "\n",
    "+ DrainageNetworkOrderMap is the input drainage network ordering map\n",
    "+ FlowDirectionMap is the input flow direction map\n",
    "+ FlowaccumulationMap is the flow accumulation map\n",
    "+ DEM is the input DEM map\n",
    "+ StreamOrderSystem = strahler | shreve, specify whether you wish to use the Strahler or the Shreve ordering system.\n",
    "+ StreamOrderValue is the stream order value (integer value > 0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c05452c-5848-4732-b4b8-c64ef7c179b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "CatchMerge_order3 = ilwis.do('MapCatchmentMergeWithStreamOrder', 'dr_net.mpr', dem_fd, dem_fa, fillsink, 'strahler', 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a6ce9ec-0c75-46aa-87cc-d4fd9d21b91f",
   "metadata": {},
   "outputs": [],
   "source": [
    "CatchMerge_order3[0].store('catchmerge_strahler3_polygon.mpa') \n",
    "CatchMerge_order3[1].store('catchmerge_strahler3_raster.mpr') \n",
    "CatchMerge_order3[2].store('catchmerge_strahler3_long_flowpath.mps')\n",
    "CatchMerge_order3[3].store('catchmerge_strahler3_segment.mps') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ee69803-6410-4424-9f4d-f9bb66dcf0e4",
   "metadata": {},
   "source": [
    "### Overland Flow Length\n",
    "\n",
    "The Overland flow length operation calculates for each pixel the overland distance towards the 'nearest' drainage according to the flow paths available in a Flow Direction map.\n",
    "\n",
    "MapOverlandFlowLength(DrainageNetworkOrderMap,FlowDiractionMap)\n",
    "\n",
    "+ DrainageNetworkOrderMap is the input drainage network ordering map\n",
    "+ FlowDirectionMap is the input flow direction map\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70910a11-3c4b-4505-8e49-856a3abbe40c",
   "metadata": {},
   "outputs": [],
   "source": [
    "OverLandFlow = ilwis.do('MapOverlandFlowLength', outDNO, dem_fd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12942391-9f2c-4b27-9748-757d93f502f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "OverLandFlow.store('overland_flow.mpr') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "896f7c9c-d8c8-4b76-a240-6a3eff43dad8",
   "metadata": {},
   "source": [
    "### -\tFlow Length to Outlet\n",
    "\n",
    "MapFlowLength2Outlet(DrainageNetworkOrderMap,FlowDiractionMap)\n",
    "+ DrainageNetworkOrderMap is the input drainage network ordering map\n",
    "+ FlowDiractionMap is the input flow direction map\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4240f843-56de-43d1-ad07-f24b8efb47d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "length2outlet = ilwis.do('MapFlowLength2Outlet', outDNO, dem_fd) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4225cfb1-b230-4c7c-b9ac-807a7f282d1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "length2outlet.store('length2outlet.mpr') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ca306b6-cff8-491a-a327-9719d263e02c",
   "metadata": {},
   "source": [
    "### Some other ilwis hydro dem processing routines examples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e68db600-bf91-4d5f-9477-677fb40e6e1a",
   "metadata": {},
   "source": [
    "#### Calculate color shaded refief map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d57cb46-d7cf-43ed-8512-7c5cd2be3087",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create shadow filters: 'code=, dimensions = 3 by 3, followed by 9 filter values and eventually added by a gain factor, here 1'\n",
    "shadow_w = 'code=3,3,-2 -1 2 -3 1 4 -2 -1 2,1'\n",
    "shadow_n = 'code=3,3,-2 -3 -2 -1 1 -1 2 4 2,1'\n",
    "shadow_nw = 'code=3,3,-3 -2 -1 -2 1 2 -1 2 4,1' # or use 'shadownw'\n",
    "#shadow_nw = 'shadownw' #when applying a standard filter available in the ilwis filter library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45fca506-75f1-4afa-aa71-0688ecb0c2a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#shadowW\n",
    "shadow_west = ilwis.do('linearrasterfilter', dem, shadow_w)\n",
    "shadow_west = ilwis.do('linearstretch',shadow_west, 5)\n",
    "shadow_west = ilwis.do('setvaluerange', shadow_west, 0, 255, 1)\n",
    "shadow_west = ilwis.do('mapcalc','iff(@1==?,0,@1)',shadow_west)\n",
    "shadow_west_2np = np.fromiter(iter(shadow_west), np.ubyte, shadow_west.size().linearSize()) \n",
    "shadow_west_2np = shadow_west_2np.reshape((shadow_west.size().ysize, shadow_west.size().xsize))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52181e26-76c8-4070-a7a9-589c2ce223a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#shadowN\n",
    "shadow_north = ilwis.do('linearrasterfilter', dem, shadow_n)\n",
    "shadow_north = ilwis.do('linearstretch',shadow_north, 5)\n",
    "shadow_north = ilwis.do('setvaluerange', shadow_north, 0, 255, 1)\n",
    "shadow_north = ilwis.do('mapcalc','iff(@1==?,0,@1)',shadow_north)\n",
    "shadow_north_2np = np.fromiter(iter(shadow_north), np.ubyte, shadow_north.size().linearSize()) \n",
    "shadow_north_2np = shadow_north_2np.reshape((shadow_north.size().ysize, shadow_north.size().xsize))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4369a2e0-8f48-4e5a-b34e-46e0b0be3d28",
   "metadata": {},
   "outputs": [],
   "source": [
    "#shadowNW\n",
    "shadow_northwest = ilwis.do('linearrasterfilter', dem, shadow_nw)\n",
    "shadow_northwest = ilwis.do('linearstretch',shadow_northwest, 5)\n",
    "shadow_northwest = ilwis.do('setvaluerange', shadow_northwest, 0, 255, 1)\n",
    "shadow_northwest = ilwis.do('mapcalc','iff(@1==?,0,@1)',shadow_northwest)\n",
    "shadow_northwest_2np = np.fromiter(iter(shadow_northwest), np.ubyte, shadow_northwest.size().linearSize()) \n",
    "shadow_northwest_2np = shadow_northwest_2np.reshape((shadow_northwest.size().ysize, shadow_northwest.size().xsize))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c3ccf80-8b8b-432d-a9c6-83c06c6d7752",
   "metadata": {},
   "outputs": [],
   "source": [
    "nc_array = np.dstack((shadow_north_2np, shadow_northwest_2np, shadow_west_2np))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb9757ae-b7db-4388-8301-85d8e2f216c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#load shape file and note extent\n",
    "gdf = gpd.read_file(work_dir + \"/dr_net1.shp\")\n",
    "print(\"Shapefile bounds:\" ,gdf.total_bounds)\n",
    "print()\n",
    "#check coordinate system\n",
    "print(gdf.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8dde54c-4dad-4356-b583-7aca66571a24",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot DEM, NW-gray scale shaded relief and color shaded relief maps with overlay of drainage vector map\n",
    "\n",
    "# common extent\n",
    "extent = [299657.5, 336437.5, 9911539.0, 9983419.0]\n",
    "\n",
    "# Create 3 subplots\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 7))\n",
    "\n",
    "#DEM background\n",
    "im1 = ax1.imshow(dem_2np, extent=extent, cmap='terrain', vmin=1500, vmax=4500)\n",
    "gdf.plot(ax=ax1,facecolor=\"none\", edgecolor=\"red\",linewidth=0.5)\n",
    "\n",
    "ax1.set_title(\"DEM\")\n",
    "ax1.set_xlabel(\"Easting (m)\")\n",
    "ax1.set_ylabel(\"Northing (m)\")\n",
    "\n",
    "cbar1 = fig.colorbar(im1, ax=ax1,shrink=0.8)\n",
    "cbar1.set_label(\"Elevation (m)\")\n",
    "\n",
    "# Shaded Relief\n",
    "im2 = ax2.imshow(shadow_northwest_2np, extent=extent, cmap='gray')\n",
    "gdf.plot(ax=ax2, color='cyan', linewidth=0.5)\n",
    "\n",
    "ax2.set_title(\"Shaded Relief\")\n",
    "ax2.set_xlabel(\"Easting (m)\")\n",
    "ax2.set_ylabel(\"Northing (m)\")\n",
    "\n",
    "# Color shaded RGB Composite\n",
    "im3 = ax3.imshow(nc_array, extent=extent,zorder=1)\n",
    "gdf.plot(ax=ax3, facecolor=\"none\", edgecolor=\"yellow\", linewidth=0.5, zorder=2)\n",
    "\n",
    "ax3.set_title(\"RGB Composite\")\n",
    "ax3.set_xlabel(\"Easting (m)\")\n",
    "ax3.set_ylabel(\"Northing (m)\")\n",
    "\n",
    "fig.subplots_adjust(\n",
    "    left=0.10,\n",
    "    right=0.90,\n",
    "    bottom=0.08,\n",
    "    top=0.92,\n",
    "    wspace=0.01\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ecfb537-478f-49ae-a5de-d9342ca052d5",
   "metadata": {},
   "source": [
    "#### Conduct the processing, also using a looping procedure to stretch the initial filter results and store the final maps as an ilwis maplist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b461ba15-95c9-4d0b-bb54-50b8c1ec9fe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#once more compute the shadow maps\n",
    "shadow_west = ilwis.do('linearrasterfilter', dem, shadow_w)\n",
    "shadow_north = ilwis.do('linearrasterfilter', dem, shadow_n)\n",
    "shadow_northwest = ilwis.do('linearrasterfilter', dem, shadow_nw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a87005b-1a84-46cc-bcb8-fe45460a0b75",
   "metadata": {},
   "outputs": [],
   "source": [
    "rcNew = ilwis.RasterCoverage()\n",
    "rcNew.setGeoReference(ilwis.GeoReference(shadow_west.coordinateSystem(), shadow_west.envelope(), shadow_west.size()))\n",
    "rcNew.setSize(ilwis.Size(rcNew.size().xsize, rcNew.size().ysize, 3))\n",
    "rcNew.setDataDef(ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1)))\n",
    "\n",
    "rcNew.addBand(0, shadow_west.begin())\n",
    "rcNew.addBand(1, shadow_northwest.begin())\n",
    "rcNew.addBand(2, shadow_north.begin())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fce58295-2c6a-4b9a-88f5-3619ad04adff",
   "metadata": {},
   "outputs": [],
   "source": [
    "#stretch the rasterbands contained in rcNew\n",
    "multiple_stretch = []\n",
    "multiple_bands = ilwis.do('selection',rcNew,\"rasterbands(0..2)\") # for specific bands use band number as \"rasterbands(0,1,2,6)\"\n",
    "ls = ilwis.do('linearstretch',multiple_bands, 5)\n",
    "mb_stretch = ilwis.do('setvaluerange', ls, 0, 255, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c08443aa-1a24-431f-b337-c6740a0e6789",
   "metadata": {},
   "outputs": [],
   "source": [
    "#write the maplist to disk\n",
    "mb_stretch.store('dem_col.mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eba3543e-b725-464e-bf9e-779b1d26bf2a",
   "metadata": {},
   "source": [
    "##### Compare your ilwis and matplotlib results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bc7284a-2a7e-4658-a7cb-2c297a186b0a",
   "metadata": {},
   "source": [
    "### Create contour lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b1846ee-b732-4a8e-92b2-81b8dd6bcbc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check the numpy array size and elevation range in DEM\n",
    "print(dem_2np.shape)\n",
    "print(np.nanmin(dem_2np))\n",
    "print(np.nanmax(dem_2np))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "185d945d-af98-4207-9611-37684ada8d2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Remove the 'ilwis' no-data value(s) represented as -1e+308 to get the actual DEM data range\n",
    "dem_adj = np.where(dem_2np == -1e308, np.nan, dem_2np)\n",
    "print(np.nanmin(dem_adj))\n",
    "print(np.nanmax(dem_adj))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a7c01d8-16ed-42b4-94a6-69da05145ea2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Use the extent as defined before\n",
    "\n",
    "x = np.linspace(extent[0], extent[1], dem_adj.shape[1])\n",
    "y = np.linspace(extent[3], extent[2], dem_adj.shape[0])\n",
    "\n",
    "plt.figure(figsize=(10, 7))\n",
    "\n",
    "# plot adjusted DEM\n",
    "plt.imshow(dem_adj, extent=extent, cmap='terrain')\n",
    "\n",
    "# create contours, interval = 250 m\n",
    "levels = np.arange(np.floor(np.nanmin(dem_adj)/250)*250, np.ceil(np.nanmax(dem_adj)/250)*250 + 250, 250)\n",
    "\n",
    "# plot contours\n",
    "cs = plt.contour(x, y, dem_adj, levels=levels, colors='black',linewidths=0.5)\n",
    "\n",
    "plt.clabel(cs, inline=True, fontsize=8, fmt='%d')\n",
    "plt.xlabel(\"Easting (m)\")\n",
    "plt.ylabel(\"Northing (m)\")\n",
    "plt.title(\"DEM with 250 m contourline interval\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65ed41a5-f37e-4576-8da9-6646add0bbba",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Extract the contour lines\n",
    "lines = []\n",
    "elevations = []\n",
    "\n",
    "for level, seglist in zip(cs.levels, cs.allsegs):\n",
    "    for seg in seglist:\n",
    "        if len(seg) > 1:\n",
    "            lines.append(LineString(seg))\n",
    "            elevations.append(level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccab52a2-29d5-46fd-93b9-4b3c919a61b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create a GeoDataFrame\n",
    "gdf_contours = gpd.GeoDataFrame(\n",
    "    {\"elevation\": elevations},\n",
    "    geometry=lines,\n",
    "    crs=\"EPSG:32737\"   # UTM Zone 37S\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "862c3166-b5a2-40fd-a2d8-e8cdaae240bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check geodataframe created\n",
    "print(gdf_contours.head())\n",
    "print(len(gdf_contours))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd3b831f-f29f-4750-9153-2e4cd2263d5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Save as shapefile\n",
    "gdf_contours.to_file(work_dir+\"/mount_kenya_contours.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "529081d0-5b8d-42ea-a036-53ef832aae0e",
   "metadata": {},
   "source": [
    "#### Derive class-coverage statistics\n",
    "\n",
    "We continue working a bit with tables derived from crossing of raster maps. First load the raster coverages that are going to be used. Note that the 'catchmerge_raster' map also has an associated 'domain' and attribute table "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c4c84d4-1d31-4095-a891-856a76794e06",
   "metadata": {},
   "outputs": [],
   "source": [
    "#selected merged catchment map upstream of outlet location\n",
    "catch_sel = ilwis.RasterCoverage('catchmerge_raster.mpr')\n",
    "\n",
    "#selected cross map - here the threshold map with 4 classes is used\n",
    "raster_sel = ilwis.RasterCoverage('thres4.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60280812-56d8-4d95-b5b6-fa1a89d70665",
   "metadata": {},
   "outputs": [],
   "source": [
    "#perform the crossing operation - here the variable threshold map is used\n",
    "ct = ilwis.do('cross',catch_sel, raster_sel, 'ignoreundef')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba5dfa56-9ab3-49db-add3-eecd3b86213e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ct.store('ct.tbt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62e508c1-2a5d-459b-8605-0932499d872d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#review table size and columns names\n",
    "print(ct.columnCount())\n",
    "print(ct.recordCount())\n",
    "print(ct.columns())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7e64922-d856-4bdc-9950-6d6171d42093",
   "metadata": {},
   "outputs": [],
   "source": [
    "#review content of table\n",
    "for i in range(ct.recordCount()):\n",
    "    rec = ct.record(i)\n",
    "    print(rec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6070503c-bad8-4c3a-8a12-7b4cb5205278",
   "metadata": {},
   "source": [
    "ILWIS also has a table calculator, called: tabcalc. Below you find an example how the data in columns of a table is used to conduct calculations and a new output column and output table is created\n",
    "\n",
    "An explanation expression below: \n",
    "+ 'ct1' = output table\n",
    "+ '(@1/@2)' = calculation expression using 2 parameters, here (a/b), parameters to be defined lateron in the expression\n",
    "+ 'test_calc' = new output column'\n",
    "+ 'ct' = input table\n",
    "+ 'yes' or 'no' = if 'yes' option to create new table', if 'no' - existing table is used\n",
    "+ column pixel_count = @1\n",
    "+ column second_raster = @2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c03e31f-580c-4d2a-9653-db67cfbfbe1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#table calculator in ilwis\n",
    "ct1 = ilwis.do('tabcalc', '(@1/@2)', ct, 'test_calc','yes', 'pixel_count', 'second_raster')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33f9efd3-5168-43db-9aab-849edb9fdaee",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ct1.columns())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "627f930c-454e-4c84-a1fa-6ce2f2c39fd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(ct1.recordCount()):\n",
    "    rec = ct1.record(i)\n",
    "    print(rec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "565be62b-08dc-4365-848c-b36a681fb389",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ilwis syntax to get map statistics:|ilwis.PropertySets.pMAX|ilwis.PropertySets.pDISTANCE|ilwis.PropertySets.pDELTA|ilwis.PropertySets.pNETTOCOUNT|ilwis.PropertySets.pCOUNT|ilwis.PropertySets.pSUM|ilwis.PropertySets.pMEAN|ilwis.PropertySets.pMEDIAN|ilwis.PropertySets.pSTDEV|ilwis.PropertySets.pHISTOGRAM\n",
    "#to get percentage per class we need the total number of pixels!\n",
    "ct_stat = ilwis.PropertySets.pSUM\n",
    "pc_stat = ct.statistics('pixel_count',ct_stat) # derive the sum for the column 'pixel_count'\n",
    "pixel_count_sum = pc_stat.prop(ilwis.PropertySets.pSUM)\n",
    "print('Sum = ', pixel_count_sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "094ff368-6a13-4209-896d-0b2a104a7b33",
   "metadata": {},
   "outputs": [],
   "source": [
    "ct2 = ilwis.do('tabcalc', '(@1/@2*100)', ct, 'area_percent','yes', 'pixel_count', pixel_count_sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a695e2c9-d067-4a5e-91af-ef3337feac92",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ct2.columnCount())\n",
    "print(ct2.recordCount())\n",
    "print(ct2.columns())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "200d1374-bde8-4eea-9061-4d09047f7123",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(ct2.recordCount()):\n",
    "    rec = ct2.record(i)\n",
    "    print(rec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b568550e-514c-477d-a3bc-1a899e00dcd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ct2.store('ct_fin.tbt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "804d6193-477f-43bf-966f-ca5e7e7e50df",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ilwis cross table to pandas dataframe\n",
    "data = {}\n",
    "for col in list(ct2.columns()):\n",
    "    data[col] = list(ct2.column(col))\n",
    "df = pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "222f36e4-63a8-4eed-8de4-2ea1e2688f6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dd3c96a-c4a6-4342-a83c-f3b1947c59a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df['pixel_count'].sum())\n",
    "print()\n",
    "print(df['area_percent'].sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ac06fda-19db-4869-b7fc-62ae00b91624",
   "metadata": {},
   "outputs": [],
   "source": [
    "#store dataframe as csv file - without index column\n",
    "df.to_csv(work_dir+'/ct_fin.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d4f8feb-323e-4483-bee3-f82c74885953",
   "metadata": {},
   "source": [
    "#### Derive aggregate statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f468164-c998-4e42-a16e-83d16f602a41",
   "metadata": {},
   "outputs": [],
   "source": [
    "#perform the crossing operation - here the variable dem_fill map with elevation values is used (with precision of 1 - meter)\n",
    "raster_sel1 = ilwis.RasterCoverage('thres4.mpr')\n",
    "raster_sel2 = ilwis.RasterCoverage('dem_fill.mpr')\n",
    "raster_sel2 = ilwis.do('setvaluerange', raster_sel2, 0, 10000, 1)#set dem to interval of 1 meter - no decimals\n",
    "ct = ilwis.do('cross',raster_sel1, raster_sel2, 'ignoreundef')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52fb3165-00f9-46b6-94e3-d96436084e23",
   "metadata": {},
   "outputs": [],
   "source": [
    "print()\n",
    "print(ct.columnCount(), ct.recordCount())\n",
    "print(ct.columns())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58845fe3-98aa-4825-9e27-84acecdf6a13",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#check the first 10 records of the cross-table\n",
    "for i in range(min(10,ct.recordCount())):\n",
    "    rec = ct.record(i)\n",
    "    print(rec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18460904-fe70-4f79-9990-4d27ea4df893",
   "metadata": {},
   "outputs": [],
   "source": [
    "ct.store('aggr.tbt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d3b1e2a-7dcf-421d-ab06-321eddbcdc7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ilwis cross table to pandas dataframe\n",
    "data = {}\n",
    "for col in list(ct.columns()):\n",
    "    data[col] = list(ct.column(col))\n",
    "df = pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e98dc53-db59-4fa9-935b-a5e9b7dd3f09",
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef0b43b6-d7fb-4ef5-912a-da327372b26a",
   "metadata": {},
   "source": [
    "##### Retrieve (weigthed) aggregate statisitics of the elevation for each of the threshold classess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d960126b-1a6c-45d5-9ada-f02f81704f1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- Weighted Functions -----\n",
    "def weighted_median(values, weights):\n",
    "    sorted_idx = np.argsort(values)\n",
    "    values_sorted = np.array(values)[sorted_idx]\n",
    "    weights_sorted = np.array(weights)[sorted_idx]\n",
    "    cum_weights = np.cumsum(weights_sorted)\n",
    "    cutoff = 0.5 * sum(weights_sorted)\n",
    "    return values_sorted[np.searchsorted(cum_weights, cutoff)]\n",
    "\n",
    "def weighted_std(values, weights):\n",
    "    average = np.average(values, weights=weights)\n",
    "    variance = np.average((values - average)**2, weights=weights)\n",
    "    return np.sqrt(variance)\n",
    "\n",
    "# ----- Unweighted Aggregations -----\n",
    "unweighted = (\n",
    "    df.groupby('first_raster')['second_raster']\n",
    "    .agg(\n",
    "        min_elevation='min',\n",
    "        max_elevation='max',\n",
    "        std_elevation='std',\n",
    "        median_elevation='median',\n",
    "        sum_elevation='sum',\n",
    "        count_elevation='count',\n",
    "        predominant_elevation=lambda x: x.mode().iloc[0] if not x.mode().empty else np.nan\n",
    "    )\n",
    "    .reset_index()\n",
    ")\n",
    "\n",
    "# ----- Weighted Average -----\n",
    "df['weighted_elevation'] = df['second_raster'] * df['pixel_count']\n",
    "weighted_avg = (\n",
    "    df.groupby('first_raster')\n",
    "    .apply(lambda g: (g['second_raster'] * g['pixel_count']).sum() / g['pixel_count'].sum())\n",
    "    .reset_index(name='weighted_avg_elevation')\n",
    ")\n",
    "\n",
    "# ----- Weighted Median & Std -----\n",
    "def group_weighted_stats(group):\n",
    "    values = group['second_raster'].values\n",
    "    weights = group['pixel_count'].values\n",
    "    return pd.Series({\n",
    "        'weighted_median_elevation': weighted_median(values, weights),\n",
    "        'weighted_std_elevation': weighted_std(values, weights)\n",
    "    })\n",
    "\n",
    "weighted_stats = df.groupby('first_raster').apply(group_weighted_stats).reset_index()\n",
    "\n",
    "# ----- Merge All Together -----\n",
    "final = (\n",
    "    unweighted\n",
    "    .merge(weighted_avg, on='first_raster')\n",
    "    .merge(weighted_stats, on='first_raster')\n",
    ")\n",
    "\n",
    "# Display or export final DataFrame\n",
    "print(final)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38cfbe0a-9dc3-40b0-8410-f3b75c81a1e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#store dataframe as csv file - without index column\n",
    "final.to_csv(work_dir+'/agg_stats.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dce6cc8-4e7c-4504-91bc-1b5d223313cf",
   "metadata": {},
   "source": [
    "#### Cumulative hypsometric curve\n",
    "A cumulative hypsometric curve, also simply called a hypsometric curve, is a graph that shows the proportion of an area (like a watershed or the Earth's surface) located at or below a certain elevation. It's a cumulative distribution function of elevations, essentially a histogram of elevation data. This curve helps geomorphologists and hydrologists understand the distribution of elevations and compare the geomorphology of different areas. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "398ac6c8-a5dc-4e8a-944a-e0adfd52a80e",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_catchment_map = ilwis.RasterCoverage('catchment.mpr') \n",
    "selected_catchment = 49 #use catchment-id 50  -note difference of '1' due to index offset\n",
    "elevation_map = raster_sel2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "016dc555-8621-41ed-91e7-eadd29113e87",
   "metadata": {},
   "outputs": [],
   "source": [
    "c_temp = ilwis.do('mapcalc','iff(@1==' + str(selected_catchment) + ',1,?)',merged_catchment_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0712d66-2af5-4da5-8784-2c23283a91c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "c_temp.store('sel_catch.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63513c1f-5ce8-4cf7-ba63-1a56f8809691",
   "metadata": {},
   "outputs": [],
   "source": [
    "ct = ilwis.do('cross',c_temp, elevation_map, 'ignoreundef')#note first output is table, second is map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a995919-f12b-4656-b61f-d1ae38a85664",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ilwis cross table to pandas dataframe\n",
    "data = {}\n",
    "for col in list(ct.columns()):\n",
    "    data[col] = list(ct.column(col))\n",
    "df = pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3532591c-053f-49b5-b180-55b0909dfaaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "#sort on elevation - from low to high elevation\n",
    "df_sorted = df.sort_values(by='second_raster')\n",
    "df_sorted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74a337ff-1d98-43ec-aa29-bd406e43a7e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_pixels = df_sorted['pixel_count'].sum()\n",
    "df_sorted['Cum_pixels'] = df_sorted['pixel_count'].cumsum()\n",
    "df_sorted['Cum_area'] = df_sorted['Cum_pixels']* total_pixels\n",
    "df_sorted['Percentage'] = (df_sorted['pixel_count'] / total_pixels) * 100\n",
    "df_sorted['Cum_perc'] = df_sorted['Percentage'].cumsum().round(3)\n",
    "\n",
    "total_percent = df_sorted['Percentage'].sum()\n",
    "\n",
    "if total_percent > 100:\n",
    "    # Find the index of the row with the max percentage\n",
    "    max_idx = df_sorted['Cum_perc'].idxmax()\n",
    "    # Reduce the max value by the excess\n",
    "    df_sorted.loc[max_idx, 'Cum_perc'] -= round(total_percent - 100, 2)\n",
    "\n",
    "#print(\"Total number of pixels:\", total_pixels)\n",
    "df_sorted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c79f72d0-fb4e-4112-84e1-7d24542ba1b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot(df_sorted['second_raster'], df_sorted['Cum_perc'], label = 'elevation-percentage')\n",
    "plt.xlabel(\"Elevation\")\n",
    "plt.ylabel(\"Percentage\")\n",
    "plt.title(\"Cumulative hypsometric curve of selected catcment\")\n",
    "plt.grid(True)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34409e89-10e5-40d3-8b75-60dcaa3e1f08",
   "metadata": {},
   "outputs": [],
   "source": [
    "#store dataframe as csv file - without index column\n",
    "df_sorted.to_csv(work_dir+'/ct_hypso.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0520a6e4-0b8f-4ce8-a4aa-7652d08cd163",
   "metadata": {},
   "source": [
    "### Plot the longest flowpath longitudinal profile \n",
    "\n",
    "In order to do this the segment file created before needs to be converted to raster, crossed with the elevation model and the cross table obtained shows per pixel the corresponding elevation. From the segment attribute table the length of the segment can be extracted and when creating the plot the elevation values can be evenly distributed over the length to create the profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc8beb2e-9d49-444c-b305-9d099d28ae0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import longest flow path vector map\n",
    "longest_flow_path = ilwis.FeatureCoverage('catchmerge_long_flowpath.mps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df25a9a2-f99f-4392-a374-3e9f75b45f44",
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert longest flow path segment to raster\n",
    "longest_flow_path_r = ilwis.do('line2raster', longest_flow_path, dem.geoReference())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac16b9e5-7944-49e0-8463-845a91ea36e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ct_lfp = ilwis.do('cross',longest_flow_path_r, dem, 'ignoreundef')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1499080a-489b-433f-9baf-909c8f8b9a43",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ilwis cross table to pandas dataframe\n",
    "data = {}\n",
    "for col in list(ct_lfp.columns()):\n",
    "    data[col] = list(ct_lfp.column(col))\n",
    "df = pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65424a22-f7b9-4a97-817b-45b67864e407",
   "metadata": {},
   "outputs": [],
   "source": [
    "#sort on elevation - from low to high elevation\n",
    "df_sorted = df.sort_values(by='second_raster')\n",
    "df_sorted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ebb4265-0fed-4545-b9b6-b5c6f7e44242",
   "metadata": {},
   "outputs": [],
   "source": [
    "#load the catchmerge_long_flowpath table to get length of longest flow path\n",
    "table = ilwis.Table(\"catchmerge_long_flowpath.tbt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a16fe11-98d5-4ded-b3b1-63b2804930b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "recCount = table.recordCount()\n",
    "colCount = table.columnCount()\n",
    "columns = table.columns()  # (listing of columns in table: 'column1','column2',...)\n",
    "\n",
    "print(recCount)\n",
    "print(colCount)\n",
    "print(columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fd130bb-7ea6-4177-a1fc-c5919da96209",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get first record of column 'Length'\n",
    "seg_length = table.cell(\"Length\", 0)\n",
    "print(seg_length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4698738c-f556-417d-88a0-7dbe0db37fe2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Elevation values\n",
    "elevation = np.sort(df['second_raster'].values)  #from low to high\n",
    "\n",
    "# Create distance x-axis\n",
    "total_length = seg_length  # meters\n",
    "\n",
    "distance = np.linspace(0, total_length, len(elevation)) #distribute the elevation values evenly along the length of segment\n",
    "\n",
    "distance_km = distance / 1000 #transform to km\n",
    "\n",
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot(distance_km, elevation, linewidth=2, label = 'flow-path elevation')\n",
    "\n",
    "plt.xlabel('Distance (km)')\n",
    "plt.ylabel('Elevation (m)')\n",
    "plt.title('Longest Flow-path Longitudinal Profile')\n",
    "\n",
    "plt.legend()\n",
    "plt.grid(True, alpha=0.3)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22784044-6871-4936-8fca-d4a69fd9cc42",
   "metadata": {},
   "source": [
    "#### Calculate DEM derived compound indexes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5ca20c8-290e-4fe1-b401-1fbe96f910f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#assign input maps\n",
    "dem = ilwis.RasterCoverage('dem_fill.mpr')\n",
    "dem_fa = ilwis.RasterCoverage('dem_fa.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2984e02d-127f-4d35-884d-5a2ce27b3895",
   "metadata": {},
   "outputs": [],
   "source": [
    "#assume metric coordinate system, if lat-lon then first transform to metric / UTM\n",
    "dx = ilwis.do('linearrasterfilter', dem, 'dfdx')\n",
    "dy = ilwis.do('linearrasterfilter', dem, 'dfdy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b05c1afa-6a09-4636-abff-f869eebd67b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the pixelsize from the image dimensions\n",
    "dimensions = dem.envelope()\n",
    "print(dimensions)\n",
    "pix_size = (dimensions.maxCorner().x - dimensions.minCorner().x) / dem.size().xsize\n",
    "print(pix_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "569dc13b-1681-4ab6-9619-057c8738a36c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate slope map in percentage\n",
    "slope_percent = ilwis.do('mapcalc','100*(sqrt(pow(@1,2)+(pow(@2,2)))/'+ str(pix_size) + ')',dx,dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60d3e6e0-beb2-47c9-b042-87190d4a2845",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_percent.store('slope_percent.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd9b8bc7-b804-47c4-9c96-5385a0f2f494",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate slope map in degree - ilwis method\n",
    "slope_degree = ilwis.do('mapcalc','180/3.141592654*atan(@1/100)',slope_percent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcc61889-f63b-493d-8fda-a85bf72dda4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_degree.store('slope_degree.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5649899-4915-4849-94bd-a0b1cd6f1af8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate slope in degree - python numpy method\n",
    "slope_percent_2np = np.fromiter(iter(slope_percent), np.float64, slope_percent.size().linearSize()) \n",
    "slope_percent_2np = slope_percent_2np.reshape((slope_percent.size().ysize, slope_percent.size().xsize))\n",
    "slope_degrees_map = np.degrees(np.arctan(slope_percent_2np / 100))\n",
    "\n",
    "grf = ilwis.GeoReference(slope_percent.coordinateSystem(), slope_percent.envelope(), slope_percent.size())\n",
    "dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 1000.0, 0.001))\n",
    "rcNew = ilwis.RasterCoverage()\n",
    "rcNew.setSize(ilwis.Size(slope_percent.size().xsize, slope_percent.size().ysize, 1))\n",
    "rcNew.setGeoReference(grf)\n",
    "rcNew.setDataDef(dfNum)\n",
    "\n",
    "rcNew.array2raster(slope_degrees_map)\n",
    "rcNew.store('slope_degree_numpy.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc29b113-6444-4da8-8195-5e3f0f0b8a7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#remove 0-slope values, change these to 0.1, else compound index calculation will report division by 0\n",
    "slp_temp = ilwis.do('mapcalc', 'iff(@1>0,@1,0.1)',rcNew)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55775f6e-2bd6-48f7-aa59-f59c53c31883",
   "metadata": {},
   "outputs": [],
   "source": [
    "slp_temp.store('slope_no0.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3ff60fa-fd7a-4470-a1f7-f3798d075d21",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate the wetness index\n",
    "wet_index = ilwis.do('mapcalc','ln((@1*' + str(pix_size * pix_size) + ') / tan(@2 * 3.141592654 / 180.0))',dem_fa,slp_temp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81a95dfa-4252-4f8e-a715-3acd8bfba332",
   "metadata": {},
   "outputs": [],
   "source": [
    "wet_index.store ('wetness_index.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faaed446-5ed6-482c-9859-a3ef1913dd22",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate the stream power index\n",
    "stream_index = ilwis.do('mapcalc','((@1*' + str(pix_size * pix_size) + ') * tan(@2 * 3.141592654 / 180.0))',dem_fa,slp_temp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "712e2fe3-08df-452f-8abf-b048c5643439",
   "metadata": {},
   "outputs": [],
   "source": [
    "stream_index.store ('stream_power_index.mpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5734cc60-00e9-41ca-a4ec-0e2ecf7eb0b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate sediment transport index\n",
    "sediment_index = ilwis.do('mapcalc','(pow(@1*' + str(pix_size * pix_size) + '/22.13, 0.6) * (pow(sin(@2 * 3.141592654 / 180.0)/0.0896,1.3)))',dem_fa,slp_temp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1716d51-88ee-4b24-a967-6bd2543952bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "sediment_index.store ('sediment_stream_index.mpr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a44cfd0d-595c-4b75-bfb1-8ea37732c054",
   "metadata": {},
   "source": [
    "### Develop a SCS dimensionless Unit Hydrograph and apply this to the merged catchment and calculate the discharge for a given rainfall event\n",
    "\n",
    "The direct runoff hydrograph resulting from one unit depth of excess rainfall (e.g., 1 cm or 1 inch), uniformly distributed over a watershed at a constant rate for a specified duration. Therefore the unit hydrograph represents the runoff response to:  1 cm of effective rainfall occurring uniformly over the watershed during one hour.\n",
    "\n",
    "The key properties are: \n",
    "\n",
    "+ Volume conservation: V=A×d\n",
    "    where:\n",
    "\n",
    "        A = watershed area\n",
    "        d = unit rainfall excess depth\n",
    "  \n",
    "+ Linearity: If rainfall excess doubles, discharge doubles.\n",
    "\n",
    "+ Superposition: Multiple rainfall pulses are handled through convolution.\n",
    "\n",
    "For a more detailed description of the Unit Hydrograph concept see: https://www.hec.usace.army.mil/confluence/hmsdocs/hmstrm/transform/scs-unit-hydrograph-model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b94c168-42db-42b4-8357-c9038e62b116",
   "metadata": {},
   "source": [
    "### Retrieve a small mountainous catchment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95e91f8d-0087-4656-b211-d871786a4142",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First import the outlet point map\n",
    "outlet = ilwis.FeatureCoverage('UH_location.mpp') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "419df910-871e-4777-bc62-4be3eb136214",
   "metadata": {},
   "outputs": [],
   "source": [
    "#perform the catchment merge\n",
    "CatchMerge_outlet = ilwis.do('MapCatchmentMergeWithOutlet','dr_net.mpr', dem_fd,dem_fa,fillsink,outlet,'no')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5f5be72-2dd5-45ba-b9fa-cb9cec80c8be",
   "metadata": {},
   "outputs": [],
   "source": [
    "#store the results\n",
    "#note different polygon and vector maps created\n",
    "CatchMerge_outlet[0].store('UH_catchmerge_polygon.mpa') \n",
    "CatchMerge_outlet[1].store('UH_catchmerge_raster.mpr') \n",
    "CatchMerge_outlet[2].store('UH_catchmerge_long_flowpath.mps')\n",
    "CatchMerge_outlet[3].store('UH_catchmerge_segments.mps')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ab7deae-9c44-47e1-86b0-565a5e716cb6",
   "metadata": {},
   "source": [
    "#### Read the required data from merged catchment table \"catchmerge_polygon.tbt\":\n",
    "\n",
    "+ Catchment area (km²) = \"CatchmentArea\"\n",
    "+ Longest flow path (m) = \"LongestFlowPathLength\"\n",
    "+ Longest drainage length (m) = \"LongestDrainageLength\"\n",
    "+ Elevation divide (m) = \"LDPUpstreamElevation\"\n",
    "+ Elevation outlet (m) = \"OutletElevation\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "244ba1f9-4b89-4da1-81fc-9e7c4d39b2ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "#load the catchmerge_long_flowpath table to get length of longest flow path\n",
    "UH_table = ilwis.Table(\"UH_catchmerge_polygon.tbt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33336476-8405-412e-9038-63f21a334fca",
   "metadata": {},
   "outputs": [],
   "source": [
    "recCount = UH_table.recordCount()\n",
    "colCount = UH_table.columnCount()\n",
    "columns = UH_table.columns()  # (listing of columns in table: 'column1','column2',...)\n",
    "\n",
    "print(recCount)\n",
    "print(colCount)\n",
    "print(columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f1beef1-79b1-40dc-8037-d6b8546da291",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get first record of required columns\n",
    "A_m2 = UH_table.cell(\"CatchmentArea\",0)\n",
    "L = UH_table.cell(\"LongestFlowPathLength\",0)\n",
    "drainage_length = UH_table.cell(\"LongestDrainageLength\",0)\n",
    "elev_divide = UH_table.cell(\"LDPUpstreamElevation\",0)\n",
    "elev_outlet = UH_table.cell(\"OutletElevation\",0)\n",
    "\n",
    "print(\"Catchment area (km2) =\",A_m2/1000000)\n",
    "print(\"Longest flow path (m) =\",L) \n",
    "print(\"Longest drainage length (m) =\",drainage_length)\n",
    "print(\"Elevation divide (m) =\",elev_divide)\n",
    "print(\"Elevation outlet (m) =\",elev_outlet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f0df004-ae8f-4cdd-bd53-17abd80f1d06",
   "metadata": {},
   "outputs": [],
   "source": [
    "H = elev_divide - elev_outlet  # elevation difference (m)\n",
    "S = H / L                      # watershed slope (m/m)\n",
    "\n",
    "print(\"elevation difference (m) =\",H)\n",
    "print(\"watershed slope (m/m) =\",S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07b00a2a-4386-404b-a1ed-7279ffdce25d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Time of Concentration (The time required for runoff from the hydraulically most distant point in the watershed to reach the outlet)\n",
    "# Tc according to Kirpich (1940) - for small watersheds\n",
    "L_feet = drainage_length/0.3048\n",
    "\n",
    "#Tc_min = 0.0195 * (L ** 0.77) * (S ** (-0.385)) # SI units\n",
    "Tc_min = 0.0078 * (L_feet ** 0.77) * (S ** (-0.385)) # U.S. customary form (feet)\n",
    "Tc_K_hr = Tc_min / 60.0\n",
    "\n",
    "# TC according to California Culverts Practice (1942):\n",
    "dr_length_miles = (L/1000)/1.6093\n",
    "#print(dr_length_miles)\n",
    "elev_dif_feet = H/0.3048\n",
    "#print(elev_dif_feet)\n",
    "\n",
    "Tc_CCP_hr = ((11.9 * dr_length_miles**3) / elev_dif_feet) ** 0.385\n",
    "\n",
    "Tc_Avg_hr = (Tc_K_hr+Tc_CCP_hr)/2\n",
    "\n",
    "\n",
    "print(\"\\nKirpich and CCP Tc approximation\")\n",
    "print(\"--------------------\")\n",
    "print(\"TC-Kirpich (hour) =\",Tc_K_hr)\n",
    "print(\"TC-CCP (hour) =\", Tc_CCP_hr)\n",
    "print(\"Tc-Average (hour) =\", Tc_Avg_hr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57bd2b4b-0c8e-4297-a932-6a4384036e01",
   "metadata": {},
   "outputs": [],
   "source": [
    "#SCS dimensionless Unit Hydrograph \n",
    "#for the coordinates / ratios of the dimensionless Unit Hydrograph see also: https://www.nohrsc.noaa.gov/technology/gis/uhg_manual.html \n",
    "\n",
    "tr = 1 #time step in hrs\n",
    "Depth = 1 #rainfall depth in cm/hr\n",
    "Tc = Tc_Avg_hr #time of concentration in hrs\n",
    "tp = 0.6 * Tc #lag time in hrs (time between the centroid of effective rainfall and the peak of the runoff hydrograph-the average watershed response)\n",
    "Tp = Depth / 2 + tp #time to peak in hrs\n",
    "Area = A_m2/1000000 #catchment area in km2\n",
    "\n",
    "Qp = 2.08 * Area / Tp #Peak discharde [m3/(s.cm)]\n",
    "\n",
    "print('Lag time (hr) =', tp)\n",
    "print('Time to peak (hr)=', Tp)\n",
    "print('Catchment area (km2)=', Area)\n",
    "print('Peak discharge (m3/sec)= ', Qp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6f2cca1-bdb1-43c9-9568-9d64f88fdbce",
   "metadata": {},
   "outputs": [],
   "source": [
    "#dimensionless UH shape\n",
    "#t/TP - time ratios (x-axis)\n",
    "t_Tp = np.array([0.00, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.20, 2.40, 2.60, 2.80, 3.00, 3.20, 3.40, 3.60, 3.80, 4.00, 4.50, 5.00])\n",
    "\n",
    "#q/Qp - discharge ratios (y-axis)\n",
    "q_Qp = np.array([0.00, 0.03, 0.10, 0.19, 0.31, 0.47, 0.66, 0.82, 0.93, 0.99, 1.00, 0.99, 0.93, 0.86, 0.78, 0.68, 0.56, 0.46, 0.39, 0.33, 0.28, 0.207, 0.147, 0.107, 0.077, 0.055, 0.04, 0.029, 0.021, 0.015, 0.011, 0.005, 0])\n",
    "\n",
    "# Plot\n",
    "plt.figure(figsize=(6, 4))\n",
    "\n",
    "plt.plot(t_Tp, q_Qp, marker='o', linewidth=2)\n",
    "\n",
    "plt.xlabel('t/Tp')\n",
    "plt.ylabel('Q/Qp')\n",
    "plt.title('Standard SCS Dimensionless Unit Hydrograph')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.xlim(0, 5)\n",
    "plt.ylim(0, 1.05)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11929b11-1404-4751-932c-256fa888b81b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate t (hrs)\n",
    "t = t_Tp*Tp\n",
    "print(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52a6678f-032d-4a6f-88fb-9e2fe9aec114",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate discharge [m3/(s.cm)]\n",
    "discharge = q_Qp * Qp\n",
    "print(discharge)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31a2a7ad-6945-4171-b267-f2db1ece5e8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#transform to 1 hr time steps (note that t as provided above is not having hourly intervals)\n",
    "t_hr = np.arange(0, np.ceil(t.max()) + 1, 1.0)\n",
    "print(len(t))\n",
    "print(t_hr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc4444ac-44e7-4985-942a-44d97c7c5556",
   "metadata": {},
   "outputs": [],
   "source": [
    "#conduct the interpolation of discharge for the correcponding hourly timesteps\n",
    "UH_1hr = np.interp(t_hr, t, discharge)\n",
    "print(len(UH_1hr))\n",
    "print(UH_1hr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d94d94ba-ac42-40bd-876b-06607b916fd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6, 4))\n",
    "\n",
    "plt.plot(t_hr, UH_1hr, marker='o', linewidth=2, markersize=5, label='1-hour Unit Hydrograph')\n",
    "\n",
    "plt.xlabel('Time (hr)')\n",
    "plt.ylabel('Discharge (m³/s)')\n",
    "plt.title('1-Hour SCS Unit Hydrograph of selected catchment')\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "872b677a-4702-4950-9596-547348601dfe",
   "metadata": {},
   "source": [
    "#### Rainfall excess. \n",
    "\n",
    "In unit hydrograph theory, rainfall excess is the portion of rainfall that actually becomes direct runoff at the watershed outlet.\n",
    "\n",
    "Not all rainfall contributes immediately to streamflow. Some is lost to:\n",
    "\n",
    "+ Infiltration into the soil\n",
    "+ Interception by vegetation\n",
    "+ Depression storage (water trapped in puddles, surface irregularities)\n",
    "+ Evaporation and evapotranspiration\n",
    "\n",
    "The rainfall remaining after these abstractions is called: Rainfall Excess=Total Rainfall−Losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e7db88d-8704-4c6b-9d93-58aa6dbceb1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rainfall excess (hourly, for 4 consecutive hours)\n",
    "rain = np.array([0.03, 1.77, 0.22, 0.02])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb882421-6d74-4344-ad92-a619eecc511f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# conduct the convolution\n",
    "Q = np.convolve(rain, UH_1hr)\n",
    "time_q = np.arange(len(Q))\n",
    "\n",
    "#show resulting table\n",
    "df = pd.DataFrame({\"Hour\": time_q,\"Discharge_m3s\": Q})\n",
    "\n",
    "print(\"\\n=========================\")\n",
    "print(\"CONVOLVED HYDROGRAPH\")\n",
    "print(\"=========================\")\n",
    "print(df)\n",
    "\n",
    "peak_idx = np.argmax(Q)\n",
    "peak_q = np.max(Q)\n",
    "\n",
    "print(\"\\nPeak discharge:\")\n",
    "print(\"Qpeak (m3/s)=\", peak_q)\n",
    "print(\"Time (hr) =\", peak_idx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01900311-4011-44ad-96e3-1cec8d88177f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the results\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8), sharex=False, gridspec_kw={\"height_ratios\": [1, 2]})\n",
    "\n",
    "\n",
    "# Rainfall Hyetograph - showing rainfall versus time\n",
    "rain_time = np.arange(1, len(rain) + 1)\n",
    "\n",
    "ax1.bar(rain_time, rain, width=0.8)\n",
    "\n",
    "ax1.invert_yaxis()\n",
    "\n",
    "ax1.set_ylabel(\"Rainfall Excess\")\n",
    "ax1.set_title(\"Rainfall Excess Hyetograph\")\n",
    "ax1.grid(True, alpha=0.3)\n",
    "\n",
    "# Runoff Hydrograph - showing streamflow/discharge versus time\n",
    "\n",
    "ax2.plot(time_q, Q, linewidth=2.5, label=\"Convolved Hydrograph\")\n",
    "\n",
    "ax2.axvline(Tc_Avg_hr, linestyle=\"--\", linewidth=2, label=\"Tc-Average (hour)\")\n",
    "ax2.scatter(peak_idx, peak_q, s=60, zorder=5,label=f\"Peak Discharge = {peak_q:.2f} m3/s\")\n",
    "\n",
    "ax2.set_xlabel(\"Time (hours)\")\n",
    "ax2.set_ylabel(\"Discharge (m3/s)\")\n",
    "ax2.set_title(\"Direct Runoff Hydrograph from Unit Hydrograph Convolution\")\n",
    "\n",
    "ax2.grid(True, alpha=0.3)\n",
    "ax2.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "363bd352-8d2a-4fba-9316-87f3dd15b406",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check: estimated time of the peak discharge\n",
    "\n",
    "#sum rainfall\n",
    "P_total_cm =  np.sum(rain) #note dominant rainfall event during the second hour\n",
    "print('Total excess rainfall (cm)=', f\"{P_total_cm:.2f}\")\n",
    "\n",
    "#midpoint of each 1-hour interval\n",
    "time_mid = np.array([0.5, 1.5, 2.5, 3.5])\n",
    "\n",
    "#rainfall centroid time\n",
    "tc = np.sum(rain * time_mid) / np.sum(rain)\n",
    "print('rainfall centroid time (hr)',f\"{tc:.2f}\")\n",
    "\n",
    "Tlag = 0.6 * tc\n",
    "print('lag time (hr)', f\"{Tlag:.2f}\")\n",
    "\n",
    "#Expected peak timing\n",
    "tpeak = tc + Tlag\n",
    "print('hydrograph peak expected around (hr):', f\"{tpeak:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "717ab577-8d46-41f2-a490-32df37ef7e39",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check: Unit Hydrograph consistency check\n",
    "\n",
    "# Unit hydrograph time step (1 hour)\n",
    "dt = 3600.0  # seconds\n",
    "\n",
    "# Expected runoff volume for 1 cm excess rainfall\n",
    "rainfall_excess_m = 0.01  # 1 cm = 0.01 m\n",
    "expected_volume = rainfall_excess_m * A_m2\n",
    "\n",
    "# Volume represented by the Unit Hydrograph\n",
    "UH_volume = np.sum(UH_1hr) * dt\n",
    "\n",
    "# Compare both\n",
    "percent_error = ((UH_volume - expected_volume) / expected_volume) * 100\n",
    "\n",
    "print(\"\\n==============================\")\n",
    "print(\"UNIT HYDROGRAPH CHECK\")\n",
    "print(\"==============================\")\n",
    "\n",
    "print(\"Catchment area (km2) =\", A_m2/1000000)\n",
    "print(\"1 cm runoff volume (m3) =\", expected_volume)\n",
    "print(\"Volume under Unit Hydrograph (m3) =\", UH_volume)\n",
    "print(\"Difference (m3) =\", UH_volume - expected_volume)\n",
    "print(\"Percent error (%) =\", f\"{percent_error:.2f}\")\n",
    "\n",
    "if abs(percent_error) < 5:\n",
    "    print(\"\\n Unit Hydrograph is consistent with a 1-cm UH.\")\n",
    "else:\n",
    "    print(\"\\n Check UH scaling or watershed area.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8cafc43-bdb0-4ea9-93b3-b954948f844d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check: Total volume check\n",
    "P_total_m =  np.sum(rain)/100\n",
    "print('Total excess rainfall (m) =', P_total_m)\n",
    "print('Catchment area (m2) =', A_m2)\n",
    "\n",
    "TotalVolume_basin = P_total_m * A_m2\n",
    "print('Total basin volume (m3) =',TotalVolume_basin)\n",
    "\n",
    "runoff_volume = np.sum(Q) * 3600\n",
    "print('Total runoff volume (m3) =', runoff_volume)\n",
    "\n",
    "# Compare both\n",
    "Vpercent_error = ((TotalVolume_basin - runoff_volume) / runoff_volume) * 100\n",
    "\n",
    "print(\"Difference (m3) =\", TotalVolume_basin - expected_volume)\n",
    "print(\"Percent error (%) =\", f\"{Vpercent_error:.2f}\")\n",
    "\n",
    "if abs(percent_error) < 5:\n",
    "    print(\"\\n Total runoff volume is consistent with total basin precipitation volume.\")\n",
    "else:\n",
    "    print(\"\\n Check  rainfall or discharge volumes used.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7bfb6f3-6dbc-4cd3-a428-2f3c2bc60473",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export to csv\n",
    "df.to_csv(\"convolved_hydrograph.csv\", index=False)"
   ]
  }
 ],
 "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
}