在terra包中实现无监督分类：K-means与聚类指标计算问题

1. 解决K-means报错问题 #

你的报错核心原因是直接将SpatRaster对象传入kmeans()——即使使用na.omit()处理，返回的仍是栅格对象，而base包的kmeans()仅支持数值矩阵或数据框输入。正确操作是先提取栅格数值并移除缺失值：

library(terra)
library(vegan)

set.seed(99)
# 提取栅格数值，移除含NA的行
clim_values <- na.omit(values(Clim))
# 运行K-means聚类
kmncluster <- kmeans(clim_values, centers = 5, iter.max = 500, nstart = 5)

2. 将聚类结果转回栅格 #

无需单独创建带单元格信息的数据框，直接基于原栅格的空间属性生成结果栅格：

# 复制原栅格的空间框架（取第一个波段即可）
cluster_raster <- Clim[[1]]
# 将聚类结果赋值给非NA的栅格单元格
cluster_raster[!is.na(Clim[[1]])] <- kmncluster$cluster

# 可视化结果
plot(cluster_raster, main = "K-means聚类结果")

3. 实现类似cascadeKM的最优簇数选择 #

直接用vegan::cascadeKM()处理提取的数值矩阵（无需将整个栅格转成数据框保存），同时计算Calinski或SSI指标：

方法1：用cascadeKM批量测试簇数 #

# 测试2-5类，基于Calinski指标
cascade_result <- cascadeKM(clim_values, 
                            inf.gr = 2, 
                            sup.gr = 5, 
                            iter = 100, 
                            criterion = "calinski")

# 查看结果（最优簇数对应最高的Calinski值）
print(cascade_result)
plot(cascade_result)

方法2：自定义循环计算指标（灵活控制） #

如果需要同时计算SSI指标，可以手动循环不同簇数：

k_range <- 2:5
calinski_scores <- numeric(length(k_range))
ssi_scores <- numeric(length(k_range))

for (idx in seq_along(k_range)) {
  k <- k_range[idx]
  km_fit <- kmeans(clim_values, centers = k, iter.max = 500, nstart = 5)
  
  # 计算Calinski-Harabasz指数
  calinski_scores[idx] <- calinski(km_fit$cluster, clim_values)
  
  # 计算Simple Structure Index（SSI）
  ssi_scores[idx] <- vegan:::ssi(km_fit$centers, clim_values, km_fit$cluster)
}

# 输出指标结果
data.frame(
  簇数 = k_range,
  Calinski指数 = calinski_scores,
  SSI指数 = ssi_scores
)

关键注意事项 #

• 始终用values(Clim)提取栅格的数值矩阵，这是连接terra与统计函数的核心步骤，无需将结果保存为单独的数据框。
• na.omit()只对数值矩阵生效，直接用于SpatRaster只会返回移除了NA单元格的栅格，而非适合统计计算的数值格式。

内容的提问来源于stack exchange，提问作者Derek Corcoran