入门示例
<p>[TOC]</p>
<p>Meteva 是一个模块化的检验python包，里面基础了各种常用的检验指标（例如ME，RMSE，TS等）的计算算法，高级用户可以组织数据并调用这些算法完成检验计算并绘制结果，但对入门用户（特别是对python不太熟练的用户）来说组织数据和绘制结果可能都比较困难。此种情况下，可学习本篇中的内容作为Meteva的入门。</p>
<pre><code class="language-python">import meteva.base as meb # 该模块用于IO和基础计算
import meteva.method as mem # 该模块基础了检验的基础算法
import meteva.product as mpd # 该模块包含了检验的工具
import numpy as np
import datetime
import copy
import pandas as pd</code></pre>
<p>采用Meteva的检验流程可以分为两个部分：</p>
<ul>
<li>part1: 收集数据</li>
<li>part2: 检验分析</li>
</ul>
<h1>收集数据</h1>
<pre><code class="language-python">###################以下开始为数据收集部分的程序
#设置关注的起始时段
time_begin= datetime.datetime(2019,7,1,8,0)
time_end = datetime.datetime(2019,8,1,8,0)
#读取站点列表，并将站点内容为缺省值，当其作为读取站点数据的参数时，如果站点文件中某个站号不存在时,返回结果中该站点保持为缺省值
station = meb.read_stadata_from_micaps3(r"H:\Example_data\ob\temp_2m\BT19070102.000")
station.iloc[:,-1] = meb.IV
##读取收集观测数据
dir_ob = r"H:\Example_data\ob\temp_2m\BTYYMMDDHH.000" # 观测数据的路径模板
sta_list = []
time0 = time_begin
time_end_ob = time_end + datetime.timedelta(hours = 72)
while time0 < time_end_ob:
path = meb.get_path(dir_ob,time0) # 基于模板和时间生成观测数据路径
#从micaps3格式文件中读取观测数据
sta = meb.read_stadata_from_micaps3(path,station = station,time = time0,dtime = 0,level = 0,data_name = "OBS",show = True)
sta_list.append(sta) # 将数据加入到列表当中
time0 += datetime.timedelta(hours = 3) # 移动时间，循环读取下一个时刻数据
ob_sta_all = meb.concat(sta_list) #将观测数据列表拼接成一个DataFrame
#读取收集ec预报数据
dir_ec = r"H:\Example_data\ecmwf\temp_2m\YYMMDD\BTYYMMDDHH.TTT" #ECMWF data path patten
sta_list =[]
time0 = time_begin
while time0 < time_end:
for dh in range(0,73,3):
path = meb.get_path(dir_ec,time0,dh) # 基于模板，时间和时效生成预报数据路径
#从micaps4格式文件中读取预报数据
grd = meb.read_griddata_from_micaps4(path,show = False)
if grd is not None:
sta = meb.interp_gs_linear(grd,station) #将数据插值到站点
meb.set_stadata_coords(sta,time = time0,dtime = dh,level = 0) #为数据添加时间、时效和层次等坐标信息
meb.set_stadata_names(sta,["ECMWF"]) #为数据添加名称
sta_list.append(sta)
time0 += datetime.timedelta(hours = 12)
ec_sta_all = meb.concat(sta_list) #将ECMWF预报数据列表拼接成一个DataFrame
#读取收集grapes预报数据
dir_grapes = r"H:\Example_data\grapes\temp_2m\YYMMDD\BTYYMMDDHH.TTT"
sta_list =[]
time0 = time_begin
while time0 < time_end:
for dh in range(0,73,3):
path = meb.get_path(dir_grapes,time0,dh) # 基于模板，时间和时效生成预报数据路径
#从micaps4格式文件中读取预报数据
grd = meb.read_griddata_from_micaps4(path,show = False)
if grd is not None:
sta = meb.interp_gs_linear(grd,station) #将数据插值到站点
meb.set_stadata_coords(sta,time = time0,dtime = dh,level = 0) #为数据添加时间、时效和层次等坐标信息
meb.set_stadata_names(sta,["GRAPES"])
sta_list.append(sta)
time0 += datetime.timedelta(hours = 12)
grapes_sta_all = meb.concat(sta_list) #将GRAPES预报数据列表拼接成一个DataFrame
#将观测和预报数据按时空坐标进行匹配,合并成一个DataFrame
sta_all = meb.combine_on_obTime_id(ob_sta_all,[ec_sta_all,grapes_sta_all],need_match_ob=True)
sta_all.to_hdf(r"H:\Example_data\sta_all.h5","df")
###################以上为数据收集部分的程序</code></pre>
<pre><code>success read from H:\Example_data\ob\temp_2m\BT19070108.000
success read from H:\Example_data\ob\temp_2m\BT19070111.000
success read from H:\Example_data\ob\temp_2m\BT19070114.000
success read from H:\Example_data\ob\temp_2m\BT19070117.000
success read from H:\Example_data\ob\temp_2m\BT19070120.000
。。。
success read from H:\Example_data\ob\temp_2m\BT19080323.000
success read from H:\Example_data\ob\temp_2m\BT19080402.000
success read from H:\Example_data\ob\temp_2m\BT19080405.000</code></pre>
<h1>检验分析</h1>
<pre><code class="language-python">#首先导入之前收集好的数据
sta_all = pd.read_hdf(r"H:\Example_data\sta_all.h5")
print(sta_all)</code></pre>
<pre><code> level time dtime id lon lat OBS ECMWF \
0 0 2019-07-01 08:00:00 0 54398 116.6 40.1 25.8 24.724
1 0 2019-07-01 08:00:00 0 54410 116.1 40.6 18.9 20.284
2 0 2019-07-01 08:00:00 0 54416 116.9 40.4 25.1 22.864
3 0 2019-07-01 08:00:00 0 54419 116.6 40.4 27.5 21.796
4 0 2019-07-01 08:00:00 0 54499 116.2 40.2 27.5 23.076
... ... ... ... ... ... ... ... ...
9187 0 2019-07-31 20:00:00 72 54410 116.1 40.6 18.8 21.556
9188 0 2019-07-31 20:00:00 72 54416 116.9 40.4 26.7 24.960
9189 0 2019-07-31 20:00:00 72 54419 116.6 40.4 27.0 22.920
9190 0 2019-07-31 20:00:00 72 54499 116.2 40.2 29.3 23.056
9191 0 2019-07-31 20:00:00 72 54412 116.6 40.7 25.4 21.640
GRAPES
0 23.844
1 21.740
2 22.704
3 23.136
4 23.736
... ...
9187 19.968
9188 23.404
9189 23.056
9190 22.568
9191 21.232
[9192 rows x 9 columns]</code></pre>
<pre><code class="language-python"># 按照不同的时效分组检验模式的平均绝对误差，并绘制成图形
#其中mem.mae是 平均绝对误差计算的算法函数名
#g="dtime" 代表按时效进行分组
#plot = "bar"表示以柱状图形式绘制检验结果
result = mpd.score(sta_all,mem.mae,g= "dtime",plot = "bar") </code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/2fd334bd9962497658f415bb4b000693" alt="" /></p>
<pre><code class="language-python"># 按照起报时间分组检验模式的平均绝对误差，并绘制成图形
#其中mem.me是 平均误差计算的算法函数名
#g="time" 代表按时间进行分组
#plot = "line"表示以折线图形式绘制检验结果
result = mpd.score(sta_all,mem.rmse,g= "time",plot = "line") </code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/7aa69b86df88afa100760bd47d988255" alt="" /></p>
<pre><code class="language-python"># 按照起报时间的时钟数分组检验模式的TS评分，并绘制成图形
#其中mem.ts是 TS评分算法函数名
#grade_list = [30,35,37,40]，是TS评分检验的阈值
#g="hour" 代表按时间时钟数进行分组
#plot = "bar"表示以折线图形式绘制检验结果
result = mpd.score(sta_all,mem.ts,grade_list = [30,35,37,40],g = "hour",plot = "bar") </code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/09e0cb22067f19ec88e66a7590942124" alt="" /></p>
<pre><code class="language-python"># 按照起报时间-预报时效联合分组对模式的平均误差进行检验，并绘制成图形
#其中mem.me是 平均绝对误差算法函数名
#x="time_dtime" 是横纵坐标
result = mpd.score_tdt(sta_all,mem.me,x_y = "time_dtime") </code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/70e94c0d805cfc79bcb75e40a51db8b1" alt="" /></p>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/b8461b35bdf12c4fb74315ad6a0b7816" alt="" /></p>
<pre><code class="language-python"># 绘制平均误差的空间分布
#其中mem.me是 平均绝对误差算法函数名
# subplot = "member"设置用于将多个模式绘制成多个子图
#ncol 用于设置subplot中子图列数
# width = 10 用户设置Fig的宽度
result = mpd.score_id(sta_all,mem.me,subplot = "member",ncol = 2,width = 10) </code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/af7f3108161e9e125541ad6094f40c8a" alt="" /></p>
<pre><code class="language-python"># 绘制观测和预报的回归关系
#其中mem.scatter_regress是回归图绘制函数
result = mpd.plot(sta_all,mem.scatter_regress) </code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=3f1b0920ff4b953cbbf6445b567ea62b" alt="" /></p>
<pre><code class="language-python"># 绘制观测和预报的频率对应关系
#其中mem.frequency_histogram是频率图绘制函数
#grade_list = np.arange(10,40,5).tolist()是频率统计时等级的设置方式
result = mpd.plot(sta_all,mem.frequency_histogram,grade_list = np.arange(10,40,5).tolist()) </code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/8a5317d09cea5d473651f041b7047a8c" alt="" /></p>