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

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

23. cdo select,name=var1,var2, … ,varN input output; So, if you want to extract sea surface temperature (variable name: sst) only from my_data.nc then,

cdo select,name=sst my_data.nc new_data.nc

24. for “gadsdf: Time increment too large for ‘minutes since’ time units attribute” 

cdo settunits,[seconds|minutes|days|months|years] input-filename output-filename

25. You can take values excepting the values greater than -2 degrees Celsius: cdo ifthen -gtc,-2 input.nc input.nc myregion.nc

Meaning: if first input.nc is greater than -2 then let it be second input.nc and it writes as myregion.nc

26  how to mask or extract data from nc files by one certain (e.g., China) region:

cdo remapnn,China.nc -selname,temp temp_2051-2110.nc temp_China_a.nc

cdo mul temp_China_a.nc -div China.nc China.nc temp_China.nc

27.  When you encounter a problem using mergetime such as

“cdf_put_vara_double : ncid = 196608 varid = 3 val0 = -42744.000000

Error (cdf_put_vara_double) : NetCDF: Numeric conversion not representable”

try  cdo -b F64 mergetime *.nc outfile.nc

利用CDO进行EOF分析

现在有气候模式模拟的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、特征值和主成分了。下面就可以绘图分析了。