[R] Soil texture triangle in R?

[Sander Oom]

Hi Jim,

Your email was classified as spam, so I missed it previously. 
Unfortunately you did not send a cc to the R-help list. I send a cc to 
the mailing list now, so all code gets archived.

Thanks for the improvements on your function! It seems that drawing the 
ternary graph and the points with the generic plot functions works well. 
Makes the function less dependent on other packages!

Will merge the two functions into one and post it back to the mailing 
list! Then the graph might be ready for the graph gallery!

Thanks,

Sander.


Jim Lemon wrote:
 > Sander Oom wrote:
 >> Hi Jim,
 >>
 >> This looks impressive! It gives me the 'background' graph. However,
 >> I'm not sure how I can use this function to plot my soil texture
 >> values! Can you explain?
 >>
 >> I would like to be able to plot my soil texture samples in the same
 >> graph as the one your function plots.
 >>
 > Yes, that's just the background figure for which you attached the link.
 > I was unsure of how you were going to use this, that is, do you just
 > place a symbol for each soil sample? Aha! Just read your latest email
 > and I think that's what you want.
 >
 > Let's say you have one or more soil samples with the proportions of
 > clay, sand and silt.
 >
 > sample1<-c(0.1,0.4,0.5)
 > sample2<-c(0.2,0.5,0.3)
 > soil.samples<-rbind(sample1,sample2)
 > colnames(soil.samples)<-c("clay","sand","silt")
 >
 > This more complicated function allows you to overlay colored points
 > representing soil samples on either an empty triangle, a gridded
 > triangle or a triangle with soil type names. It also has an optional
 > legend.
 >
 > par(ask=TRUE)
 > soil.triangle(soil.samples)
 > soil.triangle(soil.samples,show.grid=TRUE)
 > soil.triangle(soil.samples,soil.names=TRUE,legend=TRUE)
 > par(ask=FALSE)
 >
 > Jim
 >
 > ------------------------------------------------------------------------
 >
 > soil.triangle<-function(soilprop,pch=NULL,col=NULL,soil.names=FALSE,
 >  show.grid=FALSE,show.legend=FALSE) {
 >  if(missing(soilprop))
 >   stop("Usage: 
soil.triangle(soilprop,pch=NULL,col=NULL,soil.names=FALSE,show.grid=FALSE)")
 >  if(!is.matrix(soilprop))
 >   stop("soilprop must be a matrix with at least three columns and one 
row.")
 >  if(any(soilprop > 1) || any(soilprop < 0))
 >   stop("All soil proportions must be between zero and one.")
 >  if(any(abs(rowSums(soilprop)-1) > 0.01))
 >   warning("At least one set of soil proportions does not equal one.")
 >  oldpar<-par(no.readonly=TRUE)
 >  plot(0:1,type="n",axes=FALSE,xlim=c(0,1.1),ylim=c(0,1),
 >   main="Soil Triangle",xlab="",ylab="")
 >  # first draw the triangle
 >  x1<-c(0,0,0.5)
 >  sin60<-sin(pi/3)
 >  x2<-c(1,0.5,1)
 >  y1<-c(0,0,sin60)
 >  y2<-c(0,sin60,0)
 >  segments(x1,y1,x2,y2)
 >  # now the bottom internal ticks
 >  bx1<-seq(0.1,0.9,by=0.1)
 >  bx2<-bx1-0.01
 >  by1<-rep(0,9)
 >  by2<-rep(0.02*sin60,9)
 >  segments(bx1,by1,bx2,by2)
 >  text(bx1,by1-0.03,as.character(rev(seq(10,90,by=10))))
 >  # now the left internal ticks
 >  ly1<-bx1*sin60
 >  lx1<-bx1*0.5
 >  lx2<-lx1+0.02
 >  ly2<-ly1
 >  segments(lx1,ly1,lx2,ly2)
 >  text(lx1-0.03,ly1,as.character(seq(10,90,by=10)))
 >  # right internal ticks
 >  rx1<-rev(lx1+0.5-0.01)
 >  rx2<-rx1+0.01
 >  ry1<-ly1-0.02*sin60
 >  ry2<-ly2
 >  segments(rx1,ry1,rx2,ry2)
 >  if(show.grid) {
 >   segments(bx2,by2,lx1,ly1,lty=3)
 >   segments(lx2,ly2,rx2,ry2,lty=3)
 >   segments(rev(rx1),rev(ry1),bx1,by1,lty=3)
 >  }
 >  text(rx2+0.03,ry1+0.025,as.character(rev(seq(10,90,by=10))))
 >  text(0.5,0.9,"100% clay")
 >  par(xpd=TRUE)
 >  text(-0.1,0,"100% sand")
 >  text(1.1,0,"100% loam")
 >  text(0.07,0.43,"percent clay")
 >  text(0.93,0.43,"percent silt")
 >  text(0.5,-0.1,"percent sand")
 >  if(soil.names) {
 >   # boundary of clay with extensions
 >   x1<-c(0.275,0.35,0.6)
 >   x2<-c(0.4,0.79,0.7)
 >   y1<-c(0.55*sin60,0.41*sin60,0.41*sin60)
 >   y2<-c(0.285*sin60,0.41*sin60,0.6*sin60)
 >   segments(x1,y1,x2,y2)
 >   # lower bound of clay loam & silty divider
 >   x1<-c(0.4,0.68)
 >   x2<-c(0.86,0.6)
 >   y1<-c(0.285*sin60,0.285*sin60)
 >   y2<-c(0.285*sin60,0.41*sin60)
 >   segments(x1,y1,x2,y2)
 >   x1<-c(0.185,0.1,0.37)
 >   x2<-c(0.36,0.37,0.4)
 >   y1<-c(0.37*sin60,0.2*sin60,0.2*sin60)
 >   y2<-c(0.37*sin60,0.2*sin60,0.285*sin60)
 >   segments(x1,y1,x2,y2)
 >   # sand corner
 >   x1<-c(0.05,0.075)
 >   x2<-c(0.12,0.3)
 >   y1<-c(0.1*sin60,0.15*sin60)
 >   y2<-c(0,0)
 >   segments(x1,y1,x2,y2)
 >   x1<-c(0.37,0.42,0.5,0.8,0.86)
 >   x2<-c(0.42,0.54,0.65,0.86,0.94)
 >   y1<-c(0.2*sin60,0.08*sin60,0,0,0.12*sin60)
 >   y2<-c(0.08*sin60,0.08*sin60,0.285*sin60,0.12*sin60,0.12*sin60)
 >   segments(x1,y1,x2,y2)
 >   text(0.5,0.57,"Clay")
 >   text(0.7,0.49*sin60,"Silty")
 >   text(0.7,0.44*sin60,"clay")
 >   text(0.73,0.37*sin60,"Silty clay")
 >   text(0.73,0.33*sin60,"loam")
 >   text(0.5,0.35*sin60,"Clay loam")
 >   text(0.27,0.43*sin60,"Sandy")
 >   text(0.27,0.39*sin60,"clay")
 >   text(0.27,0.3*sin60,"Sandy clay")
 >   text(0.27,0.26*sin60,"loam")
 >   text(0.25,0.13*sin60,"Sandy loam")
 >   text(0.13,0.075*sin60,"Loamy")
 >   text(0.15,0.035*sin60,"sand")
 >   text(0.055,0.021,"Sand")
 >   text(0.49,0.18*sin60,"Loam")
 >   text(0.72,0.15*sin60,"Silt loam")
 >   text(0.9,0.06*sin60,"Silt")
 >  }
 >  if(is.null(pch)) pch<-1:nrow(soilprop)
 >  if(is.null(col)) col<-2:(nrow(soilprop)+1)
 >  points(1-soilprop[,2]+(soilprop[,2]-(1-soilprop[,3]))*0.5,
 >   soilprop[,1]*sin60,pch=pch,col=col)
 >  if(show.legend) {
 >   samplenames<-rownames(soilprop)
 > 
legend(0,0.8+0.05*length(samplenames),legend=samplenames,pch=pch,col=col)
 >  }
 >  par(oldpar)
 > }