Re: [R] Soil texture triangle in R?

From: Sander Oom <slist_at_oomvanlieshout.net>
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.html
Received 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