Open In Colab

Berdasarkan isu #106: menghitung evapotranspirasi

Referensi Isu:

  • Limantara, Lily M. (2018): Rekayasa Hidrologi, Edisi Revisi. Penerbit Andi Offset, Yogyakarta. (hal. 31-47).

Deskripsi Isu:

  • Menghitung Evapotranspirasi dengan tiga rumus yaitu Rumus Blaney-Criddle, Rumus Radiasi, Rumus Penman.

Strategi:

  • Urutan pengembangan dari Blaney-Criddle, Radiasi, kemudian diakhiri Penman. Beberapa tabel yang digunakan pada Penman tersedia pada rumus sebelumnya.
  • Buat Tabel lampiran yang digunakan sebagai acuan. Tabel berupa:
    • Tabel rel_P_LL, hubungan $P$ dan Letak Lintang (LL) (Untuk Indonesia $5^{\circ}$ s/d $10^{\circ} \text{LS}$).
    • Tabel cor_C_BC, angka koreksi $C$ menurut Blaney Criddle.
    • Tabel rel_T_W, hubungan $T$ dan $W$ (untuk Indonesia, El. $0-500\ \text{m}$)
    • Tabel cor_C_RAD, angka koreksi $C$ menurut Rumus Radiasi.
    • Tabel val_Rg, harga $R_\gamma$ (Untuk Indonesia $5^{\circ}\ \text{LU}$ s/d $10^{\circ}\ \text{LS}$).
    • Tabel cor_C_PEN, angka koreksi $C$ menurut Penman.
    • Tabel rel_T_PEN, hubungan $T$ dengan $\varepsilon_\gamma, W, F(T)$.
  • Pengembangan fungsi mengikuti hk89.model_NRECA().

Catatan:

  • Sampel menggunakan soal latihan yang tersedia dalam buku.

PERSIAPAN DAN DATASET

In [0]:
import pandas as pd
import numpy as np

TABEL

Nilai tabel mengikuti referensi buku. Tabel yang digunakan dalam perhitungan dibangkitkan dengan kode dibawah ini:

In [0]:
t_rel_P_LL = pd.DataFrame({
    '5 U':      [0.27]*3 + [0.28]*6 + [0.27]*3,
    '2.5 U':    [0.27]*3 + [0.28]*6 + [0.27]*3,
    '0':        [.27]*12,
    '2.5 S':    [.28]*12,
    '5 S':      [.28]*12,
    '7.5 S':    [.28]*4 + [.27]*4 + [.28]*3 + [.29]*1,
    '10 S':     [.29] + [.28]*2 + [.27] + [.26]*4 + [.27] +[.28]*2 + [.29]
}, index=range(1, 13))

t_cor_C_BC = pd.DataFrame({
    'C':    [.8, .8, .75, .7, .7, .7, .7, .75, .8, .8, .8, .8]
}, index=range(1, 13))

t_rel_T_W = pd.DataFrame({
    'suhu': np.arange(24.0, 30, 0.2),
    'W': np.arange(0.735, 0.794, .002)
})

t_cor_C_RAD = pd.DataFrame(
    {'C': [.8]*2 + [.75]*5 + [.8]*5},
    index=range(1, 13))

_data = [13.0, 14.3, 14.7, 15.0, 15.3, 15.5, 15.8, 16.1, 16.1, 
         14.0, 15.0, 15.3, 15.5, 15.7, 15.8, 16.0, 16.1, 16.0, 
         15.0, 15.5, 15.6, 15.7, 15.7, 15.6, 15.6, 15.1, 15.3, 
         15.1, 15.5, 15.3, 15.3, 15.1, 14.9, 14.7, 14.1, 14.0, 
         15.3, 14.9, 14.6, 14.4, 14.1, 13.8, 13.4, 13.1, 12.6, 
         15.0, 14.4, 14.2, 13.9, 13.9, 13.2, 12.8, 12.4, 12.6, 
         15.1, 14.6, 14.3, 14.1, 14.1, 13.4, 13.1, 12.7, 11.8, 
         15.3, 15.1, 14.9, 14.8, 14.8, 14.3, 14.0, 13.7, 12.2, 
         15.1, 15.3, 15.3, 15.3, 15.3, 15.1, 15.0, 14.9, 13.1, 
         15.7, 15.1, 15.3, 15.4, 15.4, 15.6, 15.7, 15.8, 14.6, 
         14.8, 14.5, 14.8, 15.1, 15.1, 15.5, 15.8, 16.0, 15.6, 
         14.6, 14.1, 14.4, 14.8, 14.8, 15.4, 15.7, 16.0, 16.0]
_data = np.array(_data).reshape((12, 9)).T
t_val_Rg = pd.DataFrame(_data,
                      columns=range(1, 13),
                      index=['5 LU', '4 LU', '2 LU', '0', 
                             '2 LS', '4 LS', '6 LS', '8 LS', 
                             '10 LS']).T

t_cor_C_PEN = pd.DataFrame({
    'C': [1.1]*3 + [.9]*4 + [1.1]*5,
}, index=range(1, 13))

t_rel_T_PEN = pd.DataFrame({
    'suhu': np.arange(24., 29.1, .2),
    'e_mbar': [29.85, 30.21, 30.57, 30.94, 31.31, 31.69, 32.06, 32.45, 32.83, 
               32.22, 33.62, 34.02, 34.42, 34.83, 35.25, 35.66, 36.09, 36.50, 
               36.94, 37.37, 37.81, 38.25, 38.70, 39.14, 39.61, 40.06],
    'w': np.arange(0.735, 0.786, .002),
    'f_t': [15.40, 15.45, 15.50, 15.55, 15.60, 15.65, 15.70, 15.75, 15.80, 
            15.85, 15.90, 15.94, 15.98, 16.02, 16.06, 16.10, 16.14, 16.18, 
            16.22, 16.26, 16.30, 16.34, 16.38, 16.42, 16.46, 16.5]
})

DIAGRAM

Diagram ini dibuat untuk memberikan gambaran proses perhitungan evapotranspirasi pada modul ini. Diagram dibuat oleh taruma menggunakan draw.io dengan lisensi CC-BY-SA-4.0. Diagram mengacu pada buku Rekayasa Hidrologi: Edisi Revisi (2018) oleh Prof. Dr. Ir. Lily Montarcih Limantara, M.Sc.

Diagram Rumus Blaney Criddle dan Rumus Radiasi

Dua diagram ini menggambarkan proses perhitungan pada fungsi ETo_BlaneyCriddle() dan ETo_Radiation.

Blaney Criddle and Radiation Diagram

Diagram Rumus Penman

Diagram ini menggambarkan proses perhitungan pada fungsi ETo_Penman()

Blaney Criddle and Radiation Diagram

KODE

In [0]:
def __lat_to_num(lat):
    num, lat = lat.split(' ')
    num = float(num)
    num = -num if lat.lower() == 'lu' else num
    return num

Blaney Criddle

In [0]:
def BC_ETo(c, ETo_x):
    return c * ETo_x

def BC_ETo_x(P, temp):
    return P * (.457*temp + 8.13)

def BC_find_P(latitude, month, table=t_rel_P_LL):
    m = table.loc[month].values
    x = [-5, -2.5, 0, 2.5, 5, 7.5, 10]
    return np.interp(__lat_to_num(latitude), x, m)

def BC_find_C(month, table=t_cor_C_BC, col='C'):
    return table.loc[month, col]
In [0]:
def ETo_BlaneyCriddle(df, temp_col, lat, 
                      as_df=True, report='ETo'):

    # sub_df
    data = df.loc[:, [temp_col]]
    data_array = data.values

    # info_df
    nrows = data.shape[0]

    # initialization
    (P, ETo_x, C, ETo) = (np.zeros(nrows) for _ in range(4))

    # calculation
    temp = data_array[:, 0]
    month = data.index.month.values

    for i in range(nrows):
        P[i]        = BC_find_P(lat, month[i])
        ETo_x[i]    = BC_ETo_x(P[i], temp[i])
        C[i]        = BC_find_C(month[i])
        ETo[i]      = BC_ETo(C[i], ETo_x[i])

    if report.lower() == 'full':
        results = np.stack((
            month, temp, P, ETo_x, C, ETo
        ), axis=1)
        columns_name = [
            'Month', 'Temp', 'P', 'ETo_x', 'C', 'ETo'
        ]
    elif report.lower() == 'eto':
        results = ETo
        columns_name = ['ETo']
    
    if as_df:
        return pd.DataFrame(
            data=results, index=data.index, columns=columns_name
        )
    else:
        return results

RADIASI

In [0]:
def RAD_ETo(C, ETo_x):
    return C * ETo_x

def RAD_ETo_x(w, Rs):
    return w * Rs

def RAD_find_W(temp, table=t_rel_T_W):
    t = table['suhu'].values
    w = table['W'].values
    return np.interp(temp, t, w)

def RAD_find_Rg(latitude, month, table=t_val_Rg):
    m = table.loc[month].values
    x = [-5, -4, -2, 0, 2, 4, 6, 8, 10]
    return np.interp(__lat_to_num(latitude), x, m)

def RAD_Rs(sun_duration, Rg):
    return (0.25 + 0.54 * sun_duration/100) * Rg

def RAD_find_C(month, table=t_cor_C_RAD, col='C'):
    return table.loc[month, col]
In [0]:
def ETo_Radiation(df, temp_col, sun_col, lat,
                  as_df=True, report='ETo'):

    # sub_df
    data = df.loc[:, [temp_col, sun_col]]
    data_array = data.values

    # info_df
    nrows = data.shape[0]

    # initialization
    (w, Rg, Rs, ETo_x, C, ETo) = (np.zeros(nrows) for _ in range(6))

    # calculation
    temp = data_array[:, 0]
    sun = data_array[:, 1]
    month = data.index.month.values

    for i in range(nrows):
        w[i]        = RAD_find_W(temp[i])
        Rg[i]       = RAD_find_Rg(lat, month[i])
        Rs[i]       = RAD_Rs(sun[i], Rg[i])
        ETo_x[i]    = RAD_ETo_x(w[i], Rs[i])
        C[i]        = RAD_find_C(month[i])
        ETo[i]      = RAD_ETo(C[i], ETo_x[i])


    if report.lower() == 'full':
        results = np.stack((
            month, temp, sun, w, Rg, Rs, ETo_x, C, ETo
        ), axis=1)
        columns_name = [
            'Month', 'Temp', 'Sun', 'W', 'R_G', 'R_s', 'ETo_x', 'C', 'ETo'
        ]
    elif report.lower() == 'eto':
        results = ETo
        columns_name = ['ETo']
    
    if as_df:
        return pd.DataFrame(
            data=results, index=data.index, columns=columns_name
        )
    else:
        return results

Penman

In [0]:
def PEN_ETo(C, ETo_x):
    return C * ETo_x

def PEN_ETo_x(w, Rs, Rn1, fU, e_g, e_d):
    return w * (0.75 * Rs - Rn1) + (1 - w) * fU * (e_g - e_d)

def PEN_Rs(sun_duration, RG):
    return (0.25 + 0.54 * sun_duration / 100) * RG

def PEN_Rn1(ft, fe_d, fsun):
    return ft * fe_d * fsun

def PEN_fe_d(e_d):
    return 0.34 - 0.044 * np.sqrt(e_d)

def PEN_e_d(e_d, RH):
    return e_d * RH /100

def PEN_fsun(sun_duration):
    return 0.1 + (0.9 * sun_duration / 100)

def PEN_fU(U):
    return 0.27 * (1 + 0.864 * U)

def PEN_find_from_T(temp, table=t_rel_T_PEN):
    t = table['suhu'].values
    e = table['e_mbar'].values
    w = table['w'].values
    ft = table['f_t'].values

    return (
        np.interp(temp, t, e),
        np.interp(temp, t, w),
        np.interp(temp, t, ft)
    )

def PEN_find_Rg(latitude, month, table=t_val_Rg):
    m = table.loc[month].values
    x = [-5, -4, -2, 0, 2, 4, 6, 8, 10]
    return np.interp(__lat_to_num(latitude), x, m)

def PEN_find_C(month, table=t_cor_C_PEN, col='C'):
    return table.loc[month, col]
In [0]:
def ETo_Penman(df, temp_col, humid_col, wind_col, sun_col, lat,
               as_df=True, report='ETo'):

    # sub_df
    data = df.loc[:, [temp_col, humid_col, wind_col, sun_col]]
    data_array = data.values

    # info_df
    nrows = data.shape[0]

    # initialization
    (e_g, w, ft, fed, e_d, Rg, Rs, fsun,
     fU, Rn1, C, ETo_x, ETo) = (np.zeros(nrows) for _ in range(13))

    # calculation
    temp = data_array[:, 0]
    RH = data_array[:, 1]
    wind = data_array[:, 2]
    sun = data_array[:, 3]
    month = data.index.month.values

    for i in range(nrows):
        e_g[i], w[i], ft[i] = PEN_find_from_T(temp[i])
        e_d[i]      = PEN_e_d(e_g[i], RH[i])
        fed[i]      = PEN_fe_d(e_d[i])
        Rg[i]       = PEN_find_Rg(lat, month[i])
        Rs[i]       = PEN_Rs(sun[i], Rg[i])
        fsun[i]     = PEN_fsun(sun[i])
        fU[i]       = PEN_fU(wind[i])
        Rn1[i]      = PEN_Rn1(ft[i], fed[i], fsun[i])
        C[i]        = PEN_find_C(month[i])
        ETo_x[i]    = PEN_ETo_x(w[i], Rs[i], Rn1[i], fU[i], e_g[i], e_d[i])
        ETo[i]      = PEN_ETo(C[i], ETo_x[i])


    if report.lower() == 'full':
        results = np.stack((
            month, temp, RH, wind, sun, e_g, w, ft, e_d, fed, Rg, Rs,
            fsun, fU, Rn1, C, ETo_x, ETo
        ), axis=1)
        columns_name = [
            'Month', 'Temp', 'Humidity', 'Wind', 'Sun', 'e_g', 'w', 'f_t',
            'e_d', 'f_e_d', 'R_g', 'R_s', 'f_sun', 'f_U', 'R_n1', 'C', 'ETo_x',
            'ETo'
        ]
    elif report.lower() == 'eto':
        results = ETo
        columns_name = ['ETo']
    
    if as_df:
        return pd.DataFrame(
            data=results, index=data.index, columns=columns_name
        )
    else:
        return results

FUNGSI

Dalam modul ini terdapat tiga fungsi utama yang memiliki rumus berbeda yaitu ETo_BlaneyCriddle() untuk rumus Blaney Criddle, ETo_Radiation() untuk rumus Radiasi, ETo_Penman() untuk rumus Penman.

Ketiga ini memiliki keserupaan argumen dalam fungsi, berikut argumen yang terdapat pada ketiga fungsi:

  • df: dataset objek pandas.DataFrame dengan index berupa objek DatetimeIndex atau yang serupa.
  • lat: posisi lintang lokasi dalam bentuk string seperti 5 LU/3.5 LS.
  • (opsional) as_df=True (default), keluaran berupa pandas.DataFrame. Keluaran berupa numpy.ndarray jika False.
  • (opsional) report='ETo' (default), terdapat beberapa nilai yang dapat diterima oleh argumen report:
    • full: keluaran akan menyertakan seluruh peubah yang dihitung dalam model.
    • Eto: keluaran hanya menyertakan kolom ETo (hasil akhir model).

Fungsi ETo_BlaneyCriddle()

Selain argumen posisi df, lat dan argumen opsional as_df, report, fungsi ini membutuhkan argumen temp_col yaitu nama kolom suhu.

In [0]:
sample = pd.DataFrame(
    data=[27.2,27.4,27.3,27.9,27.9,26.8,26.9,28.6,27.9,27.5,27.7,27.3],
    index=pd.date_range('20010101', periods=12, freq='MS'),
    columns=['temp_C']
)

default, as_df=False

In [11]:
ETo_BlaneyCriddle(df=sample, 
                  temp_col='temp_C', 
                  lat='7.5 LS',
                  as_df=False)
Out[11]:
array([4.6055296, 4.6260032, 4.327281 , 4.0925388, 3.9463767, 3.8513664,
       3.8600037, 4.2930405, 4.6771872, 4.63624  , 4.6567136, 4.7806152])

report='full'

In [12]:
ETo_BlaneyCriddle(df=sample, 
                  temp_col='temp_C', 
                  lat='7.5 LS',
                  report='full')
Out[12]:
Month Temp P ETo_x C ETo
2001-01-01 1.0 27.2 0.28 5.756912 0.80 4.605530
2001-02-01 2.0 27.4 0.28 5.782504 0.80 4.626003
2001-03-01 3.0 27.3 0.28 5.769708 0.75 4.327281
2001-04-01 4.0 27.9 0.28 5.846484 0.70 4.092539
2001-05-01 5.0 27.9 0.27 5.637681 0.70 3.946377
2001-06-01 6.0 26.8 0.27 5.501952 0.70 3.851366
2001-07-01 7.0 26.9 0.27 5.514291 0.70 3.860004
2001-08-01 8.0 28.6 0.27 5.724054 0.75 4.293041
2001-09-01 9.0 27.9 0.28 5.846484 0.80 4.677187
2001-10-01 10.0 27.5 0.28 5.795300 0.80 4.636240
2001-11-01 11.0 27.7 0.28 5.820892 0.80 4.656714
2001-12-01 12.0 27.3 0.29 5.975769 0.80 4.780615

Fungsi ETo_Radiation()

Selain argumen posisi df, lat dan argumen opsional as_df, report, fungsi ini membutuhkan argumen:

  • temp_col yaitu nama kolom suhu.
  • sun_col yaitu nama kolom penyinaran matahari.
In [0]:
sample_2 = pd.DataFrame({
    'temp': [25.4, 25.7, 26.1, 26.9, 26.4, 26.6, 26.2, 25.9,
             26.6, 27.9, 28.4, 27.2],
    'sun': [41.8, 41.8, 53.6, 49.1, 60.0, 63.6, 60.9, 57.3,
            61.8, 65.5, 61.8, 47.3],
}, index=pd.date_range('20010101', periods=12, freq='MS'))

default, report='full'

In [14]:
ETo_Radiation(df=sample_2, 
              temp_col='temp', 
              sun_col='sun', 
              lat='7.5 LS', 
              report='full')
Out[14]:
Month Temp Sun W R_G R_s ETo_x C ETo
2001-01-01 1.0 25.4 41.8 0.749 16.025 7.623413 5.709936 0.80 4.567949
2001-02-01 2.0 25.7 41.8 0.752 16.075 7.647199 5.750694 0.80 4.600555
2001-03-01 3.0 26.1 53.6 0.756 15.225 8.212974 6.209008 0.75 4.656756
2001-04-01 4.0 26.9 49.1 0.764 14.250 7.340745 5.608329 0.75 4.206247
2001-05-01 5.0 26.4 60.0 0.759 13.175 7.562450 5.739900 0.75 4.304925
2001-06-01 6.0 26.6 63.6 0.761 12.500 7.418000 5.645098 0.75 4.233823
2001-07-01 7.0 26.2 60.9 0.757 12.800 7.409408 5.608922 0.75 4.206691
2001-08-01 8.0 25.9 57.3 0.754 13.775 7.706010 5.810332 0.80 4.648266
2001-09-01 9.0 26.6 61.8 0.761 14.925 8.712021 6.629848 0.80 5.303878
2001-10-01 10.0 27.9 65.5 0.774 15.775 9.523368 7.371086 0.80 5.896869
2001-11-01 11.0 28.4 61.8 0.779 15.950 9.310334 7.252750 0.80 5.802200
2001-12-01 12.0 27.2 47.3 0.767 15.925 8.048813 6.173440 0.80 4.938752

Fungsi ETo_Penman()

Selain argumen posisi df dan lat; argumen opsional as_df dan report, fungsi ini membutuhkan argumen:

  • temp_col: nama kolom suhu.
  • sun_col: nama kolom penyinaran matahari.
  • humid_col: nama kolom kelembaban.
  • wind_col: nama kolom kecepatan angin.
In [0]:
sample_3 = pd.DataFrame({
    'temp': [25.4, 25.7, 26.1, 26.9, 26.4, 26.6, 26.2, 25.9,
             26.6, 27.9, 28.4, 27.2],
    'humid': [76.6, 79.6, 74.4, 78.0, 79.1, 78.8, 79.6, 79.8,
              77.4, 79.4, 77.2, 77.7],
    'wind': [2.3, 1.8, 1.9, 2.0, 2.0, 2.4, 2.5, 3.0, 3.3, 3.0, 2.0, 2.4],
    'sun': [41.8, 41.8, 53.6, 49.1, 60.0, 63.6, 60.9, 57.3, 61.8, 
            65.5, 61.8, 47.3]
}, index=pd.date_range('20010101', periods=12, freq='MS'))
In [16]:
ETo_Penman(df=sample_3, 
           temp_col='temp', 
           humid_col='humid', 
           wind_col='wind', 
           sun_col='sun', 
           lat='7.5 LS', 
           report='full')
Out[16]:
Month Temp Humidity Wind Sun e_g w f_t e_d f_e_d R_g R_s f_sun f_U R_n1 C ETo_x ETo
2001-01-01 1.0 25.4 76.6 2.3 41.8 32.450 0.749 15.750 24.85670 0.120631 16.025 7.623413 0.4762 0.806544 0.904754 1.1 5.141999 5.656198
2001-02-01 2.0 25.7 79.6 1.8 41.8 32.525 0.752 15.825 25.88990 0.116119 16.075 7.647199 0.4762 0.689904 0.875055 1.1 4.790219 5.269241
2001-03-01 3.0 26.1 74.4 1.9 53.6 33.820 0.756 15.920 25.16208 0.119288 15.225 8.212974 0.5824 0.713232 1.106015 1.1 5.327334 5.860068
2001-04-01 4.0 26.9 78.0 2.0 49.1 35.455 0.764 16.080 27.65490 0.108613 14.250 7.340745 0.5419 0.736560 0.946428 0.9 4.839053 4.355148
2001-05-01 5.0 26.4 79.1 2.0 60.0 34.420 0.759 15.980 27.22622 0.110413 13.175 7.562450 0.6400 0.736560 1.129221 0.9 4.724821 4.252339
2001-06-01 6.0 26.6 78.8 2.4 63.6 34.830 0.761 16.020 27.44604 0.109489 12.500 7.418000 0.6724 0.829872 1.179394 0.9 4.800835 4.320751
2001-07-01 7.0 26.2 79.6 2.5 60.9 34.020 0.757 15.940 27.07992 0.111031 12.800 7.409408 0.6481 0.853200 1.147031 0.9 4.777259 4.299533
2001-08-01 8.0 25.9 79.8 3.0 57.3 32.920 0.754 15.875 26.27016 0.114481 13.775 7.706010 0.6157 0.969840 1.118960 1.1 5.100576 5.610634
2001-09-01 9.0 26.6 77.4 3.3 61.8 34.830 0.761 16.020 26.95842 0.111545 14.925 8.712021 0.6562 1.039824 1.172601 1.1 6.036265 6.639892
2001-10-01 10.0 27.9 79.4 3.0 65.5 37.590 0.774 16.280 29.84646 0.099620 15.775 9.523368 0.6895 0.969840 1.118236 1.1 6.360059 6.996065
2001-11-01 11.0 28.4 77.2 2.0 61.8 38.700 0.779 16.380 29.87640 0.099499 15.950 9.310334 0.6562 0.736560 1.069471 1.1 6.042748 6.647023
2001-12-01 12.0 27.2 77.7 2.4 47.3 36.090 0.767 16.140 28.04193 0.107000 15.925 8.048813 0.5257 0.829872 0.907870 1.1 5.489920 6.038912

Changelog

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