file-glob こと k.daibaの日記 このページをアンテナに追加 RSSフィード

2012-03-19 震源地mapを作る

震源mapを作る

はじめに

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

f:id:kdaiba:20120322144129g:image

準備

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日分の地図を表示するようにしています