高空要素统计检验
<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>在对数值模式和AI大模型数据进行检验时，经常要对长序列的多层次高空数据进行统计检验。这种长序列的统计非常耗时，并且在批量数据统计时，时常会因为各种原因导致程序中断，完成统计并非易事。改进的策略时及时的将中间结果进行保存，出错后再重复启动的时候不必从头再来。本模块的目标就是实现这样的功能。
该模块可以根据设定的实况和预报数据路径等参数，读取多个层次、多个时效、一段时间、多个模式的高空要素，生成数值型检验指标计算所需的中间量。如果程序因异常中断，重新启动后，程序能够重新加载此前已经统计过的大部分结果，对其余部分进行补齐。 之所以只能加载此前的大部分结果，是因为考虑到输出频次过高会影响效率，程序默认是完成约500个单时刻单时次单层次的中间量统计时 输出一次。</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>
<pre><code class="language-python">para = {
&quot;mid_method&quot;:mem.tase, # 统计检验的总量计算函数，
&quot;grade_list&quot;:None, #二分类检验所需的等级参数
&quot;compare&quot;:None, #二分类检验所需的对比方式参数
&quot;middle_result_path&quot;: r&quot;D:\data\ssc\veri_t\tase_t_mcv_1024.h5&quot;, #检验中间量数据的存储路径
&quot;begin_time&quot;: datetime.datetime(2023,7,20,12), #统计起始时刻， 它是预报的起报时间，而不是实况实况
&quot;end_time&quot;: datetime.datetime(2023,7,20,12), #统计结束时刻，它是预报的起报时间，而不是实况实况
&quot;time_type&quot;:&quot;UT&quot;, # 最终检验结果呈现时，采用北京时还是世界时，UT代表世界时，BT代表北京时
&quot;time_step&quot;:12, #起报时间间隔
&quot;grid&quot; : meteva.base.grid([70, 140, 0.5], [20, 60, 0.5]), # 检验区域
&quot;level_list&quot;: [1000.0,975.0,950.0, 925.0,900.0,850.0, 800.0, 750,700.0,650, 600.0,550, 500.0,450, 400.0, 350,300.0, 250.0, 200.0, 150.0, 100.0, 70.0, 50.0, 30.0, 20.0, 10.0],
#检验的层次
&quot;step&quot;:2, #masker间隔（单位：°）
&quot;recover&quot;:True, #False表示保留已经收集好的中间量，补齐未收集的数据，True表示删除已有的统计结果文件，重新开始收集
&quot;cpu&quot;:1, #cpu大于1时表示会采用python的multiprocessing进行并行的处理不同时间、时效和层次的统计任务，cpu是并行采用的核心数
&quot;ob_data&quot;: { #支持多个实况数据
&quot;ERA5&quot;:{ #实况数据名称之一
&quot;dirm&quot;: r&quot;D:\data\ssc\era5\YYYYMMDD\ERA5-plevels-YYYYMMDDHH_t.nc&quot;, # 实况场数据路径
&quot;read_method&quot;: meteva.base.read_griddata_from_nc, #读取数据的函数
&quot;read_para&quot;: {}, #读取数据的函数参数
&quot;time_type&quot;: &quot;UT&quot;, # 数据文件中的时间类型，UT代表世界时
&quot;multiple&quot;: 1, # 用于单位转换的系数，例如将比湿由kg/kg 转换成g/kg时， 该参数设置为10000
}
},
&quot;fo_data&quot;: {
&quot;FV3&quot;:{ #模式名称之一
&quot;dirm&quot;: r&quot;D:\data\ssc\mcv_ncepphys\6hour_10days\post_YYYYMMDDHH0000.ctl&quot;, # 数据路径
&quot;dtime&quot;: [0,241,6], #检验时效
&quot;read_method&quot;: meteva.base.read_griddata_from_ctl, #读取数据的函数
&quot;read_para&quot;: {&quot;value_name&quot;: &quot;t&quot;, &quot;add_block_head_tail&quot;:True, &quot;dtime_dim&quot;:&quot;tdef&quot;, &quot;endian&quot;:&quot;&gt;&quot;}, #读取数据的函数参数
&quot;time_type&quot;: &quot;UT&quot;, # 预报数据时间类型是北京时，即08时起报
&quot;multiple&quot;: 1, # 用于单位转换的系数，例如将比湿由kg/kg 转换成g/kg时， 该参数设置为10000
&quot;move_fo_time&quot;:0, #是否对预报的时效进行平移，12 表示将1月1日08时的36小时预报转换成1月1日20时的24小时预报后参与对比
&quot;veri_by&quot;:&quot;ERA5&quot;, #self表示用自己的零场作为真值进行检验,
&quot;gh&quot;:{
&quot;dirm&quot;: r&quot;D:\data\ssc\mcv_ncepphys\6hour_10days\post_YYYYMMDDHH0000.ctl&quot;, # 数据路径
&quot;read_method&quot;: meteva.base.read_griddata_from_ctl, # 读取数据的函数
&quot;read_para&quot;: {&quot;value_name&quot;: &quot;h&quot;, &quot;add_block_head_tail&quot;: True, &quot;dtime_dim&quot;: &quot;tdef&quot;, &quot;endian&quot;: &quot;&gt;&quot;},
# 读取数据的函数参数
},
},
&quot;MCV&quot;: { #模式名称之
&quot;dirm&quot;: r&quot;D:\data\ssc\mcv_grapesphys\6hour_10days\other\lnmask_ozone_mcvsst_cumict5\post_YYYYMMDDHH0000.ctl&quot;, #数据路径
&quot;dtime&quot;: [0,241,6], #检验时效
&quot;read_method&quot;: meteva.base.read_griddata_from_ctl, #读取数据的函数
&quot;read_para&quot;: {&quot;value_name&quot;: &quot;t&quot;, &quot;add_block_head_tail&quot;:True, &quot;dtime_dim&quot;:&quot;tdef&quot;, &quot;endian&quot;:&quot;&gt;&quot;}, #读取数据的函数参数
&quot;time_type&quot;: &quot;UT&quot;,#预报数据时间类型是北京时，即08时起报
&quot;multiple&quot;:1, #用于单位转换的系数，例如文件数据的单位是位势米，通过乘以9.8可转换成位势
&quot;move_fo_time&quot;: 0,
&quot;veri_by&quot;: &quot;ERA5&quot;, #指定某一个实况数据作为真值进行检验,
&quot;gh&quot;:{
&quot;dirm&quot;: r&quot;D:\data\ssc\mcv_ncepphys\6hour_10days\post_YYYYMMDDHH0000.ctl&quot;, # 数据路径
&quot;read_method&quot;: meteva.base.read_griddata_from_ctl, # 读取数据的函数
&quot;read_para&quot;: {&quot;value_name&quot;: &quot;h&quot;, &quot;add_block_head_tail&quot;: True, &quot;dtime_dim&quot;: &quot;tdef&quot;, &quot;endian&quot;: &quot;&gt;&quot;},
# 读取数据的函数参数
},
},
},
&quot;zs&quot;: { #地形高度数据读取所需信息
&quot;path&quot;: r&quot;D:\data\ssc\mcv_ncepphys\6hour_10days\post_20230720120000.ctl&quot;, # 数据路径
&quot;read_method&quot;: meteva.base.read_griddata_from_ctl, # 读取数据的函数
&quot;read_para&quot;: {&quot;value_name&quot;: &quot;zs&quot;, &quot;add_block_head_tail&quot;: True, &quot;dtime_dim&quot;: &quot;tdef&quot;, &quot;endian&quot;: &quot;&gt;&quot;,&quot;dtime&quot;:0}, # 读取数据的函数参数
}
}
</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>
<li>veri_by 用于指定检验实况，每个模式可用选择相同的实况，也可以选择不同的实况。</li>
<li>para中zs关键词用于读取地形高度数据，每个模式的gh 参数用来读取每一层的位势高度数据，当位势高度低于地形高度时，对应网格点的数据样本会被设置成nan，并在最终的检验样本中被剔除。如果不需要实用地形剔除的功能，则将zs的取值设置为None即可。</li>
</ol>
<h2>调用方法</h2>
<pre><code class="language-python"># 统计rmse，me,mae 相关的中间量
mps.middle_of_score(para)</code></pre>
<pre><code>success read from D:\data\ssc\era5\20230720\ERA5-plevels-2023072012_t.nc
success read data with D:\data\ssc\mcv_ncepphys\6hour_10days\post_20230720120000.ctl
success read data with D:\data\ssc\mcv_ncepphys\6hour_10days\post_20230720120000.ctl
success read data with D:\data\ssc\mcv_grapesphys\6hour_10days\other\lnmask_ozone_mcvsst_cumict5\post_20230720120000.ctl
....
success read from D:\data\ssc\era5\20230730\ERA5-plevels-2023073012_t.nc
success read data with D:\data\ssc\mcv_ncepphys\6hour_10days\post_20230720120000.ctl
success read data with D:\data\ssc\mcv_ncepphys\6hour_10days\post_20230720120000.ctl
success read data with D:\data\ssc\mcv_grapesphys\6hour_10days\other\lnmask_ozone_mcvsst_cumict5\post_20230720120000.ctl
success read data with D:\data\ssc\mcv_ncepphys\6hour_10days\post_20230720120000.ctl
中间量统计程序运行完毕
中间结果已输出至D:\data\ssc\veri_t\tase_t_mcv_1024.h5
单进程耗时：931.0928151607513
检验计算耗时：45.604432344436646
读取ERA5耗时：806.7255029678345
读取FV3耗时：40.8194842338562
读取MCV耗时：29.385673761367798
总耗时：931.6759107112885</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> level time dtime member id T E \
0 10.0 2023-07-20 12:00:00 0 FV3 20070 0.939693 -1.038601
1 10.0 2023-07-20 12:00:00 0 FV3 20072 3.758770 -3.005474
2 10.0 2023-07-20 12:00:00 0 FV3 20074 3.758770 -2.884430
3 10.0 2023-07-20 12:00:00 0 FV3 20076 3.758770 6.605784
4 10.0 2023-07-20 12:00:00 0 FV3 20078 3.758770 1.503090
.. ... ... ... ... ... ... ...
179 1000.0 2023-07-20 12:00:00 240 MCV 60072 0.500000 2.608284
180 1000.0 2023-07-20 12:00:00 240 MCV 60078 1.507538 2.711890
181 1000.0 2023-07-20 12:00:00 240 MCV 60080 6.090306 5.337546
182 1000.0 2023-07-20 12:00:00 240 MCV 60082 6.657801 10.963189
183 1000.0 2023-07-20 12:00:00 240 MCV 60084 2.082534 1.287826
A S M
0 1.038601 1.147919 216.234548
1 4.630156 11.051691 859.313193
2 5.575627 11.684649 854.057493
3 6.605784 11.803160 851.092762
4 1.901097 1.233169 862.144487
.. ... ... ...
179 2.608284 13.606290 147.014000
180 2.740965 7.665755 449.925190
181 5.418934 6.935573 1823.526962
182 10.963189 20.224089 1983.339122
183 1.823327 2.479527 617.696259
[1446086 rows x 10 columns]</code></pre>
<pre><code class="language-python">#根据中间量结果绘制误差曲线
result = mps.score_df(df_mid,mem.rmse,s = {&quot;level&quot;:[850,300]},g = [&quot;level&quot;,&quot;member&quot;,&quot;dtime&quot;],plot =&quot;line&quot;,ncol = 2,sup_title = &quot;RMSE&quot;)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=9f53d4b50f221cdb5397731a4f978405&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">result = mps.score_xy_df(df_mid,mem.me,s = {&quot;level&quot;:[850]},g = &quot;member&quot;,ncol = 2)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=c8ec5fbbd9736d5d61dcd6189bbaba90&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">result = mps.score_yz_df(df_mid,mem.me,s={&quot;lon&quot;:90},g = &quot;member&quot;,ncol = 2)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=c7ac1056fd500c6080398feddeb10e7f&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">result = mps.score_xz_df(df_mid,mem.me,s = {&quot;lat&quot;:30},g = &quot;member&quot;,ncol = 2)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=c3ac6f10d15bc0775eb4284757095a3e&amp;file=file.png" alt="" /></p>
<p>上面是对比不同模式的rmse的实现方式。如果需要对比模式对副高的预报性能，则可以用二分类检验指标来实现</p>
<pre><code class="language-python"># 设置统计参数
para_hfmc = {
&quot;mid_method&quot;:mem.hfmc, # 统计检验的总量计算函数，
&quot;grade_list&quot;:[5880,5920], #二分类检验所需的等级参数
&quot;compare&quot;:None, #二分类检验所需的对比方式参数
&quot;middle_result_path&quot;: r&quot;H:\test_data\output\mps\hfmc_500.h5&quot;, #检验中间量数据的存储路径
&quot;begin_time&quot;: datetime.datetime(2018,8,1,0), #统计起始时刻， 它是预报的起报时间，而不是实况实况
&quot;end_time&quot;: datetime.datetime(2018,8,2,0), #统计结束时刻，它是预报的起报时间，而不是实况实况
&quot;time_type&quot;:&quot;UT&quot;, # 最终检验结果呈现时，采用北京时还是世界时，UT代表世界时，BT代表北京时
&quot;time_step&quot;:12, #起报时间间隔
&quot;grid&quot; : meteva.base.grid([70,140,0.25],[10,60,0.25]), # 检验区域
&quot;level_list&quot;: [500], #检验的层次
&quot;step&quot;:5, #masker间隔（单位：°）
&quot;recover&quot;:False, #False表示保留已经收集好的中间量，补齐未收集的数据，True表示删除已有的统计结果文件，重新开始收集
&quot;cpu&quot;:64, #cpu大于1时表示会采用python的multiprocessing进行并行的处理不同时间、时效和层次的统计任务，cpu是并行采用的核心数
&quot;ob_data&quot;: { #支持多个实况数据
&quot;CRA40&quot;:{ #实况数据名称之一
&quot;dirm&quot;: r&quot;\\10.28.16.234\data2\AI\CRA40\2018\Z\LLL\YYYYMMDDHH.nc&quot;, # 实况场数据路径
&quot;read_method&quot;: meteva.base.read_griddata_from_nc, #读取数据的函数
&quot;read_para&quot;: {}, #读取数据的函数参数
&quot;time_type&quot;: &quot;UT&quot;, # 数据文件中的时间类型，UT代表世界时
&quot;multiple&quot;: 1, # 用于单位转换的系数，例如将比湿由kg/kg 转换成g/kg时， 该参数设置为10000
}
},
&quot;fo_data&quot;: {
&quot;ECMWF&quot;:{ #模式名称之一
&quot;dirm&quot;: r&quot;\\10.28.16.234\data2\AI\ECMWF\grib\YYMMDD\YYMMDDHH.TTT.grib&quot;, # 数据路径
&quot;dtime&quot;: [12, 240, 12], #检验的时效
&quot;read_method&quot;: meteva.base.read_griddata_from_grib, #读取数据的函数
&quot;read_para&quot;: {&quot;level_type&quot;:&quot;isobaricInhPa&quot;,&quot;value_name&quot;:&quot;gh&quot;}, #读取数据的函数参数
&quot;time_type&quot;: &quot;UT&quot;, # 预报数据时间类型是北京时，即08时起报
&quot;multiple&quot;: 1, # 用于单位转换的系数，例如将比湿由kg/kg 转换成g/kg时， 该参数设置为10000
&quot;move_fo_time&quot;:12, #是否对预报的时效进行平移，12 表示将1月1日08时的36小时预报转换成1月1日20时的24小时预报后参与对比
&quot;veri_by&quot;:&quot;self&quot; #self表示用自己的零场作为真值进行检验
},
&quot;NMC1&quot;: { #模式名称之
&quot;dirm&quot;: r&quot;\\10.28.16.177\NMC_Tsinghua_reforecast\cra40\YYYYMMDDHH.nc&quot;, #数据路径
&quot;dtime&quot;: [12, 240, 12], #检验时效
&quot;read_method&quot;: meteva.base.read_griddata_from_nc, #读取数据的函数
&quot;read_para&quot;: {&quot;value_name&quot;: &quot;GPH&quot;}, #读取数据的函数参数
&quot;time_type&quot;: &quot;UT&quot;,#预报数据时间类型是北京时，即08时起报
&quot;multiple&quot;:1, #用于单位转换的系数，例如文件数据的单位是位势米，通过乘以9.8可转换成位势
&quot;move_fo_time&quot;: 0,
&quot;veri_by&quot;: &quot;CRA40&quot; #指定某一个实况数据作为真值进行检验
},
},
}
# 统计rmse，me,mae 相关的中间量
mps.middle_of_score(para_hfmc)</code></pre>
<pre><code>NMC12018年08月01日00时起报的检验中间量已存在
NMC12018年08月01日12时起报的检验中间量已存在
NMC12018年08月02日00时起报的检验中间量已存在
Waiting for all subprocesses done...
All subprocesses done.
中间量拼接程序运行完毕
拼接后的中间结果已输出至H:\test_data\output\mps\hfmc_500.h5
总耗时：8.337354898452759</code></pre>
<pre><code class="language-python">#加载中结果数据
df_mid = pd.read_hdf(para_hfmc[&quot;middle_result_path&quot;])
print(df_mid)</code></pre>
<pre><code> level time dtime member id grade H F M C
0 500 2018-08-01 12 ECMWF 30080 5880 0.0 0.0 0.0 1600.0
1 500 2018-08-01 12 ECMWF 30080 5920 0.0 0.0 0.0 1600.0
0 500 2018-08-01 12 ECMWF 20100 5880 0.0 0.0 0.0 1600.0
1 500 2018-08-01 12 ECMWF 20100 5920 0.0 0.0 0.0 1600.0
0 500 2018-08-01 12 ECMWF 40070 5880 0.0 0.0 0.0 40.0
.. ... ... ... ... ... ... ... ... ... ...
1 500 2018-08-02 240 NMC1 30070 5920 0.0 0.0 0.0 40.0
0 500 2018-08-02 240 NMC1 20090 5880 0.0 0.0 0.0 1600.0
1 500 2018-08-02 240 NMC1 20090 5920 0.0 0.0 0.0 1600.0
0 500 2018-08-02 240 NMC1 10110 5880 0.0 0.0 0.0 40.0
1 500 2018-08-02 240 NMC1 10110 5920 0.0 0.0 0.0 40.0
[10944 rows x 10 columns]</code></pre>
<pre><code class="language-python">#根据实况和预报的频次（累计格点数）的对比
result = mps.score_df(df_mid,mem.ob_fo_hc,g = [&quot;grade&quot;,&quot;member&quot;,&quot;dtime&quot;],plot =&quot;line&quot;,ncol = 2,sup_title = &quot;RMSE&quot;)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=509f4225ac3e7d4eaa37b621b588fbd2&amp;file=file.png" alt="" /></p>
<p>如果数据时段很长，实况的副高累计格点数随时效基本是不变化的，但上面的数据集只包含一小段时间的数据，因此实况的副高累计格点数随时效变化非常剧烈。从上面的结果还可以看出，两种模式对副高面积（格点数）的预报都是偏低的。</p>
<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=&#039;offscreen&#039;
python hgt_data_prepare.py
</code></pre>
<pre><code class="language-python"></code></pre>