From: Sander Oom <slist_at_oomvanlieshout.net>

Date: Sat 28 May 2005 - 21:58:27 EST

require(vcd)

)

oldpar <- par(no.readonly=TRUE)

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

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

text(0.5,0.57,"Clay", col=coltext)

# lower bound of clay loam & silty divider

segments(x1,y1,x2,y2, col=colgrid)

y1<-c(0.37*sin60,0.2*sin60,0.2*sin60) y2<-c(0.37*sin60,0.2*sin60,0.285*sin60)

}

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")

*>
*

*>
*

Date: Sat 28 May 2005 - 21:58:27 EST

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@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@oomvanlieshout.net 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.htmlReceived on Sat May 28 22:07:32 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:32:13 EST
*