[R] Soil texture triangle in R?

Sander Oom slist at oomvanlieshout.net
Fri May 27 18:33:29 CEST 2005


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




More information about the R-help mailing list