Submitted by dylan on Thu, 2006-10-12 20:40.
See attached file 'pipette_vs_granulometer.csv_.txt' at the bottom of this page.
Load Required Packages and Input Data
#the package 'plotrix' can be installed with:
# install.packages('plotrix', repos='http://cran.r-project.org', dependencies=TRUE)
# note 1: you can accomplish this through the R-GUI in windows / mac os
# note 2: on UNIX-like systems you will need to start R as superuser to
# install packages
#load required package
require(plotrix)
#read in text data: this is a CSV file: sep=","
x <- read.table('pipette_vs_granulometer.csv_.txt', header=T, sep=",")
#subset the data to include only sand, silt, clay columns
p <- data.frame(sand=x$p_sand, silt=x$p_silt, clay=x$p_clay)
g <- data.frame(sand=x$g_sand, silt=x$g_silt, clay=x$g_clay)
Sample 2D Plot: comparison between pipette and laser granulometer sand, silt, and clay values
Simple 2D plot of corelation between pipette and granulometer: for sand, silt, clay
#plot with custom settings:
# annotated axis, 0-100% range, plot symbols scaled by 0.8
plot(p$clay ~ g$clay, ylab="Pct. Clay: Pipette Method", xlab="Pct. Clay: Granulometer Method", main="Simple Plot", xlim=c(0,100), ylim=c(0,100), cex=0.6, pch=16, col=1)
#add a 1:1 corrospondance line:
abline(0,1, col="gray", lty=2)
#add silt fraction to the SAME plot:
points(p$silt ~ g$silt, cex=0.6, pch=16, col=2)
# add sand fraction to the SAME plot:
points(p$sand ~ g$sand, cex=0.6, pch=16, col=3)
#add a legend:
legend(x=2.7, y=94, legend=c('Clay','Silt','Sand'), col=1:3, pch=16, cex=0.6)
#add locally smoothing estimator lines: lowess(x,y)
lines(lowess(g$clay, p$clay), col=1, lwd=2)
lines(lowess(g$silt, p$silt), col=2, lwd=2)
lines(lowess(g$sand, p$sand), col=3, lwd=2)
Soil Texture Triangle 2: visualization of differences between the two methods
Compare the two datasets on the textural triangle
#plot soil textures on triangle
p.tri <- soil.texture(p,show.lines=T, show.names=T, col.symbol='red', pch=16, cex=0.7, main="Soil
Texture Triangle")
g.tri <- triax.points(g, col.symbol='blue', pch=16, cex=0.7)
#create dataframe of segment midpoints
mid <- data.frame(x=(p.tri$x + g.tri$x) / 2, y=(p.tri$y + g.tri$y) / 2 )
#plot a lowess function along the midpoints, ordered by x-coordinate
#what does this really mean?
lines( lowess(mid[order(mid$x), ]) , lwd=2, lty=2, col='green')
# plot overall shift: average p.x,p.y ---> avg g.x,g.y
arrows(mean(p.tri$x), mean(p.tri$y), mean(g.tri$x), mean(g.tri$y), len=.1, lwd=2)
#add a legend:
legend(.68,.79, legend=c("pipette","granulometer"), pch=16, cex=0.7, col=c('red','blue'))
# optionally: add clutter
#plot segments connecting related measurements
segments(p.tri$x, p.tri$y, g.tri$x, g.tri$y, col='gray', lty=1)
