合并多个nc数据_气象数据处理的火箭加速器—CDO
關注
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小時的，而且每個變量的分辨率與模式的數據都不對應。
分析：解決這個問題，有以下幾點要注意
觀測資料需要進行日平均
要對每個觀測資料的變量進行水平插值和垂直插值
修改觀測資料中的變量名，變為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. 汪君老師的課件
在此表示感謝。
如果覺得這個有幫助，麻煩各位小可愛點一下廣告，謝謝啦。總結
以上是生活随笔為你收集整理的合并多个nc数据_气象数据处理的火箭加速器—CDO的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Sharepoint 2007 用代码聚
- 下一篇： Cognos 云最佳实践: 调整架构提供