Re: [R] Soil texture triangle in R?

From: Sander Oom <>
Date: Sun 29 May 2005 - 05:51:46 EST

Dear R users,

Please find attached a new plot function, plot.soiltexture, to plot soil texture data on a triangular plot with an optional backdrop of the USDA soil texture classification, written by Jim Lemon and me.

I tried to write the function and documentation confirm the R conventions. However, this is a new experience for me, so any comments and suggestions are welcome!

I tried to find a suitable package for the plot function, but none are obvious.

I have approached Todd Skaggs to ask permission to include sample data from the paper:
Skaggs, T.H., L.M. Arya, P.J. Shouse, and B.P. Mohanty. 2001. Estimating particle-size distribution from limited soil texture data. Soil Sci. Soc. Am. J., 65:1038-1044.

Things still to do:

- rotate axis labels;
- rotate axis tick labels
- provide option to plot ticks inside or outside plot area.

Thus making it look like: or



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


  # 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

    segments(x1,y1,x2,y2, col=col.lines)     # lower bound and divider of clay loam and silty clay loam

    segments(x1,y1,x2,y2, col=col.lines)     # upper/lower bounds of sandy clay loam and divider with loam
    segments(x1,y1,x2,y2, col=col.lines)     # dividers sand, loamy sand and sandy loam

    # lower bounds of loam and upper bounds of silt     segments(x1,y1,x2,y2, col=col.lines)
    segments(x1,y1,x2,y2, col=col.lines)   }     

  if(soil.names) {
    # soil texture labels

    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)

  # 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


  text(bx1,by1-0.03,as.character(rev(seq(10,90,by=10))))   # left internal ticks and labels

  text(lx1-0.03,ly1,as.character(seq(10,90,by=10)))   # right internal ticks and labels


  # 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



  # plot soil texture points

  if(show.points) {
    if(is.null(pch)) pch<-1:nrow(soiltexture)
    if(is.null(col)) col<-2:(nrow(soiltexture)+1)
      soiltexture[,3]*sin60,pch=pch,col=col, bg=bg.pch)    

  # legend
  if(show.legend) {
    legend(0,0.8+0.05*length(samplenames),legend=samplenames,pch=pch,col=col)   }

} mailing list PLEASE do read the posting guide! Received on Sun May 29 06:02:10 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:13 EST