R

気象庁地震カタログ01

mapdata パッケージ

(参考)

日本の地震カタログ（1923年から現在）の解説
地震月報(カタログ編)
震源データ

震源レコードフォーマット
「オリジンタイム、震央の緯度、震央の経度、震源の深さ(km)、マグニチュード1、震源決定フラグ」を取り出し整形。

例として2014年のデータを使う

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
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
#データをダウンロード、解凍
date<-"h2014"
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<-NULL
#system.time(
#for (i in 1:length(x)){
#if (substr(x[i],53,53)=="-"){
#mag[i]=as.numeric(paste0("-0.",substr(x[i],54,54)))
#}else if (substr(x[i],53,53)=="A"){
#mag[i]=as.numeric(paste0("-1.",substr(x[i],54,54)))
#}else if (substr(x[i],53,53)=="B"){
#mag[i]=as.numeric(paste0("-2.",substr(x[i],54,54)))
#}else if (substr(x[i],53,53)=="C"){
#mag[i]=as.numeric(paste0("-3.",substr(x[i],54,54)))
#}else if (substr(x[i],53,53)==" "){
#mag[i]=NA
#}else {
#mag[i]=as.numeric(paste0(substr(x[i],53,53),".",substr(x[i],54,54)))
#}
#}
#)
#上のようにfor文を使うと遅いので以下のコードにした。
#system.time(
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)
#データフレーム作成
eq<-data.frame(time,longitude,latitude,depth,mag,flag)
#データの確認
summary(eq)
#マグニチュードのデータがあり、かつ日本周辺のデータのみ抽出
eqdata<-subset(eq,longitude>=119 & longitude<=155 & latitude>=19 & latitude<=51 & mag>=-4)
#データの確認
head(eqdata) ; tail(eqdata)
#保存する場合はコメントアウトを外す。
#save("eqdata", file=paste0("eqdata",date,".RData"))

2.分類。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#保存したデータセットを取り込む場合は以下のコメントアウトを外す。
#load("eqdatah2014.RData")
#date<-"h2014"
#M>=3、微小地震、極微小地震で色を変えるので分類する
#極微小地震
eqdata1<-subset(eqdata,mag<1)
#微小地震
eqdata2<-subset(eqdata,mag>=1 & mag<3)
#M>=3
eqdata3<-subset(eqdata,mag>=3)
#確認
summary(eqdata1) 
summary(eqdata2) 
summary(eqdata3) 
#マグニチュードと震源の深さ
depth_class <- factor(cut(eqdata$depth,breaks =c(0,60,300,max(eqdata$depth)),
include.lowest=T, labels=c("0～60km","60～300km","300km～")))
mag_class<-factor(cut(eqdata$mag,breaks =c(-4,1,3,max(eqdata$mag)),right=F),
labels=c("[ ,1)","[1,3)","[3, )"))
library(xtable)
print(xtable(table(depth_class,mag_class)),type = "html")





[ ,1) [1,3) [3, )
0～60km 59929 59217 5652
60～300km 2428 8105 727
300km～ 0 69 107

3.作図。

legendの記述にexpression関数を使った。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#png("depth01.png",width=1000,height=800)
#透過率 alphaを指定する
alpha=0.2
par(mfrow=c(1,2))
plot(eqdata1$longitude,-eqdata1$depth,cex=(eqdata1$mag+4)/20,col=rgb(1,0,0,alpha=alpha),las=1,pch=20,
xlim=range(eqdata$longitude),ylim=range(-eqdata$depth),xlab="longitude",ylab="depth")
points(eqdata2$longitude,-eqdata2$depth,cex=(eqdata2$mag+4)/20,col=rgb(0,1,0,alpha=alpha),pch=20)
points(eqdata3$longitude,-eqdata3$depth,cex=(eqdata3$mag+4)/20,col=rgb(0,0,1,alpha=alpha),pch=20)
legend("bottomright",title="Magnitude",legend = c("M<1",expression(paste(1<=M,"<3")),expression(M>=3)),  
cex=1,pch=20,col=c(rgb(1,0,0,alpha=1),rgb(0,1,0,alpha=1),rgb(0,0,1,alpha=1)),
pt.cex =2,bg ="gray", x.intersp = 1, y.intersp =1)
title(paste("震源：経度・深さ (",date,")"),cex.main=2)
plot(eqdata1$latitude,-eqdata1$depth,cex=(eqdata1$mag+4)/20,col=rgb(1,0,0,alpha=alpha),las=1,pch=20,
xlim=range(eqdata$latitude),ylim=range(-eqdata$depth),xlab="latitude",ylab="depth")
points(eqdata2$latitude,-eqdata2$depth,cex=(eqdata2$mag+4)/20,col=rgb(0,1,0,alpha=alpha),pch=20)
points(eqdata3$latitude,-eqdata3$depth,cex=(eqdata3$mag+4)/20,col=rgb(0,0,1,alpha=alpha),pch=20)
legend("bottomright",title="Magnitude",legend = c("M<1",expression(paste(1<=M,"<3")),expression(M>=3)),  
cex=1,pch=20,col=c(rgb(1,0,0,alpha=1),rgb(0,1,0,alpha=1),rgb(0,0,1,alpha=1)),
pt.cex =2,bg ="gray", x.intersp = 1, y.intersp =1)
title(paste("震源：緯度・深さ (",date,")"),cex.main=2)
par(mfrow=c(1,1))
#dev.off()

4.地図上にプロット

1
2
3
4
5
6
7
8
9
10
11
12
13
library(mapdata)
#png("eqmap01.png",width=1000,height=800)
alpha=0.2
map('japan',fill=T,col="gray90")
map('japan',col="gray85",boundary =F,add=T)
points(eqdata1$longitude,eqdata1$latitude,cex=(eqdata1$mag+4)/20,col=rgb(1,0,0,alpha=alpha),pch=20)
points(eqdata2$longitude,eqdata2$latitude,cex=(eqdata2$mag+4)/20,col=rgb(0,1,0,alpha=alpha),pch=20)
points(eqdata3$longitude,eqdata3$latitude,cex=(eqdata3$mag+4)/20,col=rgb(0,0,1,alpha=alpha),pch=20)
legend("topleft",title="Magnitude",legend = c("M<1",expression(paste(1<=M,"<3")),expression(M>=3)),  
cex=1,pch=20,col=c(rgb(1,0,0,alpha=1),rgb(0,1,0,alpha=1),rgb(0,0,1,alpha=1)),
pt.cex =2,bg ="gray", x.intersp = 1, y.intersp =1)
title(paste(date,"地震(微小地震、極微小地震を含む)"))
#dev.off()

1
2
3
4
5
6
7
8
9
10
11
12
#png("eqmap02.png",width=1000,height=800)
alpha=0.4
map('japan',xlim=c(130,139),ylim=c(30,36),fill=T,col="gray90")
map('japan',col="gray85",boundary =F,add=T)
points(eqdata1$longitude,eqdata1$latitude,cex=(eqdata1$mag+4)/20,col=rgb(1,0,0,alpha=alpha),pch=20)
points(eqdata2$longitude,eqdata2$latitude,cex=(eqdata2$mag+4)/20,col=rgb(0,1,0,alpha=alpha),pch=20)
points(eqdata3$longitude,eqdata3$latitude,cex=(eqdata3$mag+4)/20,col=rgb(0,0,1,alpha=alpha),pch=20)
legend("bottomright",title="Magnitude",legend = c("M<1",expression(paste(1<=M,"<3")),expression(M>=3)),  
cex=1,pch=20,col=c(rgb(1,0,0,alpha=1),rgb(0,1,0,alpha=1),rgb(0,0,1,alpha=1)),
pt.cex =2,bg ="gray", x.intersp = 1, y.intersp =1)
title(paste(date,"地震(微小地震、極微小地震を含む)"))
#dev.off()

プロットする大きさを大きく、M>=3はより大きく

1
2
3
4
5
6
7
8
9
10
11
12
13
#png("eqmap03.png",width=1000,height=800)
alpha=0.6
map('japan',xlim=c(132.5,135),ylim=c(34.5,35.9),fill=T,col="gray90")
map('japan',col="gray85",boundary =F,add=T)
points(eqdata1$longitude,eqdata1$latitude,cex=(eqdata1$mag+4)/10,col=rgb(1,0,0,alpha=alpha),pch=20)
points(eqdata2$longitude,eqdata2$latitude,cex=(eqdata2$mag+4)/10,col=rgb(0,1,0,alpha=alpha),pch=20)
points(eqdata3$longitude,eqdata3$latitude,cex=(eqdata3$mag+4)/5,col=rgb(0,0,1,alpha=alpha),pch=20)
legend("topleft",title="Magnitude",legend = c("M<1",expression(paste(1<=M,"<3")),expression(M>=3)),  
cex=1,pch=20,col=c(rgb(1,0,0,alpha=1),rgb(0,1,0,alpha=1),rgb(0,0,1,alpha=1)),
pt.cex =2,bg ="gray", x.intersp = 1, y.intersp =1)
title(paste(date,"地震(微小地震、極微小地震を含む)"))
box()
#dev.off()