いちからはじめる GMT その30 - 時系列グラフ

2010年 4月 1日 掲載

ある現象の時間変化を観測したら，当然，時系列グラフを描くだろう． というわけで，今回から何回かにわたり，GMT で時系列グラフを描画する方法を解説する．

まずは，サンプルデータを作る． あつかうのは潮汐である． 某所の主要四分潮成分を合成して，擬似潮汐データを作成する． ちょっと気合いを入れて awkスクリプトを書いてみた （参考文献：日本沿岸潮汐調和定数表／海上保安庁）．

BEGIN{
  
  d = julianday(year, month, day) - julianday(2000, 1, 1.5)

  t = 0
  h = 280.466449 + 0.985647360 * d
  s = 218.316656 + 13.17639647754 * d
  p = 83.353243 + 0.11140352392 * d
  n = 125.044556 - 0.05295376276 * d

  t_hour = 15
  s_hour = 0.54901652
  h_hour = 0.04106864
  p_hour = 0.00464181

  s_m2 = 2 * t_hour - 2 * s_hour + 2 * h_hour + 0 * p_hour
  s_k1 = 1 * t_hour + 0 * s_hour + 1 * h_hour + 0 * p_hour
  s_s2 = 2 * t_hour + 0 * s_hour + 0 * h_hour + 0 * p_hour
  s_o1 = 1 * t_hour - 2 * s_hour + 1 * h_hour + 0 * p_hour

  lon = -140.05
  time_zone = -9

  h_m2 = 50.9
  k_m2 = 152.5
  h_k1 = 25.7
  k_k1 = 178.9
  h_s2 = 24.8
  k_s2 = 181.3
  h_o1 = 20.0
  k_o1 = 160.0
  z0 = 120

  v_m2_0 = v0_m2(t, s, h, p) + u_m2(n) - 2 * lon + s_m2 * time_zone - k_m2
  v_k1_0 = v0_k1(t, s, h, p) + u_k1(n) - 1 * lon + s_k1 * time_zone - k_k1
  v_s2_0 = v0_s2(t, s, h, p) + u_s2(n) - 2 * lon + s_s2 * time_zone - k_s2
  v_o1_0 = v0_o1(t, s, h, p) + u_o1(n) - 1 * lon + s_o1 * time_zone - k_o1

  for(hour = 0; hour < 24; hour++){
    eta_m2 = f_m2(n) * h_m2 * cos(rad(s_m2 * hour + v_m2_0))
    eta_k1 = f_k1(n) * h_k1 * cos(rad(s_k1 * hour + v_k1_0))
    eta_s2 = f_s2(n) * h_s2 * cos(rad(s_s2 * hour + v_s2_0))
    eta_o1 = f_o1(n) * h_o1 * cos(rad(s_o1 * hour + v_o1_0))
    printf "%d-%d-%dT%02d:00:00 %f\n", year, month, day, hour, eta_m2 + eta_k1 + eta_s2 + eta_o1 + z0
  }      
}


function v0_m2(t, s, h, p){
  return 2 * t - 2 * s + 2 * h + 0 * p + 0
}

function v0_k1(t, s, h, p){
  return 1 * t + 0 * s + 1 * h + 0 * p + 90
}

function v0_s2(t, s, h, p){
  return 2 * t + 0 * s + 0 * h + 0 * p + 0
}

function v0_o1(t, s, h, p){
  return 1 * t - 2 * s + 1 * h + 0 * p + 270
}

function f_m2(n){
  tmp = 1.0004
  tmp += -0.0373 * cos (rad(n))
  tmp += 0.0002 * cos(rad(n) * 2)
  tmp += 0.0000 * cos(rad(n) * 3)
  return tmp
}

function u_m2(n){
  tmp = -2.14 * sin(rad(n))
  tmp += 0.00 * sin(rad(n) * 2)
  tmp += 0.00 * sin(rad(n) * 3)
  return tmp
}

function f_k1(n){
  tmp = 1.0060
  tmp += 0.1150 * cos(rad(n))
  tmp += -0.0088 * cos(rad(n) * 2)
  tmp += 0.0006 * cos(rad(n) * 3)
  return tmp
}

function u_k1(n){
  tmp = -8.86 * sin(rad(n))
  tmp += 0.68 * sin(rad(n) * 2)
  tmp += -0.07 * sin(rad(n) * 3)
  return tmp
}

function f_s2(n){
  return 1
}

function u_s2(n){
  return 0
}

function f_o1(n){
  tmp = 1.0089
  tmp += 0.1871 * cos(rad(n))
  tmp += -0.0147 * cos(rad(n) * 2)
  tmp += 0.0014 * cos(rad(n) * 3)
  return tmp
}

function u_o1(n){
  tmp = 10.80 * sin(rad(n))
  tmp += -1.34 * sin(rad(n) * 2)
  tmp += 0.19 * sin(rad(n) * 3)
  return tmp
}

function rad(deg){
  return deg * atan2(0, -1) / 180
}

function julianday(y, m, d){

  if(m <= 2){
    y -= 1
    m += 12
  }

  a = int(y / 100)
  b = 2 - a + int(a / 4)
  c = int(365.25 * y)
  
  return 1720994 + int(30.6001 * (m + 1)) + b + c + d + 0.5
}

長い…．100行以上ある…．気合い入りすぎ…． これを打ち込むのは大変なので， ダウンロード用のファイルを用意した（tide4.awk）．

さて，これを実行してみる．年，月，日を awk の -vオプションで与えると， その日の毎時の潮高が出力される．

$ awk -v year=2010 -v month=4 -v day=1 -f tide4.awk
2010-4-1T00:00:00 65.853609
2010-4-1T01:00:00 80.759239
2010-4-1T02:00:00 109.825872
2010-4-1T03:00:00 145.494661
2010-4-1T04:00:00 178.243813
2010-4-1T05:00:00 199.048868
2010-4-1T06:00:00 201.731466
2010-4-1T07:00:00 184.591789
2010-4-1T08:00:00 150.903318
2010-4-1T09:00:00 108.140360
2010-4-1T10:00:00 66.133637
2010-4-1T11:00:00 34.623863
2010-4-1T12:00:00 20.836895
2010-4-1T13:00:00 27.697640
2010-4-1T14:00:00 53.135288
2010-4-1T15:00:00 90.652210
2010-4-1T16:00:00 131.005283
2010-4-1T17:00:00 164.564596
2010-4-1T18:00:00 183.743269
2010-4-1T19:00:00 184.877006
2010-4-1T20:00:00 169.076418
2010-4-1T21:00:00 141.841345
2010-4-1T22:00:00 111.545634
2010-4-1T23:00:00 87.190762

1カラム目が日時で，2カラム目が潮高（cm）である． 日時は，ごらんの通り「年-月-日T時:分:秒」のフォーマットで出力している． これが，GMT であつかう日時の基本フォーマットである．

では，さっそく時系列グラフを書いてみよう． まずは，tide4.awk の出力 2日分をファイル tide.txt に保存する．

$ awk -v year=2010 -v month=4 -v day=1 -f tide4.awk > tide.txt
$ awk -v year=2010 -v month=4 -v day=2 -f tide4.awk >> tide.txt

これを psxy で描画する． 時系列データをあつかう際には，

• -J の日時データの座標に「T」をつける
• -R の日時の範囲指定は「年-月-日T時:分:秒」のフォーマットで行う

必要がある．

$ psxy tide.txt -Jx8T/0.04 -R2010-4-1T00:00:00/2010-4-3T00:00:00/0/250 > tide.eps

とりあえず，潮汐曲線だけ描いてみた．

awkスクリプトの掲載でページが長くなりすぎたので，今回はここまで．

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

(C) 2010 KDO