Precipitation Event Week

Are satellite rainfall products useful for hydrological simulations?

Design of the exercise

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:

• -GPCC Global Precipitation Climatology Centre product for the initial calibration;
• -ERA-5 reanalysis product;
• -SM2RAIN product;
• -H23
• -H64

All the data are stored in different text files and are named considering the different rainfall products.

Import the necessary python libraries

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

Plot Africa basin

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()

Let's take a look to the rainfall and discharge time series

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')

STEP1: calibration of the MILc rainfall-runoff model (lumped version of the MISDc) over a basin in Africa using the GPCC rainfall product

Model calibration with Pyswarm optimization tool (PSO)

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

• 1) Initial conditions, fraction of W_max (0-1);
• 2) Field capacity;
• 3) Exponent of drainage;
• 4) Parameter of infiltration and drainage;
• 5) Fraction of drainage verusu interflow;
• 6) Coefficient lag-time relationship;
• 7) Parameter of potential evapotranspiration;
• 8) Runoff exponent.

#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

STEP2: running the rainfall-runoff model with GPCC-calibrated parameters considering different precipitation inputs

Model run using GPCC-calibrated parameters and input precipitation from ERA-5 rainfall dataset

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

Model run using GPCC-calibrated parameters and input precipitation from SM2RAIN rainfall dataset

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

Model run using GPCC-calibrated parameters and input precipitation from H23 rainfall dataset

name6='AFRICA_H23_2011_14'
data_input6=pd.read_csv(name6+'.txt',index_col=0,header = None, names = ['P','T','Q'],sep=',',parse_dates=True)

QobsQsim,data=MILC(name6,data_input6,PAR,Ab,fig=1)

Let's control what leads to poor performance of H23 rainfall product

name1='AFRICA_GPCC_2011_14'
data_GPCC=pd.read_csv(name1+'.txt',index_col=0,header = None, names = ['P','T','Q'],sep=',',parse_dates=True)
name2='AFRICA_H23_2011_14'
data_H23=pd.read_csv(name2+'.txt',index_col=0,header = None, names = ['P','T','Q'],sep=',',parse_dates=True)


data_GPCC_H23=pd.concat([data_GPCC['P'],data_H23['P']],axis=1)
data_GPCC_H23.columns = ['GPCC', 'H23']
data_GPCC_H23.head()

Corr=data_GPCC_H23['GPCC'].corr(data_GPCC_H23['H23'])
Bias=data_GPCC_H23['H23'].mean()/data_GPCC_H23['GPCC'].mean()*100
print('Corr = ',np.round(Corr,3),' Multiplicative Bias [%]=',np.round(Bias,2))

plt.figure(figsize=(12,8))
plt.plot(data_GPCC['P'],alpha=1,color='black')
plt.title('Observed discharge')
plt.ylabel('P [mc/s]', fontsize=14)
plt.plot(data_H23['P'],alpha=0.3,color='red')
plt.ylabel('P [mm/day]', fontsize=14)
plt.title('Corr = '+str(np.round(Corr,3))+' Multiplicative Bias [%]='+str(np.round(Bias,2)))
plt.legend(('GPCC','H23'))

Corr =  0.639  Multiplicative Bias [%]= 206.08

<matplotlib.legend.Legend at 0x7fdda0696650>

So bias is the main driver of the poor model simulation...

Step 3: we might try to correct its bias and re-run the model with the bias corrected product

alpha=data_GPCC_H23['H23'].mean()/data_GPCC_H23['GPCC'].mean() # ratio of means of the two products 
data_GPCC_H23['H23']=data_GPCC_H23['H23']*(1/alpha)

Bias=data_GPCC_H23['H23'].mean()/data_GPCC_H23['GPCC'].mean()*100

plt.figure(figsize=(12,8))
plt.plot(data_GPCC['P'],alpha=1,color='black')
plt.title('Observed discharge')
plt.ylabel('P [mc/s]', fontsize=14)
plt.plot(data_GPCC_H23['H23'],alpha=0.3,color='red')
plt.ylabel('P [mm/day]', fontsize=14)
plt.title('Corr = '+str(np.round(Corr,3))+' Multiplicative Bias [%]='+str(np.round(Bias,2)))
plt.legend(('GPCC','H23'))

data_H23['P']=data_GPCC_H23['H23'].values

QobsQsim,data=MILC(name6,data_H23,PAR,Ab,fig=1)

Correcting the bias of the rainfall is effective!! But let's see if recalibration is better...

Let's try to re-calibrate...

#The objective function and pso options have been previously defined

#Call instance of PSO with bounds argument
name=name6
d_input=data_input6 #input from H23 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 H23 input

PAR6=np.loadtxt('X_opt_' + name6+'.txt') #Calibrated parameters with GPM input
QobsQsim,data=MILC(name6,data_input6,PAR6,Ab,fig=1)

2020-12-17 07:45:45,748 - pyswarms.single.global_best - INFO - Optimize for 20 iters with {'c1': 0.5, 'c2': 0.3, 'w': 0.9}
pyswarms.single.global_best: 100%|██████████|20/20, best_cost=0.319
2020-12-17 07:45:53,027 - pyswarms.single.global_best - INFO - Optimization finished | best cost: 0.31904337737128685, best pos: [6.28556493e-01 7.39698353e+02 5.86969286e+00 3.50836962e+00
 7.93635129e-01 2.77171662e+00 1.73966406e+00 8.28948869e+00]

Results are now even better...

This is true also for SM2RAIN? Does recalibration improve the performance?

#The objective function and pso options have been previously defined

#Call instance of PSO with bounds argument
name=name3
d_input=data_input3 #input from SM2RAIN 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 SM2RAIN input
PAR3=np.loadtxt('X_opt_' + name3+'.txt') #Calibrated parameters with SM2RAIN input
QobsQsim,data=MILC(name3,data_input3,PAR3,Ab,fig=1)

2020-12-17 07:45:54,710 - pyswarms.single.global_best - INFO - Optimize for 20 iters with {'c1': 0.5, 'c2': 0.3, 'w': 0.9}
pyswarms.single.global_best: 100%|██████████|20/20, best_cost=0.354
2020-12-17 07:46:01,915 - pyswarms.single.global_best - INFO - Optimization finished | best cost: 0.3540204792434174, best pos: [7.07931723e-01 4.97739029e+02 6.23243054e+00 1.20055972e+00
 4.93248821e-01 1.16880625e+00 1.12294303e+00 1.18327885e+01]

Yes, recalibration boost the performance of the model as parameters are able to take into account of the bias of the the rainfall observations

Step 4: What about integrating SM2RAIN and H23? Let's calibrate with H64 product...

#The objective function and pso options have been previously defined

name5='AFRICA_H64_2011_14'
data_input5=pd.read_csv(name5+'.txt',index_col=0,header = None, names = ['P','T','Q'],sep=',',parse_dates=True)

#Call instance of PSO with bounds argument
name=name5
d_input=data_input5 #input from H64 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 H64 input
PAR5=np.loadtxt('X_opt_' + name5+'.txt') #Calibrated parameters with GPM input
QobsQsim,data=MILC(name5,data_input5,PAR5,Ab,fig=1)

2020-12-17 07:46:03,576 - pyswarms.single.global_best - INFO - Optimize for 20 iters with {'c1': 0.5, 'c2': 0.3, 'w': 0.9}
pyswarms.single.global_best: 100%|██████████|20/20, best_cost=0.34 
2020-12-17 07:46:10,725 - pyswarms.single.global_best - INFO - Optimization finished | best cost: 0.3396252546195634, best pos: [  0.51326055 455.65202055   6.79377241   2.17929674   0.73650463
   1.48149649   1.02340757  10.39675986]

 We obtain really good results. Let's see if we can do better than ERA5

Calibration using ERA-5 rainfall dataset

#The objective function and pso options have been previously defined

#Call instance of PSO with bounds argument
name=name2
d_input=data_input2 #input from ERA-5 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 ERA-5 input
PAR2=np.loadtxt('X_opt_' + name2+'.txt') #Calibrated parameters with ERA-5 input
QobsQsim,data=MILC(name2,data_input2,PAR2,Ab,fig=1)

2020-12-17 07:46:12,457 - pyswarms.single.global_best - INFO - Optimize for 20 iters with {'c1': 0.5, 'c2': 0.3, 'w': 0.9}
pyswarms.single.global_best: 100%|██████████|20/20, best_cost=0.321
2020-12-17 07:46:19,480 - pyswarms.single.global_best - INFO - Optimization finished | best cost: 0.320746678114346, best pos: [  0.59573721 419.72524882   6.10617761   4.07172258   0.48580042
   1.61975724   1.44246621   6.15771813]

Well, in this basin H64 is even better than ERA5. So are H-SAF rainfall products useful? YES THEY ARE!

Lab done!!!