温度地形高度订正
<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.product as mpd
import numpy as np
import datetime
import copy
import matplotlib.pyplot as plt
import pandas as pd</code></pre>
<pre><code>The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload</code></pre>
<pre><code class="language-python">sta_all = pd.read_hdf(r"H:\test_data\input\mpd\Example_data\sta_all.h5")
sta_all = meb.sele_by_para(sta_all,value = [-100,100]) #选取合理取值范围的数据，即进行简单的阈值质控
print(sta_all)</code></pre>
<pre><code> level time dtime id lon lat OBS ECMWF \
0 0 2019-07-01 08:00:00 0 54398 116.6 40.1 25.8 24.724
1 0 2019-07-01 08:00:00 0 54410 116.1 40.6 18.9 20.284
2 0 2019-07-01 08:00:00 0 54416 116.9 40.4 25.1 22.864
3 0 2019-07-01 08:00:00 0 54419 116.6 40.4 27.5 21.796
4 0 2019-07-01 08:00:00 0 54499 116.2 40.2 27.5 23.076
... ... ... ... ... ... ... ... ...
9187 0 2019-07-31 20:00:00 72 54410 116.1 40.6 18.8 21.556
9188 0 2019-07-31 20:00:00 72 54416 116.9 40.4 26.7 24.960
9189 0 2019-07-31 20:00:00 72 54419 116.6 40.4 27.0 22.920
9190 0 2019-07-31 20:00:00 72 54499 116.2 40.2 29.3 23.056
9191 0 2019-07-31 20:00:00 72 54412 116.6 40.7 25.4 21.640
GRAPES
0 23.844
1 21.740
2 22.704
3 23.136
4 23.736
... ...
9187 19.968
9188 23.404
9189 23.056
9190 22.568
9191 21.232
[9192 rows x 9 columns]</code></pre>
<pre><code class="language-python">mpd.score_id(sta_all,mem.me,subplot = "member",ncol = 2)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/831d9303a9a367873dc2dff8eed46337" alt="" /></p>
<pre><code>( level time dtime id lon lat ECMWF GRAPES
0 0 2019-07-01 08:00:00 0 54398 116.6 40.1 -1.085090 -2.624936
1 0 2019-07-01 08:00:00 0 54410 116.1 40.6 2.368527 1.228972
2 0 2019-07-01 08:00:00 0 54412 116.6 40.7 -1.357502 -1.741227
3 0 2019-07-01 08:00:00 0 54416 116.9 40.4 -0.726175 -2.728313
4 0 2019-07-01 08:00:00 0 54419 116.6 40.4 -2.398982 -2.982574
5 0 2019-07-01 08:00:00 0 54499 116.2 40.2 -2.769563 -3.780637,
[None])</code></pre>
<h1>温度地形订正</h1>
<p><strong><font face="黑体" color=blue size = 3>terrain_height_correct(sta_temp, grid, sta_alt = None, member_list = None,rate = 0.6)</font></strong><br />
在基于站点观测对网格温度进行检验时，有时站点在山顶，周围格点在山脚，或者站点在山谷，周围网格在腰或山顶，此时采用双线性插值将网格温度插值到站点，会存在一个系统性的偏差。在温度检验时理应扣除该偏差，本函数功能是完成该偏差的自动扣除。
该偏差的幅度除了和高度有关，也和温度垂直递减率有关，然而目前尚没有合理的方式动态的估算温度垂直递减率，因此本算法仅仅考虑温度垂直递减率为常数的情况。
订正的算法：</p>
<ol>
<li>本程序库提供了一套0.008333°分辨率的地形高度数据文件，程序第一次启动时会自动从网站下载该地形文件，如果网络不通时，则可按运行提示下载相关文件保存至相关路径</li>
<li>程序读取地形数据，获得网格温度插值前的网格点所在位置的地形高度。</li>
<li>将站点附近4个网格点的地形高度插值到站点，并和站点实际高度进行比对，得到高度差</li>
<li>将高度差乘以递减率和原始的温度插值结果相加。</li>
</ol>
<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_temp</font></strong></td>
<td style="text-align: left;">站点形式的温度数据，它是网格温度插值到地面站点之后的结果，它可以包含多个层次，时间和时效，因此它可以是多个插值结果的拼接结果</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>grid</font></strong></td>
<td style="text-align: left;">格点的经纬度信息，它插值前的网格温度场所对应的网格</td>
</tr>
<tr>
<td style="text-align: left;"><strong>sta_alt</strong></td>
<td style="text-align: left;">站点形式的站点高度数据，它包含 level,time,dtime,id,lon,lat,alt共7列数据，如果该参数缺省值是，则自动读取系统中的国家站高度信息</td>
</tr>
<tr>
<td style="text-align: left;"><strong>member_list</strong></td>
<td style="text-align: left;">缺省是对所有数据列进行订正，否正对指定对某些列进行订正</td>
</tr>
<tr>
<td style="text-align: left;"><strong>rate</strong></td>
<td style="text-align: left;">垂直递减率，单位 ℃/100m</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;">订正后的温度数据, 不在sta_alt站点列表中的数据不会被保留</td>
</tr>
</tbody>
</table>
<p><strong>调用示例</strong></p>
<pre><code class="language-python">#先订正ec,假设ec的网格分辨率是0.25, 注意该算法属于mpd模块
grid_ec = meb.grid([70,140,0.25],[15,55,0.25])
sta_all_thc1 = mpd.terrain_height_correct(sta_all,grid = grid_ec, member_list = ["ECMWF"])
#再订正grapes,假设grapes的网格分辨率是0.05
grid_grapes = meb.grid([70,140,0.05],[15,55,0.05])
sta_all_thc2 = mpd.terrain_height_correct(sta_all_thc1,grid = grid_grapes, member_list = ["GRAPES"]) </code></pre>
<pre><code class="language-python">mpd.score_id(sta_all_thc2,mem.me,subplot = "member",ncol = 2)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitfile/sign/6ec8aebc96a77fa1978dd73760b3b39d" alt="" /></p>
<pre><code>( level time dtime id lon lat ECMWF GRAPES
0 0 2019-07-01 08:00:00 0 54398 116.6 40.1 -1.005936 -2.545783
1 0 2019-07-01 08:00:00 0 54410 116.1 40.6 -0.057341 -1.196895
2 0 2019-07-01 08:00:00 0 54412 116.6 40.7 0.291411 -0.092314
3 0 2019-07-01 08:00:00 0 54416 116.9 40.4 0.300203 -1.701935
4 0 2019-07-01 08:00:00 0 54419 116.6 40.4 0.019818 -0.563775
5 0 2019-07-01 08:00:00 0 54499 116.2 40.2 -2.374523 -3.385597,
[None])</code></pre>
<p>对比上述结果，注意左上角的54110站，该站位于一座山峰之上，周围网格点海拔更低，因此从模式直接插值到站点的温度值出现明显的偏高线性，在做地形高度订正之后ME转变为略低于0.</p>