R：使用 ggtern 的 Shepard 图

Question

我有一些使用 ggtern 包为 R 详细阐述三元图的良好经验。

不幸的是，我仍然无法找到重现 Shepard's diagram 的限制的脚本，这是一个基于沙子、淤泥和粘土百分比的沉积物纹理分类图。到目前为止，我发现的只是soil classification，这不符合我对海洋沉积物分类的兴趣。

有谁知道我在哪里可以找到或怎么做？

谢谢！

Answer 1

感谢本的帮助！！

另外，我通过电子邮件联系了 ggtern，他们在几个小时内发布了美国农业部适用于 Shepard 分类的图表，这里是 link。

从ggtern website 中提取的以下代码可以很好地再现 Shepard 的图表

#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()