模糊逻辑检验
<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>模糊逻辑检验评估方法概述</h2>
<p>在以二分类的思路考察预报时，当实况达到某种阈值，就认为实况事件发生，实况记为1，否则记为0；同样的当预报达到某个阈值时，就认为预报事件发生，预报记为1，否则记为0。如果某站次实况事件发生的同时预报事件也发生，就记命中数+1；如果某站次实况事件发生，预报事件不发生，就记漏报数+1；如果某站次实况事件不发生，预报事件发生，就记空报数+1。 在这种检验视角下，事件的发生与否是明确的，非0即1。
然而，受可预报性的限制，天气预报总是不可能那么准确，当某个站点（或格点）上预报值为1时，在使用该预报时我们真的就会认为该点未来有事件发生吗？当然不是。更进一步的，通常来说，如果一个点预报值为1，且周围有一大片区域的预报值都为1，那我们会倾向于认为该点事件发生的可能性较高，相反，如果一个点预报值为1，但周围其它点都是0，那我们会认为该点事件发生的可能性其实不高。对此，可以通过下面的示例说明。</p>
<pre><code class="language-python">grid1 = meb.grid([100, 120, 0.05], [24, 40, 0.05]) # 设置读入后的数据的额网格
path_ob = r'H:\test_data\input\mem\mode\ob\rain03\20072611.000.nc'
path_fo_03 = r'H:\test_data\input\mem\mode\ec\rain03\20072608.003.nc'
grd_ob = meb.read_griddata_from_nc(path_ob, grid=grid1, time=&quot;2020072611&quot;, dtime=0, data_name=&quot;OBS&quot;)
grd_fo_03 = meb.read_griddata_from_nc(path_fo_03, grid=grid1, time=&quot;2020072608&quot;, dtime=3, data_name=&quot;ECMWF&quot;)
meb.plot_tools.plot_2d_grid_list([grd_ob,grd_fo_03],cmap=&quot;rain_3h&quot;,ncol = 2) #绘制原始数据的内容</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=6071fd4a44971df42532a1d16cc2346e&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">#以1毫米为阈值，记达到1m的为事件发生，格点值记为1，否则记为事件不发生，格点值记为0
h_ob = meb.grid_data(grid1)
h_ob.values[grd_ob.values&gt;=1] = 1
h_fo = meb.grid_data(grid1)
h_fo.values[grd_fo_03.values&gt;=1] = 1
meb.plot_tools.plot_2d_grid_list([h_ob,h_fo],ncol = 2) #绘制01形式事件发生与否的数据</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=661d2d5d331dd0741269a07ffb11c910&amp;file=file.png" alt="" /></p>
<p>在上面的示例的预报场中江西南部个别格点上取值为1，但这些取值为1的区域非常零散，我们会倾向于认为、这些区域的事件发生可能性不会真的达到1。相反，我们会认为在湖北中部，事件发生的可能性会是1或接近1。
那如何定量的表达一个01场对应的事件发生可能性的场呢？有一种方式，就是采用平滑的算法获取一个点周围一定范围内的平均值。平滑方法实际上可以有很多种，在模糊逻辑检验中，使用的是一种基于二维傅里叶变换的窗口平滑算法，通过下面的示例我们可以看到这个算法的调用方法和平滑效果。</p>
<pre><code class="language-python">#采用模糊逻辑中的hoods2dsmooth_grd函数对01形式的场进行平滑，平滑窗口的大小由half_window_size控制，
#half_window_size=1时代表平滑窗口是以某格点为中心的3×3的窗口，该参数越大返回场越平滑
sm_ob = meteva.method.fuzzy_logic.lib.hoods2dsmooth.hoods2dsmooth_grd(h_ob,half_window_size = 8)
sm_fo = meteva.method.fuzzy_logic.lib.hoods2dsmooth.hoods2dsmooth_grd(h_fo,half_window_size =8)
meb.plot_tools.plot_2d_grid_list([sm_ob,sm_fo],ncol = 2)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=99f149aecf500d4e033282f6f5783ee7&amp;file=file.png" alt="" /></p>
<p>从上面的结果可以看到通过平滑预报场在江西南部的取值由1降到了0.2以下。在模糊逻辑算法的提出者看来，平滑后的场能反映对事件概率的预报能力。在应用中会有两个问题：</p>
<ol>
<li>平滑窗口该怎么选取？ </li>
<li>要不要对01形式的观测也做平滑？<br />
关于问题1，不同的窗口宽度代表我们检验的是不同尺度的信息，并没有一个最优的平滑窗口参数。关于问题2， 模糊逻辑检验是一个大类，其中又包含多种不同的检验算法，在有的算法中，需要对观测进行平滑，有的则不需要。 </li>
</ol>
<p>模糊逻辑的检验方法包括 joint、fuzzy、mincvr、multi_event和pagmatic几种，它们的具体介绍如下：</p>
<h2>joint方法</h2>
<p>本方法为模糊逻辑（fuzzy logic）中的一种特例，用于解决观测数据与预测结果无法在时空层面上完全一致导致的校验困难的问题，通过假设某一时间某地的预报数据可以从相邻的时空预测场中导出。通过假定领域内事件发生的频率即为该点处的概率，同时，使用联合概率来修复所有概率加起来不为1的问题，最后，计算POD，FAR，ETS，验证模型的准确性。详细可以参考 </p>
<ul>
<li>Ebert, E. E. (2008) Fuzzy verification of high resolution gridded forecasts: A review and proposed framework. Meteorol. Appl., 15, 51–64. doi:10.1002/met.25</li>
<li>Ebert EE. 2002. Fuzzy verification: giving partial credit to erroneous forecasts. In NCAR/FAA Verification Workshop: Making Verification More Meaningful, NCAR, Boulder, 30 July – 1 August 2002,</li>
<li>Damrath U. 2004. Verification against precipitation observations of a high density network – what did we learn? In International Verification Methods Workshop, Montreal, 15–17 September 2004,</li>
</ul>
<p>具体来说，在joint方法中，是对01形式的观测和预报都进行平滑操作，得到的平滑场代表事件的概率。假设一个格点上观测和预报的概率分别是pob和pfo，则该点上</p>
<ul>
<li>命中数=pob * pfo </li>
<li>漏报数=pob * (1-pfo)</li>
<li>空报数=（1-pob） * pfo</li>
<li>正确否定数=（1-pob） * (1-pfo)
全场的命中数、漏报数、空报数和正确否定数是每个格点上的相应数值的累加。基于全场的全场的命中数、漏报数、空报数和正确否定数，采用和传统二分类检验一样的计算公式计算命中率、空报率和ETS评分等。</li>
</ul>
<p>从上面的方法介绍可以看到，joint和传统的二分类检验方法的异同：</p>
<ol>
<li>在joint方法中命中数、漏报数、空报数和正确否定数可以是0-1之间的小数，而二分类检验中则只能是0或1</li>
<li>在一个点上命中数+漏报数+空报数+正确否定数=1</li>
</ol>
<p><font face="黑体" color=Blue size=3><strong>joint(grd_ob, grd_fo, half_window_sizes, thresholds = None,compare=">=", show=False)</strong></font></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>grd_ob</font></strong></td>
<td style="text-align: left;">网格数据形式的观测数据，只支持包含单一平面场的网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>grd_fo</font></strong></td>
<td style="text-align: left;">网格数据形式的预报数据，只支持包含单一平面场的网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong>half_window_sizes</strong></td>
<td style="text-align: left;">包含半边窗口的网格数的列表，例如当half_window_size = [1],时，计算区域中部的一个网格点的平均用到 3×3个格点，当half_window_size = [2]时 用到了5×5个格点求平均</td>
</tr>
<tr>
<td style="text-align: left;"><strong>thresholds</strong></td>
<td style="text-align: left;">缺省值时，程序自动提取观测和预报数据的0.05, 0.1, 0.25, 0.333, 0.5, 0.666, 0.75, 0.9, 0.95的分位值作为阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong>compare</strong></td>
<td style="text-align: left;">字符串，代表threshold的规则，可选值为 >=, >, <, <=</td>
</tr>
<tr>
<td style="text-align: left;"><strong>show</strong></td>
<td style="text-align: left;">布尔值，是否打印日志</td>
</tr>
</tbody>
</table>
<p><font face="黑体" color=green size=4><strong>返回结果内容说明</strong></font></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=green size=3>pod</font></strong></td>
<td style="text-align: left;">命中率， 二维数组m × n, m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>far</font></strong></td>
<td style="text-align: left;">空报率，二维数组m × n,m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>ets</font></strong></td>
<td style="text-align: left;">ETS评分，二维数组m × n,m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>thresholds_ob</font></strong></td>
<td style="text-align: left;">观测数据用到的阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>thresholds_fo</font></strong></td>
<td style="text-align: left;">预报数据用到的阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>window_sizes</font></strong></td>
<td style="text-align: left;">窗口的网格数</td>
</tr>
</tbody>
</table>
<hr />
<p>$$POD=\frac{hits}{hits + misses}$$</p>
<p>$$FAR=\frac{false\ alarms}{false\ alarms}$$</p>
<p>$$ETS=\frac{hits - hits<em>{random}}{hits + misses + false\ alarms - hits</em>{random}}$$</p>
<p>$$其中， hits_{random} = \frac{1}{N}(observed\ yes × forecast\ yes)$$</p>
<p>POD代表预测结果的命中率， FAR代表误报率， ETS代表公正预兆得分</p>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">grid1 = meb.grid([100, 120, 0.05], [24, 40, 0.05])
path_ob = r'H:\test_data\input\mem\mode\ob\rain03\20072611.000.nc'
path_fo_03 = r'H:\test_data\input\mem\mode\ec\rain03\20072608.003.nc'
grd_ob = meb.read_griddata_from_nc(path_ob, grid=grid1, time=&quot;2020072611&quot;, dtime=0, data_name=&quot;OBS&quot;)
grd_fo_03 = meb.read_griddata_from_nc(path_fo_03, grid=grid1, time=&quot;2020072608&quot;, dtime=3, data_name=&quot;ECMWF&quot;)
look = mem.fuzzy_logic.joint(grd_ob,grd_fo_03, half_window_sizes=[1,2,4,8])
print(look)</code></pre>
<pre><code>{'pod': array([[0.99750391, 0.99750391, 0.74922438, 0.66650474, 0.50013473,
0.7666196 , 0.73307372, 0.48760858, 0.32745843],
[0.99550244, 0.99550244, 0.74840328, 0.66591527, 0.50016503,
0.7652279 , 0.73220114, 0.4864535 , 0.32696631],
[0.99165483, 0.99165483, 0.74673918, 0.66459193, 0.5001135 ,
0.7614434 , 0.72862712, 0.48173299, 0.32407346],
[0.9840334 , 0.9840334 , 0.74292745, 0.66163482, 0.49992005,
0.74857448, 0.71400291, 0.46344696, 0.30841125]]), 'far': array([[0.00249609, 0.00249609, 0.00197296, 0.00184672, 0.00141164,
0.23335659, 0.26730626, 0.51268182, 0.67255853],
[0.00449756, 0.00449756, 0.00363236, 0.00345602, 0.00263389,
0.23483276, 0.26849615, 0.51403775, 0.67308448],
[0.00834517, 0.00834517, 0.00684681, 0.0066214 , 0.00502855,
0.23892403, 0.27268214, 0.51912355, 0.67607198],
[0.0159666 , 0.0159666 , 0.01375129, 0.01310202, 0.00976078,
0.25261148, 0.28847964, 0.5381244 , 0.69203124]]), 'ets': array([[0.1989431 , 0.1989431 , 0.00521905, 0.00375452, 0.00232012,
0.48129227, 0.47489763, 0.27431492, 0.17099312],
[0.19809664, 0.19809664, 0.00903571, 0.00644123, 0.0040718 ,
0.47915289, 0.47351741, 0.27319388, 0.17062719],
[0.19647176, 0.19647176, 0.01597475, 0.0112984 , 0.00734657,
0.47315453, 0.46814212, 0.26881332, 0.16851417],
[0.19326262, 0.19326262, 0.02703662, 0.01990166, 0.01366213,
0.45273331, 0.44701443, 0.25239795, 0.1572739 ]]), 'thresholds_ob': array([0. , 0. , 0. , 0. , 0. , 0.03456, 0.2748 ,
3.07148, 7.69968]), 'thresholds_fo': array([0.00000e+00, 0.00000e+00, 1.68000e-03, 1.60800e-02, 1.06960e-01,
4.22800e-01, 9.61760e-01, 3.61744e+00, 6.05900e+00]), 'window_sizes': [3, 5, 9, 17]}</code></pre>
<h2>fuzzy方法</h2>
<p>本方法为模糊逻辑（fuzzy logic）中的一种特例，用于解决观测数据与预测结果无法在时空层面上完全一致导致的校验困难的问题，通过假设某一时间某地的预报数据可以从相邻的时空预测场中导出。通过假定领域内事件发生的频率即为该点处的概率，同时，使用联合概率来修复所有概率加起来不为1的问题，最后，计算POD，FAR，ETS，验证模型的准确性。详细可以参考</p>
<ul>
<li>Ebert, E. E. (2008) Fuzzy verification of high resolution gridded forecasts: A review and proposed framework. Meteorol. Appl., 15, 51–64. doi:10.1002/met.25</li>
<li>Ebert EE. 2002. Fuzzy verification: giving partial credit to erroneous forecasts. In NCAR/FAA Verification Workshop: Making Verification More Meaningful, NCAR, Boulder, 30 July – 1 August 2002,</li>
<li>Damrath U. 2004. Verification against precipitation observations of a high density network – what did we learn? In International Verification Methods Workshop, Montreal, 15–17 September 2004,</li>
</ul>
<p>具体来说，在joint方法中，是对01形式的观测和预报都进行平滑操作，得到的平滑场代表事件的概率。假设一个格点上观测和预报的概率分别是pob和pfo，则该点上</p>
<ul>
<li>命中数=min(pob , pfo) # pob和pfo既要大，又要接近，命中才大</li>
<li>漏报数=min(pob , (1-pfo)) #pob越小，或pfo越大，漏报越小</li>
<li>空报数=min(（1-pob） , pfo) #pob越大，或pfo越小，空报越小</li>
<li>正确否定数=min(（1-pob） , (1-pfo))
全场的命中数、漏报数、空报数和正确否定数是每个格点上的相应数值的累加。基于全场的全场的命中数、漏报数、空报数和正确否定数，采用和传统二分类检验一样的计算公式计算命中率、空报率和ETS评分等。</li>
</ul>
<p>从上面的方法介绍可以看到，fuzzy和传统的二分类检验方法很不相同：</p>
<ol>
<li>在fuzzy方法中如果某个格点上观测和预报中有一个取值为0，那么命中数、漏报数、空报数和正确否定都为0，也就相当于该点没进入检验样本中。</li>
<li>在一个点上命中数+漏报数+空报数+正确否定数 !=1 </li>
</ol>
<p><font face="黑体" color=Blue size=3><strong>fuzzy(grd_ob, grd_fo, half_window_sizes, thresholds = None,compare=">=", show=False)</strong></font></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>grd_ob</font></strong></td>
<td style="text-align: left;">网格数据形式的观测数据，只支持包含单一平面场的网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>grd_fo</font></strong></td>
<td style="text-align: left;">网格数据形式的预报数据，只支持包含单一平面场的网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong>half_window_sizes</strong></td>
<td style="text-align: left;">包含半边窗口的网格数的列表，例如当half_window_size = [1],时，计算区域中部的一个网格点的平均用到 3×3个格点，当half_window_size = [2]时 用到了5×5个格点求平均</td>
</tr>
<tr>
<td style="text-align: left;"><strong>thresholds</strong></td>
<td style="text-align: left;">缺省值时，程序自动提取观测和预报数据的0.05, 0.1, 0.25, 0.333, 0.5, 0.666, 0.75, 0.9, 0.95的分位值作为阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong>compare</strong></td>
<td style="text-align: left;">字符串，代表threshold的规则，可选值为 >=, >, <, <=</td>
</tr>
<tr>
<td style="text-align: left;"><strong>show</strong></td>
<td style="text-align: left;">布尔值，是否打印日志</td>
</tr>
</tbody>
</table>
<p><font face="黑体" color=green size=4><strong>返回结果内容说明</strong></font></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=green size=3>pod</font></strong></td>
<td style="text-align: left;">命中率， 二维数组m × n, m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>far</font></strong></td>
<td style="text-align: left;">空报率，二维数组m × n,m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>ets</font></strong></td>
<td style="text-align: left;">ETS评分，二维数组m × n,m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>thresholds_ob</font></strong></td>
<td style="text-align: left;">观测数据用到的阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>thresholds_fo</font></strong></td>
<td style="text-align: left;">预报数据用到的阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>window_sizes</font></strong></td>
<td style="text-align: left;">窗口的网格数</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">look = mem.fuzzy_logic.fuzzy(grd_ob,grd_fo_03, half_window_sizes=[1,2,4,8])
print(look)</code></pre>
<pre><code>{'pod': array([[0.99626754, 0.99626754, 0.74892688, 0.666411 , 0.50027132,
0.76508055, 0.73177471, 0.48773751, 0.32906524],
[0.99329146, 0.99329146, 0.74786822, 0.66569862, 0.50038053,
0.76126997, 0.72861906, 0.48699495, 0.33170941],
[0.98761402, 0.98761402, 0.74558429, 0.66402402, 0.50036863,
0.7510783 , 0.7194299 , 0.48396571, 0.33689589],
[0.97653067, 0.97653067, 0.74020368, 0.65993317, 0.49989128,
0.72776125, 0.6978538 , 0.47633471, 0.34496972]]), 'far': array([[0.00373246, 0.00373246, 0.00315356, 0.00297293, 0.00227144,
0.23489584, 0.26860232, 0.5125501 , 0.67095164],
[0.00670854, 0.00670854, 0.00577589, 0.00564654, 0.00438568,
0.23878941, 0.27206381, 0.51348337, 0.66834077],
[0.01238598, 0.01238598, 0.01111825, 0.01112483, 0.00895435,
0.24926889, 0.2818086 , 0.51682913, 0.66324514],
[0.02346933, 0.02346933, 0.02292477, 0.02315101, 0.01926151,
0.27327886, 0.3043383 , 0.52499575, 0.65544413]]), 'ets': array([[0.33022831, 0.33022831, 0.00858054, 0.00539634, 0.00258548,
0.47835511, 0.47262874, 0.27402369, 0.17192953],
[0.32776382, 0.32776382, 0.01497374, 0.0093451 , 0.00447591,
0.47164127, 0.46732319, 0.27253401, 0.17344041],
[0.32313068, 0.32313068, 0.02602994, 0.01624011, 0.00775728,
0.45368044, 0.45226374, 0.26747699, 0.17615758],
[0.31417408, 0.31417408, 0.04336553, 0.02734287, 0.01322068,
0.41368231, 0.41805147, 0.25525126, 0.1792177 ]]), 'thresholds_ob': array([0. , 0. , 0. , 0. , 0. , 0.03456, 0.2748 ,
3.07148, 7.69968]), 'thresholds_fo': array([0.00000e+00, 0.00000e+00, 1.68000e-03, 1.60800e-02, 1.06960e-01,
4.22800e-01, 9.61760e-01, 3.61744e+00, 6.05900e+00]), 'window_sizes': [3, 5, 9, 17]}</code></pre>
<h2>pragmatic</h2>
<p>本方法为模糊逻辑（fuzzy logic）中的一种特例，用于解决观测数据与预测结果无法在时空层面上完全一致导致的校验困难的问题，通过假设某一时间某地的预报数据可以从相邻的时空预测场中导出，经过排序算法及核回归方法（kernel regression）获得处理后的概率预报（probabilistic postprocessed forecast），使用Brier score与Brier skill score比较处理后的预测结果与观测场数据，得出模拟结果是否可信的结果。详细可参考Theis, S. E., Hense, A. Damrath, U. (2005) Probabilistic precipitation forecasts from a deterministic model: A pragmatic approach. Meteorol. Appl., 12, 257–268.</p>
<p>具体来说，在pragmatic方法中，是只对01形式的预报都进行平滑操作，得到的平滑场代表事件的概率，观测仍然采用原始的01场，然后采用点对点概率预报检验中的bs和bss检验算法进行检验。</p>
<p><font face="黑体" color=Blue size=3><strong>pragmatic(grd_ob, grd_fo, half_window_sizes, thresholds = None,compare=">=", show=False)</strong></font></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>grd_ob</font></strong></td>
<td style="text-align: left;">网格数据形式的观测数据，只支持包含单一平面场的网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>grd_fo</font></strong></td>
<td style="text-align: left;">网格数据形式的预报数据，只支持包含单一平面场的网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong>half_window_sizes</strong></td>
<td style="text-align: left;">包含半边窗口的网格数的列表，例如当half_window_size = [1],时，计算区域中部的一个网格点的平均用到 3×3个格点，当half_window_size = [2]时 用到了5×5个格点求平均</td>
</tr>
<tr>
<td style="text-align: left;"><strong>thresholds</strong></td>
<td style="text-align: left;">缺省值时，程序自动提取观测和预报数据的0.05, 0.1, 0.25, 0.333, 0.5, 0.666, 0.75, 0.9, 0.95的分位值作为阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong>compare</strong></td>
<td style="text-align: left;">字符串，代表threshold的规则，可选值为 >=, >, <, <=</td>
</tr>
<tr>
<td style="text-align: left;"><strong>show</strong></td>
<td style="text-align: left;">布尔值，是否打印日志</td>
</tr>
</tbody>
</table>
<p><font face="黑体" color=green size=4><strong>返回结果内容说明</strong></font></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=green size=3>bs</font></strong></td>
<td style="text-align: left;">命中率， 二维数组m × n, m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>bss</font></strong></td>
<td style="text-align: left;">空报率，二维数组m × n,m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>thresholds_ob</font></strong></td>
<td style="text-align: left;">观测数据用到的阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>thresholds_fo</font></strong></td>
<td style="text-align: left;">预报数据用到的阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>window_sizes</font></strong></td>
<td style="text-align: left;">窗口的网格数</td>
</tr>
</tbody>
</table>
<hr />
<p>$$BS=\frac{1}{n}\sum_{k=1}^{n}{(y_k-o_k)^2}$$</p>
<p>$$BS代表预测结果的平方误差。其中y_k代表预测的事件发生几率，而o_k代表对应的观测事件是否发生，1代表已发生，0代表未发生$$</p>
<p>$$BSS=1-\frac{BS}{BS_{reference}}$$</p>
<p>$$BSS代表预测结果的可改进度。其中BS_{reference}为假定的最佳预测的BS值$$</p>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">look = mem.fuzzy_logic.pragmatic(grd_ob,grd_fo_03, half_window_sizes=[1,2,4,8])
print(look)</code></pre>
<pre><code>{'bs': array([[0.00124914, 0.00124914, 0.23834974, 0.31847384, 0.48419575,
0.14396243, 0.12446423, 0.09579048, 0.06321414],
[0.00225231, 0.00225231, 0.23192674, 0.30966383, 0.47418034,
0.13688835, 0.11872908, 0.09137602, 0.06050345],
[0.00418468, 0.00418468, 0.22275198, 0.29656584, 0.4591749 ,
0.12591444, 0.10976235, 0.08419843, 0.05593387],
[0.00802702, 0.00802702, 0.21162679, 0.28165156, 0.44029288,
0.11243604, 0.09831977, 0.07460869, 0.04929544]]), 'bss': array([[ -inf, -inf, -inf, -inf, -inf,
0.35284637, 0.33620107, -0.06427257, -0.33063798],
[ -inf, -inf, -inf, -inf, -inf,
0.38464646, 0.36678806, -0.01522603, -0.27357865],
[ -inf, -inf, -inf, -inf, -inf,
0.43397741, 0.41460989, 0.06452006, -0.17739053],
[ -inf, -inf, -inf, -inf, -inf,
0.49456678, 0.47563605, 0.17106607, -0.03765367]]), 'thresholds_ob': array([0. , 0. , 0. , 0. , 0. , 0.03456, 0.2748 ,
3.07148, 7.69968]), 'thresholds_fo': array([0.00000e+00, 0.00000e+00, 1.68000e-03, 1.60800e-02, 1.06960e-01,
4.22800e-01, 9.61760e-01, 3.61744e+00, 6.05900e+00]), 'window_sizes': [3, 5, 9, 17]}</code></pre>
<h2>minimum coverage</h2>
<p>本方法为模糊逻辑（fuzzy logic）中的一种特例，用于解决观测数据与预测结果无法在时空层面上完全一致导致的校验困难的问题，通过假设在关心的区域中，预测事件与实际事件都不准确，两者均有可能在区域中的任意一点发生，使用以下模型，当预测的概率超过了预测和实际区域两者的最小覆盖区域的发生概率则认为事件发生，否则则认为不发生，最后，计算POD，FAR，ETS，验证模型的准确性。详细可以参考</p>
<ul>
<li>Ebert, E. E. (2008) Fuzzy verification of high resolution gridded forecasts: A review and proposed framework. Meteorol. Appl., 15, 51–64. doi:10.1002/met.25</li>
<li>Ebert, E. E. (2009) Neighborhood verification: A strategy for rewarding close forecasts. Wea. Forecasting, 24, 1498–1510, doi:10.1175/2009WAF2222251.1.</li>
<li>Damrath U. 2004. Verification against precipitation observations of a high density network – what did we learn? In International Verification Methods Workshop, Montreal, 15–17 September 2004,</li>
</ul>
<p>具体来说，在 minimum coverage方法中，是对01形式的观测和预报都进行平滑操作，得到的平滑场代表事件的概率。进一步的，设定一个新的阈值(=1/(window_sizes*window_sizes))，将概率达到阈值的部分设置为1，否则设置为0，由此重新将观测预报转换成01场。之后基于该01场，采用二分类的检验算法计算命中率，空报率，ETS评分等。</p>
<p><font face="黑体" color=Blue size=3><strong>minimum_coverage(grd_ob, grd_fo, half_window_sizes, thresholds = None,compare=">=", show=False)</strong></font></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>grd_ob</font></strong></td>
<td style="text-align: left;">网格数据形式的观测数据，只支持包含单一平面场的网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>grd_fo</font></strong></td>
<td style="text-align: left;">网格数据形式的预报数据，只支持包含单一平面场的网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong>half_window_sizes</strong></td>
<td style="text-align: left;">包含半边窗口的网格数的列表，例如当half_window_size = [1],时，计算区域中部的一个网格点的平均用到 3×3个格点，当half_window_size = [2]时 用到了5×5个格点求平均</td>
</tr>
<tr>
<td style="text-align: left;"><strong>thresholds</strong></td>
<td style="text-align: left;">缺省值时，程序自动提取观测和预报数据的0.05, 0.1, 0.25, 0.333, 0.5, 0.666, 0.75, 0.9, 0.95的分位值作为阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong>compare</strong></td>
<td style="text-align: left;">字符串，代表threshold的规则，可选值为 >=, >, <, <=</td>
</tr>
<tr>
<td style="text-align: left;"><strong>show</strong></td>
<td style="text-align: left;">布尔值，是否打印日志</td>
</tr>
</tbody>
</table>
<p><font face="黑体" color=green size=4><strong>返回结果内容说明</strong></font></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=green size=3>pod</font></strong></td>
<td style="text-align: left;">命中率， 二维数组m × n, m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>far</font></strong></td>
<td style="text-align: left;">空报率，二维数组m × n,m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>ets</font></strong></td>
<td style="text-align: left;">ETS评分，二维数组m × n,m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>thresholds_ob</font></strong></td>
<td style="text-align: left;">观测数据用到的阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>thresholds_fo</font></strong></td>
<td style="text-align: left;">预报数据用到的阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>window_sizes</font></strong></td>
<td style="text-align: left;">窗口的网格数</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">look = mem.fuzzy_logic.minimum_coverage(grd_ob,grd_fo_03, half_window_sizes=[1,2,4,8])
print(look)</code></pre>
<pre><code>{'pod': array([[1. , 1. , 1. , 1. , 1. ,
0.7800676 , 0.75222267, 0.51466578, 0.36911366],
[1. , 1. , 1. , 1. , 1. ,
0.7802016 , 0.77063791, 0.57334242, 0.41159451],
[1. , 1. , 1. , 1. , 1. ,
0.79263808, 0.79744575, 0.64931277, 0.47954263],
[1. , 1. , 1. , 1. , 1. ,
0.82835239, 0.8404791 , 0.73663604, 0.60891739]]), 'far': array([[0. , 0. , 0.22541 , 0.30300417, 0.46852495,
0.23167798, 0.26058148, 0.47457052, 0.63962716],
[0. , 0. , 0.19647921, 0.26359335, 0.41502164,
0.22122295, 0.25628166, 0.44508543, 0.60163866],
[0. , 0. , 0.1593524 , 0.21660801, 0.34735591,
0.20161085, 0.24677072, 0.37648007, 0.53348468],
[0. , 0. , 0.11918024, 0.163229 , 0.26489073,
0.16857257, 0.21538435, 0.29262303, 0.41912331]]), 'ets': array([[ nan, nan, 0. , 0. , 0. ,
0.47895572, 0.4835546 , 0.29638933, 0.19293652],
[ nan, nan, 0. , 0. , 0. ,
0.46277451, 0.48216909, 0.32761041, 0.21767914],
[ nan, nan, 0. , 0. , 0. ,
0.45149099, 0.48315689, 0.38565621, 0.26068616],
[ nan, nan, 0. , 0. , 0. ,
0.45380024, 0.50045343, 0.45819029, 0.35024374]]), 'thresholds_ob': array([0. , 0. , 0. , 0. , 0. , 0.03456, 0.2748 ,
3.07148, 7.69968]), 'thresholds_fo': array([0.00000e+00, 0.00000e+00, 1.68000e-03, 1.60800e-02, 1.06960e-01,
4.22800e-01, 9.61760e-01, 3.61744e+00, 6.05900e+00]), 'window_sizes': [3, 5, 9, 17]}</code></pre>
<h2>multi_event</h2>
<p>本方法为模糊逻辑（fuzzy logic）中的一种特例
具体来说，在 minimum coverage方法中，是对01形式的预报都进行平滑操作，得到的平滑场代表事件的概率。进一步的，设定一个新的阈值(=1/(window_sizes*window_sizes))，将概率达到阈值的部分设置为1，否则设置为0，由此重新将预报转换成01场。之后基于原始的01观测场和重构后的预报01场，采用二分类的检验算法计算命中率，空报率，ETS评分等。</p>
<p><font face="黑体" color=Blue size=3><strong>multi_event(grd_ob, grd_fo, half_window_sizes, thresholds = None,compare=">=", show=False)</strong></font></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>grd_ob</font></strong></td>
<td style="text-align: left;">网格数据形式的观测数据，只支持包含单一平面场的网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>grd_fo</font></strong></td>
<td style="text-align: left;">网格数据形式的预报数据，只支持包含单一平面场的网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong>half_window_sizes</strong></td>
<td style="text-align: left;">包含半边窗口的网格数的列表，例如当half_window_size = [1],时，计算区域中部的一个网格点的平均用到 3×3个格点，当half_window_size = [2]时 用到了5×5个格点求平均</td>
</tr>
<tr>
<td style="text-align: left;"><strong>thresholds</strong></td>
<td style="text-align: left;">缺省值时，程序自动提取观测和预报数据的0.05, 0.1, 0.25, 0.333, 0.5, 0.666, 0.75, 0.9, 0.95的分位值作为阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong>compare</strong></td>
<td style="text-align: left;">字符串，代表threshold的规则，可选值为 >=, >, <, <=</td>
</tr>
<tr>
<td style="text-align: left;"><strong>show</strong></td>
<td style="text-align: left;">布尔值，是否打印日志</td>
</tr>
</tbody>
</table>
<p><font face="黑体" color=green size=4><strong>返回结果内容说明</strong></font></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=green size=3>pod</font></strong></td>
<td style="text-align: left;">命中率， 二维数组m × n, m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>far</font></strong></td>
<td style="text-align: left;">空报率，二维数组m × n,m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>ets</font></strong></td>
<td style="text-align: left;">ETS评分，二维数组m × n,m对应输入参数half_window_sizes的长度， n对应输入参数thresholds的长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>thresholds_ob</font></strong></td>
<td style="text-align: left;">观测数据用到的阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>thresholds_fo</font></strong></td>
<td style="text-align: left;">预报数据用到的阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>window_sizes</font></strong></td>
<td style="text-align: left;">窗口的网格数</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">look = mem.fuzzy_logic.multi_event(grd_ob,grd_fo_03, half_window_sizes=[1,2,4,8])
print(look)</code></pre>
<pre><code>{'pod': array([[1. , 1. , 1. , 1. , 1. ,
0.74057355, 0.7135801 , 0.46398404, 0.30757322],
[1. , 1. , 1. , 1. , 1. ,
0.69056978, 0.67701702, 0.4403951 , 0.28857694],
[1. , 1. , 1. , 1. , 1. ,
0.63109417, 0.62536166, 0.40062757, 0.25630056],
[1. , 1. , 1. , 1. , 1. ,
0.56259831, 0.54280055, 0.33904834, 0.21334074]]), 'f': array([[0. , 0. , 1. , 1. , 1. ,
0.10906765, 0.08081397, 0.05187094, 0.03350417],
[0. , 0. , 1. , 1. , 1. ,
0.09711249, 0.06788865, 0.04603753, 0.03135967],
[0. , 0. , 1. , 1. , 1. ,
0.08045447, 0.05344397, 0.0358927 , 0.02711777],
[0. , 0. , 1. , 1. , 1. ,
0.05666185, 0.03828494, 0.02199897, 0.02030904]]), 'hk': array([[1. , 1. , 0. , 0. , 0. ,
0.6315059 , 0.63276613, 0.41211309, 0.27406904],
[1. , 1. , 0. , 0. , 0. ,
0.5934573 , 0.60912837, 0.39435756, 0.25721727],
[1. , 1. , 0. , 0. , 0. ,
0.55063969, 0.57191769, 0.36473487, 0.22918279],
[1. , 1. , 0. , 0. , 0. ,
0.50593646, 0.50451561, 0.31704937, 0.1930317 ]]), 'thresholds_ob': array([0. , 0. , 0. , 0. , 0. , 0.03456, 0.2748 ,
3.07148, 7.69968]), 'thresholds_fo': array([0.00000e+00, 0.00000e+00, 1.68000e-03, 1.60800e-02, 1.06960e-01,
4.22800e-01, 9.61760e-01, 3.61744e+00, 6.05900e+00]), 'window_sizes': [3, 5, 9, 17]}</code></pre>