happy科研

關(guān)注

Hi～新朋友，記得點藍(lán)字關(guān)注我們喲

科研到頭禿(7)

學(xué)渣

今天碰到了一個科研難題，想想都頭大

學(xué)霸

學(xué)渣

模式運行得到了360個逐月輸出的文件，NCL處理起來慢爆了，想哭

學(xué)霸

哈哈哈，來來來，我給你介紹一下強大的氣象處理軟件cdo

╮(￣▽￣"")╭

所以今天主要來介紹一個超級強大的氣象類的數(shù)據(jù)處理運算軟件——cdo

目標(biāo)：快速實現(xiàn)數(shù)據(jù)的預(yù)處理等

劃重點

01

cdo簡介

cdo是一款氣象領(lǐng)域基于Linux處理數(shù)據(jù)十分強大的工具，是climate data operator的縮寫。它提供了600多個常見的操作，能夠?qū)?shù)據(jù)進(jìn)行快速的操作和分析，能夠很快速的處理nc、grid等常見的數(shù)據(jù)。常見的功能包括：

1、數(shù)據(jù)的提取合并(提取特定時間、空間、經(jīng)緯度等等)

2、數(shù)據(jù)的簡單運算(加減乘除、方差、均方差、和、最值、滑動均值、滑動方差、滑動最值、區(qū)域平均、區(qū)域方差、區(qū)域最值等等)

3、數(shù)據(jù)的統(tǒng)計運算(相關(guān)、線性回歸、EOF、濾波、水平插值、垂直插值等等)

4、數(shù)據(jù)的轉(zhuǎn)換(binary轉(zhuǎn)nc、HDF轉(zhuǎn)nc等等)

5、各種氣候指數(shù)的運算(極端有關(guān)的指數(shù)等等)

有很多軟件都可以處理氣象數(shù)據(jù)，常見的向Matlab、Python和NCL等，除此之外也有快速處理氣象數(shù)據(jù)的軟件如Cdo、NCO等。那么如果把Cdo與傳統(tǒng)的氣象軟件NCL做對比，它有如下的優(yōu)缺點。

優(yōu)點：

1. 數(shù)據(jù)處理的速度極快

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

....

缺點：

1.與NCL一樣都是的基于linux系統(tǒng)才能操作

2. 不能中途查看數(shù)據(jù)，而且是交互式命令，不利于查錯【但是可以把命令批量寫在bash腳本里面，然后執(zhí)行】。

3. 參考資料和官網(wǎng)信息沒有NCL豐富。

.....

但是“唯快不破”是最大的優(yōu)勢，這可以給我們數(shù)據(jù)處理節(jié)省很多的時間，特別是處理模式運行的結(jié)果。然后再結(jié)合NCL進(jìn)行后續(xù)的處理和分析，能夠在一定程度上提高效率。

劃重點

02

cdo安裝

推薦兩種安裝方法：

1.參考以下網(wǎng)站安裝，但是比較麻煩，不推薦。

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

2. 使用conda安裝，推薦

可以基于自己的linux子系統(tǒng)或者是服務(wù)器，輸入命令：

> conda install cdo

如果還沒有在電腦上安裝linux子系統(tǒng)的，可以參考這一篇推送里面的教程安裝哦(Windows下也可以安裝Pyngl和Pynio，你還不知道？)

劃重點

03

cdo功能

3.1

查看數(shù)據(jù)信息

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

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

3. diff, diffn? ? ? ? ? ? ? ## 查看兩個文件的差異

4. npar, nlevel, nyear, nmon, ndate, ntime? ?## 顯示文件中變量個數(shù)、層數(shù)、年、月等的長度

5.showformat, showcode, showname, showstdname, showlevel, showtype, showyear, showmon, showdate, showtime, showtimestamp? ? ? ? ? ? ?## 顯示格式信息，變量名，具體的時間信息等

6. pardes, griddes, zaxisdes, vct? ? ?##顯示網(wǎng)格信息

例1：查詢下載的ERSSTv4文件中的變量個數(shù)，年份和月份長度

代碼：

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文件的經(jīng)緯度網(wǎng)格信息

代碼：

cdo griddes slp.mon.mean.nc

注：這個功能在經(jīng)常與差異函數(shù)結(jié)合起來，對數(shù)據(jù)進(jìn)行插值，詳情見后面的介紹

3.2

數(shù)據(jù)切片

1.selname, sellevel, …? ? ?? ## 特定的變量和高度場的切片

2. seltimestep, seldate, selmon, …? ? ## 根據(jù)時間信息進(jìn)行篩選

3.sellonlatbox, …? ? ? ? ? ? ? ?##? 根據(jù)經(jīng)緯度信息切片

例5：根據(jù)時間信息選擇2000年1月到2002年12月的GPCP降水?dāng)?shù)據(jù)。

代碼：

步驟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是原始的數(shù)據(jù)，它的時間格式的“年-月-01”，GPCP_00_02.nc是切片后的文件

例6：對于NCEP的風(fēng)場數(shù)據(jù)，篩選1980-2010年中國地區(qū)(10-55N, 70-135E)的夏季850hPa的風(fēng)場。

分析：這個問題涉及到多步的操作，可以一次性的在cdo里面執(zhí)行，但是注意兩點：①設(shè)計到多步操作時，每個函數(shù)名稱前面加上“-”；②不同函數(shù)名稱中間用空格間隔。

代碼：

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

數(shù)據(jù)修改

1.?setname, setlevel, …? ?## 設(shè)置更改變量、高度場的屬性

2. settaxis, setcalendar, …? ## 更改時間

3. chname, …? ? ? ? ? ? ? ? ? ?## 更改變量名

4. setgrid, …, setzaxis, setgatt, … ## 更改經(jīng)緯度坐標(biāo)信息

5. invertlat, invertlev, …? ? ?## 經(jīng)緯度高度場導(dǎo)致，如1000-10hPa轉(zhuǎn)變成10-1000hPa的排列

6. maskregion, masklonlatbox, …## mask得到想要的區(qū)域

7. setclonlatbox, …

9. setmissval, setctomiss, setmisstoc, setrtomiss, … ## 設(shè)置為缺測值

例7：對例6中生成的uwnd.nc數(shù)據(jù)，把時間改為從1000-01-01開始，并且時間的間隔為1年

代碼：

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

例8：對于ERSST資料，單位是kelvin，講0~273.15范圍的值都設(shè)置為缺測。

代碼：

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

3.4

文件處理

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

2. replace? ? ? ? ? ? ? ? ? ? ? ? ## 替代數(shù)據(jù)

3. merge, mergetime? ? ? ?## 合并數(shù)據(jù)

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

例9：將兩個不同時間段的降水進(jìn)行合并

代碼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：現(xiàn)有11618個文件，從1979年1月1日到2010年12月31日，格式為 GLOBAL19870405.nc 每個文件有三個變量ABC。現(xiàn)在需要整合A變量的每年4月1號到9月30號的數(shù)據(jù)。

代碼：

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

解釋：這里用到了正則化的方法進(jìn)行篩選，[1-2]表示取1或者2，[90]表示取0或者9，[0-9]表示取0到9中的任意一個數(shù)字，*則表示滿足文件名的前面是GLOBAL[1-2][90][0-9][0-9]0[4-9]的所有后續(xù)文件。

3.5

數(shù)學(xué)運算

1. expr, exprf? ? ? ? ? ? ? ? ? ? ? ? ## 執(zhí)行語句的運算

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

3. addc, subc, mulc, divc,? ? ## 加減一個常數(shù)

4. add, sub, mul, …? ? ? ? ? ? ?## 兩個文件對應(yīng)值加減

5. monadd, monsub, monmul, mondiv

6. ymonadd,?ydayadd, yhouradd, …

9. muldpm, …

例11：將HadISST資料的單位從K轉(zhuǎn)變到℃。

代碼1：

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

代碼2：

cdo addc,-273.15 HadISST.nc?HadISST_C.nc

解釋：這expr就是代表后面講出現(xiàn)一個運算的表達(dá)式，并執(zhí)行后面的表達(dá)式，sst是原始HadISST.nc里面的變量名，經(jīng)過這個運算之后會生成新的變量名TSC，并存儲到HadISST_C.nc，這種方法可以只對特定的變量進(jìn)行操作。但是addc就很簡單粗暴，對所有的變量都減去了273.15

3.6

統(tǒng)計運算

1. ensmean, ensstd, ensmin, ensmax, enspctl, …? ## 集合平均、標(biāo)準(zhǔn)差等操作

2. fldmean; zonmean; mermean; gridboxmean, …? ? ## 空間區(qū)域平均、緯向平均等

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

4. timselmean,? …? ? ##? 時間序列上的平均

5. runmean, …? ? ? ? ?## 滑動平均

6. timmean, yearmean, monmean,? ?## 整個時間段平均、年平均、月平均

7. ymonmean, …? ? ? ## 特定月份，所有年數(shù)的平均，如1981到2010年所有1月的平均計算得到30年的1月平均值

例12：如何計算夏季平均？冬季平均？以及季節(jié)性平均的序列？(如3-5月的平均作為第一個值，6-8月的平均作為第二個值，....)

代碼1：夏季平均

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

解釋：這里的思路是先篩選得到6,7,8月的數(shù)據(jù)，然后再計算年平均，這樣就是夏季平均的數(shù)據(jù)了。

代碼2：冬季平均

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

解釋：這里不能用yearmean和selmon來操作了，因為冬季的平均跨越了兩個年份。因此這里用timselmean函數(shù)，設(shè)計三個參數(shù)，第一個3表示3個月的平均，第二個11表示跳過最初的11個月，第三個9表示間隔的時間步長，即從2月份到12月份需要間隔9個月。

代碼3：季節(jié)性平均

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

解釋：與冬季平均類似，3是表示3個月的平均，2表示跳過最初的2個月，因為這是屬于上一年的冬季，缺少第三個參數(shù)，是因為這里間隔為0。3-5為季節(jié)性平均的第一個值，6-8為季節(jié)性平均的第二個值。

例13：模式運行得到了逐月的輸出結(jié)果，從MMEAN1930-01.nc到MMEAN1990-02.nc，現(xiàn)在需要對后30年的結(jié)果求集合平均，然后進(jìn)行畫圖分析。

代碼：

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

解釋：這里同樣是用到了正則化的表達(dá)式，MMEAN19[6-8]* 表示滿足文件名為MMEAN196，MMEAN197，MMEAN198的所有文件。

3.7

專業(yè)技巧

1.?fldcor, timecor, fldcovar, timecovar? ?## 空間相關(guān)、時間序列相關(guān)、空間序列協(xié)方差、時間序列協(xié)方差

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? ? ? ? ? ? ? ? ? ? ?## 時間、高度場插值

相關(guān)插值函數(shù)的介紹：

例14：把ERSSTV4里面的SST插值成與GPCP的Pr具有相同分辨率的數(shù)據(jù)。

步驟一代碼：

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

解釋：這里用griddes主要是為了獲取precip變量的空間分辨率等信息。> 這個是重定向符號，然后把這些信息保存到txt當(dāng)中。

步驟二代碼：

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

解釋：這里是先讀取txt中的分辨率信息，然后remapbil按照這個信息對sst.mon.nc進(jìn)行插值得到sst_low.nc

劃重點

04

實例分享

例15：篩選CMIP5中的數(shù)據(jù)中1850-01-01到1899-12-31范圍內(nèi)的數(shù)據(jù)，并求平均后輸出

代碼：

如下所示，講這些代碼放在了bash腳本中，可以復(fù)制下面代碼放在新建的xx.sh文件里面，然后執(zhí)行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之間的數(shù)據(jù)，然后對時間維求平均cdo -timmean -seldate,1850-01-01,1899-12-31 $INPATH/$INFILE $OUTPATH/tas_Amon_${MODEL}_historical_r1i1p1.ncdone

例16：現(xiàn)在需要把觀測場中的五個變量(比濕、地面氣壓、溫度、風(fēng)場)替換到模式場中，然后驅(qū)動模式運行。模式場中對應(yīng)的為逐日資料，變量分別為Q, PS, T, US, VS，處理PS外都為四維的數(shù)據(jù)。觀測場的資料則是逐6小時的，而且每個變量的分辨率與模式的數(shù)據(jù)都不對應(yīng)。

分析：解決這個問題，有以下幾點要注意

觀測資料需要進(jìn)行日平均

要對每個觀測資料的變量進(jìn)行水平插值和垂直插值

修改觀測資料中的變量名，變?yōu)镼, PS, T, US, VS

代碼如下

#!/bin/shinputpath=original_data/outputpath=interpolate_data_f19vars=(Q PS T US VS)## 這里five_varible_f19.nc是模式場中整合下來的五個變量的合集## 這一步是獲取高度場的氣壓值，方便后面進(jìn)行插值## 本來原始得到的氣壓值是 [3.267 76.98 189.09]，sed的處理就是讓其變?yōu)閇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) ## 講模式場中相應(yīng)變量的經(jīng)緯度信息輸出到txt，方便對觀測資料進(jìn)行插值 cdo -griddes -select,name=${vars[i]} five_varible_f19.nc > horizontal_interplo.txt ## 因為這里PS是三維的，不需要垂直插值，所以這里判斷一下 if [ $i -ne 1 ] ## 這里intlevel是進(jìn)行垂直插值，remapbil是水平看見插值，daymean是求日平均，setname是修改變量名與模式里面一致 ## ${outputpath}/$(basename $filepath .nc).preprocess.nc 表示講修改后的數(shù)據(jù)保存到 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. 汪君老師的課件

在此表示感謝。

如果覺得這個有幫助，麻煩各位小可愛點一下廣告，謝謝啦。

總結(jié)

以上是生活随笔為你收集整理的输出nc数据_气象数据处理的火箭加速器—CDO的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。