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
