蒸散发计算：Hargreaves ET0 方程(公式、注意细节，python代码)

最近在使用日最高最低气温计算潜在蒸散发能力，从气温数据准备，到Hargreaves方程计算，也是几经波澜。写下这篇博客，一方面是为了记录，另外一方面也是想为正在做相关计算的你，减少一些无用功~
网上Hargreaves方程计算的教程基本都没写清楚，部分csdn下载包虽然收费，但那个代码真是一言难尽呀，代码混乱就算了，还没个清楚的说明文件，看的就很费劲。而且很多都是从GitHub上无脑复制过来，未经任何改动就直接拿来卖钱的。抄袭别人代码来卖钱，服务还做不好，嘴脸真是有够丑恶哈~
废话不说，咱直接开始 Hargreaves ET0 方程的介绍趴

Hargreaves ET0 方程计算公式

当太阳辐射、相对温度和（或）风速数据缺失时，或者说只有日最高最低气温时，可以选择使用Hargreaves ET0 方程来计算潜在蒸散发，其公式为
$$E{T}_0=0.0023(T_{mean}+17.8)(T_{max}+T_{min}{)}^{0.5}{R}_a$$
其中，$T_{mean},T_{max},T_{min}$分别为日平均气温、最高气温和最低气温，单位都为$C°$。其中$T_{mean}$如果没有小时气温数据来平均计算，可以用$T_{mean}=(T_{max}+T_{min})/2$来简化替代。$R_a$是天顶辐射，单位是$mmda{y}^{-1}$（这个单位需要注意一下，我就是在这被卡了好久的(哭死~））。
下面来详细介绍天顶辐射$R_a$的计算公式
$$R_a=\frac{24\times 60}{\pi }{G}_{sc}{d}_r[w_{s}sin(\phi )sin(\delta )+cos(\phi )cos(\delta )sin(w_s)]$$
公式中$R_a$单位为$MJ{m}^{-2}da{y}^{-1}$，注意这里单位和Hargreaves 方程里$R_a(mmda{y}^{-1})$的单位不一样。
$G_{sc}$为太阳常数为0.0820$(MJ{m}^{-2}mi{n}^{-1})$
$d_r$为日地间相对距离的倒数
$w_s$为太阳时角（rad）
$\phi$为地理纬度（rad）
$\delta$为太阳磁偏角（rad）
下面分辨介绍$\phi ,d_r,\delta ,w_s$的具体计算公式
首先$\phi$
$$\phi =\frac{\pi }{180}\times [十进制度数纬度]$$
为防止理解错，这里给了个例子~
日地间相对距离的倒数 $d_r$ 和太阳磁偏角$\delta$的计算见下式：
$$d_r=1+0.033cos(\frac{2\pi }{365}J)$$
$$\delta =0.409sin(\frac{2\pi }{365}J-1.39)$$
公式中： $J$是年内某天的日序数，在 1（ 1 月 1 日）至 365 或 366（ 12 月 31 日）之间。
日落时角$w_s$ 见下式：
$$w_s=arccos[-tan(\phi )tan(\delta )]$$
在算法语言里没有反余弦函数，所以日落时角${\omega }_s$ 亦可用反正切函数计算：
$$w_s=\frac{\pi }{2}-arctan[\frac{-tan(\phi )tan(\delta )}{X^{0.5}}]$$
其中：
$$X=1-[tan(\phi ){]}^2[tan(\delta ){]}^2$$
当 $X\le 0$时，$X=0.00001$
因为公式Hargreaves 方程里 $R_a$单位为$mmda{y}^{-1}$，所以需要把$R_a(MJ{m}^{-2}da{y}^{-1})$转换成$R_a(mmda{y}^{-1})$。具体转换如下：
$$R_a(mmda{y}^{-1})=0.408\times R_a(MJ{m}^{-2}da{y}^{-1})$$

Hargreaves ET0 python代码

这里给了我的Hargreaves方程代码。因为这是适应我的数据的代码，利用插值后的站点气温数据（2922,140,230）(时间,纬度,经度)计算南方[99°E~122°E, 21°N ~ 35°N]地区0.1°分辨率的蒸散发能力
在你们的程序上，应该不能直接运行。不过我给了详尽的注释，希望能给你们带来一些帮助。
先 $R_a$ 的计算

## 计算天顶辐射Ra
import h5py
import numpy as np
import pandas as pd
import math   #导入包

with h5py.File('temperature_daily.h5', 'r') as file:
    lons = file['longitude'][:]
    lats = file['latitude'][:] 
#导入进度纬度文件lons是[140]的一维数组，lats是[230]的一维数组
    
j=np.array(range(1,366),dtype=np.int32)
for ii in range(2016,2023):
    if ii % 4 == 0:
        mid = np.array(range(1,367),dtype=np.int32)
        j = np.hstack([j, mid]) 
    else:
        mid = np.array(range(1,366),dtype=np.int32)
        j = np.hstack([j, mid])
    del mid       # j 年内某天的日序数。把2015年到2022年的日序数都给标好，赋给j变量


hudu = np.empty((lats.shape[0],lons.shape[0]), dtype=np.float32)

for ii in range(lats.shape[0]):
    hudu[ii,:] = lats[ii]*math.pi/180   #hudu是φ 地理纬度，(140,230)

dr = np.empty((lats.shape[0],lons.shape[0]), dtype=np.float32)  #dr 日地间相对距离的倒数
fi = np.empty((lats.shape[0],lons.shape[0]), dtype=np.float32)  #fi 是 δ 太阳磁偏角
Ra = np.empty((2922,lats.shape[0],lons.shape[0]), dtype=np.float32) # Ra 是天顶辐射

for jj in range(2922):
    dr[:,:] = 1+0.033*math.cos(2*math.pi*j[jj]/365)
    fi[:,:] = 0.409*math.sin(2*math.pi*j[jj]/365-1.39)
    x = 1-np.tan(hudu)**2*np.tan(fi)**2 
    x = np.where(x<=0, 0.00001, x)
    ws = math.pi/2-np.arctan((-np.tan(hudu))*np.tan(fi)/(x**0.5))
    Ra[jj,:,:] = (
        24*60/math.pi*0.082*dr*
        (ws*np.sin(fi)*np.sin(hudu)+np.cos(hudu)*np.cos(fi)*np.sin(ws))
        )*0.408
    
with h5py.File('Ra.h5', 'w') as h5file:
	h5file.create_dataset('Ra', data=Ra)
	h5file.create_dataset('longitude', data=lons)
	h5file.create_dataset('latitude', data=lats)
	h5file.attrs['Description'] = 'ra 2015~2022'   #保存变量至HDF5文件

再计算ET0

import h5py
import numpy as np
import pandas as pd
import math #导入包

with h5py.File('temperature_daily.h5', 'r') as file:
    temperature_max = file['temperature_max'][:]
    temperature_min = file['temperature_min'][:]
    temperature_mean = file['station_temputure_mean'][:] #导入气温数据都是（2922，140，230）的数组

with h5py.File('Ra.h5', 'r') as file:
    Ra = file['Ra'][:] #导入气温 （2922，140，230）

et = 0.0023*(temperature_mean+17.8)*((temperature_max-temperature_min)**0.5)*Ra #ET0计算

with h5py.File('Ra.h5', 'a') as h5file:
	h5file.create_dataset('ET0', data=et) #保存至 HDF5文件

基本表现

多年空间平均

时间序列平均

参考

FAO Irrigation and Drainage Paper

https://blog.csdn.net/weixin_51956867/article/details/136112272?fromshare=blogdetail&sharetype=blogdetail&sharerId=136112272&sharerefer=PC&sharesource=m0_46211544&sharefrom=from_link