Re: [R] Soil texture triangle in R?

From: Sander Oom <slist_at_oomvanlieshout.net>
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.
http://soil.scijournals.org/cgi/content/full/65/4/1038

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:
http://soils.usda.gov/technical/manual/images/fig3-16_large.jpg or
http://soil.scijournals.org/content/vol65/issue4/images/large/1038f2.jpeg

Enjoy,

Sander.

# 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)
}



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html 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