スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Rによる3Dプロット(RGLパッケージ使用)

Rのrglライブラリでヒストグラムを作成する方法がわかった。本ページにて紹介する方法で、下図のような立体のヒストグラムを描くことができる。

iceT010_convert_20101211190016.png


カラーマップに比べ、z軸方向で表される値の高低に注意を向かわせる表現をしたい場合には優れたプロット方法であろう。


次の時系列データがあるとして、

1000000 3912.882090 325.638341 2583.000700
1000010 3938.663450 328.008606 2571.402200
1000020 4005.116200 325.475949 2529.966400
1000030 3817.944600 343.527834 2486.984100
1000040 3822.547790 330.767038 2454.971200
1000050 3836.393540 329.273641 2477.533000
1000060 3827.740110 328.972475 2393.342400
1000070 3953.690080 337.964218 2423.753900
1000080 3945.989290 337.305724 2475.190800
1000090 3897.403140 347.804742 2488.807900
1000100 3867.529670 325.866540 2447.205200
1000110 3797.758710 339.163194 2487.733500
1000120 3777.881240 344.835873 2489.342800
1000130 3790.815560 329.731020 2496.429000
1000140 3885.044590 333.818744 2473.586900
1000150 3847.940810 336.612753 2435.929300
1000160 3784.553010 338.839554 2495.632500
1000170 3779.970170 332.510576 2438.425500
1000180 3782.549810 341.432855 2428.608000
1000190 3854.629460 336.316134 2427.246100
1000200 3817.810030 331.365585 2474.987800
1000210 3887.884110 329.141570 2496.600600
1000220 3760.931430 325.840837 2470.318600
1000230 3843.553530 331.861819 2413.599900
1000240 3771.818380 356.066626 2449.877400
1000250 3867.856270 332.201465 2444.941900
             ・
             ・
             ・



上記データファイルを"histdata.dat"とする。

第2列と第3列を使い、時系列頻度を二次元座標中に表現することが目標である。
gregmiscライブラリの「hist2d」を、前回と同じく使用し、この関数から出力されるデータをrglの関数に渡して三次元プロットを行う。

まず、以下を実行する:

library(gregmisc)
data1 <- read.table("histdata.dat")
datax <- data1[,2]
datay <- data1[,4]
h2d <- hist2d(datax, datay, show=FALSE, same.scale=TRUE, nbins = c(800,800))
zzz <- h2d$counts
xxx <- h2d$x
yyy <- h2d$y


今回使う方法では、扱うデータは行列形式でなければならず、volcanoと同様のものでなければならない。すなわち、行と列で座標を表現し、要素の値でz軸方向を指定するというタイプのデータである(もしかしたらそうじゃないかもしれないが、自分は他の方法を知らない。自分も初心者なのでわかりません)。

hist2d()から得られるh2dというデータから$countsを抜き出してzzzというデータに渡しているが、これが行列形式になっている。だが、これでは範囲が広くとられているはず。そして、自分が知る限りでは、本方法で行うプロット方法にはxlim、ylim、あるいはzlimのような機能はなく、行列データで指定されるすべての範囲を漏れなくプロットすることになるため、必要に応じて行列を削る必要がある。その操作が以下である(削る範囲は問題・目的に応じて変える必要があります)。

元の行列は800×800である。これを上下左右削り、220×620にする。

zzz2 <- tail(zzz, n=380)
xxx2 <- tail(xxx, n=380)
zzz3 <- head(zzz2, n=200)
xxx3 <- head(xxx2, n=200)
zzzt <- t(zzz3)
zzzt2 <- head(zzzt, n=180)
yyy2 <- head(yyy,n=180)
zzz2 <- t(zzzt2)


データ整形が完了した。以下、rglによるプロットを行う手順を示す。

library(rgl)
zzz2[1,1] <- 2500 # 最低温の領域でピークが最も高くなる場所の頂点が2500付近なので2500と設定した。colourのみ効いて、z軸方向の範囲指定には効かない
zlim <-c(0,2500)
zlen <- zlim[2] - zlim[1] + 1
colorlut <- cm.colors(zlen) # height color lookup table
col <- colorlut[(zzz2-zlim[1]+1)] # assign colors to heights for each point
open3d()
surface3d(xxx3/150.0, yyy2/100.0, zzz2/1000.0, color=col,alpha=0.75, back="lines")
axes3d(c('x','y','z-'))
grid3d(c("x", "y+", "z"),at = NULL,col=c("gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray"), n =8,lwd = 1, lty = "solid")
title.name <- paste("T=", "1.00", sep = "")
title3d(main = title.name, sub = "", xlab = "radius of gyration", ylab = "Sigma theta", zlab = "")
rgl.bringtotop()


上のようにすることで、rglのウィンドウが出てきて、マウスで動かせるようになったと思う(本ページ後半部分は、慣れないうちに筆者が試行錯誤し苦悶していたときの履歴が書かれているのでそれも参考になると思います)

また、「rgl.quit()」でrglのウィンドウは閉じられる。

pngやeps形式で表示されている画像を保存する場合は次のようにする:

rgl.snapshot("iceT010.png")
rgl.postscript("iceT010.eps", fmt="eps" )



しかし、rglウィンドウは、マウスで操作しないと角度が調整できないので、同じ角度、倍率から異なるデータ由来のプロット画像を連続して得たい場合、その角度・倍率・位置をファイルに出力する関数があるので、それをまずいったんは使い、

save(pv, file="pv.txt", ascii=T)



その後次のようにすると同じ視点からプロットが出てくる(少しウィンドウを動かすと反映される):

load("pv.txt")
par3d(windowRect=pv$windowRect) #rglのWindowの表示位置、大きさの再現
par3d(zoom=pv$zoom) #ズームの再現
par3d(userMatrix=pv$userMatrix) #視点およびその他の再現



この状態で先ほどの画像出力のコマンドを入力すれば画像が得られることになる。


参考にした主なページ:
rglで表示方法を保存、再現する
Experience with R
メモ帳




(以下、不慣れな時期の自分のメモ)
------

Rには「volcano」というデータがサンプルとして入っているが、これを最も美麗にプロットしていると思ったのが、「RGLパッケージ」を使った3Dプロットだ。

以下のページを参考にした:
http://ja.efreedom.com/Question/1-1896419/3-D-%E3%81%AE%E8%A1%A8%E9%9D%A2%E3%83%95%E3%83%AD%E3%83%83%E3%83%88%E3%82%92%E7%AD%89%E9%AB%98%E7%B7%9A%E3%82%AA%E3%83%8F%E3%83%AC%E3%82%A4R-%E3%82%92%E4%BD%BF%E7%94%A8%E3%81%97%E3%81%A6%E3%83%95%E3%83%AD%E3%83%83%E3%83%88

install.packages("rgl")



その後は上記ページと同じコマンドを入力するだけである。

data(volcano)
z <- 2 * volcano # Exaggerate the relief
x <- 10 * (1:nrow(z)) # 10 meter spacing (S to N)
y <- 10 * (1:ncol(z)) # 10 meter spacing (E to W)
zlim <- range(y)
zlen <- zlim[2] - zlim[1] + 1
colorlut <- terrain.colors(zlen,alpha=0) # height color lookup table
col <- colorlut[ z-zlim[1]+1 ] # assign colors to heights for each point
open3d()
rgl.surface(x, y, z, color=col, alpha=0.75, back="lines")



結果、以下が出力されるが、マウスポインタで自由に回転させることができる。軸があるといいのだが。そして、下に等高線のプロットがあるとなお良い(このページで示されているFig.13.8のように)。

rgl_volcano1.png
rgl_volcano2.png


gridlineの設定は実線(solid line )とした。変更の場合はltyを編集する。パラメータとの対応はR-Source参照

以下を参考に、下図を作るまでには至った(等高線は無い方が良いと結論)。あとは実際に使うデータが手元にないと本番用の図が描けない。月曜に作成する。

参考サイト
Experience with R
メモ帳

rgl_volcano3_convert_20101204143718.png

コード(あまりに汚いので、あとで編集しておくこと)

library(rgl);
data(volcano)
dim(volcano)

z <- 2*volcano;
x <- 10*(1:nrow(volcano));
y <- 10*(1:ncol(volcano));

zlim <- range(z)
zlen <- zlim[2] - zlim[1] + 1
colorlut <- rainbow(zlen) # height color lookup table
col <- colorlut[(z-zlim[1]+1)] # assign colors to heights for each point
open3d()

title.name <- paste("free energy ", "landscape", sep = "");
surface3d(x, y, z, color=col,alpha=0.75, back="lines", main = title.name);
grid3d(c("x", "y+", "z"),at = NULL,col=c("gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray","gray"), n =8,lwd = 1, lty = "solid")

sample.name <- paste("col.", 1:ncol(volcano), sep="");
sample.label <- as.integer(seq(1, length(sample.name), length = 5));

axis3d('y+',at = y[sample.label], x[sample.label], cex = 0.3,col="#525252");
axis3d('y',at = y[sample.label], x[sample.label], cex = 0.3,col="#525252")
axis3d('z',pos=c(0, 0, NA),col="#525252")

ppm.label <- as.integer(seq(1, length(x), length = 10));
axes3d('x', at=c(x[ppm.label], 0, 0), abs(round(y[ppm.label], 2)), cex = 0.3);

title3d(main = title.name, sub = "subtitle", xlab = "x", ylab = "y", zlab = "z")
rgl.bringtotop();




rglで描画した場合、見る角度等を再現することは困難である。解決方法が下記のリンク先にて紹介されている。

rglで表示方法を保存、再現する
スポンサーサイト
Profile

zoa

Author:zoa
Metropolis (1927)

Calendar
11 | 2010/12 | 01
- - - 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 -
Labels
Admin
Previous Posts
Recent Comments
Recent Trackback
Archive
このページのトップへ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。