y_uti のブログ

統計、機械学習、自然言語処理などに興味を持つエンジニアの技術ブログです

Dynamic Time Warping による時系列データの類似度計算

台風の経路情報を題材にして、Dynamic Time Warping (DTW) を用いた時系列データの類似度の計算を試してみます。DTW は二つの時系列データの類似度を測る方法の一つで、英語版Wikipedia に簡単な説明と実装例があります。
Dynamic time warping - Wikipedia, the free encyclopedia

過去の台風の経路情報は、各国の機関によって公表されているようです。たとえば、気象庁のデータや、米軍の Joint Typhoon Warning Center (JTWC) という機関のデータが、それぞれ以下のウェブページで公表されています。

これらのデータは各機関が独自の観測によって取得したもので、同一の台風を表す情報でも少しずつ数値がずれています。また、台風と認定される基準も各々異なっており、天気予報で伝えられる「台風何号」という番号自体も機関ごとに異なる場合があるようです。2012 年に上陸した台風 17 号 (JELAWAT) を例に挙げると、気象庁と JTWC のデータはそれぞれ以下のようになっており*1気象庁では台風 17 号、JTWC では 18 号として扱われていることがわかります。

これらの台風の経路を描いてみると、それぞれ以下のようになります。左が気象庁のデータ、右が JTWC のデータです。なお描画には D3.js を利用し、世界地図のデータは Natural Earth からダウンロードしたものを加工して利用しました。

さて、以下では、JTWC のデータの中から、気象庁の T1217.pdf と同一の台風を探してみます*2。観測データに多少のずれがあっても、同じ台風であれば似たような経路情報になっているだろうと考えて、JTWC のデータの中で T1217.pdf との DTW 距離が最小になるものを見つけます。

気象庁のデータの前処理

まず、気象庁の PDF ファイルから緯度経度のデータを抽出します。今回は、pdftotext*3でテキスト情報を抽出し、その後は場当たり的なテキスト処理で何とかする方針で進めます。ファイルをダウンロードしてテキスト変換するまでは以下のとおりです。

$ wget http://www.data.jma.go.jp/fcd/yoho/data/typhoon/T1217.pdf
$ pdftotext T1217.pdf

出力されたテキストの一部分を抜粋すると、以下のようになります。表組みなどは崩れてしまい、このようにぐちゃぐちゃな出力が得られます。

km
hPa m/s
9 20 09 13.7 N 134.8 E 1010 -----15 13.8
133.4
1006 -----21 13.5
132.0
1006 -----21 03 13.3

このような状態から必要な情報を抽出する処理は、データに合わせた工夫が必要になります。このファイルは都合のよいことに、緯度経度のほかに小数点を含むフィールドが無かったので、次のように簡単に抽出できました。

$ tr ' ' '\n' <T1217.txt | grep '\.' | paste -d, - - >T1217.csv

できあがりは以下のようになります。

$ head -n 5 T1217.csv
13.7,134.8
13.8,133.4
13.5,132.0
13.3,131.5
13.1,130.9

JTWC のデータの前処理

同様に、比較対象にする JTWC のデータからも緯度経度の情報を抽出します。JTWC のウェブサイトでは年ごとに全データが圧縮ファイルにまとめられていますので*4、今回は 2012 のデータ全体を利用することにします。

$ wget http://www.usno.navy.mil/NOOC/nmfc-ph/RSS/jtwc/best_tracks/2012/2012s-bwp/bwp.2012.tar.gz
$ tar xf bwp.2012.tar.gz

展開して得られた dat ファイルは以下の形式になっています。各フィールドの定義はウェブサイトに説明がありますが、わざわざ読まなくても第 7, 8 列が緯度経度だろうことは想像できます。

$ head -n 5 bwp012012.dat
WP, 01, 2012021700,   , BEST,   0, 100N, 1168E,  20, 1007, DB,   0,    ,    0,    0,    0,    0,
WP, 01, 2012021706,   , BEST,   0,  99N, 1157E,  20, 1007, DB,   0,    ,    0,    0,    0,    0, 1008,  200,  45,   0,   0,   W,   0,    ,   0,   0,        ONE, S,
WP, 01, 2012021712,   , BEST,   0,  97N, 1144E,  20, 1007, DB,   0,    ,    0,    0,    0,    0, 1008,  200,  45,   0,   0,   W,   0,    ,   0,   0,        ONE, S,
WP, 01, 2012021718,   , BEST,   0,  98N, 1133E,  20, 1007, DB,   0,    ,    0,    0,    0,    0, 1008,  200,  45,   0,   0,   W,   0,    ,   0,   0,        ONE, S,
WP, 01, 2012021800,   , BEST,   0,  98N, 1127E,  20, 1007, DB,   0,    ,    0,    0,    0,    0, 1008,  200,  45,   0,   0,   W,   0,    ,   0,   0,        ONE, S,

次のようにして csv ファイルに変換します。JTWC のデータは、同一時刻の観測が複数行になっていることがあるようなので、最後に uniq でまとめています。なお、この処理は緯度経度が北緯、東経であることを前提にしています。この記事で利用する 2012 年のデータには南緯や西経の行が存在していないことを確認していますが、他の年のデータには含まれていました。いろいろ試してみる場合には注意してください。

$ for f in *.dat; do awk -F, '{ print $7 / 10, $8 / 10 }' $f | tr ' ' , | uniq >${f%.dat}.csv; done

こちらのできあがりは以下のようになります。

$ head -n 5 bwp012012.csv
10,116.8
9.9,115.7
9.7,114.4
9.8,113.3
9.8,112.7

DTW の計算

DTW の計算は、この記事の冒頭に示した Wikipedia のウェブページに記載されているように、簡単に実装できます。今回は PHP で書いてみました。関数 dtw で、引数 $a と $b の DTW 距離を計算します。二点間の距離は単純にユークリッド距離で求めています*5

<?php
function dtw($a, $b, $distance = 'euclid') {
    $d = array_fill(0, count($a) + 1, array_fill(0, count($b) + 1, INF));
    $d[0][0] = 0;
    for ($i = 1; $i <= count($a); ++$i) {
        for ($j = 1; $j <= count($b); ++$j) {
            $d[$i][$j] =
                min([$d[$i - 1][$j - 1], $d[$i][$j - 1], $d[$i - 1][$j]])
                + $distance($a[$i - 1], $b[$j - 1]);
        }
    }
    return $d[count($a)][count($b)];
}

function euclid($a, $b) {
    return sqrt(array_sum(array_map(function ($x, $y) {
        return pow($x - $y, 2);
    }, $a, $b)));
}

function readCsv($filename) {
    return array_map(
        function ($line) { return explode(',', $line); },
        file($filename, FILE_IGNORE_NEW_LINES));
}

echo dtw(readCsv($argv[1]), readCsv($argv[2])) . "\n";

このプログラムを次のように繰り返し実行すれば、T1217.csv と JTWC の各データ間の DTW 距離を計算できます。実行例では、計算結果を DTW 距離の昇順にソートしています。T1217.csv に最も近い系列は bwp182012.csv という結果になりました。最初に示したように、T1217.csv と bwp182012.csv は同一の台風 (JELAWAT) を表すデータであり、DTW によって正しい対応付けが得られています。

$ for f in bwp*.csv; do echo $(php dtw.php T1217.csv $f) $f; done | sort -nsk1,1
188.48919354809 bwp182012.csv
356.67621274902 bwp042012.csv
473.55617321595 bwp052012.csv

... (省略)

1777.8991470182 bwp252012.csv
1845.8456855743 bwp132012.csv
1905.2481949742 bwp022012.csv

JTWC の 2012 年のデータの中で、T1217.csv との DTW 距離が近い台風 3 つと距離が遠い台風 3 つをそれぞれ図示すると、次のようになります。赤色で示した軌跡が T1217.csv です。たしかに、DTW 距離が近い方が似たような軌跡を描いているように見えます。

*1:どちらも前述のウェブページからリンクを辿って到達できます。

*2:ここで示したように、正解は bwp182012.dat になります。

*3:CentOS 7 では、poppler-util パッケージに含まれています。

*4:年によっては圧縮ファイルが提供されていなかったり、圧縮形式が tar.gz ではなく zip だったりします。

*5:台風の経路といった規模でユークリッド距離を用いるのは不適切だと思いますが、今回は特に気にしないことにします。