图片型检验产品
<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 datetime
import copy
import matplotlib.pyplot as plt
import pandas as pd</code></pre>
<pre><code>The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload</code></pre>
<p> 在做预报检验时，数据提取和检验计算的通常是交织在一起的，在本程序中的meteva.method模块中也是如此，在meteva.method中至少计算检验中间量是和数据提取是交织在一起的。有时，我们面对一批预报和观测数据要开展检验工作时，事先并不知道后期要计算什么检验指标，这时需要反复尝试。由于数据提取和检验技术的代码交织在一起，意味着每次增删或更改检验指标时，数据提取的部分也需要重新运行一遍。另外上述问题也意味着想把数据提取部分的代码做进一步的封装难以实现或没有意义。<br />
本程序库则为上述问题提供了解决方案，基于meteva.base和meteva.method中的函数，可以将检验过程分割以下三个部分：<br />
1.数据收集(文件读取，<a href="https://www.showdoc.cc/meteva?page_id=3975602867727263">拼接合并</a>）<br />
2.<a href="https://www.showdoc.cc/meteva?page_id=3975604785954540">数据选取</a><br />
3.检验计算<br />
这为检验程序的进一步模块化提供了基础。以下结合具体的例子进行说明 </p>
<p><a href="https://www.showdoc.cc/nmc?page_id=3831227192066999">本模块测试数据集简介和更详细的数据收集代码说明</a> </p>
<pre><code class="language-python">###################以下开始为数据收集部分的程序
dir_ob = r"H:\test_data\input\mpd\ob\temp_2m\BTYYMMDDHH.000"
dir_ec = r"H:\test_data\input\mpd\ec\temp_2m\YYMMDD\BTYYMMDDHH.TTT"
dir_grapes = r"H:\test_data\input\mpd\grapes\temp_2m\YYMMDD\BTYYMMDDHH.TTT"
time0 = datetime.datetime(2019,1,1,2,0)
path = meb.get_path(dir_ob,time0)
station = meb.read_stadata_from_micaps3(path)
station.iloc[:,-1] = meb.IV
##读取收集观测数据
time_end = datetime.datetime(2020,1,1,2,0)
ob_sta_list = []
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")
ob_sta_list.append(sta)
time0 += datetime.timedelta(hours = 3)
ob_sta_all = pd.concat(ob_sta_list,axis = 0) #数据拼接
##读取收集预报数据
ec_sta_list =[]
grapes_sta_list = []
time0 = datetime.datetime(2019,1,1,8,0)
time_end = datetime.datetime(2020,1,1,8,0)
while time0 < time_end:
for dh in range(0,73,3):
#读取ec预报数据
path = meb.get_path(dir_ec,time0,dh)
grd = meb.read_griddata_from_micaps4(path)
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"])
ec_sta_list.append(sta)
#读取grapes预报数据
path = meb.get_path(dir_grapes,time0,dh)
grd = meb.read_griddata_from_micaps4(path)
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"])
grapes_sta_list.append(sta)
time0 += datetime.timedelta(hours = 12)
ec_sta_all = pd.concat(ec_sta_list,axis = 0) #数据拼接
grapes_sta_all = pd.concat(grapes_sta_list,axis = 0) #数据拼接
#数据匹配合并
sta_all = meb.combine_on_obTime(ob_sta_all,[ec_sta_all,grapes_sta_all])
sta_all = meb.not_IV(sta_all) #删除包含缺省值的样本
sta_all = meb.not_equal_to(sta_all,9999)
###################以上开始为数据收集部分的程序</code></pre>
<h1>图片型检验产品制作</h1>
<p><strong><font face="黑体" color=blue size=3>plot(sta_ob_and_fos,method,s = None,g = None,gll = None,save_dir = None,save_path = None,title =None,</strong>kwargs)</font>**<br />
根据输入的站点数据和检验方法，生成图片型检验产品，并自动生成后缀名为.png的文件批量输出到指定目录 </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;"><font face="黑体" color=blue size=5><strong>sta_ob_and_fos</strong></font></td>
<td style="text-align: left;">实况和预报合并对齐后的数据，形式为站点数据格式如上述例子中的sta_all</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5><strong>method</strong></font></td>
<td style="text-align: left;">method中的各类图片型的函数名称，例如mem.performance,mem.scatter_regress,mem.pdf_plot,mem.box_plot_continue,mem.frequency_histogram, mem.rank_histogram,mem.comprehensive_probability,mem.scatter_uv,mem.scatter_uv_error,mem.statistic_uv等函数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>s</strong></td>
<td style="text-align: left;">用于选择数据样本的字典参数，具体的参数说明可参见meb.sele_by_dict中的<a href="https://www.showdoc.cc/meteva?page_id=3975604785954540"><font face="黑体" color=red size=5>s</font></a>参数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>g</strong></td>
<td style="text-align: left;">用于分组检验的参数，具体用法可参见meb.group中的<a href="https://www.showdoc.cc/meteva?page_id=4071849185300418"><font face="黑体" color=red size=5>g</font></a>参数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>gll</strong></td>
<td style="text-align: left;">用于分组检验的参数，具体用法可参见meb.group中的<a href="https://www.showdoc.cc/meteva?page_id=4071849185300418"><font face="黑体" color=red size=5>gll</font></a>参数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>save_dir</strong></td>
<td style="text-align: left;">生成的表格文件的保存目录</td>
</tr>
<tr>
<td style="text-align: left;"><strong>save_path</strong></td>
<td style="text-align: left;">指定图片输出的文件路径，当批量生成多张图片时save_path为包含所有图片的输出路径的列表</td>
</tr>
<tr>
<td style="text-align: left;"><strong>title</strong></td>
<td style="text-align: left;">图片标题的确定分为1、全自动，2、半自动，3手动三种方式。<br>1、全自动：不设置title参数，系统自动采用title的缺省值+自动补齐的其它信息来确定每一幅图的标题；<br> 2、半自动：title 为字符串类型，系统会采用title + 自动补齐的其它信息来确定每一幅图的标题；<br>3、手动：title为一个包含多个字符串的列表，且列表的长度必须和要绘制的图的数量一致，每一幅图会依次采用列表中的字符串作为标题</td>
</tr>
<tr>
<td style="text-align: left;">**kwargs</td>
<td style="text-align: left;">检验方法 meteva.method 中的可选参数，具体用法参见下面的示例</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">mpd.plot(sta_all,mem.scatter_regress,s = {"id":54398}) #选取54398站进行检验，group_by = None</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=478c3b6aa85245ffbc00271a6981529c" alt="" /></p>
<p>从上面的整体检验可以看出grapes模式在54398站上的温度预报存在系统性偏低的情况。如果想进一步分析不同时刻起报的预报是否都存在这样的问题，可以进一步设置group_by = "hour"来检验，如下：</p>
<pre><code class="language-python">mpd.plot(sta_all,mem.scatter_regress) #输入需要检验的数据和方法，即可自动绘制相应的检验图形产品</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=2e920d7b8e1f3d35455567e7796dc9a6" alt="" /></p>
<pre><code class="language-python">mpd.plot(sta_all,mem.pdf_plot) #切换一下检验方法名，即可以完成新的检验产品的制作</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=c79a1b446c0a7a2189b24f9a9ec25d27" alt="" /></p>
<pre><code class="language-python">mpd.plot(sta_all,mem.box_plot_continue) #制作观测预报频率对比箱须图</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/2c9ca69aadf420557cd73b1ed722d8f0?showdoc=.jpg" alt="" /></p>
<pre><code class="language-python">mpd.plot(sta_all,mem.performance,grade_list = [0,10,15,20,25,30,35])
#制作二分类预报的综合表现图，grade_list 是mem.performance函数中的参数</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/150091dd6e03b6df8785bfde7cc8a6eb?showdoc=.jpg" alt="" /></p>
<pre><code class="language-python">mpd.plot(sta_all,mem.frequency_histogram,grade_list = np.arange(-20,40,5).tolist()) #制作预报观测频率</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/072b2f3a1688a2fb9560924b548eba01?showdoc=.jpg" alt="" /></p>
<pre><code class="language-python">mpd.plot(sta_all,mem.frequency_histogram,s = {"member":["ob","grapes"]},g = "hour",
grade_list = np.arange(-20,40,5).tolist()) </code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/073587b0e832ea4384f392916cd24288?showdoc=.jpg" alt="" /></p>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/baeba889436e92c756bfb59a1cb0b78a?showdoc=.jpg" alt="" /></p>
<pre><code>[8, 20]</code></pre>
<pre><code class="language-python">mpd.plot(sta_all,mem.frequency_histogram,s = {"member":["ob","grapes"]},g = "hour",
grade_list = np.arange(-20,40,5).tolist(),save_dir = r"H:\test_data\output\mpd\program\plot",show = True,
title = "预报观测频率对比") #半自动的方式设置title</code></pre>
<pre><code>检验结果已以图片形式保存至H:/test_data/output/mpd/program/plot/frequency_histogram__hour8_.png</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/6b213d541bade6f37bee998aab818e36?showdoc=.jpg" alt="" /></p>
<pre><code>检验结果已以图片形式保存至H:/test_data/output/mpd/program/plot/frequency_histogram__hour20_.png</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/9ade753b4c7c5cc8e1bd9297d601c91a?showdoc=.jpg" alt="" /></p>
<pre><code>[8, 20]</code></pre>
<pre><code class="language-python">mpd.plot(sta_all,mem.frequency_histogram,s = {"member":["ob","grapes"]},g = "hour",
grade_list = np.arange(-20,40,5).tolist(),
save_path = [r"H:\test_data\output\mpd\program\plot\08.png",r"H:\test_data\output\mpd\program\plot\20.png"],
show = True,title = ["预报观测频率对比(08时起报)","预报观测频率对比(20时起报)"]) #手动设置每一幅图片的title。 </code></pre>
<pre><code>检验结果已以图片形式保存至H:\test_data\output\mpd\program\plot\08.png</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/65af454a1c0bc037280c00354965aa86?showdoc=.jpg" alt="" /></p>
<pre><code>检验结果已以图片形式保存至H:\test_data\output\mpd\program\plot\20.png</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/e878680127ea61d323133c86e4ae4427?showdoc=.jpg" alt="" /></p>
<pre><code>[8, 20]</code></pre>
<pre><code class="language-python">#一下构造一个集合预报预报用于测试
import copy
sta_ensemble = meb.in_member_list(sta_all,["ob","ecmwf"]) #提取观测和ec预报的温度数据
sample_count = len(sta_ensemble.index)
meb.set_stadata_names(sta_ensemble,["ob",0])
for i in range(1,21):
sta_ensemble[i] = sta_ensemble[0].values + 3*np.random.randn(sample_count)
print(sta_ensemble) #打印构建的概率预报结果</code></pre>
<pre><code> level time dtime id lon lat ob 0 \
0 0 2019-01-06 20:00:00 0 54398 116.6 40.1 -6.3 -6.016
1 0 2019-01-06 20:00:00 0 54410 116.1 40.6 -7.0 -8.712
2 0 2019-01-06 20:00:00 0 54416 116.9 40.4 -9.5 -7.156
3 0 2019-01-06 20:00:00 0 54419 116.6 40.4 -7.6 -8.164
4 0 2019-01-06 20:00:00 0 54499 116.2 40.2 -5.4 -7.756
... ... ... ... ... ... ... ... ...
8336 0 2019-12-24 20:00:00 72 54410 116.1 40.6 -5.7 -6.452
8337 0 2019-12-24 20:00:00 72 54416 116.9 40.4 -6.1 -2.892
8338 0 2019-12-24 20:00:00 72 54419 116.6 40.4 -3.9 -4.068
8339 0 2019-12-24 20:00:00 72 54499 116.2 40.2 -1.7 -4.276
8340 0 2019-12-24 20:00:00 72 54412 116.6 40.7 -4.9 -6.140
1 2 ... 11 12 13 14 \
0 -7.566625 -6.437295 ... -9.558936 -4.635423 -10.903679 -4.811984
1 -12.347778 -12.479825 ... -1.083987 -4.070836 -8.878130 -9.745126
2 -10.008093 -5.485106 ... -9.679604 -5.349538 -7.903990 -16.055495
3 -4.042009 -7.960807 ... -7.751965 -4.617165 -4.992131 -5.873452
4 -6.230813 -8.777621 ... -8.992249 -6.110431 -11.436333 -9.774550
... ... ... ... ... ... ... ...
8336 -6.174789 -6.248870 ... -3.543153 -3.378253 -9.016864 -7.741782
8337 -3.518854 -2.581045 ... -2.962865 0.217085 0.092268 0.813913
8338 -4.586136 -2.664622 ... -5.660530 -3.670290 -6.224347 -3.528454
8339 -2.879676 -2.346811 ... -10.231584 -6.849797 -6.258126 -2.753724
8340 -5.290389 -4.875566 ... -1.940827 -4.807420 2.678128 -7.881496
15 16 17 18 19 20
0 -4.469150 -8.294785 -4.724681 -9.563168 -4.338762 -4.480826
1 -9.947279 -13.137376 -11.496214 -7.163931 -15.700751 -8.397545
2 -8.237473 -3.193422 -9.569797 -6.469109 -6.176549 -9.091136
3 -11.232146 -4.764005 -10.934699 -3.689846 -9.509125 -5.312376
4 -8.894650 -7.136992 -6.280799 -10.644803 -15.073189 -6.677337
... ... ... ... ... ... ...
8336 -5.520388 -7.147116 -9.735487 -1.328938 -4.930237 -8.043137
8337 -5.029152 -3.088244 -1.774121 -7.478492 -5.065802 -2.604842
8338 -8.077163 -2.243704 -2.921360 -5.901501 -10.235178 -6.571381
8339 -5.432124 -2.817576 -0.843425 -2.783871 -2.412422 0.113051
8340 -3.312898 -10.513175 -7.778807 -6.579240 -3.464236 -7.615210
[8341 rows x 28 columns]</code></pre>
<pre><code class="language-python">mpd.plot(sta_ensemble,mem.box_plot_ensemble)</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/b83739502e71eac7a0ebd71e477fee8f?showdoc=.jpg" alt="" /></p>
<pre><code class="language-python">mpd.plot(sta_ensemble,mem.rank_histogram)</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/c28fe06daa67d4a91aa52052b9bc6ccc?showdoc=.jpg" alt="" /></p>
<pre><code class="language-python">#一下构造一个概率预报用于测试
import copy
sta_ec = meb.in_member_list(sta_all,["ob","ecmwf"]) #提取观测和ec预报的温度数据
sta_ec = meb.between_value_range(sta_ec,-4,4) #测试数据中仅保留观测和预报的取值范围都在 -4 到4 之内的样本
sta_p0 =copy.deepcopy(sta_ec) #开始构建概率预报，设预报的目标是观测值大于 0 的概率
p_ob = np.zeros(len(sta_p0.index)) #原始观测值小于0，则观测样本的概率值被记为0
p_ob[sta_p0["ob"]> 0] = 1 #原始观测值大于0，则观测样本的概率值记为1
sta_p0["ob"] = p_ob[:]
p_fo_ecmwf = (sta_p0["ecmwf"].values+4)/8 #概率预报 采用 原始预报的线性函数代表，取值范围会在0-1之间。
p_fo_grapes = np.sqrt(p_fo_ecmwf) #生成grapes的概率预报，它的取值范围也会在0-1之间。
sta_p0["ecmwf"] = p_fo_ecmwf
sta_p0["grapes"] = p_fo_grapes
print(sta_p0) #打印构建的概率预报结果</code></pre>
<pre><code> level time dtime id lon lat ob ecmwf \
13 0 2019-03-02 20:00:00 0 54410 116.1 40.6 1.0 0.7530
18 0 2019-03-08 08:00:00 0 54398 116.6 40.1 1.0 0.4915
20 0 2019-03-08 08:00:00 0 54416 116.9 40.4 0.0 0.3940
21 0 2019-03-08 08:00:00 0 54419 116.6 40.4 0.0 0.3460
22 0 2019-03-08 08:00:00 0 54499 116.2 40.2 1.0 0.2590
... ... ... ... ... ... ... ... ...
8325 0 2019-12-13 20:00:00 72 54416 116.9 40.4 0.0 0.2430
8326 0 2019-12-13 20:00:00 72 54419 116.6 40.4 0.0 0.1170
8327 0 2019-12-13 20:00:00 72 54499 116.2 40.2 0.0 0.0915
8328 0 2019-12-13 20:00:00 72 54412 116.6 40.7 0.0 0.0055
8335 0 2019-12-24 20:00:00 72 54398 116.6 40.1 0.0 0.2285
grapes
13 0.867756
18 0.701071
20 0.627694
21 0.588218
22 0.508920
... ...
8325 0.492950
8326 0.342053
8327 0.302490
8328 0.074162
8335 0.478017
[860 rows x 9 columns]</code></pre>
<pre><code class="language-python">mpd.plot(sta_p0,mem.discrimination) #绘制概率预报区分能力图</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/bf4c1dc447dc57a9fd52ede3c4d15f32?showdoc=.jpg" alt="" /></p>
<pre><code class="language-python">mpd.plot(sta_p0,mem.reliability) #绘制概率预报区分能力图</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/0339a2291896927651a5b95744952ac9?showdoc=.jpg" alt="" /></p>
<pre><code class="language-python">mpd.plot(sta_p0,mem.probability.plot.roc)</code></pre>
<p><img src="https://www.showdoc.cc/server/api/common/visitfile/sign/80fea95d03e88a259486ab07d39240b6?showdoc=.jpg" alt="" /></p>
<pre><code class="language-python">mpd.plot(sta_p0,mem.comprehensive_probability)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/e11cc553b16a52d36e184ba986b0fb5a?showdoc=.jpg" alt="" /></p>
<pre><code class="language-python">ob_uv = meb.speed_angle_to_wind(ob_wind10m) #观测数据中原始数据为风速风向，则可以通过该函数转换为u，v
ec_uv = pd.read_hdf(r"H:\test_data\input\mpd\time_compair\ec_wind10m.h5","df") #读露点预报
grapes_uv = pd.read_hdf(r"H:\test_data\input\mpd\time_compair\grapes_wind10m.h5","df") #读露点预报
wind_all =meb.combine_on_obTime_id(ob_uv,[ec_uv,grapes_uv],need_match_ob=True) #合并预报观测数据，need_match_ob为缺省参数False
print(wind_all) #绘制风场检验的数据以u、v、...、u、v的形式排列，前两列默认为是观测，之后是预报 </code></pre>
<pre><code> level time dtime id lon lat u v \
0 0.0 2020-04-06 08:00:00 3 54511 119 40 1.394383 0.125288
1 0.0 2020-04-06 20:00:00 3 54511 119 40 -0.154433 -0.475553
2 0.0 2020-04-07 08:00:00 3 54511 119 40 -0.172479 3.195349
3 0.0 2020-04-07 20:00:00 3 54511 119 40 -1.628287 1.624402
4 0.0 2020-04-08 08:00:00 3 54511 119 40 -1.191180 0.145221
.. ... ... ... ... ... ... ... ...
337 0.0 2020-04-08 08:00:00 120 54511 119 40 1.165081 -0.287380
338 0.0 2020-04-06 08:00:00 126 54511 119 40 3.698046 -0.120228
339 0.0 2020-04-06 20:00:00 126 54511 119 40 0.818901 -1.824664
340 0.0 2020-04-07 08:00:00 126 54511 119 40 1.738523 1.177938
341 0.0 2020-04-07 20:00:00 126 54511 119 40 -0.738476 -0.945861
u_ecmwf v_ecmwf u_grapes v_grapes
0 0.239926 -0.012376 -0.381486 -0.779833
1 -1.186846 -0.657165 0.849711 1.086778
2 -1.340084 0.247622 -1.338557 -0.651618
3 -1.006525 1.442748 -1.816748 1.523139
4 -2.370218 0.372719 -2.036030 0.673341
.. ... ... ... ...
337 -0.346077 -3.494583 -0.338349 -4.285983
338 -0.958621 -1.002192 -0.142277 -0.728994
339 1.204885 -2.738038 0.646286 -2.432010
340 -1.205878 -2.645335 -1.821869 -1.263226
341 -0.119772 -4.037322 0.006010 -3.782761
[342 rows x 12 columns]</code></pre>
<pre><code class="language-python">mpd.plot(wind_all,mem.scatter_uv,add_randn_to_ob = 0.3)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/81d9102a93b6db4b74cc093750f7bd23" alt="" /></p>
<pre><code class="language-python">mpd.plot(wind_all,mem.scatter_uv_error)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/95ca05d00340eaa88971c90994a232af" alt="" /></p>
<pre><code class="language-python">mpd.plot(wind_all,mem.statisitic_uv)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/f8a2334a5c25e4b061dcbfb0f2c59186" alt="" /></p>
<pre><code class="language-python"></code></pre>