#Construction of Shepards Diagram
library(ggtern)
library(plyr)
#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"))
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)
)
)
polygons$PointOrder <- 1:nrow(polygons)
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
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
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()