如何创建土壤质地三元图?

How to create soil texture ternary plot?

我有一个包含土壤质地信息的数据框,我想创建一个类似于 this 的土壤质地图。

我该如何使用:

  1. 每个区域有不同的颜色并且
  2. 三角形内 类 的标签,例如 this ?

示例数据:

area <- c('S1','S2','S3','S4','S5','S6','S7','S8')
sand <- c(76.4,56.9,61.7,64.5,71,70.1,60.5,53.7)
silt<-c(9.3,23.1,23,17.4,13.5,13.4,21.1,30.6)
clay<-c(14.3, 20,15.4,18,15.5,16.6,18.4,15.7)
my_data<-data.frame(area,sand,silt,clay)

ggtern 包的基本情节:

theme_set(theme_bw())
my_data %>%
    ggtern(aes(
        x = sand,
        y = clay,
        z = silt )) +
    geom_point()
library(tidyverse)
library(ggtern)
#> Registered S3 methods overwritten by 'ggtern':
#>   method           from   
#>   grid.draw.ggplot ggplot2
#>   plot.ggplot      ggplot2
#>   print.ggplot     ggplot2
#> --
#> Remember to cite, run citation(package = 'ggtern') for further info.
#> --
#> 
#> Attaching package: 'ggtern'
#> The following objects are masked from 'package:ggplot2':
#> 
#>     aes, annotate, ggplot, ggplot_build, ggplot_gtable, ggplotGrob,
#>     ggsave, layer_data, theme_bw, theme_classic, theme_dark,
#>     theme_gray, theme_light, theme_linedraw, theme_minimal, theme_void

area <- c("S1", "S2", "S3", "S4", "S5", "S6", "S7", "S8")
sand <- c(76.4, 56.9, 61.7, 64.5, 71, 70.1, 60.5, 53.7)
silt <- c(9.3, 23.1, 23, 17.4, 13.5, 13.4, 21.1, 30.6)
clay <- c(14.3, 20, 15.4, 18, 15.5, 16.6, 18.4, 15.7)
my_data <- data.frame(area, sand, silt, clay)

theme_set(theme_bw())

my_data %>%
  ggtern(aes(
    x = sand,
    y = clay,
    z = silt,
    color = area
  )) +
  geom_point() +
  labs(
    x = "Axis sand",
    y = "Axis clay",
    z = "Axis silt"
  )

reprex package (v2.0.1)

于 2021-10-04 创建

你的数据:

library(plotrix)

soil.texture(my_data[,2:4],col.symbols=1:8,bg.symbols=1:8,point.labels=my_data$area,pch=21)
legend( x=1, 
        legend=my_data$area,
        col=1:8, 
        fill=1:8 )

这给你:

编辑

要避免标签背景,您必须查看 soil.texture 的源代码。如果这样做,您会注意到该函数使用 boxed.labels 来绘制三角形图中的文本,因此您无法通过使用函数的参数来避免这种行为。相反,您可以编写自己的函数版本,将 boxed.labels 更改为更简单的 text 函数,如下所示:

my_soil_texture <- function (soiltexture = NULL, main = "", at = seq(0.1, 0.9, 
                                                  by = 0.1), axis.labels = c("percent sand", "percent silt", 
                                                                             "percent clay"), tick.labels = list(l = seq(10, 90, by = 10), 
                                                                                                                 r = seq(10, 90, by = 10), b = seq(10, 90, by = 10)), show.names = TRUE, 
          show.lines = TRUE, col.names = "gray", bg.names = par("bg"), 
          show.grid = FALSE, col.axis = "black", col.lines = "gray", 
          col.grid = "gray", lty.grid = 3, show.legend = FALSE, label.points = FALSE, 
          point.labels = NULL, col.symbols = "black", pch = par("pch"), 
          ...) 
{
  par(xpd = TRUE)
  triax.plot(x = NULL, main = main, at = at, axis.labels = axis.labels, 
             tick.labels = tick.labels, col.axis = col.axis, show.grid = show.grid, 
             col.grid = col.grid, lty.grid = lty.grid)
  arrows(0.12, 0.41, 0.22, 0.57, length = 0.15)
  arrows(0.78, 0.57, 0.88, 0.41, length = 0.15)
  arrows(0.6, -0.1, 0.38, -0.1, length = 0.15)
  if (show.lines) {
    triax.segments <- function(h1, h3, t1, t3, col) {
      segments(1 - h1 - h3/2, h3 * sin(pi/3), 1 - t1 - 
                 t3/2, t3 * sin(pi/3), col = col)
    }
    h1 <- c(85, 70, 80, 52, 52, 50, 20, 8, 52, 45, 45, 65, 
            45, 20, 20)/100
    h3 <- c(0, 0, 20, 20, 7, 0, 0, 12, 20, 27, 27, 35, 40, 
            27, 40)/100
    t1 <- c(90, 85, 52, 52, 43, 23, 8, 0, 45, 0, 45, 45, 
            0, 20, 0)/100
    t3 <- c(10, 15, 20, 7, 7, 27, 12, 12, 27, 27, 55, 35, 
            40, 40, 60)/100
    triax.segments(h1, h3, t1, t3, col.lines)
  }
  if (show.names) {
    xpos <- c(0.5, 0.7, 0.7, 0.73, 0.73, 0.5, 0.275, 0.275, 
              0.27, 0.27, 0.25, 0.135, 0.18, 0.055, 0.49, 0.72, 
              0.9)
    ypos <- c(0.66, 0.49, 0.44, 0.36, 0.32, 0.35, 0.43, 
              0.39, 0.3, 0.26, 0.13, 0.072, 0.032, 0.024, 0.18, 
              0.15, 0.06) * sin(pi/3)
    snames <- c("clay", "silty", "clay", "silty clay", "loam", 
                "clay loam", "sandy", "clay", "sandy clay", "loam", 
                "sandy loam", "loamy", "sand", "sand", "loam", "silt loam", 
                "silt")
    text(xpos, ypos, labels=snames)# here I switched from boxed.labels(xpos, ypos, snames, border = FALSE, xpad = 0.5) to text(xpos, ypos, labels=snames)
  }
  par(xpd = FALSE)
  if (is.null(soiltexture)) 
    return(NULL)
  soilpoints <- triax.points(soiltexture, show.legend = show.legend, 
                             label.points = label.points, point.labels = point.labels, 
                             col.symbols = col.symbols, pch = pch, ...)
  invisible(soilpoints)
}

所以现在你可以只使用你自己的函数来绘制了:

my_soil_texture(my_data[,2:4],col.symbols=1:8,bg.symbols=1:8,point.labels=my_data$area,pch=21,col.grid=3)
legend( x=1, 
        legend=my_data$area,
        col=1:8, 
        fill=1:8 )

这给你:

我想强调刚刚所做的更改的一个小问题。由于您可以更改主三角形的背景颜色,因此 boxed labels() 尝试根据背景颜色确定白色或黑色文本是否更易于阅读,并相应地显示文本。因此,如果您想更改情节的背景,从 boxed.labels() 切换到 text() 可能会导致问题