The exercise is divided in four different steps:
1) We will calibrate the MILc rainfall-runoff model (lumped version of the MISDc) over a basin in South Africa using the GPCC rainfall product. Being based on rain gauge data, we will assume this product the most accurate one. For the calibration we will use the Particle Swarm Optimization algorithm, optimizing the Kling-Gupta Efficiency index (KGE).
2) We will run the MILc over the basin using the calibrated parameters from GPCC and different rainfall products as inputs: ERA5, SM2RAIN, H23;
3) We will show the potential of improving flood simuation by using simple bias correction and recalibration of H23
4) We will show the potential of improving flood simuation by using an integrated product: H64
In particular, rainfall products used will be:
All the data are stored in different text files and are named considering the different rainfall products.
from MILc_2 import * # hydrological model library Brocca et al. (2011)
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import pyswarms as ps
from pyswarms.utils.plotters import plot_cost_history
map = Basemap(llcrnrlon=-20.,llcrnrlat=-35.,urcrnrlon=60.,urcrnrlat=10.,\
rsphere=(6378137.00,6356752.3142),\
resolution='l',projection='merc')
map.drawcoastlines()
map.drawcountries()
map.fillcontinents(color = 'coral')
map.drawmapboundary()
# draw parallels
map.drawparallels(np.arange(-40,10,20),labels=[1,1,0,1])
# draw meridians
map.drawmeridians(np.arange(-20,60,30),labels=[1,1,0,1])
#lonlat of the basin
lon = 30.3318
lat = -29.0799
x,y = map(lon, lat)
map.plot(x, y, 'bo', markersize=10)
label = 'Basin'
plt.text(x, y, label,color='white',
fontsize="large", weight='heavy',
horizontalalignment='right',
verticalalignment='bottom')
plt.show()
name1='AFRICA_GPCC_2011_14'
data_input1=pd.read_csv(name1+'.txt',index_col=0,header = None, names = ['P','T','Q'],sep=',',parse_dates=True)
Ab=340 #area of the basin
data_input1.head()
P | T | Q | |
---|---|---|---|
2011-01-01 | 1.7496 | 20.32 | 3.712 |
2011-01-02 | 8.1949 | 22.47 | 5.328 |
2011-01-03 | 10.8695 | 17.85 | 5.788 |
2011-01-04 | 14.7804 | 17.37 | 6.166 |
2011-01-05 | 8.4735 | 18.07 | 8.940 |
plt.figure(figsize=(12,8))
plt.subplot(211)
plt.fill_between(data_input1.index,data_input1['Q'].values,0,color='gray',alpha=0.8)
plt.title('Observed discharge')
plt.ylabel('Q [mc/s]', fontsize=14)
plt.subplot(212)
plt.plot(data_input1['P'],alpha=0.7,color='black')
plt.ylabel('P [mm/day]', fontsize=14)
plt.title('GPCC Rainfall')
Text(0.5, 1.0, 'GPCC Rainfall')
see the documentation and installation at https://pypi.org/project/pyswarms/0.1.9/
PSO https://link.springer.com/article/10.1007/s00500-016-2474-6
We are going to calibrate 8 model parameters: https://onlinelibrary.wiley.com/doi/abs/10.1002/hyp.8042
#Define objective function
def func(PARv):
#global d_input, Ab
n_particles = PARv.shape[0]
err = np.zeros(n_particles)
for i in range(n_particles):
KGE,data=MILC(name,d_input,PARv[i],Ab,Wobs=[],K=0,fig=0)
err[i] = 1 - KGE
return err
#PSO needs defined parameters boundaries.
bnds1 = (np.array([0, 200, 1, 0,
0, 0.5, 1, 1]),
np.array([1, 1000, 10, 5,
1, 3, 2, 15]))
#PSO needs defined options
options = {'c1': 0.5, 'c2': 0.3, 'w': 0.9} #https://link.springer.com/article/10.1007/s00500-016-2474-6
#Call instance of PSO with bounds argument
name=name1 #variables renamed to be used in the objective function
d_input=data_input1 #input from GPCC becomes the input for the objective function
optimizer = ps.single.GlobalBestPSO(n_particles=20, dimensions=8, options=options, bounds=bnds1)
cost, pos = optimizer.optimize(func, 20)
#Obtain cost history from optimizer instance
plot_cost_history(cost_history=optimizer.cost_history)
plt.show()
PARn = pos
np.savetxt('X_opt_' + name+'.txt', PARn) #save parameters calibrated for GPCC input
PAR=np.loadtxt('X_opt_' + name+'.txt') #Calibrated parameters with GPCC input
QobsQsim,data=MILC(name1,data_input1,PAR,Ab,fig=1)
2020-12-17 07:45:29,263 - pyswarms.single.global_best - INFO - Optimize for 20 iters with {'c1': 0.5, 'c2': 0.3, 'w': 0.9} pyswarms.single.global_best: 0%| |0/202020-12-17 07:45:29,300 - numexpr.utils - INFO - Note: NumExpr detected 16 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8. 2020-12-17 07:45:29,301 - numexpr.utils - INFO - NumExpr defaulting to 8 threads. pyswarms.single.global_best: 100%|██████████|20/20, best_cost=0.314 2020-12-17 07:45:36,986 - pyswarms.single.global_best - INFO - Optimization finished | best cost: 0.31436389197148007, best pos: [3.91652942e-01 8.06653247e+02 7.51379437e+00 4.12671604e+00 2.79176692e-01 1.88613910e+00 1.06276983e+00 4.88578284e+00]
The figure shows, in the upper plot, relative saturation and rainfall. In the bottom plot the temporal comparison between observed and simulated river discharge is shown. The scores are displayed in the title. We obtain a satisfactory KGE of 0.73 with the calibrated parameters considering the entire calibration period
name2='AFRICA_ERA5_2011_14'
data_input2=pd.read_csv(name2+'.txt',index_col=0,header = None, names = ['P','T','Q'],sep=',',parse_dates=True)
QobsQsim,data=MILC(name2,data_input2,PAR,Ab,fig=1) #Parameters have been calibrated using GPCC
name3='AFRICA_SM2R_2011_14'
data_input3=pd.read_csv(name3+'.txt',index_col=0,header = None, names = ['P','T','Q'],sep=',',parse_dates=True)
QobsQsim,data=MILC(name3,data_input3,PAR,Ab,fig=1) #Parameters calibrated using GPCC