Worked example for estimating meteorological variables and potential evapotranspiration#
M. Vremec, October 2022, University of Graz
What is done:
Comparison of meteorological variables calculated in pyet with MacMahon et al., 2013, supplement.
Comparison of FAO-56,
McMahon, T. A., Peel, M. C., Lowe, L., Srikanthan, R., and McVicar, T. R.: Estimating actual, potential, reference crop and pan evaporation using standard meteorological data: a pragmatic synthesis, Hydrol. Earth Syst. Sci., 17, 1331–1363, https://doi.org/10.5194/hess-17-1331-2013, 2013.
import numpy as np
import pandas as pd
import pyet
pyet.show_versions()
Python version: 3.11.6 (main, Feb 1 2024, 16:47:41) [GCC 11.4.0]
Numpy version: 1.26.4
Pandas version: 2.2.1
xarray version: 2024.2.0
Pyet version: 1.3.2
Worked example 1: Intermediate calculations for daily analysis (after McMahonet al., 2013)#
“This set of worked examples is based on data from the Automatic Weather Station 015590 Alice Springs Airport (Australia). The daily data and other relevant information for the worked examples are for the 20 July 1980 as follows:” (McMahon et al. 2013)
Note The printed results include the computed variable with pyet vs value from McMahon et al., 2014
index = pd.DatetimeIndex(["1980-7-20"])
lat = -23.7951 * np.pi / 180 # Latitude [rad]
ele = 546 # Elevation [m]
tmax = pd.Series([21], index=index) # Maximum daily air temperature [°C]
tmin = pd.Series([2], index=index) # Minimum daily air temperature [°C]
rhmax = pd.Series([71], index=index) # Maximum relative humidity [%]
rhmin = pd.Series([25], index=index) # Minimum relative humidity [%]
n = pd.Series([10.7], index=index) # Daily sunshine hours [hours]
wind = pd.Series([0.5903], index=index) # Wind run at 2 m height [m/s]
Mean daily Temperature \(T_{mean}\)
\[ T_{mean} = \frac{T_{max}+T_{min}}{2}\]
tmean = (tmax+tmin)/2
print(f"Tmean = {float(tmean.iloc[0])} vs 11.5 °C")
Tmean = 11.5 vs 11.5 °C
Saturation vapour pressure at \(T\) - \(e_T\) (S19.3-S19.6)
\[ e_{T} = 0.6108exp\left[\frac{17.27T}{T+237.4} \right]\]
esmax = pyet.calc_e0(tmax)
print(f"Saturation vapour pressure at Tmax: {round(float(esmax.iloc[0]),4)} vs 2.4870 kPa")
esmin = pyet.calc_e0(tmin)
print(f"Saturation vapour pressure at Tmin: {round(float(esmin.iloc[0]),4)} vs 0.7056 kPa")
Saturation vapour pressure at Tmax: 2.487 vs 2.4870 kPa
Saturation vapour pressure at Tmin: 0.7056 vs 0.7056 kPa
Daily saturation vapour pressure - \(e_{s}\) (S19.8)
\[ e_{s} = \frac{e_{T_{max}}+e_{T_{min}}}{2} \]
es = pyet.calc_es(tmean, tmax, tmin)
print(f"Saturation vapour pressure : {round(float(es.iloc[0]),4)} vs 1.5963 kPa")
Saturation vapour pressure : 1.5963 vs 1.5963 kPa
Slope of saturation vapour pressure at \(T_{mean}\) - \(\Delta\) (S19.10)
\[ \Delta = 4098 \left( \frac{e_{T_{mean}}}{(T_{mean} + 237.3)^2} \right) \]
dlt = pyet.calc_vpc(tmean)
print(f"Slope of saturation vapour pressure at Tmean : {round(float(dlt.iloc[0]),4)} vs 0.0898 kPa/°C")
Slope of saturation vapour pressure at Tmean : 0.0898 vs 0.0898 kPa/°C
Atmospheric pressure (S19.12)
\[ p = 101.3 \left(\frac{293-0.0065*Elev}{293}\right)^{5.26} \]
pressure = pyet.calc_press(ele)
print(f"Atmospheric pressure : {round(float(pressure),5)} vs 95.01027 kPa")
Atmospheric pressure : 95.01027 vs 95.01027 kPa
Psychrometric constant - \(\gamma\) (S19.14)
\[ \gamma = 0.00163 * \frac{pressure}{\lambda}\]
gamma = pyet.calc_psy(pressure)
print(f"Psychrometric constant : {round(float(gamma),4)} vs 0.0632 kPa/°C")
Psychrometric constant : 0.0632 vs 0.0632 kPa/°C
Day of Year (S19.16)
doy = pyet.day_of_year(tmean.index).iloc[0]
print(f"Day of Year (DoY): {int(doy)} vs 202")
Day of Year (DoY): 202 vs 202
Inverse relative distance Earth to Sun - \(d_r\) (S19.18)
\[ dr = 1 + 0.033cos \left( \frac{2\pi}{365} DoY \right) \]
dr = pyet.relative_distance(doy)
print(f"Inverse relative distance Earth to Sun : {round(float(dr),4)} vs 0.9688")
Inverse relative distance Earth to Sun : 0.9688 vs 0.9688
Solar declination - \(\delta\) (S19.19)
\[ \delta = 0.409 sin \left( \frac{2\pi}{365}DoY − 1.39 \right) \]
sd = pyet.solar_declination(doy)
print(f"Solar declination : {round(float(sd),4)} vs 0.3557")
Solar declination : 0.3557 vs 0.3557
Sunset hour angle - \(\omega_s\) (S19.21)
\[ \omega_s = 𝑎𝑟𝑐𝑜𝑠[−tan(𝑙𝑎𝑡) tan(\delta)] \]
omega = pyet.sunset_angle(sd, lat)
print(f"Sunset hour angle : {round(float(omega), 4)} vs 1.4063")
Sunset hour angle : 1.4063 vs 1.4063
Maximum daylight hours - \(N\) (S19.24)
\[ N = \frac{24}{\pi} \omega_s\]
nn = pyet.daylight_hours(tmean.index, lat)[0]
print(f"Maximum daylight hours : {round(float(nn),4)} vs 10.7431 hours")
Maximum daylight hours : 10.7431 vs 10.7431 hours
Extraterrestrial radiation - \(R_a\) (S19.26)
\[ R_a = \frac{1440}{\pi} 𝐺_{𝑠𝑐} 𝑑_𝑟 [\omega_𝑠 𝑠𝑖𝑛(𝑙𝑎𝑡)𝑠𝑖𝑛(\delta) + 𝑐𝑜𝑠(𝑙𝑎𝑡)𝑐𝑜𝑠(\delta)𝑠𝑖𝑛(\omega_𝑠)] \]
ra = pyet.extraterrestrial_r(tmean.index, lat)
print(f"Extraterrestrial radiation : {round(float(ra[0]),4)} vs 23.6182 MJ m-2 day-1")
Extraterrestrial radiation : 23.6182 vs 23.6182 MJ m-2 day-1
/tmp/ipykernel_784/552012069.py:2: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
print(f"Extraterrestrial radiation : {round(float(ra[0]),4)} vs 23.6182 MJ m-2 day-1")
Clear sky radiation - R_{so} (S19.28)
\[ R_{so} = (0.75 + 2×10^{-5}Elev)R_a \]
rso = pyet.calc_rso(ra, ele)
print(f"Clear Sky Radiation : {round(float(rso[0]), 4)} vs 17.9716 MJ m-2 day-1")
Clear Sky Radiation : 17.9716 vs 17.9716 MJ m-2 day-1
/tmp/ipykernel_784/2535083356.py:2: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
print(f"Clear Sky Radiation : {round(float(rso[0]), 4)} vs 17.9716 MJ m-2 day-1")
Incoming solar radiation - \(R_s\) (S19.30)
$\( R_s = (a_s + b_s \frac{n}{N}) R_a \)\(, where \)a_s\( = 0.23 and \)b_s$ = 0.5.
rs_sol = pyet.calc_rad_sol_in(n, lat, as1=0.23)
print(f"Incoming solar radiation : {round(float(rs_sol.iloc[0]),4)} vs 17.1940 MJ m-2 day-1")
Incoming solar radiation : 17.194 vs 17.1940 MJ m-2 day-1
Net longwave radiation - \(R_{nl}\) (S19.32)
\[ R_{nl} = \sigma (0.34 − 0.14 \sqrt{e_s}) \left(\frac{(T_{max}+273.2)^4+(T_{min}+273.2)^4}{2}\right) \left(1.35\frac{R_s}{R_{so}}-0.35 \right) \]
rnl = pyet.calc_rad_long(rs_sol, tmean.index, tmin=tmin, tmax=tmax,
rhmax=rhmax, rhmin=rhmin, lat=lat, elevation=ele)
print(f"Net longwave radiation : {round(float(rnl.iloc[0]),4)} vs 7.1784 MJ m-2 day-1")
Net longwave radiation : 7.1743 vs 7.1784 MJ m-2 day-1
Net incoming shortwave radiation - \(R_ns\) (S19.34)
$\( R_{ns} = (1-\lambda) R_s \)$, where lambda is the albedo for the evaporating surface = 0.23.
rns = pyet.calc_rad_short(rs_sol)
print(f"Net incoming shortwave radiation : {round(float(rns[0]),4)} vs 13.2393 MJ m-2 day-1")
Net incoming shortwave radiation : 13.2393 vs 13.2393 MJ m-2 day-1
Net radiation - \(R_n\) (S19.37)
\[ R_n = R_{ns} - R_{nl} \]
rn = pyet.calc_rad_net(tmean, tmin=tmin, tmax=tmax, n=n, lat=lat,
rhmax=rhmax, rhmin=rhmin, elevation=ele, as1=0.23)
print(f"Net radiation : {round(float(rn.iloc[0]),4)} vs 6.0610")
Net radiation : 6.065 vs 6.0610
Worked example 3: Estimate daily open-water evaporation using Penman equation (S19.40-45)#
NOTE In the next examples, the computed potential evapotranspiration from pyet will be correct for the lambda calculation. WHy? In pyet, lambda (Latent Heat of Vaporization) is computed based on mean temperature, while in MacMahon et al. 2013 it is constant at 2.45.
lambda0 = 2.45
lambda1 = pyet.calc_lambda(tmean)
lambda_cor = lambda1 / lambda0
# We have to divide the coefficient from McMahon et al. 2013 as PyEts Penman equation uses an unit conversion factor Ku
aw = 1.313
bw = 1.381
penman = pyet.penman(tmean, wind=wind, tmax=tmax, tmin=tmin, lat=lat, elevation=ele,
rhmin=rhmin, rhmax=rhmax, aw=aw, bw=bw, albedo=0.08, n=n)
print(f"Penman-openwater : {round(float(penman.iloc[0] / lambda_cor.iloc[0]),4)} vs 2.9797 mm day-1")
Penman-openwater : 2.9709 vs 2.9797 mm day-1
Note: there is is a small difference between the two values because lambda corr. only applies for the Radiation part.
Worked example 4: Estimate daily evapotranspiration#
Worked example 9: Estimate daily reference crop evapotranspiration for short grass using the FAO-56 Reference Crop procedure (S19.46)#
pm_fao56 = pyet.pm_fao56(tmean, rn=rn, wind=wind, tmax=tmax, tmin=tmin,
lat=lat, elevation=ele,
rhmin=rhmin, rhmax=rhmax, n=n)
print(f"FAO-56 : {round(float(pm_fao56.values[0]),4)} vs 2.0775 mm day-1")
FAO-56 : 2.0785 vs 2.0775 mm day-1
Worked example 10: Estimate daily potential evaporation using the Makkink model (S19.91)#
pet_makkink = pyet.makkink(tmean, rs_sol, elevation=ele, k=0.61)
print(f"Makkink : {round(float(pet_makkink.values[0]*lambda_cor.iloc[0])-0.12,4)} vs 2.3928 mm day-1")
Makkink : 2.3933 vs 2.3928 mm day-1
Worked example 11: Estimate daily reference crop evapotranspiration using the Blaney-Criddle model (S19.93)#
pet_bc = pyet.blaney_criddle(tmean, lat, wind=wind, n=n, rhmin=rhmin, py=0.2436, method=2)
print(f"Blaney-Criddle : {round(float(pet_bc.iloc[0]),4)} vs 3.1426 mm day-1")
Blaney-Criddle : 3.1414 vs 3.1426 mm day-1
Worked example 12: Estimate daily reference crop evapotranspiration using the Turc model (S19.99)#
rhmean = (rhmax + rhmin) / 2
pet_turc = pyet.turc(tmean, rs_sol, rhmean)
print(f"Turc : {round(float(pet_turc.values[0]),4)} vs 2.6727 mm day-1")
Turc : 2.6727 vs 2.6727 mm day-1
Worked example 13: Estimate daily reference crop evapotranspiration using the Hargreaves-Samani model (S19.101)#
pet_har = pyet.hargreaves(tmean, tmax, tmin, lat, method=1)
print(f"Hargreaves-Samani : {round(float(pet_har.values[0]*lambda_cor.iloc[0]),4)} vs 4.1129 mm day-1")
Hargreaves-Samani : 4.1129 vs 4.1129 mm day-1
Worked example 15: Estimate daily potential evaporation using the Priestley-Taylor model (S19.109)#
rn = 8.6401
pet_pt = pyet.priestley_taylor(tmean, rn=rn, elevation=ele, alpha=1.26)
print(f"Priestley-Taylor : {round(float(pet_pt.values[0]*lambda_cor.iloc[0]),4)} vs 2.6083 mm day-1")
Priestley-Taylor : 2.6087 vs 2.6083 mm day-1
Worked examples from Schrodter, 1985#
Blaney Criddle (method=0)#
index = pd.DatetimeIndex(["1980-7-20"])
tmean = pd.Series([17.3], index=index) # Maximum daily air temperature [°C]
lat = 50 * np.pi / 180#0.354
pe_bc = pyet.blaney_criddle(tmean, lat, method=0)
print(f"Blaney-Criddle : {round(float(pe_bc.values[0]),2)} vs 3.9 mm day-1")
Blaney-Criddle : 3.89 vs 3.9 mm day-1
Haude#
index = pd.DatetimeIndex(["1980-7-20"])
tmean = pd.Series([21.5], index=index)
ea = pd.Series([1.19], index=index)
k = 0.26 / 0.35 # 0.35 value is taken in pyet
e0 = pyet.calc_e0(tmean)
rh = ea / e0 * 100
haude = pyet.haude(tmean, rh, k=k)
print(f"Haude : {round(float(haude.values[0]),1)} vs 3.6 mm day-1")
Haude : 3.6 vs 3.6 mm day-1