スポンサーサイト

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

R: 大量の線を重ねてプロットする(geom_line使用)

細かい折れ線を同時に大量に描きたい場合、以下の方法が便利である:

下図のような図を描くことが目標である ↓

geomline

ライブラリ関数「gglot」を使うが、今回の場合はデータの形式をいかなるものにするかが要であろう。

コマンドは以下の通り、極めて単純である:

> ggplot(data=dataR, aes(x=num, y=value, group=gen)) + ylim(c(1.1,3.1))+
+ ylab(expression(theta)) + geom_line(size=0.02,colour="#1D2117")


この「 + geom_line()」でggplotに「折れ線で図を描く」ということを命令するのだが、ggplotの方で「group=」という指定を行うと、複数の線を描いてくれる。「aes(x=num, y=value, group=gen)」としたが、データには3軸、x、y、groupを与えれば、目的の図は描ける。

次のようなデータファイルを作成する:

gen num value
0 0.500000 1.943523
0 1.500000 1.285030
0 2.500000 1.618999
0 3.500000 2.098752
     ・
     ・
     ・
1000 0.500000 2.001602
1000 1.500000 1.232749
1000 2.500000 1.615162
1000 3.500000 2.012310
     ・
     ・
     ・
2000 0.500000 2.087536
2000 1.500000 1.195257
2000 2.500000 1.680437
2000 3.500000 2.139283
     ・
     ・
     ・


名前を「data.dat」としよう。作業フォルダに保存する。

次のコマンドで、最初の行をheaderとしてRに読み込ませ、このデータをR内の「dataR」というデータに渡す。

dataR <- read.table("data.dat",header=T)


「header=T」というオプションによりheaderを読み込ませた場合、その名前をRは記録し、また、その列にある数値を、行数を記録しつつ記録してくる。なので今の場合、gen、num、valueというこちらが勝手に指定した名前の軸を持つ3次元空間の中にある、データとして入力した点すべてを記録する。

先ほどのggplotのコマンドを入力すれば、先ほどの図が出力される。ggplotが動かない、という場合は、このブログの検索フォームに「install.R」と入れて検索すると方法が示されるだろう。


以下、管理人用の解読不能ページ

------------

手順:

1._rplotファイルから適当な場合(ΘとR2gで判断する)を抜きだし、trajectory12X.cに渡す。trajectory12X.c内で"r"ファイルとなる「recorded????.dat」を書き換えればよい。

2.trajectory12X.cのT、Cageを書き換えた上でシミュレーションを実行する。

3.trajectory_theta_12X.datにthetaの系列(このページで説明したデータファイルがこれに相当)が得られるので、R作業フォルダ内の「theta.dat」に移し、なおかつheaderを「gen re??? value」として編集する。

4.このページで紹介した方法でggplotからは出力される。その後、trajectoryのプロットをまた行わなければならない。シミュレーションを行ったフォルダ内にある「trajectory_12X.dat」の名前を適当に変更し、R作業フォルダへ移す。中身を操作する必要はない。trajectoryへ作業を移行
スポンサーサイト

出発前に

この記事を閲覧するにはパスワードが必要です
パスワード入力

R: ヒストグラム上の軌道の描画

ggplot2の「geom_path」を使うと、下図のような図が描ける。ヒストグラム中で動力学的なトラジェクトリーを示したい場合に有効な手段である。

trjC03_convert_20101231181759.png

http://had.co.nz/ggplot2/geom_path.html

上図の背景は等高線プロットである。

trajectory03_convert_20101231182107.png

ggplotの、「図を重ねて次々に描ける」という機能が発揮された図だと思う。

申し訳ありませんが、たぶんこの記事は自分以外は意味がわからないかと存じます。コマンドが参考にできる程度にしか訳には立たないでしょう。私が持っているデータ、およびデータ変換方法まで書かないと意味を持たない記事になると思います。備忘録として私も書いていますし、その点、ご了承ください。

等高線プロットの描き方:
test3.Rによりhist3d*.datを変換してヒストグラムファイルとして渡した後、生成される「data3」を使ってtilesによる二次元ヒストグラムを描く要領(本ブログの過去の記事参照)でコマンドを入力する。

ggplot(data3, aes(x=data3[,1]/150.0, y=data3[,2]/100.0, z=data3[,3]/400000.0)) + stat_contour(binwidth=0.00003, size=0.2, colour = "#280B0B") + xlim(c(23,29)) + ylim(c(22.5,27.0)) + xlab(expression(R[g]^2)) + ylab(expression(paste(Sigma,theta)))


stat_contour()が、等高線プロットを指定する関数で、binwidthにより等高線の幅を決め(ちなみにmaxとminを考え、こちらから指定するある数で割った商だけ等高線を描くという方法もある)、等高線の線の太さをsizeにより決める。

ここで、このページで行った方法によりtrajectoryのファイルができているはずなので、これをRの作業フォルダに移し、これを読み込む。それぞれ「trajectory_???.dat」という名前で保存されているはずだ。

trajectory_???.datを、描きたい軌道の数だけR内のデータ配列に渡す。

datatA <- read.table("trajectory_C03A.dat")
datatD <- read.table("trajectory_C03D.dat")
datatE <- read.table("trajectory_C03E.dat")



そして、次のコマンドを実行する:

png("trj03.png",width=450,height=450)
ggplot(data3, aes(x=data3[,1]/150.0, y=data3[,2]/100.0, z=data3[,3]/400000.0)) + stat_contour(binwidth=0.00003, size=0.2, colour = "#280B0B") + xlim(c(23,29)) + ylim(c(22.5,27.0)) + xlab(expression(R[g]^2)) + ylab(expression(paste(Sigma,theta))) + geom_path(data=datatA, aes(datatA[,1],datatA[,2]),size=0.02, color=2,alpha=0.4) + geom_path(data=datatD, aes(datatD[,1],datatD[,2]),size=0.02, color=6,alpha=0.4) + geom_path(data=datatE, aes(datatE[,1],datatE[,2]),size=0.02, color=6,alpha=0.4)


おもしろいのはggplotでデータをdata3と指定している一方、軌道を描く「geom_path」というオプション関数(?)ではそれぞれ渡されているデータが異なるという点だ。

二面角

この記事を閲覧するにはパスワードが必要です
パスワード入力

夜道、イチョウの木の下での結論

この記事を閲覧するにはパスワードが必要です
パスワード入力

Linux コマンド集

データファイルの行数を調べる

$ wc データ名.dat


こうすることで、

451154 31580710 348290302 データ名.dat


などと表示される。意味は、
「このファイルは451154行のファイルで、単語が31580710語含まれていて、サイズは348290302byteです」
ということ。





epsファイルをpngへ変換する


epsで作成した画像をconvertコマンドでpng形式に変換しようとすると、図の中に白い線が入ったり、文字が汚くなったり、画像が粗くなる場合が多い。これは、変換の際に解像度が悪くなるためらしい。以下のようにすると解決された:

$ convert -density 288 -geometry 50% picture.eps picure.png






epsファイルのpng形式への一括変換

$ mogrify -density 100 -format png *.eps







tar.gzファイルの解凍


tar xvzf ファイル名.tar.gz

これが済んだら干潟に行きたい

クリスマス・イブ~クリスマスの夜、2kmある場所「C」から場所「S」まで歩くことになってしまった。「L」共用の自転車の鍵が無くなっていて、用があった場所「S」には徒歩で行かなくてはならなくなったため。一人で行った。たぶん、月曜には鍵の件で「L」内はひと騒ぎだろう。


場所「S」から戻って、聖夜も徹夜で作業である。こんなクリスマスは今年で最後にします!
案外、場所「S」と場所「C」が徒歩でも十分に往復できる距離だということを知った。徒歩30分ぐらい?



やがて朝になり、また、ジョウビタキ目当てにカメラを持って場所「C」内を歩いた。


シロハラがいた。

シロハラ25Dec2010駒場キャンパス_convert_20101225094834



アオジもいた。ただ、夜明け直後だったから画質が悪くなっている。

アオジ0

冬鳥

一昨日、場所「L」で徹夜をした。そして夜が明け、場所「C」内をしばらく歩いた。午前6時半ごろなので、外を歩く者はほとんどいない。しかし、冬鳥が木々の間を飛んでいることに気付いた。秋にはいなかった鳥である。

いったんは場所「C」を出て、買い物をしてからまた場所「L」に戻ろうとテニスコートの横を歩いていたとき、桜の木の下に見慣れない赤い小鳥が見えて、そして、別の木の上にまた留まった。アカハラにしては随分と紅くて、それに、後ろには白い斑点が遠くからでもはっきりと見えた。警戒させないように、かつ迅速に近くへ早歩きし、種類の確認に努めた。

また。赤、というか橙だろうか。紅から橙の間の色だ。遠くでよく見えない。たぶん珍しい鳥だろう、そして、渡り鳥だろうと推測。そう考えている間にも、右手にあった椿の木の中にはメジロがいっぱい。どうやら、鳥の写真を撮るには格好の季節のようだ。2ヶ月前に買った、光学26倍ズームデジタルカメラ「Nikon COOLPIX P100」が遂に初出撃する日がやってきた。

このカメラは、ほとんど生物写真を撮るために買った。被写体がなかったためにしまって置かれていたが、明日にでもこの鳥を撮りに行こうと考えた。朝7時頃、この場所に来ればまた遭遇するはず。場所「L」に戻り、インターネットで鳥の種類を調べた。どうやら「ジョウビタキ」という鳥のようだ。ただ、見えた範囲も狭いので確認できない。

Google画像検索:ジョウビタキ


昨日は決行できなかった。そして、今日になった。今日は天皇誕生日。だが、祝日とは知らず、そして1時間の電車内、普段は満員電車のはずが空いている、これにも自分は異常を感じず、自分は休みの日に、朝の7時から場所「C」に行き、ついに9時頃まで祝日だということに気づかぬままに鳥を追うことになる。

そうとは知らず、たまたま祝日で人がいないことも好都合に、ジョウビタキらしき鳥を撮影しようと場所「C」に到着した。

、と、門を入る前に、フェンスにあの赤い鳥が留まっている!!急いでカメラを鞄から取り出し、構えた。一昨日と違い、今回は距離が近い。頭が灰色であることから、ジョウビタキのオスであることを確認した。是非とも写真に撮っておきたい、が、そもそもこのカメラには慣れていないために操作にも手間取る。そうしているうち、ジョウビタキは上に飛んでいって、レンガづくりのマンションの上に留まった後、どこかへ行ってしまった。撮影失敗。


他の冬鳥でも撮ろうと思い、場所「C」内、および近くの自然公園を散策した。カメラ操作に慣れるため、近くにいたキジバトを撮ってみた。シャッタースピードはマニュアルで操作した方がいいらしい。前はオリンパスC770という機種だったが、ほとんど自動でやってくれた。今回はそうでもない。倍率が高いと、生物写真では写真の質の決定的な要因となる精細さがブレにより大きく減じられる。手ブレ補正がついているとはいえ、生物写真用にはまるで追いつかない。羽毛の一枚一枚が精微に写し出されていなければならない。マニュアルで補正すると満足のいく写真が撮れるようになった。おそらく、三脚がそのうち欲しくなることだろう。プロミレベルの望遠レンズも、金が貯まったら買ってしまうかもしれない。


キジバト





自然公園で、朽ち木の上にツグミが留まっていた。


ツグミ



ジョウビタキは早めに決着をつけたいところだ。いなくなってしまうかもしれないし。
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ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。