Berdasarkan isu #90: kalibrasi model NRECA
Referensi isu:
hidrokit.contrib.taruma.hk89
#89. (manual/notebook). Pemodelan NRECADeskripsi permasalahan:
Strategi permasalahan:
GridSearchCV
pada scikit-learn
, yang membuat parameter grid, dan melakukan pemodelan berdasarkan seluruh kombinasi yang memungkinkan dalam parameter grid.Catatan:
nama_metrik(simulasi, prediksi)
.import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
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
try:
import HydroErr as he
except ModuleNotFoundError:
!pip install HydroErr -q
import HydroErr as he
Building wheel for HydroErr (setup.py) ... done
# 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]
Saving NRECA_sample.xlsx to NRECA_sample.xlsx
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
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 |
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
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.
calibration()
¶Fungsi meminta 7 argumen posisi dan 2 argumen opsional. Argumen ini bisa dibagi menjadi beberapa bagian:
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.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.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
.Dalam notebook ini akan menggunakan metrik yang telah disediakan pada paket HydroErr.
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.
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
}
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.
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.
func_parameter = {
'df': dataset,
'precip_col': 'PRECIP',
'pet_col': 'PET',
'AREA': 1450.6e6,
'CF': 0.6,
'C': 0.25,
'as_df': False
}
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()
.
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 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:
results.head(10)
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()
best_param = _best_parameter(results, calibration_parameter)
best_param
{'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
.
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')
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
- 20191215 - 1.0.0 - Initial
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.