地震

震源地と原子力発電所03

ggplot2 , maps , mapdata , geosphere , gridExtra パッケージ

（前の記事）
震源地と原子力発電所
震源地と原子力発電所02

今回は日本国内の各原子力発電所近くで起きたマグニチュード５以上の地震をプロット。
震源データは気象庁。マグニチュード5以上のデータ。　　

原子力発電所の位置は震源地と原子力発電所 と同じ

左：1923年から日本で最初の原子力発電が行われた1963年10月26日以前の震源データで作成した図
右：1923年から2016-11-21までの震源データで作成した図　　

泊発電所

東通原子力発電所

女川原子力発電所

福島第二原子力発電所

東海第二発電所

柏崎刈羽原子力発電所

浜岡原子力発電所

志賀原子力発電所

敦賀発電所

美浜発電所

大飯発電所

高浜発電所

島根原子力発電所

伊方発電所

玄海原子力発電所

川内原子力発電所

(準備)過去記事参照
気象庁地震カタログ、震源データからマグニチュード5以上のデータを抽出。保存。eq51923_20161121.RData
原発の位置を保存。in_service.dat

データをロード。震源のプロット準備。

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
library("ggplot2")
library("maps")
library("mapdata")
library("geosphere")
library(gridExtra)
load("in_service.dat")
nujap<-in_service[in_service$Country=="Japan",]
# nrow(nujap)
# plot limits
#
load("eq51923_20161121.RData")
eqdata<-eq51923_20161121
#
#震源の深さ
# 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,]
#
# 1963-10-26以前のデータを抽出
eq1963<-eq[substring(eq$time,1,10) < "1963-10-27",]

地図作成。保存。

今回のポイントとなる関数は、

  • destPointRhumb(地図にする範囲を決めるのに使用)
  • destPoint(地図上に円を描くのに使用)
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
japan = map_data("japan")
#
#60km
d=60000
#
for (i in 1:nrow(nujap)) {
lon1<-destPointRhumb(c(nujap$longitude[i],nujap$latitude[i]), 270, d, r = 6378137)[1]
lon2<-destPointRhumb(c(nujap$longitude[i],nujap$latitude[i]), 90, d, r = 6378137)[1]
Lon.range = c(lon1,lon2)
#
lon1<-destPointRhumb(c(nujap$longitude[i],nujap$latitude[i]), 180, d, r = 6378137)[2]
lon2<-destPointRhumb(c(nujap$longitude[i],nujap$latitude[i]), 0, d, r = 6378137)[2]
Lat.range = c(lon1,lon2)
#Lon.range = c(nujap$longitude[i]-0.5, nujap$longitude[i]+0.5)
#Lat.range = c(nujap$latitude[i]-0.5, nujap$latitude[i]+0.5)
#
p0<-ggplot() +
    borders("japan",fill="gray75", size = 0.3) +
    theme_bw() +
    coord_map(xlim=Lon.range ,ylim=Lat.range) +
    labs(y="",x="")
#マグニチュード5以上
eq1=eq1963[ eq1963$mag>=5,]
p0<-p0 + geom_point(data=eq1, 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以上"))+
#プロットする大きさの範囲を指定したい場合
  scale_size(breaks=c(5,6,7,8,9),range = c(3,10)) +
# 凡例のタイトルの指定
  labs(size="Magnitude",fill="Depth",title="震源分布（マグニチュード5以上）1923 ～ 1963-10-26[データ：気象庁]") +
#凡例の順番を入れ替える
  guides(size = guide_legend(order = 1), fill = guide_legend(order = 2)) +
  theme(text=element_text(size=12, family="TakaoExMincho"))
p0<- p0 + geom_point(data=nujap,aes(x=longitude, y=latitude),shape=24,size =5, fill=rgb(0,0,0,0.8),color="black")
p0<- p0 + geom_text(data=nujap[i,],aes(x=longitude, y=latitude, label = station),vjust = 0, nudge_y = 0.02)
#30km円
circle30=as.data.frame(destPoint(c(nujap$longitude[i],nujap$latitude[i]), b=1:365, d=30000))
p0<-p0 + geom_polygon(data=circle30, aes(x=lon,y=lat), color="black",fill="red",alpha=0.2)
p0<-p0 + geom_text(aes(x=nujap$longitude[i],y=max(circle30$lat)),label="30km", nudge_y = 0.02)
#50km円
circle50=as.data.frame(destPoint(c(nujap$longitude[i],nujap$latitude[i]), b=1:365, d=50000))
p0<-p0 + geom_polygon(data=circle50, aes(x=lon,y=lat), color="black",fill="blue",alpha=0.2)
p0<-p0 + geom_text(aes(x=nujap$longitude[i],y=max(circle50$lat)),label="50km", nudge_y = 0.02)
#
p<-ggplot() +
    borders("japan",fill="gray75", size = 0.3) +
    theme_bw() +
    coord_map(xlim=Lon.range ,ylim=Lat.range) +
    labs(y="",x="")
#マグニチュード5以上
eq1=eq[ eq$mag>=5,]
p<-p + geom_point(data=eq1, 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以上"))+
#プロットする大きさの範囲を指定したい場合
  scale_size(breaks=c(5,6,7,8,9),range = c(3,10)) +
# 凡例のタイトルの指定
  labs(size="Magnitude",fill="Depth",title="震源分布（マグニチュード5以上）1923 ～ 2016-11-21[データ：気象庁]") +
#凡例の順番を入れ替える
  guides(size = guide_legend(order = 1), fill = guide_legend(order = 2)) +
  theme(text=element_text(size=12, family="TakaoExMincho"))
p<- p + geom_point(data=nujap,aes(x=longitude, y=latitude),shape=24,size =5, fill=rgb(0,0,0,0.8),color="black")
p<- p + geom_text(data=nujap[i,],aes(x=longitude, y=latitude, label = station),vjust = 0, nudge_y = 0.02)
#30km円
circle30=as.data.frame(destPoint(c(nujap$longitude[i],nujap$latitude[i]), b=1:365, d=30000))
p<-p + geom_polygon(data=circle30, aes(x=lon,y=lat), color="black",fill="red",alpha=0.2)
p<-p + geom_text(aes(x=nujap$longitude[i],y=max(circle30$lat)),label="30km", nudge_y = 0.02)
#50km円
circle50=as.data.frame(destPoint(c(nujap$longitude[i],nujap$latitude[i]), b=1:365, d=50000))
p<-p + geom_polygon(data=circle50, aes(x=lon,y=lat), color="black",fill="blue",alpha=0.2)
p<-p + geom_text(aes(x=nujap$longitude[i],y=max(circle50$lat)),label="50km", nudge_y = 0.02)
#
png(paste0("nuclear",i,".png"),width=1600,height=600)
grid.arrange(p0,p,ncol = 2)
dev.off()
}