空间栅格栈K-Means聚类异常问题求解及替代方案咨询

问题描述 #

我尝试结合不同空间数据集，寻找这些因子组合的相似模式，目前使用terra包处理空间栅格，尝试通过K-Means或BIC聚类实现目标。

已编写如下代码：

#I am stacking the different rasters together (they have the same extent and resolution)
stack1 <- c(resamp.pop, eleva.rec.tif, resamp.lc, resamp.lc.bare)
#need to make matrix to create a cluster
val <- values(stack1)
#do the cluster
k.clust <- kmeans(na.omit(val), centers = 10)

#create raster with same dimensions as stack1
krast <- stack1
values(krast) <- k.clust$cluster
plot(krast)

但聚类结果呈现异常随机模式，未按因子组合有效聚类。现寻求改进方法，或可实现空间组合数据K-Means/BIC聚类的替代包推荐。

一、现有代码的改进方法 #

• 数据标准化：不同栅格的数值量级差异会严重干扰K-Means聚类（比如人口数据可能是千级，海拔是百级，土地覆盖是分类编码）。先对每个栅格图层做标准化处理：

  # 对堆叠后的栅格每个图层做标准化
  stack_norm <- scale(stack1)
  val_norm <- values(stack_norm)
  

• 修复缺失值映射问题：原代码中na.omit(val)丢弃含NA的行后，直接赋值回栅格会导致缺失值位置被错误填充，聚类结果混乱。正确做法是保留原索引：

  val <- values(stack1)
  # 记录非NA的位置
  non_na_idx <- !is.na(rowSums(val))
  val_non_na <- val[non_na_idx, ]
  # 标准化后聚类
  k.clust <- kmeans(scale(val_non_na), centers = 10)
  # 初始化聚类结果向量，填充NA
  clust_result <- rep(NA, nrow(val))
  clust_result[non_na_idx] <- k.clust$cluster
  # 用单个图层做模板赋值更高效
  krast <- stack1[[1]]
  values(krast) <- clust_result
  plot(krast)
  

• 合理选择聚类数：固定centers=10可能不符合数据实际分布，可通过肘部法则或BIC准则确定最优聚类数：

  # 肘部法则示例：遍历聚类数，计算组内平方和
  wss <- sapply(1:15, function(k) {
    kmeans(scale(val_non_na), centers = k)$tot.withinss
  })
  # 绘制肘部图，找拐点对应的聚类数
  plot(1:15, wss, type="b", xlab="聚类数", ylab="组内平方和")
  

• 加入空间约束：普通K-Means不考虑空间位置信息，易导致结果随机。可加入空间坐标作为聚类变量，或使用空间权重约束：

  # 示例：将栅格单元格坐标加入聚类变量
  coords <- xyFromCell(stack1, 1:ncell(stack1))
  val_with_coords <- cbind(val_non_na, coords[non_na_idx, ])
  # 标准化后再聚类
  k.clust_spatial <- kmeans(scale(val_with_coords), centers = 10)
  

二、替代包推荐 #

• spatstat：支持空间点、栅格数据的聚类分析，提供多种空间约束聚类方法，适配带空间依赖的数据集。
• fpc：包含密度聚类（DBSCAN）、层次聚类等算法，支持结合空间信息的聚类评估，可验证K-Means结果的合理性。
• clustGeo：专门针对空间数据的聚类工具，能同时平衡变量相似性与空间邻近性，通过权重调节适配不同场景，非常适合因子组合+空间模式的分析需求。
• RStoolbox：封装了遥感栅格数据的聚类流程，内置标准化、聚类数选择等功能，操作简洁，适合批量处理空间数据。

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