震源を地図上にplotする

地震続いてますね。
ということで、震源データを可視化してみます。
使う道具は、perlとRです。

# 震源データをダウンロードする
震源地のデータは、気象庁からダウンロード可能です。

[気象庁震源リスト](http://www.data.jma.go.jp/svd/eqev/data/daily_map/index.html)からは、震源地とマグニチュードと地震発生日時が、日付別にアクセス可能になっています。

実際のデータは、CSV等になっておらず、日付別のhtml中に埋め込まれているので、これを切り出すscriptをperlで以下のように作成しました。

```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/60<sup>2</sup>のようにして変換します。

# Rで処理する
## タブ区切りテキストを読む
Rで、カレントフォルダにあるありったけのtxtファイルを読みこみ、順に処理します。

```R
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なようです。

```R
install.packages("ggmap")
library(ggmap)
```

震源のデータがセットされているdatを使って、ggmapで地図上に震源地をプロットします。そのために、まず、地図をGoogle Mapからダウンロードします。

今回は、日本全体の地図にするために、日本の中心地(ウソ)である多治見市を中心にズームレベル5の地図データをダウンロードします。

```R
center <- geocode('tajimi')
map <-   get_map(c(center$lon, center$lat), zoom = 5, maptype = 'roadmap')
```

次に、元のデータファイルの名前を使って、地図画像ファイルを保存するためのファイル名を生成します。

```R
fname <- gsub('.txt', '', files[i])
png.japan <- paste('japan', paste(fname, 'png', sep = '.'), sep = '/')
```

あとは、ggplot2風に、地図上に震源の位置をプロットします。ちなみに、プロットする点の大きさはmagnitudeの値を使います。

```R
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日の震源地です。

![20160701.png](https://qiita-image-store.s3.amazonaws.com/0/15553/c3065255-b55a-f7fb-8e9e-03bf40892a34.png "20160701.png")

日本に安全な場所というのは、なかなかなさそうですね。

```R
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を作ると、このようになります。
![201604.j.gif](https://qiita-image-store.s3.amazonaws.com/0/15553/58eea2c7-bebf-df18-2f3e-100d47e24e6b.gif "201604.j.gif")

## おまけ その2
震央をプロットしてから、ggmapでkernel密度推定をしてcontour図を描くこともできます。

```R
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)
```

![20160701.con.png](https://qiita-image-store.s3.amazonaws.com/0/15553/313076a4-f853-b725-a02c-b83e468c9ec0.png "20160701.con.png")

## 今回のコード

- https://github.com/fukuit/epicenter_plot