震央plot
データのダウンロード
ダウンロードできるデータは、 http://www.data.jma.go.jp/svd/eqev/data/daily_map/20160414.html のように、日付別になっているので、とりあえず、4/14〜4/19までのデータをダウンロードしておく*1。
http://www.data.jma.go.jp/svd/eqev/data/daily_map/20160414.html http://www.data.jma.go.jp/svd/eqev/data/daily_map/20160415.html http://www.data.jma.go.jp/svd/eqev/data/daily_map/20160416.html http://www.data.jma.go.jp/svd/eqev/data/daily_map/20160417.html http://www.data.jma.go.jp/svd/eqev/data/daily_map/20160418.html http://www.data.jma.go.jp/svd/eqev/data/daily_map/20160419.html
データの加工
とりあえず、いらんところを捨てるだけのperl script。
#!/usr/bin/perl use strict; use warnings; use utf8; opendir my $dh, "." or die $!; for my $infile(grep /\.html$/, readdir($dh)){ my $outfile = $infile; $outfile =~ s/html/csv/; open my $fhIn, '<:utf8', $infile or die $!; open my $fhOut, '>:utf8', $outfile or die $!; my $flg = 0; my $cnt=0; while (my $line = <$fhIn>) { chomp($line); $flg = 0 if $line =~ /<\/pre>/; if ($flg == 1){ if ($cnt==0){ print $fhOut "date,time,latitude,longitude,depth,magnitude,epicenter\n"; }else{ if (($line !~ /^\-+$/) && length($line) != 0){ my @tmp = split /\s+/, $line; my $yy=$tmp[0]+0; my $mm=$tmp[1]+0; my $dd=$tmp[2]+0; my $hm=$tmp[3]; my $ss=$tmp[4]+0.0; my $lat1; my $lat2; my $lon1; my $lon2; my $depth; my $magnitude; my $epicenter; if ($#tmp == 11){ $lat1 = $tmp[5]; $lat2 = $tmp[6]; $lon1 = $tmp[7]; $lon2 = $tmp[8]; $depth = $tmp[9]; $magnitude = $tmp[10]; $epicenter = $tmp[11]; }elsif($#tmp == 10){ if ($tmp[5] =~ /N/){ $lat1 = $tmp[5]; $lat2 = ''; $lon1 = $tmp[6]; $lon2 = $tmp[7]; }else{ $lat1 = $tmp[5]; $lat2 = $tmp[6]; $lon1 = $tmp[7]; $lon2 = ''; } $depth = $tmp[8]; $magnitude = $tmp[9]; $epicenter = $tmp[10]; }else{ $lat1 = $tmp[5]; $lat2 = ''; $lon1 = $tmp[6]; $lon2 = ''; $depth = $tmp[7]; $magnitude = $tmp[8]; $epicenter = $tmp[9]; } my $date = $yy . "-" . $mm . "-" . $dd . "," . $hm . ":" . $ss; my $latitude = $lat1 . $lat2; my $longitude =$lon1 . $lon2; $latitude =~ s/'[EN]//; $latitude =~ s/\./\:/; $latitude =~ s/°/\:/; $longitude =~ s/'[EN]//; $longitude =~ s/\./\:/; $longitude =~ s/°/\:/; print $fhOut "$date,$latitude,$longitude,$depth,$magnitude,$epicenter\n"; } } $cnt++; } $flg = 1 if $line =~ /<pre>/; } close $fhOut; close $fhIn; print "$infile => $outfile\n"; } closedir $dh; __END__
固定長データだから、split(/\s+/, ...) するのは、得策じゃなかったっぽい。変な場合分けをするはめになってしまった。*2
地図データへの重畳表示
あとは、http://d.hatena.ne.jp/fukuit/20160418/1460976967 同様にggmapする。複数csvファイルに対応させたのと、日付別に色分けするところだけ追加。
files = list.files(pattern=".csv") for( i in 1:length(files)){ d.tmp <- read.delim(files[i], header=T, sep=',') if ( i == 1 ){ dat <- d.tmp }else{ dat <- rbind(dat,d.tmp) } } 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) cols <- as.numeric(dat$date) library(ggmap) center <- geocode('tajimi') map <- get_map(c(center$lon, center$lat), zoom=5, maptype='roadmap') p<-ggmap(map)+geom_point(data=df, aes(x=x, y=y), size=magnitude, color=cols, alpha=0.2)
地震列島ですね、本当に。
