CRA
<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 numpy as np
import time
import math</code></pre>
<h1>CRA方法</h1>
<p>CRA方法是基于目标识别和刚体变换的方法，诊断要素目标（例如雨带）的位置偏差、形态偏差、强度偏差的一种空间检验方法。CRA的实现步骤包括如下部分：</p>
<ol>
<li>根据MODE算法模块中的目标识别、目标匹配和目标合并方法，获取具有匹对关系的观测和预报目标</li>
<li>根据刚体变换方法（rigider）对每个匹配上的预报目标进行平移和旋转，根据平移和旋转的场计算预报误差的不同分量。</li>
</ol>
<p><font face="黑体" color=Blue size=4><strong>craer(look_merge, stages = True,translate = True,rotate = True)</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=3>look_merge</font></strong></td>
<td style="text-align: left;">MODE算法完成目标识别匹配和合并后的输出结果</td>
</tr>
<tr>
<td style="text-align: left;"><strong>stages</strong></td>
<td style="text-align: left;">布尔值，该参数为True时，采用方案1同时计算最优的平移和旋转，否作采用方案2</td>
</tr>
<tr>
<td style="text-align: left;"><strong>translate</strong></td>
<td style="text-align: left;">stages = False时有效，当translate为True时代表会计算最优平移，否则不计算</td>
</tr>
<tr>
<td style="text-align: left;"><strong>rotate</strong></td>
<td style="text-align: left;">stages = False时有效，当rotate为True时代表会计算最优旋转，否则不计算</td>
</tr>
</tbody>
</table>
<p><strong>————————————————————————————————————————————————————————————————————————————————————————————————————————————</strong></p>
<p><font face="黑体" color=green size=4><strong>craer 返回结果是一个两层嵌套的字典，自动的第一级关键词是每个成功匹配的目标的编号，第二级关键词对应该目标的各类检验结果，内容说明如下：</strong></font> </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=green size=3>grd_ob</font></strong></td>
<td style="text-align: center;">从输入数据种继承的二维矩阵形式的网格观测场</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>grd_fo</font></strong></td>
<td style="text-align: center;">从输入数据种继承的二维矩阵形式的网格预报场</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>grd_fo_shift</font></strong></td>
<td style="text-align: center;">平移后的预报场，网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>grd_fo_shift_rotate</font></strong></td>
<td style="text-align: center;">平移+旋转后的预报场，网格数据</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>delta_lon</font></strong></td>
<td style="text-align: center;">质心经度偏差（预报-观测）</td>
<td style="text-align: left;">°</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>delta_lat</font></strong></td>
<td style="text-align: center;">质心纬度偏差（预报-观测）</td>
<td style="text-align: left;">°</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>delta_angle</font></strong></td>
<td style="text-align: center;">目标倾角偏差（预报-观测）</td>
<td style="text-align: left;">°</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>ob_mean</font></strong></td>
<td style="text-align: center;">观测场均值</td>
<td style="text-align: left;">和要素原始单位一致</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>fo_mean</font></strong></td>
<td style="text-align: center;">原始预报场均值</td>
<td style="text-align: left;">和要素原始单位一致</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>MSE_total</font></strong></td>
<td style="text-align: center;">总误差=原始预报和观测之间误差场的均方</td>
<td style="text-align: left;">要素原始单位的平方，例如当grd_ob是降水场时，单位为mm<sup>2</sup></td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>MSE_shift_left</font></strong></td>
<td style="text-align: center;">平移场残差=平移后的预报与原始预报之间误差场的均方</td>
<td style="text-align: left;">要素原始单位的平方</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>MSE_shift</font></strong></td>
<td style="text-align: center;">平移误差=MSE_total -MSE_shift_left</td>
<td style="text-align: left;">要素原始单位的平方</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>MSE_shift_rotate_left</font></strong></td>
<td style="text-align: center;">平移旋转场残差=平移旋转后的预报与观测之间误差场的均方</td>
<td style="text-align: left;">要素原始单位的平方</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>MSE_rotate</font></strong></td>
<td style="text-align: center;">旋转误差= MSE_shift_left - MSE_shift_rotate_left</td>
<td style="text-align: left;">要素原始单位的平方</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>MSE_volume</font></strong></td>
<td style="text-align: center;">强度误差 = （预报均值-观测均值）<sup>2</sup></td>
<td style="text-align: left;">要素原始单位的平方</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>MSE_pattern</font></strong></td>
<td style="text-align: center;">形态误差 = （MSE_shift_rotate_left-MSE_volume）</td>
<td style="text-align: left;">要素原始单位的平方</td>
</tr>
</tbody>
</table>
<p>在CRA算法中采用最小的能够覆盖grd_ob、grd_fo、grd_fo_shift、grd_fo_shift_rotate非0值的矩形区域内的格点数作为各项平均值以及均方差统计的样本数。
如果将预报(fo)、平移后的预报(fo_shift)、平移旋转后的预报(fo_shift_rotate)以及观测(ob)看做是相空间中的一个点，下面通过一张示意图大致说明CRA输出的不同的均方差的含义。</p>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=a9c5b312d479c7bd550c69793734b319&amp;file=file.jpg" alt="" /></p>
<p>上面的示意图中MSE_shift 并不等于矢量fo->fo_shift的长度的平方，MSE_rotate 也并不等于矢量fo_shift->fo_shift_rotate的长度的平方，因此上面的示意图是不严谨的，它只是大致示意。</p>
<p><strong>调用示例</strong> </p>
<pre><code class="language-python"># 首先读取预报和观测数据
grid1 = meb.grid([100,120,0.05],[24,40,0.05])
filename_ob = r'H:\test_data\input\mem\mode\ob\rain03\20072611.000.nc'
filename_fo = r'H:\test_data\input\mem\mode\ec\rain03\20072608.003.nc'
grd_ob =meb.read_griddata_from_nc(filename_ob,grid = grid1)
grd_fo = meb.read_griddata_from_nc(filename_fo,grid = grid1)</code></pre>
<pre><code class="language-python"># 调用MODE算法中的目标识别和匹配算法
look_featureFinder = mem.mode.feature_finder(grd_ob,grd_fo,smooth=1,threshold=5,minsize=30,compare=&quot;&gt;&quot;)
look_match = mem.mode.centmatch(look_featureFinder)
look_merge = mem.mode.merge_force(look_match)
mem.mode.plot_label(look_merge) #绘图查看目标匹配情况</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=0e5c6147beb12d47b8bf11aba6e30328&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">#调用CRA中的craer算法，对每个不进行刚体变换，以及误差分量计算
look_cra = mem.cra.craer(look_merge,stages = True,translate = True,rotate = True)
print(&quot;CRA 算法执行完毕&quot;)</code></pre>
<pre><code>CRA 算法执行完毕</code></pre>
<pre><code class="language-python"># 查看craer返回结果的第一级关键词，也就是此前成功匹配的目标的编号
print(look_cra.keys())</code></pre>
<pre><code>dict_keys([1, 2, 3, 4])</code></pre>
<pre><code class="language-python">rire = look_cra[2] #选取其中一个目标的cra检验结果
print(&quot;质心经度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lon&quot;])+&quot;°&quot;)
print(&quot;质心纬度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lat&quot;])+&quot;°&quot;)
print(&quot;目标倾角偏差（预报-观测）：&quot; +str(rire[&quot;delta_angle&quot;])+&quot;°&quot;)
print(&quot;观测均值：&quot; +str(rire[&quot;ob_mean&quot;]) +&quot;mm&quot;)
print(&quot;预报均值：&quot; +str(rire[&quot;fo_mean&quot;]) +&quot;mm&quot;)
print(&quot;总误差：&quot; +str(math.sqrt(rire[&quot;MSE_total&quot;])) +&quot;mm&quot;)
print(&quot;平移场残差：&quot; +str(math.sqrt(rire[&quot;MSE_shift_left&quot;])) +&quot;mm&quot;)
print(&quot;平移误差：&quot; +str(math.sqrt(rire[&quot;MSE_shift&quot;])) +&quot;mm&quot;)
print(&quot;平移旋转场残差：&quot; +str(math.sqrt(rire[&quot;MSE_shift_rotate_left&quot;])) +&quot;mm&quot;)
print(&quot;旋转误差：&quot; +str(math.sqrt(rire[&quot;MSE_rotate&quot;])) +&quot;mm&quot;)
print(&quot;强度误差 ：&quot; +str(math.sqrt(rire[&quot;MSE_volume&quot;])) +&quot;mm&quot;)
print(&quot;形态误差：&quot; +str(math.sqrt(rire[&quot;MSE_pattern&quot;])) +&quot;mm&quot;)</code></pre>
<pre><code>质心经度偏差（预报-观测）：-0.7214864249375318°
质心纬度偏差（预报-观测）：0.02260510632156949°
目标倾角偏差（预报-观测）：16.71712676149444°
观测均值：1.6832297411063535mm
预报均值：2.01499193890855mm
总误差：4.587094559638624mm
平移场残差：3.255281683669486mm
平移误差：3.231807181598589mm
平移旋转场残差：3.1242793978190972mm
旋转误差：0.9141865697969876mm
强度误差 ：0.33176219780219673mm
形态误差：3.1066148135464617mm</code></pre>
<pre><code class="language-python">mem.cra.plot_value(look_cra,1)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=5f1be47b7eb084b695c50868f08c34d5&amp;file=file.png" alt="" /></p>
<p><font face="黑体" color=Blue size=4><strong>operation(grd_ob,grd_fo,smooth,threshold,minsize,compare = ">=", stages = True,translate = True,rotate = True)</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=3>grd_ob</font></strong></td>
<td style="text-align: left;">二维矩阵形式的网格观测场</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=3>grd_fo</font></strong></td>
<td style="text-align: left;">二维矩阵形式的网格预报场</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=3>smooth</font></strong></td>
<td style="text-align: left;">平滑系数,采用圆盘形的卷积内核对网格数据进行平滑，smooth相当圆盘的半径</td>
<td>可以为单个整数或列表形式，当其为单个整数时，观测和预报采用相同的参数，否则相应采用不同的参数，如smooth = [5, 4]。平滑系数为0时不做平滑</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=3>threshold</font></strong></td>
<td style="text-align: left;">雨量阈值，平滑后低于阈值的部分会被置0</td>
<td>可以为单个实数或列表形式，当其为单个实数时，观测和预报采用相同的参数，否则相应采用不同的参数，，如threshold = [9,10]</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=3>minsize</font></strong></td>
<td style="text-align: left;">目标的最小面积,面积小于该阈值的目标将会被删除</td>
<td>可以为单个整数或列表形式，当其为单个整数时，观测和预报采用相同的参数，否则相应采用不同的参数，<br>如minsize = [5, 4]，此列表中的数值表示的是格点数。</td>
</tr>
<tr>
<td style="text-align: left;"><strong>compare</strong></td>
<td style="text-align: left;">比较方法，可选项包括">=",">","<=","<"</td>
<td>当关注目标是要素场大于阈值时,采用默认参数，如识别暴雨目标时选择参数">=",若关注对象是要素值小于阈值时，则可选取"<="，如根据能见度场识别大雾目标</td>
</tr>
<tr>
<td style="text-align: left;"><strong>stages</strong></td>
<td style="text-align: left;">布尔值，该参数为True时，采用方案1同时计算最优的平移和旋转，否作采用方案2</td>
</tr>
<tr>
<td style="text-align: left;"><strong>translate</strong></td>
<td style="text-align: left;">stages = False时有效，当translate为True时代表会计算最优平移，否则不计算</td>
</tr>
<tr>
<td style="text-align: left;"><strong>rotate</strong></td>
<td style="text-align: left;">stages = False时有效，当rotate为True时代表会计算最优旋转，否则不计算</td>
</tr>
</tbody>
</table>
<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=green size=3>merge</font></strong></td>
<td style="text-align: center;">look_merge,MODE算法完成目标识别匹配和合并后的输出结果</td>
<td style="text-align: left;">字典</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>cra</font></strong></td>
<td style="text-align: center;">look_cra, craer 函数的返回结果</td>
<td style="text-align: left;">字典</td>
</tr>
</tbody>
</table>
<pre><code class="language-python">result = mem.cra.operation(grd_ob,grd_fo,smooth=1,threshold=5,minsize=30,compare=&quot;&gt;&quot;)</code></pre>
<pre><code class="language-python">result[&quot;cra&quot;][1]</code></pre>
<pre><code>{'grd_fo': &lt;xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 321, lon: 401)&gt;
array([[[[[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]]]]])
Coordinates:
* member (member) &lt;U3 'FST'
* level (level) float64 0.0
* time (time) datetime64[ns] 2020-07-26T11:00:00
* dtime (dtime) int32 0
* lat (lat) float64 24.0 24.05 24.1 24.15 24.2 ... 39.85 39.9 39.95 40.0
* lon (lon) float64 100.0 100.0 100.1 100.2 ... 119.8 119.9 120.0 120.0,
'grd_ob': &lt;xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 321, lon: 401)&gt;
array([[[[[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]]]]])
Coordinates:
* member (member) &lt;U3 'OBS'
* level (level) float64 0.0
* time (time) datetime64[ns] 2020-07-26T11:00:00
* dtime (dtime) int32 0
* lat (lat) float64 24.0 24.05 24.1 24.15 24.2 ... 39.85 39.9 39.95 40.0
* lon (lon) float64 100.0 100.0 100.1 100.2 ... 119.8 119.9 120.0 120.0,
'grd_fo_shift': &lt;xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 321, lon: 401)&gt;
array([[[[[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]]]]])
Coordinates:
* member (member) &lt;U9 'FST_SHIFT'
* level (level) float64 0.0
* time (time) datetime64[ns] 2020-07-26T11:00:00
* dtime (dtime) int32 0
* lat (lat) float64 24.0 24.05 24.1 24.15 24.2 ... 39.85 39.9 39.95 40.0
* lon (lon) float64 100.0 100.0 100.1 100.2 ... 119.8 119.9 120.0 120.0,
'grd_fo_shift_rotate': &lt;xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 321, lon: 401)&gt;
array([[[[[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]]]]])
Coordinates:
* member (member) &lt;U16 'FST_SHIFT_ROTATE'
* level (level) float64 0.0
* time (time) datetime64[ns] 2020-07-26T11:00:00
* dtime (dtime) int32 0
* lat (lat) float64 24.0 24.05 24.1 24.15 24.2 ... 39.85 39.9 39.95 40.0
* lon (lon) float64 100.0 100.0 100.1 100.2 ... 119.8 119.9 120.0 120.0,
'delta_lon': -1.423332167295316,
'delta_lat': -0.571371766284003,
'delta_angle': -7.684624380852102,
'rmse_before': 3.875106574500946,
'rmse_after': 3.3201910574660793,
'ob_mean': 2.4824262196278006,
'fo_mean': 1.8856434146584844,
'MSE_total': 48.21964238147072,
'MSE_shift_left': 44.24292831839233,
'MSE_shift': 3.97671406307839,
'MSE_shift_rotate_left': 35.39833491334686,
'MSE_rotate': 8.844593405045472,
'MSE_volume': 0.35614971630704495,
'MSE_pattern': 35.042185197039814}</code></pre>
<pre><code class="language-python"></code></pre>