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

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

脚本程序:

#-----------------------------------------------------
# Author: Yaqiang Wang
# Date: 2015-9-30
# Purpose: Get meteorological data along trajectory
# Note: Sample
#-----------------------------------------------------

# 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)
# Create trajectory layer
trajLayer = trajf.trajlayer()

# 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.txt\')
outf = open(outfn, \'w\')
outf.write(\'Lon,Lat,Time,Heigh,PBLH,UWND\n\')
pblvar = \'PBLH\'
uvar = \'UWND\'
idx = 0
for tline in trajLayer.shapes():
    t = trajLayer.cellvalue(\'StartDate\', idx)
    h = trajLayer.cellvalue(\'StartHour\', idx)    
    t.replace(hour=h)
    for ps in tline.getPoints():
        lon = ps.X
        lat = ps.Y 
        z = ps.M
        pres = ps.Z
        pbl = meteof.tostation(pblvar, lon, lat, None, t)
        uwnd = meteof.tostation(uvar, lon, lat, pres, t)
        print \'lon: %.2f; lat: %.2f; time: %s; height: %.2f; PBLH: %.2f; UWND: %.2f\' % (lon, lat, t.strftime(\'%Y%m%d_%H:%M\'), z, pbl, uwnd)
        line = \'%.4f,%.4f,%s,%.2f,%.2f,%.2f\' % (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...\'

版权声明:本文为yaqiang原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://www.cnblogs.com/yaqiang/p/4852094.html