MeteoInfoLab脚本示例：水汽通量散度计算

MeteoInfo

用ncep数据计算水汽通量散度的脚本。需要air, uwnd, vwnd和rhum变量。数据是4维数据，需要固定时间维和高度维。脚本中用到几个内置的函数：cdiff, hdivg和magnitude，和GrADS中同名函数的作用一样。计算出来的水汽通量散点单位应该是：kg/(s•cm2•hPa)。

脚本程序如下：

1. print 'Open data files...'
   
2. f_air = addfile('D:/Temp/nc/air.2011.nc')
   
3. f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')
   
4. f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')
   
5. f_rhum = addfile('D:/Temp/nc/rhum.2011.nc')
   
6. 
   
7. print 'Read data array...'
   
8. tidx = 173    # Jun 23, 2011
   
9. t = f_air.gettime(tidx)
   
10. lidx = 3    # 700 hPa
    
11. air = f_air['air'][tidx,lidx,:,:]
    
12. uwnd = f_uwnd['uwnd'][tidx,lidx,:,:]
    
13. vwnd = f_vwnd['vwnd'][tidx,lidx,:,:]
    
14. rhum = f_rhum['rhum'][tidx,lidx,:,:]
    
15. 
    
16. # Calculate
    
17. print 'Calculate...'
    
18. prs = 700
    
19. g = 9.8
    
20. es = 6.112*exp(17.67*(air-273.16)/(air-29.65))
    
21. qs = 0.62197*es/(prs-0.378*es)
    
22. q = qs*rhum/100
    
23. qhdivg = hdivg(q*uwnd/g,q*vwnd/g)
    
24. 
    
25. #Plot
    
26. print 'Plot...'
    
27. axesm()
    
28. mlayer = shaperead('D:/Temp/map/country1.shp')
    
29. geoshow(mlayer, edgecolor='black')
    
30. layer = contourfm(qhdivg, cmap='grads_rainbow')
    
31. title('Water Vapor Flux Divergency (' + t.strftime('%Y-%m-%d') + ')')
    
32. colorbar(layer)
    
33. xlim(0, 360)
    
34. ylim(-90, 90)

print 'Open data files...'<br /> f_air = addfile('D:/Temp/nc/air.2011.nc')<br /> f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')<br /> f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')<br /> f_rhum = addfile('D:/Temp/nc/rhum.2011.nc')<br /> <br /> print 'Read data array...'<br /> tidx = 173    # Jun 23, 2011<br /> t = f_air.gettime(tidx)<br /> lidx = 3    # 700 hPa<br /> air = f_air['air'][tidx,lidx,:,:]<br /> uwnd = f_uwnd['uwnd'][tidx,lidx,:,:]<br /> vwnd = f_vwnd['vwnd'][tidx,lidx,:,:]<br /> rhum = f_rhum['rhum'][tidx,lidx,:,:]<br /> <br /> # Calculate<br /> print 'Calculate...'<br /> prs = 700<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 /> qhdivg = hdivg(q*uwnd/g,q*vwnd/g)<br /> <br /> #Plot<br /> print 'Plot...'<br /> axesm()<br /> mlayer = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(mlayer, edgecolor='black')<br /> layer = contourfm(qhdivg, cmap='grads_rainbow')<br /> title('Water Vapor Flux Divergency (' + t.strftime('%Y-%m-%d') + ')')<br /> colorbar(layer)<br /> xlim(0, 360)<br /> ylim(-90, 90)

结果图：

江南剑侠

王老师您好，我是MeteoInfo的初学者，python也刚刚开始学，第一次使用了这个脚本，因为我用的是FNL 1°*1°的Grib2数据，所以对脚本改动了一下，但是不知道对不对，因为画出来的水汽通量散度场有很大面积的负值异常区域，但是几个正、负值中心位置又好像是正确的。还望您指教！

print 'Open data files...'
f=addfile('F:/datafiles/fnl_20160701_00_00.grib2')

print 'Read data array...'

air = f['Temperature_isobaric'][:,25,:,:]
uwnd = f['u-component_of_wind_isobaric'][:,25,:,:]
vwnd = f['v-component_of_wind_isobaric'][:,25,:,:]
rhum = f['Relative_humidity_isobaric'][:,25,:,:]

# Calculate
print 'Calculate...'
prs = 850
g = 9.8
es = 6.112*exp(17.67*(air-273.16)/(air-29.65))
qs = 0.62197*es/(prs-0.378*es)
q = qs*rhum/100
test = cdiff(q, True)
qhdivg = hdivg(q*uwnd/g,q*vwnd/g)
qv = rhum*es/100
uv = magnitude(uwnd, vwnd)
uvq = uv*qv/(9.8*1000)

#Plot
print 'Plot...'
axesm()
mlayer = shaperead('F:/MeteoInfo/map/country1.shp')
geoshow(mlayer, linecolor='black')
#layer = contourfm(qhdivg, 20)
layer = contourfm(qhdivg, cmap='GMT_red2green')
title('Water Vapor Flux Divergency')
colorbar(layer)
xlim(80, 150)
ylim(0, 50)

MeteoInfo

江南剑侠 发表于 2016-8-18 13:19
王老师您好，我是MeteoInfo的初学者，python也刚刚开始学，第一次使用了这个脚本，因为我用的是FNL 1°*1° ...

脚本看起来没什么问题。

江南剑侠

MeteoInfo 发表于 2016-8-18 15:08
脚本看起来没什么问题。

                               
登录/注册后可看大图

                               
登录/注册后可看大图

当天长江中下游出现了区域性暴雨，在该区域上空可以看到有明显的水汽通量负值区，但是其他大部分地区也是负值区，值还比较单一，有点不正常。
一开始以为是填色的问题，所以后来换了一种填色发现还是这样

火枪裸鞋

王老师您好！想请教您一个问题：前面加载地图层都没错，但是一使用contourf就报错，错误信息是java.lang.nullpointerexception？该如何解决呢？

MeteoInfo

火枪裸鞋 发表于 2017-9-19 05:43
王老师您好！想请教您一个问题：前面加载地图层都没错，但是一使用contourf就报错，错误信息是java.lang.nu ...

建议你先仔细看看这里：http://www.meteothinker.com/docs ... /data_tutorial.html

打打基础再来写脚本。

GS晨光

王老师，850hpa lidx 应该等于几啊2 还是4

MeteoInfo

GS晨光 发表于 2017-11-28 22:39
王老师，850hpa lidx 应该等于几啊2 还是4

获取变量或者数组后可有用dimvalue方法获取某个维的具体数值。该方法只有一个参数，即维的序号。

离岸风

王老师，你的代码中时间tidx=173表示2011.06.23，是不是2011.01.01是tidx=0
那为啥高度层次lidx=3表示700hpa，是lidx起始数是1不是0，也就是lidx=1是1000hpa的意思么

瑛雩

老师您好，我有一个问题想请教一下。
我在计算水汽通量散度时，前面的计算都有结果，但是在调用函数hdivg，也就是qhdivg = hdivg(q*uwnd/g,q*vwnd/g)之后返回值为none，在这之前我输出q,uwnd,vwnd均有值，请问这是什么原因？应该怎么解决？