r - Plotting error ellipses in R

Question

I'm trying to plot a U-Pb isochron, which is basically an xy plot, with an error ellipses around each datapoint to indicate the precision of that data point. Here's an example of what I'm trying to achieve, the only difference being that I have a lot more data points:

I have symmetrical + and - errors for each data point in x and y space that I would like to use to dictate the shape of the error ellipse.

Below is the code that I have written so far. This plots up my data, but instead of error ellipses I currently have errorbars.

# 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)

How would I go about replacing the errorbars with error ellipses?

Here's a Google docs link to my data which is in .txt format.

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

Columns are UPb, UPbErrP (error on UPb at 1 sigma %), UPbErrAbs (absolute error on UPb), and then the same again for the PbPb data.

If you need me to clarify anything just let me know and I'll do my best

Answer 1

A few month ago I wrote a little function to draw ellipses to answer someone else's question. By simplifying it a little we can achieve something useful for your issue I think.

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"))