地震続いてますね。
ということで、震源データを可視化してみます。
使う道具は、perlとRです。
震源データをダウンロードする
震源地のデータは、気象庁からダウンロード可能です。
気象庁震源リストからは、震源地とマグニチュードと地震発生日時が、日付別にアクセス可能になっています。
実際のデータは、CSV等になっておらず、日付別のhtml中に埋め込まれているので、これを切り出すscriptをperlで以下のように作成しました。
#!/usr/bin/perl
use strict;
use utf8;
use LWP::Simple;
use Encode;
use DateTime;
if (@ARGV == 1){
my $today = $ARGV[0];
die if ($today !~ /[0-9]{8}/);
retrieveData($today);
}else{
for(my $i=1;$i<31; $i++){
my $dt = DateTime->now;
my $ddt= $dt->subtract( days=> $i );
my $tdate = $ddt->strftime("%Y%m%d");
retrieveData($tdate);
}
}
sub retrieveData{
my $tdate = shift;
my $outfile = $tdate . ".txt";
eval{
if (! -f $outfile){
my $baseurl = "http://www.data.jma.go.jp/svd/eqev/data/daily_map/";
my $url = $baseurl . $tdate . ".html";
my $html = get($url) or die "$tdate : $!";
my @lines = split /\n/, $html;
my $flg = 0;
my $cnt=0;
my $outline = "";
open my $fhOut, '>:encoding(utf-8)', $outfile or die $!;
for (my $l=0;$l<=$#lines;$l++){
my $line = $lines[$l];
chomp($line);
$flg = 0 if $line =~ /<\/pre>/;
if ($flg == 1){
if ($cnt==0){
$outline = "date\ttime\tlatitude\tlongitude\tdepth\tmagnitude\tepicenter\r\n";
}else{
$outline .= formatLine($line);
}
$cnt++;
}
$flg = 1 if $line =~ /<pre>/;
}
print $fhOut $outline;
close $fhOut;
print " => $outfile\n";
}
};
if($@){
print STDERR "ERROR: $@\n";
}
}
sub formatLine{
my $line = shift;
my $outline;
if (($line !~ /^\-+$/) && length($line) != 0){
my $yy = substr($line,0,4)+0;
my $mm = trim(substr($line,5,2))+0;
my $dd = trim(substr($line,8,2))+0;
my $hm = substr($line, 11,5)."";
my $ss = trim(substr($line, 17,2))+0;
my $date = sprintf("%d-%02d-%02d¥t%s:%02d", $yy, $mm, $dd, $hm, $ss);
my $lat1 = trim(substr($line,22,3))+0.0;
my $lat2 = substr($line,26,2)+0.0;
my $lat3 = substr($line,29,1)+0.0;
my $latitude = $lat1+$lat2/60+$lat3/60**2;
my $lon1 = trim(substr($line,33,3))+0.0;
my $lon2 = substr($line,37,2)+0.0;
my $lon3 = substr($line,40,1)+0.0;
my $longitude = $lon1+$lon2/60+$lon3/60**2;
my $depth = trim(substr($line,44,4));
my $magnitude = trim(substr($line, 49,7));
my $epicenter = substr($line, 58);
$outline = "$date\t$latitude\t$longitude\t$depth\t$magnitude\t$epicenter\r\n";
}
return $outline;
}
sub trim{
my $str = shift;
$str =~ s/\s+//g;
return $str;
}
__END__
事前に、libwww-perlとかDateTimeとか、perl -MCPAN -eshellして、install しておく必要があります。
実行時に引数として、YYYYMMDD形式で日付を与えると、その日のデータをダウンロードします。引数なしの場合は、今日から数えて過去31日分のデータを取得します。
取得したデータはtab区切りテキストで保存します。
緯度・経度のデータが、「36度47.5分」のような形式で記載されているので、36+47/60+5/602のようにして変換します。
Rで処理する
タブ区切りテキストを読む
Rで、カレントフォルダにあるありったけのtxtファイルを読みこみ、順に処理します。
files = list.files(pattern = ".txt")
for (i in 1:length(files)) {
fname <- gsub('.txt', '', files[i])
dat <- read.delim(files[i], header = T, sep = '\t')
longitude <- c()
latitude <- c()
magnitude <- c()
for (j in 1:nrow(dat)) {
longitude[j] <- as.numeric(as.character(dat$longitude[j]))
latitude[j] <- as.numeric(as.character(dat$latitude[j]))
if (as.character(dat$magnitude[j]) == '-') {
magnitude[j] <- 0
} else{
if (as.numeric(as.character(dat$magnitude[j])) < 0) {
magnitude[j] <- 0
} else{
magnitude[j] <- as.numeric(as.character(dat$magnitude[j]))
}
}
}
dat$longitude <- longitude
dat$latitude <- latitude
dat$magnitude <- magnitude
}
Rでread.delim()でファイルを読みこんだ時に、各列の要素がfactor形式で読み込まれてしまうことがあります。この時、安易にas.numeric()で整数値にすると、エラい目にあうので、as.character()したものをas.numeric()にしています。
また、magnictudeの値は、負の数値だったり、'-'しか入ってなかったりすることがあるので、そういう場合は0にしています。
以上で、datという変数名のdata.frameに日時と震源地とマグニチュードと深さのデータが読み込まれました。
ggmap()で地図上に震源を書き込む
ggmap()というのは、地図用に特化したggplot2のwrapperなようです。
install.packages("ggmap")
library(ggmap)
震源のデータがセットされているdatを使って、ggmapで地図上に震源地をプロットします。そのために、まず、地図をGoogle Mapからダウンロードします。
今回は、日本全体の地図にするために、日本の中心地(ウソ)である多治見市を中心にズームレベル5の地図データをダウンロードします。
center <- geocode('tajimi')
map <- get_map(c(center$lon, center$lat), zoom = 5, maptype = 'roadmap')
次に、元のデータファイルの名前を使って、地図画像ファイルを保存するためのファイル名を生成します。
fname <- gsub('.txt', '', files[i])
png.japan <- paste('japan', paste(fname, 'png', sep = '.'), sep = '/')
あとは、ggplot2風に、地図上に震源の位置をプロットします。ちなみに、プロットする点の大きさはmagnitudeの値を使います。
p <- ggmap(map)+
geom_point(data=dat, aes(x=longitude, y=latitude), size=magnitude, colour='red', alpha=0.7)+
ggtitle(fname)
ggsave(png.japan, plot=p)
というワケで、今年の7月1日の震源地です。
日本に安全な場所というのは、なかなかなさそうですね。
library(ggmap)
center <- geocode('tajimi')
map <- get_map(c(center$lon, center$lat), zoom = 5, maptype = 'roadmap')
files = list.files(pattern = ".txt")
for (i in 1:length(files)) {
dat <- read.delim(files[i], header = T, sep = '\t')
longitude <- c()
latitude <- c()
magnitude <- c()
for (j in 1:nrow(dat)) {
longitude[j] <- as.numeric(dat$longitude[j])
latitude[j] <- as.numeric(dat$latitude[j])
if (as.character(dat$magnitude[j]) == '-') {
magnitude[j] <- 0
} else{
if (as.numeric(as.character(dat$magnitude[j])) < 0) {
magnitude[j] <- 0
} else{
magnitude[j] <- as.numeric(as.character(dat$magnitude[j]))
}
}
}
dat$longitude <- longitude
dat$latitude <- latitude
dat$magnitude <- magnitude
fname <- gsub('.txt', '', files[i])
png.japan <- paste(fname, 'png', sep = '.')
p <- ggmap(map)+
geom_point(data=dat, aes(x=longitude, y=latitude), size=magnitude, colour='red', alpha=0.7)+
ggtitle(fname)+
stat_density2d(data=dat, aes(x=longitude, y=latitude), geom='polygon', alpha=0.2)
ggsave(png.japan, plot=p)
}
おまけ その1
日付別に作成した画像から、animated gifを作ると、このようになります。
おまけ その2
震央をプロットしてから、ggmapでkernel密度推定をしてcontour図を描くこともできます。
p<-ggmap(map)+geom_point(data=dat, aes(x=longitude, y=latitude), alpha=0.1)+
stat_density2d(data=dat, aes(x=longitude, y=latitude), geom='polygon', alpha=0.2)
