R

気象庁地震カタログ(アニメーション:鳥取県西部地震）

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

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

鳥取県西部地震 2000年（平成12年）10月6日13時30分18秒

2000年～2015年（年）

ggmap

2000年10月から2002年12月まで（月）

mapdata

2000年10月1日から10月31日まで（日）

mapdata

2000年10月6日13時から10月11日23時まで（時間）

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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
library(mapdata)
library(ggplot2)
#devtools::install_github("dgrtwo/gganimate")
library(gganimate)
library(xts)
library(ggmap)
#プロットする範囲
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)
#
##################
#
bd <- borders("japan", "Tottori", colour =rgb(0.8,1,0.4,alpha=0.7),fill=rgb(0.8,1,0.4,alpha=0.5)) 
#年ごとにプロットする準備
eq$period<-substring(eq$time,1,4)
#月ごとにプロットする準備
#eq$period<-substring(eq$time,1,7)
#日ごとにプロットする準備
#eq$period<-substring(eq$time,1,10)
#プロットする期間を抽出する場合
eq.xts<-xts(read.zoo(eq))
subeq<-data.frame(coredata(eq.xts["2000::2015"]),stringsAsFactors = F)
#subeq<-data.frame(coredata(eq.xts["2000-10::2002-12"]),stringsAsFactors = F)
#subeq<-data.frame(coredata(eq.xts["2000-10"]),stringsAsFactors = F)
subeq<-transform(subeq,longitude = as.numeric(longitude),latitude= as.numeric(latitude),depth= as.numeric(depth),mag= as.numeric(mag))
#マグニチュードの昇順に並べ替え
sortlist <- order(subeq$mag)
subeq<- subeq[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')
#
p<-ggmap(tottori) + 
  coord_cartesian(xlim=c(lon1[1],lon1[2]),ylim=c(lat1[1],lat1[2])) +
  bd +
  geom_point(data=subeq, aes(x=longitude, y=latitude,fill=Magnitude,size=mag, frame =period),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))+
  guides(fill=FALSE , size=FALSE)
#print(p)
#
gg_animate(p, "tottori_year.gif", interval = 2)
#gg_animate(p, "tottori_seibumon.gif", interval = 1.5)
#gg_animate(p, "tottori_seibuday.gif", interval = 1.5)

1時間ごと(コードの後半）

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
#時間ごとにプロットする準備
eq$period<-substring(eq$time,1,13)
#プロットする期間を抽出する場合
#時間単位で取り出す場合 xts(eq[,-1], as.POSIXct(eq$time)) とする
eq.xts<-xts(eq[,-1], as.POSIXct(eq$time))
#
subeq<-data.frame(coredata(eq.xts["2000-10-06 13::2000-10-11 23"]),stringsAsFactors = F)
subeq<-transform(subeq,longitude = as.numeric(longitude),latitude= as.numeric(latitude),depth= as.numeric(depth),mag= as.numeric(mag))
#マグニチュードの昇順に並べ替え
sortlist <- order(subeq$mag)
subeq<- subeq[sortlist,]
#
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=subeq, aes(x=longitude, y=latitude,fill=Magnitude,size=mag, frame =period),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_seibuhour.gif", interval = 1.5)