R

気象庁地震カタログ(最近9ヶ月 日データ)

raster , mapview パッケージ

(参考)

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

日本の活火山表 List of Active Volcanoes in Japan

まだ地震カタログに入っていない「最近9ヶ月、日データ」を地図上にプロット。
今回は背景となる地図データをraster::getDataで入手。

背景となる地図データを入手

1
2
3
4
5
6
7
library(raster)
j<-getData("alt", country = "JAPAN")
#陰影図を描く準備
slope <- terrain(j, opt = "slope")
aspect <- terrain(j, opt = "aspect")
hill <- hillShade(slope, aspect, 40, 270)
#load("kazan.RData")

例として2016年1月10日のデータを入手

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
date<-"20160110"
#ウェブスクレイピング
url<-paste0("http://www.data.jma.go.jp/svd/eqev/data/daily_map/",date,".html")
doc<-readLines(url,encoding ="UTF-8")
#必要な行を探すための準備
kensaku<-paste0(substr(date,1,4)," ",gsub("\\<0"," ",substr(date,5,6))," ",gsub("\\<0"," ",substr(date,7,8)))
#必要な行を探す
x<-doc[grep(paste0("\\<",kensaku),doc)]
#全角文字 ° が入っているので半角スペース２つに変換する
x<-gsub("°", "  ",x)
#データの整理
time<-paste0(substr(x,1,4),"-",substr(x,6,7),"-",substr(x,9,10)," ",substr(x,12,16),":",substr(x,18,21))
#南緯、西経の場合も考慮に入れると
#地震カタログとは分の有効桁が1桁異なる
latitude = as.numeric(paste0(substr(x,23,23),as.character(as.numeric(substr(x,24,25))+as.numeric(substr(x,28,29))/60+as.numeric(substr(x,31,31))/600)))
longitude = as.numeric(paste0(substr(x,34,34),as.character(as.numeric(substr(x,35,37))+as.numeric(substr(x,40,41))/60+as.numeric(substr(x,43,43))/600)))
#
depth<-as.numeric(substr(x,48,50))
mag<-as.numeric(substr(x,55,58))
eq<-data.frame(time,longitude,latitude,depth,mag)
eqdata<-subset(eq,mag>=-4)
#M>=3、微小地震、極微小地震で色を変えるので分類する
#極微小地震
eqdata1<-subset(eqdata,mag<1)
#微小地震
eqdata2<-subset(eqdata,mag>=1 & mag<3)
#M>=3(マグニチュードの昇順に並べ替え)
eqdata3<-subset(eqdata,mag>=3)
sortlist <- order(eqdata3$mag)
eqdata3<- eqdata3[sortlist,]
#確認
head(eqdata1);tail(eqdata1)
head(eqdata2);tail(eqdata2)
head(eqdata3);tail(eqdata3)

地図上にプロット

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#png("eqjapan01.png",width=1000,height=800)
alpha=0.8
plot(hill, col = grey(0:100/100), legend = F,axes=F, box=F)
plot(j, col = terrain.colors(20, alpha = 0.3), add = TRUE)
#points(kazan$longitude,kazan$latitude,pch=17,col="orange",cex=0.5)
points(eqdata1$longitude,eqdata1$latitude,cex=(eqdata1$mag+4)/5,col=rgb(1,0,0,alpha=alpha),pch=20)
points(eqdata2$longitude,eqdata2$latitude,cex=(eqdata2$mag+4)/5,col=rgb(0,1,0,alpha=alpha),pch=20)
points(eqdata3$longitude,eqdata3$latitude,cex=(eqdata3$mag+4)/5,col=rgb(0,0,1,alpha=alpha),pch=20)
#legend("topleft",title="Magnitude",legend = c("M<1",expression(paste(1<=M,"<3")),expression(M>=3),"活火山"),  
#cex=1,pch=c(20,20,20,17),col=c(rgb(1,0,0,alpha=1),rgb(0,1,0,alpha=1),rgb(0,0,1,alpha=1),"orange"),
#pt.cex =2,bg ="gray", x.intersp = 1, y.intersp =1)
legend("topleft",title="Magnitude",legend = c("M<1",expression(paste(1<=M,"<3")),expression(M>=3)),  
cex=1,pch=c(20,20,20),col=c(rgb(1,0,0,alpha=1),rgb(0,1,0,alpha=1),rgb(0,0,1,alpha=1)),
pt.cex =2,bg ="gray", x.intersp = 1, y.intersp =1)
title(paste(date,"地震  (微小地震、極微小地震を含む)"))
#dev.off()

描く範囲を狭くする場合

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
lon<-c(130.8,136)
lat<-c(33.5,36)
e <- extent(lon[1], lon[2], lat[1], lat[2])
rc <- crop(j, e)
slope1 <- terrain(rc, opt = "slope")
aspect1 <- terrain(rc, opt = "aspect")
hill1 <- hillShade(slope1, aspect1, 40, 270)
#
subeqdata<-subset(eqdata,longitude>=lon[1] & longitude<=lon[2] & latitude>=lat[1] & latitude<=lat[2] & mag>=-4)
#
#M>=3、微小地震、極微小地震で色を買えるので分類する
#極微小地震
subeqdata1<-subset(subeqdata,mag<1)
#微小地震
subeqdata2<-subset(subeqdata,mag>=1 & mag<3)
#M>=3(マグニチュードの昇順に並べ替え)
subeqdata3<-subset(subeqdata,mag>=3)
sortlist <- order(subeqdata3$mag)
subeqdata3<- subeqdata3[sortlist,]
#
head(subeqdata1);tail(subeqdata1)
head(subeqdata2);tail(subeqdata2)
head(subeqdata3);tail(subeqdata3)
#
#png("eqjapan02.png",width=1000,height=800)
alpha=0.6
plot(hill1, col = grey(0:100/100), legend = FALSE, main = "",axes=F, box=F)
plot(rc, col = terrain.colors(20, alpha = 0.3), add = TRUE)
#points(kazan$longitude,kazan$latitude,pch=17,col="orange",cex=2)
points(subeqdata1$longitude,subeqdata1$latitude,cex=(subeqdata1$mag+4)/2,col=rgb(1,0,0,alpha=alpha),pch=20)
points(subeqdata2$longitude,subeqdata2$latitude,cex=(subeqdata2$mag+4)/2,col=rgb(0,1,0,alpha=alpha),pch=20)
points(subeqdata3$longitude,subeqdata3$latitude,cex=(subeqdata3$mag+4)/2,col=rgb(0,0,1,alpha=alpha),pch=20)
#legend("topleft",title="Magnitude",legend = c("M<1",expression(paste(1<=M,"<3")),expression(M>=3),"活火山"),  
#cex=1,pch=c(20,20,20,17),col=c(rgb(1,0,0,alpha=1),rgb(0,1,0,alpha=1),rgb(0,0,1,alpha=1),"orange"),
#pt.cex =2,bg ="gray", x.intersp = 1, y.intersp =1)
legend("topleft",title="Magnitude",legend = c("M<1",expression(paste(1<=M,"<3")),expression(M>=3)),  
cex=1,pch=c(20,20,20),col=c(rgb(1,0,0,alpha=1),rgb(0,1,0,alpha=1),rgb(0,0,1,alpha=1)),
pt.cex =2,bg ="gray", x.intersp = 1, y.intersp =1)
title(paste(date,"地震  (微小地震、極微小地震を含む)"))
#dev.off()

mapviewパッケージを使ってみる

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
map1<-eqdata1 
map2<-eqdata2
map3<-eqdata3
#
library(leaflet)
library(mapview)
library(sp)
library(raster)
#
coordinates(map1) <- ~longitude+latitude
proj4string(map1) <- CRS("+init=epsg:4326")
map1<-spTransform(map1, CRS("+init=epsg:3857"))
#
coordinates(map2) <- ~longitude+latitude
proj4string(map2) <- CRS("+init=epsg:4326")
map2<-spTransform(map2, CRS("+init=epsg:3857"))
#
coordinates(map3) <- ~longitude+latitude
proj4string(map3) <- CRS("+init=epsg:4326")
map3<-spTransform(map3, CRS("+init=epsg:3857"))
#
e1<-mapView(map1, color = "red",cex=5) 
e2<-mapView(map2, color = "green",cex=5) 
e3<-mapView(map3, color = "blue",cex=5) 
#
e1 + e2 + e3
#
#ちなみに
#nrow(eqdata)
#[1] 672

結果は省略(背景となる地図はOpenStreetMap , OpenTopoMap などから選択できます。)