sqlite3

データベース03

RSQLiteパッケージ

日本の地震カタログ（1923年から現在）の解説より抜粋
検知能力に関する主な変更
1997/10/01〜 大学や関係機関の観測点を震源計算に使用開始
2000/10/01〜 大阪管内、福岡管内のHi-net観測点を震源計算に使用開始
2001/10/01〜 札幌管内、仙台管内のHi-net観測点を震源計算に使用開始
2002/10/01〜 東京管内のHi-net観測点を震源計算に使用開始
2003/03/03〜 トリガーパラメータ変更により震源決定総数が減少（検知能力に変更なし）
2003/05/01〜 振幅値の分解能の向上

（過去の記事）
気象庁地震カタログ01

気象庁地震カタログから作成したデータフレームをデータベースに INSERT

この部分はメモリを多く積んだパソコンで行った。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
load("japan1998_2014.RData")
head(eqdata)
#
library(RSQLite)
driver=dbDriver("SQLite")
dbname="japan1998_2014.db"
con=dbConnect(driver,dbname)
# データフレームをテーブルとして一括書き込み
rs=dbWriteTable(con,"japan1998_2014",eqdata,row.names=F)
# 書き込めているかテスト
dbGetQuery(con,"SELECT * FROM japan1998_2014 limit 5;")
# テーブルをデータフレームとして一括取得
#dbReadTable(con,"japan1998_2014")
#データベース内のテーブルを確認（.tables;(show tables;)に相当）
dbListTables(con)
#データベース内のデータセットをRの変数に代入
eqdata <- dbGetQuery(con, "select * from japan1998_2014")
#データベース接続を終了
dbDisconnect(con)

データフレームの取得(必要な期間のデータ)

上で作成したdb ファイルをネットブックにコピーして以下のコードを実行。

1
2
3
4
5
6
7
8
library(mapdata)
library(xts)
library(RSQLite)
#
driver=dbDriver("SQLite")
dbname="japan1998_2014.db"
con=dbConnect(driver,dbname)
#dbGetQuery(con,"SELECT * FROM japan1998_2014 limit 5;")

2003/05/01以降に起きた大きな地震の震源域(約1週間）

新潟県中越地震 2004年10月23日

  • 北緯37度17分30秒
  • 東経138度52分0秒 
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
# 2004年10月のデータを抽出
eqdata <- dbGetQuery(con, "SELECT * FROM japan1998_2014 WHERE time LIKE '2004-10-%';")
# 2004-10-23 から　2004-10-29
#eqdata <- dbGetQuery(con, "SELECT * FROM japan1998_2014 WHERE time >= '2004-10-23' AND  time <= '2004-10-30' ;")
#137.87 <= longitude <= 139.87),36.3 <= latitude <= 38.3
#eqdata <- dbGetQuery(con, "SELECT * FROM japan1998_2014 WHERE longitude >= 137.87 AND longitude <= 139.87 AND latitude >= 36.3 AND latitude  <= 38.3 ;")
#
subset(eqdata,mag==max(eqdata$mag))
#                      time longitude latitude depth mag flag
#166 2004-10-23 17:56:00.30  138.8672  37.2925 13.08 6.8    K
eqdata.xts<-as.xts(read.zoo(eqdata[,-6]))
#
period<-"2004-10-23::2004-10-29"
eqdata<-data.frame(eqdata.xts[period])
#
#震源の深さ
# depth<10km
eqdata1<-subset(eqdata,depth<10)
eqdata1$depth_rank<-" 10km 未満"
eqdata1$col<-"red"
# 10<=depth<20
eqdata2<-subset(eqdata,depth>=10 & depth<20)
eqdata2$depth_rank<-" 10km以上20km未満"
eqdata2$col<-"orange"
# 20<=depth<40
eqdata3<-subset(eqdata,depth>=20 & depth<40)
eqdata3$depth_rank<-" 20km以上40km未満"
eqdata3$col<-"gold4"
# 40<=depth<80
eqdata4<-subset(eqdata,depth>=40 & depth<80)
eqdata4$depth_rank<-" 40km以上80km未満"
eqdata4$col<-"green"
# 80<=depth<150
eqdata5<-subset(eqdata,depth>=80 & depth<150)
eqdata5$depth_rank<-" 80km以上150km未満"
eqdata5$col<-"blue"
# 150<=depth
eqdata6<-subset(eqdata,depth>=150)
eqdata6$depth_rank<-"150km以上"
eqdata6$col<-"purple"
#
eq<-rbind(eqdata1,eqdata2,eqdata3,eqdata4,eqdata5,eqdata6)
#並べ替え:マグニチュードの昇順
sortlist <- order(eq$mag)
eq<- eq[sortlist,]
##
#mag>=1の地震を抽出
subeq<-subset(eq,mag>=1)
#mag<1
mag1und<-subset(eq,mag<1)
#
par(family = "TakaoExMincho")
#png("eq20041023_02.png",width=1000,height=800)
map('japan',xlim=c(137.87,139.87),ylim=c(36.3,38.3),fill=T,col="gray85",mar = c(5,4,4,15))
map('japan',col="gray10",boundary =T,add=T)
points(mag1und$longitude,mag1und$latitude,cex=0.2,bg=adjustcolor(subeq$col,0.7),pch=21,col="black")
points(subeq$longitude,subeq$latitude,cex=(subeq$mag-0.2)/2,bg=adjustcolor(subeq$col,0.7),pch=21,col="black")
par(xpd=T)
#
cexsize=(c(1,3,5,7)-0.2)/2
legend(par("usr")[2]*1.001,par("usr")[3]*0.05+par("usr")[4]*0.95,title="Magnitude",
legend = c("M<1","M=1","M=3","M=5","M=7"),  
col=1,pch=21,pt.cex=c(0.2,cexsize),
cex =1,bg ="gray90", x.intersp = 2, y.intersp =1.8,bty = "n")
#
legend(par("usr")[2],par("usr")[3]*0.6+par("usr")[4]*0.4,title="Depth",
legend = c("D<10km","10km<=D<20km","20km<=D<40km","40km<=D<80km","80km<=D<150km","150km<=D"),  
cex=1,pch=19,col=c("red","orange","gold4","green","blue","purple"),
pt.cex =1.5,bg ="gray90", x.intersp = 1, y.intersp =1.2,bty = "n")
par(xpd=F)
title(paste0(period," 新潟県中越地震"))
title("","( 地震データ：気象庁地震カタログ )",line=2.5)
map.axes(cex.axis=0.8,las=1)
#dev.off()

岩手・宮城内陸地震 2008年6月14日

  • 北緯39度01.7分
  • 東経140度52.8分 
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
# 2008年6月のデータを抽出
eqdata <- dbGetQuery(con, "SELECT * FROM japan1998_2014 WHERE time LIKE '2008-06-%';")
subset(eqdata,mag==max(eqdata$mag))
#                       time longitude latitude depth mag flag
#4681 2008-06-14 08:43:45.36  140.8807 39.02983  7.77 7.2    K
eqdata.xts<-as.xts(read.zoo(eqdata[,-6]))
#
period<-"2008-06-14::2008-06-20"
eqdata<-data.frame(eqdata.xts[period])
#
#震源の深さ
# depth<10km
eqdata1<-subset(eqdata,depth<10)
eqdata1$depth_rank<-" 10km 未満"
eqdata1$col<-"red"
# 10<=depth<20
eqdata2<-subset(eqdata,depth>=10 & depth<20)
eqdata2$depth_rank<-" 10km以上20km未満"
eqdata2$col<-"orange"
# 20<=depth<40
eqdata3<-subset(eqdata,depth>=20 & depth<40)
eqdata3$depth_rank<-" 20km以上40km未満"
eqdata3$col<-"gold4"
# 40<=depth<80
eqdata4<-subset(eqdata,depth>=40 & depth<80)
eqdata4$depth_rank<-" 40km以上80km未満"
eqdata4$col<-"green"
# 80<=depth<150
eqdata5<-subset(eqdata,depth>=80 & depth<150)
eqdata5$depth_rank<-" 80km以上150km未満"
eqdata5$col<-"blue"
# 150<=depth
eqdata6<-subset(eqdata,depth>=150)
eqdata6$depth_rank<-"150km以上"
eqdata6$col<-"purple"
#
eq<-rbind(eqdata1,eqdata2,eqdata3,eqdata4,eqdata5,eqdata6)
#並べ替え:マグニチュードの昇順
sortlist <- order(eq$mag)
eq<- eq[sortlist,]
##
#mag>=1の地震を抽出
subeq<-subset(eq,mag>=1)
#mag<1
mag1und<-subset(eq,mag<1)
#
par(family = "TakaoExMincho")
#png("eq20080614_02.png",width=1000,height=800)
map('japan',xlim=c(139.8807,141.8807),ylim=c(38.02983,40.02983),fill=T,col="gray85",mar = c(5,4,4,15))
map('japan',col="gray10",boundary =T,add=T)
points(mag1und$longitude,mag1und$latitude,cex=0.2,bg=adjustcolor(subeq$col,0.7),pch=21,col="black")
points(subeq$longitude,subeq$latitude,cex=(subeq$mag-0.2)/2,bg=adjustcolor(subeq$col,0.7),pch=21,col="black")
par(xpd=T)
#
cexsize=(c(1,3,5,7)-0.2)/2
legend(par("usr")[2]*1.001,par("usr")[3]*0.05+par("usr")[4]*0.95,title="Magnitude",
legend = c("M<1","M=1","M=3","M=5","M=7"),  
col=1,pch=21,pt.cex=c(0.2,cexsize),
cex =1,bg ="gray90", x.intersp = 2, y.intersp =1.8,bty = "n")
#
legend(par("usr")[2],par("usr")[3]*0.6+par("usr")[4]*0.4,title="Depth",
legend = c("D<10km","10km<=D<20km","20km<=D<40km","40km<=D<80km","80km<=D<150km","150km<=D"),  
cex=1,pch=19,col=c("red","orange","gold4","green","blue","purple"),
pt.cex =1.5,bg ="gray90", x.intersp = 1, y.intersp =1.2,bty = "n")
par(xpd=F)
title(paste0(period," 岩手・宮城内陸地震"))
title("","( 地震データ：気象庁地震カタログ )",line=2.5)
map.axes(cex.axis=0.8,las=1)
#dev.off()

1
2
#データベース接続を終了
dbDisconnect(con)