happy科研

关注

Hi～新朋友，记得点蓝字关注我们哟

科研到头秃(7)

学渣

今天碰到了一个科研难题，想想都头大

学霸

学渣

模式运行得到了360个逐月输出的文件，NCL处理起来慢爆了，想哭

学霸

哈哈哈，来来来，我给你介绍一下强大的气象处理软件cdo

╮(￣▽￣"")╭

所以今天主要来介绍一个超级强大的气象类的数据处理运算软件——cdo

目标：快速实现数据的预处理等

划重点

01

cdo简介

cdo是一款气象领域基于Linux处理数据十分强大的工具，是climate data operator的缩写。它提供了600多个常见的操作，能够对数据进行快速的操作和分析，能够很快速的处理nc、grid等常见的数据。常见的功能包括：

1、数据的提取合并(提取特定时间、空间、经纬度等等)

2、数据的简单运算(加减乘除、方差、均方差、和、最值、滑动均值、滑动方差、滑动最值、区域平均、区域方差、区域最值等等)

3、数据的统计运算(相关、线性回归、EOF、滤波、水平插值、垂直插值等等)

4、数据的转换(binary转nc、HDF转nc等等)

5、各种气候指数的运算(极端有关的指数等等)

有很多软件都可以处理气象数据，常见的向Matlab、Python和NCL等，除此之外也有快速处理气象数据的软件如Cdo、NCO等。那么如果把Cdo与传统的气象软件NCL做对比，它有如下的优缺点。

优点：

1. 数据处理的速度极快

2. 文件很小，基本上不占空间

....

缺点：

1.与NCL一样都是的基于linux系统才能操作

2. 不能中途查看数据，而且是交互式命令，不利于查错【但是可以把命令批量写在bash脚本里面，然后执行】。

3. 参考资料和官网信息没有NCL丰富。

.....

但是“唯快不破”是最大的优势，这可以给我们数据处理节省很多的时间，特别是处理模式运行的结果。然后再结合NCL进行后续的处理和分析，能够在一定程度上提高效率。

划重点

02

cdo安装

推荐两种安装方法：

1.参考以下网站安装，但是比较麻烦，不推荐。

http://www.studytrails.com/blog/install-climate-data-operator-cdo-with-netcdf-grib2-and-hdf5-support/

2. 使用conda安装，推荐

可以基于自己的linux子系统或者是服务器，输入命令：

> conda install cdo

如果还没有在电脑上安装linux子系统的，可以参考这一篇推送里面的教程安装哦(Windows下也可以安装Pyngl和Pynio，你还不知道？)

划重点

03

cdo功能

3.1

查看数据信息

1. info, infon, map   ## 查看文件基本信息

2. sinfo, sinfon        ## 查看文件基本信息

3. diff, diffn              ## 查看两个文件的差异

4. npar, nlevel, nyear, nmon, ndate, ntime   ## 显示文件中变量个数、层数、年、月等的长度

5.showformat, showcode, showname, showstdname, showlevel, showtype, showyear, showmon, showdate, showtime, showtimestamp             ## 显示格式信息，变量名，具体的时间信息等

6. pardes, griddes, zaxisdes, vct     ##显示网格信息

例1：查询下载的ERSSTv4文件中的变量个数，年份和月份长度

代码：

cdo nyear ERSSTV4.sst.mon.mean.2.202001.nc

cdo nmon ERSSTV4.sst.mon.mean.2.202001.nc

cdo npar ERSSTV4.sst.mon.mean.2.202001.nc

例2：查询nc文件的基本信息

代码：

cdo sinfon precip.mon.mean.nc

例3：查询nc文件的时间范围。

代码：

cdo showdata slp.mon.mean.nc

例4：查询nc文件的经纬度网格信息

代码：

cdo griddes slp.mon.mean.nc

注：这个功能在经常与差异函数结合起来，对数据进行插值，详情见后面的介绍

3.2

数据切片

1.selname, sellevel, …       ## 特定的变量和高度场的切片

2. seltimestep, seldate, selmon, …    ## 根据时间信息进行筛选

3.sellonlatbox, …               ##  根据经纬度信息切片

例5：根据时间信息选择2000年1月到2002年12月的GPCP降水数据。

代码：

步骤1：cdo showdata xxx.nc

解释：目的是为了查看时间格式信息到底是什么，可能是1979-01-01,1979-02-01.....也可能是1979-01,1979-02等等

步骤2：

cdo seldate,2000-01-01,2002-12-01 pre.mon.nc GPCP_00_02.nc

解释：这里pre.mon.nc是原始的数据，它的时间格式的“年-月-01”，GPCP_00_02.nc是切片后的文件

例6：对于NCEP的风场数据，筛选1980-2010年中国地区(10-55N, 70-135E)的夏季850hPa的风场。

分析：这个问题涉及到多步的操作，可以一次性的在cdo里面执行，但是注意两点：①设计到多步操作时，每个函数名称前面加上“-”；②不同函数名称中间用空格间隔。

代码：

cdo -selmon,6,7,8 -seldate,1981-01-01,2010-12-01 -sellonlatbox,70,135,10,55 uwnd.mon.mean.nc uwnd.nc

解释：这里uwnd.mon.mean.nc是原始文件，uwnd.nc是处理之后生成的文件

3.3

数据修改

1. setname, setlevel, …   ## 设置更改变量、高度场的属性

2. settaxis, setcalendar, …  ## 更改时间

3. chname, …                   ## 更改变量名

4. setgrid, …, setzaxis, setgatt, … ## 更改经纬度坐标信息

5. invertlat, invertlev, …     ## 经纬度高度场导致，如1000-10hPa转变成10-1000hPa的排列

6. maskregion, masklonlatbox, …## mask得到想要的区域

7. setclonlatbox, …

9. setmissval, setctomiss, setmisstoc, setrtomiss, … ## 设置为缺测值

例7：对例6中生成的uwnd.nc数据，把时间改为从1000-01-01开始，并且时间的间隔为1年

代码：

cdo settaxis,1000-01-01,00:00:00,1year uwnd.nc uwnd_yr.nc 

例8：对于ERSST资料，单位是kelvin，讲0~273.15范围的值都设置为缺测。

代码：

cdo setrtomiss,0,273.15 ERSSTV4.nc ERSSTV4_missing.nc

3.4

文件处理

1. copy, cat                      ## 合并，多个文件合成成一个

2. replace                        ## 替代数据

3. merge, mergetime       ## 合并数据

4. splitcode, splitname, splitlevel, splithour, splitday, splitsel                             ## 文件切分

例9：将两个不同时间段的降水进行合并

代码1：

cdo mergetime precip.200001-200012.nc precip.200101-200112.nc precip.200001-200112.nc

代码2：

cdo copy precip.200001-200012.nc precip.200101-200112.nc precip.200001-200112.nc

解释：这里precip.200001-200012.nc是2000年1月到2000年12月的降水，precip.200101-200112.nc是2001年1月到2001年12月的降水。生成的文件是precip.200001-200112.nc

例10：现有11618个文件，从1979年1月1日到2010年12月31日，格式为 GLOBAL19870405.nc 每个文件有三个变量ABC。现在需要整合A变量的每年4月1号到9月30号的数据。

代码：

cdo select,name=A GLOBAL[1-2][90][0-9][0-9]0[4-9]*  sm_1979_2010.nc

解释：这里用到了正则化的方法进行筛选，[1-2]表示取1或者2，[90]表示取0或者9，[0-9]表示取0到9中的任意一个数字，*则表示满足文件名的前面是GLOBAL[1-2][90][0-9][0-9]0[4-9]的所有后续文件。

3.5

数学运算

1. expr, exprf                        ## 执行语句的运算

2. abs, int, nint, pow, sqr, sqrt, exp, ln, log10, sin, cos, tan, asin, acos, …                       ## 常见的运算

3. addc, subc, mulc, divc,    ## 加减一个常数

4. add, sub, mul, …             ## 两个文件对应值加减

5. monadd, monsub, monmul, mondiv

6. ymonadd, ydayadd, yhouradd, …

9. muldpm, …

例11：将HadISST资料的单位从K转变到℃。

代码1：

cdo expr,’TSC=sst-273.15;’ HadISST.nc HadISST_C.nc    

代码2：

cdo addc,-273.15 HadISST.nc HadISST_C.nc

解释：这expr就是代表后面讲出现一个运算的表达式，并执行后面的表达式，sst是原始HadISST.nc里面的变量名，经过这个运算之后会生成新的变量名TSC，并存储到HadISST_C.nc，这种方法可以只对特定的变量进行操作。但是addc就很简单粗暴，对所有的变量都减去了273.15

3.6

统计运算

1. ensmean, ensstd, ensmin, ensmax, enspctl, …  ## 集合平均、标准差等操作

2. fldmean; zonmean; mermean; gridboxmean, …    ## 空间区域平均、纬向平均等

3. vertmean, …        ## 垂直平均，即高度场的平均

4. timselmean,  …    ##  时间序列上的平均

5. runmean, …         ## 滑动平均

6. timmean, yearmean, monmean,   ## 整个时间段平均、年平均、月平均

7. ymonmean, …      ## 特定月份，所有年数的平均，如1981到2010年所有1月的平均计算得到30年的1月平均值

例12：如何计算夏季平均？冬季平均？以及季节性平均的序列？(如3-5月的平均作为第一个值，6-8月的平均作为第二个值，....)

代码1：夏季平均

cdo  -yearmean –selmon,6,7,8  aaa.nc  bbb.nc

解释：这里的思路是先筛选得到6,7,8月的数据，然后再计算年平均，这样就是夏季平均的数据了。

代码2：冬季平均

cdo timselmean, 3, 11, 9   aaa.nc  bbb.nc

解释：这里不能用yearmean和selmon来操作了，因为冬季的平均跨越了两个年份。因此这里用timselmean函数，设计三个参数，第一个3表示3个月的平均，第二个11表示跳过最初的11个月，第三个9表示间隔的时间步长，即从2月份到12月份需要间隔9个月。

代码3：季节性平均

cdo timselmean, 3, 2  aaa.nc  bbb.nc

解释：与冬季平均类似，3是表示3个月的平均，2表示跳过最初的2个月，因为这是属于上一年的冬季，缺少第三个参数，是因为这里间隔为0。3-5为季节性平均的第一个值，6-8为季节性平均的第二个值。

例13：模式运行得到了逐月的输出结果，从MMEAN1930-01.nc到MMEAN1990-02.nc，现在需要对后30年的结果求集合平均，然后进行画图分析。

代码：

cdo ensmean MMEAN19[6-8]*  MEAN_60_89.nc

解释：这里同样是用到了正则化的表达式，MMEAN19[6-8]* 表示满足文件名为MMEAN196，MMEAN197，MMEAN198的所有文件。

3.7

专业技巧

1. fldcor, timecor, fldcovar, timecovar   ## 空间相关、时间序列相关、空间序列协方差、时间序列协方差

2. regress, detrend, trend, subtrend     ## 回归，去趋势

3. eof, eoftime, eofspatial, eof3d, eofcoeff   ## eof计算

4. remapbil, remapbic, remapdis, remapnn, remapcon, remapcon2, remaplaf                         ## 空间插值

5. remapeta, ml2pl, ml2hl, intlevel, intlevel3d, intlevelx3d

inttime, intntime, intyear                     ## 时间、高度场插值

相关插值函数的介绍：

例14：把ERSSTV4里面的SST插值成与GPCP的Pr具有相同分辨率的数据。

步骤一代码：

cdo -griddes -select,name=precip GPCP.mon.nc > horizontal_interplo.txt

解释：这里用griddes主要是为了获取precip变量的空间分辨率等信息。> 这个是重定向符号，然后把这些信息保存到txt当中。

步骤二代码：

cdo -remapbil,horizontal_interplo.txt sst.mon.nc sst_low.nc

解释：这里是先读取txt中的分辨率信息，然后remapbil按照这个信息对sst.mon.nc进行插值得到sst_low.nc

划重点

04

实例分享

例15：筛选CMIP5中的数据中1850-01-01到1899-12-31范围内的数据，并求平均后输出

代码：

如下所示，讲这些代码放在了bash脚本中，可以复制下面代码放在新建的xx.sh文件里面，然后执行bash xx.sh即可以运行。

#!/bin/sh -eINPATH=/data                                      OUTPATH=/A1_1850-1899mkdir -p $OUTPATH## 列出所有input下面的文件for INFILE in `ls $INPATH`do## 获取每个模式的名称MODEL=`echo $INFILE | cut -d_ -f3`## 输出这个模式的名称echo $MODEL## 筛选1850-01-01到1899-12-31之间的数据，然后对时间维求平均cdo -timmean -seldate,1850-01-01,1899-12-31 $INPATH/$INFILE  $OUTPATH/tas_Amon_${MODEL}_historical_r1i1p1.ncdone

例16：现在需要把观测场中的五个变量(比湿、地面气压、温度、风场)替换到模式场中，然后驱动模式运行。模式场中对应的为逐日资料，变量分别为Q, PS, T, US, VS，处理PS外都为四维的数据。观测场的资料则是逐6小时的，而且每个变量的分辨率与模式的数据都不对应。

分析：解决这个问题，有以下几点要注意

1. 观测资料需要进行日平均

2. 要对每个观测资料的变量进行水平插值和垂直插值

3. 修改观测资料中的变量名，变为Q, PS, T, US, VS

代码如下

#!/bin/shinputpath=original_data/outputpath=interpolate_data_f19vars=(Q PS T US VS)## 这里five_varible_f19.nc是模式场中整合下来的五个变量的合集## 这一步是获取高度场的气压值，方便后面进行插值## 本来原始得到的气压值是 [3.267 76.98 189.09]，sed的处理就是让其变为[3.267,76.98,189.09]，也是方便后面插值levels=$(echo `cdo -showlevel -select,name=US five_varible_f19.nc`|sed 's/[ ][ ]*/,/g')i=0## 获取每个观测变量的完整路径for filepath in `ls $inputpath*.nc`do  ## 获得每个文件里面的变量，如sst.mon.mean.nc的文件，这里的fullname就得到sst  fullname= echo $(basename $filepath .nc)  ## 讲模式场中相应变量的经纬度信息输出到txt，方便对观测资料进行插值  cdo -griddes -select,name=${vars[i]} five_varible_f19.nc > horizontal_interplo.txt  ## 因为这里PS是三维的，不需要垂直插值，所以这里判断一下  if [ $i -ne 1 ]    ## 这里intlevel是进行垂直插值，remapbil是水平看见插值，daymean是求日平均，setname是修改变量名与模式里面一致      ## ${outputpath}/$(basename $filepath .nc).preprocess.nc  表示讲修改后的数据保存到 output下，文件名如US.preprocess.nc等    then cdo -intlevel,$levels -remapbil,horizontal_interplo.txt -daymean -setname,${vars[i]} $filepath ${outputpath}/$(basename $filepath .nc).preprocess.nc    else cdo -remapbil,horizontal_interplo.txt -daymean -setname,${vars[i]} $filepath ${outputpath}/$(basename $filepath .nc).preprocess.nc  fi  rm horizontal_interplo.txt  i=$[i+1]  ## 使得i变成i+1done

参考资料：

1.https://code.mpimet.mpg.de/projects/cdo/embedded/cdo.pdf  (cdo官方文档)

2. 汪君老师的课件

在此表示感谢。

如果觉得这个有帮助，麻烦各位小可爱点一下广告，谢谢啦。