2016-04-24

震央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)

地震列島ですね、本当に。

*1:どうでもいいけど、このhtmlファイルを作る作業って人手でやってるんだろうか。全国から上がってくるデータを機械的にDBから吐き出すだけだと思ってたんだけれど。あるいは熊本方面からのデータの吸い上げができなくて、公開が遅くなっているんだろうか。

*2:さらにこうやって作ったCSVファイルをWindows版のRでread.delim()したら文字化けした。Windows版はutf-8に対応してないのかな。

2016-04-18

カーネル密度推定など

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で公開されているよりはマシ

2016-01-20

炊飯器ケーキ(3)

Diary

懲りずに炊飯器でケーキを炊いている。

なんていうか、美味しいワケでもないんだけれど、ご飯を炊く機械でケーキができるという面白さだけの、言わば出落ちみたいなハナシだな。

今回はバナナケーキにした。作り方は、全く同じで、チョコレートやチーズの代わりにバナナを2本分つぶしたものをよく混ぜた。

味は、まあ、バナナの入ったホットケーキって、そのまんま。
焼きあがった時の甘い匂いが期待させるほどの味にならないのは、砂糖が足りない*1のかな。やはり、お菓子は甘くないとな。

*1:そもそも、ホットケーキミックスやバナナに含まれている糖分で足りるんじゃないかと思って、砂糖を入れてない。

2016-01-19

炊飯器ケーキ(2)

Diary

さて、今度は、チョコレートケーキに引き続き、チーズケーキも焼いてみた。

材料

クリームチーズ 100g
生クリーム 100mL
卵 2個
ホットケーキミックス 100g

作り方

ほぼ、チョコレートケーキと同じだ。違いは、チョコレートを湯せんにする変わりに、クリームチーズをよくこねるだけだ。

クリームチーズは、近所のスーパーで売ってるプライベートブランドのものだ。フィラデルフィアのクリームチーズに似てる。

さて、炊飯器のスイッチが切れて、中を見てみたところ、まだまだ生地が固まっていないようだったので、炊飯スイッチをもう一回押してみた。

待つこと30分くらい(忘れた)、今度はちゃんと焼きあがった(炊き上がった)。

さてお味の方は

ホットケーキミックス特有の甘いにおいがプンプンしているが、味はびっくりするくらい甘くない。砂糖を入れないといけなかったようだ。

そして、ホットケーキミックスを焼いたモノに期待されるフワフワ感もなく、ずっしりと重い。炊飯スイッチを2回押したのがよくなかったのかな。あるいは全卵をまぜるんじゃなくて、ちゃんとメレンゲを作るとかしたほうが良かっただろうか。

そして、チーズケーキなのに、チーズ味もほんのかすかにしか感じられない。

基本的にホットケーキミックスを焼いたものなので、不味いとか食べられないとかではないんだけれど、なんだか残念な結果になった。

2016-01-18

炊飯器ケーキ

Diary

以前から、炊飯器でケーキを作ることに憧れていた。スーパーで、板チョコが1枚88円だったので、それでガトーショコラを作ることにした。

材料

板チョコ2枚
生クリーム 100mL
卵 2個
ホットケーキミックス 100g

作り方

卵を冷蔵庫から出しておく。常温に戻すためだ。

板チョコは包丁で刻んで湯せんする。

卵を溶いて、少しずつ加えながら混ぜる。本当は、白身と黄身を分けて、白身はメレンゲにすべきなんだろうけれど、面倒臭い。

さらに、生クリームを少しずつ加えながら混ぜる。

生クリームがいい具合に混ざったら、ホットケーキミックスをダマにならないように混ぜる。

写真を見れば分かるように、ゴムベらみたいなのを持ってないので、カレースプーンで混ぜる。

炊飯器にクッキングシートを敷いて、その上から混ぜたヤツをドバーッと。この後、内釜を何回かトントンとテーブルに数cmの高さから落として、空気を抜く。

後は、炊飯器のスイッチオン。何にも考えずに白米を炊くモードで炊く。

炊けたら、皿に取り出して粗熱を取る。

なかなか悪くない。

さてお味の方は

完全に、ホットケーキミックスにチョコレートが加わっただけの味。うーむ。もっと、濃厚なガトーショコラっぽい味になるかと期待したんだけれど、ホットケーキミックスはホットケーキミックスだった。

とはいえ、コレを冷蔵庫で冷やすと、ちょっといい感じにチョコレートケーキになった。

次は、チーズケーキを作ってみよう。

2016-01-17

WD Elements Portableの故障→保障修理

Diary

仕事で使っているWD Elements Portableというハードディスクが故障して、保障で交換してもらったので、それを記録しておく。

2015年11月頃

2014年初ころに購入したWD Elements Portableがマウントできなくなったり、マウントできてもアクセスできなくなったりする。

これは故障の予兆だと思って、マウントできたときに、robocopyでバックアップを取ろうとしても、エラーが発生してコピーできない。ファイルサイズが数Kbのごくごく小さなファイルであっても、コピー完了までに数十分かかることもある。数時間待っても終わってないことがある。すなわち、「故障の予兆」というよりは、壊れている。

いったん、フォーマットしなおして、その際にCHKDSKを実施すると、マウントできることがあるが、その際にも新しいフォルダを作成しようとしたり、反応が返ってこなくなる。どうやら、ホントの本当に壊れてしまっているようだ。

購入時の資料等がないが、保障は効くだろうか？というのが不安で確認するのが面倒で、ぐだぐだと過ごしてしまったが意を決してサポートに連絡してみることにした。

2015/12/14 製品登録実施

WDのWebSiteにて、製品登録を実施し、故障の可能性があるということで、Trouble Shootingとしてサポートケースに登録する。

2015/12/15 メールにて返信がある

WDのサポート担当者からメールがくる。要約すると、以下のような内容だった。

Data Lifeguard Diagnosticを使って、HDDを初期化したら回復するかどうかを確認すれ

Data Lifeguard DiagnosticでのHDD初期化(WRITE ZEROS)が非常に時間のかかる作業っぽかったので、空いているWindows 7(32bit)機を探して、まずはData Lifeguard Diagnosticのソフトウェアをインストールしようとしたのだが、なぜか起動しなかったので、他に使えるPCを探すまで放置。

翌日、「メールは届きましたか？」とWDのサポートからメールが来たので、「もちろんメールは届いたし、読んだんだけれど、指示されたツールをインストールできるPCがないので、少し時間かかりそうです」と返事をしておく。

確認でき次第再度ご連絡をいただくようお願い致します。
それでは一旦このケースをクローズさせていただきます。

ということになった。

2015/12/21 Data Lifeguard Diagnosticの実施と結果の送信

WDのサポート担当者の指示通りに、Data Lifeguard Diagnosticでフォーマット後に「新しいシンプルボリュームの作成」をしようとしたところ、「新しいシンプルボリュームの作成」が終わらない。仕方ないので、「指示されたことをやろうと思ったけど、無理でした」とサポートに連絡。その際に、名前と住所を合わせて返信した。

2015/12/22 WDから「RMAの確認」というメールがくる

これはお客様のRMA(返品許可)要請の確認です。
RMA年月日 12月 22, 2015
有効期限 1月 22, 2016

お客様のRMAを迅速に処理するため、お客様から発送する製品の梱包方法についてはガイドラインを必ずご確認ください。

ということで、梱包の仕方、発行されたRMA番号を箱に明記することなどの注意書きが書かれたURLが送られてくる。このURLを印刷することで、送り状になるようになっているので、2箇所に署名と署名した日付を書き込む。

すべての製品は、運搬中の静電気（ESD）による影響や破損から保護するために必要な、安定した十分な梱包材を必要とします。必要な梱包手順:
A. ドライブを ESD （帯電防止）バッグに入れます（内蔵ドライブのみ)。
B.厚さ2インチ（約5センチ）のバブルラップまたは固定された発泡緩衝材で包みます。
C.頑丈な段ボールの箱に入れます。クリップボードは輸送中の取扱いに耐えられないため使用しないでください。段ボール箱には損傷のないことと構造的な欠陥のないことを確認してください。注意：ハードドライブを封筒で製品を発送すると、保証が無効となります。
D.商品返品用ラベルを印刷し、箱に貼って下さい。箱の外側、3側面に太字でRMA番号を明記してください。

必要な梱包手順として指示される「厚さ2インチ(約5cm)のバブルラップ」というのが意外と無理難題。プチプチをいろんなところからかき集めることになった。また、厚さ5cmのバブルラップでくるんだWD Elements Portableを入れるのにちょうどよい大きさで頑丈そうな段ボールというのもなかなか見つからないものだ。

ところで、いろいろなWebサイトで、RMAの送り先はシンガポールだと書かれているが、今回指示されたのは東京宛だった。ということで、2015/12/25には発送を完了した。

2015/12/28 RMAのサイトを見てみる

RMA Statusを見てみると、「Receiving Detail 26-Dec-2015」として、こちらから送ったハードディスクのシリアル番号が記載されている。そして、そのすぐ下の欄には、「Shipped Date 28-Dec-2015」として、別のシリアル番号が記載されている。もう、発送された？年末なのと、送り先として職場を指定したこともあって、この時点で、受取は年明けになることを予想する。
ちなみに、Fedexのトラッキング番号が記載されているが、Fedexの画面ではこの番号を入力しても表示されなかった。
なお、翌12/29には、WDから「RMA発送のお知らせ」というメールも来た。

2016/01/07

新しいハードディスクがベトナムから到着した。invoiceに書いてあるshipping dateは03JAN16ということになっているので、WDの倉庫の出庫からFedexに入庫するまでに年末年始を挟んだがために少し時間がかかる結果になったようだが、営業日換算して考えるとすると、非常に迅速に対応してもらえたようだ。