happy科研

關(guān)注

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

科研到頭禿(7)

學(xué)渣

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

學(xué)霸

學(xué)渣

模式運(yùn)行得到了360個(gè)逐月輸出的文件，NCL處理起來(lái)慢爆了，想哭

學(xué)霸

哈哈哈，來(lái)來(lái)來(lái)，我給你介紹一下強(qiáng)大的氣象處理軟件cdo

╮(￣▽￣"")╭

所以今天主要來(lái)介紹一個(gè)超級(jí)強(qiáng)大的氣象類(lèi)的數(shù)據(jù)處理運(yùn)算軟件——cdo

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

劃重點(diǎn)

01

cdo簡(jiǎn)介

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

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

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

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

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

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

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

優(yōu)點(diǎn)：

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

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

....

缺點(diǎn)：

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

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

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

.....

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

劃重點(diǎn)

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

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

劃重點(diǎn)

03

cdo功能

3.1

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

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

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

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

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

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

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

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

代碼：

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文件的時(shí)間范圍。

代碼：

cdo showdata slp.mon.mean.nc

例4：查詢nc文件的經(jīng)緯度網(wǎng)格信息

代碼：

cdo griddes slp.mon.mean.nc

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

3.2

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

1.selname, sellevel, …? ? ?? ## 特定的變量和高度場(chǎng)的切片

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

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

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

代碼：

步驟1：cdo showdata xxx.nc

解釋：目的是為了查看時(shí)間格式信息到底是什么，可能是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ù)，它的時(shí)間格式的“年-月-01”，GPCP_00_02.nc是切片后的文件

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

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

代碼：

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è)置更改變量、高度場(chǎng)的屬性

2. settaxis, setcalendar, …? ## 更改時(shí)間

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

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

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

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

7. setclonlatbox, …

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

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

代碼：

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

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

代碼：

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

3.4

文件處理

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

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

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

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

例9：將兩個(gè)不同時(shí)間段的降水進(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個(gè)文件，從1979年1月1日到2010年12月31日，格式為 GLOBAL19870405.nc 每個(gè)文件有三個(gè)變量ABC。現(xiàn)在需要整合A變量的每年4月1號(hào)到9月30號(hào)的數(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中的任意一個(gè)數(shù)字，*則表示滿足文件名的前面是GLOBAL[1-2][90][0-9][0-9]0[4-9]的所有后續(xù)文件。

3.5

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

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

2. abs, int, nint, pow, sqr, sqrt, exp, ln, log10, sin, cos, tan, asin, acos, …? ? ? ? ? ? ? ? ? ? ? ?## 常見(jiàn)的運(yùn)算

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

4. add, sub, mul, …? ? ? ? ? ? ?## 兩個(gè)文件對(duì)應(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)一個(gè)運(yùn)算的表達(dá)式，并執(zhí)行后面的表達(dá)式，sst是原始HadISST.nc里面的變量名，經(jīng)過(guò)這個(gè)運(yùn)算之后會(huì)生成新的變量名TSC，并存儲(chǔ)到HadISST_C.nc，這種方法可以只對(duì)特定的變量進(jìn)行操作。但是addc就很簡(jiǎn)單粗暴，對(duì)所有的變量都減去了273.15

3.6

統(tǒng)計(jì)運(yùn)算

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

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

3. vertmean, …? ? ? ? ## 垂直平均，即高度場(chǎng)的平均

4. timselmean,? …? ? ##? 時(shí)間序列上的平均

5. runmean, …? ? ? ? ?## 滑動(dòng)平均

6. timmean, yearmean, monmean,? ?## 整個(gè)時(shí)間段平均、年平均、月平均

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

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

代碼1：夏季平均

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

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

代碼2：冬季平均

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

解釋：這里不能用yearmean和selmon來(lái)操作了，因?yàn)槎镜钠骄缭搅藘蓚€(gè)年份。因此這里用timselmean函數(shù)，設(shè)計(jì)三個(gè)參數(shù)，第一個(gè)3表示3個(gè)月的平均，第二個(gè)11表示跳過(guò)最初的11個(gè)月，第三個(gè)9表示間隔的時(shí)間步長(zhǎng)，即從2月份到12月份需要間隔9個(gè)月。

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

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

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

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

代碼：

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

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

3.7

專(zhuān)業(yè)技巧

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

2. regress, detrend, trend, subtrend? ? ?## 回歸，去趨勢(shì)

3. eof, eoftime, eofspatial, eof3d, eofcoeff? ?## eof計(jì)算

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

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

inttime, intntime, intyear? ? ? ? ? ? ? ? ? ? ?## 時(shí)間、高度場(chǎng)插值

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

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

步驟一代碼：

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

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

步驟二代碼：

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

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

劃重點(diǎn)

04

實(shí)例分享

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

代碼：

如下所示，講這些代碼放在了bash腳本中，可以復(fù)制下面代碼放在新建的xx.sh文件里面，然后執(zhí)行bash xx.sh即可以運(yùn)行。

#!/bin/sh -eINPATH=/data OUTPATH=/A1_1850-1899mkdir -p $OUTPATH## 列出所有input下面的文件for INFILE in `ls $INPATH`do## 獲取每個(gè)模式的名稱(chēng)MODEL=`echo $INFILE | cut -d_ -f3`## 輸出這個(gè)模式的名稱(chēng)echo $MODEL## 篩選1850-01-01到1899-12-31之間的數(shù)據(jù)，然后對(duì)時(shí)間維求平均cdo -timmean -seldate,1850-01-01,1899-12-31 $INPATH/$INFILE $OUTPATH/tas_Amon_${MODEL}_historical_r1i1p1.ncdone

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

分析：解決這個(gè)問(wèn)題，有以下幾點(diǎn)要注意

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

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

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

代碼如下

#!/bin/shinputpath=original_data/outputpath=interpolate_data_f19vars=(Q PS T US VS)## 這里five_varible_f19.nc是模式場(chǎng)中整合下來(lái)的五個(gè)變量的合集## 這一步是獲取高度場(chǎng)的氣壓值，方便后面進(jìn)行插值## 本來(lái)原始得到的氣壓值是 [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## 獲取每個(gè)觀測(cè)變量的完整路徑for filepath in `ls $inputpath*.nc`do ## 獲得每個(gè)文件里面的變量，如sst.mon.mean.nc的文件，這里的fullname就得到sst fullname= echo $(basename $filepath .nc) ## 講模式場(chǎng)中相應(yīng)變量的經(jīng)緯度信息輸出到txt，方便對(duì)觀測(cè)資料進(jìn)行插值 cdo -griddes -select,name=${vars[i]} five_varible_f19.nc > horizontal_interplo.txt ## 因?yàn)檫@里PS是三維的，不需要垂直插值，所以這里判斷一下 if [ $i -ne 1 ] ## 這里intlevel是進(jìn)行垂直插值，remapbil是水平看見(jiàn)插值，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. 汪君老師的課件

在此表示感謝。

如果覺(jué)得這個(gè)有幫助，麻煩各位小可愛(ài)點(diǎn)一下廣告，謝謝啦。

總結(jié)

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

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