0

I would like to create a graph in R looking like this: enter image description here

It is a so-called triaxial ratio diagram to display ratios of plant nutrient contents. It needs a log-scale from 0.01 to 100, the axes cross at 1. I have found two scripts on this page, however, they had another purpose and don't really fit my needs. Here is one: 

get.coords <- function(a, d, x0, y0) {
a <- ifelse(a <= 90,90 - a, 450 - a)
data.frame(x = x0 + d * cos(a / 180 * pi),
y = y0+ d * sin(a / 180 * pi))
}

rotatedAxis <- function(x0, y0, a, d, symmetrical=FALSE, tickdist,  
ticklen,...) {
if(isTRUE(symmetrical)) {
axends <- get.coords(c(a, a + 180), d, x0, y0)
tick.d <- c(seq(0, d, tickdist), seq(-tickdist, -d, -tickdist))
} else {
axends <- rbind(get.coords(a, d, x0, y0), c(x0, y0))
tick.d <- seq(0, d, tickdist)
}
invisible(lapply(apply(get.coords(a, d=tick.d, x0, y0), 1, function(x) {
get.coords(a + 90, c(-ticklen, ticklen), x[1], x[2])
}), function(x) lines(x$x, x$y, ...)))
lines(axends$x, axends$y, ...)
}

plot.new()
plot.window(xlim=c(-70, 70), ylim=c(-70, 70), asp=1)
# Plot the rotated axes -original
rotatedAxis(0, 0, a=60, d=60,symmetrical=TRUE, tickdist=10, ticklen=1)
rotatedAxis(0, 0, a=120, d=60, symmetrical=TRUE, tickdist=10, ticklen=1)
rotatedAxis(0, 0, a=180, d=60, symmetrical=TRUE, tickdist=10, ticklen=1)

# Add text labels to circumference -original
text.coords <- get.coords(seq(0, 300, 60), 65, 0, 0)
text(text.coords$x, text.coords$y, c('I', 'A', 'S', 'E', 'C', 'R'))

points(0, 0, pch=21, bg=1, col=0, lwd=2, cex=2)

This code creates the rotated axis, but not the log-scales and plotting points to the diagram is only possible via xy-coordinates, not xyz.

It would be great if someone could help me, the original paper doesn't offer description of the method and I haven't found anything helpful on the internet. Many Thanks!

bVa
  • 3,839
  • 1
  • 13
  • 22
J.Doe
  • 1
  • 2

1 Answers1

0

This code should give you a plot close to what you desire. I used it with some example data to create the plot given below. The code uses base graphics, so it should be simple to modify. The functions create a plot with log axes from 0.01 to 100.

Because the plotted values are ratios, the only two of them are independent. The third ratio can be calculated from the other two, so only specifying two ratios is required.

I didn't find much on the underlying calculations for the plot either. This code assumes that the ratios for each point are read by tracing the point perpendicular to the axes. For example, the point at 10 on Ca/Mg will not be about 0.3 for each, not 1.

enter image description here

makeAxes <- function(a3, # angles of axes, assumed at 120-degree separation, value indicates positive side from 3-o'clock
                     d, # length of half-axis in log10-units, 2 => 0.01-100
                     axLab, # labels for axes
                     cen=c(0,0), # center of axis
                     ticks=T, # should tick marks be included
                     tickLocs=c((1:9)*10^-2,(1:9)*10^-1,(1:9),(1:9)*10^1,100), # location of tick marks
                     tickLen=0.05, # dimension of tick marks
                     endText=T # are the values at the end of axes plotted
){

  # making an empty plot
  plot(NA,NA
       ,ylim=c(-d,d)*1.5 + cen[1]
       ,xlim=c(-d,d)*1.5 + cen[2]
       ,bty='n'
       ,axes = F
       ,ylab=''
       ,xlab='')

  # adding axis lines and notations
  for(i in 1:3){

    # coordinates of x and y ends of axes
    # col1 = x values
    matTemp <- matrix(0,2,2)

    matTemp[,1] <- c(cen[1] + d*cos(a3[i]*pi/180),cen[1] - d*cos(a3[i]*pi/180))
    matTemp[,2] <- c(cen[2] + d*sin(a3[i]*pi/180),cen[2] - d*sin(a3[i]*pi/180))

    # plotting lines
    lines(matTemp[,1],matTemp[,2])

    # adding text for axes
    text(x=matTemp[1,1]*1.4,y=matTemp[1,2]*1.4,
         labels=axLab[i])

    if(endText){
      text(x=matTemp[,1]*1.1,
           y=matTemp[,2]*1.1,
           labels=c(10^d,10^-d),
           cex=0.8
          )
      text(cen[1],cen[2],
           lab=1,
           pos=4,
           cex=0.8
           )
    }

    # adding tick marks
    if(ticks){

      # values are the point on the axis plus a correction to make it perpendicular
      xVals1 <- cen[1] + log10(tickLocs)*cos(a3[i]*pi/180) - tickLen*sin(a3[i]*pi/180)
      yVals1 <- cen[2] + log10(tickLocs)*sin(a3[i]*pi/180) - tickLen*cos(a3[i]*pi/180)

      xVals2 <- cen[1] + log10(tickLocs)*cos(a3[i]*pi/180) + tickLen*sin(a3[i]*pi/180)
      yVals2 <- cen[2] + log10(tickLocs)*sin(a3[i]*pi/180) + tickLen*cos(a3[i]*pi/180)

      # plotting lines as tick marks
      for(j in 1:length(xVals1)){
        lines(c(xVals1[j],xVals2[j]),c(yVals2[j],yVals1[j]))
      }
    }
  }
}

addPoints <- function(df, # data frame of values, must have axis columns
                     cen=c(0,0), # center of axis
                     a3 # angles of axes, assumed at 120-degree separation, value indicates positive side from 3-o'clock
){

  # temporary data frame for calculations
  dftemp <- df

  # calculation of log values
  dftemp$axis1log <- log10(dftemp$axis1)
  dftemp$axis2log <- log10(dftemp$axis2)
  dftemp$axis3log <- log10(dftemp$axis3)

  # calculation of values along axis 1
  dftemp$axis1log_xvals <- cen[1] + dftemp$axis1log*cos(a3[1]*pi/180)
  dftemp$axis1log_yvals <- cen[2] + dftemp$axis1log*sin(a3[1]*pi/180)

  # calculation of values along axis 2
  dftemp$axis2log_xvals <- cen[1] + dftemp$axis2log*cos(a3[2]*pi/180)
  dftemp$axis2log_yvals <- cen[2] + dftemp$axis2log*sin(a3[2]*pi/180)

  # slope calcuation
  dftemp$axis1_perp <- sin((a3[1]+90)*pi/180)/cos((a3[1]+90)*pi/180)
  dftemp$axis2_perp <- sin((a3[2]+90)*pi/180)/cos((a3[2]+90)*pi/180)

  # intersection points
  # y1 = ax1 + c
  # y2 = bx2 + d
  a <- dftemp$axis1_perp
  b <- dftemp$axis2_perp

  c <- dftemp$axis1log_yvals - dftemp$axis1_perp*dftemp$axis1log_xvals
  d <- dftemp$axis2log_yvals - dftemp$axis2_perp*dftemp$axis2log_xvals

  dftemp$intersectx <- (d-c)/(a-b)
  dftemp$intersecty <- a*(d-c)/(a-b)+c

  points(dftemp$intersectx
         ,dftemp$intersecty
         ,pch=19 # plotting symbol
         ,col='gray48' # colors of points
         ,cex=1.2 # multiplier for plotting symbols
         ) 
}

# example data
labs <- c("a","b","c","d","e","f","g")
k <- c(150,  100, 3000, 3500, 10,  10,  30)
ca <- c(500, 120, 350,  100,  10,  10,  30)
mg <- c(50,  450, 220,  350,  10,  100, 60)

# conversion into data frame
dataPlot <- as.data.frame(list('labs'=labs,'k'=k,'ca'=ca,'mg'=mg)

# calcuations with data frame
dataPlot$axis1 <- dataPlot$ca/dataPlot$mg
dataPlot$axis2 <- dataPlot$k/dataPlot$ca
dataPlot$axis3 <- dataPlot$mg/dataPlot$k

# making axes
makeAxes(a3=c(90,-30,210),
         d=2,
         axLab=c("Ca/Mg", "K/Ca",  "Mg/K" ))

# adding points
addPoints(df=dataPlot,
          a3=c(90,-30,210))
Calvin
  • 1,309
  • 1
  • 14
  • 25