Nov 6 2016

マグニチュード６以上

ggplot2 , mapdata , xts パッケージ

東海・東南海・南海地震の想定震源域のデータをFuyuki Hirose’s HP で見つけたのでRで地図上に描いてみた。（プラス：マグニチュード６以上の震源）

(参考 : データの引用文献等)
関東地震，東海地域のアスペリティ，東海・東南海・南海地震の想定震源域
駿河トラフ,相模トラフ,南海トラフ

Fuyuki Hirose’s HP
Fuyuki Hirose’s HP:フィリピン海プレート

日本周辺で起きたマグニチュード６以上の地震の回数

東海・東南海・南海地震の想定震源域とマグニチュード６以上の震源

（準備）

関東地震，東海地域のアスペリティ，東海・東南海・南海地震の想定震源域
駿河トラフ,相模トラフ,南海トラフにヘッダーを加えて、csvファイルに変換した。
気象庁の地震カタログ、震源リストからデータをダウンロード、加工、日本周辺で起きたマグニチュード６以上データを抽出。

今回はコードの一部のみ

地図を描いた部分のコード

tokai=read.csv("./mapdata/tokai.csv")
tonankai=read.csv("./mapdata/tonankai.csv")
nankai=read.csv("./mapdata/nankai.csv")
# Asperities
# Kanto (Wald and Somerville, 1995, BSSA)
kasp=read.csv("./mapdata/kanto_eq.csv")
# Tokai (Matsumura, 1997, Tectono.)
tasp=read.csv("./mapdata/tokai_asperity.csv")
#
trench=read.csv("./mapdata/trench.csv")
#voldata=read.csv("./mapdata/voldata.csv")
#地図上に書く文字と書く位置
trough<-data.frame(names=c("駿河トラフ","相模トラフ","南海トラフ","フィリピン海プレート","東海","東南海","南海"),
   longitude=c(138.9,140.2,136.5,136,138,136.8,134),
   latitude=c(33.8,34.4,32.5,30.7,34.5,34,32.9))
#
library("ggplot2")
library("maps")
library("mapdata")
library("geosphere")
# plot limits
xlim = c(130,141)
ylim = c(30,37)
#
japan = map_data("japan")
#
p<-ggplot() +
    borders("japan",fill="lightgray", size = 0.3) +
    theme_bw() +
    coord_map(xlim=xlim,ylim=ylim) +
    labs(y="",x="")
#想定震源域等(geom_polygon) 
p <-p +geom_polygon(data = tokai,aes(x = longitude, y = latitude),fill="green",alpha = 0.3)
p <-p +geom_polygon(data = tonankai,aes(x = longitude, y = latitude),fill="red",alpha = 0.3)
p <-p +geom_polygon(data = nankai,aes(x = longitude, y = latitude),fill="blue",alpha = 0.3)
p <-p +geom_polygon(data = kasp,aes(x = longitude, y = latitude),fill="yellow",alpha = 0.3)
p <-p +geom_polygon(data = tasp,aes(x = longitude, y = latitude),fill="darkgreen",alpha = 0.3)
#trench(geom_line)
p <-p +geom_line(data = trench,aes(x = longitude, y = latitude),colour="gray30",linetype="dashed")
#p <-p +geom_point(data = voldata,aes(x = longitude, y = latitude),colour="red",shape=17,size=2)
#文字(geom_text)
p<-p + geom_text(data = trough,aes(x = longitude, y = latitude,label =names))
#矢印の終点を計算
arrow1<-destPointRhumb(c(140.8,33.5),90-131, 63750/3*2)
arrow2<-destPointRhumb(c(137,31.1),90-145, 105000/3*2)
#矢印(geom_path)
p<-p+ geom_path(aes(x=c(140.8,arrow1[1]),y=c(33.5,arrow1[2])),arrow=arrow(angle = 20,type = "closed"),size=1) +
geom_path(aes(x=c(137,arrow2[1]),y=c(31.1,arrow2[2])),arrow=arrow(angle = 20,type = "closed"),size=1)
#
p<-p + geom_point(data=eq, aes(x=longitude, y=latitude,fill=depth_rank,size=mag),alpha=0.8,shape=21,color="gray20")+ 
#プロットする色を指定したい場合
  scale_fill_manual(values=c("red","orange","gold4","green","blue","purple"), 
  breaks = c(" 10km 未満"," 10km以上20km未満"," 20km以上40km未満"," 40km以上80km未満"," 80km以上150km未満","150km以上"))+
# 凡例のタイトルの指定
  labs(size="Magnitude",fill="Depth",title="震源分布（マグニチュード６以上）1923 ～ 2016/10/31[データ：気象庁]") +
#凡例の順番を入れ替える
  guides(size = guide_legend(order = 1), fill = guide_legend(order = 2)) +
  theme(text=element_text(size=12, family="TakaoExMincho"))
#png("eq61923_20161031_02.png",width=800,height=600)
p
#dev.off()

時間のデータは秒まででそれ以下は切り捨てた。

地震カタログの1967_1982のデータに緯度経度が分までのデータがあったため過去の記事のコードでは緯度経度がNAになってしまった。

その部分のコード

時間

1	time<-paste0(substr(x,2,5),"-",substr(x,6,7),"-",substr(x,8,9)," ",substr(x,10,11),":",substr(x,12,13),":",substr(x,14,15))

震央の緯度、震央の経度（ ifelse(is.na( ）

1
2

latitude = as.numeric(paste0(substr(x,22,22),as.character(as.numeric(substr(x,23,24))+as.numeric(substr(x,25,26))/60+ifelse(is.na(as.numeric(substr(x,27,28))), 0,as.numeric(substr(x,27,28)) )/6000)))
longitude = as.numeric(paste0(substr(x,33,33),as.character(as.numeric(substr(x,34,36))+as.numeric(substr(x,37,38))/60+ifelse(is.na(as.numeric(substr(x,39,40))), 0,as.numeric(substr(x,39,40)) )/6000)))

続けてダウンロード

続けてダウンロード1

eq<-read.table(text="",col.names=c("time", "longitude", "latitude", "depth", "mag", "flag"))
date<-c("h1923","h1951","h1961","h1967",paste0("h",as.character(seq(1983,1996))),"h199701","h199710",paste0("h",as.character(seq(1998,2015))) )
#
for (i in 1:length(date)){
temp <- tempfile()
download.file(paste0("http://www.data.jma.go.jp/svd/eqev/data/bulletin/data/hypo/",date[i],".zip"),temp)
x<- readLines(unz(temp,date[i]))
x<-x[grep("^J", x)]

続けてダウンロード2(rbindでつなげる)

1
2
3

#データフレーム作成
eq<-rbind(eq,data.frame(time,longitude,latitude,depth,mag,flag))
}