R

気象庁地震カタログ(地震の回数）

ggmap , xts パッケージ

記事気象庁地震カタログ(鳥取県：震源の深さ）の
「10月1日から12月31日までの地震の回数」で1日ごとの地震の回数を求めるのに「table関数」を使いましたが、今回はxtsパッケージの関数を使ってみます。

地震の回数を求める範囲(薄い赤色で塗りつぶした部分）

２０１５年（前回の記事の地図を作成したコードを少し変更）

１９９８年から２０１４年

地震の回数

2015年10月1日から12月31日まで（1日ごと）

  • xts::apply.daily 関数

２０１５年（1日ごと）

  • xts::apply.daily 関数

2015年10月17日（1時間ごと）

  • xts::endpoints 関数　xts::period.apply 関数

1998年から2015年（1月ごと）

  • xts::apply.monthly 関数

コード

地図を描くコードの主な変更点

プロットする前に地図の範囲のデータを抽出した（サイズを適切にするため）

attr(yurihama,”bb”)でget_map関数でゲットした地図の範囲を得ることができる。
1
2
3
4
5
eq<-subset(eq,attr(yurihama,"bb")$ll.lon<=longitude & longitude<=attr(yurihama,"bb")$ur.lon & 
attr(yurihama,"bb")$ll.lat<=latitude &  latitude<=attr(yurihama,"bb")$ur.lat)
#
eqdata4<-subset(eqdata4,attr(yurihama,"bb")$ll.lon<=longitude & longitude<=attr(yurihama,"bb")$ur.lon & 
attr(yurihama,"bb")$ll.lat<=latitude &  latitude<=attr(yurihama,"bb")$ur.lat)

magの大きさに合わせてsizeを変更

aes(x=longitude, y=latitude,color=Magnitude) を aes(x=longitude, y=latitude,color=Magnitude,size=mag) に

コード（地震の回数集計）

  • グラフにするとき xlim で期間を指定すること
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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
#### 薄い赤色で塗りつぶした範囲の地震の頻度
#
##### データを読み込み、グラフ化する範囲のデータを抽出。xtsクラスへ
#
load("2015.RData")
eq2015<-eqdata
load("tottori1998_2014.RData")
#データをつなげる
eqdata<-rbind(eqdata[,-6],eq2015)
#重複した行を取り除く
eqdata<-unique(eqdata)
#グラフ化する範囲
lon1<-c(133.85,133.975) ; lat1<-c(35.39,35.49)
#
#期間を区切ったデータを扱いやすくするためにxtsオブジェクトにする
library(xts)
#年月日まで
#eq.xts<-as.xts(read.zoo(eqdata))
#年月日時分秒（小数点以下切捨て）まで
eq.xts<-xts(eqdata[,2:5], strptime(eqdata$time, "%Y-%m-%d %H:%M:%S"))
#
#指定した範囲のデータを抽出。極微小地震(1未満)、微小地震（1以上３未満）,M>=3 にクラス分け
subeq1.xts<-subset(eq.xts,lon1[1]<=longitude & longitude<=lon1[2] & lat1[1]<=latitude & latitude<=lat1[2] & mag<1)
subeq2.xts<-subset(eq.xts,lon1[1]<=longitude & longitude<=lon1[2] & lat1[1]<=latitude & latitude<=lat1[2] & mag>=1 & mag<3)
subeq3.xts<-subset(eq.xts,lon1[1]<=longitude & longitude<=lon1[2] & lat1[1]<=latitude & latitude<=lat1[2] & mag>=3)
#
#### データを集計してグラフにする
#
###### 2015年
#
date<-"2015"
#
count1=apply.daily(subeq1.xts[date],nrow)
count2=apply.daily(subeq2.xts[date],nrow)
count3=apply.daily(subeq3.xts[date],nrow)
#
#png("eqgg24.png",width=800,height=800)
lwd=3
par(mfrow=c(3,1),mar=c(3,4,3,2))
c1<-data.frame(time=as.Date(substring(index(count1),1,10)),frequency=coredata(count1))
plot(frequency~time, data=c1,type="h",main=paste0("極微小地震 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.Date(c("2015-01-01","2015-12-31")))
#
c2<-data.frame(time=as.Date(substring(index(count2),1,10)),frequency=coredata(count2))
plot(frequency~time, data=c2,type="h",main=paste0("微小地震 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.Date(c("2015-01-01","2015-12-31")))
#
c3<-data.frame(time=as.Date(substring(index(count3),1,10)),frequency=coredata(count3))
plot(frequency~time, data=c3,type="h",main=paste0("M>=3 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.Date(c("2015-01-01","2015-12-31")),ylim=c(0,5))
#dev.off()
#
######  2015-10-01 から2015-12-31
#
start<-"2015-10-01"  
end<-"2015-12-31"  
date<-paste0(start,"::",end)  
#
count1=apply.daily(subeq1.xts[date],nrow)
count2=apply.daily(subeq2.xts[date],nrow)
count3=apply.daily(subeq3.xts[date],nrow)
#
#png("eqgg25.png",width=800,height=800)
lwd=4
par(mfrow=c(3,1),mar=c(3,4,3,2))
c1<-data.frame(time=as.Date(substring(index(count1),1,10)),frequency=coredata(count1))
plot(frequency~time, data=c1,type="h",main=paste0("極微小地震 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.Date(c(start,end))  )
#
c2<-data.frame(time=as.Date(substring(index(count2),1,10)),frequency=coredata(count2))
plot(frequency~time, data=c2,type="h",main=paste0("微小地震 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.Date(c(start,end))  )
#
c3<-data.frame(time=as.Date(substring(index(count3),1,10)),frequency=coredata(count3))
plot(frequency~time, data=c3,type="h",main=paste0("M>=3 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.Date(c(start,end))  ,ylim=c(0,5))
#dev.off()
#
###### 2015-10-17 1時間ごとの回数
#
date<-"2015-10-17"
#
t1<-endpoints(subeq1.xts[date],on="hours")
count1<-period.apply(subeq1.xts[date],t1,nrow)
names(count1)<-"frequency"
#count1
#
t2<-endpoints(subeq2.xts[date],on="hours")
count2<-period.apply(subeq2.xts[date],t2,nrow)
names(count2)<-"frequency"
#count2
#
t3<-endpoints(subeq3.xts[date],on="hours")
count3<-period.apply(subeq3.xts[date],t3,nrow)
names(count3)<-"frequency"
#count3
#
#png("eqgg26.png",width=800,height=800)
lwd=7
par(mfrow=c(3,1),mar=c(3,4,3,2))
c1<-data.frame(time=as.POSIXct(paste0(substring(index(count1),1,13),":00:00")),coredata(count1))
plot(frequency~time, data=c1,type="h",main=paste0("極微小地震 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.POSIXct(c(paste0(date," 00:00:00"),paste0(date," 23:59:59"))) )
#
c2<-data.frame(time=as.POSIXct(paste0(substring(index(count2),1,13),":00:00")),coredata(count2))
plot(frequency~time, data=c2,type="h",main=paste0("微小地震 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.POSIXct(c(paste0(date," 00:00:00"),paste0(date," 23:59:59"))) )
#
c3<-data.frame(time=as.POSIXct(paste0(substring(index(count3),1,13),":00:00")),coredata(count3))
plot(frequency~time, data=c3,type="h",main=paste0("M>=3 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.POSIXct(c(paste0(date," 00:00:00"),paste0(date," 23:59:59"))) ,ylim=c(0,2))
#dev.off()
#
#subeq1.xts[date]
#subeq2.xts[date]
#subeq3.xts[date]
#
###### 1998_2015 の月ごとの地震の回数
#
start<-"1998-01-01"  
end<-"2015-12-31"  
date<-paste0(start,"::",end)  
#
count1=apply.monthly(subeq1.xts[date],nrow)
count2=apply.monthly(subeq2.xts[date],nrow)
count3=apply.monthly(subeq3.xts[date],nrow)
#
count3
#                   [,1]
#2001-11-20 23:13:17    1
#2015-10-19 11:14:03    4
#2015-12-14 15:01:37    1
#
#このままだと月ごとの日付がばらばら。
#
#png("eqgg27.png",width=800,height=800)
lwd=2.5
par(mfrow=c(3,1),mar=c(3,4,3,2))
#日付を01に揃える
c1<-data.frame(time=as.Date(paste0(substring(index(count1),1,7),"-01")),frequency=coredata(count1))
plot(frequency~time, data=c1,type="h",main=paste0("極微小地震 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.Date(c(start,end)) )
#
c2<-data.frame(time=as.Date(paste0(substring(index(count2),1,7),"-01")),frequency=coredata(count2))
plot(frequency~time, data=c2,type="h",main=paste0("微小地震 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.Date(c(start,end)) )
#
c3<-data.frame(time=as.Date(paste0(substring(index(count3),1,7),"-01")),frequency=coredata(count3))
plot(frequency~time, data=c3,type="h",main=paste0("M>=3 (",date,")"),lend=1,lwd=lwd,col="gray30",
xlim=as.Date(c(start,end)) ,ylim=c(0,5))
#dev.off()
#
subeq3.xts[date]
#                    longitude latitude depth mag
#2001-11-20 23:13:17  133.8995 35.41633  11.4 3.0
#2015-10-17 17:53:54  133.9117 35.43667   8.0 3.8
#2015-10-18 08:30:35  133.9050 35.43667   8.0 4.2
#2015-10-18 08:36:52  133.9117 35.43833   8.0 4.3
#2015-10-19 11:14:03  133.8967 35.43333   8.0 3.0
#2015-12-14 15:01:37  133.9100 35.44500   8.0 4.2
メモ

start<-“2015-10-01”
end<-“2015-12-31”
date<-paste0(start,”::”,end)
plot(・・・・・, xlim=as.Date(c(start,end)) )

date<-“2015-10-17”
plot(・・・・・, xlim=as.POSIXct(c(paste0(date,” 00:00:00”),paste0(date,” 23:59:59”))) )

start<-“1998-01-01”
end<-“2015-12-31”
date<-paste0(start,”::”,end)
plot(・・・・・, xlim=as.Date(c(start,end)) )