ggplot2

地震動のスペクトル ggplot2 (Hanningウィンドウ)

Rcpp , ggplot2 , ggfortify , scales , reshape2 , gridExtra パッケージ

(過去の記事)
地震動のスペクトル1(Hanningウィンドウ)
地震動のスペクトル2(Hanningウィンドウ)

前半は過去の記事と同じ。作図にggplot2 パッケージを使った。

Hanning Window の関数を定義

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
library(Rcpp)
cppFunction('
            NumericVector hanningRcpp(List x) 
{          
            NumericVector acc=x["acc"];
            int nhan=x["nhan"];
            int nd = acc.size();
            NumericVector ret(nd);
            NumericVector temp=clone(acc);
            for(int ihan=0; ihan < nhan; ihan++){
            ret[0]=0.5*(temp[0]+temp[1]); 
            ret[nd-1]=0.25*temp[nd-2]+0.5*temp[nd-1]+0.25*temp[0]; 
            for(int i=1; i < nd-1; i++){
            ret[i]=0.25*temp[i-1]+0.5*temp[i]+0.25*temp[i+1]; 
            }
            temp=clone(ret);
            }
            return ret;
}
            ')

気象庁のサイトから鳥取県西部地震の境港市東本町のデータを取り込み、処理。

windowのバンド幅 bd は1 とした。
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
#入力
#加速度　境港市東本町
seibu1<-read.csv("http://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/001006_tottori-seibu/dat/AA06E9E1.csv",skip=6,header=T)
#サンプリング周期
dt     <- 0.01
#windowのバンド幅
bd<-1
for (direction in 1:3){
if (direction==1) {d<-seibu1$NS}
if (direction==2) {d<-seibu1$EW}
if (direction==3) {d<-seibu1$UD}
    k  <- length(d)
#入力データ数が２の乗数になるように調整
    nn <- 2
    while( nn < k) nn <- nn*2
    x      <- numeric(nn)
    x[1:k] <- d[1:k]
    z      <- fft(x)
    z1     <- Re(z)/nn
    z2     <- Im(z)/nn
    nfold  <- nn/2+1
    df     <- 1.0/dt/nn
    fs1    <- numeric(nfold)
    fs2    <- numeric(nfold)
    frq    <- numeric(nfold)
    fs1    <- dt*nn*sqrt(z1[1:nfold]^2+z2[1:nfold]^2)
    ii     <- seq(0,nfold-1,by=1)
    frq    <- ii*df
#ウィンドウのバンド幅 bdのときの繰り返し回数nを計算 
hann<-(bd*(3*(length(x)*dt)/8))^2
n=ceiling(hann)
band<-round(8/3*sqrt(n)*(1/(length(x)*dt)),3)
fs2<-hanningRcpp(list(acc=fs1,nhan=n)) 
if (direction==1) {resNS<-data.frame(frq,fs1,fs2)}
if (direction==2) {resEW<-data.frame(frq,fs1,fs2)}
if (direction==3) {resUD<-data.frame(frq,fs1,fs2)}
}

ggplot2 パッケージを用いて作図

グラフの表示順は、 EW , NS , UD になっています。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
library(ggplot2)
library(ggfortify)
library(scales)
library(reshape2)
library(gridExtra)
#
#波形を見る(ts class)
seibu1.ts<-ts(seibu1,start=0,freq=100)
lim=max(abs(max(seibu1.ts)),abs(min(seibu1.ts)))
#[1] 748.396
#
#png("quakeggplot01.png")
autoplot(seibu1.ts,size=0.2)+
ylim(-lim,lim) +
ggtitle("波形（加速度）")
#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
lim=max(max(resEW$fs1),max(resNS$fs1),max(resUD$fs1))
#EW
x1<-reshape2::melt(resEW,id.vars="frq")
p1<-ggplot(data=x1,aes(frq,value,colour=variable)) + geom_line() 
p1<-p1 + scale_colour_manual(name="",labels = c(fs1= "spectrum",fs2  = "Hanning"),values = c(fs1= "tomato",fs2  = "gray20") )
p1<-p1 + scale_y_log10(labels=trans_format("log10", math_format(10^.x))) 
p1<-p1 + coord_cartesian(xlim=c(0,20), ylim=c(1e-2,lim)) + annotation_logticks(sides="l",colour="white")  
p1<-p1 + labs(x="Frequency (Hz)",y="",title="Fourier spectrum (gal * sec) EW")
#NS
x1<-reshape2::melt(resNS,id.vars="frq")
p2<-ggplot(data=x1,aes(frq,value,colour=variable)) + geom_line() 
p2<-p2 + scale_colour_manual(name="",labels = c(fs1= "spectrum",fs2  = "Hanning"),values = c(fs1= "tomato",fs2  = "gray20") )
p2<-p2 + scale_y_log10(labels=trans_format("log10", math_format(10^.x))) 
p2<-p2 + coord_cartesian(xlim=c(0,20), ylim=c(1e-2,lim)) + annotation_logticks(sides="l",colour="white")  
p2<-p2 + labs(x="Frequency (Hz)",y="",title="Fourier spectrum (gal * sec) NS")
#UD
x1<-reshape2::melt(resUD,id.vars="frq")
p3<-ggplot(data=x1,aes(frq,value,colour=variable)) + geom_line() 
p3<-p3 + scale_colour_manual(name="",labels = c(fs1= "spectrum",fs2  = "Hanning"),values = c(fs1= "tomato",fs2  = "gray20") )
p3<-p3 + scale_y_log10(labels=trans_format("log10", math_format(10^.x))) 
p3<-p3 + coord_cartesian(xlim=c(0,20), ylim=c(1e-2,lim)) + annotation_logticks(sides="l",colour="white") 
p3<-p3 + labs(x="Frequency (Hz)",y="",title="Fourier spectrum (gal * sec) UD")
#
#png("quakeggplot02.png")
grid.arrange(p1,p2,p3,ncol=1)
#dev.off()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
lim=max(max(resEW$fs2),max(resNS$fs2),max(resUD$fs2))
temp<-resEW[-1,c("frq","fs2")]
names(temp)<-c("frq","EW")
x1<-reshape2::melt(temp,id.vars="frq")
temp<-resNS[-1,c("frq","fs2")]
names(temp)<-c("frq","NS")
x1<-rbind(x1,reshape2::melt(temp,id.vars="frq"))
temp<-resUD[-1,c("frq","fs2")]
names(temp)<-c("frq","UD")
x1<-rbind(x1,reshape2::melt(temp,id.vars="frq"))
#
p1<-ggplot(x1,aes(1/frq,value,colour=variable)) + geom_line()
p1<-p1 + scale_x_log10() + scale_y_log10()
p1<-p1 + coord_cartesian(xlim=c(0.03,10), ylim=c(0.1,lim)) + annotation_logticks(sides="bl",colour="white")  
p1<-p1 + scale_colour_discrete(name="Direction")
p1<-p1 + labs(x="Period(sec)",y="",title="Fourier spectrum (gal * sec)")
#png("quakeggplot03.png")
p1
#dev.off()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
grdx<- as.vector((2:9) %o% 10^(-2:1))
grdy<- as.vector((2:9) %o% 10^(-1:3))
p1<-ggplot(x1,aes(1/frq,value,colour=variable)) + geom_line()
#grid lineの色、太さ変更する場合
#p1<-p1 + theme(panel.grid.minor = element_line(colour="white", size=0.2))
p1<-p1 + scale_x_log10(minor_breaks = grdx) + scale_y_log10(breaks =c(1,10,100),minor_breaks = grdy)
p1<-p1 + coord_cartesian(xlim=c(0.03,10), ylim=c(0.1,lim))
p1<-p1 + scale_colour_discrete(name="Direction")
p1<-p1 + labs(x="Period(sec)",y="",title="Fourier spectrum (gal * sec)")
#png("quakeggplot04.png")
p1
#dev.off()
#(x,y両軸の補助線をかかない)theme_classic の補助線をかく。
library(cowplot)
#png("quakeggplot05.png")
p1<-p1 + theme_classic()  + background_grid(major = "xy", minor = "xy")
p1
#dev.off()

cowplot パッケージ

Introduction to cowplot
Plot annotations