R

気象庁地震カタログ02

mapdata パッケージ

(参考)

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

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

1998年から2014年までのデータを整形して保存しておく。

鳥取県周辺のデータのみ抽出

1
2
3
4
5
6
7
8
9
10
11
eq<-data.frame(matrix(rep(NA,6),nrow=1))[-1,]
names(eq)<-c("time", "longitude", "latitude", "depth", "mag", "flag")
year<-seq(1998,2014,1)
for (i in 1:length(year)) {
load(paste0("eqdatah",year[i],".RData"))
#鳥取県周辺のデータのみ抽出
eqdata<-subset(eqdata,longitude>=132.5 & longitude<=135 & latitude>=34.5 & latitude<=35.9 & mag>=-4)
eq<-rbind(eq,eqdata)
}
eqdata<-eq
#save("eqdata", file="tottori1998_2014.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
#load("tottori1998_2014.RData")
#M>=3、微小地震、極微小地震で色を変えるので分類する
#極微小地震
eqdata1<-subset(eqdata,mag<1)
#微小地震
eqdata2<-subset(eqdata,mag>=1 & mag<3)
#M>=3(マグニチュードの昇順に並べ替え)
eqdata3<-subset(eqdata,mag>=3)
sortlist <- order(eqdata3$mag)
eqdata3<- eqdata3[sortlist,]
#
head(eqdata1) ; tail(eqdata1)
head(eqdata2) ; tail(eqdata2)
head(eqdata3) ; tail(eqdata3)
date<-"1998_2014"
#
depth_class <- factor(cut(eqdata$depth,breaks =c(0,10,20,30,max(eqdata$depth)),
include.lowest=T, labels=c("0～10km","10～20km","20～30km","30km～")))
mag_class<-factor(cut(eqdata$mag,breaks =c(-4,1,3,max(eqdata$mag)),right=F),
labels=c("[ ,1)","[1,3)","[3, )"))
table(depth_class,mag_class)
#png("tottorieq01.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()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(mapdata)
#png("tottorieq02.png",width=1000,height=800)
alpha=0.2
map('japan',xlim=c(132.5,135),ylim=c(34.5,35.9),fill=T,col="gray90")
map('japan',col="gray60",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,"地震(微小地震、極微小地震を含む)"))
rect(xleft=133,ybottom=35.1,xright=133.5,ytop=35.5,lwd=0.5)
box()
#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
#四角で囲まれた部分
subeqdata<-subset(eqdata,longitude>=133 & longitude<=133.5 & latitude>=35.1 & latitude<=35.5)
#M>=3、微小地震、極微小地震で色を買えるので分類する
#極微小地震
subeqdata1<-subset(subeqdata,mag<1)
#微小地震
subeqdata2<-subset(subeqdata,mag>=1 & mag<3)
#M>=3(マグニチュードの昇順に並べ替え)
subeqdata3<-subset(subeqdata,mag>=3)
sortlist <- order(subeqdata3$mag)
subeqdata3<- subeqdata3[sortlist,]
#
#四角で囲まれた部分の地震の頻度時系列
#png("tottorieq05.png",width=1000,height=800)
par(mfrow=c(3,1),mar=c(3,4,3,2))
hindo<-data.frame(table(as.Date(subeqdata1$time)))
dummy<-data.frame(Var1=seq(as.Date("1998-01-01"),as.Date("2014-12-31"),by="days"))
eq_hindo<-merge(dummy, hindo, by="Var1",all=T) 
eq_hindo[is.na(eq_hindo)]<-0
plot(eq_hindo,type="h",main="極微小地震")
#
hindo<-data.frame(table(as.Date(subeqdata2$time)))
dummy<-data.frame(Var1=seq(as.Date("1998-01-01"),as.Date("2014-12-31"),by="days"))
eq_hindo<-merge(dummy, hindo, by="Var1",all=T) 
eq_hindo[is.na(eq_hindo)]<-0
plot(eq_hindo,type="h",main="微小地震")
#
hindo<-data.frame(table(as.Date(subeqdata3$time)))
dummy<-data.frame(Var1=seq(as.Date("1998-01-01"),as.Date("2014-12-31"),by="days"))
eq_hindo<-merge(dummy, hindo, by="Var1",all=T) 
eq_hindo[is.na(eq_hindo)]<-0
plot(eq_hindo,type="h",main="M>=3")
par(mfrow=c(1,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
#
depth_class <- factor(cut(subeqdata$depth,breaks =c(0,10,20,30,max(subeqdata$depth)),
include.lowest=T, labels=c("0～10km","10～20km","20～30km","30km～")))
mag_class<-factor(cut(subeqdata$mag,breaks =c(-4,1,3,max(subeqdata$mag)),right=F),
labels=c("[ ,1)","[1,3)","[3, )"))
table(depth_class,mag_class)
#
#png("tottorieq03.png",width=1000,height=800)
#透過率 alphaを指定する
alpha=0.2
par(mfrow=c(1,2))
plot(subeqdata1$longitude,-subeqdata1$depth,cex=(subeqdata1$mag+4)/20,col=rgb(1,0,0,alpha=alpha),las=1,pch=20,
xlim=range(subeqdata$longitude),ylim=range(-subeqdata$depth),xlab="longitude",ylab="depth")
points(subeqdata2$longitude,-subeqdata2$depth,cex=(subeqdata2$mag+4)/20,col=rgb(0,1,0,alpha=alpha),pch=20)
points(subeqdata3$longitude,-subeqdata3$depth,cex=(subeqdata3$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)
rect(xleft=133.23,ybottom=-37,xright=133.29,ytop=-25,lwd=0.5)
title(paste("震源：経度・深さ (",date,")"),cex.main=2)
plot(subeqdata1$latitude,-subeqdata1$depth,cex=(subeqdata1$mag+4)/20,col=rgb(1,0,0,alpha=alpha),las=1,pch=20,
xlim=range(subeqdata$latitude),ylim=range(-subeqdata$depth),xlab="latitude",ylab="depth")
points(subeqdata2$latitude,-subeqdata2$depth,cex=(subeqdata2$mag+4)/20,col=rgb(0,1,0,alpha=alpha),pch=20)
points(subeqdata3$latitude,-subeqdata3$depth,cex=(subeqdata3$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)
rect(xleft=35.23,ybottom=-37,xright=35.31,ytop=-25,lwd=0.5)
title(paste("震源：緯度・深さ (",date,")"),cex.main=2)
par(mfrow=c(1,1))
#dev.off()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#png("tottorieq04.png",width=1000,height=800)
alpha=0.3
map('japan',xlim=c(132.5,134),ylim=c(35,35.6),fill=T,col="gray90")
map('japan',col="gray60",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,"地震(微小地震、極微小地震を含む)"))
rect(xleft=133,ybottom=35.1,xright=133.5,ytop=35.5,lwd=0.8,lty=2)
rect(xleft=133.23,ybottom=35.23,xright=133.29,ytop=35.31,lwd=2)
box()
#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
38
39
40
41
42
43
44
45
46
47
#四角で囲まれた部分(2)
subeqdata<-subset(eqdata,longitude>=133.23 & longitude<=133.29 & latitude>=35.23 & latitude<=35.31)
#M>=3、微小地震、極微小地震で色を買えるので分類する
#極微小地震
subeqdata1<-subset(subeqdata,mag<1)
#微小地震
subeqdata2<-subset(subeqdata,mag>=1 & mag<3)
#M>=3(マグニチュードの昇順に並べ替え)
subeqdata3<-subset(subeqdata,mag>=3)
sortlist <- order(subeqdata3$mag)
subeqdata3<- subeqdata3[sortlist,]
#
head(subeqdata1) ; tail(subeqdata1)
head(subeqdata2) ; tail(subeqdata2)
head(subeqdata3) ; tail(subeqdata3)
#
depth_class <- factor(cut(subeqdata$depth,breaks =c(0,10,20,30,max(subeqdata$depth)),
include.lowest=T, labels=c("0～10km","10～20km","20～30km","30km～")))
mag_class<-factor(cut(subeqdata$mag,breaks =c(-4,1,3,max(subeqdata$mag)),right=F),
labels=c("[ ,1)","[1,3)","[3, )"))
table(depth_class,mag_class)
#png("tottorieq06.png",width=1000,height=800)
#透過率 alphaを指定する
alpha=0.6
par(mfrow=c(1,2),mar=c(15,5,3,2))
plot(subeqdata1$longitude,-subeqdata1$depth,cex=(subeqdata1$mag+4)/20,col=rgb(1,0,0,alpha=alpha),las=1,pch=20,
xlim=range(subeqdata$longitude),ylim=range(-subeqdata$depth),xlab="longitude",ylab="depth")
points(subeqdata2$longitude,-subeqdata2$depth,cex=(subeqdata2$mag+4)/20,col=rgb(0,1,0,alpha=alpha),pch=20)
points(subeqdata3$longitude,-subeqdata3$depth,cex=(subeqdata3$mag+4)/20,col=rgb(0,0,1,alpha=alpha),pch=20)
par(xpd=T)
legend(x=par()$usr[1]+0.01,y=par()$usr[3]-7,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,horiz=T)
par(xpd=F)
title(paste("震源：経度・深さ (",date,")"),cex.main=2)
plot(subeqdata1$latitude,-subeqdata1$depth,cex=(subeqdata1$mag+4)/20,col=rgb(1,0,0,alpha=alpha),las=1,pch=20,
xlim=range(subeqdata$latitude),ylim=range(-subeqdata$depth),xlab="latitude",ylab="depth")
points(subeqdata2$latitude,-subeqdata2$depth,cex=(subeqdata2$mag+4)/20,col=rgb(0,1,0,alpha=alpha),pch=20)
points(subeqdata3$latitude,-subeqdata3$depth,cex=(subeqdata3$mag+4)/20,col=rgb(0,0,1,alpha=alpha),pch=20)
par(xpd=T)
legend(x=par()$usr[1]+0.01,y=par()$usr[3]-7,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,horiz=T)
par(xpd=F)
title(paste("震源：緯度・深さ (",date,")"),cex.main=2)
par(mfrow=c(1,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
38
39
40
41
#そのうち20kmより深いところで起こった地震
subeqdata<-subset(eqdata,longitude>=133.23 & longitude<=133.29 & latitude>=35.23 & latitude<=35.31 & depth>20)
#M>=3、微小地震、極微小地震で色を買えるので分類する
#極微小地震
subeqdata1<-subset(subeqdata,mag<1)
#微小地震
subeqdata2<-subset(subeqdata,mag>=1 & mag<3)
#M>=3(マグニチュードの昇順に並べ替え)
subeqdata3<-subset(subeqdata,mag>=3)
sortlist <- order(subeqdata3$mag)
subeqdata3<- subeqdata3[sortlist,]
#
head(subeqdata1) ; tail(subeqdata1)
head(subeqdata2) ; tail(subeqdata2)
head(subeqdata3) ; tail(subeqdata3)
#png("tottorieq08.png",width=1000,height=800)
#透過率 alphaを指定する
alpha=0.6
par(mfrow=c(1,2),mar=c(15,5,3,2))
plot(subeqdata1$longitude,-subeqdata1$depth,cex=(subeqdata1$mag+4)/20,col=rgb(1,0,0,alpha=alpha),las=1,pch=20,
xlim=range(subeqdata$longitude),ylim=range(-subeqdata$depth),xlab="longitude",ylab="depth")
points(subeqdata2$longitude,-subeqdata2$depth,cex=(subeqdata2$mag+4)/20,col=rgb(0,1,0,alpha=alpha),pch=20)
points(subeqdata3$longitude,-subeqdata3$depth,cex=(subeqdata3$mag+4)/20,col=rgb(0,0,1,alpha=alpha),pch=20)
par(xpd=T)
legend(x=par()$usr[1]+0.01,y=par()$usr[3]-3,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,horiz=T)
par(xpd=F)
title(paste("震源：経度・深さ (",date,")"),cex.main=2)
plot(subeqdata1$latitude,-subeqdata1$depth,cex=(subeqdata1$mag+4)/20,col=rgb(1,0,0,alpha=alpha),las=1,pch=20,
xlim=range(subeqdata$latitude),ylim=range(-subeqdata$depth),xlab="latitude",ylab="depth")
points(subeqdata2$latitude,-subeqdata2$depth,cex=(subeqdata2$mag+4)/20,col=rgb(0,1,0,alpha=alpha),pch=20)
points(subeqdata3$latitude,-subeqdata3$depth,cex=(subeqdata3$mag+4)/20,col=rgb(0,0,1,alpha=alpha),pch=20)
par(xpd=T)
legend(x=par()$usr[1]+0.01,y=par()$usr[3]-3,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,horiz=T)
par(xpd=F)
title(paste("震源：緯度・深さ (",date,")"),cex.main=2)
par(mfrow=c(1,1))
#dev.off()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#png("tottorieq07.png",width=1000,height=800)
par(mfrow=c(2,1),mar=c(3,4,3,2))
hindo<-data.frame(table(as.Date(subeqdata1$time)))
dummy<-data.frame(Var1=seq(as.Date("1998-01-01"),as.Date("2014-12-31"),by="days"))
eq_hindo<-merge(dummy, hindo, by="Var1",all=T) 
eq_hindo[is.na(eq_hindo)]<-0
plot(eq_hindo,type="h",main="極微小地震")
#
hindo<-data.frame(table(as.Date(subeqdata2$time)))
dummy<-data.frame(Var1=seq(as.Date("1998-01-01"),as.Date("2014-12-31"),by="days"))
eq_hindo<-merge(dummy, hindo, by="Var1",all=T) 
eq_hindo[is.na(eq_hindo)]<-0
plot(eq_hindo,type="h",main="微小地震")
#dev.off()