R

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

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

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

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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
library(mapdata)
#
eqdata<-read.table(text="",col.names=c("time", "longitude", "latitude", "depth", "mag"))
#ymd<-seq(as.Date("2016-01-01"),as.Date("2016-04-19"), by="1 day")
ymd<-seq(as.Date("2016-04-14"),as.Date("2016-04-19"), 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)
#
#ymd<-seq(as.Date("2016-01-01"),as.Date("2016-04-13"), by="1 day")
#
par(family = "TakaoExMincho")
#
#png("eq20160414_0419_01.png",width=1000,height=800)
map('japan',fill=T,col="gray90")
map('japan',col="gray50",boundary =F,add=T)
points(subeq$longitude,subeq$latitude,cex=(subeq$mag-0.2)/2,bg=adjustcolor(subeq$col,1),pch=21,col="black")
par(xpd=T)
#
cexsize=(c(1,3,5,7,9)-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=3","M=5","M=7","M=9"),  
col=1,pch=21,pt.cex=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)]," 地震(Magnitude >=1)"))
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
#png("eq20160414_0419_02.png",width=1000,height=800)
map('japan',xlim=c(128.5,139.5),ylim=c(30,36),fill=T,col="gray90",mar = c(5,4,4,15))
points(mag1und$longitude,mag1und$latitude,cex=0.2,bg=adjustcolor(subeq$col,1),pch=21,col="black")
points(subeq$longitude,subeq$latitude,cex=(subeq$mag-0.2)/2,bg=adjustcolor(subeq$col,1),pch=21,col="black")
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)]," 地震"))
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
#png("eq20160414_0419_03.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="gray30",boundary =T,add=T)
points(mag1und$longitude,mag1und$latitude,cex=0.2,bg=adjustcolor(subeq$col,1),pch=21,col="black")
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
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#library(mapdata)
library(ggplot2)
#devtools::install_github("dgrtwo/gganimate")
library(gganimate)
library(xts)
#
lon1=c(129.0449,132.5606)
lat1=c(31.003,33.96793)
#時間ごとにプロットする準備
eq$period<-substring(eq$time,1,13)
#プロットする期間を抽出する場合
#時間単位で取り出す場合 xts(eq[,-1], as.POSIXct(eq$time)) とする
eq.xts<-xts(eq[,-1], as.POSIXct(eq$time))
#
subeq<-data.frame(coredata(eq.xts["2016-04-14 21::2016-04-19 23"]),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, "kumamoto2016.gif", interval = 1.5)