R

気象庁震源リスト02 ( 平成28年（2016年）熊本地震 )

mapdata , xts , ggplot2 , gganimate パッケージ

（データ）
気象庁 震源リスト

(過去の記事)
気象庁震源リスト01 ( 平成28年（2016年）熊本地震 ) 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
library(mapdata)
#
eqdata<-read.table(text="",col.names=c("time", "longitude", "latitude", "depth", "mag"))
ymd<-seq(as.Date("2016-04-14"),as.Date("2016-04-28"), by="1 day")
for (i in 1:length(ymd)){
date<-gsub("-","",ymd[i])
#ウェブスクレイピング
url<-paste0("http://www.data.jma.go.jp/svd/eqev/data/daily_map/",date,".html")
#
doc<-readLines(url,encoding ="UTF-8")
kensaku<-paste0(substr(date,1,4)," ",gsub("\\<0"," ",substr(date,5,6))," ",gsub("\\<0"," ",substr(date,7,8)))
#
x<-doc[grep(paste0("\\<",kensaku),doc)]
#全角文字 ° が入っているので半角スペース２つに変換する
x<-gsub("°", "  ",x)
time<-paste0(substr(x,1,4),"-",substr(x,6,7),"-",substr(x,9,10)," ",substr(x,12,16),":",substr(x,18,21))
#
#南緯、西経の場合も考慮に入れると
latitude = as.numeric(paste0(substr(x,23,23),as.character(as.numeric(substr(x,24,25))+as.numeric(substr(x,28,29))/60+as.numeric(substr(x,31,31))/600)))
longitude = as.numeric(paste0(substr(x,34,34),as.character(as.numeric(substr(x,35,37))+as.numeric(substr(x,40,41))/60+as.numeric(substr(x,43,43))/600)))
depth<-as.numeric(substr(x,48,50))
mag<-as.numeric(substr(x,55,58))
eq<-data.frame(time,longitude,latitude,depth,mag)
eqdata<-rbind(eqdata,subset(eq,mag>=-4))
}
#
#重複した行を取り除く
eqdata<-unique(eqdata)
summary(eqdata)
#save(eqdata,file="eqdata20160414_0419.RData")
#load("eqdata20160414_0419.RData")
#
#震源の深さ
# depth<10km
eqdata1<-subset(eqdata,depth<10)
eqdata1$depth_rank<-" 10km 未満"
eqdata1$col<-"red"
# 10<=depth<20
eqdata2<-subset(eqdata,depth>=10 & depth<20)
eqdata2$depth_rank<-" 10km以上20km未満"
eqdata2$col<-"orange"
# 20<=depth<40
eqdata3<-subset(eqdata,depth>=20 & depth<40)
eqdata3$depth_rank<-" 20km以上40km未満"
eqdata3$col<-"gold4"
# 40<=depth<80
eqdata4<-subset(eqdata,depth>=40 & depth<80)
eqdata4$depth_rank<-" 40km以上80km未満"
eqdata4$col<-"green"
# 80<=depth<150
eqdata5<-subset(eqdata,depth>=80 & depth<150)
eqdata5$depth_rank<-" 80km以上150km未満"
eqdata5$col<-"blue"
# 150<=depth
eqdata6<-subset(eqdata,depth>=150)
eqdata6$depth_rank<-"150km以上"
eqdata6$col<-"purple"
#
eq<-rbind(eqdata1,eqdata2,eqdata3,eqdata4,eqdata5,eqdata6)
#並べ替え:マグニチュードの昇順
sortlist <- order(eq$mag)
eq<- eq[sortlist,]
##
#mag>=1の地震を抽出
subeq<-subset(eq,mag>=1)
#mag<1
mag1und<-subset(eq,mag<1)

地図を描く

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
par(family = "TakaoExMincho")
#png("eq20160414_0428_01.png",width=1000,height=800)
map('japan',xlim=c(130,132),ylim=c(31.8,33.8),fill=T,col="gray90",mar = c(5,4,4,15))
map('japan',col="gray20",boundary =T,add=T)
points(mag1und$longitude,mag1und$latitude,cex=0.2,bg=adjustcolor(subeq$col,1),pch=21,col=adjustcolor(subeq$col,1))
points(subeq$longitude,subeq$latitude,cex=(subeq$mag-0.2)/2,bg=adjustcolor(subeq$col,1),pch=21,col="black")
points(c(130.1894,129.8372,132.3111),c(31.83361,33.51556,33.49083),cex=2,col="black",bg=c(adjustcolor("gray30",0.8),rep(rgb(1,1,1,0.8),2)),pch=24) 
par(xpd=T)
#
cexsize=(c(1,3,5,7)-0.2)/2
legend(par("usr")[2]*1.001,par("usr")[3]*0+par("usr")[4]*1,title="Magnitude",
legend = c("M<1","M=1","M=3","M=5","M=7"),  
col=1,pch=21,pt.cex=c(0.2,cexsize),
cex =1,bg ="gray90", x.intersp = 2, y.intersp =1.65,bty = "n")
#
legend(par("usr")[2],par("usr")[3]*0.6+par("usr")[4]*0.4,title="Depth",
legend = c("D<10km","10km<=D<20km","20km<=D<40km","40km<=D<80km","80km<=D<150km","150km<=D"),  
cex=1,pch=19,col=c("red","orange","gold4","green","blue","purple"),
pt.cex =1.5,bg ="gray90", x.intersp = 1, y.intersp =1.2,bty = "n")
par(xpd=F)
title(paste0(ymd[1],"～",ymd[length(ymd)]," 地震"))
title("","( データ：気象庁 震源リスト )",line=2)
map.axes(cex.axis=0.8,las=1)
#dev.off()

範囲を狭くする

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
par(family = "TakaoExMincho")
#png("eq20160414_0428_02.png",width=1000,height=800)
map('japan',xlim=c(130.2,132),ylim=c(32.1,33.4),fill=T,col="gray90",mar = c(5,4,4,15))
map('japan',col="gray20",boundary =T,add=T)
points(mag1und$longitude,mag1und$latitude,cex=0.2,bg=adjustcolor(subeq$col,1),pch=21,col=adjustcolor(subeq$col,1))
points(subeq$longitude,subeq$latitude,cex=(subeq$mag-0.2)/2,bg=adjustcolor(subeq$col,0.8),pch=21,col="black")
points(c(130.1894,129.8372,132.3111),c(31.83361,33.51556,33.49083),cex=2,col="black",bg=c(adjustcolor("gray30",0.8),rep(rgb(1,1,1,0.8),2)),pch=24) 
par(xpd=T)
#
cexsize=(c(1,3,5,7)-0.2)/2
legend(par("usr")[2]*1.001,par("usr")[3]*0.05+par("usr")[4]*0.95,title="Magnitude",
legend = c("M<1","M=1","M=3","M=5","M=7"),  
col=1,pch=21,pt.cex=c(0.2,cexsize),
cex =1,bg ="gray90", x.intersp = 2, y.intersp =1.8,bty = "n")
#
legend(par("usr")[2],par("usr")[3]*0.6+par("usr")[4]*0.4,title="Depth",
legend = c("D<10km","10km<=D<20km","20km<=D<40km","40km<=D<80km","80km<=D<150km","150km<=D"),  
cex=1,pch=19,col=c("red","orange","gold4","green","blue","purple"),
pt.cex =1.5,bg ="gray90", x.intersp = 1, y.intersp =1.2,bty = "n")
par(xpd=F)
title(paste0(ymd[1],"～",ymd[length(ymd)]," 地震"))
map.axes(cex.axis=0.8,las=1)
#dev.off()

上の地図の範囲の地震の回数（2016-04-14～2016-04-28)

1
2
3
4
5
6
7
8
subset(eqdata,mag==max(eqdata$mag))
#5312 2016- 4-16 01:25: 5.4  130.7617 32.75333    12 7.3
#
eqdata_k<-subset(eqdata,longitude>=130.2 & longitude<=132 & latitude>=32.1 & latitude<=33.4)
#マグニチュード 1以上の地震のみ
mag_cut<-cut(eqdata_k$mag,breaks=c(1,2,3,4,5,6,6.5,7,7.5),right = FALSE, include.lowest = TRUE)
library(knitr)
kable(table(substring(eqdata_k$time,1,10),mag_cut))
[1,2) [2,3) [3,4) [4,5) [5,6) [6,6.5) [6.5,7) [7,7.5]
2016- 4-14 59 122 57 13 3 0 1 0
2016- 4-15 146 328 76 11 1 1 0 0
2016- 4-16 36 354 299 46 8 0 0 1
2016- 4-17 34 245 68 7 0 0 0 0
2016- 4-18 21 167 43 1 1 0 0 0
2016- 4-19 8 117 42 4 2 0 0 0
2016- 4-20 16 109 40 3 0 0 0 0
2016- 4-21 36 73 17 2 0 0 0 0
2016- 4-22 31 80 19 0 0 0 0 0
2016- 4-23 21 58 9 0 0 0 0 0
2016- 4-24 30 51 10 2 0 0 0 0
2016- 4-25 35 39 6 1 0 0 0 0
2016- 4-26 23 46 7 0 0 0 0 0
2016- 4-27 21 51 10 0 0 0 0 0
2016- 4-28 19 52 8 2 0 0 0 0

動画作成

1分ごと(2016-04-16 01:00::2016-04-16 01:59)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#library(mapdata)
library(ggplot2)
#devtools::install_github("dgrtwo/gganimate")
library(gganimate)
library(xts)
#
# 1分ごと
lon1=c(129.0449,132.5606)
lat1=c(31.003,33.96793)
# 1分ごとにプロットする準備
eq$period<-substring(eq$time,1,16)
#プロットする期間を抽出する場合
#時間単位で取り出す場合 xts(eq[,-1], as.POSIXct(eq$time)) とする
eq.xts<-xts(eq[,-1], as.POSIXct(eq$time))
#
subeq<-data.frame(coredata(eq.xts["2016-04-16 01:00::2016-04-16 01:59"]),stringsAsFactors = F)
#str(subeq)
#文字列になってしまったところをnumericに変換
subeq<-transform(subeq,longitude = as.numeric(longitude),latitude= as.numeric(latitude),depth= as.numeric(depth),mag= as.numeric(mag))
#マグニチュードの昇順に並べ替え
sortlist <- order(subeq$mag)
subeq<- subeq[sortlist,]
#
kyusyu <- map("japan",xlim=c(lon1[1],lon1[2]),ylim=c(lat1[1],lat1[2]), plot = FALSE,fill=T)
#
p<-ggplot(kyusyu, aes(long, lat)) +
  geom_polygon(aes(group = group),fill=rgb(0.8,0.2,0,alpha=0.7))+
  coord_cartesian(xlim=c(lon1[1],lon1[2]),ylim=c(lat1[1],lat1[2])) +
  geom_point(data=subeq, aes(x=longitude, y=latitude,fill=depth_rank,size=mag, frame =period),alpha=0.8,shape=21,color="black")+ 
#プロットする色を指定したい場合
  scale_fill_manual(values=c("red","orange","gold4","green","blue","purple"), 
  breaks = c(" 10km 未満"," 10km以上20km未満"," 20km以上40km未満"," 40km以上80km未満"," 80km以上150km未満","150km以上"))+
#プロットする大きさの範囲を指定したい場合
  scale_size(range = c(0.1,10))+
  guides(fill=FALSE , size=FALSE)
#p
#
gg_animate(p, "kumamoto2016041601.gif", interval = 1.5)

1日ごと(2016-04-14::2016-04-28)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
lon1=c(129.0449,132.5606)
lat1=c(31.003,33.96793)
# 1日ごとにプロットする準備
eq$period<-substring(eq$time,1,10)
#プロットする期間を抽出する場合
eq.xts<-xts(eq[,-1], as.POSIXct(eq$time))
#
subeq<-data.frame(coredata(eq.xts["2016-04-14::2016-04-28"]),stringsAsFactors = F)
#str(subeq)
#文字列になってしまったところをnumericに変換
subeq<-transform(subeq,longitude = as.numeric(longitude),latitude= as.numeric(latitude),depth= as.numeric(depth),mag= as.numeric(mag))
#マグニチュードの昇順に並べ替え
sortlist <- order(subeq$mag)
subeq<- subeq[sortlist,]
#
kyusyu <- map("japan",xlim=c(lon1[1],lon1[2]),ylim=c(lat1[1],lat1[2]), plot = FALSE,fill=T)
#
p<-ggplot(kyusyu, aes(long, lat)) +
  geom_polygon(aes(group = group),fill=rgb(0.8,0.2,0,alpha=0.7))+
  coord_cartesian(xlim=c(lon1[1],lon1[2]),ylim=c(lat1[1],lat1[2])) +
  geom_point(data=subeq, aes(x=longitude, y=latitude,fill=depth_rank,size=mag, frame =period),alpha=0.8,shape=21,color="black")+ 
#プロットする色を指定したい場合
  scale_fill_manual(values=c("red","orange","gold4","green","blue","purple"), 
  breaks = c(" 10km 未満"," 10km以上20km未満"," 20km以上40km未満"," 40km以上80km未満"," 80km以上150km未満","150km以上"))+
#プロットする大きさの範囲を指定したい場合
  scale_size(range = c(0.1,10))+
  guides(fill=FALSE , size=FALSE)
#p
#
gg_animate(p, "kumamoto20160414_28.gif", interval = 1.5)