ACC
<p>[TOC]</p>
<pre><code class="language-python">%matplotlib inline
%load_ext autoreload
%autoreload 2
import meteva.base as meb
import meteva.method as mem
import meteva.perspact as mps
import numpy as np
import pandas as pd
import datetime
import meteva</code></pre>
<h2>ACC(距平相关系数）概述</h2>
<p><font face="黑体" color=Black size=4>
<br>将预报场和观测场都扣除气候平均之后，再计算相关系数。<br></font></p>
<p>在应用该算法时须了解以下几点： </p>
<ol>
<li>气候平均场是基于1989-2008年的ERA-40，采用61天的窗口滑动平均得到。例如7曰31日08时的平均场，使用的是20年内每一年的7曰1日至8曰30日每日的08时的z500取平均得到，每个时刻的平均值使用到了20×61=1220个再分析数据。 </li>
<li>气候平均场以配置文件形式存放在meteva安装目录的resources 文件夹下，但考虑到文件较大，默认不随MetEva一起自动安装，但会在第一次运行ACC算法时自动下载到指定目录，如果自动下载失败，则需要手动下载文件 <a href="https://github.com/nmcdev/meteva/blob/master/meteva/resources/z500_cli.nc">https://github.com/nmcdev/meteva/blob/master/meteva/resources/z500_cli.nc</a><br />
并存放在MetEva安装目录的resources 文件夹。 </li>
<li>气候平均场中只包含每日00UTC和12UTC（（北京时08和20时）的结果，检验的对象也需是这些时刻 </li>
<li>气候平均场的空间分辨率是1.5°×1.5°，范围是全球，如果检验对象的网格和此不一致，则采用双线性插值方法将气候平均场插值到检验对象所在的网格上。 </li>
<li>气候平均场的数据有效位数为小数点后1位。 </li>
<li>气候平均场的单位时位势，而不是位势米，当检验对象时位势米时，需将检验对象乘以9.8</li>
</ol>
<h1>单时刻ACC计算</h1>
<p><font face="黑体" color=Blue size=3><strong>acc(grd_ob, grd_fo）</strong></font> </p>
<p>计算单个网格预报场和实况场之间的距平相关系数</p>
<table>
<thead>
<tr>
<th style="text-align: left;">参数</th>
<th style="text-align: center;">说明</th>
<th style="text-align: left;">备注</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>grd_ob</font></strong></td>
<td style="text-align: center;">网格观测数据</td>
<td style="text-align: left;"><a href="https://www.showdoc.com.cn/meteva?page_id=3975600815874861">网格数据格式</a></td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>grd_fo</font></strong></td>
<td style="text-align: center;">网格预报数据</td>
<td style="text-align: left;"><a href="https://www.showdoc.com.cn/meteva?page_id=3975600815874861">网格数据格式</a></td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: center;">实数</td>
<td style="text-align: left;"></td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">path_ec = r&quot;H:\test_data\input\mem\acc\ECMWF\z500\22050108.012.nc&quot;
grd_ec = meb.read_griddata_from_nc(path_ec) #读取预报场数据
print(grd_ec) </code></pre>
<pre><code>&lt;xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 241, lon: 281)&gt;
array([[[[[[5857.1, 5857.6, 5857.9, ..., 5870.1, 5870.1, 5869.6],
[5857.4, 5857.6, 5857.4, ..., 5869.6, 5869.9, 5869.6],
[5857.4, 5857.4, 5857.4, ..., 5869.4, 5869.9, 5870.1],
...,
[5397.9, 5400.1, 5402.4, ..., 5461.4, 5463.4, 5465.6],
[5390.1, 5392.1, 5394.1, ..., 5461.9, 5464.1, 5466.4],
[5382.4, 5384.4, 5386.4, ..., 5462.1, 5464.4, 5466.6]]]]]])
Coordinates:
* member (member) &lt;U5 'data0'
* level (level) float64 500.0
* time (time) datetime64[ns] 2022-05-01T08:00:00
* dtime (dtime) int32 12
* lat (lat) float64 0.0 0.25 0.5 0.75 1.0 ... 59.0 59.25 59.5 59.75 60.0
* lon (lon) float64 70.0 70.25 70.5 70.75 ... 139.2 139.5 139.8 140.0
Attributes:
dtime_units: hour</code></pre>
<pre><code class="language-python">path_ec0 = r&quot;H:\test_data\input\mem\acc\ECMWF\z500\22050120.000.nc&quot;
grd_ec0 = meb.read_griddata_from_nc(path_ec0) #读取实况场数据，以模式的零场来代表实况
print(grd_ec0)</code></pre>
<pre><code>&lt;xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 241, lon: 281)&gt;
array([[[[[[5856.4, 5856.4, 5856.7, ..., 5864.4, 5864.7, 5865.4],
[5856.7, 5856.7, 5856.4, ..., 5865.2, 5865.4, 5865.9],
[5857.2, 5856.7, 5856.7, ..., 5866.2, 5866.2, 5866.2],
...,
[5399.7, 5401.9, 5404.4, ..., 5460.2, 5462.4, 5464.4],
[5392.2, 5394.2, 5396.7, ..., 5460.4, 5462.7, 5464.9],
[5384.4, 5386.4, 5388.9, ..., 5460.7, 5462.9, 5465.2]]]]]])
Coordinates:
* member (member) &lt;U5 'data0'
* level (level) float64 500.0
* time (time) datetime64[ns] 2022-05-01T20:00:00
* dtime (dtime) int32 0
* lat (lat) float64 0.0 0.25 0.5 0.75 1.0 ... 59.0 59.25 59.5 59.75 60.0
* lon (lon) float64 70.0 70.25 70.5 70.75 ... 139.2 139.5 139.8 140.0
Attributes:
dtime_units: hour</code></pre>
<pre><code class="language-python">grd_ec.values *= 9.8 #根据数据结果判断预报和实况的单位都是位势米，因此乘以9.8
grd_ec0.values *= 9.8
acc1 = mem.acc(grd_ec0,grd_ec) #调用距平相关系数计算函数进行计算
print(acc1)</code></pre>
<pre><code>0.9985477061247351</code></pre>
<p>函数 mem.acc.acc 只能计算单个时刻预报和实况之间的距平相关系数，如果需要计算任意时段内的ACC则需要首先计算各时刻的相关系数中间量，再利用透视分析功能，计算任意时段的ACC。</p>
<h1>任意时段ACC计算</h1>
<h2>中间量统计</h2>
<p><font face="黑体" color=Blue size=3><strong>acc_middle_z500(para_acc）</strong></font> </p>
<p>为统计任意时段的acc准备中间量数据。该程序可以自动跳过已经统计过的起报时间和时效。</p>
<table>
<thead>
<tr>
<th style="text-align: left;">参数</th>
<th style="text-align: center;">说明</th>
<th style="text-align: left;">备注</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>para_acc</font></strong></td>
<td style="text-align: center;">待统计的观测和预报数据和检验参数描述</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: center;">无</td>
<td style="text-align: left;">以hdf文件格式将记录了每个起报时刻、每个时效、每个模式、每个区域内观测和预报的样本数、观测平均值、预报平均值、观测方差、预报方差、协方差的pandas.DataFrame格式数据输出到文件中</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">para_acc = {
&quot;begin_time&quot;: datetime.datetime(2022,5,1,8), #起始时间（预报）
&quot;end_time&quot;: datetime.datetime(2022,5,6,8), #结束时间（预报）
&quot;time_step&quot;:12, #起报时间间隔
&quot;middle_result_path&quot;: r&quot;H:\test_data\output\method\space\acc\tmmsss_z500.h5&quot;,
&quot;grid&quot; : meb.grid([70,140,0.25],[0,60,0.25]), # 检验区域
&quot;step&quot;:2, #masker间隔（单位：°），该参数只能设置为None或数据网格距的倍数，无需绘制ACC空间分布时设为None
&quot;ob_data&quot;: {
&quot;dir_ob&quot;:r&quot;H:\test_data\input\mem\acc\ECMWF\z500\YYMMDDHH.TTT.nc&quot;, # 实况场数据路径
&quot;read_method&quot;: meteva.base.io.read_griddata_from_nc, #实况数据读取函数
&quot;read_para&quot;: {}, #实况数据读取参数
&quot;time_type&quot;: &quot;BT&quot;, #实况数据时间类型
&quot;multiple&quot;: 9.8, #当文件数据的单位是位势米，通过乘以9.8转换成位势
},
&quot;fo_data&quot;: {
&quot;ECMWF&quot;: {
&quot;dir_fo&quot;:r&quot;H:\test_data\input\mem\acc\ECMWF\z500\YYMMDDHH.TTT.nc&quot;, #预报数据路径
&quot;dtime&quot;: [12, 72, 12], #检验时效
&quot;read_method&quot;: meteva.base.io.read_griddata_from_nc, #预报数据读取函数
&quot;read_para&quot;: {}, #预报数据读取参数
&quot;time_type&quot;: &quot;BT&quot;,#预报数据时间类型是北京时，即08时起报
&quot;multiple&quot;:9.8, #当文件数据的单位是位势米，通过乘以9.8转换成位势
},
&quot;NCEP&quot;: {
&quot;dir_fo&quot;: r&quot;H:\test_data\input\mem\acc\NCEP\z500\YYMMDDHH.TTT.nc&quot;, #预报数据路径
&quot;dtime&quot;: [12, 72, 12], #检验时效
&quot;read_method&quot;: meteva.base.io.read_griddata_from_nc, #预报数据读取函数
&quot;read_para&quot;: {},#预报数据读取参数
&quot;time_type&quot;: &quot;UT&quot;, #预报数据时间类型是世界时，即00时起报
&quot;multiple&quot;: 1, #当文件数据的单位是位势，则×1表示，不变
},
},
}</code></pre>
<pre><code class="language-python">mem.acc_middle_z500(para_acc) #执行中间量统计程序</code></pre>
<pre><code>ECMWF2022年05月01日08时起报的检验中间量已存在
success read from H:\test_data\input\mem\acc\ECMWF\z500\22050220.000.nc
H:\test_data\input\mem\acc\ECMWF\z500\22050120.024.nc not exists
ECMWF2022年05月02日08时起报的检验中间量已存在
success read from H:\test_data\input\mem\acc\ECMWF\z500\22050520.000.nc
H:\test_data\input\mem\acc\ECMWF\z500\22050220.072.nc not exists
ECMWF2022年05月03日08时起报的检验中间量已存在
ECMWF2022年05月03日20时起报的检验中间量已存在
ECMWF2022年05月04日08时起报的检验中间量已存在
ECMWF2022年05月04日20时起报的检验中间量已存在
ECMWF2022年05月05日08时起报的检验中间量已存在
success read from H:\test_data\input\mem\acc\ECMWF\z500\22050820.000.nc
H:\test_data\input\mem\acc\ECMWF\z500\22050520.072.nc not exists
ECMWF2022年05月06日08时起报的检验中间量已存在
NCEP2022年05月01日08时起报的检验中间量已存在
NCEP2022年05月01日20时起报的检验中间量已存在
NCEP2022年05月02日08时起报的检验中间量已存在
NCEP2022年05月02日20时起报的检验中间量已存在
NCEP2022年05月03日08时起报的检验中间量已存在
NCEP2022年05月03日20时起报的检验中间量已存在
NCEP2022年05月04日08时起报的检验中间量已存在
NCEP2022年05月04日20时起报的检验中间量已存在
NCEP2022年05月05日08时起报的检验中间量已存在
NCEP2022年05月05日20时起报的检验中间量已存在
NCEP2022年05月06日08时起报的检验中间量已存在
中间量统计程序运行完毕</code></pre>
<pre><code class="language-python">df_mid = pd.read_hdf(para_acc[&quot;middle_result_path&quot;]) #读取程序输出的中间量
print(df_mid) #查看中间量</code></pre>
<pre><code> level time dtime member id T MX \
0 500.0 2022-05-01 08:00:00 12 ECMWF 4096 63.901637 -100.572688
0 500.0 2022-05-01 08:00:00 12 ECMWF 4098 63.901637 -105.099722
0 500.0 2022-05-01 08:00:00 12 ECMWF 4100 63.901637 -108.437408
0 500.0 2022-05-01 08:00:00 12 ECMWF 4102 63.901637 -89.379083
0 500.0 2022-05-01 08:00:00 12 ECMWF 4104 63.901637 -78.855785
.. ... ... ... ... ... ... ...
0 500.0 2022-05-06 08:00:00 72 NCEP 6134 63.740954 42.878486
0 500.0 2022-05-06 08:00:00 72 NCEP 6136 63.740954 20.656238
0 500.0 2022-05-06 08:00:00 72 NCEP 6138 63.740954 -1.946539
0 500.0 2022-05-06 08:00:00 72 NCEP 6140 63.740954 3.170361
0 500.0 2022-05-06 08:00:00 72 NCEP 4094 63.901637 -121.509617
MY SX SY SXY
0 -81.292262 53.979844 55.923339 18.101133
0 -73.188211 33.639016 34.036308 9.101492
0 -70.538050 102.656192 30.606767 -5.461846
0 -91.783602 146.955801 180.611761 141.097535
0 -96.083767 129.294580 65.650976 63.215796
.. ... ... ... ...
0 77.622357 246.554394 7.959065 23.177088
0 81.058109 209.640094 9.880083 29.628266
0 84.229071 311.926179 11.967155 22.077891
0 90.052920 339.429201 10.162740 27.724071
0 -47.412352 94.031791 10.931914 5.440276
[143964 rows x 11 columns]</code></pre>
<h2>基于中间量计算ACC</h2>
<pre><code class="language-python">#绘制ACC随时效的变化图
score,gdict = mps.score_df(df_mid,mem.corr,g = [&quot;member&quot;,&quot;dtime&quot;],plot = &quot;line&quot;,
xlabel = &quot;时效（单位：小时）&quot;,ylabel = &quot;ACC&quot;,title = &quot;acc随时效的变化&quot;,grid = True,
height = 4,width = 7,sup_fontsize = 20)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=fa8e5ae8f260fba84b9aa3c9c79fd5e4&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">#选取部分样本，绘制ACC随时间的变化
score,gdict = mps.score_df(df_mid,mem.corr,s = {&quot;dtime&quot;:72},g = [&quot;member&quot;,&quot;time&quot;],plot = &quot;line&quot;,
xlabel = &quot;起报时间&quot;,ylabel = &quot;ACC&quot;,title = &quot;072H时效acc随时间的变化&quot;,grid = True,
height = 4,width = 7,sup_fontsize = 16)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=50afa1d23f2cedbb7188d5804b67f165&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">#绘制ACC的空间分布图，前提是para_acc中step的取值不为None
grd_score_list,gdict=mps.score_xy_df(df_mid,mem.corr,s = {&quot;dtime&quot;:72},g = [&quot;member&quot;],ncol = 2)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=65c8aff677e02f62346e13130402537b&amp;file=file.png" alt="" /></p>
<pre><code class="language-python"></code></pre>