R

気象庁地震カタログ(アニメーション）

animation , gganimate , mapdata , xts パッケージ

(参考)
Rで解析：ggplot2のプロットを簡単アニメーション「gganimate」パッケージ

鳥取県震源分布アニメーション（２０００年から２０１５年まで）

animation パッケージ
  • 赤：極微小地震(マグニチュード1未満）
  • 緑：微小地震（マグニチュード1以上３未満）
  • 青：マグニチュード３以上4未満の地震
  • 黒：マグニチュード4以上の地震

鳥取県震源分布アニメーション（マグニチュードの大きさ）

gganimate パッケージ

データ（１９９８年から2015年まで）

西日本震源分布アニメーション（震源の深さ）

gganimate パッケージ

データ（１９９８年から2015年まで）

コード　２０００年から２０１５年まで

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
library(xts)
library(mapdata)
library(animation)
#
load("2015.RData")
eq2015<-eqdata
load("tottori1998_2014.RData")
eqdata<-rbind(eqdata[,-6],eq2015)
#重複した行を取り除く
eqdata<-unique(eqdata)
#
lon1<-c(133,134.6) ; lat1<-c(35,35.7)
#
#期間を区切ったデータを扱いやすくするためにxtsオブジェクトにする
#年月日まで
eq.xts<-as.xts(read.zoo(eqdata))
#年月日時分秒（小数点以下切捨て）まで
#eq.xts<-xts(eqdata[,2:5], strptime(eqdata$time, "%Y-%m-%d %H:%M:%S"))
#
#指定した範囲のデータを抽出
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 & mag<4)
subeq4.xts<-subset(eq.xts,lon1[1]<=longitude & longitude<=lon1[2] & lat1[1]<=latitude & latitude<=lat1[2] & mag>=4)
#
time<-seq(2000,2015)
#
#
saveGIF({
for ( i in 1:16) {
date<-as.character(time[i])
map("japan",xlim=c(lon1[1],lon1[2]),ylim=c(lat1[1],lat1[2]),fill=T,col="gray90")
points(coredata(subeq1.xts[date])[,1],coredata(subeq1.xts[date])[,2],pch=20,col=rgb(1,0,0,alpha=0.8),cex=0.5)
points(coredata(subeq2.xts[date])[,1],coredata(subeq2.xts[date])[,2],pch=20,col=rgb(0,1,0,alpha=0.8),cex=0.8)
if (nrow(subeq3.xts[date])>0) {
points(coredata(subeq3.xts[date])[,1],coredata(subeq3.xts[date])[,2],pch=20,col=rgb(0,0,1,alpha=0.8),cex=coredata(subeq3.xts[date])[,4]/2)
}
#マグニチュード4以上のデータはマグニチュードの昇順にする
eq4<-data.frame(coredata(subeq4.xts[date]))
if (length(eq4$mag>0)){
sortlist <- order(eq4$mag)
eq4 <- eq4[sortlist,] 
points(eq4$longitude ,eq4$latitude,pch=21,col="white",bg=rgb(0,0,0,alpha=0.8),cex=eq4$mag/2)
}
title(paste0("震源分布 ",date),cex.main=3,font.main=1)
}
},movie.name = "eq2000_2015.gif", interval = 1.5,  ani.width = 640,ani.height = 640)

コード　マグニチュードの大きさ

今回は、ggmapではなく、mapdataを使ってみた。
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
library(mapdata)
library(ggplot2)
library(gganimate)
#プロットする範囲
#lon1<-c(133,134.6) ; lat1<-c(35,35.7)
lon1<- c(132.948, 134.7058) ; lat1<-c(34.70972, 36.14203)
#
#あらかじめ作成していた震源データを読み込む
load("2015.RData")
eq2015<-subset(eqdata,lon1[1]<=longitude & longitude<=lon1[2] & lat1[1]<=latitude & latitude<=lat1[2])
#
load("tottori1998_2014.RData")
eqdata<-subset(eqdata,lon1[1]<=longitude & longitude<=lon1[2] & lat1[1]<=latitude & latitude<=lat1[2])
eqdata<-rbind(eqdata[,-6],eq2015)
#重複した行を取り除く
eqdata<-unique(eqdata)
#
#M>=3、微小地震極微小地震で色を買えるので分類する
#極微小地震
eqdata1<-subset(eqdata,mag<1)
eqdata1$Magnitude<-"1 未満"
#微小地震 mag>=1 & mag<2
eqdata2<-subset(eqdata,mag>=1 & mag<2)
eqdata2$Magnitude<-"1以上2未満"
#微小地震 mag>=2 & mag<3
eqdata3<-subset(eqdata,mag>=2 & mag<3)
eqdata3$Magnitude<-"2以上3未満"
#M>=3 & M<4
eqdata4<-subset(eqdata,mag>=3 & mag<4)
eqdata4$Magnitude<-"3以上4未満"
#mag>=4 & mag<5
eqdata5<-subset(eqdata,mag>=4 & mag<5)
eqdata5$Magnitude<-"4以上5未満"
#(マグニチュードの昇順に並べ替え)
sortlist <- order(eqdata5$mag)
eqdata5<- eqdata5[sortlist,]
#mag>=5  & mag<6
eqdata6<-subset(eqdata,mag>=5 & mag<6)
eqdata6$Magnitude<-"5以上6未満"
#(マグニチュードの昇順に並べ替え)
sortlist <- order(eqdata6$mag)
eqdata6<- eqdata6[sortlist,]
#mag>=6
eqdata7<-subset(eqdata,mag>=6)
eqdata7$Magnitude<-"6以上"
#(マグニチュードの昇順に並べ替え)
sortlist <- order(eqdata7$mag)
eqdata7<- eqdata7[sortlist,]
#
head(eqdata1) ; tail(eqdata1)
head(eqdata2) ; tail(eqdata2)
head(eqdata3) ; tail(eqdata3)
head(eqdata4) ; tail(eqdata4)
head(eqdata5) ; tail(eqdata5)
head(eqdata6) ; tail(eqdata6)
head(eqdata7) ; tail(eqdata7)
#
#つなげる
eq<-rbind(eqdata1,eqdata2,eqdata3,eqdata4,eqdata5,eqdata6,eqdata7)
#
tottori <- map("japan",xlim=c(lon1[1],lon1[2]),ylim=c(lat1[1],lat1[2]), plot = FALSE,fill=T)
bd <- borders("japan", "Tottori", colour = rgb(0.8,0.2,0,alpha=1),fill=rgb(0.8,0.2,0,alpha=0.8))
#
p<-ggplot(tottori, aes(long, lat)) +
  geom_polygon(aes(group = group),fill=rgb(0.8,0.2,0,alpha=0.5))+
  coord_cartesian(xlim=c(lon1[1],lon1[2]),ylim=c(lat1[1],lat1[2])) +
  bd +
  geom_point(data=eq, aes(x=longitude, y=latitude,fill=Magnitude,size=mag, frame =Magnitude ),alpha=0.8,shape=21,color="darkgray")+ 
#プロットする色を指定したい場合
  scale_fill_manual(values=c("red","green","yellow","blue","gray50","gray20","black"), 
  breaks = c("1 未満","1以上2未満","2以上3未満","3以上4未満","4以上5未満","5以上6未満","6以上"))+
#プロットする大きさの範囲を指定したい場合
  scale_size(range = c(0.1,10))+
# 凡例のタイトルの指定
  labs(fill="マグニチュード",size="マグニチュード")+
#theme(legend.position="bottom")
   guides(fill=FALSE , size=FALSE)
#print(p)
#
gg_animate(p, "tottori_mag.gif", interval = 1.5)

geom_point(data=eq, aes(x=longitude, y=latitude,fill=Magnitude,size=mag, frame =Magnitude ),alpha=0.8,shape=21,color=”darkgray”)

frame =Magnitude を書き加えただけ

コード　震源の深さ

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
library(mapdata)
library(ggplot2)
library(gganimate)
#library(ggmap)
#プロットする範囲
lon1<- c(129.2554,139) ; lat1<-c(30.66521, 38.75755)
#
#あらかじめ作成していた震源データを読み込む
load("2015.RData")
eq2015<-subset(eqdata,lon1[1]-0.5<=longitude & longitude<=lon1[2]+0.5 & lat1[1]-0.5<=latitude & latitude<=lat1[2]+0.5)
#
load("japan1998_2014.RData")
eqdata<-subset(eqdata,lon1[1]-0.5<=longitude & longitude<=lon1[2]+0.5 & lat1[1]-0.5<=latitude & latitude<=lat1[2]+0.5)
eqdata<-rbind(eqdata[,-6],eq2015)
#重複した行を取り除く
eqdata<-unique(eqdata)
#
#M>=3、微小地震極微小地震で色を買えるので分類する
#極微小地震
eqdata1<-subset(eqdata,mag<1)
eqdata1$Magnitude<-"1 未満"
#微小地震 mag>=1 & mag<2
eqdata2<-subset(eqdata,mag>=1 & mag<2)
eqdata2$Magnitude<-"1以上2未満"
#微小地震 mag>=2 & mag<3
eqdata3<-subset(eqdata,mag>=2 & mag<3)
eqdata3$Magnitude<-"2以上3未満"
#M>=3 & M<4
eqdata4<-subset(eqdata,mag>=3 & mag<4)
eqdata4$Magnitude<-"3以上4未満"
#mag>=4 & mag<5
eqdata5<-subset(eqdata,mag>=4 & mag<5)
eqdata5$Magnitude<-"4以上5未満"
#(マグニチュードの昇順に並べ替え)
sortlist <- order(eqdata5$mag)
eqdata5<- eqdata5[sortlist,]
#mag>=5  & mag<6
eqdata6<-subset(eqdata,mag>=5 & mag<6)
eqdata6$Magnitude<-"5以上6未満"
#(マグニチュードの昇順に並べ替え)
sortlist <- order(eqdata6$mag)
eqdata6<- eqdata6[sortlist,]
#mag>=6
eqdata7<-subset(eqdata,mag>=6)
eqdata7$Magnitude<-"6以上"
#(マグニチュードの昇順に並べ替え)
sortlist <- order(eqdata7$mag)
eqdata7<- eqdata7[sortlist,]
#
head(eqdata1) ; tail(eqdata1)
head(eqdata2) ; tail(eqdata2)
head(eqdata3) ; tail(eqdata3)
head(eqdata4) ; tail(eqdata4)
head(eqdata5) ; tail(eqdata5)
head(eqdata6) ; tail(eqdata6)
head(eqdata7) ; tail(eqdata7)
#
#つなげる
eq<-rbind(eqdata1,eqdata2,eqdata3,eqdata4,eqdata5,eqdata6,eqdata7)
#
#震源の深さ
# depth<10km
eqdata1<-subset(eq,depth<10)
eqdata1$depth_rank<-" 10km 未満"
# 10<=depth<20
eqdata2<-subset(eq,depth>=10 & depth<20)
eqdata2$depth_rank<-" 10km以上20km未満"
# 20<=depth<40
eqdata3<-subset(eq,depth>=20 & depth<40)
eqdata3$depth_rank<-" 20km以上40km未満"
# 40<=depth<80
eqdata4<-subset(eq,depth>=40 & depth<80)
eqdata4$depth_rank<-" 40km以上80km未満"
# 80<=depth<150
eqdata5<-subset(eq,depth>=80 & depth<150)
eqdata5$depth_rank<-" 80km以上150km未満"
# 150<=depth
eqdata6<-subset(eq,depth>=150)
eqdata6$depth_rank<-"150km以上"
#
eq<-rbind(eqdata1,eqdata2,eqdata3,eqdata4,eqdata5,eqdata6)
#
jpn<-map("japan",xlim=c(lon1[1]-0.5,lon1[2]+0.5),ylim=c(lat1[1]-0.5,lat1[2]+0.5), plot = FALSE,fill=T)
#
p<-ggplot(jpn, aes(long, lat)) +
  geom_polygon(aes(group = group),fill=rgb(0.8,0.2,0,alpha=0.5))+
  coord_cartesian(xlim=c(lon1[1],lon1[2]),ylim=c(lat1[1],lat1[2])) +
  geom_point(data=eq, aes(x=longitude, y=latitude,colour=Magnitude,size=mag, frame =depth_rank ),alpha=0.8,shape=20)+ 
#プロットする色を指定したい場合
  scale_colour_manual(values=c("red","green","yellow","blue","gray50","gray20","black")) +
  scale_size(range = c(0,1))+
  guides(colour=FALSE,size=FALSE)
#
gg_animate(p, "japan_depth.gif", interval =2)