# Description
# 
# Plots soil texture data on an equilateral triangle. Provides the option to
# draw the USDA soil texture classification as a backdrop.
# 
# Usage
# 
# plot.soiltexture <- function(soiltexture, main=NULL, pch=1, col="black",
#   soil.names=TRUE, soil.lines=TRUE, show.points=TRUE,
#   show.clabels=FALSE, show.grid=FALSE, show.legend=FALSE, 
#   col.names="grey", col.lines="grey", bg.pch="white",
#   col.grid="grey", lty.grid=3)
# 
# Arguments
#
# soiltexture   a matrix with three columns (see details)
# main          title of the plot; defaults to NULL
# pch           variable or vector containing point symbols; defaults to 1
# col           variable or vector containing point colour; defaults to black
# soil.names    a logical value indicating whether the soil texture class names are printed; defaults to TRUE
# soil.lines    a logical value indicating whether the soil texture class division lines are plotted; defaults to TRUE
# show.points   a logical value indicating whether the soil sample points are plotted; defaults to TRUE
# show.clabels  a logical value indicating whether the corner labels for 100% sand, silt and clay are printed; defaults to FALSE
# show.grid     a logical value indicating whether the fixed grid is plotted; defaults to FALSE
# show.legend   a logical value indicating whether the legendis plotted; defaults to FALSE
# col.names     colour of the soil texture class names; defaults to grey
# col.lines     colour of the soil texture class division lines; defaults to grey
# bg.pch        colour of the points symbols when pch 21:25 are used; defaults to white
# col.grid      colour of the fixed grid; defaults to grey
# lty.grid      line style for the fixed grid; defaults to 3
#
# Details
# 
# The object soiltexture must be a matrix with at least three columns and one row.
# Columns should contain data for Sand, Silt, and Clay, in that order! Sand, Silt, 
# and Clay should be expressed in proportions between 0 and 1.
#
# The soil sample points' coordinates are calculated using simple trigonometry.
# Thus, the coordinates of a point P(Sand,Silt,Clay), where Sand + Silt + Clay = 1,
# are: P(1-Sand+(Sand-(1-Silt))*0.5, Clay*sin(pi/3))
#
# Author(s)
# 
# Jim Lemon
# Sander Oom
# 
# References
# 
# Soil Survey Division Staff. 1993. Soil survey manual. Soil Conservation Service.
# U.S. Department of Agriculture Handbook 18.
#
# Examples
# 
# # some triangular data
# library(MASS)
# tmp<-(Skye/100)
# # colnames choosen to be consistent with MASS-fig4.4
# colnames(tmp) <- c("Clay","Sand","Silt")
# soiltexture <- cbind(tmp$Sand,tmp$Silt,tmp$Clay)
# # the USDA backdrop in black
# plot.soiltexture(NULL, show.points=FALSE, col.names="black", col.lines="black")
# # the USDA backdrop and a fixed grid in grey
# plot.soiltexture(NULL, show.points=FALSE, show.grid=TRUE)
# # soil sample points with backdrop in grey
# plot.soiltexture(soiltexture)

plot.soiltexture <- function(soiltexture, main="Soil Texture Plot", pch=1, col="black",
  soil.names=TRUE, soil.lines=TRUE,
  show.points=TRUE, show.clabels=FALSE, show.grid=FALSE, show.legend=FALSE, 
  col.names="grey", col.lines="grey", bg.pch="white",
  col.grid="grey", lty.grid=3) {
  
  if(show.points) {
    ## error checking
    if(missing(soiltexture))
      stop("Usage: plot.soiltexture(soiltexture, pch=NULL, col=NULL, soil.names=FALSE, show.grid=FALSE)")
    if(!is.matrix(soiltexture))
      stop("Object soiltexture must be a matrix with at least three columns (for Sand, Silt, and Clay) and one row.")
    if(any(soiltexture > 1) || any(soiltexture < 0))
      stop("All soil texture proportions must be between zero and one.")
    if(any(abs(rowSums(soiltexture)-1) > 0.01))
      warning("At least one set of soil texture proportions does not equal one.")
  }
  
  oldpar<-par(no.readonly=TRUE)
  sin60<-sin(pi/3)
  par(xpd=TRUE) 

  # create empty canvas to draw
  plot(0:1,type="n",axes=FALSE,xlim=c(0,1.1),ylim=c(0,1),
    main=main,xlab="",ylab="")
  
  if(soil.lines) {
    # boundary of clay with extensions
    x1<-c(0.275,0.352,0.6)
    x2<-c(0.415,0.8,0.7)
    y1<-c(0.55*sin60,0.4*sin60,0.4*sin60)
    y2<-c(0.275*sin60,0.4*sin60,0.6*sin60)
    segments(x1,y1,x2,y2, col=col.lines)
    # lower bound and divider of clay loam and silty clay loam
    x1<-c(0.415,0.661)
    x2<-c(0.863,0.6)
    y1<-c(0.275*sin60,0.275*sin60)
    y2<-c(0.275*sin60,0.40*sin60)
    segments(x1,y1,x2,y2, col=col.lines)
    # upper/lower bounds of sandy clay loam and divider with loam
    x1<-c(0.1755,0.1,0.3775)
    x2<-c(0.375,0.3775,0.415)
    y1<-c(0.35*sin60,0.2*sin60,0.2*sin60)
    y2<-c(0.35*sin60,0.2*sin60,0.275*sin60)
    segments(x1,y1,x2,y2, col=col.lines)
    # dividers sand, loamy sand and sandy loam
    x1<-c(0.05,0.075)
    x2<-c(0.15,0.3)
    y1<-c(0.1*sin60,0.15*sin60)
    y2<-c(0,0)
    # lower bounds of loam and upper bounds of silt
    segments(x1,y1,x2,y2, col=col.lines)
    x1<-c(0.377,0.44,0.5,0.8,0.86)
    x2<-c(0.44,0.537,0.638,0.86,0.94)
    y1<-c(0.2*sin60,0.077*sin60,0,0,0.12*sin60)
    y2<-c(0.077*sin60,0.077*sin60,0.275*sin60,0.12*sin60,0.12*sin60)
    segments(x1,y1,x2,y2, col=col.lines)
  }
    
  if(soil.names) {
    # soil texture labels
    par(cex=0.8)
    text(0.5,0.57,"Clay", col=col.names)
    text(0.7,0.48*sin60,"Silty", col=col.names)
    text(0.7,0.44*sin60,"clay", col=col.names)
    text(0.72,0.35*sin60,"Silty clay", col=col.names)
    text(0.73,0.31*sin60,"loam", col=col.names)
    text(0.5,0.34*sin60,"Clay loam", col=col.names)
    text(0.28,0.43*sin60,"Sandy", col=col.names)
    text(0.27,0.39*sin60,"clay", col=col.names)
    text(0.27,0.3*sin60,"Sandy clay", col=col.names)
    text(0.27,0.26*sin60,"loam", col=col.names)
    text(0.25,0.13*sin60,"Sandy loam", col=col.names)
    text(0.13,0.075*sin60,"Loamy", col=col.names)
    text(0.17,0.035*sin60,"sand", col=col.names)
    text(0.065,0.035*sin60,"Sand", col=col.names)
    text(0.49,0.18*sin60,"Loam", col=col.names)
    text(0.72,0.14*sin60,"Silt loam", col=col.names)
    text(0.9,0.06*sin60,"Silt", col=col.names)
    par(cex=1)
  }
  
  # draw corner labels
  if(show.clabels) {
    text(0.5,0.9,"100% clay")
    text(-0.1,0,"100% sand")
    text(1.1,0,"100% loam")
  }
  
  ## the axis labels
  text(0.09,0.43,"% Clay")
  text(0.90,0.43,"% Silt")
  text(0.5,-0.1,"% Sand")
  
  # bottom internal ticks and labels
  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))))
  # left internal ticks and labels
  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 and labels
  rx1<-rev(lx1+0.5-0.01)
  rx2<-rx1+0.01
  ry1<-ly1-0.02*sin60
  ry2<-ly2
  segments(rx1,ry1,rx2,ry2)
  text(rx2+0.03,ry1+0.025,as.character(rev(seq(10,90,by=10))))

  # plot grid
  if(show.grid) {
    segments(bx2,by2,lx1,ly1, lty=lty.grid, col=col.grid)
    segments(lx2,ly2,rx2,ry2, lty=lty.grid, col=col.grid)
    segments(rev(rx1),rev(ry1),bx1,by1, lty=lty.grid, col=col.grid)
  }

  # draw triangle
  x1<-c(0,0,0.5)
  x2<-c(1,0.5,1)
  y1<-c(0,0,sin60)
  y2<-c(0,sin60,0)
  segments(x1,y1,x2,y2)
  
  # plot soil texture points
  if(show.points) {
    if(is.null(pch)) pch<-1:nrow(soiltexture)
    if(is.null(col)) col<-2:(nrow(soiltexture)+1)
    points(1-soiltexture[,1]+(soiltexture[,1]-(1-soiltexture[,2]))*0.5,
      soiltexture[,3]*sin60,pch=pch,col=col, bg=bg.pch)    
  }
        
  # legend
  if(show.legend) {
    samplenames<-rownames(soiltexture)
    legend(0,0.8+0.05*length(samplenames),legend=samplenames,pch=pch,col=col)
  }

  par(oldpar)
}
