いちからはじめる GMT その74 - blockmean

2013年 2月 1日 掲載

今回は，簡単な地理座標上のデータ処理を行うコマンド，blockmean を紹介する．

GMT にはデータの統計処理を行うツールがいくつかふくまれており，blockmean もそのひとつである． blockmean は，地理座標上の指定された大きさに区切られた区画内のデータの平均値を計算する． 波長の短いノイズを除去したいときなどに使うと便利である．

では，さっそく試してみよう． まずはサンプルデータ作りである． 西経 6度〜東経 6度，南緯 6度〜北緯 6度の範囲のランダムな座標値を 1000個作り，

$ awk 'BEGIN{for(i=0;i<1000;i++)print rand()*12-6,rand()*12-6}' > xy.txt

意味のあるデータ（以下，シグナルと呼ぶ）として，東西・南北方向の波長が 24度，波高が 1 の正弦波形（といっていいのかな？）を作った．

$ awk 'BEGIN{pi=atan2(0,-1)}{print $1,$2,cos(pi*$1/12)*cos(pi*$2/12)}' xy.txt > signal.txt

できあがったシグナル signal.txt を図にしてみる． カラーパレットには，赤〜青の polar を使った． 正の値が赤，負の値が青で，値が 0 だと白になる．

$ makecpt -Cpolar -D -Z > p.cpt
$ psxy signal.txt -Jm1 -R-6/6/-6/6 -Sc0.2 -Cp.cpt -K -Ba2 > signal.eps
$ psxy signal.txt -J -R -Sc0.2 -O >> signal.eps

シグナル

$ surface signal.txt -Gsignal.grd -R-6/6/-6/6 -I0.125
$ grdimage signal.grd -Jm1 -Cp.cpt -Ba2 -K > signal2.eps
$ grdcontour signal.grd -J -C0.1 -A0.2 -O >> signal2.eps

シグナル

ノイズとしては，波長 2度，波高 1 の正弦波形を作った．

$ awk 'BEGIN{pi=atan2(0,-1)}{print $1,$2,sin(pi*$1)*sin(pi*$2)}' xy.txt > noise.txt

ノイズ noise.txt の図は，こんな感じ．

$ psxy noise.txt -Jm1 -R-6/6/-6/6 -Sc0.2 -Cp.cpt -K -Ba2 > noise.eps
$ psxy noise.txt -J -R -Sc0.2 -O >> noise.eps

ノイズ

$ surface noise.txt -Gnoise.grd -R-6/6/-6/6 -I0.125
$ grdimage noise.grd -Jm1 -Cp.cpt -Ba2 > noise2.eps

ノイズ

さて，シグナルにノイズを加えて観測データを作成しよう． ノイズの大きさはシグナルの半分にしてみた．

$ paste signal.txt noise.txt | awk '{print $1,$2,$3+$6*0.5}' > data.txt

できあがった観測データ data.txt を図にするとこんな感じ．

$ psxy data.txt -Jm1 -R-6/6/-6/6 -Sc0.2 -Cp.cpt -K -Ba2 > data.eps
$ psxy data.txt -J -R -Sc0.2 -O >> data.eps

観測データ

$ surface data.txt -Gdata.grd -R-6/6/-6/6 -I0.125
$ grdimage data.grd -Jm1 -Cp.cpt -Ba2 -K > data2.eps
$ grdcontour data.grd -J -C0.1 -O >> data2.eps

観測データ

これで，ようやく準備ができた． やっと blockmean の出番である． ノイズを除去するためには，ノイズより波長の大きい範囲の値を平均すればよい． 今回は 3度×3度に区切った区画内の値を平均する．

$ blockmean data.txt -I3 -R-6/6/-6/6 > mean.txt

できあがった平均データ mean.txt を見てみると，

$ wc -l mean.txt
      25 mean.txt
$ head mean.txt
-5.32571722222	5.33426	0.0204766
-3.18644454545	5.25311181818	0.130819613636
-0.0782298807692	5.19502230769	0.206328084615
2.96212352941	5.22756764706	0.199016077941
5.43345642857	5.24136	0.121436121429
-5.25707521739	2.9417	0.197093694261
-2.97982703704	2.93283425926	0.491710814815
-0.194725390371	3.08259414286	0.652666714286
2.9566595	2.832195	0.472925925
5.264784	2.9765244	0.1445304028

あれ？なんか経緯度の値がごちゃごちゃしているな． どうなっているのか，図を描いてチェックしてみよう．

$ psxy mean.txt -Jm1 -R-6/6/-6/6 -Sc0.2 -Cp.cpt -K -Ba2g2 > mean.eps
$ psxy mean.txt -J -R -Sc0.2 -O >> mean.eps

平均データ

blockmean はデータの値を平均するだけでなく，平均に使った値の座標の経緯度も平均しているようである． なるほど，確かに理にかなっている．

最後に，平均データを図化してみる．

$ surface mean.txt -Gmean.grd -R-6/6/-6/6 -I0.125
$ grdimage mean.grd -Jm1 -Cp.cpt -Ba2 -K > mean2.eps
$ grdcontour mean.grd -J -C0.1 -A0.2 -O >> mean2.eps

平均データ

だいたい，元のシグナルに近いデータになったんじゃないかな．

ちなみに，単純に 3度間隔の格子点付近のデータを抽出するだけだと（今回は AWKスクリプト grid.awk を作って使った），

$ cat grid.awk
{
  x[NR] = $1
  y[NR] = $2
  z[NR] = $3
}
END{
  for(x0 = -6; x0 <= 6; x0 += 3){
    for(y0 = -6; y0 <= 6; y0 += 3){
      i_grid = 1
      r2_grid = (x[1] - x0) ^ 2 + (y[1] - y0) ^ 2

      for(i = 2; i <= NR; i++){
        r2 = (x[i] - x0) ^ 2 + (y[i] - y0) ^ 2
        if(r2 < r2_grid){
          i_grid = i
          r2_grid = r2
        }
      }

      print x[i_grid], y[i_grid], z[i_grid]
    }
  }
}
$ awk -f grid.awk data.txt > grid.txt

格子点データ grid.txt は，ちょっとかたよった，ゆがんだものになってしまう．

$ psxy grid.txt -Jm1 -R-6/6/-6/6 -Sc0.2 -Cp.cpt -K -Ba2g2 > grid.eps
$ psxy grid.txt -J -R -Sc0.2 -O >> grid.eps

格子データ

$ surface grid.txt -Ggrid.grd -R-6/6/-6/6 -I0.125
$ grdimage grid.grd -Jm1 -Cp.cpt -Ba2 -K > grid2.eps
$ grdcontour grid.grd -J -C0.1 -A0.2 -O >> grid2.eps

格子データ

[ KDOホームにもどる] [ いちからはじめる GMT インデックスにもどる ]

(C) 2013 KDO