MeteoInfoLab脚本示例：垂直螺旋度

MeteoInfo

参考此贴（http://bbs.06climate.com/forum.php?mod=viewthread&tid=2396）做出的计算垂直螺旋度（单层）的脚本。正确与否未经验证。

脚本程序：

1. print 'Open data files...'
   
2. f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')
   
3. f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')
   
4. f_omega = addfile('D:/Temp/nc/omega.2011.nc')
   
5. 
   
6. print 'Calculate vertical helicity...'
   
7. tidx = 173    # Jun 23, 2011
   
8. t = f_uwnd.gettime(tidx)
   
9. level = '700'
   
10. lat = '15:55'
    
11. lon = '70:135'
    
12. uwnd = f_uwnd['uwnd'][tidx,level,lat,lon][::-1,:]
    
13. vwnd = f_vwnd['vwnd'][tidx,level,lat,lon][::-1,:]
    
14. omega = f_omega['omega'][tidx,level,lat,lon][::-1,:]
    
15. wd = hcurl(uwnd, vwnd)
    
16. lx = -(wd*omega*10.)/12.64*1e6
    
17. 
    
18. print 'Plot...'
    
19. axesm()
    
20. mlayer = shaperead('D:/Temp/map/country1.shp')
    
21. geoshow(mlayer, edgecolor='black')
    
22. layer = contourfm(lx, 20)
    
23. title('Vertical helicity (' + t.strftime('%Y-%m-%d') + ')')
    
24. colorbar(layer)

print 'Open data files...'<br /> f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')<br /> f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')<br /> f_omega = addfile('D:/Temp/nc/omega.2011.nc')<br /> <br /> print 'Calculate vertical helicity...'<br /> tidx = 173    # Jun 23, 2011<br /> t = f_uwnd.gettime(tidx)<br /> level = '700'<br /> lat = '15:55'<br /> lon = '70:135'<br /> uwnd = f_uwnd['uwnd'][tidx,level,lat,lon][::-1,:]<br /> vwnd = f_vwnd['vwnd'][tidx,level,lat,lon][::-1,:]<br /> omega = f_omega['omega'][tidx,level,lat,lon][::-1,:]<br /> wd = hcurl(uwnd, vwnd)<br /> lx = -(wd*omega*10.)/12.64*1e6<br /> <br /> print 'Plot...'<br /> axesm()<br /> mlayer = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(mlayer, edgecolor='black')<br /> layer = contourfm(lx, 20)<br /> title('Vertical helicity (' + t.strftime('%Y-%m-%d') + ')')<br /> colorbar(layer)

rceclx

MeteoInfo 发表于 2015-10-9 10:47
更新了MeteoInfoLab使其能进行多次螺旋度等的运算，下面的例子演示了计算多层垂直螺旋度并绘制垂直分布图。 ...

有点“高大上”的感觉

MeteoInfo

更新了MeteoInfoLab使其能进行多次螺旋度等的运算，下面的例子演示了计算多层垂直螺旋度并绘制垂直分布图。lx1是北纬40度的剖面数据，y轴为气压，在追踪等值线前先用p2h函数转为高度，或者用其它的方式（等距或者等气压，需要确保y轴的值是递增的）。

脚本程序：

1. print 'Open data files...'
   
2. f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')
   
3. f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')
   
4. f_omega = addfile('D:/Temp/nc/omega.2011.nc')
   
5. 
   
6. print 'Calculate vertical helicity...'
   
7. tidx = 173    # Jun 23, 2011
   
8. t = f_uwnd.gettime(tidx)
   
9. level = '1000:100'
   
10. lat = '15:55'
    
11. lon = '70:135'
    
12. uwnd = f_uwnd['uwnd'][tidx,level,lat,lon][:,::-1,:]
    
13. vwnd = f_vwnd['vwnd'][tidx,level,lat,lon][:,::-1,:]
    
14. omega = f_omega['omega'][tidx,level,lat,lon][:,::-1,:]
    
15. wd = hcurl(uwnd, vwnd)
    
16. lx = -(wd*omega*10.)/12.64*1e6
    
17. lx1 = lx[:,'40',:]
    
18. lev1 = lx1.dimvalue(0)
    
19. #lev2 = 1000 - lev1
    
20. lev2 = meteo.p2h(lev1)
    
21. levels = []
    
22. for j in range(0, len(lev1)):
    
23.     levels.append('%i' % lev1[j])
    
24. lx1.setdimvalue(0, lev2)
    
25. 
    
26. print 'Plot...'
    
27. layer = contourf(lx1, 20)
    
28. title('Vertical helicity (' + t.strftime('%Y-%m-%d') + ')')
    
29. yticks(lev2, levels)
    
30. xlabel('Longitude')
    
31. ylabel('Pressure (hPa)')
    
32. colorbar(layer)

print 'Open data files...'<br /> f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')<br /> f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')<br /> f_omega = addfile('D:/Temp/nc/omega.2011.nc')<br /> <br /> print 'Calculate vertical helicity...'<br /> tidx = 173    # Jun 23, 2011<br /> t = f_uwnd.gettime(tidx)<br /> level = '1000:100'<br /> lat = '15:55'<br /> lon = '70:135'<br /> uwnd = f_uwnd['uwnd'][tidx,level,lat,lon][:,::-1,:]<br /> vwnd = f_vwnd['vwnd'][tidx,level,lat,lon][:,::-1,:]<br /> omega = f_omega['omega'][tidx,level,lat,lon][:,::-1,:]<br /> wd = hcurl(uwnd, vwnd)<br /> lx = -(wd*omega*10.)/12.64*1e6<br /> lx1 = lx[:,'40',:]<br /> lev1 = lx1.dimvalue(0)<br /> #lev2 = 1000 - lev1<br /> lev2 = meteo.p2h(lev1)<br /> levels = []<br /> for j in range(0, len(lev1)):<br />     levels.append('%i' % lev1[j])<br /> lx1.setdimvalue(0, lev2)<br /> <br /> print 'Plot...'<br /> layer = contourf(lx1, 20)<br /> title('Vertical helicity (' + t.strftime('%Y-%m-%d') + ')')<br /> yticks(lev2, levels)<br /> xlabel('Longitude')<br /> ylabel('Pressure (hPa)')<br /> colorbar(layer)

有事秦天

王老师，您好，如果我需要使用fnl资料绘制一个500-200hPa沿东经90°的涡度垂直剖面该怎么画呢，尝试了一下午都没能成功，特向王老师求教，附上数据 。
fnl_20140816_12_00.grib2 (13.94 MB, 下载次数: 11)

2017-4-23 19:08 上传

点击文件名下载附件

半颗麦芽糖

王老师，看了上一楼的问题，我尝试参照您上面给的脚本修改画图，但是不知道为什么报错了，过程中有一次成功出图，但是修改了一些参数以后又报错了。脚本和报错如下，请王老师指点，谢谢！

1. 
   
2. <P> print 'Open data files...'
   
3. fn = 'F:/tuer/fnl20140816/M08/fnl_20140816_12_00.grib2'
   
4. f = addfile(fn)
   
5. print 'Calculate vertical helicity...'
   
6. tidx = 0    # Jun 23, 2011
   
7. level = [50000, 20000]
   
8. lat = [25,40]
   
9. lon = [75,105]
   
10. u = f['u-component_of_wind_isobaric'][tidx,level,lat,lon][:,::-1,:]
    
11. v = f['v-component_of_wind_isobaric'][tidx,level,lat,lon][:,::-1,:]
    
12. vort = hcurl(u, v)
    
13. wd1 = vort[:,:,[90]]
    
14. lev1 = wd1.dimvalue(0)/100
    
15. #lev2 = 1000 - lev1
    
16. lev2 = p2h(lev1)
    
17. levels = []
    
18. for i in range(0, len(lev1)):
    
19.     levels.append('%i' % lev1[i])
    
20. wd1.setdimvalue(0, lev2)
    
21. print 'Plot...'
    
22. layer = contourf(wd1,20)
    
23. #title('Vertical helicity (' + t.strftime('%Y-%m-%d') + ')')
    
24. yticks(lev2, levels)
    
25. xlabel('Latitude')
    
26. ylabel('Pressure (hPa)')
    
27. colorbar(layer)
    
28. </P>

复制代码

报错：Traceback (most recent call last):
  File "<iostream>", line 21, in <module>
  File "F:\11\MeteoInfo\pylib\mipylib\miplot.py", line 3489, in contourf
    igraphic = GraphicFactory.createContourPolygons(gdata.data, ls, smooth)
at wContour.Contour.traceIsoline_UndefData(Contour.java:2067)
at wContour.Contour.isoline_UndefData(Contour.java:2405)
at wContour.Contour.createContourLines_UndefData(Contour.java:551)
at wContour.Contour.tracingContourLines(Contour.java:38)
at org.meteoinfo.drawing.ContourDraw.tracingContourLines(ContourDraw.java:81)
at org.meteoinfo.chart.plot.GraphicFactory.createContourPolygons(GraphicFactory.java:560)
at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
at sun.reflect.NativeMethodAccessorImpl.invoke(Unknown Source)
at sun.reflect.DelegatingMethodAccessorImpl.invoke(Unknown Source)
at java.lang.reflect.Method.invoke(Unknown Source)
java.lang.ArrayIndexOutOfBoundsException: java.lang.ArrayIndexOutOfBoundsException: -1

MeteoInfo

有事秦天 发表于 2017-4-23 19:09
王老师，您好，如果我需要使用fnl资料绘制一个500-200hPa沿东经90°的涡度垂直剖面该怎么画呢，尝试了一下 ...

参考下面的脚本：

1. fn = 'D:/Temp/grib/fnl_20140816_12_00.grib2'
   
2. f = addfile(fn)
   
3. print 'Calculate vertical helicity...'
   
4. tidx = 0
   
5. level = '50000:20000'
   
6. lat = '25:40'
   
7. lon = '75:105'
   
8. u = f['u-component_of_wind_isobaric'][tidx,level,lat,lon][::-1,:,:]
   
9. v = f['v-component_of_wind_isobaric'][tidx,level,lat,lon][::-1,:,:]
   
10. vort = hcurl(u, v)
    
11. wd1 = vort[:,:,'90']
    
12. lev1 = wd1.dimvalue(0)/100
    
13. lev1 = lev1[::-1]
    
14. lev2 = meteo.p2h(lev1)
    
15. levels = []
    
16. for j in range(0, len(lev1)):
    
17.     levels.append('%i' % lev1[j])
    
18. wd1.setdimvalue(0, lev2)
    
19. print 'Plot...'
    
20. layer = contourf(wd1,20)
    
21. title(u'Vertical helicity (Longitude: 90°E)')
    
22. yticks(lev2, levels)
    
23. xlabel('Latitude')
    
24. ylabel('Pressure (hPa)')
    
25. colorbar(layer)

fn = 'D:/Temp/grib/fnl_20140816_12_00.grib2'<br /> f = addfile(fn)<br /> print 'Calculate vertical helicity...'<br /> tidx = 0<br /> level = '50000:20000'<br /> lat = '25:40'<br /> lon = '75:105'<br /> u = f['u-component_of_wind_isobaric'][tidx,level,lat,lon][::-1,:,:]<br /> v = f['v-component_of_wind_isobaric'][tidx,level,lat,lon][::-1,:,:]<br /> vort = hcurl(u, v)<br /> wd1 = vort[:,:,'90']<br /> lev1 = wd1.dimvalue(0)/100<br /> lev1 = lev1[::-1]<br /> lev2 = meteo.p2h(lev1)<br /> levels = []<br /> for j in range(0, len(lev1)):<br />     levels.append('%i' % lev1[j])<br /> wd1.setdimvalue(0, lev2)<br /> print 'Plot...'<br /> layer = contourf(wd1,20)<br /> title(u'Vertical helicity (Longitude: 90°E)')<br /> yticks(lev2, levels)<br /> xlabel('Latitude')<br /> ylabel('Pressure (hPa)')<br /> colorbar(layer)

MeteoInfo

半颗麦芽糖 发表于 2017-4-23 22:28
王老师，看了上一楼的问题，我尝试参照您上面给的脚本修改画图，但是不知道为什么报错了，过程中有一次成功 ...

较新的MeteoInfoLab对于维数据值用字符串来表示，列表只能表示序号。

半颗麦芽糖

MeteoInfo 发表于 2017-4-24 10:03
较新的MeteoInfoLab对于维数据值用字符串来表示，列表只能表示序号。

感谢王老师回复，查看了一下自己用的还是1.3.9版本的，更新了以后照着您提供的脚本出图成功了，看来确实是这些时间关注得少了，还得继续加油！再次感谢王老师！

半颗麦芽糖

MeteoInfo 发表于 2017-4-24 10:03
较新的MeteoInfoLab对于维数据值用字符串来表示，列表只能表示序号。

王老师，您好，我将这个fnl资料放到超算上，使用milab.sh画相同的图，结果报错，报错如下，是不是milab版本的问题呢？

MeteoInfo

半颗麦芽糖 发表于 2017-4-24 14:36
王老师，您好，我将这个fnl资料放到超算上，使用milab.sh画相同的图，结果报错，报错如下，是不是milab版 ...

很有可能是版本问题