Skip to content

Soil texture triangle in R?

8 messages · Sander Oom

Sander Oom
# 
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.
Sander Oom
# 
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)
 > }
Sander Oom
# 
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))
 > }
 >
Sander Oom
# 
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:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Sander Oom
    

    

      
      #
    

  


  

      
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:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Sander Oom
    

    

      
      #
    

  


  

      
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:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Sander Oom
    

    

      
      #
    

  


  

      
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.055,0.021,"Sand")
 >   text(0.49,0.18*sin60,"Loam")
 >   text(0.72,0.15*sin60,"Silt loam")
 >   text(0.9,0.06*sin60,"Silt")
 >  }
 >  if(is.null(pch)) pch<-1:nrow(soilprop)
 >  if(is.null(col)) col<-2:(nrow(soilprop)+1)
 >  points(1-soilprop[,2]+(soilprop[,2]-(1-soilprop[,3]))*0.5,
 >   soilprop[,1]*sin60,pch=pch,col=col)
 >  if(show.legend) {
 >   samplenames<-rownames(soilprop)
 > 
legend(0,0.8+0.05*length(samplenames),legend=samplenames,pch=pch,col=col)
 >  }
 >  par(oldpar)
 > }

  
  


  


      

    

    

      
      

        


  

        

      Sander Oom
    

    

      
      #
    

  


  

      
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.
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: plot.soiltexture.R
Url: https://stat.ethz.ch/pipermail/r-help/attachments/20050528/5e09fa86/plot.soiltexture.pl