MeteoInfoLab脚本示例：计算湿位涡

MeteoInfo

利用温度、相对湿度、风场U/V分量数据计算湿位涡。示例脚本简要说明：

1、读取数组。脚本中使用的是ARL格式的4维格点数据，数据是全球1*1度的，为了节省内存，读取4维数组时限定了经纬度范围。数组读取后可以通过数组的dimlen()方法来获取某个维的长度，用dimvalue()方法获取某个维的值的1维数组。该数据集中相对湿度（RELH）变量的高度层只有21层，而温度（TEMP）和U/V（UWND/VWND）变量的高度层有23层，在做数组运算时要求公式中所有数组的各维长度都一致，因此在读取温度和U/V数组时高度维只读到第21层和相对湿度数组一致。

2、计算涡度。利用hcurl()函数和U/V数组计算涡度。

3、构造高度层4维数组。利用dimvalue()方法获取了高度层的1维数组lev，但之后的计算中其它数组都是4维，为了匹配利用zeros()函数构造一个高度层的4维数组prs（所有值均为0）。此时的prs数组是MIArray类型，并没有各维的值的信息，因此创建一个DimArray对象加入维的详细信息。然后对高度层循环给prs数组赋值。

4、计算假相当位温。论坛里有很多帖子讨论假相当位温的计算，这里就不详细说了，公式参考GrADS脚本改写。需要注意的是计算出来的假相当位温eqt/thse是4维的。

5、计算湿位涡。湿位涡的计算也参考了GrADS脚本，由于涉及到上下层的差值，只能一层层计算。先构造一个湿位涡的4维数组mpv，再对时间维和高度维进行循环，在循环中计算2维（Y/X）湿位涡，再赋值给mpv。循环结束后4维mpv数组计算完成。

6、绘图检查计算结果。从4维mpv数组中固定时间和高度维取出2维数组（乘以10的6次方）绘制平面等值线填色图，可以看出计算结果的数量级是对的，具体结果正确与否未进行严格验证。

脚本程序：

1. # Calculate moisture potential vorticity
   
2. # Set working directory
   
3. trajDir = 'D:/Temp/HYSPLIT'
   
4. meteoDir = r'U:\data\ARL\2015'
   
5. 
   
6. # Open meteorological data file
   
7. print 'Open meteorological data file...'
   
8. meteofn = os.path.join(meteoDir, 'gdas1.mar15.w5')
   
9. print 'Meteorological file: ' + meteofn
   
10. meteof = addfile(meteofn)
    
11. 
    
12. # Read data
    
13. print 'Read data...'
    
14. latlim = [10,60]
    
15. lonlim = [60,140]
    
16. rh = meteof['RELH'][:,:,latlim,lonlim]
    
17. nx = rh.dimlen(3)
    
18. ny = rh.dimlen(2)
    
19. nz = rh.dimlen(1)
    
20. nt = rh.dimlen(0)
    
21. lat = rh.dimvalue(2)
    
22. lev = rh.dimvalue(1)
    
23. t0 = meteof['TEMP'][:,:nz-1,latlim,lonlim]
    
24. uwnd = meteof['UWND'][:,:nz-1,latlim,lonlim]
    
25. vwnd = meteof['VWND'][:,:nz-1,latlim,lonlim]
    
26. vort = hcurl(uwnd, vwnd)
    
27. prs = zeros([nt,nz,ny,nx])
    
28. prs = dim_array(prs, rh.dims)
    
29. for j in range(nz):
    
30.     prs[:,j,:,:] = lev[j]
    
31. 
    
32. # Calculate pseudo-equivalent potential temperature
    
33. print 'Clalulate pseudo-equivalent potential temperature...'
    
34. es = 6.1078*exp(17.2693882*(t0-273.16)/(t0-35.86))
    
35. qq = rh*(0.62197*es/(prs-0.378*es))/100.
    
36. e = prs*qq/(0.62197+qq)+1e-10
    
37. tlcl = 55.0+2840.0/(3.5*log(t0)-log(e)-4.805)
    
38. theta = t0*pow((1000/prs),(0.2854*(1.0-0.28*qq)))
    
39. eqt = theta*exp(((3376./tlcl)-2.54)*qq*(1.0+0.81*qq))
    
40. thse = eqt-273.15
    
41. 
    
42. # Calculate moisture potential vorticity
    
43. print 'Calculate moisture potential vorticity...'
    
44. mpv = zeros([nt,nz,ny,nx], dtype='double')
    
45. mpv = dim_array(mpv, rh.dims)
    
46. mpv.setdimvalue(1, lev[1:nz-1])
    
47. for t in range(nt):
    
48.     tt = meteof.gettime(t)
    
49.     print tt.strftime('%Y-%m-%d %H:00')
    
50.     for z in range(1, nz-1):
    
51.         #print '\tLevel: %i' % z
    
52.         f = zeros([ny,nx])
    
53.         f1 = 2*7.292*sin(lat*3.14159/180.0)*0.00001
    
54.         for i in range(nx):
    
55.             f[:,i] = f1
    
56.         g = 9.8
    
57.         dp = 100*(lev[z-1]-lev[z+1])
    
58.         deqt = eqt[t,z-1,:,:]-eqt[t,z+1,:,:]
    
59.         du = uwnd[t,z-1,:,:]-uwnd[t,z+1,:,:]
    
60.         dv = vwnd[t,z-1,:,:]-vwnd[t,z+1,:,:]
    
61.         dx1 = 2.0*6370949.0*cos(lat*3.14159/180.0)*3.14159/180.0
    
62.         dx = zeros([ny,nx])
    
63.         for i in range(nx):
    
64.             dx[:,i] = dx1
    
65.         dy = 2.0*6370949.0*3.14159/180.0
    
66.         dtx = cdiff(eqt[t,z,:,:], 1)
    
67.         dty = cdiff(eqt[t,z,:,:], 0)
    
68.         pv1 = -g*(vort[t,z,:,:]+f)*deqt/dp  
    
69.         pv2 = g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))
    
70.         pv = pv1+pv2
    
71.         mpv[t,z-1,:,:] = pv
    
72. 
    
73. #Plot
    
74. axesm()
    
75. lworld = shaperead('D:/Temp/Map/country1.shp')
    
76. geoshow(lworld, edgecolor='k')
    
77. t = 0
    
78. tt = meteof.gettime(t)
    
79. z = 5
    
80. clevs = arange(-3,3.1,0.5)
    
81. layer = contourfm(mpv[t,z,:,:]*1e6, clevs)
    
82. colorbar(layer)
    
83. title('Moisture potential vorticity (%i hPa)\n' % lev[z] + \
    
84.     tt.strftime('%Y-%m-%d %H:00'))
    
85. 
    
86. print 'Finish...'

# Calculate moisture potential vorticity<br /> # Set working directory<br /> trajDir = 'D:/Temp/HYSPLIT'<br /> meteoDir = r'U:\data\ARL\2015'<br /> <br /> # Open meteorological data file<br /> print 'Open meteorological data file...'<br /> meteofn = os.path.join(meteoDir, 'gdas1.mar15.w5')<br /> print 'Meteorological file: ' + meteofn<br /> meteof = addfile(meteofn)<br /> <br /> # Read data<br /> print 'Read data...'<br /> latlim = [10,60]<br /> lonlim = [60,140]<br /> rh = meteof['RELH'][:,:,latlim,lonlim]<br /> nx = rh.dimlen(3)<br /> ny = rh.dimlen(2)<br /> nz = rh.dimlen(1)<br /> nt = rh.dimlen(0)<br /> lat = rh.dimvalue(2)<br /> lev = rh.dimvalue(1)<br /> t0 = meteof['TEMP'][:,:nz-1,latlim,lonlim]<br /> uwnd = meteof['UWND'][:,:nz-1,latlim,lonlim]<br /> vwnd = meteof['VWND'][:,:nz-1,latlim,lonlim]<br /> vort = hcurl(uwnd, vwnd)<br /> prs = zeros([nt,nz,ny,nx])<br /> prs = dim_array(prs, rh.dims)<br /> for j in range(nz):<br />     prs[:,j,:,:] = lev[j]<br /> <br /> # Calculate pseudo-equivalent potential temperature<br /> print 'Clalulate pseudo-equivalent potential temperature...'<br /> es = 6.1078*exp(17.2693882*(t0-273.16)/(t0-35.86))<br /> qq = rh*(0.62197*es/(prs-0.378*es))/100.<br /> e = prs*qq/(0.62197+qq)+1e-10<br /> tlcl = 55.0+2840.0/(3.5*log(t0)-log(e)-4.805)<br /> theta = t0*pow((1000/prs),(0.2854*(1.0-0.28*qq)))<br /> eqt = theta*exp(((3376./tlcl)-2.54)*qq*(1.0+0.81*qq))<br /> thse = eqt-273.15<br /> <br /> # Calculate moisture potential vorticity<br /> print 'Calculate moisture potential vorticity...'<br /> mpv = zeros([nt,nz,ny,nx], dtype='double')<br /> mpv = dim_array(mpv, rh.dims)<br /> mpv.setdimvalue(1, lev[1:nz-1])<br /> for t in range(nt):<br />     tt = meteof.gettime(t)<br />     print tt.strftime('%Y-%m-%d %H:00')<br />     for z in range(1, nz-1):<br />         #print '\tLevel: %i' % z<br />         f = zeros([ny,nx])<br />         f1 = 2*7.292*sin(lat*3.14159/180.0)*0.00001<br />         for i in range(nx):<br />             f[:,i] = f1<br />         g = 9.8<br />         dp = 100*(lev[z-1]-lev[z+1]) <br />         deqt = eqt[t,z-1,:,:]-eqt[t,z+1,:,:]<br />         du = uwnd[t,z-1,:,:]-uwnd[t,z+1,:,:]<br />         dv = vwnd[t,z-1,:,:]-vwnd[t,z+1,:,:]<br />         dx1 = 2.0*6370949.0*cos(lat*3.14159/180.0)*3.14159/180.0<br />         dx = zeros([ny,nx])<br />         for i in range(nx):<br />             dx[:,i] = dx1<br />         dy = 2.0*6370949.0*3.14159/180.0<br />         dtx = cdiff(eqt[t,z,:,:], 1)<br />         dty = cdiff(eqt[t,z,:,:], 0)<br />         pv1 = -g*(vort[t,z,:,:]+f)*deqt/dp  <br />         pv2 = g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))<br />         pv = pv1+pv2<br />         mpv[t,z-1,:,:] = pv<br /> <br /> #Plot<br /> axesm()<br /> lworld = shaperead('D:/Temp/Map/country1.shp')<br /> geoshow(lworld, edgecolor='k')<br /> t = 0<br /> tt = meteof.gettime(t)<br /> z = 5<br /> clevs = arange(-3,3.1,0.5)<br /> layer = contourfm(mpv[t,z,:,:]*1e6, clevs)<br /> colorbar(layer)<br /> title('Moisture potential vorticity (%i hPa)\n' % lev[z] + \<br />     tt.strftime('%Y-%m-%d %H:00'))<br /> <br /> print 'Finish...'

绘图结果：

良辰

王老师，如果后向轨迹的时效比较长，比如后向10天的，一个数据资料集就不能满足计算需求，该如何添加那？学生这样meteof = addfile(meteofn1,meteofn2,meteofn3,meteofn4)，结果rh = meteof['RELH'][:,:,latlim,lonlim]报错，错误信息为TypeError: 'NoneType' object is unsubscriptable。王老师如有时间，希望指点。

KIMO23

  感謝老師的分享

良辰

王老师，源代码（只修改了资料路径）没进行其他修改的时候运行，出现如下错误代码！已经下载更新了meteoinfolab软件，1.3.5版本！
Open meteorological data file...
Meteorological file: D:\Temp\arl\gdas1.mar15.w5
Read data...
Traceback (most recent call last):
  File "<iostream>", line 28, in <module>
  File "D:\Temp\MeteoInfo\pylib\mipylib\dimarray.py", line 170, in __setitem__
    r = ArrayMath.setSection(self.array.array, ranges, value)
        at ucar.ma2.ArrayDouble.getDouble(ArrayDouble.java:208)
        at ucar.ma2.ArrayDouble.getObject(ArrayDouble.java:232)
        at org.meteoinfo.data.ArrayMath.setSection(ArrayMath.java:1505)
        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: 21

良辰

错误行是：
     prs[:,i,:,:] = lev

MeteoInfo

良辰 发表于 2016-8-18 18:18
王老师，源代码（只修改了资料路径）没进行其他修改的时候运行，出现如下错误代码！已经下载更新了meteoinf ...

错误信息提示第28行出错，28行是什么语句？

MeteoInfo

MeteoInfo 发表于 2016-8-18 18:21
错误信息提示第28行出错，28行是什么语句？

看看lev数组的长度是多少（len(lev)）？应该和nz一样。

良辰

prs[:,i,:,:] = lev

良辰

MeteoInfo 发表于 2016-8-18 18:27
看看lev数组的长度是多少（len(lev)）？应该和nz一样。

查看了是一样的，都是21

MeteoInfo

良辰 发表于 2016-8-18 19:21
查看了是一样的，都是21

for j in range(nz):
    prs[:,j,:,:] = lev[j]

之前用的 i 做循环，论坛会自动把中括号 i 当成斜体格式处理，用 j 会避免此问题（论坛的bug）。

良辰

MeteoInfo 发表于 2016-8-18 21:51
for j in range(nz):
    prs[:,j,:,:] = lev[j]

谢谢王老师，这么不容易发现的bug都能找出来，万分感谢！我进一步测试，辛苦了老师！