MeteoInfoLab脚本示例：计算温度平流

MeteoInfo

需要温度和风场U/V分量格点数据，计算中主要用到cdiff函数（结果已用GrADS验证一致）。

脚本程序：

1. print 'Open data files...'f_air = addfile('D:/Temp/nc/air.2011.nc')
   
2. f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')
   
3. f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')
   
4. 
   
5. print 'Read data array...'
   
6. tidx = 173    # Jun 23, 2011
   
7. t = f_air.gettime(tidx)
   
8. lidx = 0    # 1000 hPa
   
9. air = f_air['air'][tidx,lidx,:,:]
   
10. uwnd = f_uwnd['uwnd'][tidx,lidx,:,:]
    
11. vwnd = f_vwnd['vwnd'][tidx,lidx,:,:]
    
12. lon = f_air['lon'][:]
    
13. lat = f_air['lat'][:]
    
14. lon, lat = meshgrid(lon, lat)
    
15. 
    
16. # Calculate
    
17. print 'Calculate...'
    
18. dtx = cdiff(air, 1)
    
19. dty = cdiff(air, 0)
    
20. dx = cdiff(lon, 1) * pi / 180
    
21. dy = cdiff(lat, 0) * pi / 180
    
22. tadv = -1*((uwnd*dtx)/(cos(lat*pi/180)*dx)+vwnd*dty/dy)/6.37e6
    
23. 
    
24. #Plot
    
25. print 'Plot...'
    
26. axesm()
    
27. geoshow('country', edgecolor='black')
    
28. layer = contourfm(tadv, cmap='grads_rainbow')
    
29. title('Temperature advection (' + t.strftime('%Y-%m-%d') + ')')
    
30. colorbar(layer)
    
31. xlim(0, 360)
    
32. ylim(-90, 90)
    
33. xticks(arange(0, 361, 30))
    
34. yticks(arange(-90, 91, 30))

print 'Open data files...'f_air = addfile('D:/Temp/nc/air.2011.nc')<br /> f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')<br /> f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')<br /> <br /> print 'Read data array...'<br /> tidx = 173    # Jun 23, 2011<br /> t = f_air.gettime(tidx)<br /> lidx = 0    # 1000 hPa<br /> air = f_air['air'][tidx,lidx,:,:]<br /> uwnd = f_uwnd['uwnd'][tidx,lidx,:,:]<br /> vwnd = f_vwnd['vwnd'][tidx,lidx,:,:]<br /> lon = f_air['lon'][:]<br /> lat = f_air['lat'][:]<br /> lon, lat = meshgrid(lon, lat)<br /> <br /> # Calculate<br /> print 'Calculate...'<br /> dtx = cdiff(air, 1)<br /> dty = cdiff(air, 0)<br /> dx = cdiff(lon, 1) * pi / 180<br /> dy = cdiff(lat, 0) * pi / 180<br /> tadv = -1*((uwnd*dtx)/(cos(lat*pi/180)*dx)+vwnd*dty/dy)/6.37e6<br /> <br /> #Plot<br /> print 'Plot...'<br /> axesm()<br /> geoshow('country', edgecolor='black')<br /> layer = contourfm(tadv, cmap='grads_rainbow')<br /> title('Temperature advection (' + t.strftime('%Y-%m-%d') + ')')<br /> colorbar(layer)<br /> xlim(0, 360)<br /> ylim(-90, 90)<br /> xticks(arange(0, 361, 30))<br /> yticks(arange(-90, 91, 30))

学而无术

很给力，很强大

river

使用grads计算一下，得到的结果如下

色标不一样，今天有些晚了，等过几天闲下来弄个一样的色标的

MeteoInfo

river 发表于 2015-7-14 21:17
使用grads计算一下，得到的结果如下
看着是一样的的，不过t=173的时候是2011年6月22日

谢谢。GrADS里t是从1开始的，MeteoInfo里是从0开始的，所以你用GrADS验证的话要用t=174

river

MeteoInfo 发表于 2015-7-14 21:21
谢谢。GrADS里t是从1开始的，MeteoInfo里是从0开始的，所以你用GrADS验证的话要用t=174

不客气。过几天再加上平滑和一样的色标再看看，总觉得有些不一样呢。

MeteoInfo

river 发表于 2015-7-14 21:34
不客气。过几天再加上平滑和一样的色标再看看，总觉得有些不一样呢。

嗯，看起来相当不一样呀，数据用的一样的吗？

river

MeteoInfo 发表于 2015-7-14 21:45
嗯，看起来相当不一样呀，数据用的一样的吗？

数据是一样的啊，2011年的日值数据。今天上班，过几天我再弄一下，更新一下图和脚本吧······

river

这次附上脚本：'reinit'
'set grid off'
'set grads off'
'sdfopen G:\ncep\2011\air.2011.nc'
'sdfopen G:\ncep\2011\uwnd.2011.nc'
'sdfopen G:\ncep\2011\vwnd.2011.nc'
'set t 174'
*'pi=3.1416'
'dtx = cdiff(air, 'x')'
'dty = cdiff(air, 'y')'
'dx = cdiff(lon, x) * 3.1416/ 180'
'dy = cdiff(lat, y) * 3.1416/ 180'
'tadv = -1*((uwnd.2*dtx)/(cos(lat*3.1416/180)*dx)+vwnd.3*dty/dy)/6.37e6'
'run G:\GrADS-NCL-colormaker-v1.1\output\11colors.gs'
'set gxout shaded'
'd tadv'
'run cbar_interp.gs 1 1 0'
'draw title Temporature advection(2011-06-23)'
'gxprint g:/adtv.png white'
;


最终出图如下：

MeteoInfo

river 发表于 2015-7-19 15:35
这次附上脚本：'reinit'
'set grid off'
'set grads off'

谢谢！找到问题所在了，已在1楼更新。

1、高度层是第一层（1000 hPa），所以 lidx = 0。（之前为：idx = 1）

2、更重要的是此数据集中Y轴是反向的，因此在取 lat 时应该让数组反向：lat = f_air['lat'][::-1]。（之前为：lat = f_air['lat'][:]）

river

MeteoInfo 发表于 2015-7-19 16:44
谢谢！找到问题所在了，已在1楼更新。

1、高度层是第一层（1000 hPa），所以 lidx = 0。（之前为：idx ...

支持王老师，赶快出更多更好地教程，我已近准备学习一下Python