模型模拟值一般为连续的值，实测值为点，且一般带有误差
# 导入必要的包，定义文件读取函数
import os
import datetime
import pandas as pd
import numpy as np
from pcse.fileinput import CABOFileReader # 读取CABO文件（作物土壤文件）
from pcse.fileinput import YAMLAgroManagementReader # 管理文件读取
from pcse.models import Wofost71_WLP_FD, Wofost71_PP # 导入模型，Wofost71_PP潜在模型
from pcse.base import ParameterProvider # 综合作物、土壤、位点参数
from pcse.fileinput import ExcelWeatherDataProvider # 读取气象文件
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()
def soil_data(s, c, om, b=1.5):
s = s # sand %
c = c # clay %
om = om # organic matter %
b = b # bulk g/cm3
# 永久萎蔫点 PWP
pwp = 1.14*(-0.00024*s+0.00487*c+0.006*om+0.00005*s*om-0.00013*c*om+0.0000068*s*c+0.031)-0.02
# 田间持水量 FC
o33 = -0.00251*s+0.00195*c+0.011*om+0.00006*s*om-0.00027*c*om+0.0000452*s*c+0.299
fc = o33+1.283*o33*o33-0.374*o33-0.015
# 饱和含水量-FC
s_33 = 1.636*(0.00278*s+0.00034*c+0.022*om-0.00018*s*om-0.00027*c*om-0.0000584*s*c+0.078)-0.107
# 饱和含水量 容重矫正后
sat = fc+s_33-0.00097*s+0.043
# 校正后田间持水量
pn = 2.65*(1-(fc+s_33-0.00097*s+0.043))
df = b/1.5
if df<0.9:
df=0.9
elif df>1.2:
df=1.2
pdf = pn*df
fc_v = fc-0.2*(fc+s_33-0.00097*s+0.043-(1-pdf/2.65))
# 校正后饱和含水量
sat_v = 1-(pdf/2.65)
return [pwp, fc_v, sat_v]
def st_loc(num):
"""
选取气象数据的时候标准化输入为0.5的分辨率
"""
import math
if num-math.floor(num) <=0.25 or num-math.floor(num)>=0.75:
result = round(num)
elif 0.25 < num-math.floor(num) < 0.75:
result = math.floor(num) + 0.5
return result
# 设置绘图地点
sites = ['内黄', '许昌', '洛宁', '罗山']
sites_e = ['NH', 'XC', 'LN', 'LS']
annote = ['a', 'b', 'c', 'd']
data_dir = r'C:\Users\Administrator\Desktop\model calibration' # 模型校准文件夹
weather_dir = r'C:\Users\Administrator\Desktop\2020气象数据' # 气象数据路径
variety = 207# 品种
data_obs = pd.read_excel('实测值.xlsx', sheet_name='实测值') # 实测值
tagp_rmse = [] # bioass rmse
twso_sim = [] # yield simulation
twso = [] # yield observation
twso_std = [] # yield standard deviation
fig = plt.figure(figsize=(10,14)) # 建立图表
for i in range(4):
# 运行模型
site = sites[i]
cropdata = CABOFileReader(os.path.join(data_dir,'%s.CAB'%site)) # 读取作物文件
soildata = CABOFileReader(os.path.join(data_dir,'EC3.NEW')) # 土壤文件
sitedata = {'SSMAX' : 0.,
'IFUNRN' : 0,
'NOTINF' : 0,
'SSI' : 0,
'WAV' : 20,
'SMLIM' : 0.03,
'RDMSOL' : 120}
parameters = ParameterProvider(cropdata=cropdata, soildata=soildata, sitedata=sitedata) # 参数集合
data_base_info = pd.read_excel('18县生物量LAI整理12.27.xlsx', sheet_name='基本情况')
sub_data = data_base_info.loc[data_base_info['试验地']==site,
['经度', '维度', '品种', '播量（计算）', '播种期',
'sa', 'cl', 'bd', 'som', '灌溉', 'WAV']]
# 基本数据获取
lon = sub_data.loc[sub_data['品种']==variety, ['经度']] # 读取经度
lat = sub_data.loc[sub_data['品种']==variety, ['维度']] # 读取纬度
SAND = sub_data.loc[sub_data['品种']==variety, ['sa']].iloc[0,0] # 砂粒
CLAY = sub_data.loc[sub_data['品种']==variety, ['cl']].iloc[0,0] # 黏粒
OM = sub_data.loc[sub_data['品种']==variety, ['som']].iloc[0,0] # 有机质
BULK_DENSITY = sub_data.loc[sub_data['品种']==variety, ['bd']].iloc[0,0] # 容重
sow_date = sub_data.loc[sub_data['品种']==variety, ['播种期']].iloc[0,0] # 播期
irrigation = sub_data.loc[sub_data['品种']==variety, ['灌溉']].iloc[0,0] # 灌溉条件
# 数据替换
parameters['TDWI'] = sub_data.loc[sub_data['品种']==variety, ['播量（计算）']].iloc[0,0] # 播量
parameters['SMW'] = soil_data(SAND, CLAY, OM, BULK_DENSITY)[0] # 萎蔫点
parameters['SMFCF'] = soil_data(SAND, CLAY, OM, BULK_DENSITY)[1] # 田间持水量
parameters['SM0'] = soil_data(SAND, CLAY, OM, BULK_DENSITY)[2] # 饱和含水量
parameters['WAV'] = sub_data.loc[sub_data['品种']==variety, ['WAV']].iloc[0,0]
# 气象数据
weatherdataprovider = ExcelWeatherDataProvider(os.path.join(weather_dir,'NASA天气文件lat={0:.1f},lon={1:.1f}.xlsx'.
format(st_loc(lat.iloc[0,0]), st_loc(lon.iloc[0,0]))))
agromanagement = YAMLAgroManagementReader(os.path.join(data_dir,'%s.agro'%site)) # 管理文件读取
wf = Wofost71_WLP_FD(parameters, weatherdataprovider, agromanagement) # 定义模型
wf.run_till_terminate() # 运行模型直到终止
output=pd.DataFrame(wf.get_output()).set_index('day')
# 实测值
sub_obs = data_obs.loc[data_obs['试验地']==site, ['品种', 'day', 'LAI', '生物量kg/ha',
'穗生物量kg/ha', '生物量标准差', '穗标准差']]
sub_obs_103 = sub_obs.loc[sub_obs['品种']==variety].set_index('day')
# ----------------- 绘图 ----------------------------------------------------
# 建子图
axes = plt.subplot(3,2,i+1) # 建立子图坐标域 三行两列 画前四个图
# 单位转换 Kg-Mg
sub_obs_103['生物量kg/ha'] = sub_obs_103['生物量kg/ha']/1000
sub_obs_103['穗生物量kg/ha'] = sub_obs_103['穗生物量kg/ha']/1000
sub_obs_103['生物量标准差'] = sub_obs_103['生物量标准差']/1000
# 绘制误差图errorbar
axes.errorbar(sub_obs_103.index, sub_obs_103['生物量kg/ha'], yerr=sub_obs_103['生物量标准差'], fmt=".", ms=10,capsize=5, label = 'Observation')
output['TAGP'] = output['TAGP']/1000
output['TWSO'] = output['TWSO']/1000
# 生物量模拟图Simulation
output.TAGP.plot(ax=axes, label = 'Simulation', lw=3)
output['TAGP'].dropna(inplace=True)
output['TWSO'].fillna(method='pad', inplace=True) # 填充后面的空值为最后一个非空值，即最终产量
# 计算模拟与实测值的差与RMSE
df_differences_tagp = sub_obs_103['生物量kg/ha'] - output['TAGP']
tagp_rmse.append(np.sqrt(np.mean(df_differences_tagp**2)))
# yield obs
twso.append(sub_obs_103['穗生物量kg/ha'][-1])
# yield sim
twso_sim.append(output['TWSO'][-1])
twso_std.append(sub_obs_103['穗标准差'][-1]/1000)
# 计算R方
r2 = 1 - np.sum((df_differences_tagp)**2) / np.sum((sub_obs_103['生物量kg/ha'] - np.mean(sub_obs_103['生物量kg/ha']))**2)
# 设置字体格式
font = {'family' : 'Helvetica',
'weight' : 'normal',
'size' : 16,
}
axes.set_xlabel('') # 去除x轴标题
# 融合坐标轴 （下次考虑使用subplots的sharex,sharey）
if i == 1 or i == 3:
axes.set_ylabel('') # 去除y轴标题
# 去除坐标轴标签
plt.tick_params(axis='y', which='both', labelleft=False,top=False, left=False, labelbottom=False)
else:
axes.set_ylabel('Aboveground biomass (Mg/ha)', font)
if i == 0 or i == 1:
plt.tick_params(axis='x',which='both', top=False, bottom=False, labelbottom=False)
else:
plt.xticks(fontproperties = 'Helvetica' ,size = 14,rotation = 30)
# 标注R和RMSE,\it表示斜体
axes.annotate("$\it{RMSE}$ = %.2f Mg/ha" % tagp_rmse[i], xy=('2019-10-20', 10), family='Helvetica', weight='normal', fontsize=16)
plt.text('2019-10-20', 17.8, annote[i], family='Helvetica', weight='normal', fontsize=17)
axes.annotate("$\it{R^2}$ = %.2f" % r2, xy=('2019-10-20', 12), family='Helvetica', weight='normal', fontsize=16)
# 设置坐标轴字体格式
plt.tick_params(axis='y', labelsize=15)
labels = axes.get_xticklabels() + axes.get_yticklabels()
[label.set_fontname('Helvetica') for label in labels]
plt.ylim(0, 19)
plt.legend(loc=(0.02, 0.72), prop=font)
# fig.subplots_adjust(wspace=0.05,hspace=0.05)
axes = plt.subplot(3,1,3) # 绘制第三列大图
df = pd.DataFrame({'Observation': twso, 'Simulation': twso_sim})
df.plot(kind='bar', ax=axes, yerr = [twso_std, [0,0,0,0]], error_kw = {'capsize' :4}, alpha=1)
plt.xticks([0,1,2,3],sites_e,fontproperties = 'Helvetica',rotation = 0, size=18)
axes.set_ylabel('Yield (Mg/ha)', font)
plt.text(-0.4, 7.8, 'e', family='Helvetica', weight='normal', fontsize=17)
labels = axes.get_xticklabels() + axes.get_yticklabels()
[label.set_fontname('Helvetica') for label in labels]
plt.tick_params(axis='y',labelsize=15) # 可以指定是x或y坐标轴参数
plt.legend(prop=font)
plt.tight_layout()
plt.savefig('simulation.tif',bbox_inches='tight', dpi=600) # bbox_inches防止保存丢失信息
结果
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
