Open In Colab

Berdasarkan isu #89: pemodelan NRECA

Referensi isu:

  • Crawford, N. H., & Thurin, S. M. (1981). Hydrologic estimates for small hydroelectric projects. 1800 Massachusetts Avenue NW, Washington, DC 20036: National Rural Electric Cooperative Association (NRECA).
  • Megariansyah, Taruma S. (2015): Kajian Penerapan Model NRECA di Bendung Pamarayan, Skripsi Program Sarjana, Universitas Katolik Parahyangan.

Deskripsi Permasalahan:

  • Memperoleh nilai debit dari model NRECA dengan parameter seperti Moist Storage MSTOR ($mm$), Groundwater Storage GSTOR ($mm$), PSUB, GWF, Crop Factor CF, C, Area AREA ($m^2$),

Strategi Penyelesaian:

  • Perhitungan dan penamaan peubah mengikuti makalah Crawford & Thurin (1981).
  • Membuat fungsi untuk setiap perhitungan yang digunakan dalam model (dengan ditandai imbuhan _ dalam penamaan fungsi).
  • Membuat fungsi model_NRECA() dengan kejelasan parameter sehingga bisa dikembangkan lebih lanjut untuk kalibrasi model atau mencari parameter terbaik.

Catatan:

  • Dataset yang digunakan diambil dari Megariansyah (2015).
  • Input data berupa observasi curah hujan PREP ($mm$) dan evapotranspirasi PET ($mm$) setiap bulannya.

PERSIAPAN DAN DATASET

In [0]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
In [0]:
from google.colab import files
try:
    pd.read_excel('NRECA_sample.xlsx')
except:
    dataset_path = list(files.upload().keys())[0]
In [19]:
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 2 columns):
PRECIP    120 non-null float64
PET       120 non-null float64
dtypes: float64(2)
memory usage: 2.8 KB
Out[19]:
PRECIP PET
1999-01-01 507.000000 142.63
1999-02-01 374.228918 128.84
1999-03-01 211.762683 138.19
1999-04-01 219.793874 138.32
1999-05-01 132.121255 125.36

KODE

In [0]:
# FUNGSI NRECA
def _STORAT(STORAGE, NOMINAL):
    return STORAGE/NOMINAL

def _STORAGE(STORAGE, DELSTOR):
    return STORAGE + DELSTOR

def _PRERAT(PRECIP, PET):
    return PRECIP/PET

def _ETRAT(STORAT, PRERAT):
    VAL = (STORAT / 2) + ((1 - (STORAT / 2)) * PRERAT) 
    return VAL if VAL < 1 else 1

def _AET(ETRAT, PET, CF=1):
    return ETRAT * PET * CF

def _WATBAL(PRECIP, AET):
    return PRECIP - AET

def _EXMRAT(WATBAL, STORAT):
    if WATBAL < 0:
        return 0
    if STORAT > 1:
        return 1 - (0.5 * (2 - STORAT)**2)
    else:
        return 0.5 * (STORAT**2)

def _EXMST(EXMRAT, WATBAL):
    return EXMRAT * WATBAL

def _DELSTOR(WATBAL, EXMST):
    return WATBAL - EXMST

def _GWRECH(PSUB, EXMST):
    return PSUB * EXMST

def _GWSTOR2(GWSTOR_1, GWRECH):
    return GWSTOR_1 + GWRECH

def _GWSTOR1(GWSTOR_2, GWFLOW):
    return GWSTOR_2 - GWFLOW

def _GWFLOW(GWRAT, GWSTOR_2):
    return GWRAT * GWSTOR_2

def _DFLOW(EXMST, GWRECH):
    return EXMST - GWRECH

def _FLOW(GWFLOW, DFLOW):
    return GWFLOW + DFLOW

def _DISCHARGE(FLOW, AREA, DAYS):
    return (FLOW/1000) * AREA / (DAYS * 24 * 60 * 60)

def _NOMINAL(C, PRECIP_MEAN_ANNUAL):
    return 100 + C*PRECIP_MEAN_ANNUAL
In [0]:
def model_NRECA(df, precip_col, pet_col,
                MSTOR, GSTOR, PSUB, GWF, CF, C, AREA,
                as_df=True, report='discharge'):

    # sub_df
    data = df.loc[:, [precip_col, pet_col]]

    # info df
    nrows = data.shape[0]

    # initialization
    storage, storat, prerat = (np.zeros(nrows) for _ in range(3))
    etrat, aet, watbal = (np.zeros(nrows) for _ in range(3))
    exmst, exmrat, delstor, gwrech = (np.zeros(nrows) for _ in range(4))
    gwstor1, gwstor2, gwflow = (np.zeros(nrows) for _ in range(3))
    dflow, flow, discharge = (np.zeros(nrows) for _ in range(3))

    # calculation
    precip_mean_annual = (data[precip_col].groupby(by=data.index.year)
                                          .sum()
                                          .mean())
    nominal = _NOMINAL(C, precip_mean_annual)

    days = data.index.days_in_month
    precip = data.iloc[:, 0].values
    pet = data.iloc[:, 1].values

    for i in range(nrows):

        if i != 0:
            storage[i] = _STORAGE(storage[i - 1], delstor[i - 1])
        else:
            storage[i] = _STORAGE(MSTOR, 0)

        storat[i] = _STORAT(storage[i], nominal)
        prerat[i] = _PRERAT(precip[i], pet[i])
        etrat[i] = _ETRAT(storat[i], prerat[i])
        aet[i] = _AET(etrat[i], pet[i], CF=CF)
        watbal[i] = _WATBAL(precip[i], aet[i])
        exmrat[i] = _EXMRAT(watbal[i], storat[i])
        exmst[i] = _EXMST(exmrat[i], watbal[i])
        delstor[i] = _DELSTOR(watbal[i], exmst[i])
        gwrech[i] = _GWRECH(PSUB, exmst[i])

        if i != 0:
            gwstor1[i] = _GWSTOR1(gwstor2[i - 1], gwflow[i - 1])
        else:
            gwstor1[i] = _GWSTOR1(GSTOR, 0)

        gwstor2[i] = _GWSTOR2(gwstor1[i], gwrech[i])
        gwflow[i] = _GWFLOW(GWF, gwstor2[i])
        dflow[i] = _DFLOW(exmst[i], gwrech[i])
        flow[i] = _FLOW(gwflow[i], dflow[i])
        discharge[i] = _DISCHARGE(flow[i], AREA, days[i])

    # results
    if report.lower() == 'full':
        results = np.stack((
            days, precip, pet, storage, storat, prerat, etrat, aet, watbal,
            exmrat, delstor, gwrech, gwstor1, gwstor2, gwflow, dflow, flow,
            discharge
        ), axis=1)
        columns_name = [
            'DAYS', 'PRECIP', 'PET', 'STORAGE', 'STORAT', 'PRERAT', 'ETRAT',
            'AET', 'WATBAL', 'EXMRAT', 'DELSTOR', 'GWRECH', 'GWSTOR1',
            'GWSTOR2', 'GWFLOW', 'DFLOW', 'FLOW', 'DISCHARGE'
        ]
    elif report.lower() == 'partial':
        results = np.stack((
            precip, pet, storage, gwstor2, flow, discharge), axis=1)
        columns_name = ['PRECIP', 'PET', 'STORAGE',
                        'GWSTOR2', 'FLOW', 'DISCHARGE']
    elif report.lower() == 'flow':
        results = flow
        columns_name = ['FLOW']
    elif report.lower() == 'discharge':
        results = discharge
        columns_name = ['DISCHARGE']
    else:
        raise ValueError(
            str(report) + ' not identified. ' +
            'Use full / partial / flow / discharge.'
        )

    if as_df:
        return pd.DataFrame(
            data=results, index=data.index, columns=columns_name
        )
    else:
        return results

FUNGSI

Fungsi model_NRECA()

Fungsi model_NRECA() membangkitkan nilai debit menggunakan model yang dikembangkan oleh Crawford & Thurin (1981). Untuk mengetahui masing-masing batasan parameter yang digunakan dalam model ini, harap melihat pada makalah terkait. Dalam notebook ini akan lebih fokus dalam penggunaan fungsi.

Fungsi memiliki 10 argumen yang harus diisi, dan 2 argumen opsional. Dibagi menjadi 3 bagian argumen yaitu:

  • Argumen untuk dataset:
    • df: dataset dalam bentuk pandas.DataFrame.
    • precip_col: nama kolom untuk curah hujan (presipitasi).
    • pet_col: nama kolom untuk evapotranspirasi potensial.
  • Argumen untuk parameter model:
    • MSTOR: Moist Storage ($mm$).
    • GSTOR: Groundwater Storage ($mm$).
    • PSUB: Percent Discharge to Subbase. ($0.3-0.8$)
    • GWF: Groundwater Fraction. ($0.2-0.9$)
    • CF: Crop Factor.
    • C: Koefisien karakteristik DAS. ($\{0.2,0.25\}$)
    • AREA: Luas DAS. ($m^2$)
  • Argumen opsional untuk fungsi:
    • as_df: True (default), keluaran berupa pandas.DataFrame. Keluaran berupa numpy.ndarray jika False.
    • report: discharge (default). Terdapat beberapa nilai untuk report:
      • full: keluaran akan menyertakan seluruh peubah yang dihitung dalam model seperti WATBAL, ETRAT, dll.
      • partial: keluaran hanya menyertakan kolom PRECIP, PET, STORAGE, GWSTOR2, FLOW, dan DISCHARGE.
      • flow: keluaran hanya menyertakan kolom FLOW.
      • discharge: keluaran hanya menyertakan kolom DISCHARGE.

Berikut informasi satuan keluaran model:

  • PRECIP, PET, STORAGE, AET, WATBAL, EXMST, DELSTOR, GWRECH, GWSTOR1, GWSTOR2, GWFLOW, DFLOW, FLOW bersatuan $mm/bulan$.
  • STORAGE, STORAT, PRERAT, ETRAT, EXMRAT merupakan rasio.

as_df=True, report='discharge' (default)

Keluaran berupa pandas.DataFrame.

In [22]:
model_NRECA(df=dataset, precip_col='PRECIP', pet_col='PET',
            MSTOR=1000, GSTOR=100, PSUB=0.4, GWF=0.2, CF=0.6, C=0.25,
            AREA=1450.6e6).head()
Out[22]:
DISCHARGE
1999-01-01 130.474329
1999-02-01 124.872153
1999-03-01 66.438083
1999-04-01 70.787505
1999-05-01 41.998656

Sebagian argumen bisa disimpan dalam bentuk dictionary.

In [23]:
parameter = {
    'MSTOR': 1000,
    'GSTOR': 100,
    'PSUB': 0.4,
    'GWF':  0.2,
    'CF': 0.6,
    'C': 0.25,
    'AREA': 1450.6e6
}

model_NRECA(df=dataset, precip_col='PRECIP', pet_col='PET', **parameter).head()
Out[23]:
DISCHARGE
1999-01-01 130.474329
1999-02-01 124.872153
1999-03-01 66.438083
1999-04-01 70.787505
1999-05-01 41.998656

as_df=False

Keluaran berupa numpy.array

In [24]:
model_NRECA(df=dataset, precip_col='PRECIP', pet_col='PET', **parameter,
            as_df=False)
Out[24]:
array([130.47432853, 124.87215318,  66.43808331,  70.78750455,
        41.99865614,  28.45640939,  24.28274088,  19.7163445 ,
        12.90107429,  44.81356496,  53.93315909,  90.06341788,
        95.31267461, 145.23779804, 117.57476262,  78.75685187,
        91.24644372,  57.71641073,  37.27870151,  22.07183021,
        24.97489079,  22.58939002,  35.79315297,  26.23750108,
       138.0338475 , 235.17042934, 116.11823068,  87.70033456,
       104.97934937, 109.51542424, 116.47710193,  34.51540744,
        71.42737681,  33.81058852,  34.87405532,  36.58594762,
       110.40885462, 112.85425887,  66.53789581,  97.93354457,
        35.04151568,  22.39154574,  47.55438496,  16.65090467,
        13.76474786,  10.65657899, 119.65349912,  57.66475933,
       137.39388673, 199.12246399, 159.01388739, 172.29628498,
       133.74151565,  44.77420904,  34.66390377,  27.73112302,
        36.38175359, 120.31395471,  77.61700498, 126.38336198,
       114.21457768, 123.66605913, 174.79252429, 195.67904295,
        99.25002079,  61.22475824,  71.16772612,  31.42972227,
        57.45303916,  44.65959205,  50.75185708,  91.32761479,
       141.91241982, 103.15499701,  96.64076894,  63.34595895,
        91.80413108, 133.77947175,  81.64626021,  31.2233851 ,
        56.84185425,  53.0918533 , 111.44781119, 124.64417381,
       133.68443001,  82.44572213, 103.95773659,  78.36551337,
        70.09140577,  28.12253459,  21.77228484,  17.41782788,
        14.39873771,  11.14740984,  26.2002207 ,  71.83216778,
        77.60770993,  94.77473292, 112.44031817,  75.38748734,
        52.32065854,  67.14854195,  35.46514257,  18.59214774,
        15.3695088 ,  21.20863532,  10.74189578,  49.26924462,
        91.25050052, 170.94322505,  77.98879732,  81.89122336,
        59.30935231,  27.5965865 ,  19.29362528,  28.03710337,
        13.98514355,  65.98636155,  69.40576331,  54.21634462])

report

report='full'

In [25]:
model_NRECA(df=dataset, precip_col='PRECIP', pet_col='PET', **parameter,
            report='full').head()
Out[25]:
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 1000.000000 1.323067 3.554652 1.0 85.578 421.422000 0.770881 96.555980 129.946408 100.000000 229.946408 45.989282 194.919612 240.908894 130.474329
1999-02-01 28.0 374.228918 128.84 1096.555980 1.450817 2.904602 1.0 77.304 296.924918 0.849199 44.776648 100.859308 183.957127 284.816435 56.963287 151.288962 208.252249 124.872153
1999-03-01 31.0 211.762683 138.19 1141.332627 1.510059 1.532402 1.0 82.914 128.848683 0.879979 15.464559 45.353650 227.853148 273.206798 54.641360 68.030474 122.671834 66.438083
1999-04-01 30.0 219.793874 138.32 1156.797186 1.530520 1.589025 1.0 82.992 136.801874 0.889794 15.076374 48.690200 218.565438 267.255638 53.451128 73.035300 126.486428 70.787505
1999-05-01 31.0 132.121255 125.36 1171.873560 1.550467 1.053935 1.0 75.216 56.905255 0.898960 5.749712 20.462217 213.804510 234.266727 46.853345 30.693325 77.546671 41.998656

report='partial'

In [26]:
model_NRECA(df=dataset, precip_col='PRECIP', pet_col='PET', **parameter,
            report='partial').head()
Out[26]:
PRECIP PET STORAGE GWSTOR2 FLOW DISCHARGE
1999-01-01 507.000000 142.63 1000.000000 229.946408 240.908894 130.474329
1999-02-01 374.228918 128.84 1096.555980 284.816435 208.252249 124.872153
1999-03-01 211.762683 138.19 1141.332627 273.206798 122.671834 66.438083
1999-04-01 219.793874 138.32 1156.797186 267.255638 126.486428 70.787505
1999-05-01 132.121255 125.36 1171.873560 234.266727 77.546671 41.998656

report='flow'

In [27]:
model_NRECA(df=dataset, precip_col='PRECIP', pet_col='PET', **parameter,
            report='flow').head()
Out[27]:
FLOW
1999-01-01 240.908894
1999-02-01 208.252249
1999-03-01 122.671834
1999-04-01 126.486428
1999-05-01 77.546671

report='discharge' (default)

In [28]:
model_NRECA(df=dataset, precip_col='PRECIP', pet_col='PET', **parameter,
            report='discharge').head()
Out[28]:
DISCHARGE
1999-01-01 130.474329
1999-02-01 124.872153
1999-03-01 66.438083
1999-04-01 70.787505
1999-05-01 41.998656

Changelog

- 20191214 - 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.