MeteoInfoLab脚本示例：获取气团轨迹每个节点的气象数据

MeteoInfo · 发表于 2015-10-2 09:56:21

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由 MeteoInfo 于 2023-10-16 22:28 编辑

读取HYSPLIT输出的轨迹数据文件和相应时间的气象数据文件，生成轨迹图层，循环每条轨迹的节点，读出该节点的经度、纬度、气压、时间，通过对气象数据插值获得该节点的气象数据。

脚本程序：

# Set working directory
trajDir = 'D:/Temp/HYSPLIT'
meteoDir = 'D:/Temp/arl'
# Open trjactory data file
print 'Open trajectory data file ...'
trajfn = os.path.join(trajDir, 'traj_20090731')
print 'Trajectory file: ' + trajfn
trajf = addfile_hytraj(trajfn)
# Read coordinates
lons = trajf['lon'][:]
lats = trajf['lat'][:]
press = trajf['PRESSURE'][:]
heights = trajf['height'][:]
tt = trajf['time'][:]
ntraj, np = lons.shape
# Open meteorological data file
print 'Open meteorological data file...'
meteofn = os.path.join(meteoDir, 'gdas1.jul09.w5')
print 'Meteorological file: ' + meteofn
meteof = addfile(meteofn)
# Get meteorological data along trajectory
print 'Get meteorological data along trajectory...'
outfn = os.path.join(trajDir, 'pblh_traj-1.txt')
outf = open(outfn, 'w')
outf.write('TrajID,Lon,Lat,Time,Heigh,PBLH,UWND\n')
pbldata = meteof['PBLH'][:]
udata = meteof['UWND'][:]
idx = 0
for i in range(ntraj):
for j in range(np):
lon = lons[i,j]
lat = lats[i,j]
pres = press[i,j]
z = heights[i,j]
t = tt[i,j]
pbl = pbldata.interpn([t, lat, lon])
uwnd = udata.interpn([t, pres, lat, lon])
t = miutil.num2date(t)
print 'TrajID: %i; lon: %.2f; lat: %.2f; time: %s; height: %.2f; PBLH: %.2f; UWND: %.2f' \
% (i, lon, lat, t.strftime('%Y%m%d_%H:%M'), z, pbl, uwnd)
line = '%i,%.4f,%.4f,%s,%.2f,%.2f,%.2f' % (i,lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pbl,uwnd)
outf.write(line + '\n')
t = t + datetime.timedelta(hours=-1)
idx += 1
outf.close()
print 'Finish...'

# Set working directory 
trajDir = 'D:/Temp/HYSPLIT' 
meteoDir = 'D:/Temp/arl' 
 
# Open trjactory data file 
print 'Open trajectory data file ...' 
trajfn = os.path.join(trajDir, 'traj_20090731') 
print 'Trajectory file: ' + trajfn 
trajf = addfile_hytraj(trajfn) 
# Read coordinates 
lons = trajf['lon'][:] 
lats = trajf['lat'][:] 
press = trajf['PRESSURE'][:] 
heights = trajf['height'][:] 
tt = trajf['time'][:] 
ntraj, np = lons.shape 
 
# Open meteorological data file 
print 'Open meteorological data file...' 
meteofn = os.path.join(meteoDir, 'gdas1.jul09.w5') 
print 'Meteorological file: ' + meteofn 
meteof = addfile(meteofn) 
 
# Get meteorological data along trajectory 
print 'Get meteorological data along trajectory...' 
outfn = os.path.join(trajDir, 'pblh_traj-1.txt') 
outf = open(outfn, 'w') 
outf.write('TrajID,Lon,Lat,Time,Heigh,PBLH,UWND\n') 
pbldata = meteof['PBLH'][:] 
udata = meteof['UWND'][:] 
idx = 0 
for i in range(ntraj): 
    for j in range(np): 
        lon = lons[i,j] 
        lat = lats[i,j] 
        pres = press[i,j] 
        z = heights[i,j] 
        t = tt[i,j] 
        pbl = pbldata.interpn([t, lat, lon]) 
        uwnd = udata.interpn([t, pres, lat, lon]) 
        t = miutil.num2date(t) 
        print 'TrajID: %i; lon: %.2f; lat: %.2f; time: %s; height: %.2f; PBLH: %.2f; UWND: %.2f' \ 
            % (i, lon, lat, t.strftime('%Y%m%d_%H:%M'), z, pbl, uwnd) 
        line = '%i,%.4f,%.4f,%s,%.2f,%.2f,%.2f' % (i,lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pbl,uwnd) 
        outf.write(line + '\n') 
        t = t + datetime.timedelta(hours=-1) 
    idx += 1 
 
outf.close() 
print 'Finish...'

二爷名声在外 · 发表于 2015-10-2 10:07:16

王老师假期还奋斗在科研岗位上，真心为您点个赞！

MeteoInfo · 发表于 2015-10-2 10:21:30

二爷名声在外发表于 2015-10-2 10:07
王老师假期还奋斗在科研岗位上，真心为您点个赞！

只能说是命苦呀

二爷名声在外 · 发表于 2015-10-2 17:48:03

MeteoInfo 发表于 2015-10-2 10:21
只能说是命苦呀

哈哈，我这条加班狗陪着王老师。。汪汪汪。。。

mm_cat · 发表于 2015-10-7 16:40:25

这个对我正有用，谢谢王老师！

mm_cat · 发表于 2015-10-12 18:13:53

本帖最后由 mm_cat 于 2015-10-12 18:22 编辑

王老师，想向您请教个问题。我的轨迹时间从19日一直到27日，需要两个GDAS文件，代码已经修改好，但是发现21日21、22、23三个小时的数据为空，用print meteof查询了GDAS文件信息后发现它的数据在每一个时间段最后一天都只到21时（实际上21时是无数据的，空缺3个时次）。因为模拟出来的轨迹是逐小时的，不想其中有时间点缺值，想问下王老师有没有什么办法可以将它补上，谢谢了！

查询到的信息：
File Name: C:/UserFiles/mty/thesis/data/GDAS\gdas1.jun15.w4
File Start Time: 2015-06-22 00:00
File End Time: 2015-06-28 21:00

File Name: C:/UserFiles/mty/thesis/data/GDAS\gdas1.jun15.w3
File Start Time: 2015-06-15 00:00
File End Time: 2015-06-21 21:00

用22日01时到21日18时轨迹测试Idle中结果如下：
Open trajectory data file ...
Trajectory file: C:/hysplit4/working/matrix/divide27\2701.txt
Open meteorological data file...
Meteorological file: C:/UserFiles/mty/thesis/data/GDAS\gdas1.jun15.w4
Get meteorological data along trajectory...
lon: 75.33; lat: 25.01; time: 20150622_01:00; height: 5698.20; PBLH: 630.49; TEMP: 404.58; RELH: -272.06
lon: 75.33; lat: 25.15; time: 20150622_00:00; height: 5703.70; PBLH: 506.52; TEMP: 409.52; RELH: -266.07
lon: 75.33; lat: 25.28; time: 20150621_23:00; height: 5724.50; PBLH: -9999.00; TEMP: -9999.00; RELH: -9999.00
lon: 75.35; lat: 25.39; time: 20150621_22:00; height: 5749.90; PBLH: -9999.00; TEMP: -9999.00; RELH: -9999.00
lon: 75.39; lat: 25.48; time: 20150621_21:00; height: 5778.90; PBLH: -9999.00; TEMP: -9999.00; RELH: -9999.00
lon: 75.43; lat: 25.57; time: 20150621_20:00; height: 5809.20; PBLH: 128.51; TEMP: 311.28; RELH: 33.58
lon: 75.48; lat: 25.66; time: 20150621_19:00; height: 5836.20; PBLH: 131.33; TEMP: 311.74; RELH: 31.72
lon: 75.54; lat: 25.75; time: 20150621_18:00; height: 5855.90; PBLH: 130.10; TEMP: 312.19; RELH: 30.18
Finish...

MeteoInfo · 发表于 2015-10-12 20:04:55

mm_cat 发表于 2015-10-12 18:13
王老师，想向您请教个问题。我的轨迹时间从19日一直到27日，需要两个GDAS文件，代码已经修改好，但是发现21 ...

目前没想到什么好办法，笨办法是记录下来两个文件之间空缺的几个时次的点的位置（x, y, z），获取前一个文件最后一个时次和后一个文件第一个时次相应的插值数据，然后再对时间进行线性插值获取最终数据。

mm_cat · 发表于 2015-10-13 12:40:38

MeteoInfo 发表于 2015-10-12 20:04
目前没想到什么好办法，笨办法是记录下来两个文件之间空缺的几个时次的点的位置（x, y, z），获取前一个 ...

好的，谢谢王老师！

还发现了一个问题，读出来的数据感觉不对，您看我上面截的测试数据中温度达到400+k，相对湿度出现了负数。
我用您上面的脚本，只修改气象数据和轨迹数据的路径进行测试，气象数据用的和您一样的，为gdas1.jul09.w5，轨迹我按着您截出来的图自行修改了一下，就只用了截图上显示出来的7行，UWND因为和轨迹数据中最后一列的气压相关结果不同我可以理解，可是PBLH只和经纬度、时间、高度相关，这个和您截图上的结果也有一点不同。我感觉我自己提取数据出现的错误值可能和这个有关系，但具体什么原因不知道，您能帮我分析一下么。

修改后的模拟结果，PBLH列也不一样
Open trajectory data file ...
Trajectory file: C:\hysplit4\working\matrix\divide27\2703.txt
Open meteorological data file...
Meteorological file: C:\UserFiles\mty\thesis\data\GDAS\gdas1.jul09.w5
Get meteorological data along trajectory...
lon: 114.10; lat: 35.77; time: 20090729_12:00; height: 635.90; PBLH: 234.01; UWND: -5.16
lon: 114.16; lat: 35.58; time: 20090729_11:00; height: 639.70; PBLH: 665.71; UWND: -1.65
lon: 114.24; lat: 35.42; time: 20090729_10:00; height: 640.40; PBLH: 1142.35; UWND: 1.05
lon: 114.34; lat: 35.29; time: 20090729_09:00; height: 639.10; PBLH: 1635.59; UWND: 2.86
lon: 114.44; lat: 35.17; time: 20090729_08:00; height: 642.90; PBLH: 1690.60; UWND: 4.19
lon: 114.52; lat: 35.06; time: 20090729_07:00; height: 658.30; PBLH: 1708.06; UWND: 5.26
lon: 114.59; lat: 34.95; time: 20090729_06:00; height: 685.80; PBLH: 1699.95; UWND: 6.07
Finish...

我使用的轨迹数据（最后气压那一列我是改了用来测试是否对UWND有影响的，这个数据对PBLH这一项应该没有影响）：
   4    1
GDAS    9    7    8    0    0
GDAS    9    7    9    0    0
GDAS    9    7 22    0    0
GDAS    9    7 29    0    0
   1 BACKWARD OMEGA
   9    7 29 12 35.770  114.100 500.0
   1 PRESSURE
   1    1 09    7 29 12    0    0    0.0 35.770  114.100 635.9 900.0
   1    1 09    7 29 11    0    1 -1.0 35.580  114.160 639.7 900.0
   1    1 09    7 29 10    0    2 -2.0 35.420  114.240 640.4 900.0
   1    1 09    7 29    9    0    3 -3.0 35.290  114.340 639.1 900.0
   1    1 09    7 29    8    0    2 -4.0 35.170  114.440 642.9 900.0
   1    1 09    7 29    7    0    1 -5.0 35.060  114.520 658.3 900.0
   1    1 09    7 29    7    0    0 -6.0 34.950  114.590 685.8 900.0

MeteoInfo · 发表于 2015-10-13 23:30:52

mm_cat 发表于 2015-10-13 12:40
好的，谢谢王老师！

还发现了一个问题，读出来的数据感觉不对，您看我上面截的测试数据中温度达到400+ ...

你用我脚本里用的轨迹数据问试试。

traj_20090731 (13.74 KB, 下载次数: 46)

mm_cat · 发表于 2015-10-14 16:58:23

本帖最后由 mm_cat 于 2015-10-14 17:01 编辑

MeteoInfo 发表于 2015-10-13 23:30
你用我脚本里用的轨迹数据问试试。

王老师，我试过了。第一次倒数第二行的风速值和您的有所不同，所以我又运行了四遍，然后对5次结果进行了比较，没有哪两次是完全相同的，每两个文件都会出现几个不同的数据。不过区别都在UWND上，PBLH都是一样的。

		自动登录	找回密码
密码			立即注册

MeteoInfoLab脚本示例：获取气团轨迹每个节点的气象数据

登录后查看更多精彩内容~

相关帖子