スポンサーサイト

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

2次元ヒストグラムの作成(ggplot2使用)

時系列データを取得したとき、2次元のヒストグラムを描きたい場合がある、が、gnuplotは不格好だと感じたため、別の方法を探すことにした。
結果、Rでヒストグラムを作るという判断が妥当と考えたので、方法を探したのだが、「カラーバー」を表示したいのだけど、それに困る。contour()関数を改造して使い二次元プロットすることもできるが(方法はこちらで紹介されている)、カラーバーが何となく不細工(特に線の入り方)と思ったから、また更に別の方法を考えることにした。

Rでは高度なグラフィックスライブラリが多数開発されている。「ggplot2」というライブラリはその一つだ。今回はこれを使用する。

参考:「ランゴリアーズが来るぞ・・・!!

次のような時系列データがあって、これの第2列と第3列を使い、両列を縦・横軸にとったヒストグラムを描きたい(第1列は時間ステップ、第4列は気にしないでいい)。

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
             ・
             ・
             ・


(データ中の数字はあとで割ったりする)

上記データを編集する。gregmiscライブラリ「hist2d()」を使うが、この関数は時系列中に出てきた最大値と最小値を注目する二変数に関して記憶する。なのでこちらから時系列データを改変し、上限と下限を設定してわざとデータに組み込んでおくと、わざとプロットした範囲のヒストグラム(つまり、時系列頻度0の領域)にすべて色が付き便利なのである。

C言語でそのようなプログラムを書き、時系列データを編集しておく(これ大事)。

そして次を実行する(source()で実行した方がよい)

library(gregmisc)
library(lattice)
install.packages("gdata")
install.packages("gtools")
install.packages("gmodels")
install.packages("gplots")
install.packages("gregmisc")
install.packages("ggplot2")
library(gregmisc)
library(lattice)
library(ggplot2)


↑ これは、なぜかはわからないが、2回続けて実行しないといけない。「install.R」とファイルにして、source("install.R")を2回実行する。ライブラリ「ggplot2」もしくは「lattice」が影響しているものと思われる(latticeって要るんだったっけ?)。

先ほどの時系列データファイルの名前を「hist2d_long_4_T49.dat」と名付けているが(40万行のデータである)、R、ggplot2が扱えるデータに次で改変する。

data1 <- read.table("hist2d_long_4_T49.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
dataX <- matrix(nrow=640800,ncol=3)
for(i in 1:800){
for(j in 1:800){
dataX[i+(800*j),1] <- xxx[i]
dataX[i+(800*j),2] <- yyy[j]
if(zzz[i,j] == 0)
dataX[i+(800*j),3] <- 1
else{
if(zzz[i,j] >= 250)
dataX[i+(800*j),3] <- 250
else
dataX[i+(800*j),3] <- zzz[i,j]
}
}
}
data3 <- as.data.frame(dataX)


「gregmisc」で使えるようになるhist2d()でヒストグラムの行列データにして、そのあと250以上の頻度で時系列中に現れた格子に関しては250を代入、また、頻度0の格子は、このままではその部分の格子の色が背景になってしまうから、1を代入する操作をしている。

その後、手動で以下を実行する:

png("T008.png")
ggplot(data3, aes(x=data3[,1]/150.0, y=data3[,2]/100.0, fill=data3[,3]/400000.0)) + geom_tile() + xlim(c(23,33)) + ylim(c(21.5,28.5))+ opts(title="T=0.80") + xlab(expression(R[g]^2)) + ylab(expression(paste(Sigma,theta))) + scale_fill_gradientn("frequency",limits = c(0, 0.000625) , colours = topo.colors(200))
dev.off()


難点:
なぜかはわからないが、ggplotを使うために使ったライブラリ「ggplot2」を使うと、source()関数が使えなくなり、大量のヒストグラムを描こうとすると、一部手動でやらなければならないために労力・時間が掛かる(ヒストグラム描画の部分だけ手動でやればよい。先のソースコードの、後半の赤い文字の部分は手動で行う(Rの反応を待って一回一回入力する。何かいい方法あったら教えて))

これでも綺麗なヒストグラムにはなるのだが、色が不満である。ので、色にこだわることにした。
colorsにわたす色の配列を変えればよい。しかしグラデーションの麗美なカラーパレットをゼロから作ることは凡人には困難である。
Color Palette Generatorというサイトがあるから、これを使った。また、細微な色の調整はRGB Color Gradient Makerで行った。
パレットカラーを作り、配列に渡す。

sun <- c("#110000","#660011","#bb0000","#C93211","#D76423","#F0BB43","#ffee55","#FFF5A3","#FFFCE8")
sea <- c("#002046","#0066ff","#11aadd","#44ccee","#00ddff","#77ffff","#bbeeee","#ffffee")


「sun」と書いたのは、先のColor Palette Generatorで、太陽の衛星写真から生成したカラーパレットを参考に編集したもの。「sea」は海底から見た海面。

そしてできあがったヒストグラムが下図。
2dhistgram_sun2dhistgram_sea

これなら満足です!

他にもカラーパレットを作成してみた。

alien <- c("#000000","#447744","#559944","#55dd33","#bbcc55","#ffff77","#ffffff")
darkmoon <- c("#000011","#336699","#5599bb","#aaccdd","#EBEDD4","#F3F3E7","#eeffff")


「alien」は、映画「Alien」のポスターから。
C85_4_convert_20101201222454.pngdarkmoon


data1 <- read.table("hist2dCNMD0_T12.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
dataX <- matrix(nrow=640800,ncol=3)
for(i in 1:800){
for(j in 1:800){
dataX[i+(800*j),1] <- xxx[i]
dataX[i+(800*j),2] <- yyy[j]
if(zzz[i,j] == 0)
dataX[i+(800*j),3] <- 1
else{
if(zzz[i,j] >= 250)
dataX[i+(800*j),3] <- 250
else
dataX[i+(800*j),3] <- zzz[i,j]
}
}
}
data3 <- as.data.frame(dataX)


> source("test3.R")
> png("C85_12.png")
> ggplot(data3, aes(x=data3[,1]/150.0, y=data3[,2]/100.0, fill=data3[,3]/400000.0)) + geom_tile() + xlim(c(20,32)) + ylim(c(21.0,28.9)) + xlab(expression(R[g]^2)) + ylab(expression(paste(Sigma,theta))) + scale_fill_gradientn("frequency",limits = c(0, 0.000625) , colours = alien) + opts(title="CageSize=8.5, T=0.65")
> dev.off()
null device
1
>





下図のように、白い線が入ってしまう場合の対処法:
線が入る方向と垂直な方の軸の範囲をズラしてみると解決する場合がある。
(具体的には、「 + xlim(c(20,32))」 のところを 「 + xlim(c(20,33))」 とするなど)

erroralienalien015

> source("hist2d.R")
> png("REMDC88_120.png")
> ggplot(data3, aes(x=data3[,1]/150.0, y=data3[,2]/100.0, fill=data3[,3]/400000.0)) + geom_tile() + xlim(c(23.0,32.0)) + ylim(c(21.5,28.5)) + xlab(expression(R[g]^2)) + ylab(expression(paste(Sigma,theta))) + scale_fill_gradientn("frequency",limits = c(0, 0.000625) , colours = alien) + opts(title="CageSize=8.8, T=1.20")
> dev.off()

スポンサーサイト

コメントの投稿


非公開コメント

コメント

Profile

zoa

Author:zoa
Metropolis (1927)

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