提取矢量点位对应5×5窗口栅格像元值

锐多宝 · 公众号 · · 2025-03-07 23:47

正文

有学员自己写了一部分代码，运行不出结果，让我改改，下面是学员代码和修改过程，供大家学习参考。

任务要求

我有一个excle表格，里边包含一组带有编号、X、Y地理坐标（WGS84）的坐标点，另外给定一个含有空间参考坐标系和空间分辨率的栅格数据，想编写一个R语言代码，提取以表格中的坐标点为中心所对应的栅格数据5x5窗口中的像元值，然后存放到表格中相应的坐标点后。

学员代码

# 加载所需的包
library(readxl)
library(sf)
library(terra)
library(writexl)

# 读取Excel表格中的坐标点
coords <- read_excel("I:/2024年课题/地面样地/KT-1.xlsx")  # 替换为实际文件路径

# 将坐标点转换为sf对象（WGS84坐标系）
pts_sf <- st_as_sf(coords, coords = c("X", "Y"), crs = 4326)

# 读取栅格数据
raster <- rast("I:/2024年课题/无人机光谱指数/SAVI_KT1.tif")  # 替换为实际文件路径

# 将点转换到栅格数据的坐标系
pts_proj <- st_transform(pts_sf, crs(raster))

# 提取转换后的坐标点坐标
pts_coords <- st_coordinates(pts_proj)

# 获取每个点对应的中心像元行列号
cells_center <- cellFromXY(raster, pts_coords)
rows_center <- rowFromCell(raster, cells_center)
cols_center <- colFromCell(raster, cells_center)

# 获取栅格的行列总数
n_rows <- nrow(raster)
n_cols <- ncol(raster)

# 初始化存储结果的矩阵（每个点25个像元值）
values_matrix <- matrix(nrow = nrow(coords), ncol = 25)

# 遍历每个坐标点
for (i in 1:nrow(coords)) {
  r <- rows_center[i]
  c <- cols_center[i]

# 生成5x5窗口的行列号范围
  rows_window <- (r - 2):(r + 2)
  cols_window <- (c - 2):(c + 2)

# 生成所有行列组合并按行优先排序
  grid <- expand.grid(row = rows_window, col = cols_window)
  grid <- grid[order(grid$row, grid$col), ]

# 检查行列号是否在栅格范围内
  valid <- grid$row >= 1 & grid$row <= n_rows & grid$col >= 1 & grid$col <= n_cols

# 提取有效的像元值
  cells_valid <- cellFromRowCol(raster, grid$row[valid], grid$col[valid])
  vals_valid <- raster[cells_valid]

# 填充到结果矩阵（无效位置为NA）
  vals <- rep(NA, 25)
  vals[valid] <- vals_valid
  values_matrix[i, ] <- vals
}

# 将结果添加到原始数据框
colnames(values_matrix) <- paste0("Cell_", 1:25)
result_df <- cbind(coords, values_matrix)

# 保存结果到Excel文件
write_xlsx(result_df, "output.xlsx")  # 替换为期望的输出路径