python计算penman蒸散（FAO-56）简化地表净辐射计算过程

penman FAO-56版本公式如下:

 

FAO-56官方文件

        本次计算利用aquacropeto包来简化地表净辐射的计算过程，使用较为容易获取的最高温、最低温、风速、平均相对湿度、日照时间、气压六项气象数据以及海拔、纬度、日序数（日序数计算请另行查找资料）即可完成计算。

       需说明的是实际水汽压的计算在此采用了平均相对湿度来进行，但根据FAO官方文件来看，这并不作为首选，其他数据如最低相对湿度、最高相对湿度以及露点温度资料可获取时建议采用其他计算方式！

import aquacropeto as pyeto
import pandas as pd
import math
data = pd.read_csv('D:/station/data_ET0.csv')
df = pd.DataFrame(data)
#计算净短波辐射
solar_dec=pyeto.sol_dec(df['index']) #输入日序数 计算太阳赤纬
latitude = pyeto.deg2rad(df['纬度']) #纬度转换单位
sunset_hour_angle=pyeto.sunset_hour_angle(latitude,solar_dec)  #输入纬度、太阳赤纬，计算日落时角
inv_rel_dist_earth_sun = pyeto.inv_rel_dist_earth_sun(df['index']) #输入日序数计算地球太阳之间的逆相对距离
Ra = pyeto.et_rad(latitude,solar_dec,sunset_hour_angle,inv_rel_dist_earth_sun) #计算天文辐射
N= pyeto.daylight_hours(sunset_hour_angle) #输入日落时角，计算白天长度
n = df['日照时数']
Rs = (0.25+(0.5*n/N))*Ra #as bs均取值为FAO默认值
Rns = 0.77*Rs  #反射率取了FAO默认值
#计算净长波辐射
Tmax = pyeto.celsius2kelvin(df['最高温']) #转换单位为开尔文 仅在计算净长波辐射时使用！！！
Tmin = pyeto.celsius2kelvin(df['最低温']) #转换单位为开尔文 仅在计算净长波辐射时使用！！！
cs_rad = pyeto.cs_rad(df['测量海拔'],Ra) #通过海拔、与天文辐射计算晴空辐射
def e0(T):
    e = 0.6108*(math.e**((17.27*T/(T+237.3))))
    return e
es = 0.5*( e0(df['最高温'])+e0(df['最低温']) )
ea = 0.01*df['平均相对湿度']*es #计算实际水汽压
#计算净长波辐射
Rnl = pyeto.net_out_lw_rad(Tmin,Tmax,Rs,cs_rad,ea)
#计算净辐射
Rn = pyeto.net_rad(Rns,Rnl)

#计算平均温度 必须为最高温与最低温的平均
tmean = (df['最高温']+df['最低温'])/2
#计算饱和水汽压曲线斜率
delta = pyeto.delta_svp(tmean)
#计算湿度计常数
r =0.665*0.001*df['平均气压'] #输入气压单位应为kpa!
#将10米风速转换为2m风速 若为2m风速数据可忽略
u2 = df['平均风速']*4.87/math.log(672.58)

#penman FAO-56
ET0_1 = 0.408 * delta * (Rn)
ET0_2 = r * 900* u2 * (es-ea)/(tmean+273)
ET0_3 = delta + r * (1 + 0.34*u2)
ET0 = (ET0_1 + ET0_2)/ET0_3