Elevation data processing for Hydrology

Notebook prepared by Ben Maathuis, Lichun Wang and Bas Retsios. ITC-University of Twente, Enschede. The Netherlands

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.

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.

# if relevant update ilwispy by uncommenting the line below, use at least ilwispy version of 20260708
#!pip install ilwis --upgrade

Load libraries

import os
import ilwis
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import LineString

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

#check your ilwispy version used
ilwis.version()

'1.0 build 20260708'

work_dir = os.getcwd()+'/DEM_processing'

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

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

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

cd_folder = os.getcwd() #check current directory folder
print(cd_folder)

d:\jupyter\notebook_scripts\Special\DEM_processing

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

d:\jupyter\notebook_scripts\Special\DEM_processing/DEM_processing

#import elevation model
dem = ilwis.RasterCoverage('Thiba_DEM_UTM.tif')

#retrieve meta data from imported map
print(dem.size())
coordSys = dem.coordinateSystem()
print(coordSys.toWKT())
print()
print(dem.envelope())

Size(1226, 2396, 1)
PROJCS["thiba_dem_utm.tif",GEOCS["thiba_dem_utm.tif",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Transverse_Mercator"],PARAMETER[,No]PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],PARAMETER["scale",0.9996],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",39],UNIT[meter,1.0]]

299657.500000 9911539.000000 336437.500000 9983419.000000

#retrieve some basic statistics from imported map
stats_dem = dem.statistics(ilwis.PropertySets.pHISTOGRAM)
#print(stats_dem.histogram())
dem_min = stats_dem[ilwis.PropertySets.pMIN] # minimum value on the map
dem_max = stats_dem[ilwis.PropertySets.pMAX] # maximum value on the map
print(dem_min)
print(dem_max)

1047.5528564453125
5158.880859375

#transform imported map from ilwis format into a numpy array using the (ilwispy) pixel itereator
dem_2np = np.fromiter(iter(dem), np.float64, dem.size().linearSize()) 
#now we overwrite the initial variable created
dem_2np = dem_2np.reshape((dem.size().ysize, dem.size().xsize))

#display the numpy array created in the previous step using matplotlib 
fig = plt.figure(figsize =(10, 7))
plt.imshow(dem_2np, interpolation='none', vmin=dem_min, vmax=dem_max, cmap='jet')

plt.axis('on')
plt.colorbar(shrink=0.45)
plt.title('DEM - south of Mount Kenya');

Fill Sinks

MapFillSinks(inputraster, fill|cut)

• Inputraster is the input DEM map
• fill|cut is a method indicate to remove local depressions or cut terrain from a DEM raster

#first we apply the cut routine followed by the fill
cut = ilwis.do('MapFillSinks', dem, 'cut')
fillsink = ilwis.do('MapFillSinks', cut, 'fill')

#store your results obtained, check the maps produced using the ilwis desktop version
fillsink.store('dem_fill.mpr')

Flow Direction

MapFlowDirection(inputraster, slope|height, yes|no)

• Inputraster is the input raster map (sink-free DEM)
• slope|height is a method to indicate how flow direction should be calculated (the steepst slope or the smallest height)
• yes|no is an option for the Parallel drainage correction algorithm

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

dem_fd.store('dem_fd.mpr')

Flow Accumulation

MapFlowAccumulation(inputraster)

• Inputraster is the input flow direction map

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

dem_fa.store('dem_fa.mpr')

Internal relief map calculation

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)

• internalrelief(inputraster,filtersize)
• Inputraster is the input DEM map
• Filtersize is an odd value to define the filter size

int_relief = ilwis.do('InternalRelief', fillsink, 5)

int_relief.store ('int_relief.mpr') 

Variable threshold computation

• variablethreshold(inputraster,filtersize,NrOfClasses,UpperBounds_and_ThresholdVals)
• inputraster is the input DEM map
• filtersize is a value for the size of the window
• NrOfClasses is a value for the number of classes
• UpperBounds_and_ThresholdVals is a string in which, for each class at a time, the upper boundary values and threshold values are specified.

#create threshold raster, 5 is the filtersize, 4 are the number of threshold classes used
thres4 = ilwis.do('VariableThreshold', dem, 5, 4, '10, 15000, 25, 10000, 60, 5000, 1000, 1500')

thres4.store('thres4.mpr')

Drainage Network extraction(1) - multiple thresholds

ExtractDrainageUseThresholdRaster(inputraster1,inputraster2,inputraster3)

• This command is used to extract the drainage according to the input threshold map
• Inputraster1is the flow accumulation map
• Inputraster2 is the threshold map
• Inputraster3 is the flow direction map

#drainage raster extraction using raster threshold map (using multiple flow accumulation thresholds)
dr_thres4 = ilwis.do('ExtractDrainageUseThresholdRaster', dem_fa, thres4, dem_fd)

dr_thres4.store('drain_thres4.mpr')

Drainage Network extraction(2) - single threshold

ExtractDrainageUseThresholdValue(FlowAccumulationRaster,thresholdval)

• FlowAccumulationRaster is the flow accumulation map
• Thresholdval is a threshold value (integer > 0)

#drainage raster map creation using single flow accumulation threshold
dr_single7500 = ilwis.do('ExtractDrainageUseThresholdValue', dem_fa, 7500)

dr_single7500.store('dr_7500.mpr')

Drainage Network Ordering

DrainageNetworkOrdering(DEMmap,FlowDiractionmap,DrainageNetworkmap,MinimumDrainageLength)

• DEMmap is the input DEM map
• FlowDirectionmap is the input Flow Direction map
• DrainageNetworkmap is the input drainage network map
• MinimumDrainageLength is a value for the minimum length (in m)

This operation will generate the output drainage network raster and segment map. The output attribute table will be stored automatically.

outTable, outSegment, outDNO = ilwis.do('DrainageNetworkOrdering',fillsink, dem_fd, dr_thres4, 500)

outDNO.store('dr_net.mpr')
outSegment.store('dr_net.mps')

#store segments as shape file
outSegment.store('dr_net1.shp', 'ESRI Shapefile', 'gdal') 

Catchment extraction

MapCatchmentExtraction(DrainageNetworkOrderingMap,FlowDirectionMap)

• DrainageNetworkOrderingMap is the input drainage network ordering map
• FlowDirectionMap is the input flow direction map

It produces the output catchment raster and polygon map. The attribute table will be stored automaticcaly.

catchm = ilwis.do('MapCatchmentExtraction', outDNO, dem_fd) 

#note the layers for polygons and raster
catchm[0].store('catchment.mpa') 
catchm[1].store('catchment.mpr')

Catchment merge - using outlet location

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.

MapCatchmentMergeWithOutlet (DrainageNetworkOrderMap,FlowDiractionMap,FlowaccumulationMap,DEM,OutletPointMap,IncludeUndefinedPixels)

• DrainageNetworkOrderMap is the input drainage network oedering raster map
• FlowDirectionMap is the input flow direction map
• DEM is the input DEM map
• OutletPointMap is the input point map that contains the outlet location
• IncludeUndefinedPixels – yes|no, specify whether include the undefined pixels in the DEM into a catchment

# First import the outlet point map
outlet = ilwis.FeatureCoverage('outlet_location.mpp') 

# check imported point map
for i in range(outlet.attributeCount()):
    print(str(outlet[i]))  

Name: Name Columnindex: 0 Domain name: Text domain, Domain type: TextDomain
Name: coverage_key Columnindex: 1 Domain name: outlet_location.mpp, Domain type: ItemDomain, Range: indexedidentifierrange:pnt|1

CatchMerge_outlet = ilwis.do('MapCatchmentMergeWithOutlet','dr_net.mpr', dem_fd,dem_fa,fillsink,outlet,'no')

#note different polygon and vector maps created
CatchMerge_outlet[0].store('catchmerge_polygon.mpa') 
CatchMerge_outlet[1].store('catchmerge_raster.mpr') 
CatchMerge_outlet[2].store('catchmerge_long_flowpath.mps')
CatchMerge_outlet[3].store('catchmerge_segments.mps')

Catchment merge - using stream order method

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.

MapCatchmentMergeWithStreamOrder(DrainageNetworkOrderMap,FlowDirectionMap,FlowaccumulationMap,DEM,StreamOrderSystem,StreamOrderValue,ExtractOriginalOrder)

• DrainageNetworkOrderMap is the input drainage network ordering map
• FlowDirectionMap is the input flow direction map
• FlowaccumulationMap is the flow accumulation map
• DEM is the input DEM map
• StreamOrderSystem = strahler | shreve, specify whether you wish to use the Strahler or the Shreve ordering system.
• StreamOrderValue is the stream order value (integer value > 0);

CatchMerge_order3 = ilwis.do('MapCatchmentMergeWithStreamOrder', 'dr_net.mpr', dem_fd, dem_fa, fillsink, 'strahler', 3)

CatchMerge_order3[0].store('catchmerge_strahler3_polygon.mpa') 
CatchMerge_order3[1].store('catchmerge_strahler3_raster.mpr') 
CatchMerge_order3[2].store('catchmerge_strahler3_long_flowpath.mps')
CatchMerge_order3[3].store('catchmerge_strahler3_segment.mps') 

Overland Flow Length

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.

MapOverlandFlowLength(DrainageNetworkOrderMap,FlowDiractionMap)

• DrainageNetworkOrderMap is the input drainage network ordering map
• FlowDirectionMap is the input flow direction map

OverLandFlow = ilwis.do('MapOverlandFlowLength', outDNO, dem_fd)

OverLandFlow.store('overland_flow.mpr') 

- Flow Length to Outlet

MapFlowLength2Outlet(DrainageNetworkOrderMap,FlowDiractionMap)

• DrainageNetworkOrderMap is the input drainage network ordering map
• FlowDiractionMap is the input flow direction map

length2outlet = ilwis.do('MapFlowLength2Outlet', outDNO, dem_fd) 

length2outlet.store('length2outlet.mpr') 

Some other ilwis hydro dem processing routines examples

Calculate color shaded refief map

#create shadow filters: 'code=, dimensions = 3 by 3, followed by 9 filter values and eventually added by a gain factor, here 1'
shadow_w = 'code=3,3,-2 -1 2 -3 1 4 -2 -1 2,1'
shadow_n = 'code=3,3,-2 -3 -2 -1 1 -1 2 4 2,1'
shadow_nw = 'code=3,3,-3 -2 -1 -2 1 2 -1 2 4,1' # or use 'shadownw'
#shadow_nw = 'shadownw' #when applying a standard filter available in the ilwis filter library

#shadowW
shadow_west = ilwis.do('linearrasterfilter', dem, shadow_w)
shadow_west = ilwis.do('linearstretch',shadow_west, 5)
shadow_west = ilwis.do('setvaluerange', shadow_west, 0, 255, 1)
shadow_west = ilwis.do('mapcalc','iff(@1==?,0,@1)',shadow_west)
shadow_west_2np = np.fromiter(iter(shadow_west), np.ubyte, shadow_west.size().linearSize()) 
shadow_west_2np = shadow_west_2np.reshape((shadow_west.size().ysize, shadow_west.size().xsize))

#shadowN
shadow_north = ilwis.do('linearrasterfilter', dem, shadow_n)
shadow_north = ilwis.do('linearstretch',shadow_north, 5)
shadow_north = ilwis.do('setvaluerange', shadow_north, 0, 255, 1)
shadow_north = ilwis.do('mapcalc','iff(@1==?,0,@1)',shadow_north)
shadow_north_2np = np.fromiter(iter(shadow_north), np.ubyte, shadow_north.size().linearSize()) 
shadow_north_2np = shadow_north_2np.reshape((shadow_north.size().ysize, shadow_north.size().xsize))

#shadowNW
shadow_northwest = ilwis.do('linearrasterfilter', dem, shadow_nw)
shadow_northwest = ilwis.do('linearstretch',shadow_northwest, 5)
shadow_northwest = ilwis.do('setvaluerange', shadow_northwest, 0, 255, 1)
shadow_northwest = ilwis.do('mapcalc','iff(@1==?,0,@1)',shadow_northwest)
shadow_northwest_2np = np.fromiter(iter(shadow_northwest), np.ubyte, shadow_northwest.size().linearSize()) 
shadow_northwest_2np = shadow_northwest_2np.reshape((shadow_northwest.size().ysize, shadow_northwest.size().xsize))

nc_array = np.dstack((shadow_north_2np, shadow_northwest_2np, shadow_west_2np))

#load shape file and note extent
gdf = gpd.read_file(work_dir + "/dr_net1.shp")
print("Shapefile bounds:" ,gdf.total_bounds)
print()
#check coordinate system
print(gdf.crs)

Shapefile bounds: [ 299702.5 9912874.   336362.5 9983374. ]

PROJCS["unknown",GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",39],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

# plot DEM, NW-gray scale shaded relief and color shaded relief maps with overlay of drainage vector map

# common extent
extent = [299657.5, 336437.5, 9911539.0, 9983419.0]

# Create 3 subplots
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 7))

#DEM background
im1 = ax1.imshow(dem_2np, extent=extent, cmap='terrain', vmin=1500, vmax=4500)
gdf.plot(ax=ax1,facecolor="none", edgecolor="red",linewidth=0.5)

ax1.set_title("DEM")
ax1.set_xlabel("Easting (m)")
ax1.set_ylabel("Northing (m)")

cbar1 = fig.colorbar(im1, ax=ax1,shrink=0.8)
cbar1.set_label("Elevation (m)")

# Shaded Relief
im2 = ax2.imshow(shadow_northwest_2np, extent=extent, cmap='gray')
gdf.plot(ax=ax2, color='cyan', linewidth=0.5)

ax2.set_title("Shaded Relief")
ax2.set_xlabel("Easting (m)")
ax2.set_ylabel("Northing (m)")

# Color shaded RGB Composite
im3 = ax3.imshow(nc_array, extent=extent,zorder=1)
gdf.plot(ax=ax3, facecolor="none", edgecolor="yellow", linewidth=0.5, zorder=2)

ax3.set_title("RGB Composite")
ax3.set_xlabel("Easting (m)")
ax3.set_ylabel("Northing (m)")

fig.subplots_adjust(
    left=0.10,
    right=0.90,
    bottom=0.08,
    top=0.92,
    wspace=0.01
)
plt.show()

Conduct the processing, also using a looping procedure to stretch the initial filter results and store the final maps as an ilwis maplist

#once more compute the shadow maps
shadow_west = ilwis.do('linearrasterfilter', dem, shadow_w)
shadow_north = ilwis.do('linearrasterfilter', dem, shadow_n)
shadow_northwest = ilwis.do('linearrasterfilter', dem, shadow_nw)

rcNew = ilwis.RasterCoverage()
rcNew.setGeoReference(ilwis.GeoReference(shadow_west.coordinateSystem(), shadow_west.envelope(), shadow_west.size()))
rcNew.setSize(ilwis.Size(rcNew.size().xsize, rcNew.size().ysize, 3))
rcNew.setDataDef(ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 1)))

rcNew.addBand(0, shadow_west.begin())
rcNew.addBand(1, shadow_northwest.begin())
rcNew.addBand(2, shadow_north.begin())

#stretch the rasterbands contained in rcNew
multiple_stretch = []
multiple_bands = ilwis.do('selection',rcNew,"rasterbands(0..2)") # for specific bands use band number as "rasterbands(0,1,2,6)"
ls = ilwis.do('linearstretch',multiple_bands, 5)
mb_stretch = ilwis.do('setvaluerange', ls, 0, 255, 1)

#write the maplist to disk
mb_stretch.store('dem_col.mpl')

Compare your ilwis and matplotlib results

Create contour lines

#check the numpy array size and elevation range in DEM
print(dem_2np.shape)
print(np.nanmin(dem_2np))
print(np.nanmax(dem_2np))

(2396, 1226)
-1e+308
5158.880859375

#Remove the 'ilwis' no-data value(s) represented as -1e+308 to get the actual DEM data range
dem_adj = np.where(dem_2np == -1e308, np.nan, dem_2np)
print(np.nanmin(dem_adj))
print(np.nanmax(dem_adj))

1047.5528564453125
5158.880859375

#Use the extent as defined before

x = np.linspace(extent[0], extent[1], dem_adj.shape[1])
y = np.linspace(extent[3], extent[2], dem_adj.shape[0])

plt.figure(figsize=(10, 7))

# plot adjusted DEM
plt.imshow(dem_adj, extent=extent, cmap='terrain')

# create contours, interval = 250 m
levels = np.arange(np.floor(np.nanmin(dem_adj)/250)*250, np.ceil(np.nanmax(dem_adj)/250)*250 + 250, 250)

# plot contours
cs = plt.contour(x, y, dem_adj, levels=levels, colors='black',linewidths=0.5)

plt.clabel(cs, inline=True, fontsize=8, fmt='%d')
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.title("DEM with 250 m contourline interval")

plt.show()

#Extract the contour lines
lines = []
elevations = []

for level, seglist in zip(cs.levels, cs.allsegs):
    for seg in seglist:
        if len(seg) > 1:
            lines.append(LineString(seg))
            elevations.append(level)

#Create a GeoDataFrame
gdf_contours = gpd.GeoDataFrame(
    {"elevation": elevations},
    geometry=lines,
    crs="EPSG:32737"   # UTM Zone 37S
)

#Check geodataframe created
print(gdf_contours.head())
print(len(gdf_contours))

   elevation                                           geometry
0     1250.0  LINESTRING (336407.476 9938194.587, 336386.742...
1     1250.0  LINESTRING (301087.775 9931133.207, 301098.676...
2     1250.0  LINESTRING (336407.476 9937943.243, 336377.451...
3     1250.0  LINESTRING (299657.5 9924743.295, 299662.378 9...
4     1250.0  LINESTRING (334155.639 9940445.145, 334151.468...
615

#Save as shapefile
gdf_contours.to_file(work_dir+"/mount_kenya_contours.shp")

Derive class-coverage statistics

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

#selected merged catchment map upstream of outlet location
catch_sel = ilwis.RasterCoverage('catchmerge_raster.mpr')

#selected cross map - here the threshold map with 4 classes is used
raster_sel = ilwis.RasterCoverage('thres4.mpr')

#perform the crossing operation - here the variable threshold map is used
ct = ilwis.do('cross',catch_sel, raster_sel, 'ignoreundef')

ct.store('ct.tbt')

#review table size and columns names
print(ct.columnCount())
print(ct.recordCount())
print(ct.columns())

4
4
('first_raster', 'second_raster', 'pixel_count', 'pixel_area')

#review content of table
for i in range(ct.recordCount()):
    rec = ct.record(i)
    print(rec)

(0, 1500.0, 51163, 46046700.0)
(0, 10000.0, 593457, 534111300.0)
(0, 5000.0, 376328, 338695200.0)
(0, 15000.0, 711811, 640629900.0)

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

An explanation expression below:

• 'ct1' = output table
• '(@1/@2)' = calculation expression using 2 parameters, here (a/b), parameters to be defined lateron in the expression
• 'test_calc' = new output column'
• 'ct' = input table
• 'yes' or 'no' = if 'yes' option to create new table', if 'no' - existing table is used
• column pixel_count = @1
• column second_raster = @2

#table calculator in ilwis
ct1 = ilwis.do('tabcalc', '(@1/@2)', ct, 'test_calc','yes', 'pixel_count', 'second_raster')

print(ct1.columns())

('first_raster', 'second_raster', 'pixel_count', 'pixel_area', 'test_calc')

for i in range(ct1.recordCount()):
    rec = ct1.record(i)
    print(rec)

(0, 1500.0, 51163, 46046700.0, 34.108666666666664)
(0, 10000.0, 593457, 534111300.0, 59.3457)
(0, 5000.0, 376328, 338695200.0, 75.2656)
(0, 15000.0, 711811, 640629900.0, 47.45406666666667)

#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
#to get percentage per class we need the total number of pixels!
ct_stat = ilwis.PropertySets.pSUM
pc_stat = ct.statistics('pixel_count',ct_stat) # derive the sum for the column 'pixel_count'
pixel_count_sum = pc_stat.prop(ilwis.PropertySets.pSUM)
print('Sum = ', pixel_count_sum)

Sum =  1732759.0

ct2 = ilwis.do('tabcalc', '(@1/@2*100)', ct, 'area_percent','yes', 'pixel_count', pixel_count_sum)

print(ct2.columnCount())
print(ct2.recordCount())
print(ct2.columns())

5
4
('first_raster', 'second_raster', 'pixel_count', 'pixel_area', 'area_percent')

for i in range(ct2.recordCount()):
    rec = ct2.record(i)
    print(rec)

(0, 1500.0, 51163, 46046700.0, 2.9526899009037035)
(0, 10000.0, 593457, 534111300.0, 34.2492522041438)
(0, 5000.0, 376328, 338695200.0, 21.718427086513472)
(0, 15000.0, 711811, 640629900.0, 41.07963080843903)

ct2.store('ct_fin.tbt')

#ilwis cross table to pandas dataframe
data = {}
for col in list(ct2.columns()):
    data[col] = list(ct2.column(col))
df = pd.DataFrame(data)

df

	first_raster	second_raster	pixel_count	pixel_area	area_percent
0	0	1500.0	51163	46046700.0	2.952690
1	0	10000.0	593457	534111300.0	34.249252
2	0	5000.0	376328	338695200.0	21.718427
3	0	15000.0	711811	640629900.0	41.079631

print(df['pixel_count'].sum())
print()
print(df['area_percent'].sum())

1732759

100.0

#store dataframe as csv file - without index column
df.to_csv(work_dir+'/ct_fin.csv', index=False)

Derive aggregate statistics

#perform the crossing operation - here the variable dem_fill map with elevation values is used (with precision of 1 - meter)
raster_sel1 = ilwis.RasterCoverage('thres4.mpr')
raster_sel2 = ilwis.RasterCoverage('dem_fill.mpr')
raster_sel2 = ilwis.do('setvaluerange', raster_sel2, 0, 10000, 1)#set dem to interval of 1 meter - no decimals
ct = ilwis.do('cross',raster_sel1, raster_sel2, 'ignoreundef')

print()
print(ct.columnCount(), ct.recordCount())
print(ct.columns())

4 12987
('first_raster', 'second_raster', 'pixel_count', 'pixel_area')

#check the first 10 records of the cross-table
for i in range(min(10,ct.recordCount())):
    rec = ct.record(i)
    print(rec)

(5000.0, 3755.0, 110, 99000.0)
(1500.0, 3508.0, 33, 29700.0)
(5000.0, 3647.0, 147, 132300.0)
(15000.0, 1179.0, 2789, 2510100.0)
(10000.0, 3129.0, 40, 36000.0)
(10000.0, 2942.0, 97, 87300.0)
(10000.0, 3283.0, 9, 8100.0)
(5000.0, 2928.0, 302, 271800.0)
(10000.0, 2980.0, 130, 117000.0)
(5000.0, 3805.0, 115, 103500.0)

ct.store('aggr.tbt')

#ilwis cross table to pandas dataframe
data = {}
for col in list(ct.columns()):
    data[col] = list(ct.column(col))
df = pd.DataFrame(data)

df

	first_raster	second_raster	pixel_count	pixel_area
0	5000.0	3755.0	110	99000.0
1	1500.0	3508.0	33	29700.0
2	5000.0	3647.0	147	132300.0
3	15000.0	1179.0	2789	2510100.0
4	10000.0	3129.0	40	36000.0
...	...	...	...	...
12982	1500.0	1085.0	1	900.0
12983	1500.0	1119.0	1	900.0
12984	1500.0	1057.0	1	900.0
12985	1500.0	1102.0	1	900.0
12986	5000.0	1056.0	5	4500.0

12987 rows × 4 columns

Retrieve (weigthed) aggregate statisitics of the elevation for each of the threshold classess

# ----- Weighted Functions -----
def weighted_median(values, weights):
    sorted_idx = np.argsort(values)
    values_sorted = np.array(values)[sorted_idx]
    weights_sorted = np.array(weights)[sorted_idx]
    cum_weights = np.cumsum(weights_sorted)
    cutoff = 0.5 * sum(weights_sorted)
    return values_sorted[np.searchsorted(cum_weights, cutoff)]

def weighted_std(values, weights):
    average = np.average(values, weights=weights)
    variance = np.average((values - average)**2, weights=weights)
    return np.sqrt(variance)

# ----- Unweighted Aggregations -----
unweighted = (
    df.groupby('first_raster')['second_raster']
    .agg(
        min_elevation='min',
        max_elevation='max',
        std_elevation='std',
        median_elevation='median',
        sum_elevation='sum',
        count_elevation='count',
        predominant_elevation=lambda x: x.mode().iloc[0] if not x.mode().empty else np.nan
    )
    .reset_index()
)

# ----- Weighted Average -----
df['weighted_elevation'] = df['second_raster'] * df['pixel_count']
weighted_avg = (
    df.groupby('first_raster')
    .apply(lambda g: (g['second_raster'] * g['pixel_count']).sum() / g['pixel_count'].sum())
    .reset_index(name='weighted_avg_elevation')
)

# ----- Weighted Median & Std -----
def group_weighted_stats(group):
    values = group['second_raster'].values
    weights = group['pixel_count'].values
    return pd.Series({
        'weighted_median_elevation': weighted_median(values, weights),
        'weighted_std_elevation': weighted_std(values, weights)
    })

weighted_stats = df.groupby('first_raster').apply(group_weighted_stats).reset_index()

# ----- Merge All Together -----
final = (
    unweighted
    .merge(weighted_avg, on='first_raster')
    .merge(weighted_stats, on='first_raster')
)

# Display or export final DataFrame
print(final)

   first_raster  min_elevation  max_elevation  std_elevation  \
0        1500.0         1057.0         5159.0    1108.993848   
1        5000.0         1055.0         4858.0    1071.421406   
2       10000.0         1048.0         4530.0     989.809924   
3       15000.0         1048.0         4367.0     708.753414   

   median_elevation  sum_elevation  count_elevation  predominant_elevation  \
0            3146.0     12017849.0             3821                 1057.0   
1            2905.0     10757621.0             3701                 1055.0   
2            2758.5      9443224.0             3422                 1048.0   
3            2069.0      4349414.0             2043                 1048.0   

   weighted_avg_elevation  weighted_median_elevation  weighted_std_elevation  
0             3238.222600                     3164.0              839.230736  
1             2455.404833                     2359.0              759.897871  
2             1810.032307                     1678.0              561.946523  
3             1258.761122                     1166.0              249.407797  

#store dataframe as csv file - without index column
final.to_csv(work_dir+'/agg_stats.csv', index=False)

Cumulative hypsometric curve

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.

merged_catchment_map = ilwis.RasterCoverage('catchment.mpr') 
selected_catchment = 49 #use catchment-id 50  -note difference of '1' due to index offset
elevation_map = raster_sel2

c_temp = ilwis.do('mapcalc','iff(@1==' + str(selected_catchment) + ',1,?)',merged_catchment_map)

c_temp.store('sel_catch.mpr')

ct = ilwis.do('cross',c_temp, elevation_map, 'ignoreundef')#note first output is table, second is map

#ilwis cross table to pandas dataframe
data = {}
for col in list(ct.columns()):
    data[col] = list(ct.column(col))
df = pd.DataFrame(data)

#sort on elevation - from low to high elevation
df_sorted = df.sort_values(by='second_raster')
df_sorted

	first_raster	second_raster	pixel_count	pixel_area
76	1.0	1767.0	8	7200.0
1303	1.0	1768.0	3	2700.0
1302	1.0	1769.0	28	25200.0
662	1.0	1770.0	30	27000.0
214	1.0	1771.0	29	26100.0
...	...	...	...	...
680	1.0	3083.0	1	900.0
518	1.0	3084.0	3	2700.0
138	1.0	3087.0	2	1800.0
676	1.0	3089.0	1	900.0
674	1.0	3093.0	2	1800.0

1304 rows × 4 columns

total_pixels = df_sorted['pixel_count'].sum()
df_sorted['Cum_pixels'] = df_sorted['pixel_count'].cumsum()
df_sorted['Cum_area'] = df_sorted['Cum_pixels']* total_pixels
df_sorted['Percentage'] = (df_sorted['pixel_count'] / total_pixels) * 100
df_sorted['Cum_perc'] = df_sorted['Percentage'].cumsum().round(3)

total_percent = df_sorted['Percentage'].sum()

if total_percent > 100:
    # Find the index of the row with the max percentage
    max_idx = df_sorted['Cum_perc'].idxmax()
    # Reduce the max value by the excess
    df_sorted.loc[max_idx, 'Cum_perc'] -= round(total_percent - 100, 2)

#print("Total number of pixels:", total_pixels)
df_sorted

	first_raster	second_raster	pixel_count	pixel_area	Cum_pixels	Cum_area	Percentage	Cum_perc
76	1.0	1767.0	8	7200.0	8	265712	0.024086	0.024
1303	1.0	1768.0	3	2700.0	11	365354	0.009032	0.033
1302	1.0	1769.0	28	25200.0	39	1295346	0.084302	0.117
662	1.0	1770.0	30	27000.0	69	2291766	0.090323	0.208
214	1.0	1771.0	29	26100.0	98	3254972	0.087313	0.295
...	...	...	...	...	...	...	...	...
680	1.0	3083.0	1	900.0	33206	1102904084	0.003011	99.976
518	1.0	3084.0	3	2700.0	33209	1103003726	0.009032	99.985
138	1.0	3087.0	2	1800.0	33211	1103070154	0.006022	99.991
676	1.0	3089.0	1	900.0	33212	1103103368	0.003011	99.994
674	1.0	3093.0	2	1800.0	33214	1103169796	0.006022	100.000

1304 rows × 8 columns

plt.figure(figsize=(12, 5))
plt.plot(df_sorted['second_raster'], df_sorted['Cum_perc'], label = 'elevation-percentage')
plt.xlabel("Elevation")
plt.ylabel("Percentage")
plt.title("Cumulative hypsometric curve of selected catcment")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show();

#store dataframe as csv file - without index column
df_sorted.to_csv(work_dir+'/ct_hypso.csv', index=False)

Plot the longest flowpath longitudinal profile

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.

#import longest flow path vector map
longest_flow_path = ilwis.FeatureCoverage('catchmerge_long_flowpath.mps')

#convert longest flow path segment to raster
longest_flow_path_r = ilwis.do('line2raster', longest_flow_path, dem.geoReference())

ct_lfp = ilwis.do('cross',longest_flow_path_r, dem, 'ignoreundef')

#ilwis cross table to pandas dataframe
data = {}
for col in list(ct_lfp.columns()):
    data[col] = list(ct_lfp.column(col))
df = pd.DataFrame(data)

#sort on elevation - from low to high elevation
df_sorted = df.sort_values(by='second_raster')
df_sorted

	first_raster	second_raster	pixel_count	pixel_area
1592	0	1049.895142	1	900.0
406	0	1050.338257	1	900.0
2826	0	1050.743896	2	1800.0
198	0	1050.864136	2	1800.0
1344	0	1051.656738	1	900.0
...	...	...	...	...
381	0	4512.178223	1	900.0
382	0	4528.774414	1	900.0
371	0	4546.852051	1	900.0
62	0	4563.362793	1	900.0
377	0	4584.410156	1	900.0

2828 rows × 4 columns

#load the catchmerge_long_flowpath table to get length of longest flow path
table = ilwis.Table("catchmerge_long_flowpath.tbt")

recCount = table.recordCount()
colCount = table.columnCount()
columns = table.columns()  # (listing of columns in table: 'column1','column2',...)

print(recCount)
print(colCount)
print(columns)

1
6
('UpstreamCoord', 'DownstreamCoord', 'Length', 'StraightLength', 'Sinuosity', 'catchmerge_long_flowpath')

#get first record of column 'Length'
seg_length = table.cell("Length", 0)
print(seg_length)

103022.412

# Elevation values
elevation = np.sort(df['second_raster'].values)  #from low to high

# Create distance x-axis
total_length = seg_length  # meters

distance = np.linspace(0, total_length, len(elevation)) #distribute the elevation values evenly along the length of segment

distance_km = distance / 1000 #transform to km

plt.figure(figsize=(12, 5))
plt.plot(distance_km, elevation, linewidth=2, label = 'flow-path elevation')

plt.xlabel('Distance (km)')
plt.ylabel('Elevation (m)')
plt.title('Longest Flow-path Longitudinal Profile')

plt.legend()
plt.grid(True, alpha=0.3)

plt.show()

Calculate DEM derived compound indexes

#assign input maps
dem = ilwis.RasterCoverage('dem_fill.mpr')
dem_fa = ilwis.RasterCoverage('dem_fa.mpr')

#assume metric coordinate system, if lat-lon then first transform to metric / UTM
dx = ilwis.do('linearrasterfilter', dem, 'dfdx')
dy = ilwis.do('linearrasterfilter', dem, 'dfdy')

#get the pixelsize from the image dimensions
dimensions = dem.envelope()
print(dimensions)
pix_size = (dimensions.maxCorner().x - dimensions.minCorner().x) / dem.size().xsize
print(pix_size)

299657.500000 9911539.000000 336437.500000 9983419.000000
30.0

#calculate slope map in percentage
slope_percent = ilwis.do('mapcalc','100*(sqrt(pow(@1,2)+(pow(@2,2)))/'+ str(pix_size) + ')',dx,dy)

slope_percent.store('slope_percent.mpr')

#calculate slope map in degree - ilwis method
slope_degree = ilwis.do('mapcalc','180/3.141592654*atan(@1/100)',slope_percent)

slope_degree.store('slope_degree.mpr')

#calculate slope in degree - python numpy method
slope_percent_2np = np.fromiter(iter(slope_percent), np.float64, slope_percent.size().linearSize()) 
slope_percent_2np = slope_percent_2np.reshape((slope_percent.size().ysize, slope_percent.size().xsize))
slope_degrees_map = np.degrees(np.arctan(slope_percent_2np / 100))

grf = ilwis.GeoReference(slope_percent.coordinateSystem(), slope_percent.envelope(), slope_percent.size())
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 1000.0, 0.001))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(slope_percent.size().xsize, slope_percent.size().ysize, 1))
rcNew.setGeoReference(grf)
rcNew.setDataDef(dfNum)

rcNew.array2raster(slope_degrees_map)
rcNew.store('slope_degree_numpy.mpr')

#remove 0-slope values, change these to 0.1, else compound index calculation will report division by 0
slp_temp = ilwis.do('mapcalc', 'iff(@1>0,@1,0.1)',rcNew)

slp_temp.store('slope_no0.mpr')

#calculate the wetness index
wet_index = ilwis.do('mapcalc','ln((@1*' + str(pix_size * pix_size) + ') / tan(@2 * 3.141592654 / 180.0))',dem_fa,slp_temp)

wet_index.store ('wetness_index.mpr')

#calculate the stream power index
stream_index = ilwis.do('mapcalc','((@1*' + str(pix_size * pix_size) + ') * tan(@2 * 3.141592654 / 180.0))',dem_fa,slp_temp)

stream_index.store ('stream_power_index.mpr')

#calculate sediment transport index
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)

sediment_index.store ('sediment_stream_index.mpr')

Develop a SCS dimensionless Unit Hydrograph and apply this to the merged catchment and calculate the discharge for a given rainfall event

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.

The key properties are:

• Volume conservation: V=A×d where:

    A = watershed area
    d = unit rainfall excess depth

• Linearity: If rainfall excess doubles, discharge doubles.

• Superposition: Multiple rainfall pulses are handled through convolution.

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

Retrieve a small mountainous catchment

# First import the outlet point map
outlet = ilwis.FeatureCoverage('UH_location.mpp') 

#perform the catchment merge
CatchMerge_outlet = ilwis.do('MapCatchmentMergeWithOutlet','dr_net.mpr', dem_fd,dem_fa,fillsink,outlet,'no')

#store the results
#note different polygon and vector maps created
CatchMerge_outlet[0].store('UH_catchmerge_polygon.mpa') 
CatchMerge_outlet[1].store('UH_catchmerge_raster.mpr') 
CatchMerge_outlet[2].store('UH_catchmerge_long_flowpath.mps')
CatchMerge_outlet[3].store('UH_catchmerge_segments.mps')

Read the required data from merged catchment table "catchmerge_polygon.tbt":

• Catchment area (km²) = "CatchmentArea"
• Longest flow path (m) = "LongestFlowPathLength"
• Longest drainage length (m) = "LongestDrainageLength"
• Elevation divide (m) = "LDPUpstreamElevation"
• Elevation outlet (m) = "OutletElevation"

#load the catchmerge_long_flowpath table to get length of longest flow path
UH_table = ilwis.Table("UH_catchmerge_polygon.tbt")

recCount = UH_table.recordCount()
colCount = UH_table.columnCount()
columns = UH_table.columns()  # (listing of columns in table: 'column1','column2',...)

print(recCount)
print(colCount)
print(columns)

1
20
('DrainageID', 'UpstreamLinkCatchment', 'DownstreamLinkCatchment', 'DrainageLen', 'TotalDrainageLength', 'DrainageDensity', 'LongestFlowPathLength', 'LongestDrainageLength', 'CenterDrainage', 'OutletCoord', 'OutletElevation', 'LFPUpstreamCoord', 'LFPUpstreamElevation', 'LDPUpstreamCoord', 'LDPUpstreamElevation', 'CenterCatchment', 'Perimeter', 'CatchmentArea', 'TotalUpstreamArea', 'UH_catchmerge_polygon')

#get first record of required columns
A_m2 = UH_table.cell("CatchmentArea",0)
L = UH_table.cell("LongestFlowPathLength",0)
drainage_length = UH_table.cell("LongestDrainageLength",0)
elev_divide = UH_table.cell("LDPUpstreamElevation",0)
elev_outlet = UH_table.cell("OutletElevation",0)

print("Catchment area (km2) =",A_m2/1000000)
print("Longest flow path (m) =",L) 
print("Longest drainage length (m) =",drainage_length)
print("Elevation divide (m) =",elev_divide)
print("Elevation outlet (m) =",elev_outlet)

Catchment area (km2) = 129.696974999
Longest flow path (m) = 29571.961
Longest drainage length (m) = 26874.609
Elevation divide (m) = 2874.308
Elevation outlet (m) = 1549.637

H = elev_divide - elev_outlet  # elevation difference (m)
S = H / L                      # watershed slope (m/m)

print("elevation difference (m) =",H)
print("watershed slope (m/m) =",S)

elevation difference (m) = 1324.671
watershed slope (m/m) = 0.044794831157798434

# Time of Concentration (The time required for runoff from the hydraulically most distant point in the watershed to reach the outlet)
# Tc according to Kirpich (1940) - for small watersheds
L_feet = drainage_length/0.3048

#Tc_min = 0.0195 * (L ** 0.77) * (S ** (-0.385)) # SI units
Tc_min = 0.0078 * (L_feet ** 0.77) * (S ** (-0.385)) # U.S. customary form (feet)
Tc_K_hr = Tc_min / 60.0

# TC according to California Culverts Practice (1942):
dr_length_miles = (L/1000)/1.6093
#print(dr_length_miles)
elev_dif_feet = H/0.3048
#print(elev_dif_feet)

Tc_CCP_hr = ((11.9 * dr_length_miles**3) / elev_dif_feet) ** 0.385

Tc_Avg_hr = (Tc_K_hr+Tc_CCP_hr)/2


print("\nKirpich and CCP Tc approximation")
print("--------------------")
print("TC-Kirpich (hour) =",Tc_K_hr)
print("TC-CCP (hour) =", Tc_CCP_hr)
print("Tc-Average (hour) =", Tc_Avg_hr)

Kirpich and CCP Tc approximation
--------------------
TC-Kirpich (hour) = 2.761356931327822
TC-CCP (hour) = 2.9759151598128573
Tc-Average (hour) = 2.86863604557034

#SCS dimensionless Unit Hydrograph 
#for the coordinates / ratios of the dimensionless Unit Hydrograph see also: https://www.nohrsc.noaa.gov/technology/gis/uhg_manual.html 

tr = 1 #time step in hrs
Depth = 1 #rainfall depth in cm/hr
Tc = Tc_Avg_hr #time of concentration in hrs
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)
Tp = Depth / 2 + tp #time to peak in hrs
Area = A_m2/1000000 #catchment area in km2

Qp = 2.08 * Area / Tp #Peak discharde [m3/(s.cm)]

print('Lag time (hr) =', tp)
print('Time to peak (hr)=', Tp)
print('Catchment area (km2)=', Area)
print('Peak discharge (m3/sec)= ', Qp)

Lag time (hr) = 1.7211816273422038
Time to peak (hr)= 2.2211816273422036
Catchment area (km2)= 129.696974999
Peak discharge (m3/sec)=  121.45324122850681

#dimensionless UH shape
#t/TP - time ratios (x-axis)
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])

#q/Qp - discharge ratios (y-axis)
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])

# Plot
plt.figure(figsize=(6, 4))

plt.plot(t_Tp, q_Qp, marker='o', linewidth=2)

plt.xlabel('t/Tp')
plt.ylabel('Q/Qp')
plt.title('Standard SCS Dimensionless Unit Hydrograph')
plt.grid(True)

plt.xlim(0, 5)
plt.ylim(0, 1.05)

plt.tight_layout()
plt.show()

#calculate t (hrs)
t = t_Tp*Tp
print(t)

[ 0.          0.22211816  0.44423633  0.66635449  0.88847265  1.11059081
  1.33270898  1.55482714  1.7769453   1.99906346  2.22118163  2.44329979
  2.66541795  2.88753612  3.10965428  3.33177244  3.5538906   3.77600877
  3.99812693  4.22024509  4.44236325  4.88659958  5.33083591  5.77507223
  6.21930856  6.66354488  7.10778121  7.55201753  7.99625386  8.44049018
  8.88472651  9.99531732 11.10590814]

#calculate discharge [m3/(s.cm)]
discharge = q_Qp * Qp
print(discharge)

[  0.           3.64359724  12.14532412  23.07611583  37.65050478
  57.08302338  80.15913921  99.59165781 112.95151434 120.23870882
 121.45324123 120.23870882 112.95151434 104.44978746  94.73352816
  82.58820404  68.01381509  55.86849097  47.36676408  40.07956961
  34.00690754  25.14082093  17.85362646  12.99549681   9.35189957
   6.67992827   4.85812965   3.522144     2.55051807   1.82179862
   1.33598565   0.60726621   0.        ]

#transform to 1 hr time steps (note that t as provided above is not having hourly intervals)
t_hr = np.arange(0, np.ceil(t.max()) + 1, 1.0)
print(len(t))
print(t_hr)

33
[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12.]

#conduct the interpolation of discharge for the correcponding hourly timesteps
UH_1hr = np.interp(t_hr, t, discharge)
print(len(UH_1hr))
print(UH_1hr)

13
[0.00000000e+00 4.74077316e+01 1.20243830e+02 9.95302066e+01
 4.73053129e+01 2.32806153e+01 1.11506540e+01 5.30013685e+00
 2.54437294e+00 1.26034840e+00 6.04705739e-01 5.79101066e-02
 0.00000000e+00]

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

plt.plot(t_hr, UH_1hr, marker='o', linewidth=2, markersize=5, label='1-hour Unit Hydrograph')

plt.xlabel('Time (hr)')
plt.ylabel('Discharge (m³/s)')
plt.title('1-Hour SCS Unit Hydrograph of selected catchment')
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

Rainfall excess.

In unit hydrograph theory, rainfall excess is the portion of rainfall that actually becomes direct runoff at the watershed outlet.

Not all rainfall contributes immediately to streamflow. Some is lost to:

• Infiltration into the soil
• Interception by vegetation
• Depression storage (water trapped in puddles, surface irregularities)
• Evaporation and evapotranspiration

The rainfall remaining after these abstractions is called: Rainfall Excess=Total Rainfall−Losses

# Rainfall excess (hourly, for 4 consecutive hours)
rain = np.array([0.03, 1.77, 0.22, 0.02])

# conduct the convolution
Q = np.convolve(rain, UH_1hr)
time_q = np.arange(len(Q))

#show resulting table
df = pd.DataFrame({"Hour": time_q,"Discharge_m3s": Q})

print("\n=========================")
print("CONVOLVED HYDROGRAPH")
print("=========================")
print(df)

peak_idx = np.argmax(Q)
peak_q = np.max(Q)

print("\nPeak discharge:")
print("Qpeak (m3/s)=", peak_q)
print("Time (hr) =", peak_idx)

=========================
CONVOLVED HYDROGRAPH
=========================
    Hour  Discharge_m3s
0      0       0.000000
1      1       1.422232
2      2      87.519000
3      3     226.247186
4      4     204.989422
5      5     108.730344
6      6      53.938982
7      7      25.963503
8      8      12.376330
9      9       5.930394
10    10       2.914723
11    11       1.400231
12    12       0.260743
13    13       0.024834
14    14       0.001158
15    15       0.000000

Peak discharge:
Qpeak (m3/s)= 226.24718579907903
Time (hr) = 3

# plot the results

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8), sharex=False, gridspec_kw={"height_ratios": [1, 2]})


# Rainfall Hyetograph - showing rainfall versus time
rain_time = np.arange(1, len(rain) + 1)

ax1.bar(rain_time, rain, width=0.8)

ax1.invert_yaxis()

ax1.set_ylabel("Rainfall Excess")
ax1.set_title("Rainfall Excess Hyetograph")
ax1.grid(True, alpha=0.3)

# Runoff Hydrograph - showing streamflow/discharge versus time

ax2.plot(time_q, Q, linewidth=2.5, label="Convolved Hydrograph")

ax2.axvline(Tc_Avg_hr, linestyle="--", linewidth=2, label="Tc-Average (hour)")
ax2.scatter(peak_idx, peak_q, s=60, zorder=5,label=f"Peak Discharge = {peak_q:.2f} m3/s")

ax2.set_xlabel("Time (hours)")
ax2.set_ylabel("Discharge (m3/s)")
ax2.set_title("Direct Runoff Hydrograph from Unit Hydrograph Convolution")

ax2.grid(True, alpha=0.3)
ax2.legend()

plt.tight_layout()
plt.show()

#Check: estimated time of the peak discharge

#sum rainfall
P_total_cm =  np.sum(rain) #note dominant rainfall event during the second hour
print('Total excess rainfall (cm)=', f"{P_total_cm:.2f}")

#midpoint of each 1-hour interval
time_mid = np.array([0.5, 1.5, 2.5, 3.5])

#rainfall centroid time
tc = np.sum(rain * time_mid) / np.sum(rain)
print('rainfall centroid time (hr)',f"{tc:.2f}")

Tlag = 0.6 * tc
print('lag time (hr)', f"{Tlag:.2f}")

#Expected peak timing
tpeak = tc + Tlag
print('hydrograph peak expected around (hr):', f"{tpeak:.2f}")

Total excess rainfall (cm)= 2.04
rainfall centroid time (hr) 1.61
lag time (hr) 0.97
hydrograph peak expected around (hr): 2.58

#Check: Unit Hydrograph consistency check

# Unit hydrograph time step (1 hour)
dt = 3600.0  # seconds

# Expected runoff volume for 1 cm excess rainfall
rainfall_excess_m = 0.01  # 1 cm = 0.01 m
expected_volume = rainfall_excess_m * A_m2

# Volume represented by the Unit Hydrograph
UH_volume = np.sum(UH_1hr) * dt

# Compare both
percent_error = ((UH_volume - expected_volume) / expected_volume) * 100

print("\n==============================")
print("UNIT HYDROGRAPH CHECK")
print("==============================")

print("Catchment area (km2) =", A_m2/1000000)
print("1 cm runoff volume (m3) =", expected_volume)
print("Volume under Unit Hydrograph (m3) =", UH_volume)
print("Difference (m3) =", UH_volume - expected_volume)
print("Percent error (%) =", f"{percent_error:.2f}")

if abs(percent_error) < 5:
    print("\n Unit Hydrograph is consistent with a 1-cm UH.")
else:
    print("\n Check UH scaling or watershed area.")

==============================
UNIT HYDROGRAPH CHECK
==============================
Catchment area (km2) = 129.696974999
1 cm runoff volume (m3) = 1296969.74999
Volume under Unit Hydrograph (m3) = 1291268.966760582
Difference (m3) = -5700.783229417866
Percent error (%) = -0.44

 Unit Hydrograph is consistent with a 1-cm UH.

#check: Total volume check
P_total_m =  np.sum(rain)/100
print('Total excess rainfall (m) =', P_total_m)
print('Catchment area (m2) =', A_m2)

TotalVolume_basin = P_total_m * A_m2
print('Total basin volume (m3) =',TotalVolume_basin)

runoff_volume = np.sum(Q) * 3600
print('Total runoff volume (m3) =', runoff_volume)

# Compare both
Vpercent_error = ((TotalVolume_basin - runoff_volume) / runoff_volume) * 100

print("Difference (m3) =", TotalVolume_basin - expected_volume)
print("Percent error (%) =", f"{Vpercent_error:.2f}")

if abs(percent_error) < 5:
    print("\n Total runoff volume is consistent with total basin precipitation volume.")
else:
    print("\n Check  rainfall or discharge volumes used.")

Total excess rainfall (m) = 0.0204
Catchment area (m2) = 129696974.999
Total basin volume (m3) = 2645818.2899796003
Total runoff volume (m3) = 2634188.692191587
Difference (m3) = 1348848.5399896004
Percent error (%) = 0.44

 Total runoff volume is consistent with total basin precipitation volume.

# Export to csv
df.to_csv("convolved_hydrograph.csv", index=False)