R语言数字土壤制图：如何创建多土壤属性变异函数并执行克里金？

1. 环境准备与数据导入 #

首先安装并加载空间分析所需的核心包，同时注意你的CSV数据必须包含空间坐标列（如经度lon、纬度lat或投影后的XY坐标），否则无法开展空间分析。

# 首次运行安装所需包
install.packages(c("gstat", "sp", "raster", "ggplot2"))

# 加载包
library(gstat)
library(sp)
library(raster)
library(ggplot2)

# 导入CSV数据，替换为你的文件路径
soil_data <- read.csv("your_soil_data.csv", stringsAsFactors = FALSE)

# 将数据转换为空间点对象（假设坐标列是lon和lat，根据实际修改）
coordinates(soil_data) <- ~ lon + lat
# 若使用投影坐标（如UTM），设置对应坐标系，示例：
# proj4string(soil_data) <- CRS("+proj=utm +zone=50 +datum=WGS84")

2. 变异函数拟合 #

针对SAND、SILT、CLAY、OM四个属性，分别完成探索性分析与理论模型拟合。

2.1 探索性变异函数分析 #

先查看数据的空间自相关性特征：

# 以SAND为例，其他属性替换列名即可
sand_variogram <- variogram(SAND ~ 1, data = soil_data)
# 可视化探索性变异函数
plot(sand_variogram)

2.2 拟合理论变异函数 #

根据探索性图的趋势，选择合适的理论模型（球状、指数、高斯等）进行拟合：

# 拟合球状模型（若效果差可替换为"Exp"指数、"Gau"高斯）
sand_vgm_fit <- fit.variogram(sand_variogram, model = vgm("Sph"))
# 查看拟合参数
print(sand_vgm_fit)
# 叠加拟合模型到探索性图验证效果
plot(sand_variogram, model = sand_vgm_fit)

重复上述代码，将列名替换为SILT、CLAY、OM，完成四个属性的变异函数拟合。

3. 创建插值网格 #

生成覆盖研究区域的规则网格，作为克里金插值的目标载体：

# 基于样本数据范围创建网格，n为网格点数，按需调整
grid <- makegrid(soil_data, n = 10000)
colnames(grid) <- c("lon", "lat")
# 转换为空间网格对象，保持与样本数据一致的坐标系
coordinates(grid) <- ~ lon + lat
proj4string(grid) <- proj4string(soil_data)
# 转换为适合插值的像素格式
grid <- as(grid, "SpatialPixelsDataFrame")

4. 克里金插值计算 #

使用拟合好的变异函数执行普通克里金插值：

# SAND插值
sand_krige <- krige(SAND ~ 1, soil_data, grid, model = sand_vgm_fit)
# 转换为栅格对象便于制图
sand_raster <- raster(sand_krige)

# 重复处理其他属性
silt_vgm_fit <- fit.variogram(variogram(SILT ~1, soil_data), vgm("Sph"))
silt_krige <- krige(SILT ~1, soil_data, grid, model = silt_vgm_fit)
silt_raster <- raster(silt_krige)

clay_vgm_fit <- fit.variogram(variogram(CLAY ~1, soil_data), vgm("Sph"))
clay_krige <- krige(CLAY ~1, soil_data, grid, model = clay_vgm_fit)
clay_raster <- raster(clay_krige)

om_vgm_fit <- fit.variogram(variogram(OM ~1, soil_data), vgm("Sph"))
om_krige <- krige(OM ~1, soil_data, grid, model = om_vgm_fit)
om_raster <- raster(om_krige)

5. 制作土壤空间分布图 #

用ggplot2或raster包完成可视化，也可导出栅格文件：

# 以SAND为例，用ggplot2绘制分布图
sand_df <- as.data.frame(sand_raster, xy = TRUE)
ggplot(sand_df, aes(x = x, y = y, fill = layer)) +
  geom_raster() +
  scale_fill_viridis_c(option = "plasma", name = "Sand Content (%)") +
  coord_equal() +
  theme_minimal() +
  labs(x = "Longitude", y = "Latitude", title = "Spatial Distribution of Sand")

# 保存栅格结果到本地（可选）
writeRaster(sand_raster, "sand_distribution.tif", overwrite = TRUE)

关键注意事项 #

• 坐标系统：若使用经纬度数据，建议先转换为投影坐标（如UTM），避免距离计算的偏差。
• 模型优化：若默认模型拟合效果差，可尝试调整vgm()的参数，或通过krige.cv()做交叉验证评估插值精度：

# SAND的交叉验证示例
sand_cv <- krige.cv(SAND ~1, soil_data, model = sand_vgm_fit)
# 查看误差统计（均值误差、均方根误差等）
summary(sand_cv$residuals)

内容的提问来源于stack exchange，提问作者Dimitris K