Re: [R] Soil texture triangle in R?
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.namesa logical value indicating whether the soil texture class names
are printed; defaults to TRUE
# soil.linesa 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.pchcolour 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 empt
Re: [R] Soil texture triangle in R?
Apologies, removed a crucial line by accident! Here the working version
of the function, with added control over the colours used.
Enjoy!
plot.soiltexture <- function(x,pch,col,colgrid,coltext) {
## triangle plots:
## triangle.plot {ade4}
## triplot{klaR}
## ternaryplot {vcd}
require(vcd)
ternaryplot(x,
grid=FALSE,
dimnames.position = "none",
pch=pch, col=col,
scale=1, main=NULL,
prop.size=FALSE
)
oldpar <- par(no.readonly=TRUE)
ticlength <- 0.01
sin60<-sin(pi/3)
## now the bottom internal ticks
x1<-seq(0.1,0.9,by=0.1)
x2<-x1
y1<-rep(0,9)
y2<-rep(ticlength,9)
segments(x1,y1,x2,y2)
text(x1,y1-0.03,as.character(rev(seq(10,90,by=10 #, cex=0.8)
## now the left internal ticks
y1<-x1*sin60
x1<-x1*0.5
x2<-x1+ticlength*sin60
y2<-y1-ticlength*0.5
segments(x1,y1,x2,y2)
text(x1-0.03,y1+0.015,as.character(seq(10,90,by=10)))
## now the right internal ticks
x1<-rev(x1+0.5-ticlength*sin60)
x2<-x1+ticlength*sin60
segments(x1,y2,x2,y1)
text(x2+0.03,y1+0.015,as.character(rev(seq(10,90,by=10
## the labels at the corners
par(xpd=TRUE)
# 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")
# boundary of clay with extensions
x1<-c(0.275,0.355,0.6)
x2<-c(0.415,0.8,0.7)
y1<-c(0.55*sin60,0.4*sin60,0.4*sin60)
y2<-c(0.285*sin60,0.4*sin60,0.6*sin60)
segments(x1,y1,x2,y2, col=colgrid)
text(0.5,0.57,"Clay", col=coltext)
# lower bound of clay loam & silty divider
x1<-c(0.415,0.66)
x2<-c(0.856,0.6)
y1<-c(0.285*sin60,0.285*sin60)
y2<-c(0.285*sin60,0.40*sin60)
segments(x1,y1,x2,y2, col=colgrid)
text(0.7,0.49*sin60,"Silty", col=coltext)
text(0.7,0.44*sin60,"clay", col=coltext)
text(0.72,0.36*sin60,"Silty clay", col=coltext)
text(0.73,0.32*sin60,"loam", col=coltext)
text(0.5,0.35*sin60,"Clay loam", col=coltext)
x1<-c(0.185,0.1,0.37)
x2<-c(0.37,0.37,0.415)
y1<-c(0.37*sin60,0.2*sin60,0.2*sin60)
y2<-c(0.37*sin60,0.2*sin60,0.285*sin60)
segments(x1,y1,x2,y2, col=colgrid)
text(0.28,0.43*sin60,"Sandy", col=coltext)
text(0.27,0.39*sin60,"clay", col=coltext)
text(0.27,0.3*sin60,"Sandy clay", col=coltext)
text(0.27,0.26*sin60,"loam", col=coltext)
# sand corner
x1<-c(0.05,0.075)
x2<-c(0.15,0.3)
y1<-c(0.1*sin60,0.15*sin60)
y2<-c(0,0)
segments(x1,y1,x2,y2, col=colgrid)
text(0.25,0.13*sin60,"Sandy loam", col=coltext)
text(0.14,0.07*sin60,"Loamy", col=coltext)
text(0.18,0.03*sin60,"sand", col=coltext)
text(0.06,0.021,"Sand", col=coltext)
x1<-c(0.37,0.435,0.5,0.8,0.86)
x2<-c(0.435,0.537,0.64,0.86,0.94)
y1<-c(0.2*sin60,0.08*sin60,0,0,0.12*sin60)
y2<-c(0.08*sin60,0.08*sin60,0.285*sin60,0.12*sin60,0.12*sin60)
segments(x1,y1,x2,y2, col=colgrid)
text(0.49,0.18*sin60,"Loam", col=coltext)
text(0.72,0.15*sin60,"Silt loam", col=coltext)
text(0.9,0.06*sin60,"Silt", col=coltext)
par(oldpar)
}
tmp <- array(dim=c(10,3))
tmp[,2] <- abs(rnorm(10)*20)
tmp[,3] <- abs(rnorm(10)*10)
tmp[,1] <- 100-tmp[,2]-tmp[,3]
col <- rep("black",10)
pch <- rep(1, 10)
plot.soiltexture(tmp,pch,col=NULL,colgrid="black", coltext="black")
plot.soiltexture(tmp,pch,col="black",colgrid="grey", coltext="white")
Sander Oom wrote:
Cleaned up the class divisions and created a full function.
Still to do:
- rotate axis labels;
- correct the partial covering of the bottom tick labels;
- rotate ticks in order to simplify viewing the graph.
See:
http://soil.scijournals.org/content/vol65/issue4/images/large/1038f2.jpeg
Wonder whether triangle.plot{ade4} will give more flexibility!?
Anyway, hopefully the result so far is useful for other people.
Cheers,
Sander.
plot.soiltexture <- function(x,pch,col) {
## triangle plots:
## triangle.plot {ade4}
## triplot{klaR}
## ternaryplot {vcd}
require(vcd)
require(Zelig)
ternaryplot(x,
grid=FALSE,
dimnames.position = "none",
pch=pch, col=col,
scale=1, main=NULL,
prop.size=FALSE
)
oldpar <- par(no.readonly=TRUE)
ticlength <- 0.01
## now the bottom internal ticks
x1<-seq(0.1,0.9,by=0.1)
x2<-x1
y1<-rep(0,9)
y2<-rep(ticlength,9)
segments(x1,y1,x2,y2)
text(x1,y1-0.03,as.character(rev(seq(10,90,by=10 #, cex=0.8)
## now the left internal ticks
y1<-x1*sin60
x1<-x1*0.5
x2<-x1+ticlength*sin60
y2<-y1-ticlength*0.5
segments(x1,y1,x2,y2)
text(x1-0.03,y1+0.015,as.character(seq(10,90,by=10)))
## now the right internal ticks
x1<-rev(x1+0.5-ticlength*sin60)
x2<-x1+ticlength*sin60
segments(x1,y2,x2,y1)
text(x2+0.03,y1+0.015,as.character(rev(seq(10,90,by=10
## the labels at the corners
par(xpd=TRUE)
# 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")
# boundary of clay with extensions
x1<-c(0.275,0.355,0.6)
x2<-c(0
Re: [R] Soil texture triangle in R?
Hi Jim,
Your email was classified as spam, so I missed it previously.
Unfortunately you did not send a cc to the R-help list. I send a cc to
the mailing list now, so all code gets archived.
Thanks for the improvements on your function! It seems that drawing the
ternary graph and the points with the generic plot functions works well.
Makes the function less dependent on other packages!
Will merge the two functions into one and post it back to the mailing
list! Then the graph might be ready for the graph gallery!
Thanks,
Sander.
Jim Lemon wrote:
> Sander Oom wrote:
>> Hi Jim,
>>
>> This looks impressive! It gives me the 'background' graph. However,
>> I'm not sure how I can use this function to plot my soil texture
>> values! Can you explain?
>>
>> I would like to be able to plot my soil texture samples in the same
>> graph as the one your function plots.
>>
> Yes, that's just the background figure for which you attached the link.
> I was unsure of how you were going to use this, that is, do you just
> place a symbol for each soil sample? Aha! Just read your latest email
> and I think that's what you want.
>
> Let's say you have one or more soil samples with the proportions of
> clay, sand and silt.
>
> sample1<-c(0.1,0.4,0.5)
> sample2<-c(0.2,0.5,0.3)
> soil.samples<-rbind(sample1,sample2)
> colnames(soil.samples)<-c("clay","sand","silt")
>
> This more complicated function allows you to overlay colored points
> representing soil samples on either an empty triangle, a gridded
> triangle or a triangle with soil type names. It also has an optional
> legend.
>
> par(ask=TRUE)
> soil.triangle(soil.samples)
> soil.triangle(soil.samples,show.grid=TRUE)
> soil.triangle(soil.samples,soil.names=TRUE,legend=TRUE)
> par(ask=FALSE)
>
> Jim
>
>
>
> soil.triangle<-function(soilprop,pch=NULL,col=NULL,soil.names=FALSE,
> show.grid=FALSE,show.legend=FALSE) {
> if(missing(soilprop))
> stop("Usage:
soil.triangle(soilprop,pch=NULL,col=NULL,soil.names=FALSE,show.grid=FALSE)")
> if(!is.matrix(soilprop))
> stop("soilprop must be a matrix with at least three columns and one
row.")
> if(any(soilprop > 1) || any(soilprop < 0))
> stop("All soil proportions must be between zero and one.")
> if(any(abs(rowSums(soilprop)-1) > 0.01))
> warning("At least one set of soil proportions does not equal one.")
> oldpar<-par(no.readonly=TRUE)
> plot(0:1,type="n",axes=FALSE,xlim=c(0,1.1),ylim=c(0,1),
> main="Soil Triangle",xlab="",ylab="")
> # first draw the triangle
> x1<-c(0,0,0.5)
> sin60<-sin(pi/3)
> x2<-c(1,0.5,1)
> y1<-c(0,0,sin60)
> y2<-c(0,sin60,0)
> segments(x1,y1,x2,y2)
> # now the bottom internal ticks
> 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
> # now the left internal ticks
> 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
> rx1<-rev(lx1+0.5-0.01)
> rx2<-rx1+0.01
> ry1<-ly1-0.02*sin60
> ry2<-ly2
> segments(rx1,ry1,rx2,ry2)
> if(show.grid) {
> segments(bx2,by2,lx1,ly1,lty=3)
> segments(lx2,ly2,rx2,ry2,lty=3)
> segments(rev(rx1),rev(ry1),bx1,by1,lty=3)
> }
> text(rx2+0.03,ry1+0.025,as.character(rev(seq(10,90,by=10
> text(0.5,0.9,"100% clay")
> par(xpd=TRUE)
> text(-0.1,0,"100% sand")
> text(1.1,0,"100% loam")
> text(0.07,0.43,"percent clay")
> text(0.93,0.43,"percent silt")
> text(0.5,-0.1,"percent sand")
> if(soil.names) {
> # boundary of clay with extensions
> x1<-c(0.275,0.35,0.6)
> x2<-c(0.4,0.79,0.7)
> y1<-c(0.55*sin60,0.41*sin60,0.41*sin60)
> y2<-c(0.285*sin60,0.41*sin60,0.6*sin60)
> segments(x1,y1,x2,y2)
> # lower bound of clay loam & silty divider
> x1<-c(0.4,0.68)
> x2<-c(0.86,0.6)
> y1<-c(0.285*sin60,0.285*sin60)
> y2<-c(0.285*sin60,0.41*sin60)
> segments(x1,y1,x2,y2)
> x1<-c(0.185,0.1,0.37)
> x2<-c(0.36,0.37,0.4)
> y1<-c(0.37*sin60,0.2*sin60,0.2*sin60)
> y2<-c(0.37*sin60,0.2*sin60,0.285*sin60)
> segments(x1,y1,x2,y2)
> # sand corner
> x1<-c(0.05,0.075)
> x2<-c(0.12,0.3)
> y1<-c(0.1*sin60,0.15*sin60)
> y2<-c(0,0)
> segments(x1,y1,x2,y2)
> x1<-c(0.37,0.42,0.5,0.8,0.86)
> x2<-c(0.42,0.54,0.65,0.86,0.94)
> y1<-c(0.2*sin60,0.08*sin60,0,0,0.12*sin60)
> y2<-c(0.08*sin60,0.08*sin60,0.285*sin60,0.12*sin60,0.12*sin60)
> segments(x1,y1,x2,y2)
> text(0.5,0.57,"Clay")
> text(0.7,0.49*sin60,"Silty")
> text(0.7,0.44*sin60,"clay")
> text(0.73,0.37*sin60,"Silty clay")
> text(0.73,0.33*sin60,"loam")
> text(0.5,0.35*sin60,"Clay loam")
> text(0.27,0.43*sin60,"Sandy")
> text(0.27,0.39*sin60,"clay")
> text(0.27,0.3*sin60,"Sandy clay")
> text(0.27,0.26*sin60,"loam")
> text(0.25,0.13*sin60,"Sandy loam")
> text(0.13,0.075*sin60,"Loamy")
> text(0.15,0.035*sin60,"sand")
> text(0.0
Re: [R] Soil texture triangle in R?
Cleaned up the class divisions and created a full function.
Still to do:
- rotate axis labels;
- correct the partial covering of the bottom tick labels;
- rotate ticks in order to simplify viewing the graph.
See:
http://soil.scijournals.org/content/vol65/issue4/images/large/1038f2.jpeg
Wonder whether triangle.plot{ade4} will give more flexibility!?
Anyway, hopefully the result so far is useful for other people.
Cheers,
Sander.
plot.soiltexture <- function(x,pch,col) {
## triangle plots:
## triangle.plot {ade4}
## triplot{klaR}
## ternaryplot {vcd}
require(vcd)
require(Zelig)
ternaryplot(x,
grid=FALSE,
dimnames.position = "none",
pch=pch, col=col,
scale=1, main=NULL,
prop.size=FALSE
)
oldpar <- par(no.readonly=TRUE)
ticlength <- 0.01
## now the bottom internal ticks
x1<-seq(0.1,0.9,by=0.1)
x2<-x1
y1<-rep(0,9)
y2<-rep(ticlength,9)
segments(x1,y1,x2,y2)
text(x1,y1-0.03,as.character(rev(seq(10,90,by=10 #, cex=0.8)
## now the left internal ticks
y1<-x1*sin60
x1<-x1*0.5
x2<-x1+ticlength*sin60
y2<-y1-ticlength*0.5
segments(x1,y1,x2,y2)
text(x1-0.03,y1+0.015,as.character(seq(10,90,by=10)))
## now the right internal ticks
x1<-rev(x1+0.5-ticlength*sin60)
x2<-x1+ticlength*sin60
segments(x1,y2,x2,y1)
text(x2+0.03,y1+0.015,as.character(rev(seq(10,90,by=10
## the labels at the corners
par(xpd=TRUE)
# 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")
# boundary of clay with extensions
x1<-c(0.275,0.355,0.6)
x2<-c(0.415,0.8,0.7)
y1<-c(0.55*sin60,0.4*sin60,0.4*sin60)
y2<-c(0.285*sin60,0.4*sin60,0.6*sin60)
segments(x1,y1,x2,y2, col="grey")
text(0.5,0.57,"Clay", col="grey")
# lower bound of clay loam & silty divider
x1<-c(0.415,0.66)
x2<-c(0.856,0.6)
y1<-c(0.285*sin60,0.285*sin60)
y2<-c(0.285*sin60,0.40*sin60)
segments(x1,y1,x2,y2, col="grey")
text(0.7,0.49*sin60,"Silty", col="grey")
text(0.7,0.44*sin60,"clay", col="grey")
text(0.72,0.36*sin60,"Silty clay", col="grey")
text(0.73,0.32*sin60,"loam", col="grey")
text(0.5,0.35*sin60,"Clay loam", col="grey")
x1<-c(0.185,0.1,0.37)
x2<-c(0.37,0.37,0.415)
y1<-c(0.37*sin60,0.2*sin60,0.2*sin60)
y2<-c(0.37*sin60,0.2*sin60,0.285*sin60)
segments(x1,y1,x2,y2, col="grey")
text(0.28,0.43*sin60,"Sandy", col="grey")
text(0.27,0.39*sin60,"clay", col="grey")
text(0.27,0.3*sin60,"Sandy clay", col="grey")
text(0.27,0.26*sin60,"loam", col="grey")
# sand corner
x1<-c(0.05,0.075)
x2<-c(0.15,0.3)
y1<-c(0.1*sin60,0.15*sin60)
y2<-c(0,0)
segments(x1,y1,x2,y2, col="grey")
text(0.25,0.13*sin60,"Sandy loam", col="grey")
text(0.14,0.07*sin60,"Loamy", col="grey")
text(0.18,0.03*sin60,"sand", col="grey")
text(0.06,0.021,"Sand", col="grey")
x1<-c(0.37,0.435,0.5,0.8,0.86)
x2<-c(0.435,0.537,0.64,0.86,0.94)
y1<-c(0.2*sin60,0.08*sin60,0,0,0.12*sin60)
y2<-c(0.08*sin60,0.08*sin60,0.285*sin60,0.12*sin60,0.12*sin60)
segments(x1,y1,x2,y2, col="grey")
text(0.49,0.18*sin60,"Loam", col="grey")
text(0.72,0.15*sin60,"Silt loam", col="grey")
text(0.9,0.06*sin60,"Silt", col="grey")
ternarypoints(x, pch = pch, col = col)
par(oldpar)
}
tmp <- array(dim=c(10,3))
tmp[,2] <- abs(rnorm(10)*20)
tmp[,3] <- abs(rnorm(10)*10)
tmp[,1] <- 100-tmp[,2]-tmp[,3]
col <- rep("black",10)
pch <- rep(1, 10)
plot.soiltexture(tmp,pch,col="black")
Sander Oom wrote:
Right,
Got the data points plotted on top of the soil texture background,
thanks to Jim and ternaryplot{vcd}! See code below.
Now there is some fine tuning to do, as it should really look like this
graph:
http://soil.scijournals.org/content/vol65/issue4/images/large/1038f2.jpeg
Things to do:
- rotate axis labels;
- correct small errors in class divisions;
- correct the partial covering of the bottom tick labels;
- rotate ticks in order to simplify viewing the graph.
Any help still appreciated!
Cheers,
Sander.
soil.triangle <- function() {
oldpar <- par(no.readonly=TRUE)
## now the bottom internal ticks
x1<-seq(0.1,0.9,by=0.1)
x2<-x1
y1<-rep(0,9)
y2<-rep(-0.02,9)
segments(x1,y1,x2,y2)
text(x1,y1-0.03,as.character(rev(seq(10,90,by=10 #, cex=0.8)
## now the left internal ticks
y1<-x1*sin60
x1<-x1*0.5
x2<-x1+0.02*sin60
y2<-y1-0.02*0.5
segments(x1,y1,x2,y2)
text(x1-0.03,y1+0.015,as.character(seq(10,90,by=10)))
## now the right internal ticks
x1<-rev(x1+0.5-0.02*sin60)
x2<-x1+0.02*sin60
segments(x1,y2,x2,y1)
text(x2+0.03,y1+0.015,as.character(rev(seq(10,90,by=10
## the labels at the corners
par(xpd=TRUE)
# text(0.5,0.9,"100% clay")
# text(-0.1,0,"100% sand")
# text(1.1,0,"100% loam")
text(0.09,0.43,"% Clay")
text(0.90,0.43,"% Silt")
text(0.5,-0.1,"% Sand")
# boundary of clay with extensions
x1<-c(0.275,0.35,0.6)
x2<-c(0.4,0.79,0.7)
y1<-c(0.55*sin60,0.41*sin60,0.41*sin60)
Re: [R] Soil texture triangle in R?
Right,
Got the data points plotted on top of the soil texture background,
thanks to Jim and ternaryplot{vcd}! See code below.
Now there is some fine tuning to do, as it should really look like this
graph:
http://soil.scijournals.org/content/vol65/issue4/images/large/1038f2.jpeg
Things to do:
- rotate axis labels;
- correct small errors in class divisions;
- correct the partial covering of the bottom tick labels;
- rotate ticks in order to simplify viewing the graph.
Any help still appreciated!
Cheers,
Sander.
soil.triangle <- function() {
oldpar <- par(no.readonly=TRUE)
## now the bottom internal ticks
x1<-seq(0.1,0.9,by=0.1)
x2<-x1
y1<-rep(0,9)
y2<-rep(-0.02,9)
segments(x1,y1,x2,y2)
text(x1,y1-0.03,as.character(rev(seq(10,90,by=10 #, cex=0.8)
## now the left internal ticks
y1<-x1*sin60
x1<-x1*0.5
x2<-x1+0.02*sin60
y2<-y1-0.02*0.5
segments(x1,y1,x2,y2)
text(x1-0.03,y1+0.015,as.character(seq(10,90,by=10)))
## now the right internal ticks
x1<-rev(x1+0.5-0.02*sin60)
x2<-x1+0.02*sin60
segments(x1,y2,x2,y1)
text(x2+0.03,y1+0.015,as.character(rev(seq(10,90,by=10
## the labels at the corners
par(xpd=TRUE)
# text(0.5,0.9,"100% clay")
# text(-0.1,0,"100% sand")
# text(1.1,0,"100% loam")
text(0.09,0.43,"% Clay")
text(0.90,0.43,"% Silt")
text(0.5,-0.1,"% Sand")
# boundary of clay with extensions
x1<-c(0.275,0.35,0.6)
x2<-c(0.4,0.79,0.7)
y1<-c(0.55*sin60,0.41*sin60,0.41*sin60)
y2<-c(0.285*sin60,0.41*sin60,0.6*sin60)
segments(x1,y1,x2,y2, col="grey")
text(0.5,0.57,"Clay", col="grey")
# lower bound of clay loam & silty divider
x1<-c(0.4,0.68)
x2<-c(0.86,0.6)
y1<-c(0.285*sin60,0.285*sin60)
y2<-c(0.285*sin60,0.41*sin60)
segments(x1,y1,x2,y2, col="grey")
text(0.7,0.49*sin60,"Silty", col="grey")
text(0.7,0.44*sin60,"clay", col="grey")
text(0.73,0.37*sin60,"Silty clay", col="grey")
text(0.73,0.33*sin60,"loam", col="grey")
text(0.5,0.35*sin60,"Clay loam", col="grey")
x1<-c(0.185,0.1,0.37)
x2<-c(0.36,0.37,0.4)
y1<-c(0.37*sin60,0.2*sin60,0.2*sin60)
y2<-c(0.37*sin60,0.2*sin60,0.285*sin60)
segments(x1,y1,x2,y2, col="grey")
text(0.27,0.43*sin60,"Sandy", col="grey")
text(0.27,0.39*sin60,"clay", col="grey")
text(0.27,0.3*sin60,"Sandy clay", col="grey")
text(0.27,0.26*sin60,"loam", col="grey")
# sand corner
x1<-c(0.05,0.075)
x2<-c(0.12,0.3)
y1<-c(0.1*sin60,0.15*sin60)
y2<-c(0,0)
segments(x1,y1,x2,y2, col="grey")
text(0.25,0.13*sin60,"Sandy loam", col="grey")
text(0.13,0.075*sin60,"Loamy", col="grey")
text(0.15,0.035*sin60,"sand", col="grey")
text(0.055,0.021,"Sand", col="grey")
x1<-c(0.37,0.42,0.5,0.8,0.86)
x2<-c(0.42,0.54,0.65,0.86,0.94)
y1<-c(0.2*sin60,0.08*sin60,0,0,0.12*sin60)
y2<-c(0.08*sin60,0.08*sin60,0.285*sin60,0.12*sin60,0.12*sin60)
segments(x1,y1,x2,y2, col="grey")
text(0.49,0.18*sin60,"Loam", col="grey")
text(0.72,0.15*sin60,"Silt loam", col="grey")
text(0.9,0.06*sin60,"Silt", col="grey")
par(oldpar)
}
tmp <- array(dim=c(10,3))
tmp[,1] <- abs(rnorm(10)*20)
tmp[,2] <- abs(rnorm(10)*10)
tmp[,3] <- 100-tmp[,1]-tmp[,2]
tmp
library(vcd)
## Mark groups
ternaryplot(tmp,
grid=FALSE,
dimnames.position = "none",
pch=1, col="black",
scale=1, main=NULL,
prop.size=FALSE,
)
soil.triangle()
Sander Oom wrote:
Hi Jim,
This looks impressive! It gives me the 'background' graph. However, I'm
not sure how I can use this function to plot my soil texture values! Can
you explain?
I would like to be able to plot my soil texture samples in the same
graph as the one your function plots.
Maybe I should try to figure out how to replicate your code inside a
ternaryplot{vcd} call.
Cheers,
Sander.
Jim Lemon wrote:
> Sander Oom wrote:
>> Dear R users,
>>
>> has anybody made an attempt to create the soil texture triangle graph
>> in R? For an example see here:
>>
>> http://www.teachingkate.org/images/soiltria.gif
>>
>> I would like to get the lines in black and texture labels in gray to
>> allow for plotting my texture results on top.
>>
>> Any examples or suggestions are very welcome!
>>
> It's not too hard to write a plot function to do this, but I'm not sure
> whether this will be what you want. Anyway, try it out.
>
> Jim
>
>
>
> soil.triangle<-function() {
> oldpar<-par(no.readonly=TRUE)
> plot(0:1,type="n",axes=FALSE,xlim=c(0,1.1),ylim=c(0,1),
> main="Soil Triangle",xlab="",ylab="")
> # first draw the triangle
> x1<-c(0,0,0.5)
> sin60<-sin(pi/3)
> x2<-c(1,0.5,1)
> y1<-c(0,0,sin60)
> y2<-c(0,sin60,0)
> segments(x1,y1,x2,y2)
> # now the bottom internal ticks
> x1<-seq(0.1,0.9,by=0.1)
> x2<-x1
> y1<-rep(0,9)
> y2<-rep(0.02,9)
> segments(x1,y1,x2,y2)
> text(x1,y1-0.03,as.character(rev(seq(10,90,by=10
> # now the left internal ticks
> y1<-x1*sin60
> x1<-x1*0.5
> x2<-x1+0.02*sin60
> y2<-y1-0.02*0.5
> segments(x1,y1,x2
Re: [R] Soil texture triangle in R?
Hi Jari,
I assume this has been superseded by the ternaryplot{vcd} function!?
Thanks,
Sander.
Jari Oksanen wrote:
> Sander,
>
> Just a quick note before I go to the field.
>
> I attach a tri.R file for drawing ternary plots. The base function was
> posted to R News someday. One thing that I added was option to plot
> nothing (type="n") plus (invisible) return of plotting coordinates. This
> means that you can take the coordinates for drawing segments (in
> original values), feed them through this functions, and you get
> translated coordinates to use ordinary segments or lines commands to
> overlay your lines into an existing ternary plot.
>
> We used this in an Applied Vegetation Science paper (Hellström as the
> first author) to overlay arrows onto ternary plots.
>
> cheers, jari oksanen
>
>
>
>
> "tri" <-
> function(a, f, m, symb = 2, grid = F, ...)
> {
> ta <- paste(substitute(a))
> tf <- paste(substitute(f))
> tm <- paste(substitute(m))
>
> tot <- 100/(a + f +m)
> b <- f * tot
> y <- b * .878
> x <- m * tot + b/2
> par(pty = "s")
> oldcol <- par("col")
> plot(x, y, axes = F, xlab = "", ylab = "", xlim = c(-10, 110), ylim
> = c(-10, 110), type = "n", ...)
> points(x,y,pch=symb)
> par(col = oldcol)
> trigrid(grid)
> text(-5, -5, ta)
> text(105, -5, tm)
> text(50, 93, tf)
> par(pty = "m")
> invisible(cbind(x,y))
> }
>
__
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
Re: [R] Soil texture triangle in R?
Hi Jim,
This looks impressive! It gives me the 'background' graph. However, I'm
not sure how I can use this function to plot my soil texture values! Can
you explain?
I would like to be able to plot my soil texture samples in the same
graph as the one your function plots.
Maybe I should try to figure out how to replicate your code inside a
ternaryplot{vcd} call.
Cheers,
Sander.
Jim Lemon wrote:
> Sander Oom wrote:
>> Dear R users,
>>
>> has anybody made an attempt to create the soil texture triangle graph
>> in R? For an example see here:
>>
>> http://www.teachingkate.org/images/soiltria.gif
>>
>> I would like to get the lines in black and texture labels in gray to
>> allow for plotting my texture results on top.
>>
>> Any examples or suggestions are very welcome!
>>
> It's not too hard to write a plot function to do this, but I'm not sure
> whether this will be what you want. Anyway, try it out.
>
> Jim
>
>
>
> soil.triangle<-function() {
> oldpar<-par(no.readonly=TRUE)
> plot(0:1,type="n",axes=FALSE,xlim=c(0,1.1),ylim=c(0,1),
> main="Soil Triangle",xlab="",ylab="")
> # first draw the triangle
> x1<-c(0,0,0.5)
> sin60<-sin(pi/3)
> x2<-c(1,0.5,1)
> y1<-c(0,0,sin60)
> y2<-c(0,sin60,0)
> segments(x1,y1,x2,y2)
> # now the bottom internal ticks
> x1<-seq(0.1,0.9,by=0.1)
> x2<-x1
> y1<-rep(0,9)
> y2<-rep(0.02,9)
> segments(x1,y1,x2,y2)
> text(x1,y1-0.03,as.character(rev(seq(10,90,by=10
> # now the left internal ticks
> y1<-x1*sin60
> x1<-x1*0.5
> x2<-x1+0.02*sin60
> y2<-y1-0.02*0.5
> segments(x1,y1,x2,y2)
> text(x1-0.03,y1+0.015,as.character(seq(10,90,by=10)))
> x1<-rev(x1+0.5-0.02*sin60)
> x2<-x1+0.02*sin60
> segments(x1,y2,x2,y1)
> text(x2+0.03,y1+0.015,as.character(rev(seq(10,90,by=10
> text(0.5,0.9,"100% clay")
> par(xpd=TRUE)
> text(-0.1,0,"100% sand")
> text(1.1,0,"100% loam")
> text(0.07,0.43,"percent clay")
> text(0.93,0.43,"percent silt")
> text(0.5,-0.1,"percent sand")
> # boundary of clay with extensions
> x1<-c(0.275,0.35,0.6)
> x2<-c(0.4,0.79,0.7)
> y1<-c(0.55*sin60,0.41*sin60,0.41*sin60)
> y2<-c(0.285*sin60,0.41*sin60,0.6*sin60)
> segments(x1,y1,x2,y2)
> text(0.5,0.57,"Clay")
> # lower bound of clay loam & silty divider
> x1<-c(0.4,0.68)
> x2<-c(0.86,0.6)
> y1<-c(0.285*sin60,0.285*sin60)
> y2<-c(0.285*sin60,0.41*sin60)
> segments(x1,y1,x2,y2)
> text(0.7,0.49*sin60,"Silty")
> text(0.7,0.44*sin60,"clay")
> text(0.73,0.37*sin60,"Silty clay")
> text(0.73,0.33*sin60,"loam")
> text(0.5,0.35*sin60,"Clay loam")
> x1<-c(0.185,0.1,0.37)
> x2<-c(0.36,0.37,0.4)
> y1<-c(0.37*sin60,0.2*sin60,0.2*sin60)
> y2<-c(0.37*sin60,0.2*sin60,0.285*sin60)
> segments(x1,y1,x2,y2)
> text(0.27,0.43*sin60,"Sandy")
> text(0.27,0.39*sin60,"clay")
> text(0.27,0.3*sin60,"Sandy clay")
> text(0.27,0.26*sin60,"loam")
> # sand corner
> x1<-c(0.05,0.075)
> x2<-c(0.12,0.3)
> y1<-c(0.1*sin60,0.15*sin60)
> y2<-c(0,0)
> segments(x1,y1,x2,y2)
> text(0.25,0.13*sin60,"Sandy loam")
> text(0.13,0.075*sin60,"Loamy")
> text(0.15,0.035*sin60,"sand")
> text(0.055,0.021,"Sand")
> x1<-c(0.37,0.42,0.5,0.8,0.86)
> x2<-c(0.42,0.54,0.65,0.86,0.94)
> y1<-c(0.2*sin60,0.08*sin60,0,0,0.12*sin60)
> y2<-c(0.08*sin60,0.08*sin60,0.285*sin60,0.12*sin60,0.12*sin60)
> segments(x1,y1,x2,y2)
> text(0.49,0.18*sin60,"Loam")
> text(0.72,0.15*sin60,"Silt loam")
> text(0.9,0.06*sin60,"Silt")
> 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
[R] Soil texture triangle in R?
Dear R users, has anybody made an attempt to create the soil texture triangle graph in R? For an example see here: http://www.teachingkate.org/images/soiltria.gif I would like to get the lines in black and texture labels in gray to allow for plotting my texture results on top. Any examples or suggestions are very welcome! Thanks in advance, Sander. -- Dr Sander P. Oom Animal, Plant and Environmental Sciences, University of the Witwatersrand Private Bag 3, Wits 2050, South Africa Tel (work) +27 (0)11 717 64 04 Tel (home) +27 (0)18 297 44 51 Fax +27 (0)18 299 24 64 Email [EMAIL PROTECTED] Web www.oomvanlieshout.net/sander __ 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
