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

MeteoInfo

增加了trapz函数用来对多维数据的某一维进行积分，该函数的文档见此网页：http://www.meteothink.org/docs/m ... numeric.minum.trapz。利用此函数可以对多维气象数据进行垂直方向的积分，比如计算整层水汽通量散度。

示例脚本如下。结果未经验证，欢迎讨论。

1. print 'Open data files...'
   
2. datadir = 'D:/Temp/nc'
   
3. f_air = addfile(datadir + '/air.2011.nc')
   
4. f_uwnd = addfile(datadir + '/uwnd.2011.nc')
   
5. f_vwnd = addfile(datadir + '/vwnd.2011.nc')
   
6. f_rhum = addfile(datadir + '/rhum.2011.nc')
   
7. 
   
8. print 'Read data array...'
   
9. v_rhum = f_rhum['rhum']
   
10. levels = v_rhum.dimvalue(1)
    
11. zn = len(levels)
    
12. tidx = 173    # Jun 23, 2011
    
13. t = f_air.gettime(tidx)
    
14. air = f_air['air'][tidx,:zn,::-1,:]
    
15. uwnd = f_uwnd['uwnd'][tidx,:zn,::-1,:]
    
16. vwnd = f_vwnd['vwnd'][tidx,:zn,::-1,:]
    
17. rhum = f_rhum['rhum'][tidx,:zn,::-1,:]
    
18. 
    
19. # Calculate
    
20. print 'Calculate...'
    
21. zn,yn,xn = uwnd.shape
    
22. prs = zeros((zn,yn,xn))
    
23. prs = dim_array(prs, air.dims)
    
24. zidx = 0
    
25. for pr in levels:
    
26.     prs[zidx,:,:] = pr
    
27.     zidx += 1
    
28. g = 9.8
    
29. es = 6.112*exp(17.67*(air-273.16)/(air-29.65))
    
30. qs = 0.62197*es/(prs-0.378*es)
    
31. q = qs*rhum/100  
    
32. qu = q*uwnd/g
    
33. qv = q*vwnd/g
    
34. qhdivg = hdivg(qu, qv)   
    
35. vqhdivg = trapz(qhdivg, levels)
    
36. 
    
37. #Plot
    
38. print 'Plot...'
    
39. axesm()
    
40. geoshow('country', edgecolor='black')
    
41. layer = contourfm(vqhdivg, 20)
    
42. title('Vertical Integrated Water Vapor Flux Divergency (' + t.strftime('%Y-%m-%d') + ')')
    
43. colorbar(layer)
    
44. xlim(0, 360)
    
45. ylim(-90, 90)

print 'Open data files...'<br /> datadir = 'D:/Temp/nc'<br /> f_air = addfile(datadir + '/air.2011.nc')<br /> f_uwnd = addfile(datadir + '/uwnd.2011.nc')<br /> f_vwnd = addfile(datadir + '/vwnd.2011.nc')<br /> f_rhum = addfile(datadir + '/rhum.2011.nc')<br /> <br /> print 'Read data array...'<br /> v_rhum = f_rhum['rhum']<br /> levels = v_rhum.dimvalue(1)<br /> zn = len(levels)<br /> tidx = 173    # Jun 23, 2011<br /> t = f_air.gettime(tidx)<br /> air = f_air['air'][tidx,:zn,::-1,:]<br /> uwnd = f_uwnd['uwnd'][tidx,:zn,::-1,:]<br /> vwnd = f_vwnd['vwnd'][tidx,:zn,::-1,:]<br /> rhum = f_rhum['rhum'][tidx,:zn,::-1,:]<br /> <br /> # Calculate<br /> print 'Calculate...'<br /> zn,yn,xn = uwnd.shape<br /> prs = zeros((zn,yn,xn))<br /> prs = dim_array(prs, air.dims)<br /> zidx = 0<br /> for pr in levels:<br />     prs[zidx,:,:] = pr<br />     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 /> q = qs*rhum/100  <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 /> geoshow('country', edgecolor='black')<br /> layer = contourfm(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)

kmdqlb

王老师再造一个“MATLAB+GrADS”

黑色青春

王老师给力

↘禹禹獨行…

楼主你好，运行你的脚本提示出错，这是什么原因？
>>> run script...
Open data files...
Read data array...
Calculate...
Traceback (most recent call last):
  File "<iostream>", line 31, in <module>
  File "D:\MeteoInfo\pylib\mipylib\dimarray.py", line 188, in __sub__
    r = DimArray(self.array.__sub__(other.array), self.dims, self.fill_value, self.proj)
  File "D:\MeteoInfo\pylib\mipylib\miarray.py", line 148, in __sub__
    r = MIArray(ArrayMath.sub(self.array, other.array))
java.lang.ArrayIndexOutOfBoundsException

java.lang.ArrayIndexOutOfBoundsException: java.lang.ArrayIndexOutOfBoundsException

MeteoInfo

↘禹禹獨行… 发表于 2016-6-23 11:27
楼主你好，运行你的脚本提示出错，这是什么原因？
>>> run script...
Open data files...

1、首先确保你用的是最新版的MeteoInfo
2、脚本程序你是否有改动？改动的什么地方？
3、所用的数据文件是不是有变化？

↘禹禹獨行…

MeteoInfo 发表于 2016-6-23 11:34
1、首先确保你用的是最新版的MeteoInfo
2、脚本程序你是否有改动？改动的什么地方？
3、所用的数据文件 ...

是最新版的MeteoInfo

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

复制代码

脚本程序只是改动了对应路径，使用的是ncep的数据。

MeteoInfo

↘禹禹獨行… 发表于 2016-6-23 12:45
是最新版的MeteoInfo

脚本程序只是改动了对应路径，使用的是ncep的数据。

你用的日数据还是每日4次的数据？

↘禹禹獨行…

MeteoInfo 发表于 2016-6-23 14:56
你用的日数据还是每日4次的数据？

数据是逐日的

MeteoInfo

↘禹禹獨行… 发表于 2016-6-24 08:48
数据是逐日的

之前的脚本有bug，已经更新，你再试试更新后的脚本。

MeteoInfo

↘禹禹獨行… 发表于 2016-6-24 08:48
数据是逐日的

找到出错的原因了，我之前下载的rhum数据只有8层，新的数据都是17层。所以用新的数据需要修改相应的代码。