早崎トップ 研究(気候気象) 研究(大気汚染) データリスト Linux Tips Mac Tips

CDO 使用方法 (hysk)

CDO とは

Climate Data Operatos の略. NetCDF 3/4, GRIB 1/2 形式のデータファイルをコマンドで操作できる. 作図機能はないが,解析の初期段階でのデータ操作に効果的なコマンド(オペレータ)が 非常に豊富.

セットアップ手順は CDO セットアップ(hysk)参照.

CDO 使用方法・実行例

入力・出力 NetCDF ファイル名をそれぞれ input.nc & output.nc とする.

コマンドラインでの使い方はこんな感じ:

cdo [オペレータ] input.nc

オペレータの種類によって利用するオプションが異なる. オペレータは非常に多種多様(400個以上!)なので, 全部解説したら長大なマニュアルになってしまう (See CDO User's Guide (html); PDF). 私自身が使用してみたものから,便利そうな一部機能のみ記載.

なお,下記の記述中で,個別項目(セクション番号など)への直接リンクは,閲覧当時のもの. マニュアル追加や修正により,参照アドレスが変化するため,リンク切れを起こす事が多い.
CDO User's Guide (html) の目次部分から該当箇所を探す(検索する)方が見付けやすいと思われる.

 使えそうなオペレータ

ファイルの中身を見る,必要なデータを切り出す,データの格子間隔を揃える,など.

 情報表示

2.1.1 Information and simple statistics: info, infon, map

  • infon
    (記録データの順番,時刻,鉛直レベル,格子数,未定義値格子数,最小・平均・最大値,変数名)
    $ cdo infon ua_day_MRI-CGCM3_historical_r1i1p1_19800101-19801231.nc
        -1 :       Date     Time   Level Gridsize    Miss :     Minimum        Mean     Maximum : Parameter name
         1 : 1980-01-01 12:00:00  100000    51200   26841 :     -22.207    -0.86837      21.827 : ua            
         2 : 1980-01-01 12:00:00   85000    51200    5212 :     -27.068      1.5460      30.109 : ua            
         3 : 1980-01-01 12:00:00   70000    51200    2308 :     -21.053      3.7809      36.637 : ua
    (中略)
      2926 : 1980-12-31 12:00:00   10000    51200       0 :     -21.816      11.593      53.787 : ua            
      2927 : 1980-12-31 12:00:00    5000    51200       0 :     -14.295      6.6742      59.838 : ua            
      2928 : 1980-12-31 12:00:00    1000    51200       0 :     -28.478      4.8907      98.844 : ua
    
    (入力データ例: CMIP5, historical, MRI-CGCM3, r1i1p1, daily mean u, 1980年)
    

2.1.2 Short information: sinfo, sinfon

  • sinfo (short information)
    (格子点,変数,時刻などの概要表示)
    $ cdo sinfo ua_day_MRI-CGCM3_historical_r1i1p1_19800101-19801231.nc
       File format : netCDF
        -1 : Institut Source   Ttype    Levels Num    Points Num Dtype : Parameter name
         1 : unknown  MRI-CGCM3 instant       8   1     51200   1  F32  : ua            
       Grid coordinates :
         1 : gaussian                 : points=51200 (320x160)  np=80
                                  lon : 0 to 358.875 by 1.125 degrees_east  circular
                                  lat : -89.1415 to 89.1415 degrees_north
                            available : cellbounds
       Vertical coordinates :
         1 : pressure                 : levels=8
                                 plev : 100000 to 1000 Pa
       Time coordinate :  366 steps
         RefTime =  1850-01-01 00:00:00  Units = days  Calendar = standard  Bounds = true
      YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss
      1980-01-01 12:00:00  1980-01-02 12:00:00  1980-01-03 12:00:00  1980-01-04 12:00:00
      1980-01-05 12:00:00  1980-01-06 12:00:00  1980-01-07 12:00:00  1980-01-08 12:00:00
    (中略)
      1980-12-22 12:00:00  1980-12-23 12:00:00  1980-12-24 12:00:00  1980-12-25 12:00:00
      1980-12-26 12:00:00  1980-12-27 12:00:00  1980-12-28 12:00:00  1980-12-29 12:00:00
      1980-12-30 12:00:00  1980-12-31 12:00:00
    

2.1.3 Compare two datasets: diff, diffn

2.1.6 Dataset description: pardes, griddes, zaxisdes, vct

 抜き出し(変数,鉛直層の指定)

2.3.1 Select fields: select, delete

  • select の例: select,level=レベルの値
    $ cdo select,level=85000  ifile01.nc ifile02.nc ifile03.nc  ofile.nc
    
    時間方向に関して3つのファイルに分割記録してある 850 hPa面データだけを切り出し,1つのファイルに出力.
    気圧面は Pa で記録されているから,850ではなく 85000 で指定する.
    入力ファイルの並び順を間違えないように注意.
    

 抜き出し(時間)

2.3.3 Select timesteps: seltimestep, seltime, selhour, selday, selmon, selyear, selseas, seldate, selsmon

  • seltimestep, selmon
    $ cdo seltimestep,9 input.nc  (9番目のデータのみ切り出し)
    $ cdo selmon,3 input.nc  (3月だけを切り出し)
    

 格子変換・内挿

2.12.1 SCRIP grid interpolation: remapbil, remapbic, remapdis, remapnn, remapcon, remapcon2, remaplaf

  • remapbil
    $ cdo remapbil,r288x145  input_Gauss.nc  ofile_1.25deg.nc
    (入力ファイル (全球,ガウス格子) を 1.25度間隔の等緯度経度格子に変換)
    
    $ cdo remapbil,n63  input_linear.nc  ofile_T63.nc
    (入力ファイル (全球,等緯度経度格子) を全波数63のガウス格子に変換)
    

CMIP5 のデータは,モデルごとに格子間隔が異なる. 再解析データなどと比較したり,モデル間の相互比較をする際には, 格子間隔が同一である方がラク. コマンド一発で変換できるので,とても便利.

 GrADS ctl 作成

2.15.1 GRADSDES - GrADS data descriptor file: gradsdes

  • gradsdes
$ cdo gradsdes  input.nc
(入力ファイルの情報から GrADS ctl ファイルを自動作成)

NetCDF だけでなく,GRIB1/2 でもできるらしい(GRIB1,GRIB2での動作は未確認; 2015-06-20).

CMIP5のオリジナルファイルでは正常作成された. しかし,remapbil や select などで抽出・格子変換したファイルでは, ctl ファイルに記述ミスが発生した(cdo v1.6.4 で確認; 2015-06-20). 原因特定できてないが,cdo の最新版では解決されている可能性あり(記述時点で v1.6.9). See also GrADS memo (hysk).

CDO サンプルスクリプト

 1気圧面を抜き出し,等緯度経度格子に変換

daily mean u の3次元データファイル(鉛直8層)から850hPa面だけを抽出.
さらに格子系をガウス格子(320x160)から等緯度経度(2.5度,144x73)に変換.

(最も単純な実行例)
$ cdo reampbil,r144x73 -select,level=85000  input.nc  output.nc

下記は,膨大なファイル数を処理するための雛形ファイル. これに for/while などを追加してファイルパスやファイル名を入れ替えていくように 『肉付け』すればよい.

#!/bin/bash
# CDO sample: Choose 850-hPa level, 
# then grid system converted into liner lat/lon grid with 2.5 degree interval.
#
#   1st release: 2015-06-20
#
#   M. Hayasaki (NIES, Japan)

### Input data CMIP5, historical, MRI-CGCM3, RIP=r1i1p1, daily mean U (ua), Year=1980
# Original input file has eight pressure levels (1000, 850, 700, 500, 250, 100, 50, and 10 hPa).
var_cmip=ua
temporal=day
model=MRI-CGCM3
exp_name=historical
rip=r1i1p1
period="19800101-19801231"

# Select a target pressure level
level=850   # unit: hPa
  lev_4digit=`printf %4.4d ${level}`
  level_SI=`echo "${level} * 100" | bc -q`  # unit: Pa

# Set a grid system (2.5 degree interval, 73 grids in latitude, 144 grids in logitude
  nx=144
  ny=73
  remapbil_opt="r${nx}x${ny}"

var_lev=${var_cmip}${lev_4digit}    # ex. "ua0850"

ipath=$HOME/dw/tmp/cdo_test
ifile=${ipath}/${var_cmip}_${temporal}_${model}_${exp_name}_${rip}_${period}.nc
# ex. ua_day_MRI-CGCM3_historical_r1i1p1_19800101-19801231.nc

opath=$HOME/dw/tmp/cdo_test
ofile01=${opath}/hoge_${var_lev}.nc
ofile02=${opath}/hoge_${var_lev}_parallel.nc    # test output for using openmp threads

## Case 1: Single thread
time cdo      remapbil,${remapbil_opt} -select,level=${level_SI} $ifile   $ofile01
## Case 2: Use 4 threads
#time cdo -P 4 remapbil,${remapbil_opt} -select,level=${level_SI} $ifile   $ofile02

exit

# [History, personal memo]
# program path = $HOME/cmip/regrid/310_select_1lev_regrid_day_3D.sh

CMIP5のモデル群のように,解像度が異なる多数のモデルの計算結果の解析(統計操作,相互比較など)をするには, 解析の初期段階で格子系を統一するのがラク. cdo では,複数の操作をパイプを通して連続処理できる (See 1.2.4 Operator chaining; リンク切れしてたらCDO User's Guide (html) よりセクション番号や見出しで見つけよ).

中間ファイルを作成する必要がなくなるため,ディスク領域の節約になる. また,実行するコマンドの行数も減るので,スクリプト全体の見通しが良くなり,メンテナンス性が向上する.

並列化(複数スレッドでの処理,-P オプション後にスレッド数を指定)も可能だが, オペレーションの種類によってはあまり恩恵がない. 実行回数が莫大であるならば,本実行の前に少ない実行数でテストし,並列化の恩恵がありそうか否かを確認しておくべし.

CDO 関連リンク

  • CDO 本家サイト
    • source のダウンロードは,右上部の Files より.Windows 用のバイナリパッケージも利用可能らしいが,私は未検証.
    • CDO User's Guide (html): 各種オペレータの解説がある. どんなことが出来るのか,全体を見渡すときにも有効. オペレータの綴りを見れば,およそどんなことが出来るのかは推測できる(checked 2014-02-18).

変更履歴

更新日 内容
2015-06-20 リンク切れの修正, gradsdes の追加, サンプルスクリプト(気圧面抜き出し&linear latlon grid 変換)の追加.
2014-06-25 記述開始