ARL资料计算涡度输出NaN，为什么呀？

良辰

meteofn = os.path.join(meteoDir, 'gdas1.mar15.w5')
print 'Meteorological file: ' + meteofn
meteof = addfile(meteofn)        
u = meteof['UWND'][0,0,:,:]
v = meteof['VWND'][0,0,:,:]
vort = hcurl(u, v)
print vort

MeteoInfo

涡度计算需要四周都有数据，因此格点场的四边均为NaN

良辰

MeteoInfo 发表于 2016-8-4 09:07
涡度计算需要四周都有数据，因此格点场的四边均为NaN

恩，知道了王老师！那我想把我计算出来的涡度，Get meteorological data along trajectory。该如何设置那？例如温度是temp = meteof.tostation(tvar, lon, lat, pres, t)。

1. #-----------------------------------------------------
   
2. # Author: Yaqiang Wang
   
3. # Date: 2015-9-30
   
4. # Purpose: Get meteorological data along trajectory
   
5. # Note: Sample
   
6. #-----------------------------------------------------
   
7. # Set working directory
   
8. trajDir = 'D:/Temp/HYSPLIT'
   
9. meteoDir = 'D:/Temp/arl'
   
10. # Open trjactory data file
    
11. print 'Open trajectory data file ...'
    
12. trajfn = os.path.join(trajDir, 'traj_20150331')
    
13. print 'Trajectory file: ' + trajfn
    
14. trajf = addfile_hytraj(trajfn)
    
15. # Create trajectory layer
    
16. trajLayer = trajf.trajlayer()
    
17. # Open meteorological data file
    
18. print 'Open meteorological data file...'
    
19. meteofn = os.path.join(meteoDir, 'gdas1.mar15.w5')
    
20. print 'Meteorological file: ' + meteofn
    
21. meteof = addfile(meteofn)        
    
22. u = meteof['UWND'][0,0,:,:]
    
23. v = meteof['VWND'][0,0,:,:]
    
24. vort = hcurl(u, v)
    
25. #axesm()
    
26. #mlayer = shaperead('D:\Temp\MeteoInfo\map/country1.shp')
    
27. #geoshow(mlayer)
    
28. #layer = contourfm(vort, 20)
    
29. #layer = contourfm(divg, 20)
    
30. #title('Vorticity')
    
31. #title('Divergence')
    
32. #colorbar(layer, orientation='horizontal', aspect=50)
    
33. #print u;print v;print vort
    
34. # Get meteorological data along trajectory
    
35. print 'Get meteorological data along trajectory...'
    
36. outfn = os.path.join(trajDir, 'traj.txt')
    
37. outf = open(outfn, 'w')
    
38. outf.write('Lon,Lat,Time,Heigh,pres,PBLH,UWND,VWND,WWND,HGTS,TEMP,RELH,CAPE,vort\n')
    
39. pblvar = 'PBLH'
    
40. uvar = 'UWND'
    
41. vvar = 'VWND'
    
42. wvar = 'WWND'
    
43. hvar = 'HGTS'
    
44. tvar = 'TEMP'
    
45. rhvar = 'RELH'
    
46. capevar = 'CAPE'
    
47. vortvar = 'vort'
    
48. #divg = hdivg(uvar, vvar)
    
49. idx = 0
    
50. for tline in trajLayer.shapes():
    
51.     t = trajLayer.cellvalue('Date', idx)
    
52.     h = trajLayer.cellvalue('Hour', idx)   
    
53.     t.replace(hour=h)
    
54.     for ps in tline.getPoints():
    
55.         lon = ps.X
    
56.         lat = ps.Y
    
57.         z = ps.M
    
58.         pres = ps.Z
    
59.         pbl = meteof.tostation(pblvar, lon, lat, None, t)
    
60.         uwnd = meteof.tostation(uvar, lon, lat, pres, t)
    
61.         vwnd = meteof.tostation(vvar, lon, lat, pres, t)
    
62.         wwnd = meteof.tostation(wvar, lon, lat, pres, t)
    
63.         hgt = meteof.tostation(hvar, lon, lat, pres, t)
    
64.         temp = meteof.tostation(tvar, lon, lat, pres, t)
    
65.         rh = meteof.tostation(rhvar, lon, lat, pres, t)
    
66.         cape = meteof.tostation(capevar, lon, lat, None, t)
    
67. #     vort = meteof.tostation(vortvar, lon, lat, pres, t)   这个地方如何写？？？？？？？？？
    
68.        print 'lon: %.3f; lat: %.3f; time: %s; height: %.1f; pres：%.1f; PBLH: %.1f; UWND: %.1f;VWND: %.1f;WWND: %.5f;HGTS: %.1f;TEMP: %.1f;RELH: %.1f;CAPE: %.1f;vort:%.5f' % (lon, lat, t.strftime('%Y%m%d_%H:%M'), z, pres, pbl, uwnd, vwnd, wwnd, hgt, temp, rh, cape, vort)
    
69.         line = '%.3f,%.3f,%s,%.1f,%.1f,%.1f,%.1f,%.1f,%.5f,%.1f,%.1f,%.1f,%.1f,%.5f' % (lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pres,pbl,uwnd,vwnd,wwnd,hgt,temp,rh,cape,vort)</P>
    
70. <P>        outf.write(line + '\n')
    
71.         t = t + datetime.timedelta(hours=-1)
    
72.     idx += 1
    
73. outf.close()
    
74. print 'Finish...' </P>

复制代码