R

アメリカ地質調査所 ( 平成28年（2016年）熊本地震 )

RCurl , mapdata , xts パッケージ

(参考)
Search Earthquake Archives
API Documentation - Earthquake Catalog

(注意)時差は修正していません。グラフに日時があってもその点を考慮願います。

マグニチュード4以上

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
95
96
97
98
99
100
101
102
103
104
105
library(RCurl)
library(mapdata)
library(xts)
#
#lon<-c(122.9,155)
#lat<-c(20,46)
lon=c(128,139)
lat=c(30,36)
#取り出すマグニチュードが小さいとダウンロード可のイベント上限超えでエラーになる
#その場合は範囲をしぼる
#日本のデータはマグニチュード4未満はないみたいです
minmagnitude<-4
#
#一括処理(USGS:地震データダウンロード。保存。)
#震源データをダウンロードする年
start<-2016 ; end<-2016
library(RCurl)
for (i in start:end){
year<-as.character(i)
#format=csvとする
url<-paste0("http://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&starttime=",as.numeric(year),"-01-01&endtime=",as.numeric(year)+1,
"-01-01&minmagnitude=",minmagnitude,"&minlatitude=",lat[1],"&maxlatitude=",lat[2],"&minlongitude=",lon[1],"&maxlongitude=",lon[2],"&eventtype=earthquake")
doc <- getURL(url)
filename<-paste0("USGS",year,"_",lon[1],"_",lon[2],"_",lat[1],"_",lat[2],".csv")
out <- file(filename, "w") # 書き込みモードで開く
writeLines(doc,out,sep = "\n") 
close(out)   
}
#
#eq
eq<-read.table(text="",col.names=c("time","latitude","longitude","depth","mag"))
for (i in start:end){
year<-as.character(i)
filename<-paste0("USGS",year,"_",lon[1],"_",lon[2],"_",lat[1],"_",lat[2],".csv")
tmp<-read.csv(filename)
#eqdata <- subset(tmp, select=c(time,latitude,longitude,depth,mag,nst,gap,dmin,rms))
tmp <- subset(tmp, select=c(time,latitude,longitude,depth,mag))
eq<-rbind(eq,tmp)
}
#
#save(eq,file="eq2016.RData") 
#
#USGS:地震データ
eqdata <- subset(eq, select=c(time,longitude,latitude,depth,mag))
#head(eqdata)
eqdata$time<-gsub("T"," ",substring(eqdata[,1],1,19))
eq.xts<-as.xts(read.zoo(eqdata))
ymd<-"2016-04-14::2016-04-16"
eqdata<-data.frame(coredata(eq.xts[ymd]))
#
#震源の深さ
# 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)
subeq<- eq[sortlist,]
##
#
#九州
par(family = "TakaoExMincho")
#png("ks20160414_0416.png",width=1000,height=800)
map('japan',xlim=c(128,133),ylim=c(30,35),fill=T,col="gray95",mar = c(5,4,4,15))
map('japan',col="gray30",boundary =T,add=T)
points(subeq$longitude,subeq$latitude,cex=(subeq$mag-0.2)/2,col=adjustcolor(subeq$col,0.8),pch=19)
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(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=3","M=5","M=7","M=9"),  
col=1,pch=19,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," 地震 ( 地震データ：USGS )"))
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
#熊本県周辺
par(family = "TakaoExMincho")
#png("kumamoto20160414_16.png",width=1000,height=800)
map('japan',xlim=c(130,132),ylim=c(32,33.5),fill=T,col="gray95",mar = c(5,4,4,15))
map('japan',col="gray30",boundary =T,add=T)
points(subeq$longitude,subeq$latitude,cex=(subeq$mag-0.2)/2,col=adjustcolor(subeq$col,0.95),pch=19)
par(xpd=T)
#
cexsize=(c(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=3","M=5","M=7","M=9"),  
col=1,pch=19,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," 地震 ( 地震データ：USGS )"))
map.axes(cex.axis=0.8,las=1)
#dev.off()

比較のために作成。但し、プロットしたデータの期間がかなり違います。

兵庫県南部地震（発生から約一ヶ月間）

岩手宮城内陸地震（発生から約一ヶ月間）

震源のデータ(１９００年 ～）

地震(マグニチュード6以上)

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
library(RCurl)
library(mapdata)
library(xts)
#
lon<-c(122.9,155)
lat<-c(20,46)
minmagnitude<-6
#
#一括処理(USGS:地震データダウンロード。保存。)
#震源データをダウンロードする年
start<-1900 ; end<-2016
library(RCurl)
for (i in start:end){
year<-as.character(i)
#format=csvとする
url<-paste0("http://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&starttime=",as.numeric(year),"-01-01&endtime=",as.numeric(year)+1,
"-01-01&minmagnitude=",minmagnitude,"&minlatitude=",lat[1],"&maxlatitude=",lat[2],"&minlongitude=",lon[1],"&maxlongitude=",lon[2],"&eventtype=earthquake")
doc <- getURL(url)
filename<-paste0("USGS",year,"_",lon[1],"_",lon[2],"_",lat[1],"_",lat[2],".csv")
out <- file(filename, "w") # 書き込みモードで開く
writeLines(doc,out,sep = "\n") 
close(out)   
}
#データをつなぐ
#eq
eq<-read.table(text="",col.names=c("time","latitude","longitude","depth","mag"))
for (i in start:end){
year<-as.character(i)
filename<-paste0("USGS",year,"_",lon[1],"_",lon[2],"_",lat[1],"_",lat[2],".csv")
tmp<-read.csv(filename)
#eqdata <- subset(tmp, select=c(time,latitude,longitude,depth,mag,nst,gap,dmin,rms))
tmp <- subset(tmp, select=c(time,latitude,longitude,depth,mag))
eq<-rbind(eq,tmp)
}
#USGS:地震データ
eqdata <- subset(eq, select=c(time,longitude,latitude,depth,mag))
#head(eqdata)
eqdata$time<-gsub("T"," ",substring(eqdata[,1],1,19))
eq.xts<-as.xts(read.zoo(eqdata))
ymd<-"1900-01-01::2016-12-31"
eqdata<-data.frame(coredata(eq.xts[ymd]))
#
#震源の深さ
# 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)
subeq<- eq[sortlist,]
##
#日本
#png("jpn1900.png",width=1000,height=800)
map('japan',xlim=c(122.9,155),ylim=c(20,46),fill=T,col="gray95",mar = c(5,4,4,15))
points(subeq$longitude,subeq$latitude,cex=(subeq$mag-5)/2,col=adjustcolor(subeq$col,0.8),pch=19)
par(xpd=T)
#
cexsize=(c(6,7,8,9)-5)/2
legend(par("usr")[2]*1.001,par("usr")[3]*0.05+par("usr")[4]*0.95,title="Magnitude",
legend = c("M=6","M=7","M=8","M=9"),  
col=1,pch=19,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("マグニチュード６以上の地震 ( 1900 〜 2016.04.16  地震データ：USGS ) "))
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("ks1900.png",width=1000,height=800)
map('japan',xlim=c(128,133),ylim=c(30,35),fill=T,col="gray95",mar = c(5,4,4,15))
map('japan',col="gray30",boundary =T,add=T)
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) 
points(subeq$longitude,subeq$latitude,cex=(subeq$mag-3)/2,col=adjustcolor(subeq$col,0.8),pch=19)
par(xpd=T)
#
cexsize=(c(6,7,8,9)-3)/2
legend(par("usr")[2]*1.001,par("usr")[3]*0.05+par("usr")[4]*0.95,title="Magnitude",
legend = c("M=6","M=7","M=8","M=9"),  
col=1,pch=19,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("マグニチュード６以上の地震 ( 1900 〜 2016.04.16  地震データ：USGS ) "))
map.axes(cex.axis=0.8,las=1)
#dev.off()