rigider(刚体变换)
<p>[TOC]</p>
<h1>刚体变换方法概述</h1>
<p><font face="黑体" color=Black size=5>rigider方法是对于一个空间场计算其最优刚体变换参数 </font> </p>
<pre><code class="language-python">%matplotlib inline
%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib.pyplot as plt
import meteva.method as mem
import meteva.base as meb
import math</code></pre>
<h1>简单刚体变换</h1>
<p><font face="黑体" color=Blue size=3><strong>rigid_simple(grd_ob,grd_fo)</strong></font> </p>
<p>简单刚体变换的计算方法：</p>
<ul>
<li>步骤1： 计算观测场的质心位置和倾角</li>
<li>步骤2： 计算预报场的质心位置和倾角</li>
<li>步骤3： 根据质心位置和倾角差，对预报场进行平移和旋转，之后重新插值到原始网格上</li>
<li>步骤4： 返回位置和倾角偏差以及变换后的预报场。</li>
</ul>
<table>
<thead>
<tr>
<th style="text-align: left;">参数</th>
<th style="text-align: center;">说明</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>
</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>
</tr>
</tbody>
</table>
<p><strong>——————————————————————————————————————————————</strong>
<font face="黑体" color=green size=4><strong>返回结果内容说明</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>
</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>
</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>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>rmse_before</font></strong></td>
<td style="text-align: center;">刚体变换前，预报和观测场的均方根误差，单位和原始场一致</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>rmse_after</font></strong></td>
<td style="text-align: center;">刚体变换后，预报和观测场的均方根误差，单位和原始场一致</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python"># 生成方形的测试数据
grid_ob = meb.grid([100,130,1],[15,45,1],gtime = [&quot;2021080408&quot;],dtime_list = [0],member_list=[&quot;OBS&quot;])
x = np.zeros((31, 31))
x[8:10, 10:17] = 10
grd_ob = meb.grid_data(grid_ob,x)
grid_fo = meb.grid([100,130,1],[15,45,1],gtime = [&quot;2021080308&quot;],dtime_list = [24],member_list=[&quot;ECMWF&quot;])
y = np.zeros((31, 31))
for i in range(y.shape[0]):
for j in range(y.shape[1]):
a = (i - 15) * math.sin(math.pi/4) - (j - 15) * math.cos(math.pi/4)
b = (i - 15) * math.sin(math.pi/4) + (j - 15) * math.cos(math.pi/4)
if abs(a) &lt; 2 and abs(b)&lt;5:
y[i,j] = 5
grd_fo = meb.grid_data(grid_fo,y)
#调用刚体变换函数
rire = mem.rigider.rigid_simple(grd_ob, grd_fo) </code></pre>
<pre><code class="language-python">print(&quot;质心经度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lon&quot;]))
print(&quot;质心纬度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lat&quot;]))
print(&quot;质心角度偏差（预报-观测）：&quot; +str(rire[&quot;delta_angle&quot;]))
print(&quot;变换前rmse：&quot; +str(rire[&quot;rmse_before&quot;]))
print(&quot;变换后rmse：&quot; +str(rire[&quot;rmse_after&quot;]))</code></pre>
<pre><code>质心经度偏差（预报-观测）：2.0
质心纬度偏差（预报-观测）：6.5
质心角度偏差（预报-观测）：-45.00000076761726
变换前rmse：1.5554275420956378
变换后rmse：0.9174346851228268</code></pre>
<h1>绘制变换前后要素场</h1>
<p><font face="黑体" color=Blue size=3><strong>plot_value(rire,cmap="rain_1h",clevs = None)</strong></font> </p>
<p>绘制变换前后要素场，方便查看变换是否合理</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>rire</font></strong></td>
<td style="text-align: left;">目标识别，匹配或合并函数所得结果</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>cmap</font></strong></td>
<td style="text-align: left;">网格预报数据绘图colorbar的cmap参数，字符或者matplotlib 的cm对象。 缺省值”rainbow”对应彩虹色方案。 当cmap为字符时支持两种情况，一种是matplotlib 能够识别的字符串，例如”rainbow”,”bwr”等等，具体可以参考https://matplotlib.org/tutorials/colors/colormaps.html， 另外则是meteva定义的一些和预报检验有关的一些配色方案，具体可参考meb.def_cmap_clevs函数说明</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>clevs</font></strong></td>
<td style="text-align: left;">网格数据绘图colorbar的等级参数， 取值为列表或一维nump数值</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>result</font></strong></td>
<td style="text-align: left;">无</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">mem.rigider.plot_value(rire)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/62791292b80eac274e9fa5da321ac4c6" alt="" /></p>
<p>从上面的结果可以看出刚体变换只能移动位置和角度，对雨带长宽和强度没有调整。</p>
<h1>最优刚体变换</h1>
<p><font face="黑体" color=Blue size=3><strong>rigid_optimal(grd_ob,grd_fo,stages=True, translate=True, rotate=True, show=False)</strong></font> </p>
<p>最优刚体变换的计算方法：</p>
<ul>
<li>首先通过计算观测和预报场的质心位置和倾角，得出两者的距离误差和角度误差</li>
<li>接下来，以上述位置和倾角偏差作为初始值，通过最优化算法计算处使得转换后观测预报之间误差最小的变换方法，最优化的方法又分两种：
<ul>
<li>方案1，同时计算最优平移距离和最优旋转角度，</li>
<li>方案2，分别计算最优平移距离和最优旋转角度，也可以通过参数控制是否只做平移，或只做旋转。</li>
</ul></li>
</ul>
<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>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>
<tr>
<td style="text-align: left;"><strong>show</strong></td>
<td style="text-align: left;">布尔值，是否打印日志</td>
</tr>
</tbody>
</table>
<p><strong>—————————————————————————————————————————————</strong></p>
<p><font face="黑体" color=green size=4><strong>返回结果内容说明</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>
</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>
</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>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>rmse_before</font></strong></td>
<td style="text-align: center;">刚体变换前，预报和观测场的均方根误差，单位和原始场一致</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>rmse_after</font></strong></td>
<td style="text-align: center;">刚体变换后，预报和观测场的均方根误差，单位和原始场一致</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">rire = mem.rigider.rigid_optimal(grd_ob, grd_fo)
print(&quot;质心经度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lon&quot;]))
print(&quot;质心纬度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lat&quot;]))
print(&quot;质心角度偏差（预报-观测）：&quot; +str(rire[&quot;delta_angle&quot;]))
print(&quot;变换前rmse：&quot; +str(rire[&quot;rmse_before&quot;]))
print(&quot;变换后rmse：&quot; +str(rire[&quot;rmse_after&quot;]))
mem.rigider.plot_value(rire)</code></pre>
<pre><code>质心经度偏差（预报-观测）：2.0
质心纬度偏差（预报-观测）：6.5
质心角度偏差（预报-观测）：-45.00000076761726
变换前rmse：1.5554275420956378
变换后rmse：0.9174346851228268</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/7b094d9125b58c986c9c3da606047bac" alt="" /></p>
<pre><code class="language-python">rire = mem.rigider.rigid_optimal(grd_ob, grd_fo,stages=False,rotate=False)
print(&quot;质心经度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lon&quot;]))
print(&quot;质心纬度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lat&quot;]))
print(&quot;质心角度偏差（预报-观测）：&quot; +str(rire[&quot;delta_angle&quot;]))
print(&quot;变换前rmse：&quot; +str(rire[&quot;rmse_before&quot;]))
print(&quot;变换后rmse：&quot; +str(rire[&quot;rmse_after&quot;]))
mem.rigider.plot_value(rire)</code></pre>
<pre><code>质心经度偏差（预报-观测）：2.4181300313541025
质心纬度偏差（预报-观测）：6.581870001061175
质心角度偏差（预报-观测）：0
变换前rmse：1.5554275420956378
变换后rmse：1.1038059967760732</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/caabc55f1282a8813ebcd7a66c729862" alt="" /></p>
<pre><code class="language-python">rire = mem.rigider.rigid_optimal(grd_ob, grd_fo,stages=False,translate=False)
print(&quot;质心经度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lon&quot;]))
print(&quot;质心纬度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lat&quot;]))
print(&quot;质心角度偏差（预报-观测）：&quot; +str(rire[&quot;delta_angle&quot;]))
print(&quot;变换前rmse：&quot; +str(rire[&quot;rmse_before&quot;]))
print(&quot;变换后rmse：&quot; +str(rire[&quot;rmse_after&quot;]))
mem.rigider.plot_value(rire)</code></pre>
<pre><code>质心经度偏差（预报-观测）：0
质心纬度偏差（预报-观测）：0
质心角度偏差（预报-观测）：0
变换前rmse：1.5554275420956378
变换后rmse：1.5086187503344448</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/c8faaf19fcf9e78ec856d30a075c7d32" alt="" /></p>
<pre><code class="language-python">#基于实际降水数据的测试
grid0 = meb.grid([100,115,0.25],[25,35,0.25])
file_name_ob = r'H:\test_data\input\mem\rigider/ob/20070111.000.nc' #读取观测场
grd_ob_rain = meb.read_griddata_from_nc(file_name_ob,grid = grid0) # 读取预报场
file_name_fo = r'H:\test_data\input\mem\rigider/fo/20070108.003.nc'
grd_fo_rain = meb.read_griddata_from_nc(file_name_fo,grid = grid0)
rire = mem.rigider.rigid_optimal(grd_ob_rain, grd_fo_rain) #基于最优刚体变换
print(&quot;质心经度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lon&quot;]))
print(&quot;质心纬度偏差（预报-观测）：&quot; +str(rire[&quot;delta_lat&quot;]))
print(&quot;质心角度偏差（预报-观测）：&quot; +str(rire[&quot;delta_angle&quot;]))
print(&quot;变换前rmse：&quot; +str(rire[&quot;rmse_before&quot;]))
print(&quot;变换后rmse：&quot; +str(rire[&quot;rmse_after&quot;]))
mem.rigider.plot_value(rire) #绘制刚体变换前后结果对比图</code></pre>
<pre><code>质心经度偏差（预报-观测）：-1.6612719002560123
质心纬度偏差（预报-观测）：0.24661751308289914
质心角度偏差（预报-观测）：3.4687474644104235
变换前rmse：4.956656473575831
变换后rmse：4.750256269016443</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/b222d990498e8584c3593ea45c79e3e2" alt="" /></p>