如何用R语言从比重计+筛法粒径分析数据计算土壤粒组含量？

前提假设 #

假设你的150份样品数据已整理为结构化表格（如CSV格式），包含三列：

• sample_id：样品唯一标识
• d：粒径（单位：μm）
• p_pass：小于对应粒径的颗粒占比（%）

代码实现 #

# 加载数据处理包
library(tidyverse)

# 读取数据（替换为你的文件路径）
soil_data <- read_csv("your_soil_data.csv")

# 批量计算每个样品的粒级含量
texture_results <- soil_data %>%
  # 按样品分组处理
  group_by(sample_id) %>%
  # 确保粒径按升序排列（插值的必要前提）
  arrange(d) %>%
  # 线性插值计算2μm和50μm对应的通过率
  mutate(
    p_2um = approx(x = d, y = p_pass, xout = 2)$y,
    p_50um = approx(x = d, y = p_pass, xout = 50)$y
  ) %>%
  # 保留每组唯一计算结果
  distinct(sample_id, p_2um, p_50um) %>%
  # 计算各粒级含量
  mutate(
    clay = p_2um,          # 黏粒(<2μm)含量
    silt = p_50um - p_2um, # 粉粒(2μm-50μm)含量
    sand = 100 - p_50um    # 砂粒(>50μm)含量
  ) %>%
  ungroup()

# 查看结果
print(texture_results)

# 导出结果到CSV
write_csv(texture_results, "soil_texture_results.csv")

关键说明 #

• 若粒径曲线非线性程度较高，可将approx()替换为spline()函数（仅需修改代码中的函数名），用样条插值提升精度。
• 需保证每个样品的粒径数据覆盖2μm和50μm范围（至少包含一个小于2μm的点和一个大于50μm的点），否则插值会返回NA值。
• 若样品为单独文件（每个样品对应一个CSV），可批量读取处理：

# 批量读取单个样品文件
file_paths <- list.files(path = "your_sample_folder", pattern = "*.csv", full.names = TRUE)

# 循环处理所有文件
texture_results <- map_dfr(file_paths, function(file) {
  sample_data <- read_csv(file)
  # 从文件名提取样品ID，可根据实际命名规则调整
  sample_id <- str_remove(basename(file), ".csv")
  
  # 插值计算目标粒径通过率
  p_2um <- approx(x = sample_data$d, y = sample_data$p_pass, xout = 2)$y
  p_50um <- approx(x = sample_data$d, y = sample_data$p_pass, xout = 50)$y
  
  # 组装结果
  tibble(
    sample_id = sample_id,
    clay = p_2um,
    silt = p_50um - p_2um,
    sand = 100 - p_50um
  )
})

# 导出批量处理结果
write_csv(texture_results, "batch_soil_texture_results.csv")

内容的提问来源于stack exchange，提问作者Thomas Dupuis