365 lines
14 KiB
Python
365 lines
14 KiB
Python
from datetime import datetime, timedelta
|
||
from pyeto import fao,convert
|
||
from pyeto.fao import svp_from_t, avp_from_rhmean, delta_svp, psy_const,atm_pressure
|
||
import pandas as pd
|
||
#1、计算日参考蒸散量
|
||
def calculateET0(tmean, tmax, tmin, rh_mean, wind_speed, solar_rad, elevation, lat, doy):
|
||
# 计算饱和水汽压和实际水汽压
|
||
svp_max = svp_from_t(tmax)
|
||
svp_min = svp_from_t(tmin)
|
||
svp_mean =(svp_max+svp_min)/2
|
||
avp = avp_from_rhmean(svp_max, svp_min, rh_mean)
|
||
|
||
# 计算净辐射(简化版)
|
||
# 将纬度从度转换为弧度
|
||
lat_rad = convert.deg2rad(lat)
|
||
sol_dec = fao.sol_dec(doy)
|
||
sha = fao.sunset_hour_angle(lat_rad, sol_dec)
|
||
ird = fao.inv_rel_dist_earth_sun(doy)
|
||
et_rad = fao.et_rad(lat_rad, sol_dec, sha, ird)
|
||
cs_rad = fao.cs_rad(elevation, et_rad)
|
||
nisr = fao.net_in_sol_rad(solar_rad)
|
||
nolr = fao.net_out_lw_rad(tmin, tmax, solar_rad, cs_rad, avp)
|
||
rn = fao.net_rad(nisr, nolr)
|
||
deltasvp= delta_svp(tmean)
|
||
psy = psy_const(atm_pressure(elevation))
|
||
# 计算ET₀
|
||
et0 = fao.fao56_penman_monteith(
|
||
net_rad=rn,
|
||
t=tmean,
|
||
ws=wind_speed,
|
||
svp=svp_mean,
|
||
avp=avp,
|
||
delta_svp=deltasvp,
|
||
psy=psy
|
||
)
|
||
return et0
|
||
def cal_et0_List(planting_date,monitor_date,date_params):
|
||
et0_list = []
|
||
current_date = planting_date
|
||
while current_date < monitor_date:
|
||
result_params=date_params[current_date.strftime("%Y-%m-%d")]
|
||
et0=calculateET0(result_params['tmean'], result_params['tmax'], result_params['tmin'], result_params['rh_mean'], result_params['wind_speed'],
|
||
result_params['solar_rad'], result_params['elevation'],result_params['lat'], result_params['doy'])
|
||
et0_list.append(et0)
|
||
current_date += timedelta(days=1) # 增加一天
|
||
return et0_list
|
||
|
||
|
||
#计算实际每天作物蒸散量
|
||
# 1. 定义甘草生育阶段与作物系数（参考苜蓿）
|
||
# growth_stages = {
|
||
# "initial": {"duration_days": 20, "kc": 0.40}, # 生长初期
|
||
# "development": {"duration_days": 40, "kc": 0.90}, #生长快速期
|
||
# "mid": {"duration_days": 20, "kc": 1.1}, # 生长缓慢期
|
||
# "growth_slowly": {"duration_days": 30, "kc": 0.6}, # 生长缓慢期
|
||
# "dormant_period": {"duration_days": 150, "kc": 0.40}, # 营养储备期
|
||
# "fanqing_year2":{"duration_days": 30, "kc": 0.35},#返青期
|
||
# "development_year2":{"duration_days": 90, "kc": 0.9}, #第二年快速生长期
|
||
# "maturity _year2": {"duration_days": 40, "kc": 0.6}, # 第二年成熟期
|
||
# }
|
||
growth_stages_roots1 = {
|
||
'initial': {
|
||
'root_length': {5: 0.8, 10: 0.2}, # 主根长：80%概率在0-15cm，平均约8cm
|
||
'duration_days': 23, # 6月播种后1-2周
|
||
'kc': 0.40, # 作物系数
|
||
'time_range': '2025-06-17至2025-07-10',
|
||
'hight':0.03,#m
|
||
'k_rain':0.25
|
||
},
|
||
'development': {
|
||
'root_length': {20: 0.6, 45: 0.4}, # 主根长：60%概率在20cm，40%概率在45cm，平均约30cm
|
||
'duration_days': 40, # 发芽后-8月20日左右（约2个月）
|
||
'kc': 0.9, # 作物系数
|
||
'time_range': '2025-07-10至2025-08-20',
|
||
'hight':0.15,
|
||
'k_rain':0.25
|
||
},
|
||
'mid': {
|
||
'root_length': {30: 0.5, 80: 0.5}, # 主根长：50%概率在30cm，50%概率在80cm，平均约55cm
|
||
'duration_days': 30, # 8月20日-9月下旬
|
||
'kc': 1.1, # 作物系数
|
||
'time_range': '2025-08-21至2025-09-20',
|
||
'hight':0.45,
|
||
'soil_evaporation_depth': 0.65,
|
||
'k_rain':0.4
|
||
},
|
||
'growth_slowly': {
|
||
'root_length':{40: 0.4,60:0.4, 80: 0.2}, # 根据数据估算平均主根长
|
||
'duration_days': 30, # 9月下旬-10月下旬
|
||
'kc': 0.6,
|
||
'time_range': '2025-09-21至2025-10-20',
|
||
'hight':0.6
|
||
},
|
||
'dormant_period': {
|
||
'root_length': {40: 0.3,60:0.4, 80: 0.2}, # 冬季休眠，主根基本停止生长
|
||
'duration_days': 160, # 11月-次年3月
|
||
'kc': 0.4,
|
||
'time_range': '2025-10-21至2026-03-30',
|
||
'hight':0.8,
|
||
'k_rain':0.5
|
||
},
|
||
'fanqing_year2': {
|
||
'root_length': {40: 0.2,60:0.6, 80: 0.2}, # 返青期主根长变化不大
|
||
'duration_days': 30, # 3月下旬-4月下旬
|
||
'kc': 0.35,
|
||
'time_range': '2026-04-01至2026-04-30',
|
||
'hight':0.5,
|
||
'k_rain':0.5
|
||
},
|
||
'development_year2': {
|
||
'root_length': {40: 0.1,60:0.7, 80: 0.2}, # 根据数据估算平均主根长
|
||
'duration_days': 90, # 4月下旬-7月下旬
|
||
'kc': 0.9,
|
||
'time_range': '2026-05-01至2026-07-30',
|
||
'hight':0.45,
|
||
'k_rain':0.65
|
||
},
|
||
'maturity_year2': {
|
||
'root_length':{40: 0.1,60:0.6, 80: 0.3}, # 根据数据估算平均主根长
|
||
'duration_days': 40, # 7月下旬-8月下旬
|
||
'kc': 0.6,
|
||
'time_range': '2026-08-01至2026-09-10',
|
||
'hight':0.65,
|
||
'k_rain':0.5
|
||
}
|
||
}
|
||
# growth_stages_roots = {
|
||
# 'initial': {
|
||
# 'root_length': {5: 0.8, 10: 0.2}, # 主根长：80%概率在0-15cm，平均约8cm
|
||
# 'duration_days': 20, # 6月播种后1-2周
|
||
# 'kc': 0.40, # 作物系数
|
||
# 'time_range': '2025-07-08至2025-07-28',
|
||
# 'hight':0.03,#m
|
||
# 'k_rain':0.25
|
||
# },
|
||
# 'development': {
|
||
# 'root_length': {20: 0.6, 45: 0.4}, # 主根长：60%概率在20cm，40%概率在45cm，平均约30cm
|
||
# 'duration_days': 42, # 发芽后-8月20日左右（约2个月）
|
||
# 'kc': 0.9, # 作物系数
|
||
# 'time_range': '2025-07-29至2025-09-10',
|
||
# 'hight':0.15,
|
||
# 'k_rain':0.25
|
||
# },
|
||
# 'mid': {
|
||
# 'root_length': {30: 0.5, 80: 0.5}, # 主根长：50%概率在30cm，50%概率在80cm，平均约55cm
|
||
# 'duration_days': 30, # 8月20日-9月下旬
|
||
# 'kc': 1.1, # 作物系数
|
||
# 'time_range': '2025-09-11至2025-10-10',
|
||
# 'hight':0.45,
|
||
# 'soil_evaporation_depth': 0.65,
|
||
# 'k_rain':0.4
|
||
# },
|
||
# 'growth_slowly': {
|
||
# 'root_length':{40: 0.4,60:0.4, 80: 0.2}, # 根据数据估算平均主根长
|
||
# 'duration_days': 31, # 9月下旬-10月下旬
|
||
# 'kc': 0.6,
|
||
# 'time_range': '2025-10-11至2025-11-10',
|
||
# 'hight':0.6
|
||
# },
|
||
# 'dormant_period': {
|
||
# 'root_length': {40: 0.3,60:0.4, 80: 0.2}, # 冬季休眠，主根基本停止生长
|
||
# 'duration_days': 140, # 11月-次年3月
|
||
# 'kc': 0.4,
|
||
# 'time_range': '2025-11-11至2026-03-31',
|
||
# 'hight':0.8,
|
||
# 'k_rain':0.5
|
||
# },
|
||
# 'fanqing_year2': {
|
||
# 'root_length': {40: 0.2,60:0.6, 80: 0.2}, # 返青期主根长变化不大
|
||
# 'duration_days': 30, # 3月下旬-4月下旬
|
||
# 'kc': 0.35,
|
||
# 'time_range': '2026-04-01至2026-04-30',
|
||
# 'hight':0.5,
|
||
# 'k_rain':0.5
|
||
# },
|
||
# 'development_year2': {
|
||
# 'root_length': {40: 0.1,60:0.7, 80: 0.2}, # 根据数据估算平均主根长
|
||
# 'duration_days': 90, # 4月下旬-7月下旬
|
||
# 'kc': 0.9,
|
||
# 'time_range': '2026-05-01至2026-07-30',
|
||
# 'hight':0.45,
|
||
# 'k_rain':0.65
|
||
# },
|
||
# 'maturity_year2': {
|
||
# 'root_length':{40: 0.1,60:0.6, 80: 0.3}, # 根据数据估算平均主根长
|
||
# 'duration_days': 40, # 7月下旬-8月下旬
|
||
# 'kc': 0.6,
|
||
# 'time_range': '2026-08-01至2026-09-10',
|
||
# 'hight':0.65,
|
||
# 'k_rain':0.5
|
||
# }
|
||
# }
|
||
# ---------------------------
|
||
# 2. 计算每日Kc值（根据生育阶段动态分配）
|
||
# ---------------------------
|
||
def get_daily_kc(growth_stages, total_days=423):
|
||
daily_kc = []
|
||
for stage in growth_stages.values():
|
||
daily_kc.extend([stage["kc"]] * stage["duration_days"])
|
||
# 如果总天数超出定义，用末期Kc填充剩余天数
|
||
if len(daily_kc) < total_days:
|
||
daily_kc.extend([growth_stages["maturity_year2"]["kc"]] * (total_days - len(daily_kc)))
|
||
return daily_kc[:total_days] # 截断至总天数
|
||
|
||
# ---------------------------
|
||
# 3. 计算每日实际蒸散量ETc(得到历史列表)
|
||
# ---------------------------
|
||
def calculate_etc_List(et0_daily, planting_date,monitor_date, kc_daily):
|
||
"""
|
||
计算每日ETc（实际蒸散量）
|
||
:param eto_daily: 每日参考蒸散量列表（单位：cm/d）
|
||
:param planting_date: 播种日期，格式"YYYY-MM-DD"
|
||
:param kc_daily: 每日Kc值列表
|
||
:return: 日期列表、ETc列表
|
||
"""
|
||
start_date = planting_date
|
||
delta=monitor_date-start_date
|
||
days = delta.days
|
||
dates = [start_date + timedelta(days=i) for i in range(days)]
|
||
date_strings = [date.strftime("%Y-%m-%d") for date in dates]
|
||
kc_daily = kc_daily[:days]
|
||
etc = [eto * kc for eto, kc in zip(et0_daily, kc_daily)]
|
||
return etc,date_strings,kc_daily
|
||
#基于作物蒸散和土壤蒸发计算实际蒸散量
|
||
def calculate_etc_List_new(et0_daily, planting_date,monitor_date, kc_daily,ke_daily):
|
||
"""
|
||
计算每日ETc（实际蒸散量）
|
||
:param eto_daily: 每日参考蒸散量列表（单位：cm/d）
|
||
:param planting_date: 播种日期，格式"YYYY-MM-DD"
|
||
:param kc_daily: 每日Kc值列表
|
||
:return: 日期列表、ETc列表
|
||
"""
|
||
start_date = planting_date
|
||
delta=monitor_date-start_date
|
||
days = delta.days
|
||
dates = [start_date + timedelta(days=i) for i in range(days)]
|
||
date_strings = [date.strftime("%Y-%m-%d") for date in dates]
|
||
kc_daily = kc_daily[:days]
|
||
etc = [eto * (kc+ke) for eto, kc,ke in zip(et0_daily, kc_daily,ke_daily)]
|
||
return etc,date_strings,kc_daily
|
||
def calculate_etc_dayily(et0_daily, planting_date,monitor_date, kc_daily):
|
||
"""
|
||
计算每日ETc（实际蒸散量）
|
||
:param eto_daily: 每日参考蒸散量列表（单位：cm/d）
|
||
:param planting_date: 播种日期，格式"YYYY-MM-DD"
|
||
:param kc_daily: 每日Kc值列表
|
||
:return: 日期列表、ETc列表
|
||
"""
|
||
start_date = planting_date
|
||
delta=monitor_date-start_date
|
||
days = delta.days
|
||
# dates = [start_date + timedelta(days=i) for i in range(len(kc_daily))]
|
||
kc_daily = kc_daily[days]
|
||
etc = et0_daily*kc_daily
|
||
return etc
|
||
def export_etc_et0_date(et0_list,etc_list,dates,kc_list,ke_daily,soil_data,rain_daily,irrigate_daily,dk_name):
|
||
df = pd.DataFrame({
|
||
"Date": dates,
|
||
"ET0": [float(f"{x:.3f}") for x in et0_list],
|
||
"KC": kc_list,
|
||
"KE":ke_daily,
|
||
"ETc": [float(f"{x:.3f}") for x in etc_list],
|
||
"rain_effective":rain_daily,
|
||
"I_effective":irrigate_daily,
|
||
"VWC_fc":0.2884,
|
||
"VWC_wp":0.1086,
|
||
"dkbm":dk_name
|
||
})
|
||
df['VWC_SUR'] = 0.0
|
||
df['VWC_MID'] = 0.0
|
||
df['VWC_BOT'] = 0.0
|
||
for i in range(len(df)):
|
||
date=df.at[i, 'Date']
|
||
df.at[i, 'VWC_SUR']=round(float(soil_data["SOIL_MOISTURE_SURFACE"][date]),3)
|
||
df.at[i, 'VWC_MID'] = round(float(soil_data["SOIL_MOISTURE_MIDDLE"][date]),3)
|
||
df.at[i, 'VWC_BOT'] = round(float(soil_data["SOIL_MOISTURE_BOTTOM"][date]),3)
|
||
return df
|
||
|
||
|
||
def estimate_mad_from_et(max_et):
|
||
|
||
# 根据最大日蒸散量查表确定MAD值
|
||
if max_et < 2:
|
||
mad = 0.68
|
||
elif 2 <= max_et < 3:
|
||
mad = 0.58
|
||
elif 3 <= max_et < 4:
|
||
mad = 0.48
|
||
elif 4 <= max_et < 5:
|
||
mad = 0.4
|
||
elif 5 <= max_et < 6:
|
||
mad = 0.35
|
||
elif 6 <= max_et < 7:
|
||
mad = 0.33
|
||
elif 7 <= max_et < 8:
|
||
mad = 0.28
|
||
elif 8 <= max_et < 9:
|
||
mad = 0.25
|
||
else: # max_et >= 8
|
||
mad = 0.23
|
||
return mad
|
||
def calculate_weighted_soil_parameters(root_distribution, soil_layers):
|
||
# 初始化累计变量
|
||
sum_fc = 0.0
|
||
sum_wp = 0.0
|
||
# 处理每个根系分布层
|
||
for root_depth, root_ratio in root_distribution['root_length'].items():
|
||
# 找到包含该深度的土壤层
|
||
for layer in soil_layers:
|
||
min_depth, max_depth = layer['depth_range']
|
||
if min_depth <= root_depth < max_depth:
|
||
# 计算该层贡献 (简单比例分配)
|
||
layer_thickness = max_depth - min_depth
|
||
depth_in_layer = min(root_depth, max_depth) - min_depth
|
||
contribution_ratio = depth_in_layer / layer_thickness * root_ratio
|
||
sum_fc += layer['theta_33'] * contribution_ratio
|
||
sum_wp += layer['theta_1500'] * contribution_ratio
|
||
break
|
||
|
||
return sum_fc, sum_wp
|
||
def cal_stagebyDate(monitor_date,zwlx_growth_stages):
|
||
for stage, details in zwlx_growth_stages.items():
|
||
time_range_str = details['time_range']
|
||
if type(monitor_date) is str:
|
||
monitor_date= datetime.strptime(monitor_date, '%Y-%m-%d')
|
||
if '-' in time_range_str:
|
||
start_str, end_str = time_range_str.split('至')
|
||
start_date = datetime.strptime(start_str, '%Y-%m-%d')
|
||
end_date = datetime.strptime(end_str, '%Y-%m-%d')
|
||
if start_date <= monitor_date <= end_date:
|
||
return stage,details['root_length'],details['hight'],details['k_rain']
|
||
raise ValueError(f"监测日期 {monitor_date.strftime('%Y-%m-%d')} 不在任何生育阶段范围内")
|
||
#根据日期计算主根长度
|
||
def calculate_average_root_length(root_length_dist):
|
||
return sum(length * prob for length, prob in root_length_dist.items())
|
||
def cal_soil_fc_wp():
|
||
# 土壤分层数据 (33kPa和1500kPa对应的含水率)
|
||
soil_data = [
|
||
{'depth_range': (0, 5), 'theta_33': 0.296, 'theta_1500': 0.103},
|
||
{'depth_range': (5, 15), 'theta_33': 0.284, 'theta_1500': 0.104},
|
||
{'depth_range': (15, 30), 'theta_33': 0.286, 'theta_1500': 0.105},
|
||
{'depth_range': (30, 60), 'theta_33': 0.290, 'theta_1500': 0.114},
|
||
{'depth_range': (60, 100), 'theta_33': 0.286, 'theta_1500': 0.117},
|
||
{'depth_range': (100, 200), 'theta_33': 0.284, 'theta_1500': 0.114}
|
||
]
|
||
# theta_fc, theta_wp = calculate_weighted_soil_parameters(
|
||
# growth_stages_roots[stage],
|
||
# soil_data
|
||
# )
|
||
theta_fc=0.2884
|
||
theta_wp=0.1086
|
||
return theta_fc,theta_wp
|
||
|
||
if __name__ == "__main__":
|
||
# #cal_et0_List(datetime(2026, 6, 18),datetime(2026, 6, 22))
|
||
et0_daily=[23,45,67,89]
|
||
# etc,dates,kc_daily=calculate_etc_List(et0_daily, datetime(2026, 6, 18),datetime(2026, 6, 22), get_daily_kc(growth_stages_roots))
|
||
# calculate_etc_dayily(45, datetime(2026, 6, 18), datetime(2026, 6, 22), get_daily_kc(growth_stages))
|
||
# df=export_etc_et0_date(et0_daily,etc,dates,kc_daily)
|
||
# print(df)
|
||
# cal_soil_fc_wp('initial')
|
||
# stage,details = cal_stagebyDate(datetime(2025, 7, 18))
|
||
# print(details)
|
||
# rd=calculate_average_root_length(details)
|
||
# print(rd) |