Nov 20 2016

震源地と原子力発電所02

raster , marmap , rgdal , xts パッケージ

今回は日本に絞って図示。
震源データは気象庁。マグニチュード6.5以上のデータ。

日本で最初の原子力発電が行われた1963年10月26日以前の震源データで作成した図

load("in_service.dat")
#load("under_construction.dat")
#in_service[in_service$Country=="Japan",]
#under_construction[under_construction$Country=="Japan",]
#
load("eq61923_20161031.RData")
eqdata<-eq6
#震源の深さ
# 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,]

library(marmap) 
library(rgdal)
library(raster)
Lon.range = c(122, 154)
Lat.range = c(20, 46)
# 最初のみネットにアクセス　keep=TRUE で保存する
#dat<-getNOAA.bathy(Lon.range[1],Lon.range[2],Lat.range[1],Lat.range[2],res=2,keep=TRUE)
# ２度めからは保存したデータを読み込む
dat<-read.bathy("marmap_coord_122;20;154;46_res_2.csv", header=TRUE) 
r1<-marmap::as.raster(dat)
#
oce<-colorRampPalette(c(rgb(0,0,0),rgb(0,5/255,25/255),rgb(0,10/255,50/255),rgb(0,80/255,125/255),rgb(0,150/255,200/255),
  rgb(86/255,197/255,184/255),rgb(172/255,245/255,168/255),rgb(211/255,250/255,211/255),rgb(1,1,1)))
#
land1 <- colorRampPalette(c("#467832","#786432"))
land2 <- colorRampPalette(c("#786432","#927E3C"))
land3 <- colorRampPalette(c("#927E3C","#C6B250"))
land4 <- colorRampPalette(c("#C6B250","#FAE664"))
land5 <- colorRampPalette(c("#FAE664","#FAEA7E"))
breakpoints <- c(seq(-10000,0,10000/100),1  ,seq(50,500,50),seq(550,1000,50),seq(1100,2000,100),seq(2100,3000,100),seq(3500,8000,500))
colors <- c(oce(100),rgb(70/255,120/255,50/255),land1(10),land2(10),land3(10),land4(10),land5(10))
#
#png("japan65_01.png",width=800,height=600)
par(family="TakaoExMincho") 
plot(r1,breaks=breakpoints,col=colors, axis.args=list(at=c(-8000,-4000,0,4000,8000), labels=c(-8000,-4000,0,4000,8000)),
   legend.args=list(text="  height (m)", side=3, font=2, line=1.5, cex=0.8),las=1,axes=F,box=F)
par(xpd=T)
axes1<-seq(125,150,5)
axes2<-seq(20,45,5)
text(axes1,Lat.range[1],as.character(axes1),pos=1)
text(Lon.range[1],axes2,as.character(axes2),pos=2)
rect(Lon.range[1],Lat.range[1],Lon.range[2],Lat.range[2])
par(xpd=F)
#
#マグニチュード6.5以上
eq1=eq[eq$longitude>=122 & eq$longitude<=154 & eq$latitude>=20 & eq$latitude<=46 & eq$mag>=6.5,]
points(x=eq1$longitude, y=eq1$latitude,cex=eq1$mag/6,pch=21,bg=eq1$col,col="gray20")
#
in_service1=in_service[in_service$longitude>=122 & in_service$longitude<=154 & in_service$latitude>=20 & in_service$latitude<=46,]
points(in_service1$longitude,in_service1$latitude,pch=24,cex =2, bg=rgb(0,0,0,0.2),col="black")
#points(under_construction$longitude,under_construction$latitude,pch=24,cex =1, bg=rgb(0,0,0,0.2),col="black")
title("震源分布（マグニチュード6.5以上）1923 ～ 2016/10/31[データ：気象庁]","(地図上の三角 : 原子力発電所)")
#dev.off()

#library(xts)
#eq.xts<-xts(read.zoo(eq))
#eq1963.xts<-eq.xts["::1963-10-26"]
#eq1963<-data.frame(time=index(eq1963.xts), #longitude=as.numeric(coredata(eq1963.xts)[,1]), 
#  latitude=as.numeric(coredata(eq1963.xts)[,2]), #depth=as.numeric(coredata(eq1963.xts)[,3]),
#  mag=as.numeric(coredata(eq1963.xts)[,4]), #depth_rank=coredata(eq1963.xts)[,5],
#  col=coredata(eq1963.xts)[,6])
#eq1963$depth_rank<-as.character(eq1963$depth_rank)
#eq1963$col<-as.character(eq1963$col)

1963-10-26以前のデータを抽出するだけなので

1	eq1963<-eq[substring(eq$time,1,10) < "1963-10-27",]

#png("japan65_02.png",width=800,height=600)
par(family="TakaoExMincho") 
plot(r1,breaks=breakpoints,col=colors, axis.args=list(at=c(-8000,-4000,0,4000,8000), labels=c(-8000,-4000,0,4000,8000)),
   legend.args=list(text="  height (m)", side=3, font=2, line=1.5, cex=0.8),las=1,axes=F,box=F)
par(xpd=T)
axes1<-seq(125,150,5)
axes2<-seq(20,45,5)
text(axes1,Lat.range[1],as.character(axes1),pos=1)
text(Lon.range[1],axes2,as.character(axes2),pos=2)
rect(Lon.range[1],Lat.range[1],Lon.range[2],Lat.range[2])
par(xpd=F)
#
eq1=eq1963[eq1963$longitude>=122 & eq1963$longitude<=154 & eq1963$latitude>=20 & eq1963$latitude<=46 & eq1963$mag>=6.5,]
points(x=eq1$longitude, y=eq1$latitude,cex=eq1$mag/6,pch=21,bg=eq1$col,col="gray20")
#
in_service1=in_service[in_service$longitude>=122 & in_service$longitude<=154 & in_service$latitude>=20 & in_service$latitude<=46,]
points(in_service1$longitude,in_service1$latitude,pch=24,cex =2, bg=rgb(0,0,0,0.2),col="black")
#
title("震源分布（マグニチュード6.5以上）1923 ～ 1963/10/26 [データ：気象庁]","(地図上の三角 : 原子力発電所)")
#dev.off()

mapdata と　ggplot2 を使ったシンプルな地図

日本で最初の原子力発電が行われた1963年10月26日以前の震源データで作成した図

library("ggplot2")
library("maps")
library("mapdata")
#library("geosphere")
# plot limits
Lon.range = c(122, 154)
Lat.range = c(20, 46)
#
japan = map_data("japan")
#
p<-ggplot() +
    borders("japan",fill="lightgray", size = 0.3) +
    theme_bw() +
    coord_map(xlim=Lon.range ,ylim=Lat.range) +
    labs(y="",x="")
#マグニチュード6.5以上
eq1=eq[ eq$mag>=6.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(6.5,7,8,9),range = c(1,4)) +
# 凡例のタイトルの指定
  labs(size="Magnitude",fill="Depth",title="震源分布（マグニチュード6.5以上）1923 ～ 2016/10/31[データ：気象庁]") +
#凡例の順番を入れ替える
  guides(size = guide_legend(order = 1), fill = guide_legend(order = 2)) +
  theme(text=element_text(size=12, family="TakaoExMincho"))
#png("japan65_03.png",width=800,height=600)
 in_service1=in_service[in_service$Country=="Japan",]
 p + geom_point(data=in_service1,aes(x=longitude, y=latitude),shape=24,size =5, fill=rgb(0,0,0,0.3),color="black")
#dev.off()

p<-ggplot() +
    borders("japan",fill="lightgray", size = 0.3) +
    theme_bw() +
    coord_map(xlim=Lon.range ,ylim=Lat.range) +
    labs(y="",x="")
#マグニチュード6.5以上
eq1=eq1963[ eq1963$mag>=6.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(6.5,7,8,9),range = c(1,4)) +
# 凡例のタイトルの指定
  labs(size="Magnitude",fill="Depth",title="震源分布（マグニチュード6.5以上）1923 ～ 1963/10/26[データ：気象庁]") +
#凡例の順番を入れ替える
  guides(size = guide_legend(order = 1), fill = guide_legend(order = 2)) +
  theme(text=element_text(size=12, family="TakaoExMincho"))
#png("japan65_04.png",width=800,height=600)
 in_service1=in_service[in_service$Country=="Japan",]
 p + geom_point(data=in_service1,aes(x=longitude, y=latitude),shape=24,size =5, fill=rgb(0,0,0,0.3),color="black")
#dev.off()

raster , marmap , rgdal , xts パッケージ

日本で最初の原子力発電が行われた1963年10月26日以前の震源データで作成した図

1963-10-26以前のデータを抽出するだけなので

mapdata と ggplot2 を使ったシンプルな地図

日本で最初の原子力発電が行われた1963年10月26日以前の震源データで作成した図

mapdata と　ggplot2 を使ったシンプルな地図