地面要素统计检验
<p>[TOC]</p>
<pre><code class="language-python">%matplotlib inline
%load_ext autoreload
%autoreload 2
import pandas as pd
import meteva.method as mem
import meteva.base as meb
import meteva.product as mpd
import meteva.perspact as mps # 透视分析模块
import datetime
import meteva
import numpy as np
import os
import copy</code></pre>
<h1>地面要素站点统计检验</h1>
<p>在需要对长时间序列的预报数据进行检验时，如果数据量很大，且预报数据的读取非常耗时，也可以采用并行的方法来提高效率</p>
<table>
<thead>
<tr>
<th style="text-align: left;">参数</th>
<th style="text-align: left;">说明</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>para</font></strong></td>
<td style="text-align: left;">统计的参数</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;">统计量中间结果输出到文件中，函数不返回结果</td>
</tr>
</tbody>
</table>
<p><strong> 参数示例和说明</strong> </p>
<h2>降水预报检验实例</h2>
<pre><code class="language-python">para= {
&quot;mid_method&quot;: meteva.method.hfmc, # 统计检验的总量计算函数，
&quot;grade_list&quot;: [0.1,10,25,50,100], # 二分类检验所需的等级参数
&quot;compare&quot;: &quot;&gt;=&quot;, # 二分类检验所需的对比方式参数
&quot;begin_time&quot;:datetime.datetime(2024,9,1,8), #时段开始时刻(基于起报时间)
&quot;end_time&quot;:datetime.datetime(2024,9,3,20), #时段结束时刻（基于起报时间）
&quot;time_type&quot;:&quot;UT&quot;, # 最终检验结果呈现时，采用北京时还是世界时，UT代表世界时，BT代表北京时
&quot;time_step&quot;: 12, # 起报时间间隔
&quot;middle_result_path&quot;: r&quot;H:\test_data\output\mps\hfmc.h5&quot;, #检验中间量数据的存储路径
&quot;station_file&quot;:meteva.base.station_国家站, #检验站点表的存储路径
&quot;defalut_value&quot;:meteva.base.IV,
&quot;interp&quot;: meteva.base.interp_gs_nearest,
&quot;recover&quot;:False, #False表示保留已经收集好的中间量，补齐未收集的数据，True表示删除已有的统计结果文件，重新开始收集
&quot;cpu&quot;:2,
&quot;ob_data&quot;:{
&quot;dir_ob&quot;: r&quot;I:\rc\RAIN\rain0\YYYYMMDDHH.000&quot;, # 实况场数据路径
&quot;hour&quot;:None,
&quot;read_method&quot;: meteva.base.io.read_stadata_from_micaps3, #读取数据的函数
&quot;read_para&quot;: {}, #读取数据的函数参数
&quot;reasonable_value&quot;: [0, 1000], #合理的观测值的取值范围，超出范围观测将被过滤掉
&quot;operation&quot;:meteva.base.fun.sum_of_sta, #实况数据读取后处理函数
&quot;operation_para&quot;: {&quot;used_coords&quot;: [&quot;time&quot;], &quot;span&quot;: 24}, #实况数据读取后处理参数
&quot;time_type&quot;: &quot;BT&quot;, # 数据文件中的时间类型，UT代表世界时
},
&quot;fo_data&quot;:{
&quot;ECMWF&quot;: {
&quot;dir_fo&quot;: r&quot;S:\data\grid\ECMWF_HR\APCP\YYYYMMDD\YYMMDDHH.TTT.nc&quot;, # 数据路径
&quot;hour&quot;: [8, 20, 12], #一天中的起报时长，第一个时次，最后一个时次，时次间隔
&quot;dtime&quot;:[0,48,12], #检验的时效
&quot;read_method&quot;: meteva.base.io.read_griddata_from_nc, #读取数据的函数
&quot;read_para&quot;: {}, #读取数据的函数参数
&quot;reasonable_value&quot;: [0, 1000], #合理的预报值的取值范围，超出范围观测将被过滤掉
&quot;operation&quot;: meteva.base.fun.change, #预报数据读取后处理函数
&quot;operation_para&quot;: {&quot;used_coords&quot;: &quot;dtime&quot;, &quot;delta&quot;: 24}, #预报数据读取后处理参数
&quot;time_type&quot;: &quot;BT&quot;, # 预报数据时间类型是北京时，即08时起报
&quot;move_fo_time&quot;: 0 #是否对预报的时效进行平移，12 表示将1月1日08时的36小时预报转换成1月1日20时的24小时预报后参与对比
},
&quot;SCMOC&quot;: {
&quot;dir_fo&quot;: r&quot;S:\data\grid\NWFD_SCMOC\RAIN03\YYYYMMDD\YYMMDDHH.TTT.nc&quot;, # 数据路径
&quot;hour&quot;: [8, 20, 12], #一天中的起报时长，第一个时次，最后一个时次，时次间隔
&quot;dtime&quot;:[3,48,3], #检验的时效
&quot;read_method&quot;: meteva.base.io.read_griddata_from_nc, #读取数据的函数
&quot;read_para&quot;: {}, #读取数据的函数参数
&quot;reasonable_value&quot;: [0, 1000], #合理的预报值的取值范围，超出范围观测将被过滤掉
&quot;operation&quot;: meteva.base.fun.sum_of_sta, # #预报数据读取后处理函数，sum_of_sta 表示在将逐3小时格点降水插值到站点后，需要将3小时降水累加成24小时降水
&quot;operation_para&quot;: {&quot;used_coords&quot;: [&quot;dtime&quot;], &quot;span&quot;: 24}, # #预报数据读取后处理参数
&quot;time_type&quot;: &quot;BT&quot;, # 预报数据时间类型是北京时，即08时起报
&quot;move_fo_time&quot;: 0 #是否对预报的时效进行平移，12 表示将1月1日08时的36小时预报转换成1月1日20时的24小时预报后参与对比
},
},
&quot;output_dir&quot;:None #观测站点合并数据的输出路径，设置为None时不输出收集数据的中间结果
}
</code></pre>
<h2>关于para的补充说明：</h2>
<ol>
<li>mid_method 参数常用的选项包括 mem.tase,mem.tmmsss,mem.hfmc等。如果要统计平均误差、平均绝对误差、均方根误差就选mem.tase；如果要计算二分类的指标，如TS、BIAS等就选择mem.hfmc； 如果要选相关系数就选mem.tmmsss，大部分情况下高空要素一般要计算距平相关系数，而非普通的相关系数。在meteva中计算距平相关系数有专门的ACC模块。</li>
<li>grade_list 和compare参数只有当mid_method = mem.hfmc时才需要设置，它们是二分类检验要用到的参数</li>
<li>统计形成的中间量数据文件中包含但未报仅包含begin_time 至end_time 时段内的数据，如果recover=False，它还会包含上一次运行留下来的统计数据。</li>
<li>time_step参数一般用于起报时间的间隔，例如全球模式一般是12小时间隔报一次。在调试阶段，time_step也可以设置得大一些，例如1200， 这样程序会跳跃着读取部分日期的数据进行统计，用户可以根据部分数据的统计结果来计算评分和绘图，如果检验结果合理，在重新将dtime 设置为12进行统计。避免花费很长时间完成统计后，结果是错误的，则又要重新返工，浪费大量时间。 对于非常长的时间序列，time_step 也未必一定要根据实际的预报起报时间间隔来进行设置，例如时间间隔是12小时，那将time_step 设置成60小时，也可以得到很有代表性的统计结果。此处只所以设置为60而不是48或72，是考虑到后面两种设置方法不能遍历到00和12时两种起报时间。</li>
<li>read_method对应的读取数据函数应该包含参数time,level,dtime，本模块的程序会自动循环采用不同的time,dtime,level参数调用read_method。read_para 中则不必包含time,dtime,level相关的内容，因为本程序会自动设置这些参数。</li>
</ol>
<h2>调用方法</h2>
<pre><code class="language-python"># 统计rmse，me,mae 相关的中间量
mps.middle_of_score_surface(para)</code></pre>
<pre><code>Waiting for all subprocesses done...
All subprocesses done.
中间量拼接程序运行完毕
拼接后的中间结果已输出至H:\test_data\output\mps\hfmc.h5
总耗时：66.90529608726501</code></pre>
<pre><code class="language-python">#加载中结果数据
df_mid = pd.read_hdf(para[&quot;middle_result_path&quot;])
print(df_mid)</code></pre>
<pre><code> time dtime member grade H F M C
0 2024-09-01 08:00:00 24 ECMWF 0.1 659.0 1075.0 41.0 629.0
1 2024-09-01 08:00:00 24 ECMWF 10.0 34.0 159.0 97.0 2114.0
2 2024-09-01 08:00:00 24 ECMWF 25.0 3.0 10.0 34.0 2357.0
3 2024-09-01 08:00:00 24 ECMWF 50.0 0.0 0.0 9.0 2395.0
4 2024-09-01 08:00:00 24 ECMWF 100.0 0.0 0.0 0.0 2404.0
.. ... ... ... ... ... ... ... ...
0 2024-09-03 20:00:00 48 SCMOC 0.1 674.0 221.0 318.0 1192.0
1 2024-09-03 20:00:00 48 SCMOC 10.0 90.0 141.0 116.0 2058.0
2 2024-09-03 20:00:00 48 SCMOC 25.0 10.0 41.0 35.0 2319.0
3 2024-09-03 20:00:00 48 SCMOC 50.0 0.0 4.0 10.0 2391.0
4 2024-09-03 20:00:00 48 SCMOC 100.0 0.0 0.0 0.0 2405.0
[180 rows x 8 columns]</code></pre>
<pre><code class="language-python">#根据中间量结果绘制误差曲线
result = mps.score_df(df_mid,mem.ts,g = [&quot;grade&quot;,&quot;member&quot;,&quot;dtime&quot;],plot =&quot;line&quot;,ncol = 3,sup_title = &quot;TS&quot;)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=6382dbf030757fd95e2457975b4a2a65&amp;file=file.png" alt="" /></p>
<p>上面是对比不同模式的ts的实现方式。如果需要对比模式对副高的预报性能，则可以用二分类检验指标来实现</p>
<pre><code class="language-python">#根据中间量结果绘制误差曲线
result = mps.score_df(df_mid,mem.bias,g = [&quot;grade&quot;,&quot;member&quot;,&quot;dtime&quot;],plot =&quot;line&quot;,ncol = 3,sup_title = &quot;BIAS&quot;)
</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=737836e0963c1e2953cfa45ac598208a&amp;file=file.png" alt="" /></p>
<h2>风向风速预报检验实例</h2>
<pre><code class="language-python">
def read_ob_wind(path):
speed = meteva.base.read_stadata_from_micaps1_2_8(path,column=meteva.base.m2_element_column.风速)
angle = meteva.base.read_stadata_from_micaps1_2_8(path,column=meteva.base.m2_element_column.风向)
if speed is not None and angle is not None:
speed = meteva.base.sele_by_para(speed,value=[0,100])
wind = meteva.base.speed_angle_to_wind(speed,angle)
return wind
else:
return None
para_wind= {
&quot;base_on&quot;: &quot;foTime&quot;, # 程序运行时段范围是基于起报时间
&quot;mid_method&quot;: meteva.method.na_uv, # 统计检验的总量计算函数，
&quot;grade_list&quot;: None, # 二分类检验所需的等级参数
&quot;compare&quot;: None, # 二分类检验所需的对比方式参数
&quot;begin_time&quot;:datetime.datetime(2024,9,1,8),
&quot;end_time&quot;:datetime.datetime(2024,9,3,20),
&quot;time_type&quot;:&quot;BT&quot;, # 最终检验结果呈现时，采用北京时还是世界时，UT代表世界时，BT代表北京时
&quot;time_step&quot;: 12, # 起报时间间隔
&quot;middle_result_path&quot;: r&quot;H:\test_data\output\mps\na.h5&quot;, #检验中间量数据的存储路径
&quot;station_file&quot;:meteva.base.station_国家站, #检验站点表的存储路径
&quot;defalut_value&quot;:meteva.base.IV,
&quot;interp&quot;: meteva.base.interp_gs_nearest,
&quot;recover&quot;:False, #False表示保留已经收集好的中间量，补齐未收集的数据，True表示删除已有的统计结果文件，重新开始收集
&quot;cpu&quot;:6,
&quot;ob_data&quot;:{
&quot;dir_ob&quot;: r&quot;\\10.20.22.27\rc\WIND\w2mi0\YYYYMMDDHH.000&quot;, # 实况场数据路径
&quot;hour&quot;:None,
&quot;read_method&quot;: read_ob_wind, #读取数据的函数
&quot;read_para&quot;: {}, #读取数据的函数参数
&quot;reasonable_value&quot;: None,
&quot;operation&quot;:None, #实况数据读取后处理函数
&quot;operation_para&quot;: {}, #实况数据读取后处理参数
&quot;time_type&quot;: &quot;BT&quot;, # 数据文件中的时间类型，UT代表世界时
},
&quot;fo_data&quot;:{
&quot;ECMWF&quot;: {
&quot;dir_fo&quot;: r&quot;S:\data\grid\ECMWF_HR\WIND_10M\YYYYMMDD\YYMMDDHH.TTT.nc&quot;, # 数据路径
&quot;hour&quot;: [8, 20, 12],
&quot;dtime&quot;:[12,240,12], #检验的时效
&quot;read_method&quot;: meteva.base.io.read_griddata_from_nc, #读取数据的函数
&quot;read_para&quot;: {}, #读取数据的函数参数
&quot;reasonable_value&quot;: None,
&quot;operation&quot;: None, #预报数据读取后处理函数
&quot;operation_para&quot;: {}, #预报数据读取后处理参数
&quot;time_type&quot;: &quot;BT&quot;, # 预报数据时间类型是北京时，即08时起报
&quot;move_fo_time&quot;: 0 #是否对预报的时效进行平移，12 表示将1月1日08时的36小时预报转换成1月1日20时的24小时预报后参与对比
},
&quot;SCMOC&quot;: {
&quot;dir_fo&quot;: r&quot;S:\data\grid\NWFD_SCMOC\WIND\10M_ABOVE_GROUND\YYYYMMDD\YYMMDDHH.TTT.nc&quot;,
&quot;hour&quot;: [8, 20, 12],
&quot;dtime&quot;:[12,240,12],
&quot;read_method&quot;: meteva.base.io.read_griddata_from_nc,
&quot;read_para&quot;: {},
&quot;reasonable_value&quot;: None,
&quot;operation&quot;: None,
&quot;operation_para&quot;: {},
&quot;time_type&quot;: &quot;BT&quot;,
&quot;move_fo_time&quot;: 0
},
},
&quot;output_dir&quot;:None #观测站点合并数据的输出路径，设置为None时不输出收集数据的中间结果
}</code></pre>
<pre><code class="language-python"># 统计rmse，me,mae 相关的中间量
mps.middle_of_score_surface(para_wind)</code></pre>
<pre><code>ECMWF2024年09月01日08时起报的检验中间量已存在
ECMWF2024年09月01日20时起报的检验中间量已存在
ECMWF2024年09月02日08时起报的检验中间量已存在
ECMWF2024年09月02日20时起报的检验中间量已存在
ECMWF2024年09月03日08时起报的检验中间量已存在
ECMWF2024年09月03日20时起报的检验中间量已存在
SCMOC2024年09月01日08时起报的检验中间量已存在
SCMOC2024年09月01日20时起报的检验中间量已存在
SCMOC2024年09月02日08时起报的检验中间量已存在
SCMOC2024年09月02日20时起报的检验中间量已存在
SCMOC2024年09月03日08时起报的检验中间量已存在
SCMOC2024年09月03日20时起报的检验中间量已存在
[]
未指定需要并行的参数列表
Waiting for all subprocesses done...
All subprocesses done.
中间量拼接程序运行完毕
拼接后的中间结果已输出至H:\test_data\output\mps\na.h5
总耗时：0.03599262237548828</code></pre>
<pre><code class="language-python">#加载中结果数据
df_mid = pd.read_hdf(para_wind[&quot;middle_result_path&quot;])
print(df_mid)</code></pre>
<pre><code> time dtime member T ACZ ACD ACS
0 2024-09-01 08:00:00 12 ECMWF 2407.0 322.0 717.0 1049.0
0 2024-09-01 08:00:00 12 SCMOC 2407.0 393.0 706.0 1354.0
0 2024-09-01 08:00:00 24 ECMWF 2408.0 347.0 770.0 1189.0
0 2024-09-01 08:00:00 24 SCMOC 2408.0 441.0 774.0 1482.0
0 2024-09-01 08:00:00 36 ECMWF 2408.0 282.0 748.0 975.0
.. ... ... ... ... ... ... ...
0 2024-09-03 20:00:00 216 SCMOC 2408.0 266.0 541.0 1168.0
0 2024-09-03 20:00:00 228 ECMWF 2408.0 140.0 515.0 780.0
0 2024-09-03 20:00:00 228 SCMOC 2408.0 247.0 504.0 1228.0
0 2024-09-03 20:00:00 240 ECMWF 2408.0 133.0 479.0 724.0
0 2024-09-03 20:00:00 240 SCMOC 2408.0 285.0 571.0 1169.0
[240 rows x 7 columns]</code></pre>
<pre><code class="language-python">#根据中间量结果绘制误差曲线
result = mps.score_df(df_mid,mem.acd_uv,g = [&quot;member&quot;,&quot;dtime&quot;],plot =&quot;line&quot;,ncol = 3,title = &quot;风向预报准确率&quot;)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=a11a4a0f8b81d6141769abbe850dfb94&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">#根据中间量结果绘制误差曲线
result = mps.score_df(df_mid,mem.acs_uv,g = [&quot;member&quot;,&quot;dtime&quot;],plot =&quot;line&quot;,ncol = 3,title = &quot;风速预报准确率&quot;)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=ef6ebe89c01ea2cc6a584a9fe196a2fb&amp;file=file.png" alt="" /></p>
<pre><code class="language-python"></code></pre>
<h1>并行任务下sbatch脚本</h1>
<p>在大型机上如果使用并行方式运行mps.middle_of_score程序（即参数cpu大于1）时务必采用sbatch方式提交作业，否则可能导致登录节点卡顿。
下面是提交作业的sbatch脚本，其中:</p>
<ul>
<li>参数N=1表示只采用一个计算几点（其中包含多个cpu核心）。</li>
<li>mem=64G表示计算所需的内存。如果统计中国附近区域中尺度模式一整年的数据，每天2个起报时次，每个起报时次10个预报时效，每个时效20层数据，para["step"]设置为5，mem需要设置在16G以上。如果统计全球模式1整年的数据，每天2个起报时次，每个起报时间10个预报时效，每个时效有20个层次，para["step"]参数设置为5，mem需要设置在128以上。</li>
<li>参数p normal 表示用并行方式执行。</li>
</ul>
<pre><code class="language-python">#!/bin/bash
#SBATCH --comment GRAPES
#SBATCH --wckey=105-02
#SBATCH -J GRAPES
#SBATCH -N 1
#SBATCH --mem=64G
#SBATCH --ntasks-per-node=64
#SBATCH -p normal
#SBATCH -t 72:00:00
#SBATCH -o ./out.print.%j.txt
#SBATCH -e ./out.err.%j.txt
#set runtime environment variables
unset SLURM_MEM_PER_CPU
#some shell commands
ulimit -c unlimited
ulimit -s unlimited
#export QT_DEBUG_PLUGINS=1
export QT_QPA_PLATFORM='offscreen'
python hgt_data_prepare.py
</code></pre>
<pre><code class="language-python"></code></pre>