Get meteorological data along trajectoryΒΆ

Read trajectory endpoint data file and create trajectory polyline layer. Then read endpoint from the layer and interpolate the meteorological data to the endpoint location using interpn function of the data array.

# 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-1.txt')
outf = open(outfn, 'w')
pbldata = meteof['PBLH'][None]
udata = meteof['UWND'][None]
idx = 0
for tline in trajLayer.shapes():
    t = trajLayer.cellvalue('Date', idx)
    h = trajLayer.cellvalue('Hour', idx)
    for ps in tline.getPoints():
        lon = ps.X
        lat = ps.Y
        pres = ps.M
        z = ps.Z
        pbl = pbldata.interpn([t, lat, lon])
        uwnd = udata.interpn([t, pres, lat, lon])
        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

print 'Finish...'