Tools

OpenGrADS

配合GrADS 2.0的一些新特性,使GrADS的使用越发的简单与方便。2.0的新特性包括:

o HDF4 and HDF5 interfaces ENABLED
o GeoTIFF and KML/TIFF output ENABLED
o KML contour output ENABLED
o Shapefile interface ENABLED
o NetCDF output ENABLED现在支持shapefile ,那么中国地图的问题就可以简单的解决,不用再用纠结的cnworld 等,
只需用一个好点的shapefile就好了。支持KML输出,那么和 Google earth或者map等整合也比较容易了。支持输出NetCDF,那么保存中间和最终结果数据就更方便了。
==================================================而OpenGrADS的扩展就更多了
如(q udx):
———– —————————————  ————————–
 ftest       F-test                                   f_bjt@^libbjt.gex
 ttest       T-test                                   f_bjt@^libbjt.gex
 tfit        Point  linear regression                 f_bjt@^libbjt.gex
 fit         Global linear regression                 f_bjt@^libbjt.gex
 tcorr2      Time correlation                         f_bjt@^libbjt.gex
 tregr2      Point linear regression                  f_bjt@^libbjt.gex

 smth2d      Shuman smoother/de-smoother              f_smth2d@^libmf.gex
 xy2rt       cartesian -> cylindrical coord           f_xy2rt@^libmf.gex
 uv2trw      Find radial/tangential velocity          f_uv2trw@^libmf.gex
 re2         General interpolator (regrid2)           f_re2@^libmf.gex
 esmrf       Saturation vapor pressure (old MRF)      f_esmrf@^libmf.gex
 linreg      Linear regression: mx + b                f_linreg@^libmf.gex

 re          General interpolator                     ffre@^re.gex
 re_         General interpolator (verbose)           ffre_dbg@^re.gex
 regrid2     regrid2-like wrapper for re()            f_gsudf@^gsudf.gex
 geos        Interpolation to GEOS-5 Grids            f_gsudf@^gsudf.gex
 reimg       Interpolation for Image Generation’      f_gsudf@^gsudf.gex
——————————————————————————
简单的例子,要将一个分辨率的数据插值到另一个分辨率只需用 re 扩展函数就行了。
如 ts 是 2×2 水平分辨率的变量,要插值到 1×1的tsh,则可以简单的
define  tsh=re(ts, 1)
如果要把插值好的数据保存,那么只需
set sdfwrite  tsh.nc
sdfwrite tsh

GrADS/PyGrADS:

pygrads installation

On Windows you can get the pre-compiled basemap package from:
http://sourceforge.net/projects/matplotlib/files/matplotlib-toolkits/
basemap-1.0.1/
Get the package that matches the EPD Free version. Currently this is
 Python 2.7, so get basemap-1.0.1.win32-py2.7.exe
<http://sourceforge.net/projects/matplotlib/files/matplotlib-toolkits/
basemap-1.0.1/basemap-1.0.1.win32-py2.7.exe/download>
After that, install pygrads proper:
pygrads-1.1.b5.win32.exe<http://sourceforge.net/projects/opengrads/
files/python-grads/1.1.b5/pygrads-1.1.b5.win32.exe/download>
And of course, make sure you have the opengrads Win32 superpack 
installed:grads-2.0.a9.oga.1-win32_superpack.exe
<http://sourceforge.net/projects/opengrads/files/grads2-windows/
2.0.a9.oga.1/grads-2.0.a9.oga.1-win32_superpack.exe/download>
PyGrADS should be available on your [Start] menu.

Several examples:
----------------------------------------------------------------------
使用PyGrADS保存和导出变量:import grads
ga = grads.GaNum(Bin='grads',Echo=False,Window=False)  # start grads
fh = ga.open("http://monsoondata.org:9090/dods/model") # open file
ts = ga.expr('ts')                       # export a variable to python
print "   Lon     Lat      Ts"
print "-------- -------- --------"
for j in range(fh.ny):
    for i in range(fh.nx):
        print "%8.3f %8.3f %8.2f"%(ts.grid.lon[i],ts.grid.lat[j],ts[j,i])

输出如下:
Lon      Lat      Ts
-------- -------- --------
   0.000  -90.000   258.49
   5.000  -90.000   258.49
  10.000  -90.000   258.49
  15.000  -90.000   258.49
  20.000  -90.000   258.49
  25.000  -90.000   258.49
  ...
 230.000  -18.000   300.49
 235.000  -18.000   299.99
 240.000  -18.000   299.49
 245.000  -18.000   298.99
 ...
 --------------------------------------------------------------------------
利用PyGrADS绘图
view plain
  1. from grads import *
  2. from sys   import stdout

测试PyGrADS是否正确安装:

  1. try:
  2.     ga = GrADS(Verb=1, Echo=False, Port=False, Window=False,Opts=
  3.     “-c ‘q config'”)
  4.     print(“>>> OK <<< start GrADS”)
  5. except:
  6.     print(“>>> NOT OK <<< cannot start GrADS”)

尝试打开一个控制文件,打印文件中的信息:

  1. try:
  2.     fh = ga.open(“../data/model.ctl”)
  3.     print(fh.title)
  4.     print(‘              File Id: ‘, fh.fid)
  5.     print(‘         Dataset type: ‘, fh.type)
  6.     print(‘    No. of Time steps: ‘, fh.nt)
  7.     print(‘    No. of Longitudes: ‘, fh.nx)
  8.     print(‘    No. of  Latitudes: ‘, fh.ny)
  9.     print(‘    No. of     Levels: ‘, fh.nz)
  10.     print(‘      Variable  names: ‘, fh.vars)
  11.     print(‘      Variable levels: ‘, fh.var_levs)
  12.     print(‘      Variable titles: ‘, fh.var_titles)
  13.     print(“>>> OK <<< open CTL file”)
  14. except:
  15.     print(“>>> NOT OK <<< cannot open CTL file”)

读取NetCDF文件时,会自动调用sdfopen:

  1. try:
  2.     fh = ga.open(“../data/model.nc”)
  3.     print(fh.title)
  4.     print (‘              File Id: ‘, fh.fid)
  5.     print (‘         Dataset type: ‘, fh.type)
  6.     print (‘    No. of Time steps: ‘, fh.nt)
  7.     print (‘    No. of Longitudes: ‘, fh.nx)
  8.     print (‘    No. of  Latitudes: ‘, fh.ny)
  9.     print (‘    No. of     Levels: ‘, fh.nz)
  10.     print (‘      Variable  names: ‘, fh.vars)
  11.     print (‘      Variable levels: ‘, fh.var_levs)
  12.     print (‘      Variable titles: ‘, fh.var_titles)
  13.     print (“>>> OK <<< open NetCDF file”)
  14. except:
  15.     print(“>>> NOT OK <<< cannot open NetCDF file”)

查询维度信息(query dims):

  1. try:
  2.     qh = ga.query(‘dims’)
  3.     print(‘Current dimensional state: ‘)
  4.     print (‘   X is ‘+qh.x_state+‘   Lon = ‘,qh.lon,‘  X = ‘,qh.x)
  5.     print(‘   Y is ‘+qh.y_state+‘   Lat = ‘,qh.lat,‘  Y = ‘,qh.y)
  6.     print(‘   Z is ‘+qh.z_state+‘   Lev = ‘,qh.lev,‘  Z = ‘,qh.z)
  7.     print(‘   T is ‘+qh.t_state+‘  Time = ‘,qh.time,‘  T = ‘,qh.t)
  8.     print(“>>> OK <<< query dimensions”)
  9. except:
  10.     print(“>>> NOT OK <<< cannot query dimensions”)

查询文件信息(query file): 

  1. try:
  2.     qh = ga.query(‘file’)
  3.     print(‘Current file state: ‘)
  4.     print(‘         Title: ‘, qh.title)
  5.     print(‘       File Id: ‘, qh.fid)
  6.     print(‘   Description: ‘, qh.desc)
  7.     print(‘        Binary: ‘, qh.bin)
  8.     print ‘   Variable  Num levels  Description’)
  9.     for (var, nlevels0,nlevels1, desc) in qh.var_info:
  10.          print‘   ‘+var.rjust(8)+‘  ‘+str(nlevels0).rjust(10)+str(nlevels1).rjust(10)+‘  ‘+desc)
  11.     print(“>>> OK <<< query file”)
  12. except:
  13.     print(“>>> NOT OK <<< cannot query file”)

测试捕获grads输出(Line和Word):

  1. try:
  2.     ga.cmd(“q config”)
  3.  
  4.     print(“————————————————————–“)
  5.     print ”            Captured GrADS output: Line interface”
  6.     print(“————————————————————–“)
  7.     for i in range(1,ga.nLines):
  8.         print(ga.rline(i))
  9.  
  10.     print(“————————————————————–“)
  11.     print ”            Captured GrADS output: Word interface”
  12.     print(“————————————————————–“)
  13.     for i in range(1,ga.nLines):
  14.         for j in range(1,20):     # 20 is an over estimate, but i is OK
  15.             stdout.write(ga.rword(i,j)+‘ ‘)
  16.         stdout.write(‘\n’)
  17.  
  18.     print(“>>> OK <<< rline()/rword() completes”)
  19. except:
  20.     print(“>>> NOT OK <<< rline()/rword() fails”)

测试和Python交换变量,即可以将grads空间中的变量导入到python工作空间中,反之亦然:

  1. try:
  2.     ts = ga.exp(‘ts’)
  3.     print(“Ts in Kelvins: “, ts.min(), ts.max())
  4.     ts = ts – 273
  5.     print(“Ts in Celsius: “, ts.min(), ts.max())
  6.     ga.imp(‘tc’,ts)
  7.     tc = ga.exp(‘tc’)
  8.     print(“Tc in Celsius: “, tc.min(), tc.max())
  9.     print(“>>> OK <<< exp()/imp() completes”)
  10. except:
  11.     print(“>>> NOT OK <<< exp()/imp() fails”)
  1. from grads.gacore import GaCore
  2. ga = GaCore(Bin=‘grads’# 如果GrADS的版本小于2.0,可以使用ga = GaCore(Bin=’gradsnc’)
  3. fh = ga.open(“http://monsoondata.org:9090/dods/model”)
  4. ga(“display ps”)


实际上,在GrADS大于2.0的版本中,可以更简洁一点:

  1. import grads #导入库
  2. ga=grads.GrADS() #使用GrADS客户端程序,要求GrADS版本大于2.0.0
  3. fh = ga.open(“http://monsoondata.org:9090/dods/model”#读取文件
  4. ga(“display ps”)   #显示图形
  5. ga(“printim D:/x.png white”#保存图形,将背景设为白色

 

绘制的图形如下:

如果你安装了Numpy,那么可以很方便地在python和GrADS客户端之间传递数据:

  1. from grads.ganum import GaNum
  2. ga = GaNum()
  3. ts = ga.exp(“ts”)  # export variable ts from GrADS
  4. ts = ts – 273      # convert ts to Celsius
  5. ga.imp(“tc”,ts)    # send the NumPy array to GrADS
  6. ga(“display tc”)   # display the just imported variable

如果你还安装了MatPlotLib,那么绘制图形将非常轻松方便:

  1. from pylab import contourf
  2. contourf(ts.grid.lon,ts.grid.lat,tc)

如果你还安装了Matplotlib/Basemap,那么PyGrADS还提供了galab类,容许你使用basemap绘图:

  1. from grads.galab import GaLab
  2. ga = GaLab()
  3. fh = ga.open(“http://monsoondata.org:9090/dods/model”#读取文件
  4. ga.blue_marble(‘on’)
  5. ga(“set lon -180 180”)
  6. ga.contour(‘ua’)
  7. title(‘Zonal Wind’)

绘制的图形如下:

  1. import grads
  2. from pylab import title, savefig
  3. ga = grads.GrADS()
  4. ga.open(‘model.ctl’)
  5. ga.contour(‘ts’)
  6. title(‘Surface Temperature’)
  7. savefig(‘ts.png’)

如果将GrADS、MatPlotLib、Pylab结合起来威力更加强大。在GrADS中做统计分析是很不方便的,但是Python中却非常轻松。

例如计算eof(empirical ortoghonal functions):

  1. ga.open(“slp.nc”)
  2. ga(“set t 1 41”)
  3. v, d, pc = ga.eof(‘slp’)

执行批处理:

  1. try:
  2.     ga(“””
  3.            set lat 30 60
  4.            set lon -80 -50
  5.            set t 1 3
  6.            query dims
  7.        “””)
  8. print(“>>> OK <<< successfuly ran several commands at one”)
  9. except:
  10. print(“>>> NOT OK <<< could not run several commands at once”)

CDO

数据提取:几个简单的例子

假设:提取若干年的逐日冬季资料(这个对大家来说是个比较常见的问题)
困难:冬季跨年且存在闰年的问题
情景:1948-2013年逐日高度场资料:hgt.1948.nc, hgt.1949.nc , …. , hgt.2013.nc
其他解决办法:NCL、Fortran、GrADS可以批量读取并且去掉闰年,或者有条件判断,脚本较长
cdo解决办法如下:
1.将所有年份合并到一起:(参见UsersGuide P25)
cdo  cat  hgt.*.nc  hgt.1948-2013.nc
生成的hgt.1948-2013.nc为1948-2013年所有时次的整合nc文件。
2.提取所有年份的12,1,2月:(参见UsersGuide P36)
cdo  selmon,1,2,12  hgt.1948-2013.nc  hgt.1948-2013.1212.nc
生成的hgt.1948-2013.1212.nc为1948-2013年所有12-2月的整合nc文件。
3.由于冬季的第一年是从1948年12月开始,那么1948年的1,2月以及2013年的12月份(假设数据到2013年12月份)需要剔除:(参见UsersGuide P36)
cdo  seldate,1948-12-01,2013-02-28  hgt.1948-2013.1212.nc  hgt.1948-2012.DJF.nc
生成的hgt.1948-2012.DJF.nc 即为1948/1949至2012/2013的逐日冬季资料。
4. 如果提取djfmam平均值的时间序列,假设你有d,j,f,m,a,m6个月的时间序列数据,则可用timselmean命令,具体如下:cdo timselmean,6 ifile ofile.

5. For adding two nc files: cdo add ifile1 ifile2 ofile  %where ifile 1 and ifile2 are the input files to be added and ofile is the output file from CDO.

6. For subtracting two nc files: cdo sub ifile1 ifile2 ofile   %where ifile 1 and ifile2 are the input files to be added and ofile is the output file from CDO.
7.For converting an nc file to grib format: cdo -f grb copy filename.nc filename.grb
8. creating the ctl file for the converted grib file.
cdo gradsdes2 filename.grb  %which will create a descriptor file for grads to read.
9. selecting months from an nc file.
cdo selmon,1,2,3 max2013.nc max_jan_mar_2013.nc  
10. Finding out daily maximum and daily minimum of a parameter in nc file.
cdo daymax 2016.nc max2016.nc
cdo daymin 2016.nc min2016.nc
11. concatenating and merging several files.
cdo -b F64 -f nc2 mergetime 2012.nc 2013.nc 2014.nc 2015.nc 2016.nc out.nc    
or cdo cat 2012.nc 2013.nc 2014.nc 2015.nc 2016.nc out.nc
12. Selecting a variable from a file and writing into another file.
cdo -selname,mx2t 2017.nc max.nc   %where mx2t is the variable name to be selected from the file 2017.nc and the output file max.nc only contains the mx2t.
13. Finding out the spatial trend of a variable from a file….
cdo trend ifile.nc trend1.nc trend2.nc   %where ifile is the input file and trend1.nc is the mean file and trend2.nc is the trend file.
14. Splitting a file into two time frames
For example to split the data from 1950/2005 we have to write command:    cdo selyear, 1950/2005 gdp_histsoc_0p5deg_annual_1861-2005.nc gdp_1950_2005.nc
15. For selecting a season from monthly data..
cdo select,season=JJAS dtr_mon_1950_2013.nc summer_dtr.nc
16. To calculate the climatological monthly mean of multi years, I used
cdo -f nc ymonmean input.nc output.nccdo fldcor file1 file2 output.nc
17. For spatial correlation
cdo fldcor file1 file2 output.nc
18. For temporal correlation
cdo timcor dtrw.nc gdp.nc dtrw_gdp_scorr.nc
19.  selecting a season
cdo select,season=DJF input.nc output.nc
20. regrid using cdo: If you already have a netCDF file which has data on your target grid, you can also extract the grid definition from that file:
cdo griddes my.nc > mygrid.txt
Then you can interpolate data. Assume that we do a bilinear interpolation. This is done with the remapbil operator via: cdo remapbil,mygrid INPUT_FILE.nc OUTPUT_FILE.nc
21.selecting a season
If the coordinate variables are missing, then they can be created via the setgrid operator. If we want to add a grid definition to the source file then we need the source file's grid definition in the text file (see 'define grid').

cdo setgrid,mySourceGridDef INFILE_NO_COORDS.nc INFILE_WI_COORDS.nc
22.change celendar type: cdo setcalendar,standard ifile ofile

22. change name of variable:   cdo chname, v1,v2 ifile.nc ofile.nc

EOF 分析实例:学习和使用CDO

现在有气候模式模拟的3500-4499年的月平均海面温度(tsw)、海面气压(slp)、风应力的纬向和经向分量(ustrw、vstrw)数据,分别存放在NetCDF文件tsw.nc、slp.nc、ustrw.nc、vstrw.nc中。以海面温度tsw为例进行北太平洋海域100年(3500-3600)的EOF分析。可以用命令
cdo sinfon tsw.nc
察看文件的信息。sinfon运算符返回文件的简短信息,并已适合阅读的方式呈现。类似的还有info, infon, sinfo等运算符。
1.因为cdo能同时对一个文件中的所有变量分别执行同样的操作,所以首先将所有变量合并到一个文件中以方便操作。在linux终端中输入命令
cdo merge tsw.nc slp.nc ustrw.nc vstrw.nc all.nc
其中最后一个文件是输出文件,前面的都是输入文件。merge运算符是将不同文件合并起来,要求这些文件中的不同变量定义在相同的时间轴上。如果文件中的变量相同而时间不同,则应该用copy运算符。
2.选取要研究的时间和空间范围内的数据:
cdo selyear,3500/3600 -sellonlatbox,110,-70,-20,60 all.nc all_sub.nc
sellonlatbox运算符的作用是从输入文件中选取一定的经纬度区域,写入到输出文件中。上面选取了100E:70W,20S:60N的区域,即赤道太平洋和北太平洋。
在运算符前面加短线-,表示将该运算符的输出作为前一个运算符的输入,即实现了类似与linux shell中管道的作用。如在linux shell中man ls | less表示将man ls的输出结果导入到less中作为输入,而cdo中如下命令
cdo sub -dayavg ifle2 -timavg ifile1 ofile
相当于
cdo timavg ifile1 tmp1
cdo dayavg ifile2 tmp2
cdo sub tmp2 tmp1 ofile
rm tmp1 tmp2
即将dayavg和timavg的输出结果都作为前面第一个不加-的运算符(sub)的输入。这样既避免了不必要的文件操作,又能利用cdo内部的并行操作特性从而加快运算速度。
于是,上面的命令就是将sellonlatbox的输出结果,即提取出北太平洋区域后的文件,作为selyear运算符的输入文件。相当于
cdo sellonlatbox,110,-70,-20,60 all.nc tmp.nc
cdo selyear,3500/3600 tmp.nc sub.nc
rm tmp
运算符后面的逗号引导运算符的参数列表。例如上面的sellonlatbox运算符接受4个参数,分别表示选取区域的西、东、南、北四个方向的边界。cdo中可以用负数表示西经的经度。
selyear运算符用于从输入文件中选取指定年份的数据,selyear,3500/3600 tmp sub.nc表示从tmp.nc中选取3500-3600年写入sub.nc中。3500/3600的写法在cdo中表示一个整数范围,表示从3500到3600,间隔为1.而3500/3600/10则表示每隔10年选取一年。
3. 将陆地点设为无效数据:
标准的NetCDF文件中的变量应该具有一个属性“_FillValue”,表示缺测或无效数据的值,用以识别陆地或无效数据。但这里所用的数据文件不太标准,没有这个属性。数据中陆地点上的变量值被设为了0.而进行EOF分析时0是会作为数据考虑在内的,如果提出所有=0的点,又会把确实等于零的有效数据也去掉,所以最好给数据加上无效值。
这里我们有ETOPO全球高程水深数据etopo120.nc,利用它来给数据加上陆地标记。但这个水深数据的网格跟前面用的数据不同,因此还需要先插值到数据所用的网格上。
首先生成sub.nc数据的网格描述文件:
cdo griddes sub.nc > subgrid.txt
griddes运算符生成网格描述信息,然后用linux shell的重定向特性将网格描述写入到文件中。然后进行插值:
cdo remapbil,subgrid.txt etopo120.nc etopo_interp.nc
remapbil运算符用双线性(bilinear)算法进行插值,逗号后面的参数是插值的目标网格的描述文件。除了双线性插值,cdo还有双立方插值remapbic、反距离加权插值remapdis、临近点插值remapnn、最大面积比插值remaplaf等一系列插值运算符,还能进行时间插值。
现在就可以利用插值后的etopo数据来识别陆地:
cdo ifthen -ltc,0 etopo_interp.nc sub.nc sub_masked.nc
这里再次用了管道技术。ltc,c ifile ofile 运算符用来比较输入文件ifile是否小于常量c,小于c的点赋值为1,否则为0.即
{ 1 if i(t,x)<c
o(t,x)={ 0 if i(t,x)>=c
{miss if i(t,x)=miss
其中i(t,x)表示输入文件中定义在时间t和空间x上的变量,o(t,x)表示输入变量,miss表示无效值。
ifthen运算符用于有条件地选取输入文件中的某些点写入输出文件。ifthen ifile1 ifile2 ofile表示从ifile2中选取ifile1的对应位置不为零的点,其他点设为无效值,即
o(t,x)={i2(t,x) if i1(t,x)!=0
{ miss if i1(t,x)=0
于是前面的命令从sub.nc中选取了水深小于0的点写入sub_masked.nc,而水深>=0的点则被设为了无效值。
4. 计算异常(annomally)
EOF分析要求分析的对象是异常时间序列,即减掉平均值的时间序列。用cdo可以如下实现:
cdo sub sub_masked.nc -timmean sub_masked.nc annomally.nc
其中timmean计算时间平均,sub运算符将两个文件相减。
5. 计算EOF
用cdo中专门的eof运算符进行eof计算:
cdo eof,10 annomally.nc eigenvals.nc eofs.nc
输出文件eigenvals.nc存放所有的特征值(已经按降序排序),eofs.nc存放特征向量,即EOF,参数10表示只保存前10个EOF。eofs.nc中不同的时间步就是不同的EOF,所以这里eofs.nc中有10个时间步。
cdo eofcoeff eofs.nc annomally.nc pc
这是计算EOF的时间系数,即主成分(PC)。输入文件是刚才生成的eofs.nc和原异常时间序列文件annomally.nc。因为eofs.nc中保存了前10个EOF,这里将会生成10个文件,每个文件都以pc开头,然后后面跟一个数字,表示第几个主成分。
到此为止我们已经得到了每个变量的eof、特征值和主成分了。下面就可以绘图分析了。