シェルスクリプトとawkによるデータ解析

2006/01/10 筑波大学農林工学系 西田顕郎

はじめに

テキスト形式の数値データの処理・解析には、awk+シェルスクリプトが最強である。誇張ではない。これまでエクセルやC言語、FORTRANなどしか知らなかった人には、ぜひawk+シェルスクリプトの世界を体験して欲しい。気象データなど、大量のテキスト形式数値データをいじる人には、人生の転換点と言えるほどの衝撃が待っている。

なお、awk+シェルスクリプトの使用環境は、当然ながらUNIXシェルである。Windowsユーザーは、Cygwinで擬似的にUNIXシェル環境を構築してください。

参考になるページ

こことか
こことか


awkとは?

awkは、テキストデータの1行づつに対して処理をする、プログラミング言語(スクリプト言語)である。コンパイルは必要ない。UNIXのコマンドラインの出力に、パイプでつないで実行できる。

awkは列を抜き出す!

例1: ファイルのリスト(ls -lの出力)をawkで加工してみよう。たとえばls -lの出力が下記のようであったとする:
$ ls -l
-rw-r--r--    1 nishida  nishida       802 Dec 27 15:08 IJ_Prefs.txt
-rw-r--r--    1 nishida  nishida    238197 Jan  6 15:22 ame_master.pdf
-rw-r--r--    1 nishida  nishida    816064 Jan  6 15:47 wg3ts.pdf
-rw-r--r--    1 nishida  nishida   2527691 Jan  6 15:21 tebiki.pdf
-rw-r--r--    1 nishida  nishida   5019689 Jan  6 15:36 haichi.pdf
このとき、ファイルサイズとファイル名のみを抜き出して表示してみよう。
$ ls -l | awk '{print $5,$NF}'
802 IJ_Prefs.txt
238197 ame_master.pdf
816064 wg3ts.pdf
2527691 tebiki.pdf
5019689 haichi.pdf
ここで、$5は「5列目のデータ」の意味。NFは、各行のフィールド数(列の数)。$NFは、「NF列目のデータ」、つまり、その行の最後の列のデータ、ということ。

awkは計算もする!

例2: 例1の出力で、ファイルサイズを1024で割る、つまりキロバイト単位で表示する。例1で打ったコマンドを、キーボードの「上向き矢印」のキーで読み出して、ちょっとだけ変更すればよい(/1024を書き加える):
$ ls -l | awk '{print $5/1024,$NF}' 0.783203 IJ_Prefs.txt 232.614 ame_master.pdf 796.938 wg3ts.pdf 2468.45 tebiki.pdf 4902.04 haichi.pdf

awkは条件分岐もできる!

例3: 例1の出力で、1000バイト以上のファイルについてのみ処理する:
$ ls -l | awk '1000<=$5{print $5,$NF}'
238197 ame_master.pdf
816064 wg3ts.pdf
2527691 tebiki.pdf
5019689 haichi.pdf
これらの例を見れば、awkの文法は、なんとなく直観でわかるだろう。実際、awkは極めて直観的な言語である。文法はCに似ているが、Cよりもはるかにフレキシブルで直観的だ。ifもforも使える。変数も配列も使えるが、それらの初期宣言は不要だ。これだけわかれば、もうawkのポテンシャルを発揮できる。以下、実際のデータをいじってみよう。

データ列のならべかえ

以下のデータは、5つの生態観測サイトの、経度・緯度・名前を順に並べたものである:

137.42 36.14 Takayama
141.51 42.73 Tomakomai
138.8  35.45 Fuji-Yoshida
130.71 33.13 Kahoku
140.10 36.10 TERC
これが、site.txtという名前のテキストファイルになっているとしよう。

これを、あるアプリケーションに渡すために、名前・緯度・経度の順に並べ替える必要があるなら、このようにする:

$ cat site.txt | awk '{print $3,$2,$1}'
Takayama 36.14 137.42
Tomakomai 42.73 141.51
Fuji-Yoshida 35.45 138.8
Kahoku 33.13 130.71
TERC 36.10 140.10
極めて簡単である。「並べ替えくらいエクセルでもカットペーストでできる」と思う人、あなたは甘い。エクセルで扱うことができるのは、65536行までのデータである。しかし、awkなら、データが数万行にわたっても問題なく走る。もっと大切なことは、後述するが、awkの処理はシェルスクリプトに書き込むことができる、ということである。

この出力を、アルファベット順に並べたいとする。その場合、UNIXのsortコマンドと組み合わせて、

	
$ cat site.txt | awk '{print $3,$2,$1}' | sort
Fuji-Yoshida 35.45 138.8
Kahoku 33.13 130.71
TERC 36.10 140.10
Takayama 36.14 137.42
Tomakomai 42.73 141.51

大規模データの処理

ここでは、実際にフィールドで観測されたデータを解析してみよう。サンプルは、有名な森林生態観測サイトである、岐阜県高山フラックスサイトで測られた日射(光合成有効放射;Photosynthetic Active Radiation; PAR)のデータである。このリンクを右クリックなどで保存し(ファイル名はsample.txt)、中を覗いてみよう:
$ cat sample.txt
2005,170,600,1.058,.044,-.033
2005,170,605,.991,.06,.067
2005,170,610,1.022,.033,.087
2005,170,615,1.062,.02,-.017
2005,170,620,1.313,.044,.013
2005,170,625,1.186,.03,-.023
......
各行に、各年,日付,時刻,センサ1,センサ2,センサ3の順にデータが入っている。日付は、1月1日を1とし、12月31日を365とする、年間通算の日付け(Day of Year またはJulian Dayともいう)、時刻は、1または2桁で時、2桁で分をあらわす(たとえば605は午前6時05分)。各列はコンマで区切られている。

センサは1から3まであり、センサ1は森の樹冠の上に突き出たタワーの上で、空からやってくる光を観測している。センサ2, 3は、林床に設置されており、森の樹冠を透過した光を、それぞれ別の場所で観測している。センサ1〜センサ3の出力は、それぞれセンサの電圧(mV)であらわされており、これを実際の物理的な値(この場合はμmol/m2/sec)に変換するには、それぞれのセンサに固有の係数をかけてやる必要がある。

まず、センサの出力を物理値に変換しよう。別途の較正作業の結果から、センサ1〜センサ3の変換係数は、以下のようにわかっている:

センサ1  センサ2  センサ3
293.73   306.82   290.84
これらの値を、各センサの出力(mV)にかけてやれば、センサの電圧出力は物理値に変換される。それにはawkをこのように使う:
$ cat sample.txt | awk 'BEGIN{FS=","}{print $1,$2,$3, $4*293.73, $5*306.82, $6*290.84}'
2005 170 600 310.766 13.5001 -9.59772
2005 170 605 291.086 18.4092 19.4863
2005 170 610 300.192 10.1251 25.3031
2005 170 615 311.941 6.1364 -4.94428
2005 170 620 385.667 13.5001 3.78092
......
ここで、BEGIN{FS=","}というのは、「最初にまず、列のさかいめ(field separator,略してFS)はコンマで表す」という意味。これを宣言することで、各行の列をコンマで区切ることが指示される(ちなみにデフォルトではFSはスペースかタブである)。その次の{}内の意味は自明だろう。

次に、時刻を"605"のような時分表示ではなく、1日の通算の分で表示しよう。というのも、時分表示では、繰り上がりにギャップが出来るからである(例えば659の次が660でなくて700になる):

$ cat sample.txt | awk 'BEGIN{FS=","}{min=$3%100; hour=($3-min)/100; print $1,$2,min+hour*60, $4*293.73, $5*306.82, $6*290.84}'
2005 170 360 310.766 13.5001 -9.59772
2005 170 365 291.086 18.4092 19.4863
2005 170 370 300.192 10.1251 25.3031
2005 170 375 311.941 6.1364 -4.94428
2005 170 380 385.667 13.5001 3.78092
さらに、時刻を1日単位で日付けにあわせてしまおう。これは、複数の日にまたがるデータを統一的に扱うためである:
$ cat sample.txt | awk 'BEGIN{FS=","}{min=$3%100; hour=($3-min)/100; print $1,$2+(min+hour*60)/1440, $4*293.73, $5*306.82, $6*290.84}'
2005 170.25 310.766 13.5001 -9.59772
2005 170.253 291.086 18.4092 19.4863
2005 170.257 300.192 10.1251 25.3031
2005 170.26 311.941 6.1364 -4.94428
2005 170.264 385.667 13.5001 3.78092
この出力を、"conv.txt"という新規ファイルに入れよう(リダイレクト):
$ cat sample.txt | awk 'BEGIN{FS=","}{min=$3%100; hour=($3-min)/100; print $1,$2+(min+hour*60)/1440, $4*293.73, $5*306.82, $6*290.84}' > conv.txt
このconv.txtをもとに、gnuplotという、コマンドライン対応型のグラフソフト(言いかえれば、グラフを描くためのスクリプト言語)で、センサ1,2,3の観測値の変動のグラフを描こう:
$ gnuplot
gnuplot> plot "conv.txt" u 2:3 w l, "conv.txt" u 2:4 w l, "conv.txt" u 2:5 w l 

gnuplotから抜けるには、

gnuplot> exit
さて、この一連の処理は、viなどのエディタを使って、以下のようなreadme.shというファイル名のシェルスクリプトにまとめてしまうことができる:
#!/bin/sh
cat sample.txt | \
awk 'BEGIN{FS=","}
     {min=$3%60; hour=($3-min)/100; print $1,$2+(min+hour*60)/1440, $4*293.73, $5*306.82, $6*290.84}' > conv.txt

gnuplot <<EOF
set terminal png
set output "Takayama_PAR.png"
set title "Takayama PAR"
set yrange [0:2500]
set ylabel "PAR [umol/m2/sec]"
set xlabel "DOY (2005)"
plot "conv.txt" u 2:3 title "above canopy" w l, "conv.txt" u 2:4 title "under canopy" w l, "conv.txt" u 2:5 title "under canopy" w l
EOF
これをviなどのエディタで作ってreadme.shという名前で保存して、

$ chmod +x readme.sh
$ ./readme.sh
とすれば、一連の処理がぜんぶ自動的に実行され、以下のようなグラフがTakayama_PAR.pngというファイル名でできあがる:

このように処理をシェルスクリプトにすることは、以下のようなメリットがある:

1. 1次データから最終結果(グラフ)までの処理工程を全て完全に記録でき、必要ならば簡単に再現できる。
2. 同様のデータについて、同様の処理をするときに使いまわすことができる。
3. 処理の一部を変更して全てをやり直すのも容易である。たとえばこの例では、センサの出力を物理値に変換する係数を変更する、などの場合である。
4. 処理のノウハウを公開したり、他の研究者と共有することが容易である。

ここで述べた1の事項、つまり1次データから最終結果(グラフ)までの処理工程を全て完全に記録するということは、一般にはあまり意識されていないが、実は科学研究では極めて重要なことである。科学者・技術者の倫理として、なんらかの研究・解析をやった結果のデータを世間に発表する前に、その正しさを十分にチェックしなければならないし、発表したあとも、常にそのデータ処理の細部に関して、だれかから問い合わせがあったときに対応できなくてはならない。そのためには、最終的に発表するデータが、どのようなデータから、どのような工程を経て加工されたか、たどれなければならない。これを、「処理のトレーサビリティ」と呼ぼう。処理のトレーサビリティを確保するためには、データに施した全ての処理を記録する必要がある。シェルスクリプトはそれを実現するのだ。そのようなことを、シェルスクリプトを使わずに、たとえば表計算ソフトなどのグラフィカルインターフェースに基づく処理で行うのは大変である。


統計量の計算

上の例で、こんどはセンサー1の出力について、各日の日積算値を求めてみよう:

$ cat conv.txt | awk 'BEGIN{D=170;a=0;t=170.0}D+1<=$2{print D, a; a=0; D++}{a=a+$3*($2-t)*60*60*24/1000000; t=$2}'
170 46.912
171 44.421
172 47.185
173 27.6938
174 23.9657
175 43.4051
176 55.7856
177 38.7133
178 4.91674
179 4.5035
単位は、mol/dayである。この考え方を使えば、時間ごとの積算や平均値も簡単に計算できる。

もっと簡単な例として、日射(PAR)の、森林樹冠透過率を求めてみよう。センサ2やセンサ3の出力を、センサ1の出力で割れば良い:

$ cat conv.txt | awk '{print $1,$2,$4/$3,$5/$3}'
2005 170.25 0.0434414 -0.0308841
2005 170.253 0.0632432 0.0669434
2005 170.257 0.0337287 0.0842897
2005 170.26 0.0196717 -0.01585
...............
2005 170.735 -0.0335411 -0.317942
2005 170.753 0 0
2005 170.757 0 0
awk: コマンドライン:1: (FILENAME=- FNR=147) 致命的: ゼロでの除算を行いました。
上のようなエラーが出てしまった。これは、センサ1の出力がゼロのとき、ゼロが分母になるような割り算をすることになってしまうからである。これを回避するには、以下のようにすればよい:
$ cat conv.txt | awk '$3>0{print $1,$2,$4/$3,$5/$3}'
2005 170.25 0.0434414 -0.0308841
2005 170.253 0.0632432 0.0669434
2005 170.257 0.0337287 0.0842897
2005 170.26 0.0196717 -0.01585
................
2005 180.725 0.113284 0.0655461
2005 180.728 0.0378628 0.0844489
2005 180.732 0 0.106089
2005 180.735 0.175067 0.0912721
これを、前と同様のグラフにまで描くスクリプトにすると、以下のようになる:
#!/bin/sh
cat conv.txt | awk '$3>0{print $1,$2,$4/$3,$5/$3}' > trans.txt
gnuplot << EOF
set terminal png
set output "Takayama_PAR_trans.png"
set title "Takayama PAR transmittance"
set yrange [0:1]
set ylabel "PAR transmittance"
set xlabel "DOY (2005)"
plot "trans.txt" u 2:3 title "sensor 2" w l, "trans.txt" u 2:4 title "sensor 3" w l
EOF

ファイルを並べてくっつけるには?

たとえば、

001 2
002 3
003 1
004 5
005 6

というデータを含む、data1.txtというファイルと、

1
5
4
4
5

というデータを含む、data2.txtというファイルのデータを、横に並べて

001 2 1
002 3 5
003 1 4
004 5 4
005 6 5

というふうなファイルにしたい、ということがよくある。そういうときには、

$ paste data1.txt data2.txt

とすればよい。標準出力に書き出されるから、ファイルにリダイレクトしてもいいし、
パイプでawkなどに流し込んでもよい。3つ以上のファイルのデータをいっぺんに
並べることもできる:

$ paste data1.txt data2.txt data3.txt

ファイル間でデータの並びに不整合があるときは...joinコマンドが便利! 双方のファイルに
キーとなる列(時刻とか)があると、それを合致させるように行どうしをくっつける。


awkメモ
外から引数を与えるには:

awk '......' x=1 y=$y filename
とか、
cat filename | awk '......' x=1 y=$y

注意: 外から与える引数は、BEGIN{}の中には反映されない。

逆三角関数: atan2(y,x) ... y/xのアークタンジェント