周期構造の可視化
概要
ここでは，地点データを用いた時空間特製の可視化で表示している図について，データの処理から 可視化までを説明していく．データの説明については， アメダスの気象データの説明ページを見てもらいたい． まずは，可視化については，3次元で行うことを考える．次に，アメダスのデータの中から3種類のデータを選ぶ． ここでは「降雨量・平均気温・平均風速」，「降雨量・最高気温・最低気温」の２パターンを採用した． まず，各データの周期性を推定するために スプライン回帰を用いる． 次に，観測期間で一年の各日についてどのくらいばらつきがあるかを可視化するため，各日のデータだけを抜き出して 主成分分析を 用いて3次元に落とし込み，固有ベクトルを用いて各日のばらつきを楕円で可視化する． また，スプライン回帰によって推定した周期構造について，曲率(curvature)と捩率(torsion)による色付けによる比較も行う．
周期スプラインによる周期性の推定
ここからは，概要で述べたうち全観測所における「降雨量・平均気温・平均風速」を扱っていく．ここでは年周期に興味を絞り，データのうち一年の各日で平均をとり，各日の平均値を それぞれのデータごとで周期スプラインで回帰する．つまり，降雨量，平均気温，平均風速の一年の各日の平均値をそれぞれ\( {\bar y^{(rain)}_i}, {\bar{y}^{(temp)}_i}, {\bar{y}^{(wind)}_i} \)， \(t=1,\ldots,366\)として\(t\)で回帰する．ここで1月1日を1とし，閏日を含め12月31日を366とした．数式で表すとそれぞれ
という関係となる．以下にRのコードとスプラインによって回帰された結果を図１示す．
## amedas.html に使う図(日本地図と観測所以外)を作成する
data_path <- "データが保存されているフォルダへのパス"
## packages =======
library(tidyverse)
library(data.table)
#===================
## import =========
rain <- data.table::fread(file.path(data_path, "Amedas.csv"))
aveTemp <- data.table::fread(file.path(data_path, "aveTemp.csv"))/10
aveVelo <- data.table::fread(file.path(data_path, "aveVelo.csv"))/10
maxTemp <- data.table::fread(file.path(data_path, "maxTemp.csv"))/10
minTemp <- data.table::fread(file.path(data_path, "minTemp.csv"))/10
maxVelo <- data.table::fread(file.path(data_path, "maxVelo.csv"))
Data_list <- list(rain=rain, aveTemp=aveTemp, maxTemp=maxTemp,
minTemp=minTemp, aveVelo=aveVelo, maxVelo=maxVelo)
## 日毎の平均をとる
res <- Data_List[c(1,4,6)] %>% lapply(function(L1){
L1[,-c(1:3)] %>% rowMeans(na.rm=TRUE)
}) %>% data.frame()
res <- data.frame(Data_List[[1]][,1:3], res)
#----------------
#----------------
## 一年の各日で平均をとり周期スプラインで関数推定を行う
## smoothing_calcは著者自作の関数 githubよりダウンロード可能
smooth_res <- smoothing_calc(Data=res, colnum = 4:6, type="spline", leap=TRUE)
my_col <- ranbow(366)
v_lines <- cumsum(c(1,31,28,31,30,31,30,31,31,30,31,30,31))
par(mfrow=c(3,1), mar=c(2,5,1,1))
plot(smooth_res[,3], type="b", ylab="rain", cex.lab=1.5, lwd=1.5, col=my_col, xaxt="n")
abline(v=v_lines, lty="dashed", col="grey")
axis(side=1, at=v_lines, labels=1:12)
plot(smooth_res[,4], type="b", ylab="temp", cex.lab=1.5, lwd=1.5, col=my_col, xaxt="n")
abline(v=v_lines, lty="dashed", col="grey")
axis(side=1, at=v_lines, labels=1:12)
plot(smooth_res[,5], type="b", ylab="wind", cex.lab=1.5, lwd=1.5, col=my_col, xaxt="n")
abline(v=v_lines, lty="dashed", col="grey")
axis(side=1, at=v_lines, labels=1:12)
ここで，1~366の違いがわかりやすいようにRのrainbow()関数によってでカラーリングしている．
各日の分散の推定
次に，一年の各日のデータのばらつきを可視化する．手順としては，各日における3つの項目のデータを抽出し，主成分分析を行い，その第1,第2主成分ベクトルを短軸と長軸に持ち， 中心が先ほどのスプラインの推定の結果の点に持つ楕円を描く．これの処理についても著者自作の関数によって行なっているが，実際にはRではデフォルトで組み込まれている prcomp()関数を利用している．
## 各日の分散計算
x1A <- PCA_each_day(Data=res, colnum=4:6, leap=TRUE)
描画
最後に，計算結果を描画する．描画にはrglパッケージのplot3dを利用する．ただし，下記の二つの関数はどちらも著者自作の関数である ので，ダウンロードしてから利用してほしい．Rのコードと結果を示そう（図２）．
## 描画に必要な要素を計算する関数
calcA <- calc.smooth(smoothing_calc_result=xA, PCA_each_day=x1A, num=v_lines, pcol="day", tms.prm = 0.5)
## 描画する関数
plot.smooth(calc.smooth=calcA, num=num, col=calcA$mycol, addAxis=FALSE, size=2,
xlab="Rain", ylab="Ave Temp", zlab="Ave Wind", PCA=FALSE)
曲率と捩率
最後に，3次元の輪に対する曲率と捩率について紹介する．いま，3次元空間の曲線を\(x_t, y_t, z_t, \ (t=1,\ldots,T)\)と表し，各点における 曲率\(\kappa_t\)と捩率\(\tau_t\)を次のように定義する．
曲率\(\kappa_t\)は，局所的な曲がり具合を表し，曲率半径の逆数であるので非負である．曲率が大きいほど急な曲がり方をしており，小さいほど緩やかなカーブを 描く．捩率はどのくらいねじれているかを表し，曲率とは異なり負の値もとる．捩率\(\tau_t\)とは曲線と平面との関係を表しており，曲線が平面に埋め込まれている場合は 0を取る．つまり2次元上の曲線に対して捩率は常に0である．捩率が0でない場合，曲線が平面から逸れていく状況を表しており，値が大きいほどその具合は顕著になる． また，捩率の正負は，右手系の座標においてそれぞれ左回りの螺旋，右回りの螺旋を描くかを意味している． 以下の図では，曲率については対数を取ったものを赤のグラデーションで色をつけている．赤が濃くなるほど曲率が大きいことを示す． 捩率については，正の値が赤，0を黒，負の値を青としてグラデーションをつけている．Rのコードと結果を図３に示す．
## 曲率と捩率
#-------------------------------
# curvatureのlogスケールのカラーリング
curv_logcolRule <- function(x, colorsp=c("red2","red1","red"), colorsm=c("red3", "red4", "black")){
# colorRamp for plus value
color_rulep <- colorRamp(colorsp)
# colorRamp for minus value
color_rulem <- colorRamp(colorsm)
# -1~1にスケール変換
tmp1 <- max(abs(x))
x2 <- x/tmp1
plus_num <- which(x2 >= 0)
minus_num <- which(x2 < 0)
res <- c()
res[plus_num] <- rgb(color_rulep(x2[plus_num])/255)
res[minus_num] <- rgb(color_rulem(-x2[minus_num])/255)
return(res)
}
#-------------------------------
#-------------------------------
# torsionのスケールのカラーリング
tors_colRule <- function(x, colorsp=c("red2","red1","red","red3", "red4", "black"), colorsm=c("royalblue1","royalblue2","royalblue","royalblue3","royalblue4","black")){
plus_num <- which(x >= 0)
minus_num <- which(x < 0)
log_x <- my_log(x)
lxp <- log_x[plus_num]
lxm <- log_x[minus_num]
res <- c()
res[plus_num] <- curv_logcolRule(lxp, colorsp=colorsp[1:3], colorsm=colorsp[4:6])
res[minus_num] <- curv_logcolRule(lxm, colorsp=colorsm[1:3], colorsm=colorsm[4:6])
return(res)
}
#-------------------------------
## 曲率と捩率の計算
res_CT <- mycalc_3dCT(smooth_res[,-c(1,2)])
col_curv <- curv_logcolRule(log(res_CT$curvature))
col_tors <- tors_logcolRule(res_CT$torsion)
laymat <- matrix(c(1,2,3,4), 2, 2, byrow=TRUE)
open3d()
layout3d(laymat, sharedMouse = TRUE, heights = c(1,10))
text3d(0,0,0, texts="curvation"); next3d()
text3d(0,0,0, texts="torsion"); next3d()
plot3d(smooth_res[,3:5], col=col_curv, xlab="rain", ylab="temp", zlab="wind")
plot3d(smooth_res[,3:5], col=col_tors, xlab="rain", ylab="temp", zlab="wind")