请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
楼主: MeteoInfo

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

[复制链接]

新浪微博达人勋

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

谢谢楼主,可以正常出图的。最近再学习使用这个软件,楼主太赞了!!!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-8-21 10:28:13 | 显示全部楼层
王老师,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的官网和别的地方都找不到

密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-8-21 12:16:26 | 显示全部楼层
江南剑侠 发表于 2016-8-21 10:28
王老师,MeteoInfoLab中的很多函数或语法看不懂,像这个程序中的
levels = v_rhum.dimvalue(1)   中的“di ...

目前MeteoInfoLab的文档还很不完善,不过所有的函数都可以在源代码中看到,MeteoInfo -> pylib -> mipylib目录中。注意不要随意修改这些代码。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-8-22 14:11:21 | 显示全部楼层
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. print 'Read data array...'

  4. f_air = f['Temperature_isobaric'][:,:,:,:]
  5. f_uwnd = f['u-component_of_wind_isobaric'][:,:,:,:]
  6. f_vwnd = f['v-component_of_wind_isobaric'][:,:,:,:]
  7. f_rhum = f['Relative_humidity_isobaric'][:,:,:,:]
  8. v_rhum = f_rhum['Relative_humidity_isobaric']

  9. levels = v_rhum.dimvalue(1)

  10. t = f_air.gettime(tidx)
  11. uwnd = f_uwnd['uwnd'][:,:,:,:]
  12. vwnd = f_vwnd['vwnd'][:,:,:,:]

  13. # Calculate
  14. print 'Calculate...'
  15. zn = uwnd.dimlen(0)
  16. yn = uwnd.dimlen(1)
  17. xn = uwnd.dimlen(2)
  18. q = zeros((zn,yn,xn))
  19. q = DimArray(q, uwnd.dims, uwnd.fill_value, uwnd.proj)
  20. zidx = 0
  21. for prs in uwnd.dimvalue(0):
  22. air = f_air['Temperature_isobaric'][:,zidx,:,:]
  23. rhum = f_rhum['Relative_humidity_isobaric'][:,zidx,:,:]
  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. #Plot
  35. print 'Plot...'
  36. axesm()
  37. mlayer = shaperead(’F:/MeteoInfo/map/country1.shp‘)
  38. geoshow(mlayer, edgecolor='black')
  39. layer = contourfm(vqhdivg, 20)
  40. title('Vertical Integrated Water Vapor Flux Divergency ')
  41. colorbar(layer)
  42. xlim(0, 360)
  43. 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'
密码修改失败请联系微信:mofangbao

新浪微博达人勋

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

v_rhum = f_rhum['Relative_humidity_isobaric']
改为:
v_rhum = f['Relative_humidity_isobaric']
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-8-22 14:28:05 | 显示全部楼层
MeteoInfo 发表于 2016-8-22 14:14
v_rhum = f_rhum['Relative_humidity_isobaric']
改为:
v_rhum = f['Relative_humidity_isobaric']

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

                               
登录/注册后可看大图
密码修改失败请联系微信:mofangbao

新浪微博达人勋

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

估计是MeteoInfo软件有bug,我抽空看看。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-8-22 15:09:51 | 显示全部楼层
MeteoInfo 发表于 2016-8-22 15:03
估计是MeteoInfo软件有bug,我抽空看看。

麻烦了,会不会因为FNL资料和NC不同?FNL的rhum垂直方向有31层,和这个有关么?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-8-22 15:35:19 | 显示全部楼层
江南剑侠 发表于 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. print 'Read data array...'
  4. v_rhum = f['Relative_humidity_isobaric']
  5. levels = v_rhum.dimvalue(1)
  6. tidx = 0    # July 2, 2016
  7. t = f.gettime(tidx)
  8. uwnd = f['u-component_of_wind_isobaric'][tidx,:,:,:]
  9. vwnd = f['v-component_of_wind_isobaric'][tidx,:,:,:]

  10. # Calculate
  11. print 'Calculate...'
  12. zn = uwnd.dimlen(0)
  13. yn = uwnd.dimlen(1)
  14. xn = uwnd.dimlen(2)
  15. q = zeros((zn,yn,xn))
  16. q = DimArray(q, uwnd.dims, uwnd.fill_value, uwnd.proj)
  17. zidx = 0
  18. for prs in uwnd.dimvalue(0):
  19.     prs = prs * 0.001
  20.     air = f['Temperature_isobaric'][tidx,zidx,::-1,:]
  21.     rhum = f['Relative_humidity_isobaric'][tidx,zidx,::-1,:]
  22.     g = 9.8
  23.     es = 6.112*exp(17.67*(air-273.16)/(air-29.65))
  24.     qs = 0.62197*es/(prs-0.378*es)
  25.     qq = qs*rhum/100
  26.     q[zidx,:,:] = qq
  27.     zidx += 1
  28. qu = q*uwnd/g
  29. qv = q*vwnd/g
  30. qhdivg = hdivg(qu, qv)   
  31. vqhdivg = trapz(qhdivg, levels)

  32. #Plot
  33. print 'Plot...'
  34. axesm()
  35. mlayer = shaperead('D:/Temp/map/country1.shp')
  36. geoshow(mlayer, edgecolor='black')
  37. layer = imshowm(vqhdivg, 20)
  38. title('Vertical Integrated Water Vapor Flux Divergency (' + t.strftime('%Y-%m-%d') + ')')
  39. colorbar(layer)
  40. xlim(0, 360)
  41. ylim(-90, 90)


hdivg_trapz.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-8-22 15:48:35 | 显示全部楼层
本帖最后由 江南剑侠 于 2016-8-22 15:50 编辑
MeteoInfo 发表于 2016-8-22 15:35
应该是等值线填色算法的bug,将contourfm改为imshowm,



                               
登录/注册后可看大图
可以了!看来确实是这个原因呢
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表