MeteoInfoLab脚本示例：整层水汽通量散度

↘禹禹獨行…

MeteoInfo 发表于 2016-6-24 09:30
找到出错的原因了，我之前下载的rhum数据只有8层，新的数据都是17层。所以用新的数据需要修改相应的代码 ...

谢谢楼主，可以正常出图的。最近再学习使用这个软件，楼主太赞了！！！

江南剑侠

王老师，MeteoInfoLab中的很多函数或语法看不懂，像这个程序中的
levels = v_rhum.dimvalue(1)   中的“dimvalue”，
zn = uwnd.dimlen(0)  中的dimlen，
20行q = zeros((zn,yn,xn))，21行又q = DimArray(q, uwnd.dims, uwnd.fill_value, uwnd.proj)，DimArray、 uwnd.dims、uwnd.fill_value、 uwnd.proj又是什么意思？
这些函数哪里能查到呢？MeteoInfo的官网和别的地方都找不到

MeteoInfo

江南剑侠 发表于 2016-8-21 10:28
王老师，MeteoInfoLab中的很多函数或语法看不懂，像这个程序中的
levels = v_rhum.dimvalue(1)   中的“di ...

目前MeteoInfoLab的文档还很不完善，不过所有的函数都可以在源代码中看到，MeteoInfo -> pylib -> mipylib目录中。注意不要随意修改这些代码。

江南剑侠

MeteoInfo 发表于 2016-8-21 12:16
目前MeteoInfoLab的文档还很不完善，不过所有的函数都可以在源代码中看到，MeteoInfo -> pylib -> mipyli ...

1. print 'Open data files...'
   
2. f=addfile('F:/datafiles/fnl_20160702_00_00.grib2')
   
3. 
   
4. 
   
5. print 'Read data array...'
   
6. 
   
7. f_air = f['Temperature_isobaric'][:,:,:,:]
   
8. f_uwnd = f['u-component_of_wind_isobaric'][:,:,:,:]
   
9. f_vwnd = f['v-component_of_wind_isobaric'][:,:,:,:]
   
10. f_rhum = f['Relative_humidity_isobaric'][:,:,:,:]
    
11. v_rhum = f_rhum['Relative_humidity_isobaric']
    
12. 
    
13. levels = v_rhum.dimvalue(1)
    
14. 
    
15. t = f_air.gettime(tidx)
    
16. uwnd = f_uwnd['uwnd'][:,:,:,:]
    
17. vwnd = f_vwnd['vwnd'][:,:,:,:]
    
18. 
    
19. # Calculate
    
20. print 'Calculate...'
    
21. zn = uwnd.dimlen(0)
    
22. yn = uwnd.dimlen(1)
    
23. xn = uwnd.dimlen(2)
    
24. q = zeros((zn,yn,xn))
    
25. q = DimArray(q, uwnd.dims, uwnd.fill_value, uwnd.proj)
    
26. zidx = 0
    
27. for prs in uwnd.dimvalue(0):
    
28. air = f_air['Temperature_isobaric'][:,zidx,:,:]
    
29. rhum = f_rhum['Relative_humidity_isobaric'][:,zidx,:,:]
    
30. g = 9.8
    
31. es = 6.112*exp(17.67*(air-273.16)/(air-29.65))
    
32. qs = 0.62197*es/(prs-0.378*es)
    
33. qq = qs*rhum/100
    
34. q[zidx,:,:] = qq
    
35. zidx += 1
    
36. qu = q*uwnd/g
    
37. qv = q*vwnd/g
    
38. qhdivg = hdivg(qu, qv)
    
39. vqhdivg = trapz(qhdivg, levels)
    
40. 
    
41. #Plot
    
42. print 'Plot...'
    
43. axesm()
    
44. mlayer = shaperead(’F:/MeteoInfo/map/country1.shp‘)
    
45. geoshow(mlayer, edgecolor='black')
    
46. layer = contourfm(vqhdivg, 20)
    
47. title('Vertical Integrated Water Vapor Flux Divergency ')
    
48. colorbar(layer)
    
49. xlim(0, 360)
    
50. ylim(-90, 90)

复制代码

我使用的是NCEP的FNL的1°*1°GRIB2每天4次数据，但是运行总是出错，报错的地方不太明白，还请王老师指教
>>> run script...
Open data files...
Read data array...
indices must be 3 dimensions!
Traceback (most recent call last):
  File "<iostream>", line 13, in <module>
AttributeError: 'NoneType' object has no attribute 'dimvalue'

MeteoInfo

江南剑侠 发表于 2016-8-22 14:11
我使用的是NCEP的FNL的1°*1°GRIB2每天4次数据，但是运行总是出错，报错的地方不太明白，还请王老师指 ...

v_rhum = f_rhum['Relative_humidity_isobaric']
改为：
v_rhum = f['Relative_humidity_isobaric']

江南剑侠

MeteoInfo 发表于 2016-8-22 14:14
v_rhum = f_rhum['Relative_humidity_isobaric']
改为：
v_rhum = f['Relative_humidity_isobaric']

明白了
现在可以出图，但是和我之前画的水汽通量散度一样，有很多负值的异常大值区，我也在找原因呢

                               
登录/注册后可看大图

MeteoInfo

江南剑侠 发表于 2016-8-22 14:28
明白了
现在可以出图，但是和我之前画的水汽通量散度一样，有很多负值的异常大值区，我也在找原因呢

估计是MeteoInfo软件有bug，我抽空看看。

江南剑侠

MeteoInfo 发表于 2016-8-22 15:03
估计是MeteoInfo软件有bug，我抽空看看。

麻烦了，会不会因为FNL资料和NC不同？FNL的rhum垂直方向有31层，和这个有关么？

MeteoInfo

江南剑侠 发表于 2016-8-22 15:09
麻烦了，会不会因为FNL资料和NC不同？FNL的rhum垂直方向有31层，和这个有关么？

应该是等值线填色算法的bug，将contourfm改为imshowm，

1. print 'Open data files...'
   
2. f = addfile('D:/Temp/grib/fnl_20160702_00_00.grib2')
   
3. 
   
4. print 'Read data array...'
   
5. v_rhum = f['Relative_humidity_isobaric']
   
6. levels = v_rhum.dimvalue(1)
   
7. tidx = 0    # July 2, 2016
   
8. t = f.gettime(tidx)
   
9. uwnd = f['u-component_of_wind_isobaric'][tidx,:,:,:]
   
10. vwnd = f['v-component_of_wind_isobaric'][tidx,:,:,:]
    
11. 
    
12. # Calculate
    
13. print 'Calculate...'
    
14. zn = uwnd.dimlen(0)
    
15. yn = uwnd.dimlen(1)
    
16. xn = uwnd.dimlen(2)
    
17. q = zeros((zn,yn,xn))
    
18. q = DimArray(q, uwnd.dims, uwnd.fill_value, uwnd.proj)
    
19. zidx = 0
    
20. for prs in uwnd.dimvalue(0):
    
21.     prs = prs * 0.001
    
22.     air = f['Temperature_isobaric'][tidx,zidx,::-1,:]
    
23.     rhum = f['Relative_humidity_isobaric'][tidx,zidx,::-1,:]
    
24.     g = 9.8
    
25.     es = 6.112*exp(17.67*(air-273.16)/(air-29.65))
    
26.     qs = 0.62197*es/(prs-0.378*es)
    
27.     qq = qs*rhum/100
    
28.     q[zidx,:,:] = qq
    
29.     zidx += 1
    
30. qu = q*uwnd/g
    
31. qv = q*vwnd/g
    
32. qhdivg = hdivg(qu, qv)   
    
33. vqhdivg = trapz(qhdivg, levels)
    
34. 
    
35. #Plot
    
36. print 'Plot...'
    
37. axesm()
    
38. mlayer = shaperead('D:/Temp/map/country1.shp')
    
39. geoshow(mlayer, edgecolor='black')
    
40. layer = imshowm(vqhdivg, 20)
    
41. title('Vertical Integrated Water Vapor Flux Divergency (' + t.strftime('%Y-%m-%d') + ')')
    
42. colorbar(layer)
    
43. xlim(0, 360)
    
44. ylim(-90, 90)

print 'Open data files...'<br /> f = addfile('D:/Temp/grib/fnl_20160702_00_00.grib2')<br /> <br /> print 'Read data array...'<br /> v_rhum = f['Relative_humidity_isobaric']<br /> levels = v_rhum.dimvalue(1)<br /> tidx = 0    # July 2, 2016<br /> t = f.gettime(tidx)<br /> uwnd = f['u-component_of_wind_isobaric'][tidx,:,:,:]<br /> vwnd = f['v-component_of_wind_isobaric'][tidx,:,:,:]<br /> <br /> # Calculate<br /> print 'Calculate...'<br /> zn = uwnd.dimlen(0)<br /> yn = uwnd.dimlen(1)<br /> xn = uwnd.dimlen(2)<br /> q = zeros((zn,yn,xn))<br /> q = DimArray(q, uwnd.dims, uwnd.fill_value, uwnd.proj)<br /> zidx = 0<br /> for prs in uwnd.dimvalue(0):<br />     prs = prs * 0.001<br />     air = f['Temperature_isobaric'][tidx,zidx,::-1,:]<br />     rhum = f['Relative_humidity_isobaric'][tidx,zidx,::-1,:]<br />     g = 9.8<br />     es = 6.112*exp(17.67*(air-273.16)/(air-29.65))<br />     qs = 0.62197*es/(prs-0.378*es)<br />     qq = qs*rhum/100<br />     q[zidx,:,:] = qq<br />     zidx += 1<br /> qu = q*uwnd/g<br /> qv = q*vwnd/g<br /> qhdivg = hdivg(qu, qv)    <br /> vqhdivg = trapz(qhdivg, levels)<br /> <br /> #Plot<br /> print 'Plot...'<br /> axesm()<br /> mlayer = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(mlayer, edgecolor='black')<br /> layer = imshowm(vqhdivg, 20)<br /> title('Vertical Integrated Water Vapor Flux Divergency (' + t.strftime('%Y-%m-%d') + ')')<br /> colorbar(layer)<br /> xlim(0, 360)<br /> ylim(-90, 90)

江南剑侠

MeteoInfo 发表于 2016-8-22 15:35
应该是等值线填色算法的bug，将contourfm改为imshowm，

                               
登录/注册后可看大图

可以了！看来确实是这个原因呢