5

基本的にxyプロットであるU-Pbアイソクロンをプロットしようとしています。各データポイントの周りにエラー楕円があり、そのデータポイントの精度を示しています。これは私が達成しようとしているものの例です。唯一の違いは、より多くのデータポイントがあることです:

プロット例

エラー楕円の形状を決定するために使用したい x および y 空間の各データ ポイントに対称的な + および - エラーがあります。

以下は、これまでに書いたコードです。これは私のデータをプロットしますが、エラー楕円の代わりに現在エラーバーがあります。

# Select data file.
fname<-file.choose()
# Rename imported data file.
upbiso <- read.table(file=fname, header= TRUE)
# Attaches data file as default data source.
attach(upbiso)
# Reads and display column headers in console.
names(upbiso)
# Sets margins around plot.
par(mar=c(5,5,4,2))
# Plots 238U/206Pb vs 207Pb/206Pb with superscript notation for isotopic masses, xlim and ylim sets min and max limit for axes, las rotates y axis labels, cex.lab increases the size of the axis labels.
plot(UPb, PbPb, xlab = expression({}^238*"U/"*{}^206*"Pb"), ylab = expression({}^207*"Pb/"*{}^206*"Pb"), xlim = c(0,2500),ylim = c(0, 1), las=1, cex.lab = 1.5)
# Opens Oceanographic package to run errorbars function and runs 1 sigma percent error bars for x-y co-ordinates.
library(oce)
errorbars(UPb,PbPb,UPbErrP,PbPbErrP, percent=T)

エラーバーをエラー楕円に置き換えるにはどうすればよいですか?

これは、.txt 形式のデータへの Google ドキュメント リンクです。

https://docs.google.com/file/d/0B75P9iT4wwlTNUpQand2WVRWV2s/edit?usp=sharing

列はUPb、UPbErrP(1 シグマ % での UPb の誤差)、UPbErrAbs(UPb の絶対誤差)、そしてPbPbデータについても同じです。

何か明確にする必要がある場合は、お知らせください。最善を尽くします

4

1 に答える 1

7

数か月前に、楕円を描画して他の人の質問に答える小さな関数を書きました。少し単純化することで、あなたの問題に役立つ何かを達成できると思います。

ellipse <- function(a,b,xc,yc,...){
    # a is the length of the axis parallel to the x-axis
    # b is the length of the axis parallel to the y-axis
    # xc and yc are the coordinates of the center of the ellipse
    # ... are any arguments that can be passed to function lines
    t <- seq(0, 2*pi, by=pi/100)
    xt <- xc + a*cos(t)
    yt <- yc + b*sin(t)
    lines(xt,yt,...)
    }

plot(UPb, PbPb, pch=19,
     xlab = expression({}^238*"U/"*{}^206*"Pb"), ylab = expression({}^207*"Pb/"*{}^206*"Pb"), 
     xlim = c(0,2500),ylim = c(0, 1), las=1, cex.lab = 1.5)

apply(upbiso, 1, 
      function(x)ellipse(a=x[2]*x[1]/100, b=x[5]*x[4]/100, 
                               xc=x[1], yc=x[4], col="red"))

ここに画像の説明を入力

于 2013-03-14T14:21:32.710 に答える