在R中为数据框坐标计算TPI时遇resolution函数报错问题

报错原因分析 #

你遇到的could not find function "resolution"错误主要有两个原因：

• 包冲突与函数命名差异：resolution()是terra包的函数，但代码同时加载了raster和terra，且未明确dem的对象类型。如果dem是raster包的RasterLayer，则应使用res()而非resolution()。
• 函数环境问题：focal函数内部的匿名函数无法正确访问外部包的函数，尤其是未提前加载或未明确命名空间的情况下。

此外，代码还存在几个逻辑错误，会导致TPI计算结果不正确：

• 窗口尺寸错误：你指定的是5×5和10×10窗口，但代码中win_size <- scale *5会生成25×25和50×50的窗口，完全不符合需求。
• 误用点高程数据：在focal函数中使用了提前提取的elev（177个点的高程），而非当前窗口中心单元格的高程，这会导致TPI计算逻辑完全错误。

修正后的代码 #

以下是基于terra包（替代老旧的raster/sp包）的修正代码，解决了所有问题：

library(terra)
library(haven)

# 构建坐标数据框
df <- data.frame(lon = coords$longitude, lat = coords$latitude)

# 转换为terra的空间点对象（替代SpatialPointsDataFrame）
pts <- vect(df, geom = c("lon", "lat"), crs = "EPSG:4326")

# 提取兴趣点的高程值
elev <- extract(dem, pts)

# 定义TPI计算的窗口尺度（5×5和10×10）
scales <- c(5, 10)

# 预计算DEM分辨率（避免在focal内部调用函数引发环境问题）
dem_res <- resolution(dem)[1] # 假设DEM为方形单元格，取x方向分辨率

# 循环计算不同尺度的TPI
tpi_list <- list()
for (scale in scales) {
  # 创建指定尺寸的窗口矩阵
  win <- matrix(1, nrow = scale, ncol = scale)
  
  # 计算TPI：中心单元格高程减去窗口平均高程，再除以分辨率
  tpi <- focal(dem, w = win, fun = function(x) {
    if (scale %% 2 == 1) {
      # 奇数窗口：取单个中心单元格的高程
      central_idx <- (length(x) + 1) %/% 2
      central_val <- x[central_idx]
      # 若需排除中心单元格计算平均，使用mean(x[-central_idx])
      tpi_val <- (central_val - mean(x)) / dem_res
    } else {
      # 偶数窗口：取中心4个单元格的平均高程（偶数窗口无单一中心，可根据需求调整）
      half <- scale %/% 2
      center_pos <- c((half*scale)+half, (half*scale)+half+1, 
                      ((half+1)*scale)+half, ((half+1)*scale)+half+1)
      central_val <- mean(x[center_pos])
      # 若需排除中心4格计算平均，使用mean(x[-center_pos])
      tpi_val <- (central_val - mean(x)) / dem_res
    }
    return(tpi_val)
  })
  
  # 提取兴趣点的TPI值
  tpi_vals <- extract(tpi, pts)
  
  # 将结果存入列表
  tpi_list[[as.character(scale)]] <- tpi_vals
}

# 合并不同尺度的TPI结果为数据框
tpi_df <- do.call(cbind, tpi_list)
rownames(tpi_df) <- rownames(df)
colnames(tpi_df) <- paste0("TPI_", scales)

关键修正说明 #

1. 统一使用terra包：terra是raster/sp的替代包，功能更完善且避免了包冲突问题。
2. 预计算分辨率：提前获取DEM分辨率并存储为变量，避免在focal内部调用函数引发环境错误。
3. 修正窗口尺寸：直接使用指定的5和10作为窗口大小，符合你的研究需求。
4. 正确计算TPI：在focal函数内部从窗口数据中获取中心单元格的高程，而非误用提前提取的点高程数据。
5. 处理偶数窗口：针对10×10的偶数窗口，通过取中心4个单元格的平均值作为中心高程（可根据研究需求调整逻辑）。

额外提示 #

• 如果你的dem是raster包的RasterLayer，可先转换为terra对象：dem <- rast(dem)。
• 若需计算"周围单元格"（排除中心）的平均高程，将代码中的mean(x)替换为mean(x[-central_idx])（奇数窗口）或mean(x[-center_pos])（偶数窗口）。
• 学术研究中TPI通常使用奇数尺寸的窗口（如5×5、9×9），若你并非特意使用偶数窗口，建议将scales改为c(5,9)或c(5,11)。

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