I had some good experience elaborating ternary plots using ggtern package for R.
Unfortunately I still haven't been able to find a script for reproducing the limits from Shepard's diagram, a diagram for textural classification of sediments based on percentages of Sand, Silt and Clay. So far all I found was soil classification which doesn't fit my interest of classifying marine sediments.
Does any one know where I can find or how to do it?
Thanks!
Thanks for the help Ben!!
Aditionally I contacted ggtern trhough e-mail and they posted within hours an adaptation of USDA diagram for Shepard classification, here is the link.
The following code, extracted from ggtern website, produces a really nice reproduction of Shepard's diagram
#Build a library of points, left to right, top to bottom...
points <- data.frame(
rbind(c( 1,1.000,0.000,0.000),
c( 2,0.750,0.250,0.000),
c( 3,0.750,0.125,0.125),
c( 4,0.750,0.000,0.250),
c( 5,0.600,0.200,0.200),
c( 6,0.500,0.500,0.000),
c( 7,0.500,0.000,0.500),
c( 8,0.400,0.400,0.200),
c( 9,0.400,0.200,0.400),
c(10,0.250,0.750,0.000),
c(11,0.250,0.000,0.750),
c(12,0.200,0.600,0.200),
c(13,0.200,0.400,0.400),
c(14,0.200,0.200,0.600),
c(15,0.125,0.750,0.125),
c(16,0.125,0.125,0.750),
c(17,0.000,1.000,0.000),
c(18,0.000,0.750,0.250),
c(19,0.000,0.500,0.500),
c(20,0.000,0.250,0.750),
c(21,0.000,0.000,1.000)
)
)
colnames(points) = c("IDPoint","T","L","R")
#Give each Polygon a number
polygon.labels <- data.frame(
Label=c("Clay",
"Sandy Clay",
"Silty Clay",
"Sand + Silt + Clay",
"Clayey Sand",
"Clayey Silt",
"Sand",
"Silty Sand",
"Sandy Silt",
"Silt"))
#Assign each label an index
polygon.labels$IDLabel=1:nrow(polygon.labels)
#Create a map of polygons to points
polygons <- data.frame(
rbind(c(1,1),c(1,2),c(1,4),
c(2,6),c(2,2),c(2,3),c(2,5),c(2,8),
c(3,3),c(3,4),c(3,7),c(3,9),c(3,5),
c(4,5),c(4,14),c(4,12),
c(5,6),c(5,8),c(5,12),c(5,15),c(5,10),
c(6,7),c(6,11),c(6,16),c(6,14),c(6,9),
c(7,17),c(7,10),c(7,18),
c(8,15),c(8,12),c(8,13),c(8,19),c(8,18),
c(9,13),c(9,14),c(9,16),c(9,20),c(9,19),
c(10,11),c(10,21),c(10,20)
)
)
#IMPORTANT FOR CORRECT ORDERING.
polygons$PointOrder <- 1:nrow(polygons)
#Rename the columns
colnames(polygons) = c("IDLabel","IDPoint","PointOrder")
#Merge the three sets together to create a master set.
df <- merge(polygons,points)
df <- merge(df,polygon.labels)
df <- df[order(df$PointOrder),]
#Determine the Labels Data library(plyr)
Labs = ddply(df,"Label",function(x){c(c(mean(x$T),mean(x$L),mean(x$R)))})
colnames(Labs) = c("Label","T","L","R")
#Build the final plot
library(ggtern)
base <- ggtern(data=df,aes(L,T,R)) +
geom_polygon(aes(fill=Label,group=Label),color="black",alpha=0.25) +
geom_text(data=Labs,aes(label=Label),size=4,color="black") +
theme_bw() +
custom_percent("Percent") +
labs(title="Shepard Sediment Classification Diagram",
fill = "Classification",
T="Clay",
L="Sand",
R="Silt")
print(base) #to console
#Render to file.
png("plot.png",width=800,height=600)
print(base)
dev.off()