カーネル密度推定など

R Diary

ちょっと仕事で必要があって、ggplotを使ったりしたカーネル密度推定と、それを二次元座標系で図示化するのをやってたんだけれど、仕事じゃなくても二次元の地図上にマッピングして見てみたくなるものがあった。

データをダウンロードする

気象庁のサイトで、震源の情報を探してみたら、機械可読性がそんなに高くないデータ*1が見つかった。
http://www.data.jma.go.jp/svd/eqev/data/daily_map/20160414.html
ひょっとすると、APIとかあってダウンロードできるのかもしれない。

このデータをまずExcelにコピペする。



2016  4 14 00:00 56.1  32°28.8'N 131°58.1'E   18     0.7  日向灘                     

2016  4 14 00:06 32.2  42°44.0'N 145°14.9'E   45     1.5  釧路沖

...
固定長データになっているので、Excel上でセルに分割する。時分秒のデータは、空白を「:」に置換したあと、分と秒の間隔が空白2つになっている箇所があるので、「::」を「:」に置換したあと、セルの書式を「ユーザー定義」の「hh:mm:ss.0」にしておく。

緯度・経度は、「°」と「.」を「:」に置換した上で、「'E」と「'N」を削除する。

これを、csv形式で保存する(例: 20160414.csv)。

Rで地図上に震央をマッピングする

地図は、ggmapを使うと恐ろしく簡単にGoogleMapからデータを取ってきて、ソレに書き込むことができる。

dat <- read.delim('20160414.csv', head=T, sep=',')
longitude <- c()
latitude <- c()
magnitude <- c()
for(i in 1:nrow(dat)){
# 「時:分:秒」形式の緯度・経度を10進数に変換
  a <- strsplit(as.character(dat$longitude[i]), ':')[[1]]
  longitude[i] <- as.numeric(a[1])+as.numeric(a[2])/60+as.numeric(a[3])/60/60

  a <- strsplit(as.character(dat$latitude[i]), ':')[[1]]
  latitude[i] <- as.numeric(a[1])+as.numeric(a[2])/60+as.numeric(a[3])/60/60

# magnitudeは「-」だったり負の数があったりするので、そういうのは0ってことにする
  if ( as.character(dat$magnitude[i]) == '-' ) {
    magnitude[i] <- 0
  }else{
     if ( as.numeric(as.character(dat$magnitude[i])) < 0 ) {
        magnitude[i] <- 0
     }else{
        magnitude[i] <- as.numeric(as.character(dat$magnitude[i]))
     }
  }
}
df <- data.frame(x=longitude, y=latitude)

library(ggmap)
center <- geocode('kumamoto')
map <- get_map(c(center$lon, center$lat), zoom=10, maptype='roadmap')
p<-ggmap(map)+geom_point(data=df, aes(x=x, y=y), size=magnitude, color='red')

便利な世の中だ。

ggplot2がベースなので、こんなことをすると震度を観測した地点に対してカーネル密度推定ができる。

p+stat_density_2d(df, aes(x=x,y=y)

日本全国を対象にするときには、25年くらい前に首都移転論があったときに、「東京から東濃へ」というキャンペーンを張ってた多治見市に敬意を表して、多治見の座標を使う。zoomレベルは5。

center <- geocode('tajimi')
ggmap((map <- get_map(c(center$lon, center$lat), zoom=5, maptype='roadmap')))+
    geom_point(data=df, aes(x=x, y=y), size=magnitude/2, color='red')+
    stat_density2d(data=df, aes(x=df$x, y=df$y), geom='polygon', alpha=0.2)

4/14の24時間だけでも、意外といろんなところで有感地震は起きてる。

*1:PDFで公開されているよりはマシ