地震

平成28年10月21日鳥取県中部の地震01

mapdata パッケージ

震源リストに２０１６年１０月２１日のデータがあがってないので２０１５年から２０１６年１０月２０日までの震央分布図を作成した。

前の記事
気象庁地震カタログ01
気象庁震源リスト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
#データをダウンロード、解凍
date<-"h2015"
#date<-"h201601"
#date<-"h201602"
#date<-"h201603"
temp <- tempfile()
download.file(paste0("http://www.data.jma.go.jp/svd/eqev/data/bulletin/data/hypo/",date,".zip"),temp)
x<- readLines(unz(temp,date))
unlink(temp)
#必要な項目を取り出し整形
#オリジンタイム
#02-05	I4	西暦	オリジンタイムの西暦
#06-07	I2	月	オリジンタイムの月
#08-09	l2	日	オリジンタイムの日
#10-11	I2	時	オリジンタイムの時
#12-13	I2	分	オリジンタイムの分
#14-17	F4.2	秒	オリジンタイムの秒
time<-paste0(substr(x,2,5),"-",substr(x,6,7),"-",substr(x,8,9)," ",substr(x,10,11),":",substr(x,12,13),":",substr(x,14,15),".",substr(x,16,17))
#震央の緯度、震央の経度、震源の深さ(km)
#震央の緯度
#22-24	I3	緯度(度)	震央の緯度(度)
#25-28	F4.2	緯度(分)	震央の緯度(分)
#震央の経度
#33-36	I4	経度(度)	震央の経度(度)
#37-40	F4.2	経度(分)	震央の経度(分)
#震源の深さ(km)
#45-49	F5.2	深さ(km)	深さフリーの条件で計算した時の震源の深さ(km)
#南緯、西経の場合も考慮に入れると
latitude = as.numeric(paste0(substr(x,22,22),as.character(as.numeric(substr(x,23,24))+as.numeric(substr(x,25,26))/60+as.numeric(substr(x,27,28))/6000)))
longitude = as.numeric(paste0(substr(x,33,33),as.character(as.numeric(substr(x,34,36))+as.numeric(substr(x,37,38))/60+as.numeric(substr(x,39,40))/6000)))
depth<-as.numeric(substr(x,45,49))/100
#マグニチュード
#53-54	F2.1	マグニチュード1	
#気象庁マグニチュード、気象庁ＣＭＴのモーメントマグニチュードまたはＵＳＧＳ等が計算した実体波マグニチュード
#0未満の場合は以下のように表記する
#-0.1: -1   -0.9: -9   -1.0: A0
#-1.9: A9   -2.0: B0   -3.0: C0
mag<-as.numeric(gsub(" . ","NA",
	paste0(gsub("C", "-3",
		gsub("B", "-2",
			gsub("A", "-1",
				gsub("-", "-0",substr(x,53,53))))),
					".",substr(x,54,54))))
#96	A1	震源決定フラグ	
#K: 気象庁震源
#S: 参考震源
#Kは決定精度が良いもので、防災機関へは原則としてこれのみを表示した分布図のみを提供しており、
#Sは決定精度が悪いもので、必要に応じて参考にするためのもの
flag<-substr(x,96,96)
#データフレーム作成(決定精度が良い K のみ)
eq<-subset(data.frame(time,longitude,latitude,depth,mag,flag),flag=="K")
#データの確認
summary(eq)
eq<-eq[,-6]
#マグニチュードのデータがあり、かつ日本周辺のデータのみ抽出
eqdata<-subset(eq,longitude>=119 & longitude<=155 & latitude>=19 & latitude<=51 & mag>=-4)
#データの確認
head(eqdata) ; tail(eqdata)
#保存する場合はコメントアウトを外す。
#save("eqdata", file=paste0("eqdata",date,".RData"))

震源リスト

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
eqdata<-read.table(text="",col.names=c("time", "longitude", "latitude", "depth", "mag"))
ymd<-seq(as.Date("2016-04-01"),as.Date("2016-10-20"), 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="eqdata20160401_1020.RData")
#load("eqdata20160401_1020.RData")
#
load("eqdatah2015.RData")
eq<-eqdata
load("eqdatah201601.RData")
eq<-rbind(eq,eqdata)
load("eqdatah201602.RData")
eq<-rbind(eq,eqdata)
load("eqdatah201603.RData")
eq<-rbind(eq,eqdata)
load("eqdata20160401_1020.RData")
eq<-rbind(eq,eqdata)
#
eqdata<-eq
rm(eq)
summary(eqdata)

震央分布図作成

準備

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)
#
#震源の深さ
# 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
25
26
ymd<-"2015-01-01 ～ 2016-10-20"
#
par(family = "TakaoExMincho")
#
#png("eq2015_20161020.png",width=1000,height=800)
map('japan',xlim=c(132,136),ylim=c(34,36),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,col=adjustcolor(subeq$col,1),pch=19)
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+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," 地震"))
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
#png("eq2015_20161020_02.png",width=1000,height=800)
map('japan',xlim=c(133,134.7),ylim=c(35,35.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,col=adjustcolor(subeq$col,1),pch=19)
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.0001,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," 地震"))
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
#png("eq2015_20161020_03.png",width=1000,height=800)
map('japan',xlim=c(133.6,134.1),ylim=c(35.25,35.55),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,col=adjustcolor(subeq$col,1),pch=19)
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,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," 地震"))
title("","( データ：気象庁 地震月報(カタログ編) & 震源リスト )",line=2)
map.axes(cex.axis=0.8,las=1)
#dev.off()