Help us understand the problem. What is going on with this article?

震源を地図上にplotする

More than 3 years have passed since last update.

地震続いてますね。
ということで、震源データを可視化してみます。
使う道具は、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日の震源地です。

20160701.png

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

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

おまけ その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)

20160701.con.png

今回のコード

fukuit
最近、事務系の職場に異動したので、職業プログラマではなくなりました。
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away