Open In Colab

Berdasarkan isu #90: kalibrasi model NRECA

Referensi isu:

  • hidrokit.contrib.taruma.hk89 #89. (manual/notebook). Pemodelan NRECA
  • Roberts, W., Williams, G., Jackson, E., Nelson, E., Ames, D., 2018. Hydrostats: A Python Package for Characterizing Errors between Observed and Predicted Time Series. Hydrology 5(4) 66, doi:10.3390/hydrology5040066

Deskripsi permasalahan:

  • Mencari parameter model NRECA terbaik berdasarkan metrik yang digunakan (RMSE, NSE, MAE, dll)

Strategi permasalahan:

  • Mengikuti ide GridSearchCV pada scikit-learn, yang membuat parameter grid, dan melakukan pemodelan berdasarkan seluruh kombinasi yang memungkinkan dalam parameter grid.
  • Input dalam kalibrasi berupa fungsi model (NRECA) dan metrik. Sehingga membuat kebebasan untuk menggunakan model atau metrik yang dikembangkan sendiri.

Catatan:

  • Fungsi metrik bisa dikembangkan sendiri dengan bentuk nama_metrik(simulasi, prediksi).
  • Dalam notebook ini, metrik hidrologi menggunakan paket HydroErr yang telah dikembangkan oleh BYU-Hydroinformatics (github).

PERSIAPAN DAN DATASET

In [0]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
In [2]:
try:
    import hidrokit
except ModuleNotFoundError:
    !pip install git+https://github.com/taruma/hidrokit.git@latest -q
    import hidrokit
print(f'hidrokit version: {hidrokit.__version__}')
  Building wheel for hidrokit (setup.py) ... done
hidrokit version: 0.3.5-beta.2
In [3]:
try:
    import HydroErr as he
except ModuleNotFoundError:
    !pip install HydroErr -q
    import HydroErr as he
  Building wheel for HydroErr (setup.py) ... done
In [5]:
# LOAD DATASET
try:
    pd.read_excel('NRECA_sample.xlsx')
    dataset_path = 'NRECA_sample.xlsx'
except:
    from google.colab import files
    dataset_path = list(files.upload().keys())[0]
Upload widget is only available when the cell has been executed in the current browser session. Please rerun this cell to enable.
Saving NRECA_sample.xlsx to NRECA_sample.xlsx
In [6]:
dataset = pd.read_excel(dataset_path, header=0, index_col=0, parse_dates=True)
dataset.info()
dataset.head()
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 120 entries, 1999-01-01 to 2008-12-01
Data columns (total 3 columns):
PRECIP    120 non-null float64
PET       120 non-null float64
OBS       120 non-null float64
dtypes: float64(3)
memory usage: 3.8 KB
Out[6]:
PRECIP PET OBS
1999-01-01 507.000000 142.63 259.00
1999-02-01 374.228918 128.84 164.22
1999-03-01 211.762683 138.19 86.26
1999-04-01 219.793874 138.32 60.21
1999-05-01 132.121255 125.36 64.82

KODE

In [0]:
from itertools import product

#ref: sklearn.model_selection.ParameterGrid
def _parameter_grid(parameter):
    items = parameter.items()
    keys, values = zip(*items)
    for combination in product(*values):
        grid = dict(zip(keys, combination))
        yield grid

def _best_parameter(results, calibration_parameter):
    key = list(calibration_parameter.keys())
    return dict(zip(key, results.iloc[0][key].values))

def calibration(observed, func, calibration_parameter, func_parameter,
                metrics, met_names, met_sort, met_min=True, observed_func=None):
    
    metrics = metrics if isinstance(metrics, (list, tuple)) else [metrics]
    met_names = (met_names if isinstance(met_names, (list)) 
                           else [met_names])

    param_grid = list(_parameter_grid(calibration_parameter))
    n_param = len(param_grid)
    print('N = {}'.format(n_param))

    observed = (
        observed_func(observed) if observed_func is not None else observed
    )

    results = []

    print('PROGRESS 0 [-x--xx--x-] 100')
    print('---------> [', end='')

    for i, p in enumerate(param_grid, start=1):
        simulated = func(**p, **func_parameter)
        met_res = [m(simulated, observed) for m in metrics]
        results.append(
            list(p.values()) + met_res
        )
        if (i % (n_param // 10)) == 0:
            print('=', end='')
    
    print('] DONE')

    columns_name = list(p.keys()) + met_names

    results = (pd.DataFrame(results, columns=columns_name)
                 .sort_values(by=met_sort, ascending=met_min)
                 .reset_index(drop=True))

    return results

FUNGSI

Fungsi calibration()

Fungsi calibration() mengiterasi seluruh kombinasi parameter yang ada dan memeriksa hasilnya dengan metrik yang diberikan, kemudian menyusunnya dan mengurutkan berdasarkan metrik pilihan dan disajikan dalam bentuk pandas.DataFrame.

Fungsi ini tidak terbatas dengan fungsi model model_NRECA(), fungsi ini dirancang untuk menerima input fungsi model yang dikembangkan sendiri dengan syarat fungsi tersebut memberi keluaran dalam bentuk yang diminta oleh fungsi metrik.

Argumen calibration()

Fungsi meminta 7 argumen posisi dan 2 argumen opsional. Argumen ini bisa dibagi menjadi beberapa bagian:

Argumen fungsi

  • observed=, argumen ini meminta data observasi atau data yang akan dibandingkan hasilnya dengan nilai dari pemodelan.
  • observed_func= (default=None), argumen opsional ini membungkus nilai observed yang disediakan agar nilai observed menyesuaikan dengan jenis masukan yang diminta oleh fungsi metrik.

Argumen model

  • func=, argumen ini meminta fungsi model. Contoh: untuk NRECA maka func=model_NRECA.
  • calibration_parameter=, argumen ini meminta dictionary daftar argumen yang akan digunakan sebagai input func dan diiterasi seluruh kombinasi dari parameter yang disediakan.
  • func_parameter=, argumen ini serupa dengan calibration_parameter, akan tetapi argumen ini meminta dictionary daftar input func yang tidak akan diiterasi.

Argumen metrik

  • metrics=, argumen ini meminta list fungsi atau fungsi yang digunakan untuk mengevaluasi nilai simulasi dengan observasi. Contoh: metrik NSE metrics=he.NSE atau dalam bentuk list metrics=[he.NSE, he.rmse].
  • met_names=, argumen ini meminta list nama dari fungsi metrik.
  • met_sort=, argumen ini meminta nama metrik yang akan diurutkan.
  • met_min= (default=True), mengurutkan kolom met_sort dari kecil ke besar atau besar ke kecil jika False.

PENGGUNAAN

Dalam notebook ini akan menggunakan metrik yang telah disediakan pada paket HydroErr.

Parameter metrik

Fungsi calibration() membutuhkan argumen untuk parameter metrik berupa metrics, met_names, dan met_sort; argumen opsional met_min.

Metrik dapat dikembangkan sendiri dengan bentuk nama_metrik(simulasi, observasi). Untuk penggunaan praktis, telah tersedia paket HydroErr oleh BYU-Hydroinformatics yang mempersiapkan daftar metrik dalam hidrologi. Daftar metrik yang tersedia bisa baca disini.

Dalam notebook ini akan diberikan contoh penggunaan fungsi metrik sendiri dan menggunakan metrik dari HydroErr. Metrik yang akan digunakan antara lain: nama_metrik (metrik buatan sendiri), he.rmse root mean square error, he.r_squared coefficient of determination, he.nse Nash-Sutcliffe Efficiency, he.kge_2012 Kling-Gupta efficiency (2012).

Hasil kalibrasi akan diurutkan berdasarkan NSE. Karena performa NSE dinyatakan lebih baik jika nilainya mendekati maksimum (1), maka argumen met_min=False agar pengurutan nilai dari besar ke kecil.

Fungsi metrik hanya menerima input dalam bentuk numpy.array atau list. Sehingga keluaran dari model harus disesuaikan dengan permintaan metrik.

Isian untuk metrics dan met_names harus berupa list. Panjang (list) met_names dan metrics harus sama.

In [0]:
import HydroErr as he

def nama_metrik(simulasi, observasi):
    return np.mean(simulasi)/np.mean(observasi)

# Dalam bentuk dictionary agar lebih mudah dibaca
metrics_parameter = {
    'metrics': [nama_metrik, he.rmse, he.r_squared, he.nse, he.kge_2012],
    'met_names': ['NAMA_METRIK', 'RMSE', 'R2', 'NSE', 'KGE_2012'],
    'met_sort': 'NSE',
    'met_min': False
}

Parameter Model NRECA

Diketahui bahwa model_NRECA() membutuhkan 10 argumen dan 2 argumen opsional. Argumen tersebut dapat dibagi menjadi dua kategori yaitu calibration_parameter dan func_parameter.

calibration_parameter

Dalam model NRECA, ingin mencari tahu kombinasi parameter terbaik dari MSTOR, GSTOR, PSUB, GWF. Isian setiap key (values) harus berupa list / iterable.

In [0]:
calibration_parameter = {
    'MSTOR': [1000, 1100, 1200], # dalam bentuk list
    'GSTOR': (100, 110, 120), # dalam bentuk tuple
    'PSUB': np.linspace(0.3, 0.8, 31), # menggunakan numpy.linspace
    'GWF': np.arange(0.2, 0.91, 0.1), # menggunakan numpy.arange
}

func_parameter

Argumen model yang diperlukan oleh model_NRECA() disertakan dalam func_parameter. Argumen as_df=False memastikan hasil keluaran model berbentuk numpy.ndarray karena permintaan dari fungsi metrik.

In [0]:
func_parameter = {
    'df': dataset,
    'precip_col': 'PRECIP',
    'pet_col': 'PET',
    
    'AREA': 1450.6e6,
    'CF': 0.6,
    'C': 0.25,

    'as_df': False
}

Kalibrasi

Nilai observasi disimpan pada dataset kolom OBS. Karena fungsi metrik meminta input dalam bentuk numpy.array maka nilai observasi harus diubah sebelum disertakan ke dalam fungsi calibration().

In [11]:
from hidrokit.contrib.taruma.hk89 import model_NRECA

observed_values = dataset.loc[:, 'OBS'].values

results = calibration(observed_values,
                      model_NRECA, calibration_parameter, func_parameter,
                      **metrics_parameter)
N = 2232
PROGRESS 0 [-x--xx--x-] 100
---------> [==========] DONE

Hasil Kalibrasi

Hasil kalibrasi disimpan di results dalam bentuk pandas.DataFrame yang telah diurutkan berdasarkan nilai NSE. Berikut nilai 10 NSE terbaik beserta parameter yang dikalibrasi disertai metrik yang lain:

In [12]:
results.head(10)
Out[12]:
MSTOR GSTOR PSUB GWF NAMA_METRIK RMSE R2 NSE KGE_2012
0 1200 120 0.400000 0.2 1.047834 36.735131 0.553446 0.526348 0.701803
1 1200 120 0.416667 0.2 1.047347 36.739909 0.549800 0.526225 0.693170
2 1200 120 0.383333 0.2 1.048322 36.753441 0.556862 0.525876 0.709968
3 1200 120 0.433333 0.2 1.046860 36.767768 0.545910 0.525506 0.684108
4 1200 110 0.400000 0.2 1.047183 36.771415 0.552507 0.525412 0.701492
5 1200 110 0.416667 0.2 1.046695 36.777133 0.548837 0.525264 0.692858
6 1200 110 0.383333 0.2 1.047670 36.788763 0.555946 0.524964 0.709658
7 1200 120 0.366667 0.2 1.048809 36.794806 0.560062 0.524808 0.717623
8 1200 110 0.433333 0.2 1.046208 36.805908 0.544922 0.524521 0.683794
9 1200 100 0.400000 0.2 1.046531 36.808446 0.551555 0.524455 0.701176

Mengambil nilai parameter terbaik (baris pertama) dengan fungsi _best_parameter()

In [13]:
best_param = _best_parameter(results, calibration_parameter)
best_param
Out[13]:
{'GSTOR': 120.0, 'GWF': 0.2, 'MSTOR': 1200.0, 'PSUB': 0.4}

Memodelkan dengan parameter terbaik yang diperoleh dari hasil kalibrasi. Karena pada func_parameter argumen as_df=False maka dibuat dictionary baru yang memberikan argumen as_df=True.

In [14]:
func_parameter_df = {
    'df': dataset,
    'precip_col': 'PRECIP',
    'pet_col': 'PET',
    
    'AREA': 1450.6e6,
    'CF': 0.6,
    'C': 0.25,

    'as_df': True
}

# results with best parameter
model_NRECA(**func_parameter_df, 
            **best_param, report='full')
Out[14]:
DAYS PRECIP PET STORAGE STORAT PRERAT ETRAT AET WATBAL EXMRAT DELSTOR GWRECH GWSTOR1 GWSTOR2 GWFLOW DFLOW FLOW DISCHARGE
1999-01-01 31.0 507.000000 142.63 1200.000000 1.587680 3.554652 1.000000 85.578000 421.422000 0.914996 35.822537 154.239785 120.000000 274.239785 54.847957 231.359678 286.207635 155.007764
1999-02-01 28.0 374.228918 128.84 1235.822537 1.635075 2.904602 1.000000 77.304000 296.924918 0.933415 19.770735 110.861673 219.391828 330.253501 66.050700 166.292510 232.343210 139.317568
1999-03-01 31.0 211.762683 138.19 1255.593272 1.661233 1.532402 1.000000 82.914000 128.848683 0.942619 7.393517 48.582067 264.202801 312.784868 62.556974 72.873100 135.430074 73.347844
1999-04-01 30.0 219.793874 138.32 1262.986789 1.671016 1.589025 1.000000 82.992000 136.801874 0.945885 7.403086 51.759515 250.227894 301.987409 60.397482 77.639273 138.036755 77.251588
1999-05-01 31.0 132.121255 125.36 1270.389874 1.680810 1.053935 1.000000 75.216000 56.905255 0.949059 2.898812 21.602577 241.589928 263.192505 52.638501 32.403866 85.042367 46.058265
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2008-08-01 31.0 117.425578 137.29 1192.596156 1.577884 0.855310 0.969462 79.858473 37.567105 0.910909 3.346888 13.688086 142.509338 156.197425 31.239485 20.532130 51.771615 28.039092
2008-09-01 30.0 66.624410 161.50 1195.943044 1.582312 0.412535 0.877312 85.011488 -18.387078 0.000000 -18.387078 -0.000000 124.957940 124.957940 24.991588 0.000000 24.991588 13.986419
2008-10-01 31.0 265.947354 166.58 1177.555966 1.557985 1.596514 1.000000 99.948000 165.999354 0.902311 16.216257 59.913239 99.966352 159.879591 31.975918 89.869858 121.845776 65.990697
2008-11-01 30.0 249.151252 150.56 1193.772224 1.579440 1.654830 1.000000 90.336000 158.815252 0.911565 14.044878 57.908150 127.903672 185.811822 37.162364 86.862225 124.024589 69.409749
2008-12-01 31.0 199.088447 144.15 1207.817102 1.598022 1.381120 1.000000 86.490000 112.598447 0.919207 9.097166 41.400512 148.649458 190.049970 38.009994 62.100769 100.110763 54.219188

120 rows × 18 columns

Changelog

- 20191215 - 1.0.0 - Initial

Copyright © 2019 Taruma Sakti Megariansyah

Source code in this notebook is licensed under a MIT License. Data in this notebook is licensed under a Creative Common Attribution 4.0 International.