几何特征检验
<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
import math</code></pre>
<h1>几何特征检验</h1>
<p>计算和对比观测预报落区的面积、数目和周长等特征</p>
<p>参考文献</p>
<ul>
<li>AghaKouchak, A., Nasrollahi, N., Li, J., Imam, B. and Sorooshian, S. (2011) Geometrical characterization of precipitation patterns. J. Hydrometeorology, 12, 274–285, doi:10.1175/2010JHM1298.1.</li>
</ul>
<p>以下构建两组简单的落区的示例，已对函数功能和效果进行说明</p>
<pre><code class="language-python">grid1 = meb.grid([100, 120, 0.05], [24, 40, 0.05])
dat_ob1 = np.zeros((grid1.nlat,grid1.nlon))
for j in range(grid1.nlat):
for i in range(grid1.nlon):
dat_ob1[j,i] =20* math.exp(-0.001*(i-200)**2 - 0.001*(j-200)**2)
grd_ob1 = meb.grid_data(grid1,dat_ob1)
grid1 = meb.grid([100, 120, 0.05], [24, 40, 0.05])
dat_fo1 = np.zeros((grid1.nlat,grid1.nlon))
for j in range(grid1.nlat):
for i in range(grid1.nlon):
dat_fo1[j,i] =10* math.exp(-0.001*(i-240)**2 - 0.0003*(j-240)**2)
dat_fo1[j,i] +=10* math.exp(-0.001*(i-170)**2 - 0.0003*(j-170)**2)
grd_fo1= meb.grid_data(grid1,dat_fo1)
look_ff1 = mem.mode.feature_finder(grd_ob1,grd_fo1,smooth = 0,threshold = 5,minsize = 1)
mem.mode.plot_label(look_ff1)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=420d1750c7252b21de1e80f373c8cf9e&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">grid2 = meb.grid([100, 120, 0.05], [24, 40, 0.05])
dat_ob2 = np.zeros((grid2.nlat,grid2.nlon))
for j in range(grid2.nlat):
for i in range(grid2.nlon):
dat_ob2[j,i] =20* math.exp(-0.001*(i-200)**2 - 0.001*(j-200)**2)
grd_ob2 = meb.grid_data(grid2,dat_ob2)
grid2 = meb.grid([100, 120, 0.05], [24, 40, 0.05])
dat_fo2 = np.zeros((grid2.nlat,grid2.nlon))
for j in range(grid2.nlat):
for i in range(grid2.nlon):
dat_fo2[j,i] =10* math.exp(-0.001*(i-230)**2 - 0.0003*(j-230)**2)
dat_fo2[j,i] +=10* math.exp(-0.001*(i-170)**2 - 0.0003*(j-170)**2)
grd_fo2= meb.grid_data(grid2,dat_fo2)
look_ff2 = mem.mode.feature_finder(grd_ob2,grd_fo2,smooth = 0,threshold = 5,minsize = 1)
mem.mode.plot_label(look_ff2)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=fa718ae6559c28f0688a9dd231f40eb7&amp;file=file.png" alt="" /></p>
<h2>面积指数</h2>
<p>Aindex = A/Aconvex, A 代表图形面积, and Aconvex 代表图形凸包面积， 其值在0到1之间，越接近1则图形越规则，相反则越离散</p>
<p><font face="黑体" color=Blue size=3><strong>area_index(grd_ob,grd_fo,thresholds)</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><font face="黑体" color=Blue size=5>thresholds</font></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>area_index</font></strong></td>
<td style="text-align: left;">面积指数</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>total_area</font></strong></td>
<td style="text-align: left;">面积和</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>total_convex_area</font></strong></td>
<td style="text-align: left;">最大凸包面积</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">result = mem.space.area_index(grd_ob1,grd_fo1,thresholds=5)
print(result)</code></pre>
<pre><code>{'area_index': [0.9661953159522757, 0.961511708787387], 'total_area': [4373, 8294], 'total_convex_area': [4526.0, 8626.0]}</code></pre>
<pre><code class="language-python">result = mem.space.area_index(grd_ob2,grd_fo2,thresholds=5)
print(result)</code></pre>
<pre><code>{'area_index': [0.9661953159522757, 0.8596507352941176], 'total_area': [4373, 9353], 'total_convex_area': [4526.0, 10880.0]}</code></pre>
<h2>数量指数</h2>
<p>count_index = 1 - (NC - 1)/(sqrt(NP) + NC), NP 代表图像中非0的格点数， NC 代表连通域数量， 其值在0到1之间，越接近1则图形连通程度越高，相反则越少</p>
<p><font face="黑体" color=Blue size=3><strong>count_index(grd_ob,grd_fo,thresholds)</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><font face="黑体" color=Blue size=5>thresholds</font></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>count_index</font></strong></td>
<td style="text-align: left;">数量指数</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>count_nonzero</font></strong></td>
<td style="text-align: left;">达到阈值的格点数</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>count_connect</font></strong></td>
<td style="text-align: left;">联通的目标数</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">result = mem.space.count_index(grd_ob1,grd_fo1,thresholds=5)
print(result)</code></pre>
<pre><code>{'count_index': [1.0151220357808053, 1.0], 'count_nonzero': [4373, 8294], 'count_connect': [1, 2]}</code></pre>
<pre><code class="language-python">result = mem.space.count_index(grd_ob2,grd_fo2,thresholds=5)
print(result)</code></pre>
<pre><code>{'count_index': [1.0151220357808053, 1.010340095094156], 'count_nonzero': [4373, 9353], 'count_connect': [1, 1]}</code></pre>
<h2>周长指数</h2>
<p>Sindex = Pmin/P, n代表值为正数的格点数量, 如果 floor(sqrt(n)) = sqrt(n)，Pmin = 4<em>sqrt(n) ，否则 Pmin = 2 </em> floor(2*sqrt(n)+1). P 是非0格点的周长， 其值在0到1之间，越接近1则图形越接近正方形。因为网格数据是离散的，所以计算周长是是围绕格点的折线。</p>
<p><font face="黑体" color=Blue size=3><strong>surround_index(grd_ob,grd_fo,thresholds)</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><font face="黑体" color=Blue size=5>thresholds</font></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>surround_index</font></strong></td>
<td style="text-align: left;">周长指数</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>perimeter_min</font></strong></td>
<td style="text-align: left;">同样面积的落区的最小可能周长</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=green size=3>perimeter</font></strong></td>
<td style="text-align: left;">落区的实际周长</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">result = mem.space.surround_index(grd_ob1,grd_fo1,thresholds=5)
print(result)</code></pre>
<pre><code>{'surround_index': [0.8866666666666667, 0.6310344827586207], 'perimeter_min': [266, 366], 'perimeter': [300, 580]}</code></pre>
<pre><code class="language-python">result = mem.space.surround_index(grd_ob2,grd_fo2,thresholds=5)
print(result)</code></pre>
<pre><code>{'surround_index': [0.8866666666666667, 0.7185185185185186], 'perimeter_min': [266, 388], 'perimeter': [300, 540]}</code></pre>
<p>以下是实际示例中的应用效果</p>
<pre><code class="language-python">path_ob = r'H:\test_data\input\mem\mode\ob\rain03\20072611.000.nc'
path_fo = r'H:\test_data\input\mem\mode\ec\rain03\20072608.003.nc'
grd_ob3 = meb.read_griddata_from_nc(path_ob, grid=grid1, time=&quot;2020072611&quot;, dtime=0, data_name=&quot;OBS&quot;)
grd_fo3 = meb.read_griddata_from_nc(path_fo, grid=grid1, time=&quot;2020072608&quot;, dtime=3, data_name=&quot;ECMWF&quot;)
look_ff = mem.mode.feature_finder(grd_ob3,grd_fo3,smooth = 0,threshold = 5,minsize = 1)
mem.mode.plot_label(look_ff)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=e4b4ccb07710c58f5cde8a9234e94067&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">result = mem.space.area_index(grd_ob3,grd_fo3,thresholds=5)
print(result)</code></pre>
<pre><code>{'area_index': [0.6111037934668072, 0.663979336255479], 'total_area': [9279, 8483], 'total_convex_area': [15184.0, 12776.0]}</code></pre>
<pre><code class="language-python">result = mem.space.count_index(grd_ob3,grd_fo3,thresholds=5)
print(result)</code></pre>
<pre><code>{'count_index': [0.8742449827620725, 0.8159560866765792], 'count_nonzero': [9279, 8483], 'count_connect': [16, 23]}</code></pre>
<pre><code class="language-python">result = mem.space.surround_index(grd_ob3,grd_fo3,thresholds=5)
print(result)</code></pre>
<pre><code>{'surround_index': [0.2769010043041607, 0.32062391681109187], 'perimeter_min': [386, 370], 'perimeter': [1394, 1154]}</code></pre>