スポンサーサイト

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

Rで生産構造図を書く方法

層別刈取りしたデータを生産構造図として図に出力する際に、統計ソフトRを使って書く方法を備忘録的にここに記す。別にexcelで書いてもいいんだろうけど、カッコ悪いのと相対照度のデータを上乗せするときにうまくいかないのでRで書いてみた。私は同じラボの先輩の勧めで統計や作図は基本的に全部Rでかいているが、Rを覚えるとかなり便利だ。なんたって無料だからね。

準備として、タブ区切りのtxt形式で保存した乾物(data_DW.txt)と相対照度(data_light.txt)を準備しよう。それらのデータセットをRの作業フォルダに放り込む。あとは以下のプログラムの赤文太文字の群落の一番高いところを指定してやって(以下では200cm)、プログラムを動かせば生産構造図の出来上がり。

##################
### initial setting ###
##################
par(oma =c(1,5,1,5)) #デバイス領域内で,作図領域の外に余白をとる.
par(mfrow=c(1,2)) #画面1行2列に分割する。
par(mar =c(5,0,5,0)) #複数の図間の余白を設定
par(mgp=c(1.0,0.1,0)) #軸からラベル、軸から目盛、軸から軸線までのマージンを設定
par(mgp = c(3, 1.2, 0)) #余白の使い方.説明,ラベル,軸の位置を行で指定.
par(family="Helvetica")

####################
###data preparation ###
####################
H <- 200 #群落の一番高いところを指定する。

data_DW <- read.table("data_DW.txt", header=TRUE, row.names=NULL, sep="\t") #DWのCSVファイルを読み込む
data_light <- read.table("data_light.txt", header=TRUE, row.names=NULL, sep="\t")#相対照度のCSVファイルを読み込む

data_a <- data_DW[,2] #葉身の乾物データ
data_b <- data_DW[,3] #茎及び葉鞘の乾物データ
data_c <- data_DW[,4] #枯死部の乾物データ
data_d <- data_DW[,5] #その他(雑草とか)の乾物データ
data_f <- rbind(data_b, data_c, data_d) #非同化組織の乾物データの複合

hight_light<-data_light[,1] #相対照度の高さ
light <- data_light[,2] #相対照度

result <- nls(hight_light ~ a + b*light^3 + c*log(light), start=list( a=0.1, b=0.1, c=0.1 )) #x^3 + log(x)による非線形の近似曲線
summary(result)
par <- result$m$getPars() #resultのmに格納されているa,b,cの係数を呼び出す

##########################
###start plotting (left side) ###
##########################
barplot(data_a, xlab = expression(paste("assimilatory organ[g DW/",m^{2},"]")),horiz=T, space=0, names.arg=hight,xlim=c(max(data_a)+20,0),axisnames=F, col="green") #axisnames=Fでy軸ラベルを削除

par(new=T); #グラフを重ね合わせます
plot(light, hight_light,ylim=c(0,H), xlim=c(1,100),xlab="",ylab="",axes=F, pch=1,ann=F) #相対照度のグラフを描きます

mtext("relative light intensity[%]",side=3,line=3) #2個目のx軸の名前を表示、sideは方向でaxisと同じ、lineは軸と""内の文字の間隔
mtext("hight [cm]",side=2,line=3) #y軸の名前を表示、sideは方向でaxisと同じ、lineは軸と""内の文字の間隔.lasは文字の向きを指定。

axis(3, las=1)#y軸目盛を縦書き
axis(2, las=2)#y軸目盛を縦書き

result #ちゃんと回帰式が乗っているか確認のため
par(new=T); #グラフを重ね合わせます
curve(par[1] + par[2]*x^3 + par[3]*log(x), xlim=c(0,100), ylim=c(0, H), axes=F, ann=F, col=2)#近似曲線の書き込み

box() #各座標軸の付け根を閉じます

###########################
###start plotting (right side) ###
###########################
barplot(data_f, xlab = expression(paste("nonassimilatory organ[g DW/",m^{2},"]")), horiz=T, space=0, names.arg=hight, xlim=c(0,max(data_f[1,]+data_f[2,]+data_f[3,])), axisnames=F) #非同化組織の棒グラフ導入



テスト用に以下のサンプル・データセットを使うといい。
サンプルデータ

サンプルのデータセットを使うと以下の図が出力されるはずだ。
result_production_structure.jpg
スポンサーサイト

コメントの投稿

非公開コメント

プロフィール
岐阜大学応用生物科学部の助教です。2016年までは日本学術振興会特別研究員DC1で、中国雲南省に住んでいました。中国旅行記というのはその頃から書き始めたブログです。研究者プロフィールはこちらのホームページを御覧ください。今までは中国雲南省と京都を行ったり来たりでしたが、2017年4月からは基本的に岐阜にいます。

tianzhonggui

Author:tianzhonggui
FC2ブログへようこそ!

最新記事
最新コメント
最新トラックバック
月別アーカイブ
カテゴリ
ブログカウンター
検索フォーム
RSSリンクの表示
リンク
ブロとも申請フォーム

この人とブロともになる

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