5

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:

example plot

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

Shawn Mehan
  • 4,513
  • 9
  • 31
  • 51
cjms85
  • 175
  • 4
  • 2
    Hi Chris, welcome to SO. The format in which you have posted your question requires time and effort on others behalf just to *get your data in*. Instead, focus on a small reproducible example that recreates your problem. Pasting the output of `dput( head( mydata ) )` as a code block will help. See [this post](http://stackoverflow.com/q/5963269/1478381) for how to make a great reproducible example... – Simon O'Hanlon Mar 14 '13 at 14:03
  • 1
    Hi Simon, thanks for the heads up about etiquette. As you can tell I'm completely new to R and stackoverflow, so I tried to do my best to explain myself and give examples and a data set. I'll remember the dput command and the link you posted for the future though! – cjms85 Mar 14 '13 at 14:31
  • No problem. Look how quick it was to be answered - a properly formatted question helps no end. Welcome to the Joy of SO! – Simon O'Hanlon Mar 14 '13 at 14:33

1 Answers1

7

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

enter image description here

Community
  • 1
  • 1
plannapus
  • 18,529
  • 4
  • 72
  • 94
  • Hi! Your code looks great! Do you know if there is a fancy way (like with ggplot2 or Tableu styles) to colour the areas inside the ellipses? – edwineveningfall Aug 27 '15 at 11:39
  • I think you can do this simply by replacing `lines` by `polygon` in the definition of function `ellipse`. – plannapus Aug 27 '15 at 11:45
  • I found a nice answer here: http://stackoverflow.com/questions/19408977/how-to-adjust-the-point-size-to-the-scale-of-the-plot-in-ggplot2 – edwineveningfall Sep 01 '15 at 11:27