読者です 読者をやめる 読者になる 読者になる

GeoJackass

ちゃらんぽらんの絶対領域は、是を頑なに堅持

R graphics

超高機能統計処理ソフトRを用いて、お絵描をします。

まずは、適当なデータセットを用意します。

x <- c(10, 7, 9, 3, 10, 18, 20, 6, 12, 16)
y <- c(20, 10, 10, 36, 17, 14, 30, 12, 7, 16)

とりあえず描画してみましょう。

plot(x)
plot(y)
plot(x, y)
  • グラフオプションが大量に用意されています。type=""で引き渡します。
plot(x, y, type="o")
plot(x, y, type="l")
plot(x, y, type="h")
plot(x, y, type="s")
plot(x, y, type="S")
plot(x, y, type="n")
plot(x, y, type="p")
  • 通常のグラフであれば "l", "h" ぐらいでいいのではないでしょうか。"o"は時系列解析の時に時点数が少ない場合は使った記憶があります。

  • Rは何も指定しない場合は、現在のアクティブな作図ウィンドウに対して、グラフッィクスを上書きしてしまいます。そこで、グラフを書くときに、新たなウィンドウを作成する必要があります。windowsの場合は

win.graph()

次に、現在どの作図ウィンドウを利用しているのかを調べたい場合に使うのが

dev.list()

n番目の作図ウィンドウの消去を行うのが

dev.off(n)

になります。

グラフィックスそのものを終了する場合には

graphics.off()

としてください。

作図したグラフを、後で呼び出したい場合には

sample1 <- recordPlot() #この段階でメモリ上に記録される。
sample1 #これで呼び出しの実行を行う。

作図したグラフをpdfで書き出したい場合には

dev.copy(pdf, file="sample.pdf")

ITOKAWAデータを使用して、意味のあるお絵かきをします。

Hayabusa Project Science Data Archive : JAXA

上記リンクから、Gaskell 形状モデル 25143 Itokawa Plate Model, 49,152 facets Textを選択してダウンロードしてください。 ダウンロードしたデータは、.gzの形式になっていますので、適当なソフトウェアで展開してください。 展開したデータを、改行コードがLFを読めるエディタで読んでみてください。

f:id:GeoJackass:20140404001711p:plain

左から順に1行目にindex、2,3,4行目をそれぞれx, y, zとして読みだせばいいことが分かりますね。

次に、作業用ディレクトリを展開した先に変更してください。 パッケージタブから、あらかじめrglを読めるように準備してください。

itokawa <- read.table("itokawa_f0049152.txt", skip=1, nrow=25350)

x <- itokawa[,2]
y <- itokawa[,3]
z <- itokawa[,4]
//package rglを呼び出す
library(rgl)
//lim指定を行う
plot3d(x = x, y = y, z = z, xlim=c(-0.4, 0.4), ylim=c(-0.4, 0.4), zlim=c(-0.4, 0.4))

f:id:GeoJackass:20140404002536p:plain

出来ました。
ちなみにrglで作図したitokawaは3dモデルで描画していますので、マウスでグルグル回して、好きなアングルを選択して下さい。

C-SODA/ISAS/JAXAのデータとPDS/NASAのデータをビジュアル的に比較する

Asteroid Data Sets これデータ激重たい(69Mb)ので、ダウンロードに30分くらいかかります。

先ほど描画した時に使用した、itokawaのデータはC-SODA/ISAS/JAXAが管理しているデータなのですが、PDS/NASAが管理しているデータもあります。 PDS/NASAのファイルを展開するとHAY_A_AMICA_5_ITOKAWASHAPE_V1_0が出てきます。 その中に\data\vertexがあるので、そこを作業用ディレクトリに変更しましょう。

次に、2つのデータを重ね書きして、同じかどうか検証しましょう。ちなみに、結果は超微妙に異なります。

  • 重ね書きの際に使うコマンドは par(new=T)
  • 塗りつぶしの色を指定する場合は、 col="指定するカラー"
#JAXAの管理するitokawaデータを読み込ませる
itokawa <- read.table("itokawa_f0049152.txt", skip=1, nrow=25350, col="RED")
x <- itokawa[,2]
y <- itokawa[,3]

#NASAの管理するitokawaデータを読み込ませる
itokawa_pds <- read.csv("ver64p.tab", skip=1, nrow=25350)
x1 <- itokawa[,2]
y1 <- itokawa[,3]

plot(x, y)
par(new=T)
plot(x1, y1)

f:id:GeoJackass:20140404012806p:plain

広告を非表示にする