插值
<p>[TOC]</p>
<pre><code class="language-python">%matplotlib inline
%load_ext autoreload
%autoreload 2
import meteva.base as meb
import numpy as np</code></pre>
<h1>格点—>站点 临近点插值</h1>
<p><strong><font face="黑体" color=blue size = 5>interp_gs_nearest(grd,sta,used_coords = "xy")</font></strong><br />
将水平网格数据插值到平面的离散站点上，对于每一个站点，其取值设置为周围四个网格点中最近的一个格点的取值。 </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</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a></td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>sta</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a></td>
</tr>
<tr>
<td style="text-align: left;"><strong>used_coords</strong></td>
<td style="text-align: left;">插值操作使用的维度，缺省情况下插值操作只在水平方向实现，此时返回的站点数据中id、lon、lat三列的取值会和sta一致，经纬度超出网格范围的站点将被删除，而level,time,dtime参数采用grd里的坐标信息。目前，该函数仅支持该参数缺省时的功能，该参数为其它选项时对应的功能待完善中。</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a>,其站点列表由sta处于grd网格范围内的站点列表确定</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">grid0 = meb.grid([100,102,1],[20,22,1],gtime=[&quot;2019013108&quot;],dtime_list = [24],level_list = [500,700],member_list = [&quot;GRAPES&quot;])
x= np.arange(3)
y= np.arange(3)
dat = np.array(np.meshgrid(x,y))
grd = meb.grid_data(grid0,dat) #根据网格信息和numpy数组生成网格数
print(grd)</code></pre>
<pre><code>&lt;xarray.DataArray 'data0' (member: 1, level: 2, time: 1, dtime: 1, lat: 3, lon: 3)&gt;
array([[[[[[0, 1, 2],
[0, 1, 2],
[0, 1, 2]]]],
[[[[0, 0, 0],
[1, 1, 1],
[2, 2, 2]]]]]])
Coordinates:
* member (member) &lt;U6 'GRAPES'
* level (level) int32 500 700
* time (time) datetime64[ns] 2019-01-31T08:00:00
* dtime (dtime) int32 24
* lat (lat) int32 20 21 22
* lon (lon) int32 100 101 102</code></pre>
<pre><code class="language-python">station = meb.read_station(meb.station_国家站)
meb.set_stadata_coords(station,level = 600)
print(station)</code></pre>
<pre><code> level time dtime id lon lat data0
0 600 2099-01-01 08:00:00 0 50136 122.52 52.97 0
1 600 2099-01-01 08:00:00 0 50137 122.37 53.47 0
2 600 2099-01-01 08:00:00 0 50246 124.72 52.35 0
3 600 2099-01-01 08:00:00 0 50247 123.57 52.03 0
4 600 2099-01-01 08:00:00 0 50349 124.40 51.67 0
... ... ... ... ... ... ... ...
2406 600 2099-01-01 08:00:00 0 59945 109.70 18.65 0
2407 600 2099-01-01 08:00:00 0 59948 109.58 18.22 0
2408 600 2099-01-01 08:00:00 0 59951 110.33 18.80 0
2409 600 2099-01-01 08:00:00 0 59954 110.03 18.55 0
2410 600 2099-01-01 08:00:00 0 59981 112.33 16.83 0
[2411 rows x 7 columns]</code></pre>
<pre><code class="language-python">sta = meb.interp_gs_nearest(grd,station)
print(sta) #used_coords缺省时，返回数据的level,time和dtime和grd里取值一致,由于level有两层，因此插值结果中每个id会出现两次 </code></pre>
<pre><code> level time dtime id lon lat GRAPES
1271 500 2019-01-31 08:00:00 24 56958 100.42 21.92 0
1272 500 2019-01-31 08:00:00 24 56959 100.78 22.00 0
1277 500 2019-01-31 08:00:00 24 56969 101.58 21.47 1
1271 700 2019-01-31 08:00:00 24 56958 100.42 21.92 1
1272 700 2019-01-31 08:00:00 24 56959 100.78 22.00 2
1277 700 2019-01-31 08:00:00 24 56969 101.58 21.47 1</code></pre>
<h1>格点—>站点 双线性插值</h1>
<p><strong><font face="黑体" color=blue size = 5>interp_gs_linear(grd,sta,used_coords ="xy")</font></strong><br />
将水平网格数据插值到平面的离散站点上，对于每一个站点，其取值利用周围四个网格点进行双线性插值获得 </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</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a></td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>sta</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a></td>
</tr>
<tr>
<td style="text-align: left;"><strong>used_coords</strong></td>
<td style="text-align: left;">插值操作使用的维度，该参数包括如下选项： 当该参数为"xy"时，插值操作只在水平方向实现，此时返回的站点数据中id、lon、lat三列的取值会和sta一致，经纬度超出网格范围的站点将被删除，而level,time,dtime参数采用grd里的坐标信息；当该参数为"xyz"时，插值操作只在lon,lat,level三个空间维度实现，此时返回的站点数据中id、lon、lat,level四列的取值会和sta一致，经纬度超出网格范围的站点将被删除，而time,dtime参数采用grd里的坐标信息；</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a>,其站点列表由sta处于grd网格范围内的站点列表确定</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">sta = meb.interp_gs_linear(grd,station) #在二维平面进行插值，level维度会
print(sta) #used_coords缺省时，返回数据的level,time和dtime和grd里取值一致,由于level有两层，因此插值结果中每个id会出现两次 </code></pre>
<pre><code> level time dtime id lon lat GRAPES
1271 500 2019-01-31 08:00:00 24 56958 100.42 21.92 0.42
1272 500 2019-01-31 08:00:00 24 56959 100.78 22.00 0.78
1277 500 2019-01-31 08:00:00 24 56969 101.58 21.47 1.58
1271 700 2019-01-31 08:00:00 24 56958 100.42 21.92 1.92
1272 700 2019-01-31 08:00:00 24 56959 100.78 22.00 2.00
1277 700 2019-01-31 08:00:00 24 56969 101.58 21.47 1.47</code></pre>
<pre><code class="language-python">sta = meb.interp_gs_linear(grd,station,used_coords = &quot;xyz&quot;)
print(sta) #used_coords缺省时，返回数据的level,time和dtime和grd里取值一致,由于level有两层，因此插值结果中每个id会出现两次 </code></pre>
<pre><code> level time dtime id lon lat GRAPES
1271 600 2019-01-31 08:00:00 24 56958 100.42 21.92 1.170
1272 600 2019-01-31 08:00:00 24 56959 100.78 22.00 1.390
1277 600 2019-01-31 08:00:00 24 56969 101.58 21.47 1.525</code></pre>
<pre><code class="language-python">sta = meb.interp_gs_linear(grd,station,used_coords = &quot;xydt&quot;)
print(sta) #used_coords缺省时，返回数据的level,time和dtime和grd里取值一致,由于level有两层，因此插值结果中每个id会出现两次 </code></pre>
<pre><code>dtime维度size = 1,无法开展dtime维度插值
None</code></pre>
<h1>格点—>站点 双三次插值</h1>
<p><strong><font face="黑体" color=blue size = 5>interp_gs_cubic(grd,sta,used_coords = "xy")</font></strong><br />
将水平网格数据插值到平面的离散站点上，对于每一个站点，其取值利用周围九个网格点进行双三次插值获得 </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</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a></td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>sta</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a></td>
</tr>
<tr>
<td style="text-align: left;"><strong>used_coords</strong></td>
<td style="text-align: left;">插值操作使用的维度，缺省情况下插值操作只在水平方向实现，此时返回的站点数据中id、lon、lat三列的取值会和sta一致，经纬度超出网格范围的站点将被删除，而level,time,dtime参数采用grd里的坐标信息。目前，该函数仅支持该参数缺省时的功能，该参数为其它选项时对应的功能待完善中。</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a>,其站点列表由sta处于grd网格范围内的站点列表确定</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">sta =meb.interp_gs_cubic(grd,station)
print(sta)</code></pre>
<pre><code> level time dtime id lon lat GRAPES
0 500 2019-01-31 08:00:00 24 56958 100.42 21.92 0.355852
1 500 2019-01-31 08:00:00 24 56959 100.78 22.00 0.745108
2 500 2019-01-31 08:00:00 24 56969 101.58 21.47 1.644148
3 700 2019-01-31 08:00:00 24 56958 100.42 21.92 1.943552
4 700 2019-01-31 08:00:00 24 56959 100.78 22.00 2.000000
5 700 2019-01-31 08:00:00 24 56969 101.58 21.47 1.531029</code></pre>
<h1>格点—>格点 双线性插值</h1>
<p><strong><font face="黑体" color=blue size = 5>interp_gg_linear(grd,grid,used_coords = "xy",outer_value = None)</font></strong><br />
将水平网格数据插值到平面的离散站点上，对于每一个站点，其取值利用周围四个网格点进行双线性插值获得 </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</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a></td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>grid</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格信息类变量</a>,插值的目标网格。</td>
</tr>
<tr>
<td style="text-align: left;"><strong>used_coords</strong></td>
<td style="text-align: left;">插值操作使用的维度，缺省情况下插值操作只在水平方向实现，此时返回的网格数据中lon、lat维度的范围和间距由grid确定，其它维度信息由grd确定</td>
</tr>
<tr>
<td style="text-align: left;"><strong>outer_value</strong></td>
<td style="text-align: left;">当目标网格范围大于原始数据网格时，设置超出部分的取值为outer_value。如果存在网格超出的情况，就必须设置该参数，否则返回None。</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a>，网格范围由grid参数确定</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">grid1 = meb.grid([100,102,0.5],[20,22,0.5])
grd1 = meb.interp_gg_linear(grd,grid1)
print(grd1)</code></pre>
<pre><code>&lt;xarray.DataArray 'data0' (member: 1, level: 2, time: 1, dtime: 1, lat: 5, lon: 5)&gt;
array([[[[[[0. , 0.5, 1. , 1.5, 2. ],
[0. , 0.5, 1. , 1.5, 2. ],
[0. , 0.5, 1. , 1.5, 2. ],
[0. , 0.5, 1. , 1.5, 2. ],
[0. , 0.5, 1. , 1.5, 2. ]]]],
[[[[0. , 0. , 0. , 0. , 0. ],
[0.5, 0.5, 0.5, 0.5, 0.5],
[1. , 1. , 1. , 1. , 1. ],
[1.5, 1.5, 1.5, 1.5, 1.5],
[2. , 2. , 2. , 2. , 2. ]]]]]])
Coordinates:
* member (member) &lt;U6 'GRAPES'
* level (level) int32 500 700
* time (time) datetime64[ns] 2019-01-31T08:00:00
* dtime (dtime) int32 24
* lat (lat) float64 20.0 20.5 21.0 21.5 22.0
* lon (lon) float64 100.0 100.5 101.0 101.5 102.0</code></pre>
<pre><code class="language-python">grid1 = meb.grid([100,103,0.5],[20,23,0.5])
grd1 = meb.interp_gg_linear(grd,grid1)
print(grd1)</code></pre>
<pre><code>当目标网格超出数据网格时，outer_value参数必须赋值
None</code></pre>
<pre><code class="language-python">grid1 = meb.grid([100,103,0.5],[20,23,0.5])
grd1 = meb.interp_gg_linear(grd,grid1,outer_value = -1)
print(grd1)</code></pre>
<pre><code>&lt;xarray.DataArray 'data0' (member: 1, level: 2, time: 1, dtime: 1, lat: 7, lon: 7)&gt;
array([[[[[[ 0. , 0.5, 1. , 1.5, 2. , -1. , -1. ],
[ 0. , 0.5, 1. , 1.5, 2. , -1. , -1. ],
[ 0. , 0.5, 1. , 1.5, 2. , -1. , -1. ],
[ 0. , 0.5, 1. , 1.5, 2. , -1. , -1. ],
[ 0. , 0.5, 1. , 1.5, 2. , -1. , -1. ],
[-1. , -1. , -1. , -1. , -1. , -1. , -1. ],
[-1. , -1. , -1. , -1. , -1. , -1. , -1. ]]]],
[[[[ 0. , 0. , 0. , 0. , 0. , -1. , -1. ],
[ 0.5, 0.5, 0.5, 0.5, 0.5, -1. , -1. ],
[ 1. , 1. , 1. , 1. , 1. , -1. , -1. ],
[ 1.5, 1.5, 1.5, 1.5, 1.5, -1. , -1. ],
[ 2. , 2. , 2. , 2. , 2. , -1. , -1. ],
[-1. , -1. , -1. , -1. , -1. , -1. , -1. ],
[-1. , -1. , -1. , -1. , -1. , -1. , -1. ]]]]]])
Coordinates:
* member (member) &lt;U6 'GRAPES'
* level (level) int32 500 700
* time (time) datetime64[ns] 2019-01-31T08:00:00
* dtime (dtime) int32 24
* lat (lat) float64 20.0 20.5 21.0 21.5 22.0 22.5 23.0
* lon (lon) float64 100.0 100.5 101.0 101.5 102.0 102.5 103.0</code></pre>
<h1>站点—>格点 反距离权重插值</h1>
<p><strong><font face="黑体" color=blue size = 5>interp_sg_idw(sta, grid, background=None, effectR=1000, nearNum=8,decrease = 2):</font></strong>
将站点数据插值到水平网格上 </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>sta</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a>，（暂时只支持sta中只包含单个层次时间时效数据列的情况）</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>grid</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a>，插值的目标网格。</td>
</tr>
<tr>
<td style="text-align: left;"><strong>background</strong></td>
<td style="text-align: left;">插值背景场，在远离站点的区域内将采用background的值，如果background为None，则取为0</td>
</tr>
<tr>
<td style="text-align: left;"><strong>effectR</strong></td>
<td style="text-align: left;">最大的插值半径，单位km</td>
</tr>
<tr>
<td style="text-align: left;"><strong>nearNum</strong></td>
<td style="text-align: left;">插值选择的临近站点的个数, nearNum =1时即为临近点插值</td>
</tr>
<tr>
<td style="text-align: left;"><strong>decrease</strong></td>
<td style="text-align: left;">站点权重随距离幂次衰减，即 站点权重 = 1 /(距离 ** decrease) ,其中decrease是幂函数的指数部分参数</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a>，网格范围由grid参数确定</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">grd2 = meb.interp_sg_idw(sta,grid1,nearNum = 2)
print(grd2)</code></pre>
<pre><code>&lt;xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 7, lon: 7)&gt;
array([[[[[[1.1497020000000016, 1.5875887499999985,
1.5875887499999988, 1.5875887499999985,
1.5875887499999988, 1.5875887499999988,
1.5875887499999988],
[1.1497020000000016, 1.5875887499999985,
1.5875887499999988, 1.5875887499999985,
1.5875887499999985, 1.5875887499999985,
1.5875887499999988],
[1.1497020000000016, 1.1497020000000016,
1.587588749999999, 1.5875887499999988,
1.5875887499999985, 1.5875887499999985,
1.5875887499999988],
[1.1497020000000016, 1.1497020000000016,
1.3725540000000007, 1.5875887499999988,
1.5875887499999985, 1.5875887499999985,
1.5875887499999985],
[1.1497020000000016, 1.1497020000000016,
1.3725540000000007, 1.5875887499999985,
1.5875887499999988, 1.5875887499999985,
1.5875887499999985],
[1.1497020000000016, 1.3725540000000005,
1.3725540000000007, 1.3725540000000007,
1.5875887499999985, 1.5875887499999988,
1.5875887499999988],
[1.1497020000000016, 1.3725540000000007,
1.3725540000000007, 1.3725540000000007,
1.3725540000000007, 1.5875887499999985,
1.5875887499999985]]]]]], dtype=object)
Coordinates:
* member (member) &lt;U6 'GRAPES'
* level (level) int32 500
* time (time) datetime64[ns] 2019-01-31T08:00:00
* dtime (dtime) int32 24
* lat (lat) float64 20.0 20.5 21.0 21.5 22.0 22.5 23.0
* lon (lon) float64 100.0 100.5 101.0 101.5 102.0 102.5 103.0</code></pre>
<pre><code class="language-python">station0 = meb.read_station(meb.station_国家站)
rain24 = meb.read_stadata_from_micaps3(r&quot;H:\test_data\input\meb\rain24.m3.txt&quot;,station = station0)
grid0 = meb.grid([70,140,0.1],[15,55,0.1])
import time
start = time.time()
grd = meb.interp_sg_idw(rain24,grid = grid0)
used_time = time.time() - start
print(&quot;反距离权重插值耗时： &quot; + str(used_time) +&quot;秒&quot;)
meb.plot_tools.contourf_2d_grid(grd)</code></pre>
<pre><code>反距离权重插值耗时： 0.6313118934631348秒</code></pre>
<p><img src="https://www.showdoc.cc/server/api/attachment/visitfile/sign/c20a351da919079bfdb94a15ba329baa?showdoc=.jpg" alt="" /></p>
<h1>站点—>格点 cressman插值</h1>
<p><strong><font face="黑体" color=blue size = 5>interp_sg_cressman(sta0, grid, r_list,background=None , nearNum=100):</font></strong>
将站点数据插值到水平网格上 </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>sta</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a>，（暂时只支持sta中只包含单个层次时间时效数据列的情况）</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>grid</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a>，插值的目标网格。</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>r_list</font></strong></td>
<td style="text-align: left;">插值的影响半径序列，由大到小排列，单位为km</td>
</tr>
<tr>
<td style="text-align: left;"><strong>background</strong></td>
<td style="text-align: left;">插值背景场，在距离最近的站点的距离超过r_list[0]的那些格点将采用background网格数中对应的网格值，如果background为None，则取为0</td>
</tr>
<tr>
<td style="text-align: left;"><strong>nearNum</strong></td>
<td style="text-align: left;">插值时所用的最大临近点个数。例如当 r_list = [5000,100]时，则在第一轮插值时基本上做插值时全国的站点（例如2500个）都会用到一个格点的权重计算当中，但这样计算效率会很低，当nearNum = 100时,即使用50000km的半径扫描到2500个站点，实际插值时仅保留100个离格点最近的站点，由此提高计算效率，当半径收缩后，半径内的站点数小于nearNum时，则所有半径范围内的站点都会参与插值</td>
</tr>
<tr>
<td style="text-align: left;"><strong>decrease</strong></td>
<td style="text-align: left;">站点权重随距离幂次衰减，即 站点权重 = 1 /(距离 ** decrease) ,其中decrease是幂函数的指数部分参数</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a>，网格范围由grid参数确定</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">start = time.time()
grd1 = meb.interp_sg_cressman(rain24,grid = grid0,r_list = [1000,200,100,50],nearNum = 100)
used_time = time.time() - start
print(&quot;cressman插值耗时： &quot; + str(used_time) +&quot;秒&quot;)
meb.plot_tools.contourf_2d_grid(grd1)</code></pre>
<pre><code>cressman插值耗时： 6.56514310836792秒</code></pre>
<p><img src="https://www.showdoc.cc/server/api/attachment/visitfile/sign/26ab731a3e75ee0c83e98ab1ae7bc06a?showdoc=.jpg" alt="" /></p>
<pre><code class="language-python">start = time.time()
grd1 = meb.interp_sg_cressman(rain24,grid = grid0,r_list = [1000,200,100,50],nearNum =200)
used_time = time.time() - start
print(&quot;cressman插值耗时： &quot; + str(used_time) +&quot;秒&quot;)
meb.plot_tools.contourf_2d_grid(grd1)</code></pre>
<pre><code>cressman插值耗时： 13.275934934616089秒</code></pre>
<p><img src="https://www.showdoc.cc/server/api/attachment/visitfile/sign/52bbf137fed588c62e60f4bfd697bf2c?showdoc=.jpg" alt="" /></p>
<p>通过上面的例子可以看出，在基于全国2500个站进行降水插值时，nearNum参数从100增大到200对国内的插值效果影响甚微，但耗时增加了近一倍。 </p>
<h1>站点—>站点 反距离权重插值</h1>
<p><strong><font face="黑体" color=blue size = 5>interp_ss_idw(sta, station, background=None, effectR=1000, nearNum=8):</font></strong>
将站点数据插值到另一套站点上 </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>sta</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a>，（暂时只支持sta中只包含单个层次时间时效数据列的情况）</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>station</font></strong></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a>，插值的目标站点列表</td>
</tr>
<tr>
<td style="text-align: left;"><strong>background</strong></td>
<td style="text-align: left;">插值背景场，在远离站点的区域内将采用background的值，如果background为None，则取为0</td>
</tr>
<tr>
<td style="text-align: left;"><strong>effectR</strong></td>
<td style="text-align: left;">最大的插值半径，单位km</td>
</tr>
<tr>
<td style="text-align: left;"><strong>nearNum</strong></td>
<td style="text-align: left;">插值选择的临近站点的个数, nearNum ==1时即为临近点插值</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600580125986">站点数据</a>，站点列表和station完全一致</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">sta_g = meb.trans_grd_to_sta(grd2) #将网格数据转换成站点形式
station_500 = meb.sele_by_para(sta,level = 500)
sta2= meb.interp_ss_idw(sta_g,station = station_500,nearNum = 4) #从站点数据sta_g 插值到 站点列表sta上
print(sta2)</code></pre>
<pre><code> level time dtime id lon lat GRAPES
0 500 2019-01-31 08:00:00 24 56958 100.42 21.92 1.157298
1 500 2019-01-31 08:00:00 24 56959 100.78 22.00 1.300249
2 500 2019-01-31 08:00:00 24 56969 101.58 21.47 1.587589</code></pre>
<h1>模式面插值到等压面</h1>
<p><strong><font face="黑体" color=blue size = 5>interp_zlevel_to_plevel_linear(pressure,element_interp,level_list,with_nan = False)</font></strong> </p>
<p>采用线性插值将模式面的气象要素场插值到等压面。对于每个等压面格点，根据它的气压值p<em>找到模式面的气压值比它大的最近层气压值p1和比它小的最近层的气压值p0， 则等压面格点的插值结果 v</em> = (v1 <em> (p_ - p0) + v2 </em> (p1 - p_)) / （p1-p0)， 其中 v1和v0 分别是两各模式层格点的被插值前的要素值。</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>pressure</font></strong></td>
<td style="text-align: left;">模式面的气压数据，<a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a></td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>element_interp</font></strong></td>
<td style="text-align: left;">模式面的待插值场数据，网格和pressure需完全一致</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>level_list</font></strong></td>
<td style="text-align: left;">需要插值到的气压层，气压的单位需要和pressure变量中的单位一致</td>
</tr>
<tr>
<td style="text-align: left;">with_nan</font></td>
<td style="text-align: left;">如果某个等压面层的某个区域或全部低于所有的模式面的气压或高于所有模式面的气压时，是否以nan的方式返回。如果with_nan = True，则以nan填充，如果with_nan = False,则以最底层或最顶层的数值填充</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;"><a href="https://www.showdoc.cc/meteva?page_id=3975600815874861">网格数据</a>，在level维度坐标和 level_list一致，在member，time,dtime,lat,lon维度和pressure完全一致。</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<p>下面是一个模式面数据的ctl文件的内容，下面以它为例对模式面到等压面算法的插值函数的用法进行说明
dset ^./gradsmodelvar_20230720120000.dat
options sequential big_endian
title MCV output
undef -99999.9
xdef 720 linear 0.000000000000000E+000 0.500000000000000<br />
ydef 361 linear -90.0000000000000 0.500000000000000<br />
zdef 60 levels
12.0611851562500<br />
37.6261804687500<br />
67.4102105034722<br />
104.011816927083<br />
149.855041406250<br />
207.193758940972<br />
278.116011197917<br />
364.548339843750<br />
468.260119878472<br />
590.867892968750<br />
733.839700781250<br />
898.499418315972<br />
1086.03108723958<br />
1297.48324921875<br />
1533.77327925347<br />
1795.69171901042<br />
2083.90661015625<br />
2398.96782769097<br />
2741.31141328125<br />
3111.26390859375<br />
3509.04668862847<br />
3934.78029505208<br />
4388.48876953125<br />
4870.10398706597<br />
5379.46998932292<br />
5916.34731796875<br />
6480.41734800347<br />
7071.28662109375<br />
7688.49117890625<br />
8331.50089644097<br />
8999.72381536458<br />
9692.51047734375<br />
10409.1582573785<br />
11148.9156971354<br />
11910.9868382813<br />
12694.5355558160<br />
13498.6898914063<br />
14322.5463867188<br />
15165.1744167535<br />
16025.6205231771<br />
16902.9127476563<br />
17796.0649651910<br />
18704.0812174479<br />
19625.9600460937<br />
20560.6988261285<br />
21507.2980992188<br />
22464.7659070313<br />
23432.1221245660<br />
24408.4027934896<br />
25392.6644554688<br />
26383.9884855035<br />
27381.4854252604<br />
28384.2993164063<br />
29391.6120339410<br />
30402.6476195313<br />
31416.6766148438<br />
32433.0203948785<br />
33451.0555013021<br />
34470.2179757813<br />
35490.0076933160<br />
tdef 73 linear 12Z20Jul2023 60mn
vars 41
rho 60 0 density
u 60 0 u wind
v 60 0 v wind
w 60 0 w wind
wzet 60 0 wzet wind
t 60 0 temperature
p 60 0 pressure
zs 0 0 terrain
tsfc 0 0 surface temperature
hprime 0 0 topographic sdv in m
semis 0 0 surface emissivity
sfcdlw 0 0 total sky surface downward lw flux
sfcdsw 0 0 total sky surface downward sw flux
sfcnsw 0 0 total sky surface net sw flux
snoalb 0 0 snow albedo
sfalb 0 0 surface albedo
tcldcov 0 0 total cloud cover
hcldcov 0 0 high cloud cover
mcldcov 0 0 middle cloud cover
lcldcov 0 0 low cloud cover
u10m 0 0 u 10m
v10m 0 0 v 10m
t2m 0 0 T 2m
q2m 0 0 q 2m
stsoil1 0 0 soil temp at level 1
stsoil2 0 0 soil temp at level 2
stsoil3 0 0 soil temp at level 3
stsoil4 0 0 soil temp at level 4
smsoil1 0 0 soil moisture at level 1
smsoil2 0 0 soil moisture at level 2
smsoil3 0 0 soil moisture at level 3
smsoil4 0 0 soil moisture at level 4
q1 60 0 moisture variable 1
q2 60 0 moisture variable 2
q3 60 0 moisture variable 3
q4 60 0 moisture variable 4
q5 60 0 moisture variable 5
q6 60 0 moisture variable 6
q7 60 0 moisture variable 7
q8 60 0 moisture variable 8
q9 60 0 moisture variable 9
endvars</p>
<pre><code class="language-python">ctl_path = r&quot;D:\data\ssc\mcv_grapesphys\liu1\model_20230720120000.ctl&quot;
data_path = r&quot;D:\data\ssc\mcv_grapesphys\liu1\gradsmodelvar_20230720120000.dat&quot;
grid0 = meb.grid([70,140,0.25],[15,55,0.25])</code></pre>
<pre><code class="language-python">pressure = meb.read_griddata_from_ctl(ctl_path,data_path, value_name =&quot;p&quot;, add_block_head_tail =True, dtime_dim = &quot;tdef&quot;,
endian=&quot;&gt;&quot;,dtime = 0,grid = grid0)
print(pressure)</code></pre>
<pre><code>&lt;xarray.DataArray 'data0' (member: 1, level: 60, time: 1, dtime: 1, lat: 161, lon: 281)&gt;
array([[[[[[100369.5 , 100372.03515625, 100374.5703125 , ...,
100994.359375 , 101002.70703125, 101011.0546875 ],
[100350.875 , 100352.02929688, 100353.18359375, ...,
100994.125 , 101004.90820312, 101015.69140625],
[100332.25 , 100332.0234375 , 100331.796875 , ...,
100993.890625 , 101007.109375 , 101020.328125 ],
...,
[ 98543.2890625 , 98565.60546875, 98587.921875 , ...,
100410.515625 , 100308.50390625, 100206.4921875 ],
[ 98516.2890625 , 98527.078125 , 98537.8671875 , ...,
100459.359375 , 100417.23046875, 100375.1015625 ],
[ 98489.2890625 , 98488.55078125, 98487.8125 , ...,
100508.203125 , 100525.95703125, 100543.7109375 ]]]],
[[[[100081.3671875 , 100083.90625 , 100086.4453125 , ...,
100704.0234375 , 100712.71484375, 100721.40625 ],
[100062.80078125, 100063.93359375, 100065.06640625, ...,
100704.16015625, 100715.16210938, 100726.1640625 ],
...
[ 684.64434814, 684.6782074 , 684.71206665, ...,
687.08734131, 687.1736145 , 687.2598877 ],
[ 685.14941406, 685.19223022, 685.23504639, ...,
687.44866943, 687.54232788, 687.63598633]]]],
[[[[ 523.00280762, 522.33575439, 521.66870117, ...,
522.12145996, 522.04290771, 521.96435547],
[ 523.44897461, 522.85507202, 522.26116943, ...,
522.61947632, 522.32125854, 522.02304077],
[ 523.8951416 , 523.37438965, 522.8536377 , ...,
523.11749268, 522.59960938, 522.08172607],
...,
[ 592.82427979, 592.80987549, 592.79547119, ...,
595.36773682, 595.57461548, 595.78149414],
[ 593.27096558, 593.26535034, 593.25973511, ...,
595.62716675, 595.84564209, 596.06411743],
[ 593.71765137, 593.7208252 , 593.72399902, ...,
595.88659668, 596.1166687 , 596.34674072]]]]]])
Coordinates:
* member (member) int32 0
* level (level) float64 12.06 37.63 67.41 ... 3.345e+04 3.447e+04 3.549e+04
* time (time) datetime64[ns] 2023-07-20T12:00:00
* dtime (dtime) int32 0
* lat (lat) float64 15.0 15.25 15.5 15.75 16.0 ... 54.25 54.5 54.75 55.0
* lon (lon) float64 70.0 70.25 70.5 70.75 ... 139.2 139.5 139.8 140.0
Attributes:
dtime_units: hour</code></pre>
<pre><code class="language-python">temp = meb.read_griddata_from_ctl(ctl_path,data_path, value_name =&quot;t&quot;, add_block_head_tail =True, dtime_dim = &quot;tdef&quot;, endian=&quot;&gt;&quot;,dtime = 0,grid = grid0)
print(temp)</code></pre>
<pre><code>&lt;xarray.DataArray 'data0' (member: 1, level: 60, time: 1, dtime: 1, lat: 161, lon: 281)&gt;
array([[[[[[300.78582764, 300.74124146, 300.69665527, ...,
300.28912354, 300.65727234, 301.02542114],
[300.74005127, 300.70143127, 300.66281128, ...,
300.7311554 , 300.93395233, 301.13674927],
[300.6942749 , 300.66162109, 300.62896729, ...,
301.17318726, 301.21063232, 301.24807739],
...,
[297.47692871, 297.40687561, 297.33682251, ...,
287.26333618, 286.89437866, 286.52542114],
[297.60624695, 297.62017059, 297.63409424, ...,
286.30445862, 286.06224823, 285.82003784],
[297.73556519, 297.83346558, 297.93136597, ...,
285.34558105, 285.2301178 , 285.11465454]]]],
[[[[300.52371216, 300.49238586, 300.46105957, ...,
300.04776001, 300.42070007, 300.79364014],
[300.47706604, 300.45231628, 300.42756653, ...,
300.47973633, 300.69236755, 300.90499878],
...
[242.32012177, 242.12757874, 241.93503571, ...,
243.59209442, 244.05978775, 244.52748108],
[242.41151428, 242.21092987, 242.01034546, ...,
243.66577148, 244.07575989, 244.48574829]]]],
[[[[239.11909485, 240.07873535, 241.03837585, ...,
232.61509705, 233.34049988, 234.06590271],
[238.73391724, 239.77864838, 240.82337952, ...,
233.19911194, 233.8210144 , 234.44291687],
[238.34873962, 239.4785614 , 240.60838318, ...,
233.78312683, 234.30152893, 234.81993103],
...,
[244.53910828, 244.45166016, 244.36421204, ...,
245.2043457 , 245.52687836, 245.84941101],
[244.58402252, 244.4955864 , 244.40715027, ...,
245.1865387 , 245.49009705, 245.7936554 ],
[244.62893677, 244.53951263, 244.4500885 , ...,
245.16873169, 245.45331573, 245.73789978]]]]]])
Coordinates:
* member (member) int32 0
* level (level) float64 12.06 37.63 67.41 ... 3.345e+04 3.447e+04 3.549e+04
* time (time) datetime64[ns] 2023-07-20T12:00:00
* dtime (dtime) int32 0
* lat (lat) float64 15.0 15.25 15.5 15.75 16.0 ... 54.25 54.5 54.75 55.0
* lon (lon) float64 70.0 70.25 70.5 70.75 ... 139.2 139.5 139.8 140.0
Attributes:
dtime_units: hour</code></pre>
<pre><code class="language-python">temp_p = meb.interp_zlevel_to_plevel_linear(pressure,temp,level_list=[85000,50000],with_nan=True)
meb.contourf_xy(temp_p,subplot = &quot;level&quot;,ncol = 2)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=0d59eab35c305cecb484e6d1471d5771&amp;file=file.png" alt="" /></p>
<pre><code class="language-python"></code></pre>
<pre><code class="language-python"></code></pre>