R

気象庁地震カタログ(鳥取県：震源の深さ）

ggmap , grid パッケージ（タイトルと凡例）

(参考)

気象庁 震源リスト(最近9ヶ月、日データ)
これより過去の震源については地震月報(カタログ編)

鳥取県の震源分布

薄い赤の部分と薄い青の部分の震源の深さ

米子市周辺(薄い赤の部分)

鳥取県中部から鳥取市周辺(薄い青の部分)

拡大した震源分布図

コード

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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
library(ggmap)
library(mapproj)
library(grid)
#あらかじめ作成していた震源データを読み込む
load("2015.RData")
eq2015<-eqdata
load("tottori1998_2014.RData")
eqdata<-rbind(eqdata[,-6],eq2015)
#
#
#M>=3、微小地震極微小地震で色を買えるので分類する
#極微小地震
eqdata1<-subset(eqdata,mag<1)
eqdata1$Magnitude<-"1 未満"
#微小地震
eqdata2<-subset(eqdata,mag>=1 & mag<3)
eqdata2$Magnitude<-"1以上3未満"
#M>=3
eqdata3<-subset(eqdata,mag>=3)
eqdata3$Magnitude<-"3 以上"
#
head(eqdata1) ; tail(eqdata1)
head(eqdata2) ; tail(eqdata2)
head(eqdata3) ; tail(eqdata3)
#
#つなげる(ggmapで凡例を付けるため)
eq<-rbind(eqdata1,eqdata2,eqdata3)
#mag>=4(mag>=4の地震はマグニチュードの大きさでプロット)
eqdata4<-subset(eqdata,mag>=4)
#(マグニチュードの昇順に並べ替え)
sortlist <- order(eqdata4$mag)
eqdata4<- eqdata4[sortlist,]
#
#
############### 鳥取県の震源分布図
#
#倉吉を中心にzoom=9を指定する。
center <- geocode('kurayoshi')
#maptype = c(“terrain”, “satellite”, “roadmap”, “hybrid”, “toner”, “watercolor”)
tottori <- get_map(c(center$lon, center$lat), zoom = 9, maptype = 'terrain')
#
lon1<-c(133,133.5) ; lat1<-c(35.1,35.5)
lon2<-c(133.5,134.65) ; lat2<-c(35.2,35.65)
#
#
#png("eqgg10.png",width=800,height=800)
alpha=0.8 ; size=1
ggmap(tottori) + 
geom_point(data=eq, aes(x=longitude, y=latitude,color=Magnitude),alpha=alpha,size=size)+ 
geom_point(data=eqdata4, aes(x=longitude, y=latitude,size=mag),alpha=1,shape=21,color="gray",fill="black")+ 
geom_point(aes(x=134.1844,y=35.4775),shape="Ｘ",color="black",size=7)+ 
geom_rect(aes(xmin=lon1[1], xmax=lon1[2], ymin=lat1[1],ymax=lat1[2]), fill="red", alpha=0.05)+
geom_rect(aes(xmin=lon2[1], xmax=lon2[2], ymin=lat2[1],ymax=lat2[2]), fill="blue", alpha=0.05)+
#プロットする色を指定したい場合
scale_colour_manual(values=c("red","green","blue"), breaks = c("1 未満","1以上3未満","3 以上"))+
#プロットする大きさの範囲を指定したい場合（マグニチュード 4以上）
scale_size(range = c(6,10))+
# 凡例のタイトルの指定
labs(colour="マグニチュード",size="マグニチュード 4以上")+
#凡例の順番を入れ替える
guides(size = guide_legend(order = 2), colour = guide_legend(order = 1)) +
ggtitle("震源分布 1998_2015 (Ｘ印:1943年鳥取地震)")
#dev.off()
#
#
############### 西部(薄い赤色)
#
center <- geocode('yonago')
#maptype = c(“terrain”, “satellite”, “roadmap”, “hybrid”, “toner”, “watercolor”)
yonago <- get_map(c(center$lon, center$lat), zoom = 10, maptype = 'terrain')
#png("eqgg11.png",width=800,height=800)
alpha=0.8 ; size=1
ggmap(yonago) + 
geom_point(data=eq, aes(x=longitude, y=latitude,color=Magnitude),alpha=alpha,size=size)+ 
geom_point(data=eqdata4, aes(x=longitude, y=latitude,size=mag),alpha=1,shape=21,color="gray",fill="black")+ 
#プロットする色を指定したい場合
scale_colour_manual(values=c("red","green","blue"), breaks = c("1 未満","1以上3未満","3 以上"))+
#プロットする大きさの範囲を指定したい場合（マグニチュード 4以上）
scale_size(range = c(6,10))+
# 凡例のタイトルの指定
labs(colour="マグニチュード",size="マグニチュード 4以上")+
#凡例の順番を入れ替える
guides(size = guide_legend(order = 2), colour = guide_legend(order = 1)) +
ggtitle("震源分布 1998_2015")
#dev.off()
#
subeqdata<-subset(eq,longitude>=lon1[1] & longitude<=lon1[2] & latitude>=lat1[1] & latitude<=lat1[2])
subeqdata4<-subset(eqdata4,longitude>=lon1[1] & longitude<=lon1[2] & latitude>=lat1[1] & latitude<=lat1[2])
#
alpha=0.8 ; size=1
p1<-ggplot()+
geom_point(data=subeqdata,aes(x=longitude,y=-depth,color=Magnitude),alpha=alpha,size=size)+
geom_point(data=subeqdata4, aes(x=longitude, y=-depth,size=mag),alpha=1,shape=21,color="gray",fill="black")+ 
#geom_point(aes(x=134.1844,y=0),shape="Ｘ",color="black",size=7)+ 
#プロットする色を指定したい場合
scale_colour_manual(values=c("red","green","blue"), breaks = c("1 未満","1以上3未満","3 以上"))+
#プロットする大きさの範囲を指定したい場合（マグニチュード 4以上）
scale_size(range = c(6,10))+
# 凡例のタイトルの指定
labs(colour="マグニチュード",size="マグニチュード 4以上",x="Longitude",y="Depth")+
#凡例の順番を入れ替える
guides(size = guide_legend(order = 2), colour = guide_legend(order = 1)) +
ggtitle("震源：経度・深さ 1998_2015")
#
p2<-ggplot()+
geom_point(data=subeqdata,aes(x=latitude,y=-depth,color=Magnitude),alpha=alpha,size=size)+
geom_point(data=subeqdata4, aes(x=latitude, y=-depth,size=mag),alpha=1,shape=21,color="gray",fill="black")+ 
#geom_point(aes(x=35.4775,y=0),shape="Ｘ",color="black",size=7)+ 
#プロットする色を指定したい場合
scale_colour_manual(values=c("red","green","blue"), breaks = c("1 未満","1以上3未満","3 以上"))+
#プロットする大きさの範囲を指定したい場合（マグニチュード 4以上）
scale_size(range = c(6,10))+
# 凡例のタイトルの指定
labs(colour="マグニチュード",size="マグニチュード 4以上",x="Latitude",y="Depth")+
#凡例の順番を入れ替える
guides(size = guide_legend(order = 2), colour = guide_legend(order = 1)) +
ggtitle("震源：緯度・深さ 1998_2015")
#
#png("eqgg12.png",width=800,height=800)
grid.newpage() #空の画面を作る
pushViewport(viewport(layout=grid.layout(2, 1))) #画面を区切る
print(p1, vp=viewport(layout.pos.row=1)  ) #1行目
print(p2, vp=viewport(layout.pos.row=2)  ) #2行目
#dev.off()
#
#
############### 中部から鳥取市周辺(薄い青)
#
#lon=134.1,lat=35.4を中心にzoom=10を指定
kurayoshi <- get_map(c(134.1, 35.4), zoom = 10, maptype = 'terrain')
#
#png("eqgg13.png",width=800,height=800)
alpha=0.8 ; size=1
ggmap(kurayoshi) + 
geom_point(data=eq, aes(x=longitude, y=latitude,color=Magnitude),alpha=alpha,size=size)+ 
geom_point(data=eqdata4, aes(x=longitude, y=latitude,size=mag),alpha=1,shape=21,color="gray",fill="black")+ 
geom_point(aes(x=134.1844,y=35.4775),shape="Ｘ",color="black",size=7)+ 
#プロットする色を指定したい場合
scale_colour_manual(values=c("red","green","blue"), breaks = c("1 未満","1以上3未満","3 以上"))+
#プロットする大きさの範囲を指定したい場合（マグニチュード 4以上）
scale_size(range = c(6,10))+
# 凡例のタイトルの指定
labs(colour="マグニチュード",size="マグニチュード 4以上")+
#凡例の順番を入れ替える
guides(size = guide_legend(order = 2), colour = guide_legend(order = 1)) +
ggtitle("震源分布 1998_2015 (Ｘ印:1943年鳥取地震)")
#dev.off()
#
#
subeqdata<-subset(eq,longitude>=lon2[1] & longitude<=lon2[2] & latitude>=lat2[1] & latitude<=lat2[2])
subeqdata4<-subset(eqdata4,longitude>=lon2[1] & longitude<=lon2[2] & latitude>=lat2[1] & latitude<=lat2[2])
#
alpha=0.8 ; size=1
p1<-ggplot()+
geom_point(data=subeqdata,aes(x=longitude,y=-depth,color=Magnitude),alpha=alpha,size=size)+
geom_point(data=subeqdata4, aes(x=longitude, y=-depth,size=mag),alpha=1,shape=21,color="gray",fill="black")+ 
geom_point(aes(x=134.1844,y=0),shape="Ｘ",color="black",size=7)+ 
#プロットする色を指定したい場合
scale_colour_manual(values=c("red","green","blue"), breaks = c("1 未満","1以上3未満","3 以上"))+
#プロットする大きさの範囲を指定したい場合（マグニチュード 4以上）
scale_size(range = c(6,10))+
# 凡例のタイトルの指定
labs(colour="マグニチュード",size="マグニチュード 4以上",x="Longitude",y="Depth")+
#凡例の順番を入れ替える
guides(size = guide_legend(order = 2), colour = guide_legend(order = 1)) +
ggtitle("震源：経度・深さ 1998_2015 (Ｘ印:1943年鳥取地震)")
#
p2<-ggplot()+
geom_point(data=subeqdata,aes(x=latitude,y=-depth,color=Magnitude),alpha=alpha,size=size)+
geom_point(data=subeqdata4, aes(x=latitude, y=-depth,size=mag),alpha=1,shape=21,color="gray",fill="black")+ 
geom_point(aes(x=35.4775,y=0),shape="Ｘ",color="black",size=7)+ 
#プロットする色を指定したい場合
scale_colour_manual(values=c("red","green","blue"), breaks = c("1 未満","1以上3未満","3 以上"))+
#プロットする大きさの範囲を指定したい場合（マグニチュード 4以上）
scale_size(range = c(6,10))+
# 凡例のタイトルの指定
labs(colour="マグニチュード",size="マグニチュード 4以上",x="Latitude",y="Depth")+
#凡例の順番を入れ替える
guides(size = guide_legend(order = 2), colour = guide_legend(order = 1)) +
ggtitle("震源：緯度・深さ 1998_2015 (Ｘ印:1943年鳥取地震)")
#
#png("eqgg14.png",width=800,height=800)
grid.newpage() #空の画面を作る
pushViewport(viewport(layout=grid.layout(2, 1))) #画面を区切る
print(p1, vp=viewport(layout.pos.row=1)  ) #1行目
print(p2, vp=viewport(layout.pos.row=2)  ) #2行目
#dev.off()

2015年には湯梨浜町でマグニチュード４以上の地震がありました。

赤い部分の10月1日から12月31日までの地震の回数

コード

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
library(ggmap)
load("2015.RData")
#
#M>=3、微小地震極微小地震で色を買えるので分類する
#極微小地震
eqdata1<-subset(eqdata,mag<1)
eqdata1$Magnitude<-"1 未満"
#微小地震
eqdata2<-subset(eqdata,mag>=1 & mag<3)
eqdata2$Magnitude<-"1以上3未満"
#M>=3
eqdata3<-subset(eqdata,mag>=3)
eqdata3$Magnitude<-"3 以上"
#
head(eqdata1) ; tail(eqdata1)
head(eqdata2) ; tail(eqdata2)
head(eqdata3) ; tail(eqdata3)
#
#つなげる(ggmapで凡例を付けるため)
eq<-rbind(eqdata1,eqdata2,eqdata3)
#mag>=4(mag>=4の地震はマグニチュードの大きさでプロット)
eqdata4<-subset(eqdata,mag>=4)
#(マグニチュードの昇順に並べ替え)
sortlist <- order(eqdata4$mag)
eqdata4<- eqdata4[sortlist,]
#
lon1<-c(133.85,133.975) ; lat1<-c(35.39,35.49)
yurihama <- get_map(c(133.905,35.437), zoom = 12, maptype = 'terrain')
#
#png("eqgg20.png",width=800,height=800)
alpha=0.8 ; size=1
ggmap(yurihama)+
geom_point(data=eq, aes(x=longitude, y=latitude,color=Magnitude),alpha=alpha,size=size)+ 
geom_point(data=eqdata4, aes(x=longitude, y=latitude,size=mag),alpha=1,shape=21,color="gray",fill="black")+ 
geom_rect(aes(xmin=lon1[1], xmax=lon1[2], ymin=lat1[1],ymax=lat1[2]), fill="red", alpha=0.04)+
#プロットする色を指定したい場合
scale_colour_manual(values=c("red","green","blue"), breaks = c("1 未満","1以上3未満","3 以上"))+
#プロットする大きさの範囲を指定したい場合（マグニチュード 4以上）
scale_size(range = c(6,10))+
# 凡例のタイトルの指定
labs(colour="マグニチュード",size="マグニチュード 4以上")+
#凡例の順番を入れ替える
guides(size = guide_legend(order = 2), colour = guide_legend(order = 1)) +
ggtitle("震源分布 2015 ")
#dev.off()
#
#期間を区切ったデータを扱いやすくするために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"))
#
start<-"2015-10-01"
end<-"2015-12-31"
date<-paste0(start,"::",end)
#
#class(index(eq.xts[date]))
#[1] "POSIXlt" "POSIXt" 
#POSIXltクラスのままだとマージするとき1日ずれる場合がある as.character(index(eq.xts[date]))とする
subeqdata<-data.frame(time=as.character(index(eq.xts[date])),coredata(eq.xts[date]))
#
subeqdata1<-subset(subeqdata,lon1[1]<=longitude & longitude<=lon1[2] & lat1[1]<=latitude & latitude<=lat1[2] & mag<1)
subeqdata2<-subset(subeqdata,lon1[1]<=longitude & longitude<=lon1[2] & lat1[1]<=latitude & latitude<=lat1[2] & mag>=1 & mag<3)
subeqdata3<-subset(subeqdata,lon1[1]<=longitude & longitude<=lon1[2] & lat1[1]<=latitude & latitude<=lat1[2] & mag>=3)
#
#四角で囲まれた部分の極微小地震
#png("hindoyurihama2015.png",width=800,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(start),as.Date(end),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(start),as.Date(end),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(start),as.Date(end),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()

覚え書き

lenovo S100 (メモリを2GBに増設) OS:ZORIN での話です。

data.frame と data.table の処理速度

data.frame
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#load("tottori1998_2014.RData")
system.time(
{
	#M>=3、微小地震、極微小地震で色を買えるので分類する
	#極微小地震
	eqdata1<-subset(eqdata,mag<1)
	eqdata1$Magnitude<-"1 未満"
	#微小地震
	eqdata2<-subset(eqdata,mag>=1 & mag<3)
	eqdata2$Magnitude<-"1以上3未満"
	#M>=3
	eqdata3<-subset(eqdata,mag>=3)
	eqdata3$Magnitude<-"3 以上"
	#つなげる
	eq<-rbind(eqdata1,eqdata2,eqdata3)
}
)

ユーザ システム 経過
21.864 0.784 22.716

data.table
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
library(data.table)
system.time(
{
	eqdata.dt<-setDT(eqdata)
	setkey(eqdata.dt, "mag")
	eqdata1<-eqdata.dt[mag<1,]
	eqdata1$Magnitude<-"1 未満"
	eqdata2<-eqdata.dt[mag>=1 & mag<3,]
	eqdata2$Magnitude<-"1以上3未満"
	eqdata3<-eqdata.dt[mag>=3,]
	eqdata3$Magnitude<-"3 以上"
	l = list(eqdata1,eqdata2,eqdata3)
	# bind correctly by names
	eq<-rbindlist(l, use.names=TRUE)
	eq<-setDF(eq)
}
)
キー設定なし

ユーザ システム 経過
16.288 0.420 16.757

キー設定あり

ユーザ システム 経過
15.836 0.276 16.159

save いろいろ

save(save,saveRDS,saves::saves)
load(load,readRDS,saves::loads)

The Saves Package

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
### save & load
#save(eqdata,file="tottori1998_2014.RData")
#保存前の名前でしか読み込めない
#load("tottori1998_2014.RData")
#
### saveRDS & readRDS
#saveRDS(eqdata,file="tottori1998_2014.RData")
#好きな名前を付けることができる
#eqdata<-readRDS("tottori1998_2014.RData")
#
### saves & loads(saves パッケージ)
#library(saves)
#ファイル名を指定しない場合：ファイル名は「保存前の名前.RDatas」 となる
#saves(eqdata, overwrite = T)
#読み込む変数が指定できる（指定しないと読み込めない）
#eqdata<-loads("eqdata",variables =c("time" ,"longitude" ,"latitude" ,"depth" ,"mag","flag" ,"Magnitude"), to.data.frame=T)
#