爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 6109|回复: 5

MeteoInfoLab脚本示例:GrADS转netCDF

[复制链接]

新浪微博达人勋

发表于 2016-4-15 16:27:00 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 MeteoInfo 于 2017-4-3 11:11 编辑

convert2nc(infn, outfn)函数用来将GrADS数据转为netCDF格式数据(目前只支持等经纬度投影数据),第一个参数为欲转换的GrADS数据文件名,第二个参数为转换后的netCDF数据文件名。

脚本程序:
  1. infn = 'D:/Temp/grads/model.ctl'    #Input GrADS file
  2. outfn = 'D:/Temp/grads/model.nc'    #Output netCDF file
  3. convert2nc(infn, outfn)


也可以不用convert2nc函数,更加灵活的转换:
  1. #--------------------------------------------------------        
  2. # Author: Yaqiang Wang                                          
  3. # Date: 2016-4-12                                               
  4. # Purpose: Convert GrADS data to netCDF data  
  5. # Note: Sample                                                   
  6. #-----------------------------------------------------------

  7. infn = 'D:/Temp/grads/model.ctl'    #Input GrADS file
  8. outfn = 'D:/Temp/grads/model.nc'    #Output netCDF file

  9. #Open GrADS file
  10. f = addfile(infn)

  11. #New netCDF file
  12. ncfile = addfile(outfn, 'c')

  13. #Add dimensions
  14. dims = []
  15. for dim in f.dimensions():
  16.     dims.append(ncfile.adddim(dim.getShortName(), dim.getLength()))
  17. xdim = f.finddim('X')
  18. ydim = f.finddim('Y')
  19. tdim = f.finddim('T')
  20. xnum = xdim.getLength()
  21. ynum = ydim.getLength()
  22. tnum = tdim.getLength()

  23. #Add global attributes
  24. ncfile.addgroupattr('Conventions', 'CF-1.6')
  25. for attr in f.attributes():
  26.     ncfile.addgroupattr(attr.getName(), attr.getValues())

  27. #Add dimension variables
  28. dimvars = []
  29. for dim in dims:
  30.     dname = dim.getShortName()
  31.     if dname == 'T':
  32.         var = ncfile.addvar('time', 'int', [dim])
  33.         var.addattr('units', 'hours since 1900-01-01 00:00:0.0')
  34.         var.addattr('long_name', 'Time')
  35.         var.addattr('standard_name', 'time')
  36.         var.addattr('axis', dname)
  37.         tvar = var
  38.     elif dname == 'Z':
  39.         var = ncfile.addvar('level', 'float', [dim])
  40.         var.addattr('axis', dname)
  41.     else:
  42.         var = ncfile.addvar(dim.getShortName(), 'float', [dim])
  43.         if 'Z' in dname:
  44.             var.addattr('axis', 'Z')
  45.         else:
  46.             var.addattr('axis', dname)
  47.     dimvars.append(var)

  48. #Add variables
  49. variables = []
  50. for var in f.variables():   
  51.     print 'Variable: ' + var.getShortName()
  52.     vdims = []
  53.     for vdim in var.getDimensions():
  54.         for dim in dims:
  55.             if vdim.getShortName() == dim.getShortName():
  56.                 vdims.append(dim)
  57.     #print vdims
  58.     nvar = ncfile.addvar(var.getShortName(), var.getDataType(), vdims)
  59.     nvar.addattr('fill_value', -9999.0)
  60.     for attr in var.getAttributes():
  61.         nvar.addattr(attr.getName(), attr.getValues())
  62.     variables.append(nvar)

  63. #Create netCDF file
  64. ncfile.create()

  65. #Write variable data
  66. for dimvar, dim in zip(dimvars, f.dimensions()):
  67.     if dim.getShortName() != 'T':
  68.         ncfile.write(dimvar, array(dim.getDimValue()))

  69. sst = datetime.datetime(1900,1,1)
  70. for t in range(0, tnum):
  71.     st = f.gettime(t)
  72.     print st.strftime('%Y-%m-%d %H:00')
  73.     hours = (st - sst).total_seconds() // 3600
  74.     origin = [t]
  75.     ncfile.write(tvar, array([hours]), origin=origin)
  76.     for var in variables:
  77.         print 'Variable: ' + var.name
  78.         if var.ndim == 3:
  79.             data = f[str(var.name)][t,:,:]   
  80.             data[data==nan] = -9999.0        
  81.             origin = [t, 0, 0]
  82.             shape = [1, ynum, xnum]
  83.             data = data.reshape(shape)
  84.             ncfile.write(var, data, origin=origin)
  85.         else:
  86.             znum = var.dims[1].getLength()
  87.             for z in range(0, znum):
  88.                 data = f[str(var.name)][t,z,:,:]
  89.                 data[data==nan] = -9999.0
  90.                 origin = [t, z, 0, 0]
  91.                 shape = [1, 1, ynum, xnum]
  92.                 data = data.reshape(shape)
  93.                 ncfile.write(var, data, origin=origin)

  94. #Close netCDF file
  95. ncfile.close()

  96. print 'Finished!'


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

新浪微博达人勋

 成长值: 19710
发表于 2016-4-15 16:47:48 | 显示全部楼层
这里的GrADS数据应该是一般的二进制文件吧?用ctl描述的grib应该不包括在内?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-4-15 19:22:41 | 显示全部楼层
兰溪之水 发表于 2016-4-15 16:47
这里的GrADS数据应该是一般的二进制文件吧?用ctl描述的grib应该不包括在内?

GRIB数据转换用这个:MeteoInfoLab脚本示例:GRIB转netCDF
http://bbs.06climate.com/forum.p ... 055&fromuid=106
(出处: 气象家园)
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-5-27 19:06:32 | 显示全部楼层
王老师,
假设如下脚本
fn=W_NAFP_C_ECMF_20160520174829_P_C1D05201200052012001.grib2
f = addfile(fn)
接下来怎么获取这个grib2文件中时间的个数?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-5-28 00:02:03 | 显示全部楼层
my251394667 发表于 2016-5-27 19:06
王老师,
假设如下脚本
fn=W_NAFP_C_ECMF_20160520174829_P_C1D05201200052012001.grib2

tn = f.timenum()
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2020-11-23 17:26:17 | 显示全部楼层
老师,里面不是显示可以打开grads格式资料吗,为什么要转呢,我刚试着直接打开grads格式文件,好像出不来,是要转了才可以打开吗
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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