1. 目的

この資料では，GNSS 観測点の変位時系列データを読み込み，

東西・南北・上下の 3 成分を同時にプロットします。
2011 年 3 月 11 日の地震を基準線として表示します。
地震後の東西成分に対して対数関数フィッティングを行います。

2. 作業ディレクトリの設定とパッケージ読み込み

ここでは，

フォルダを作業ディレクトリに設定します。
必要なパッケージをすべて読み込みます。

library(this.path)
setwd(this.dir())

library(lubridate)
library(ggplot2)

3. データの読み込みと日付処理

ここでは，GNSS データファイル 950296-rx.dat を読み込みます。

1列目：日付
2列目：東−西変位 [m]
3列目：北−南変位 [m]
4列目：上−下変位 [m]

その後，日付を ymd() で Date 型に変換し，データフレームにまとめます。

nel=read.table("./950296-rx.dat", header=T)
   #1列目が日付, 2列目が東-西[m], 3列目が北-南[m], 4列目が上-下[m]
da=ymd(nel[,1])
   #年月日

edat=data.frame(da, nel[,2], nel[,3], nel[,4])

4. GNSS変位時系列の描画（3成分同時）

東西・南北・上下の 3 成分を，同一の図に重ねてプロットします。
不透明度（alpha）を小さくすることで，長期時系列の分布を可視化しています。

p2 = ggplot()

p2 = p2 + geom_point(data=edat, aes(x=edat[,1],y=edat[,2]), col="blue", alpha = 0.1, size=0.1) #alphaが不透明度
p2 = p2 + geom_point(data=edat, aes(x=edat[,1],y=edat[,3]), col="red", alpha = 0.1, size=0.1) 
p2 = p2 + geom_point(data=edat, aes(x=edat[,1],y=edat[,4]), col="black", alpha = 0.1, size=0.1) 

p2 = p2 + ylab("Relative displacement [m]") + xlab("Time [year]") 
p2 = p2 + scale_y_continuous(breaks=seq(-0.16,0.16,0.04)) 
p2 = p2 + coord_cartesian(ylim=c(-0.16,0.16))
p2 = p2 + scale_x_date(breaks="3 year", date_labels = "%Y")
p2

5. 2011年3月11日の表示と凡例の追加、図の体裁

2011年3月11日の位置に縦線を引き，
さらに東西・南北・上下成分の凡例を手動で描き込みます。
白背景・枠線つきの体裁を設定します。

p2 = p2 + geom_vline(xintercept =ymd("2011-03-11"), linewidth=0.5, linetype=3, colour="black") 

##凡例を手動プロット
h1=data.frame(ymd(19960701), 0.15)
h2=data.frame(ymd(19960701), 0.13)
h3=data.frame(ymd(19960701), 0.11)

p2 = p2 + geom_point( data=h1, aes(x=h1[,1], y=h1[,2]), col="blue", size=0.5)
p2 = p2 + geom_point( data=h2, aes(x=h2[,1], y=h2[,2]), col="red", size=0.5)
p2 = p2 + geom_point( data=h3, aes(x=h3[,1], y=h3[,2]), col="black", size=0.5)

p2 = p2 + annotate("text", x=ymd(19970701), y=0.15, label="East-West", size=3, hjust=0)
p2 = p2 + annotate("text", x=ymd(19970701), y=0.13, label="North-South", size=3, hjust=0)
p2 = p2 + annotate("text", x=ymd(19970701), y=0.11, label="Up-Down", size=3, hjust=0)

p2 = p2 + theme( panel.background = element_rect(fill="white", colour = "black", linetype = "solid"),
                 panel.grid = element_blank())   #白地の普通のグラフを書く体裁
p2

6. PNG保存

その後 PNG ファイルとして保存します。
図のサイズや解像度を自由に指定できます。

ggsave(p2,filename="gnss.png",width=4,height=3,dpi=400)

7. 対数関数による余効変動フィッティング

地震後（2011-03-11 以降）の東西成分に対して，
対数関数

\[ x(t) = a\log(t-t_0)+b \]

を仮定してフィッティングを行います。
観測データに対して，対数フィット曲線を重ねて描画します。

#関数fitting
library(minpack.lm)

after=subset(edat, edat[,1]> ymd(20110311))
xa=after[,1]; ya=after[,2]

logfit=nlsLM( ya ~ a*log(as.numeric(xa-ymd(20110311)))+b, start = c(a=1,b=1))  # start= は推定の初期値

a2=summary(logfit)$parameters[[1]]
b2=summary(logfit)$parameters[[2]]

result=data.frame(xa, a2*log(as.numeric(xa-ymd(20110311)))+b2)

p2 = p2 + geom_line( data=result, aes(x=result[,1], y=result[,2]), col="green", size=0.5)
p2

8. png保存

最終図を PNG として保存します。

ggsave(p2,filename="gnss_logfit.png",width=4,height=3,dpi=400)

9. まとめ

GNSS 変位の 3 成分を同時に可視化しました。
2011年東北地方太平洋沖地震の余効変動を対数関数で近似し、図示しました。
可視化 → フィッティング → 再可視化の一連の流れを確認できました。

ggplot2演習: GNSS変位時系列と対数フィッティングの可視化（HTML版）

Y. Mitsui

2025-12-06