awkによる2次元格子点データの処理

条件によるデータの抽出

Unixなど大部分のOSの上では、awkの短いプログラムをコマンドラインから

  awk 'プログラム' 入力ファイル
のような形で実行することができる。

awkのプログラムは

  パタン { 実行文 }
のような構造からできているが、パタンだけを書くと、実行文として printが仮定され、入力ファイルの行と同じものが標準出力に書き出される。

したがって、あるパタンにあてはまる行だけを取り出すためには、

  awk 'パタン' 入力ファイル > 出力ファイル
のようにすればよい。

たとえば、in.xyvの中に、 データが欠けている地点が、V = -9000として書かれているとし、 その行を除いて有効なデータだけをout.xyvとすることは、

  awk '$3+0 != -9000' in.xyv > out.xyv
でできる。

【!=は「等しくない」である。 なお、==が「等しい」、単独の=が代入記号である。 「等しい」のつもりで単独の=を使うと、 一見正常に動作するが結果が期待と違うことがあるので注意。】

【$3のあとにつけた「+0」は、 -9000と-9000.0と-9000.000などを区別したくないので、 必ず「数値として扱う」ようにさせるための処理である。】

逆に、データの欠けているところ(といっても、実際は-9000の値が入れてある ところにすぎず、in.xyvに出てこない位置は含まれない)を missing.xyvに書き出すには、

  awk '$3+0 == -9000' in.xyv > missing.xyv
とする。

また、in.xyvのうちで、 値が0以上30以下のところをnormal.xyv、 それ以外のところをabnormal.xyvに出したいとすれば、

  awk '($3+0 >= 0) && ($3+0 <= 30)' in.xyv > normal.xyv
  awk '($3+0 < 0) || ($3+0 > 30)' in.xyv > abnormal.xyv
とすればよい。(&&と||は、論理学でいうandとorを表わす。 単独の&や|は別の意味になるので注意。) ただし、この手順では、in.xyvを2回読むことになる。 awkプログラムを組めば、 1回読んで両方に分けて書き出すこともできる。 また、awkプログラムの中でif / else構文を使えば、 判断が1回ですむ。 たとえば、例文として用意したcheck.awkを、
awk -f check.awk in.xyv
として実行すればよい(この場合、標準出力には何も出てこないはず)。

判断を2回する方式だと、 うっかりしてnormalとabnormalの両方に書き出す場合や どちらにも書き出さない場合を作ってしまうことがあるので、 そのようなまちがいを少なくするためには、判断は1回にまとめたほうがよい。

$1, $2の値を使って、緯度・経度の範囲も抽出の条件にすることができる。


1つのファイルを読んで処理

awkの標準的な使い方は、

  awk -f プログラムファイル 入力ファイル
である。ここで入力ファイルは複数ならべることもできるが、その場合、まずひとつの ファイルを最後まで読み、続いて次のファイルに進む。 ここでは入力ファイルが1個の場合に限って述べる。 結果をファイルに出したい場合は、 プログラム中でファイル名を指定する方法と、 プログラムでは標準出力に書くことにしておいて Unixのシェルのredirection機能を使う方法がある。 後者の方法をとる場合、
  awk -f プログラムファイル 入力ファイル > 出力ファイル
のような使い方をする。


[換算]

precip.xyvには、年降水量(単位: mm)がはいっており、 これを、これだけの水を蒸発させるのに必要なエネルギーフラックス密度 (単位: W / m2)に換算したいとする。 ただし、1年の日数は365.25日とみなす。 また、水の質量あたりの蒸発の潜熱は定数とみなし、2.5×106 J / kg とする。 液体の水の密度も定数 1 g / cm3とみなすので、 1 mmをそのまま1 kg / m2とみてよい。 入力ファイルでは、 値の欠けているところには -9000がはいっており、これはそのまま書き出したい とする。

この場合の処理プログラムは convert.awkのようになる。 各行ごとに値が変わらない数値の設定は、 入力ファイルを読み始める前に実行されるBEGIN部にまとめてみた。 なお、結果を縦にそろうように出力したい場合は、printの代わりに printfを使って、たとえば次のように書式指定する (BEGIN部は省略した)。 awkのprintfの書式指定は 基本的にC言語のprintfと同様だが、 変数の型の区別を厳密にしなくてよい。

(convert.awkのmain部分の書きかえ)

#main
{
  v = $3+0
  if(v == vmiss) { printf "%10.3f%10.3f%10.3f\n", $1, $2, $3 }
  else           { printf "%10.3f%10.3f%10.3f\n", $1, $2, factor*v }
}


[緯線上の東西平均]

気候データを扱う場合、同じ緯度についてすべての経度での値を平均して みることが多い。 緯度経度格子の場合、緯度が違うと格子点の代表する面積が違うが、 同じ緯度では同じなので、平均する際に格子点ごとに重みを変える必要はない。

たとえば、temp.xyvの東西全経度平均を作るには、 例文zonalmean.awkを使って、

  awk -f zonalmean.awk temp.xyv | sort -n > temp.zonalmean
とすれば、出力ファイルtemp.zonalmeanに、 第1列に緯度、第2列に平均の結果、第3列に平均に使われた格子点の個数が 書き出される。 ただし、入力ファイルでは、値が欠けているところは-9000で示されている とする。出力ファイルでは値が定義されない行は書き出されない。 このawkプログラムでは、出力結果が緯度の値の順にならぶ保証がないので、 結果をsortで処理して順序を整えた。

【このプログラムzonalmean.awkでは、 sumとcountの2つの配列の各要素の値が、 たしこみを始める前に0になっていることを仮定している。 awkの場合、これは自動的に満たされているが、多くの言語では、 配列を宣言し、 内容の各要素に0を入れる初期設定をしておく必要がある。 初期設定部分をawkで書いてみると、たとえば次のようになる。 ただし、これは、格子が2.5度間隔(極も格子)の場合である。

BEGIN {
  for(lat = -90; lat <= 90; lat += 2.5)
  { sum[lat] = 0
    count[lat] = 0
  }
} 


[集計とヒストグラム]

整数値で表わせる定性的データについて (awkで扱う場合、実際の値は整数である必要はないが)、 どの値がいくつの格子点に出現するかを数えるプログラム count.awkを書いてみる。 やはり、awkでは変数の初期値として0が仮定されることを使っている。 これも、結果をsort -nで処理したほうがよいだろう。

緯度経度格子の場合、 それぞれの格子点が代表する面積が緯度によって違う。 個数でなく面積を集計したい場合もある。 プログラムarea.awkは、 緯度経度2.5度格子、ただし極そのものには値がなく、 極から格子間隔の半分だけずれたところに値がある場合のものである。 BEGIN部で、各格子点の面積を緯度ごとにあらかじめ計算している。 この部分は、なるべく謎の定数を使わないように、 1回しか使わない値にも変数名をつけて記述した。 出力する面積の単位は、109 m2つまり 1000 km2とした。

変数が連続量(実数)の場合、それを区間に分けて度数分布を集計する。 区間の幅vintで集計する例を histogram.awkとして示す。 ただし、集計や出力の段階で両端の区間を特別扱いしないですむように、 -big未満やbig以上の値は データ中に絶対に出てこないと仮定した。

区間ごとに格子点数ではなく面積を集計することも、 すでに述べた考え方を組み合わせればできる。


複数のファイルを読んで処理: 格子点のならびが同じ場合

複数のファイルを読んで計算する場合を考える。 もし、入力ファイルの何行めがどこの格子点であるかが、 どのファイルでも共通ならば、 データを配列にためずに、 入力ファイルをそれぞれ1行ずつ並列に読んで処理することができる。 ただし、awkでコマンド行に入力ファイルとして書くものは、 ひとつ読み終わってから次を読むという順次処理しかできないので、 getline関数で入力ファイル名を指定して読む方法をとり、 コマンドラインから変数の値としてファイル名をプログラムに渡すことにする。


[3か月合計]

pjan.xyv, pfeb.xyv, pmar.xyvに それぞれ1月、2月、3月の降水量がはいっているとする。 これを読んで、3か月合計降水量を計算して書き出すプログラム add3.awk を作ってみる。

これを使うには、次のようなシェルスクリプトadd3.sh

  awk -f add3.awk -v f1=$1 -v f2=$2 $3 > $4
を作っておき
  add3.sh pjan.xyv pfeb.xyv pmar.xyv pjfm.xyv
で実行するとよい。


[年平均]

同様に、 毎月のデータvjan.xyv, ..., vdec.xyvから 年平均を計算するプログラム annmean.awk を作ってみる(ただしまだ動作確認していない)。 月ごとの日数の違いも考慮し、2月は28.25日として扱う。 全部の月の値があるとき(便宜上、合計日数360日以上という条件で表わす)だけ 計算結果を、そうでなければデータの欠けた印の値を書き出すことにする。


複数のファイルを読んで処理: 格子点のならびが同じと限らない場合

データ相互に座標の位置が同じであってもその出てくる順序が同じでない場合は、 あらかじめ(BEGIN部でgetlineにより) 一方のデータを配列に読みこんでおき、 他方のデータを1行ずつ読み、その座標値を配列のものとつき合わせて 処理していく方法をとる。

この場合、次のような注意が必要である。


[2つのデータの値の差をとる]

同じ変数のデータセットで、 ヨーロッパで作られたveuro.xyvと アメリカで作られたvamer.xyvがあるとする。 両者の値の差を書き出すプログラム diff.awk を書いてみる。

使い方は

  awk -f diff.awk -v f1=veuro.xyv vamer.xyv > vdiff.xyv
のようにする(vamer - veuroが書き出される)。


[海陸別の東西平均]

前の緯度別東西平均の例を、 海陸別の東西平均を計算するように修正してみたのが landseazonal.awkである。 landsea.xyvには、陸ならば1、海ならば0がはいっているとする。 各緯度ごとに、陸の平均と海の平均を書き出すようにする。 ただし、前の例と違って、緯度方向の格子は2.5度間隔(極も格子点)に限る ことにする。


[領域(面積)平均]

前の面積集計の例文を修正し、 type.xyvにはいっている土地被覆の種類ごとに、 別の.xyvファイルにはいっている数量の平均を計算して書き出す プログラム typearea.awk を作ってみる。 格子間隔は2.5度と仮定したが、 dlon, dlatさえ変えれば他の間隔に対応できる。 結果はsort -nで土地被覆コード番号順にして 見たほうがよいだろう。


[クロス集計と2元ヒストグラム・散布図]

ふたつの定性的変数 v1v2があるとき、 変数の組み合わせ(v1,v2)ごとに格子点の数を集計する プログラム cross.awk を作ってみる。 使い方は、

  awk -f cross.awk -v f1=v1.xyv v2.xyv > cross.txt
のような形になる。 ただし、このままでは、縦横とも、変数の順序はたまたまawk処理系が 見つけた順になり、番号の大小とも個数の大小とも一致しない。 実用的には、書き出す作業にはいる前に、 利用目的に合わせて行や列のならべかえをしたほうがよい。

2つの定量的データがある場合は、2元ヒストグラムを作ることができる。 1元ヒストグラムの例文を改造することによって プログラムしてみたのが histv1v2.awkである。 使い方は、

  awk -f histv1v2.awk -v f1=v1.xyv v2.xyv > histv1v2.txt
のようにする。 1と2が違うだけで基本的に同じことのくりかえしが多いので、 awkのユーザー定義関数の機能を使うと、もっとすっきり書けると思う。

2つの定量的データの組み合わせについては、散布図も重要な方法である。 散布図を作る道具として、gnuplotを使うとすれば、 まず2つのファイルを読んで 「経度、緯度、変数1、変数2」の4列からなる形に編集し、 その3列めと4列めを横軸・縦軸にとって作図すればよい。 この編集作業をするプログラムは、 4cols.awk のように簡単な形に書ける。 使い方は、

  awk -f 4cols.awk -v f1=v1.xyv v2.xyv > v1v2.txt
そして、gnuplotの中で
  plot "v1v2.txt" using 3:4
のようにする。


[配列の値が定義された要素だけを処理する]

1次元配列aの要素a[i]が定義されている場合にだけ 何かの処理をするには、

  for ( i in a ) { ... }
のような構文を使えばよい。ただし、 要素の順序は処理系まかせになり、iの大小とは関係ない。

2次元配列aの要素a[lon, lat]については、 awkのバージョンによっては、

  for ( (lon,lat) in a ) { ... }
の構文が使えるが、これがエラーになって動かない処理系もある。 その場合は、次のような構文を使えばよい。
  for ( str in a )
  { 
    split(str, tmparray, SUBSEP)
    lon = tmparray[1]+0
    lat = tmparray[2]+0
    ... a[lon,lat] に関する処理 ...
  }


2000-08-23; 2002-04-08; 2004-06-22
増田 耕一 (MASUDA Kooiti)

[目次へ]