GSIF包buffer.dist()函数报'subscript out of bounds'错误求修复

云上老周

2026-8-5

问题描述

我想要使用Tomislav Hengl等人（2018）开发的GSIF包中的buffer.dist()函数，该包自2019年未更新且已从CRAN下架。根据作者说明，该函数可生成每个采样点的欧氏距离栅格砖，作为随机森林模型的空间协变量以解决预测中的空间自相关问题。由于该包已无法获取，我无法知晓预期结果。

我从CRAN归档库下载了GSIF v0.5-5（2019-01-04）版本，手动将函数加载到R工作区（R版本4.2.1，运行于macOS Big Sur 11.6）。按照教程流程，我加载gstat包的meuse示例数据，转换为SpatialPointsDataFrame和SpatialPixelsDataFrame后调用函数：

# 手动加载GSIF环境（从CRAN仓库手动下载）
source("AAAA.R") # 需要先加载
# 手动加载buffer.dist()函数
source("buffer.dist.R")

# 加载依赖包
library(sp)
library(plotKML)
library(raster)
library(gstat)

## 按照教程流程操作
# 加载gstat包中的示例数据
data(meuse, echo = FALSE)
data(meuse.grid)

# 转换为SpatialPoints对象（buffer.dist()的输入数据要求）
meuse.sp <- SpatialPointsDataFrame(meuse[1:2], meuse[3:14], proj4string = CRS('+init=epsg:4326'))
meuse.grid.spdf <- SpatialPixelsDataFrame(meuse.grid[1:2], meuse.grid[6], proj4string = CRS('+init=epsg:4326'))

# 为每个单独的点生成缓冲距离
grid.dist0 <- buffer.dist(meuse.sp["zinc"], 
                          meuse.grid.spdf[1],
                          as.factor(1:nrow(meuse.sp)))

运行后出现报错：Error in x@coords[i, , drop = FALSE] : subscript out of bounds。手动排查发现错误出现在buffer.dist()函数的s <- s[predictionDomain@grid.index,]行。该函数代码如下：

setMethod("buffer.dist", signature(observations = "SpatialPointsDataFrame", predictionDomain = "SpatialPixelsDataFrame"), function(observations, predictionDomain, classes, width, ...){
  if(missing(width)){ width <- sqrt(areaSpatialGrid(predictionDomain)) }
  if(!length(classes)==length(observations)){ stop("Length of 'observations' and 'classes' does not match.") }
  ## 移除没有对应点的类别：
  xg = summary(classes, maxsum=length(levels(classes)))
  selg.levs = attr(xg, "names")[xg > 0]
  if(length(selg.levs)<length(levels(classes))){
    fclasses <- as.factor(classes)
    fclasses[which(!fclasses %in% selg.levs)] <- NA
    classes <- droplevels(fclasses)
  }
  ## 生成缓冲距离
  s <- list(NULL)
  for(i in 1:length(levels(classes))){
    s[[i]] <- raster::distance(rasterize(observations[which(classes==levels(classes)[i]),1]@coords, y=raster(predictionDomain)), width=width, ...)
  }
  s <- s[sapply(s, function(x){!is.null(x)})]
  s <- brick(s)
  s <- as(s, "SpatialPixelsDataFrame")
  s <- s[predictionDomain@grid.index,]
  return(s)
})

修复建议

错误根源

报错的核心原因是：将栅格砖转为SpatialPixelsDataFrame后，其内部的网格索引与输入的predictionDomain的grid.index不匹配，导致下标越界。这是因为新版sp/raster包对空间对象的结构处理逻辑，和旧版GSIF包开发时的逻辑存在差异。

具体修复方案

1. 替换索引匹配逻辑（最直接）

把函数中的s <- s[predictionDomain@grid.index,]替换为sp包的标准空间子集化语法，自动对齐空间网格：

s <- s[predictionDomain, ]

这种方式会基于空间范围和网格参数自动匹配数据，完全规避旧索引值不匹配的问题。

2. 确保栅格转换时的参数对齐

在将brick转为SpatialPixelsDataFrame前，强制匹配predictionDomain的栅格参数，避免分辨率或投影的细微差异：

s <- brick(s)
# 强制对齐predictionDomain的栅格分辨率、范围和投影
s <- projectRaster(s, raster(predictionDomain))
s <- as(s, "SpatialPixelsDataFrame")
# 使用标准方式子集化
s <- s[predictionDomain, ]

3. 简化循环内的栅格化代码

原循环中直接提取@coords的写法容易出错，改为直接使用空间对象：

# 替换循环内的代码
for(i in 1:length(levels(classes))){
  obs_subset <- observations[classes == levels(classes)[i], ]
  s[[i]] <- raster::distance(rasterize(obs_subset, y=raster(predictionDomain)), width=width, ...)
}

直接操作SpatialPointsDataFrame对象，减少底层结构调用的出错风险。

4. 兼容新版R的投影语法

原代码中CRS('+init=epsg:4326')在新版R中会触发警告，建议替换为兼容的写法：

# 替换投影定义
proj4string(meuse.sp) <- CRS(SRS_string = "EPSG:4326")
proj4string(meuse.grid.spdf) <- CRS(SRS_string = "EPSG:4326")

虽然这不是直接报错原因，但能避免潜在的空间对象兼容性问题。

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

GSIF包buffer.dist()函数报'subscript out of bounds'错误求修复

错误根源 #

具体修复方案 #

1. 替换索引匹配逻辑（最直接） #

2. 确保栅格转换时的参数对齐 #

3. 简化循环内的栅格化代码 #