短时强降水检验
<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.product as mpd
import numpy as np
import pandas as pd
import datetime</code></pre>
<p><strong> 为了说明短时强降水预报检验的程序功能，以下先将检验所需的数逐一进行准备</strong> </p>
<pre><code class="language-python"># 首先读取一个站点列表文件，在具体的业务中，需根据某种具体的检验规范来确定相应文件
# 此处只是为说明算法逻辑因此，仅以一个包含2411个国家站点文件为例进行叙述
station_file = r"H:\test_data\station_pta.txt" # 指定站点文件，micaps3格式
station = meb.read_station(station_file) #读取站点文件
#将站点数据列设置为999999，由此在从具体的数据文件读取数时，数据文件中缺省的站点仍然会以999999的形式存在
station.iloc[:,-1] = meb.IV
print(station) </code></pre>
<pre><code> level time dtime id lon lat data0
0 0 2099-01-01 08:00:00 0 50136 122.52 52.97 999999
1 0 2099-01-01 08:00:00 0 50137 122.37 53.47 999999
2 0 2099-01-01 08:00:00 0 50246 124.72 52.35 999999
3 0 2099-01-01 08:00:00 0 50247 123.57 52.03 999999
4 0 2099-01-01 08:00:00 0 50349 124.40 51.67 999999
... ... ... ... ... ... ... ...
2406 0 2099-01-01 08:00:00 0 59945 109.70 18.65 999999
2407 0 2099-01-01 08:00:00 0 59948 109.58 18.22 999999
2408 0 2099-01-01 08:00:00 0 59951 110.33 18.80 999999
2409 0 2099-01-01 08:00:00 0 59954 110.03 18.55 999999
2410 0 2099-01-01 08:00:00 0 59981 112.33 16.83 999999
[2411 rows x 7 columns]</code></pre>
<pre><code class="language-python">#因为在短时强降水检验规范中，不同省份对应的降水阈值标准不一样， 在程序中为每个站点设置不同的阈值
#为此首先需要区分每个站点具体属于什么省份，这个数据可以字典形式表示。
#在具体的业务中，需根据某种具体的检验规范来确定相应站号省份， 此处仅以2411个国家站作为示例
dict0 = meb.station_id_name_dict
dict1 = {}
id_list = station["id"].values.tolist()
for id in dict0.keys():
if id in id_list:
str0 = dict0[id]
strss = str0.split("_")
str1 = strss[0]
dict1[id] = str1
print(dict1)</code></pre>
<pre><code>{50136: '黑龙江', 50137: '黑龙江', 50246: '黑龙江', 50247: '黑龙江', 50349: '黑龙江', 50353: '黑龙江', 50425: '内蒙古', 50431: '内蒙古', 50434: '内蒙古', 50442: '黑龙江', 50445: '内蒙古', 50468: '黑龙江', 50514: '内蒙古', 50524: '内蒙古',
。。。,
59847: '海南', 59848: '海南', 59849: '海南', 59851: '海南', 59854: '海南', 59855: '海南', 59856: '海南', 59940: '海南', 59941: '海南', 59945: '海南', 59948: '海南', 59951: '海南', 59954: '海南', 59981: '海南'}</code></pre>
<p>我们算法的测试数据包括：</p>
<ol>
<li>2020年8月18日08时 至19日08时的逐小时降水观测；</li>
<li>两份2020年8月18日08时和20时起报的未来逐12小时预报，分别代表中央台和省台的预报。 </li>
</ol>
<pre><code class="language-python">#设置数据的时间段
time_start = datetime.datetime(2020,8,18,8,0)
time_end = datetime.datetime(2020,8,19,8,0)
###读取2020年8月18日08时 至19日08时的逐小时降水观测
dir_ob = r"Y:\qpe\autorun\data_all\001h\YYYYMMDDHH.000.001h"
sta_list = []
time0 = time_start
while time0 <= time_end:
path = meb.get_path(dir_ob,time0)
sta = meb.read_stadata_from_micaps3(path,station = station,time = time0,dtime = 0,level = 0,data_name = "ob")
sta_list.append(sta)
time0 += datetime.timedelta(hours = 1)
ob_sta_all = pd.concat(sta_list,axis = 0) #数据拼接
print(ob_sta_all)</code></pre>
<pre><code> level time dtime id lon lat ob
0 0 2020-08-18 08:00:00 0 50136 122.52 52.97 0.0
1 0 2020-08-18 08:00:00 0 50137 122.37 53.47 0.0
2 0 2020-08-18 08:00:00 0 50246 124.72 52.35 0.0
3 0 2020-08-18 08:00:00 0 50247 123.57 52.03 0.2
4 0 2020-08-18 08:00:00 0 50349 124.40 51.67 1.0
... ... ... ... ... ... ... ...
2406 0 2020-08-19 08:00:00 0 59945 109.70 18.65 0.0
2407 0 2020-08-19 08:00:00 0 59948 109.58 18.22 0.0
2408 0 2020-08-19 08:00:00 0 59951 110.33 18.80 0.0
2409 0 2020-08-19 08:00:00 0 59954 110.03 18.55 0.0
2410 0 2020-08-19 08:00:00 0 59981 112.33 16.83 0.5
[60275 rows x 7 columns]</code></pre>
<pre><code class="language-python"># 以某种网格预报数据代表中央台的预报
dir1 = r"O:\data\grid\NWFD_SCMOC_1H\RAIN01\YYYYMMDD\YYMMDDHH.TTT.nc"
sta_list = []
time0 = time_start
while time0 < time_end:
for dh in range(1, 13):
path = meb.get_path(dir1, time0, dh)
grd = meb.read_griddata_from_nc(path,time=time0,dtime = dh,level = 0,data_name="nmc")
if grd is not None:
sta = meb.interp_gs_linear(grd, station)
sta_list.append(sta)
time0 += datetime.timedelta(hours=12)
sta_nmc = pd.concat(sta_list, axis=0) # 数据拼接
#设置阈值，将预报值大于阈值的站点设置为1，否则设置为0
threshold = 2 # 考虑到当天数据降水量大于20mm/h的并不多，为了展示方便此处故意将阈值随便设置为一个较低的值
sta_nmc.iloc[:,-1] = 0 + (sta_nmc.iloc[:,-1] >= threshold) # 将预报转换成0-1形式
print(sta_nmc)</code></pre>
<pre><code> level time dtime id lon lat nmc
0 0 2020-08-18 08:00:00 1 50136 122.52 52.97 0
1 0 2020-08-18 08:00:00 1 50137 122.37 53.47 0
2 0 2020-08-18 08:00:00 1 50246 124.72 52.35 0
3 0 2020-08-18 08:00:00 1 50247 123.57 52.03 0
4 0 2020-08-18 08:00:00 1 50349 124.40 51.67 0
... ... ... ... ... ... ... ...
2406 0 2020-08-18 20:00:00 12 59945 109.70 18.65 0
2407 0 2020-08-18 20:00:00 12 59948 109.58 18.22 0
2408 0 2020-08-18 20:00:00 12 59951 110.33 18.80 0
2409 0 2020-08-18 20:00:00 12 59954 110.03 18.55 0
2410 0 2020-08-18 20:00:00 12 59981 112.33 16.83 0
[57864 rows x 7 columns]</code></pre>
<pre><code class="language-python"># 以另一种网格预报数据代表省台的预报
dir1 = r"O:\data\grid\GMOSRR\ROLLING_UPDATE\RAIN01\YYYYMMDD\YYMMDDHH.TTT.nc"
sta_list = []
time0 = time_start
while time0 < time_end:
for dh in range(1, 13):
path = meb.get_path(dir1, time0, dh)
grd = meb.read_griddata_from_nc(path,time=time0,dtime = dh,level = 0,data_name="nmc")
if grd is not None:
sta = meb.interp_gs_linear(grd, station)
sta_list.append(sta)
time0 += datetime.timedelta(hours=12)
sta_pro = pd.concat(sta_list, axis=0) # 数据拼接
#设置阈值，将预报值大于阈值的站点设置为1，否则设置为0
threshold = 3 # 考虑到当天数据降水量大于20mm/h的并不多，为了展示方便此处故意将阈值随便设置为一个较低的值
sta_pro.iloc[:,-1] = 0 + (sta_pro.iloc[:,-1] >= threshold) # 将预报转换成0-1形式</code></pre>
<pre><code class="language-python">#读取省会城市单点预报，原则上从其它文件中读取，但为了简化问题，此处从网格预报插值到站点的结果中直接选取
pcapital_id = [54511] #此处选取一个54511站作为示例
nmc_pcapital = meb.sele_by_para(sta_nmc,id = pcapital_id)
print(nmc_pcapital)</code></pre>
<pre><code> level time dtime id lon lat nmc
779 0 2020-08-18 08:00:00 1 54511 116.47 39.8 0
779 0 2020-08-18 08:00:00 2 54511 116.47 39.8 0
779 0 2020-08-18 08:00:00 3 54511 116.47 39.8 0
779 0 2020-08-18 08:00:00 4 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 5 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 6 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 7 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 8 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 9 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 10 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 11 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 12 54511 116.47 39.8 1
779 0 2020-08-18 20:00:00 1 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 2 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 3 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 4 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 5 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 6 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 7 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 8 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 9 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 10 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 11 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 12 54511 116.47 39.8 0</code></pre>
<pre><code class="language-python">pro_pcapital = meb.sele_by_para(sta_pro,id = pcapital_id)
print(pro_pcapital)</code></pre>
<pre><code> level time dtime id lon lat nmc
779 0 2020-08-18 08:00:00 1 54511 116.47 39.8 0
779 0 2020-08-18 08:00:00 2 54511 116.47 39.8 0
779 0 2020-08-18 08:00:00 3 54511 116.47 39.8 0
779 0 2020-08-18 08:00:00 4 54511 116.47 39.8 0
779 0 2020-08-18 08:00:00 5 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 6 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 7 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 8 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 9 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 10 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 11 54511 116.47 39.8 1
779 0 2020-08-18 08:00:00 12 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 1 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 2 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 3 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 4 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 5 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 6 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 7 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 8 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 9 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 10 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 11 54511 116.47 39.8 0
779 0 2020-08-18 20:00:00 12 54511 116.47 39.8 0</code></pre>
<h1>短时强降水检验模块（2020_1版）</h1>
<p><font face="黑体" color=blue size = 5><strong>mpd.short_term_heavy_rainfall.edition_2020_1(sta_ob,sta_fo_list,pcapital_fo_list, id_provice_dict)</strong></font><br />
基于对2020年3月制定的《智能预报技术方法竞赛检验方案》中关于短时强降水检验部分的理解开发的检验功能，在使用前请阅读该模块前请先仔细阅读其有关<a href="https://www.showdoc.com.cn/meteva?page_id=5229500776612325">概述</a>。<br />
在算法上它调用到了<a href="https://www.showdoc.com.cn/meteva?page_id=5225225007079073">邻域检验</a>方面的功能。具体代码请参考GitHub上 <a href="https://github.com/nmcdev/meteva/blob/master/meteva/product/regulation/short_term_heavy_rainfall.py">nmcdev/meteva项目</a> </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>sta_ob</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.com.cn/meteva?page_id=3975600580125986">站点数据</a>格式的观测数据，内容为每个站点的降水量，其中可以包含多个时刻的站点数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>sta_fo_list</font></strong></td>
<td style="text-align: left;">预报数据列表，其中第0个元素为中央台的预报插值到站点得到的站点数据的合集，第1个元素为省台预报对应结果.其中可以包含多个起报时间和预报时效的数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>pcapital_fo_list</font></strong></td>
<td style="text-align: left;">省会单点预报数据列表，其中第0个元素为中央台的预报插值到站点得到的站点数据的合集，第1个元素为省台预报对应结果.其中可以包含多个起报时间和预报时效的数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>id_province_dict</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;">包含4个元素的元组,分别是网格预报命中率技巧，网格预报空报率技巧，省会单站命中率技巧，省会单站空报率技巧，其中每个元素为长度12的列表，分别对应12个预报时效</td>
</tr>
</tbody>
</table>
<p><strong>调用示例：</strong> </p>
<pre><code class="language-python">#基于前面的设置的测试数据，测试该函数的调用方式并展示结果
spo, sfa, spo_pcapital, sfa_pcapital = mpd.short_term_heavy_rainfall.edition_2020_1(ob_sta_all,
[sta_nmc,sta_pro],
[nmc_pcapital,pro_pcapital],id_province_dict=dict1)
print("网格预报命中率技巧：")
print(spo)
print()
print("网格预报空报率技巧：")
print(sfa)
print()
print("省会单站预报命中率技巧：")
print(spo_pcapital)
print()
print("省会单站空报率技巧：")
print(sfa_pcapital)</code></pre>
<pre><code>网格预报命中率技巧：
[-0.063 -0.655 -0.394 -0.198 -0.176 -0.201 -0.169 -0.197 -0.176 -0.201
-0.117 -0.107]
网格预报空报率技巧：
[ 0.15 0.311 0.171 0.201 0.093 0.04 0.015 0.004 -0.002 -0.032
-0.011 -0.03 ]
省会单站预报命中率技巧：
[-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.]
省会单站空报率技巧：
[ 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. -1.]</code></pre>
<pre><code class="language-python"></code></pre>