震源地mapを作る

はじめに

Rを使えば震源地がどう移り変わって行くか可視化できることを知ったので作ってみました.Rでスクリプトを書いてみたのは初めてなのでわからないことだらけだったので,調べたことを自分のメモがわりにまとめてみます.ソースと使い方はgithubに載せています.

準備

Rは対話モードで使うことが多いと思うのですが,今回はスクリプトとしてまとめました.スクリプトとして動かすさいにはrscriptを使います.さらに,今回はグラフ作成にggplot2,地図データとしてmaps,時系列データの取扱にxtsを使います.zooも読み込んでいますが,xtsを読み込むと自動で読み込むので,必要ないです.後で修正しないと.ライブラリは
対話モードで"install.packages()"メソッドを使って読み込むのが楽でしょう.例えば

> install.packages("ggplot2")

と実行します.そうするとどこのミラーからライブラリを落とすか聞いてくるので,適当に選べばインストールしてくれます.毎回聞かれるのが面倒な時は~/.Rprofileに

options(repos="http://cran.cnr.Berkeley.edu")

と記述すると,ここに記載したurlから落としてくれます.なんだか国内ミラーの調子が悪かったのでberkeleyを選んだ時の例です.

#!/usr/bin/env rscript
library(ggplot2)
library(maps)
library(zoo)
library(xts)

normDepthとnormMagnitudeは関数で,それぞれ震源地の深さと規模をグルーピングしています.図示するときの利便性を確保するため,というか,震源地に深さや規模に応じた点を描く際,毎日同じ基準で描くために使ったトリックです.別の方法でmax/minを指定すれば同じ規模なら同じ大きさの点が打てるかもしれないですが,やりかたがわかりませんでした.以下はnormMagnitudeの関数定義です.

normMagnitude <- function(x){
     if (x > 5) return(5)
     if (x > 4) return(4)
     if (x > 3) return(3)
     if (x > 2) return(2)
     1
}

データ読み込み

ここの箇所がデータ変形の肝になります.まず,csv形式の震源地データをネットから取得してdata.frame形式で格納します

catalog <- read.csv(
     "http://earthquake.usgs.gov/earthquakes/catalogs/eqs7day-M2.5.txt"
)

続いて正規表現を使ってjapan regionのデータだけを抽出し,data.frame形式で格納します

japan <- catalog[grep("Japan", catalog$Region),]

Datetime列はUTC時刻表記なのでこれをJSTに変換します.これが最初のはまりどころでした.色々やってみた結果60*60*9を足すとJST表記に変わることが偶然判ったんだけど,そんなことがどこに書いてあるのか見つけられませんでした.

time <- strptime(japan$Datetime, "%A, %B %d, %Y %H:%M:%S UTC") + 60*60*9

深さと規模を正規化して,JST表記の時刻と緯度経度を組み合わせてdata.frame形式で格納します

Depth <- sapply(japan$Depth, normDepth)
Magnitude <- sapply(japan$Magnitude, normMagnitude)
quake <- cbind(time, japan[5:6], Depth, Magnitude)

ここまでで元データが完成しました.

図の下準備

ここは定型処理です.mapsというライブラリは米国をターゲットにしているようで日本地図は概要しかありません.他にもっと詳細なデータを持っているライブラリがあるようですが,今回は調べていません.以下は東経120度から150度,北緯25度から50度までの地図データを抽出しています.

xydata <- as.data.frame(
     map("world", xlim = c(120, 150), ylim = c(25, 50), plot = FALSE)
     [c("x", "y")]
)
map <- ggplot(xydata, aes(x, y))
map <- map + geom_path()

この部分だけを対話モードで実行することもできます.その際は

> print(map)

と実行すれば該当エリアの地図を表示することができます.

日々のデータに分離

時系列データを扱うのはxtsが得意なのでdata.frameから変換します.直接変換はできないようなので,zooを経由します.これもノウハウです.一旦xts形式に変換すると,splitを使って日ごとのリストに分割することができます.リストの各要素に対してplotdataという関数を適用して図を作成しています.

quake.xts <- as.xts(read.zoo(quake))
days <- split(quake.xts, "days")
for ( day in days ) plotdata(day)

plotdata(day)

この部分はきれいなスクリプトがかけていません.indexというグローバル変数でファイル名を作っているところがいまいち.img%03d.pngという名前でファイルを作成しているのはffmegで動画を作るためです.

index <- 1

plotdataは日ごとに分割されたxts形式のデータです.ggplot2はxts形式を扱えないので,data.frameに戻しています.coredata()というのがそのためのメソッドで,これを使うと各要素がfactor型になるので,必要な項目に関してはas.numeric()を使って数字型に戻しています.

plotdata <- function(x){
      j <- data.frame(coredata(as.numeric(x$Lat)),
          coredata(as.numeric(x$Lon)),
          coredata(x$Magnitude),
          coredata(x$Depth)
     )

これでdata.frame型に変換できましたが,列名がごちゃごちゃしているので,すっきりさせます.

     colnames(j) <- c("Lat", "Lon", "Magnitude", "Depth")

震源地をプロットします.サイズと色をfactorで指定しているのは,その日にどんな地震があったとしても,同じ規模なら同じサイズ,同じ深さなら同じ色で描くためのノウハウです.

     map.point <- map + geom_point(
          data = j,
          aes(x = Lon, y = Lat,
               size = factor(Magnitude, c(1, 2, 3, 4, 5, 6, 7, 8, 9)),
               colour = factor(Depth, c(0, 10, 20, 30, 50, 100, 200, 700))
          )
     )

説明表記をすっきりさせて

     map.point <- map.point + scale_colour_discrete("Depth")
     map.point <- map.point + scale_size_discrete("Magnitude")

ファイルに落とします

     png(filename = paste(sprintf("data/img%03d", index), "png", sep = "."))
     index <<- index + 1
     print(map.point)
     dev.off()
}

これでできあがり.dataディレクトリ配下に地震データが... あれ,6個しかない.なぜだ?調べなきゃ.後はdataディレクトリに移動して

ffmpeg -r 1 -f image2 -i img%03d.png  -r 1 quakemapjp.mp4

を実行すれば震源地の推移を表す動画ができあがります.

おわりに

実行すると色々ワーニングがでたり,上に書いたように1週間分のデータのはずなのに6個しか出力してなかったりまだまだですが,とりあえず動いたのでまとめてみました.随時github上で修正していくつもりです.って,Rネタに興味を持ってくれる人ってどれくらいいるんだろうかなぁ.それはそれとして,mp4に変換するとはてダにあげられないのに気づいたので,gifアニメにできるかどうかも試してみます.

追加

スクリプトで作成したpngファイルをimagemagickについてくるconvertコマンドを使ってanimation GIFに変換しました.コマンドは

convert -delay 50 img*.png quake.gif

このようになります.

さらに追加

convertでdelayオプションに渡すパラメータを間違ってました.修正して画像も変更しました.秒あたり2日分の地図を表示するようにしています.