R

気象庁地震カタログ(密度分布）

ggmap , spatstat パッケージ

(参考)

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

(注意)

日本の地震カタログ（1923年から現在）の解説
「検知能力に関する主な変更」

鳥取県の震源分布

密度分布を調べる範囲は下の地図の薄い赤で塗りつぶした部分

密度分布

（注意）

  • 年ごとの比較をするために密度に対する色を統一(しかし、観測地点や検知能力は一定ではないのでその点も考慮してください。)
  • 色分けは2015年の色分けを参考に設定した。(2015年に比較して大きい密度の箇所は着色されていない)
  • 2015年のデータだけは日々のデータから取り出したので震央の位置精度が一桁落ちる。

2000年

2001年

2002年

2003年

2004年

2005年

2006年

2007年

2008年

2009年

2010年

2011年

2012年

2013年

2014年

2015年

地震の発生数(極微小地震、微小地震を含む)

2000年

2001年

2002年

2003年

2004年

2005年

2006年

2007年

2008年

2009年

2010年

2011年

2012年

2013年

2014年

2015年

2000_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
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
load("2015.RData")
eq2015<-subset(eqdata,longitude>=132.5 & longitude<=135 & latitude>=34.5 & latitude<=35.9 & mag>=-4)
load("tottori1998_2014.RData")
#
eqdata<-rbind(eqdata[,-6],eq2015)
#
library(ggmap)
library(mapproj)
#
#
#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(大きい地震はマグニチュードの大きさでプロット)
eqdata4<-subset(eqdata,mag>=4 )
sortlist <- order(eqdata4$mag)
eqdata4<- eqdata4[sortlist,]
#
##### 鳥取県震源分布図 
#
lon1<-c(133,134.65)
lat1<-c(35.1,35.65)
#
#倉吉を中心に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')
#
#png("eqgg000.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=10)+ 
#geom_rect(aes(xmin=lon1[1], xmax=lon1[2], ymin=lat1[1],ymax=lat1[2]),fill="red", 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()
#
###### 密度分布等
#
library(spatstat)
#期間を区切ったデータを扱いやすくするために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"))
#
#例として2000年
date<-"2000"
#
subeqdata<-data.frame(time=index(eq.xts[date]),coredata(eq.xts[date]))
#
subeq<-subset(subeqdata,lon1[1]<=longitude & longitude<=lon1[2] & lat1[1]<=latitude & latitude<=lat1[2])
px<-subeq$longitude
py<-subeq$latitude
pnt <- ppp(px, py, c(lon1[1],lon1[2]), c(lat1[1],lat1[2]))
#データ点の個数をカウント
Q <- quadratcount(pnt, nx = 33, ny = 11)
#
#密度分布
#
#png(paste0("density",date,"_1.png"),width=1000,height=400)
plot(density(pnt),col=topo.colors(100),axes=T,main=paste0("density(",date,")"))
contour(density(pnt),lwd=0.5,add=T)
plot(pnt,pch=20,cex=0.3,col="gray50",add=T)
#dev.off()
#
#年ごとの比較をするために密度に対する色を統一
#色分けは2015年の色分けに近い。
#2015年に比較して大きい密度の箇所は着色されていない
#
tc <- colourmap(topo.colors(100), breaks=seq(0,30000,length=101))
#png(paste0("density",date,"_2.png"),width=1000,height=400)
plot(density(pnt),col=tc,axes=T,main=paste0("density(",date,") 色付け範囲固定"))
contour(density(pnt),lwd=0.5,add=T)
plot(pnt,pch=20,cex=0.3,col="gray50",add=T)
#dev.off()
#
#png(paste0("density",date,"_3.png"),width=1000,height=400)
contour(density(pnt),axes=T,col=rgb(0,0,1,alpha=0.5),xlab="lon",ylab="lat",main=paste0("地震の発生数(",date,")"))
#plot(pnt,cols=rgb(0,0,0,alpha=0.3),cex=0.2,add=T)
plot(pnt,cols=rgb(0,0,0,alpha=0.1),cex=0.2,add=T)
#plot(Q, add = TRUE, cex =0.8,col="red")
plot(Q, add = TRUE, cex =0.5,col="red")
#dev.off()