MeteoInfoLab脚本示例：随高度分布的平均场

MeteoInfo

气象分析中经常会用到某个变量在某个矩形范围内沿经度或纬度方向平均场的垂直分布。这里有两个例子，第一个例子是“正矩形”，第二个例子是“斜矩形”（用到格点场的旋转，参考此例：http://bbs.06climate.com/forum.p ... 1110&extra=page%3D1）。

脚本1：

1. f = addfile('D:/Temp/nc/uwnd.2011.nc')
   
2. data = f['uwnd'][0,:,'10:40','70:130']
   
3. data = mean(data, axis=1)
   
4. lev1 = data.dimvalue(0)
   
5. #lev2 = 1000 - lev1
   
6. lev2 = meteo.pressure_to_height_std(lev1)
   
7. levels = []
   
8. for j in range(0, len(lev1)):
   
9.     levels.append('%i' % lev1[j])
   
10. data.setdimvalue(0, lev2)
    
11. #Plot
    
12. layer = contourf(data, 20)
    
13. colorbar(layer)
    
14. yticks(lev2, levels)
    
15. xlabel('Longitude')
    
16. ylabel('Pressure (hPa)')
    
17. title('Latitude average uwnd profile')

f = addfile('D:/Temp/nc/uwnd.2011.nc')<br /> data = f['uwnd'][0,:,'10:40','70:130']<br /> data = mean(data, axis=1)<br /> lev1 = data.dimvalue(0)<br /> #lev2 = 1000 - lev1<br /> lev2 = meteo.pressure_to_height_std(lev1)<br /> levels = []<br /> for j in range(0, len(lev1)):<br />     levels.append('%i' % lev1[j])<br /> data.setdimvalue(0, lev2)<br /> #Plot<br /> layer = contourf(data, 20)<br /> colorbar(layer)<br /> yticks(lev2, levels)<br /> xlabel('Longitude')<br /> ylabel('Pressure (hPa)')<br /> title('Latitude average uwnd profile')

脚本2：

1. #Read data array
   
2. f = addfile('D:/Temp/nc/uwnd.2011.nc')
   
3. data = f['uwnd'][0,:,'0:70','50:150']
   
4. 
   
5. #Get rotated grid
   
6. lat = arange(10, 41, 5)
   
7. lon = arange(70, 131, 5)
   
8. xn = len(lon)
   
9. yn = len(lat)
   
10. lon, lat = meshgrid(lon, lat)
    
11. angle = 20 * pi / 180
    
12. rotate = array([[cos(angle),sin(angle)], [-sin(angle),cos(angle)]])
    
13. n = lon.size
    
14. xy = zeros([n, 2])
    
15. xy[:,0] = lon - 70
    
16. xy[:,1] = lat - 10
    
17. xy = dot(xy, rotate)    #Matrix multiplication
    
18. xy[:,0] = xy[:,0] + 70
    
19. xy[:,1] = xy[:,1] + 10
    
20. lon[:,:] = xy[:,0]
    
21. lat[:,:] = xy[:,1]
    
22. 
    
23. #Project data to the grid
    
24. ndata = data.project(lon, lat)
    
25. 
    
26. #Average data along tilt latitude
    
27. londata = mean(ndata, axis=1)
    
28. latdata = mean(ndata, axis=2)
    
29. lev1 = data.dimvalue(0)
    
30. lev2 = meteo.pressure_to_height_std(lev1)
    
31. levels = []
    
32. for j in range(0, len(lev1)):
    
33.     levels.append('%i' % lev1[j])
    
34. 
    
35. #Plot
    
36. subplot(2,1,1)
    
37. layer = contourf(lon[0,:], lev2, londata, 20)
    
38. colorbar(layer)
    
39. yticks(lev2, levels)
    
40. xlabel('Longitude')
    
41. ylabel('Pressure (hPa)')
    
42. title('Tilt latitude average uwnd profile')
    
43. 
    
44. subplot(2,1,2)
    
45. layer = contourf(lat[0,:], lev2, latdata, 20)
    
46. colorbar(layer)
    
47. yticks(lev2, levels)
    
48. xlabel('Latitude')
    
49. ylabel('Pressure (hPa)')
    
50. title('Tilt longitude average uwnd profile')

#Read data array<br /> f = addfile('D:/Temp/nc/uwnd.2011.nc')<br /> data = f['uwnd'][0,:,'0:70','50:150']<br /> <br /> #Get rotated grid<br /> lat = arange(10, 41, 5)<br /> lon = arange(70, 131, 5)<br /> xn = len(lon)<br /> yn = len(lat)<br /> lon, lat = meshgrid(lon, lat)<br /> angle = 20 * pi / 180<br /> rotate = array([[cos(angle),sin(angle)], [-sin(angle),cos(angle)]])<br /> n = lon.size<br /> xy = zeros([n, 2])<br /> xy[:,0] = lon - 70<br /> xy[:,1] = lat - 10<br /> xy = dot(xy, rotate)    #Matrix multiplication<br /> xy[:,0] = xy[:,0] + 70<br /> xy[:,1] = xy[:,1] + 10<br /> lon[:,:] = xy[:,0]<br /> lat[:,:] = xy[:,1]<br /> <br /> #Project data to the grid<br /> ndata = data.project(lon, lat)<br /> <br /> #Average data along tilt latitude<br /> londata = mean(ndata, axis=1)<br /> latdata = mean(ndata, axis=2)<br /> lev1 = data.dimvalue(0)<br /> lev2 = meteo.pressure_to_height_std(lev1)<br /> levels = []<br /> for j in range(0, len(lev1)):<br />     levels.append('%i' % lev1[j])<br /> <br /> #Plot<br /> subplot(2,1,1)<br /> layer = contourf(lon[0,:], lev2, londata, 20)<br /> colorbar(layer)<br /> yticks(lev2, levels)<br /> xlabel('Longitude')<br /> ylabel('Pressure (hPa)')<br /> title('Tilt latitude average uwnd profile')<br /> <br /> subplot(2,1,2)<br /> layer = contourf(lat[0,:], lev2, latdata, 20)<br /> colorbar(layer)<br /> yticks(lev2, levels)<br /> xlabel('Latitude')<br /> ylabel('Pressure (hPa)')<br /> title('Tilt longitude average uwnd profile')

xyan88

多谢楼主共享

风格而才

王老师您好，麻烦我想问一下有没有关于时间剖面（nc数据）的例子，类似于micaps中ec数据的时间剖面图。

Finnick

王老师您好，我初学MI,我想问一下王老师，就是这种剖面图，我想左侧Y轴显示特定的层次高度，王老师上面的图70 50 10hPa我想在Y轴不显示出来，把100hPa显示出来，我该这么设置啊

独孤酒见

版上的牛人们啊，请教问题啊。
#Get rotated grid
lat = arange(10, 41, 5)
lon = arange(70, 131, 5)
xn = len(lon)
yn = len(lat)
lon, lat = meshgrid(lon, lat)
angle = 20 * pi / 180
请教两个问题。
王老师给的脚本里面的这几行中，我把arange(70, 131, 5)中的参数改为2，即arange(70, 131, 2)。能出pressure_longitude图，后面的pressure_latitude图出不来，报错，说是下标“溢出”。

第二个问题，angle = 20 * pi / 180。这句话里面的这个20是啥意思啊。

多谢老师们。