気象庁ウェブサイトの「過去の気象データ ダウンロード」の機能をつかって、 長崎、大阪、東京の3地点の、2020年の1年間の、1時間ごとの気圧のデータをダウンロードした。 [気象観測データ (1地点) のグラフ作成と基本的な統計処理] と同様だが、日付を、そこでは年・月・日にわけて別々の列にしたのに対して、 ここでは「日付リテラル」のままダウンロードした。
データをすなおにPandasのデータフレームに読みこんでみる。 日付は、Python の datetime パッケージ の datetime 型に変換する。 (気象庁の「日付リテラル」は datetime パッケージの書きだす形とすこしちがうのだが、 この程度のちがいは自動処理されるようだ。) 気圧の値は浮動小数点数として読む。
読みこんだデータを単純に時系列グラフにしてみる。 グラフの凡例をまだつけていないが、 オレンジ色=長崎、緑=大阪、紫=東京である。
ScyPy の signal パッケージの welch という関数で、 Welch の方法によって時系列のパワースペクトルをもとめることができる。 これを、長崎の気圧に適用してみた。 データフレームの気圧の列を、NumPy の配列につめかえてあたえた。 データフレームの日付の列はつかっておらず、 sampling frequency (welchのマニュアルぺージでの変数名 fs) をあたえる。 fs の単位は本来 Hz (つまり 1/秒) だが、 ここのデータの時間間隔は1時間で、見たい周期は1日前後なので、 fs の単位を「1/日」としてつかうことにし、24という数値をあたえた。
単純に適用したのでは、意味のある結果が出なかった。 この気圧データには欠損があり、 csvファイル上では null string (コンマどうしのあいだに文字がない) になっている。 これを補間してやったら、パワースペクトルらしい結果が出た。 補間は、Pandas の interpolate メソッドを利用し、線形補間によった。
結果をみると、明確なピークが、周期が数日 (未確認だが7日ぐらいか) のところと、 1/2日 (グラフの横軸の「2」)のところにある。
scipy.signal の coherence で、Welchの方法による2つの時系列のコヒーレンスがえられる。 また、csd で、クロススペクトルがえられる。 このクロススペクトルは複素数なので、その絶対値と偏角を計算し、それぞれグラフにした。 (偏角の計算は、cmath パッケージの phase 関数によった。 結果の数値の範囲は [-π, π] である。 この偏角を「位相」としたが、クロススペクトルの位相の定義にあっているか、確認が必要である。)
この例で、わたしは、温帯低気圧をともなう偏西風波動が、 周期が数日の振動数帯でコヒーレンスがあり、位相が西のほうが早く東のほうが遅い、 という形でみられるだろうと期待していた。 結果にそれらしいものはあるが、期待どおりになっているか、確認の途上である。
位相は π が 半周期なので、0.8 は 約 8分の1 周期ということになる。 長崎と東京との距離は約 1000 km なので、この位相差が波のあらわれならば、 波長が約 8000 km ということになるが、 温帯低気圧をともなう偏西風波動にしては、波長が長すぎるようだ。 温帯低気圧をともなう偏西風波動と、長崎と東京で同位相になるような別の振動現象とが まざった結果だとすれば説明できるのかもしれない。(まだよくわからない。)