notebook.community

Edit and run 





In [1]:



    


import numpy as np
import pandas as pd
import scipy.optimize as so
import sys
sys.path.append('../models')


import matplotlib.pyplot as plt
%matplotlib inline



    











In [2]:



    


data = pd.read_csv('../data/hbv_s_data.csv', index_col=0, parse_dates=True)



    











In [3]:



    


def calibrate(data, simulation_routine='HBV-s', objective_function='NS', method='DE'):
    '''
    Inputs:
    1. data - meteorological forcing data packed in pandas dataframe:
        'Temp' - daily temperature (Celsium degrees)
        'Prec' - daily precipitation (mm/day)
        'Evap' - daily potential evapotranspiration (mm/day)
    2.  'Qobs' - daily river runoff (mm\day)
    3. simulation routines (default='HBV')
    'HBV'
    'HBV-s'
    'GR4J'
    'GR4J-Cema-Neige'
    'SIMHYD'
    'SIMHYD-Cema-Neige'
    4. objective function for optimization:
        'NS'   - Nash-Sutcliffe model efficiency criterion
    5. optimization method:
        'DE' - differential evolution (global optimization routine)
    Output:
        list of optimal model parameters
    '''
    ### data preparation
    Qobs = data['Qobs']
    Qm = Qobs.mean()
    
    ### initialize the model used for parameters optimization
    if simulation_routine == 'HBV':
        import hbv
        model = hbv.simulation
        bnds = hbv.bounds()
        
    elif simulation_routine == 'HBV-s':
        import hbv_s
        model = hbv_s.simulation
        bnds = hbv_s.bounds()
        
    elif simulation_routine == 'GR4J':
        import gr4j
        model = gr4j.simulation
        bnds = gr4j.bounds()
        
    elif simulation_routine == 'GR4J-Cema-Neige':
        import gr4j_cemaneige
        model = gr4j_cemaneige.simulation
        bnds = gr4j_cemaneige.bounds()
        
    elif simulation_routine == 'SIMHYD':
        import simhyd
        model = simhyd.simulation
        bnds = simhyd.bounds()
        
    elif simulation_routine == 'SIMHYD-Cema-Neige':
        import simhyd_cemaneige
        model = simhyd_cemaneige.simulation
        bnds = simhyd_cemaneige.bounds()
        
    else:
        print("Incorrect simulation routine name, try one of:\
        'HBV', 'HBV-s', 'GR4J', 'GR4J-Cema-Neige', 'SIMHYD', 'SIMHYD-Cema-Neige' ")
    
    ### initialize objective function for optimization
    def obj_func_calc(params):
            # simulate hydrograph
            Qsim = model(data, params)
            # calculate objective function value
            return ((Qobs-Qsim)**2).sum()/((Qobs-Qm)**2).sum()
    
    if objective_function == 'NS':
        pass
    else:
        print("Incorrect objective function name, only 'NS' is available ")
    
    ### initialize optimization algorithm
    if method == 'DE':
        optimizer = so.differential_evolution
    else:
        print("Incorrect optimization method name, only 'DE' is available ")
    
    result = optimizer(obj_func_calc, bnds, maxiter=5, polish=False, disp=True)
    
    opt_param = result.x   
    
    return opt_param



    











In [4]:



    


calibrate(data)



    


















    






differential_evolution step 1: f(x)= 0.178606
differential_evolution step 2: f(x)= 0.178606
differential_evolution step 3: f(x)= 0.17535
differential_evolution step 4: f(x)= 0.17535
differential_evolution step 5: f(x)= 0.17535










    
Out[4]:








array([  3.95942878e+00,   2.40262977e+02,  -3.00000055e-01,
         1.82313905e+00,   8.07447193e-01,   3.58194271e-02,
         7.41046183e-02,   1.00000000e+00,   1.00000000e+00])

















In [4]:



    


calibrate(data, 'HBV')



    


















    






differential_evolution step 1: f(x)= 0.39603
differential_evolution step 2: f(x)= 0.39603
differential_evolution step 3: f(x)= 0.371808
differential_evolution step 4: f(x)= 0.360193
differential_evolution step 5: f(x)= 0.306383










    
Out[4]:








array([  1.42162095e+00,   1.99878780e-01,   2.56654524e+02,
         2.63302986e-01,   4.53993848e-02,   8.05797820e-02,
         6.43914664e-01,   2.67625928e+00,   8.59815986e-01,
         4.67813984e+02,   9.34958026e-01,  -5.88223092e-01,
         1.60717951e+00,   8.41218969e-01,   8.65562258e-02,
         3.68449293e-02])

















In [6]:



    


calibrate(data, 'GR4J-Cema-Neige')



    


















    






differential_evolution step 1: f(x)= 0.326109
differential_evolution step 2: f(x)= 0.29128
differential_evolution step 3: f(x)= 0.29128
differential_evolution step 4: f(x)= 0.29128
differential_evolution step 5: f(x)= 0.29128










    
Out[6]:








array([  5.43173141e+02,   1.00107790e-01,   1.83903799e+02,
         3.44666414e+00,   6.90110293e-01,   4.78345871e+00])

















In [7]:



    


calibrate(data, 'SIMHYD-Cema-Neige')



    


















    






differential_evolution step 1: f(x)= 0.40939
differential_evolution step 2: f(x)= 0.287975
differential_evolution step 3: f(x)= 0.287975
differential_evolution step 4: f(x)= 0.287975
differential_evolution step 5: f(x)= 0.287975










    
Out[7]:








array([  1.98586482e+01,   2.99931418e+02,   1.62807039e+00,
         7.89643878e+02,   6.67188456e-02,   4.97919324e-01,
         5.62709510e-02,   5.08991978e-01,   3.84080597e+00,
         1.58011502e-01,   5.79110124e-01,   3.24948487e+00])

















In [ ]:



    






    











In [ ]:



    






    











In [ ]:



    






    











In [ ]:



    


import pandas as pd
import sys
sys.path.append('../tools/')
from wfdei_to_lumped_dataframe import dataframe_construction
from metrics import NS
#import hbv_s



    











In [ ]:



    


data = dataframe_construction('../data/pur_scheme.csv')
#data['Qsim'] = hbv_s.simulation(data, [1, 260, -0.03, 1.3, 0.61, 0.03, 0.08, 1, 1])
obs = pd.read_csv('../data/pur_observations.csv', index_col=0, parse_dates=True,
                      squeeze=True, header=None, names=['Date', 'Qobs'])
data = pd.concat([data, obs], axis=1)
#data_for_obs = data.ix[obs.index, ['Qsim', 'Qobs']].dropna()



    











In [ ]:



    


data = data['1979':]
NS(data_for_obs['Qobs']['1982':'1985'], data_for_obs['Qsim']['1982':'1985'])
NS(data['Qobs']['1982':'1985'], data['Qsim']['1982':'1985']) 




In [ ]:



    


# now it is a data from gridded pur river
calibrate(data)



    











In [ ]:



    


"""
[  5.14239023e+00,   3.59334992e+02,   2.39564956e+00,
         9.05581073e+00,   9.73719068e-01,   5.58885084e-02,
         8.07197632e-02,   1.00000000e+00,   1.00000000e+00])
"""



    











In [ ]:



    






    











In [ ]:



    






    











In [ ]:



    






    











In [ ]:



    


# step by step
import hbv_s
model = hbv_s.simulation
bounds = hbv_s.bounds()

Qobs = data['Qobs']
Qm = Qobs.mean()
def obj_func_calc(params):
            # simulate hydrograph
            Qsim = model(data, params)
            # calculate objective function value
            return ((Qobs-Qsim)**2).sum()/((Qobs-Qm)**2).sum()

optimizer = so.differential_evolution

result = optimizer(obj_func_calc, bounds, maxiter=5, disp=True)

plt.plot(range(len(data.Qobs)), data.Qobs, 'b', 
         range(len(data.Qobs)), model(data, result.x), 'r')



    











In [ ]:

Content source: hydrogo/MORS

Similar notebooks:

notebook.community | gallery | about