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.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)
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,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 at 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
-------------------------------------------- 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 sander at oomvanlieshout.net Web www.oomvanlieshout.net/sander