有学员自己写了一部分代码，运行不出结果，让我改改，下面是学员代码和修改过程，供大家学习参考。
任务要求
我有一个excle表格，里边包含一组带有编号、X、Y地理坐标（WGS84）的坐标点，另外给定一个含有空间参考坐标系和空间分辨率的栅格数据，想编写一个R语言代码，提取以表格中的坐标点为中心所对应的栅格数据5x5窗口中的像元值，然后存放到表格中相应的坐标点后。
学员代码
# 加载所需的包
library(readxl)
library(sf)
library(terra)
library(writexl)
# 读取Excel表格中的坐标点
coords "I:/2024年课题/地面样地/KT-1.xlsx") # 替换为实际文件路径
# 将坐标点转换为sf对象（WGS84坐标系）
pts_sf "X", "Y"), crs = 4326)
# 读取栅格数据
raster "I:/2024年课题/无人机光谱指数/SAVI_KT1.tif") # 替换为实际文件路径
# 将点转换到栅格数据的坐标系
pts_proj
# 提取转换后的坐标点坐标
pts_coords
# 获取每个点对应的中心像元行列号
cells_center
rows_center
cols_center
# 获取栅格的行列总数
n_rows
n_cols
# 初始化存储结果的矩阵（每个点25个像元值）
values_matrix
# 遍历每个坐标点
for (i in 1:nrow(coords)) {
r
c
# 生成5x5窗口的行列号范围
rows_window
cols_window
# 生成所有行列组合并按行优先排序
grid
grid $row, grid$col), ]
# 检查行列号是否在栅格范围内
valid $row >= 1 & grid$row <= n_rows & grid$col >= 1 & grid$col <= n_cols
# 提取有效的像元值
cells_valid $row[valid], grid$col[valid])
vals_valid
# 填充到结果矩阵（无效位置为NA）
vals
vals[valid]
values_matrix[i, ]
}
# 将结果添加到原始数据框
colnames(values_matrix) "Cell_", 1:25)
result_df
# 保存结果到Excel文件
write_xlsx(result_df, "output.xlsx") # 替换为期望的输出路径
代码点评
上面的代码无法正常执行，是因为出现了一个典型问题，R包的混用问题。
terra
包和
sf
包不能混用！
-
terra
包有完整的空间分析功能，包括矢量和栅格
-
矢量：
vect()
函数，生成SpatVector
-
栅格：
rast()
函数，生成SpatRaster
要想提取矢量对应的栅格值，注意矢量一定要使用
vect()
函数，读取后获得的SpatVector，不能使用
sf
包，
sf
包生成的sf对象，terra包是无法处理的，无法和SpatRaster进行叠置分析。
代码修改
发现了上面的问题，我就把代码进行了修改，把多余的包去掉，使用
terra
包函数来做，结合DeepSeek辅助编程。
改了一部分代码，把空间分析完全使用
terra
包来写，然后把情况完整的描述给DeepSeek，让它帮我完成代码：
# 加载所需的包
library(readxl)
library(terra)
library(openxlsx)
# 读取Excel表格中的坐标点
coords "./KT1.xlsx") # 替换为实际文件路径
coords
# 将坐标点转换为sf对象（WGS84坐标系）
pts_sf "X", "Y"
), crs = "EPSG:4326")
pts_sf
plot(pts_sf)
# 读取栅格数据
raster "./NDVI_KT1.tif") # 替换为实际文件路径
raster
plot(raster)
# 将点转换到栅格数据的坐标系
rast_4326 = project(raster, "EPSG:4326")
rast_4326
已有代码如上所示，下面有一个任务，提取pts_sf点位对应的rast_4326栅格数据5x5窗口中的像元值，生成一个新的数据框，coords_rast_df，里面有ID，X,Y，还有对应的25个像元值的值，V1,V2一直到V25，请补全代码
对话过程
对话过程太长了，放到公众号里不合适，感兴趣的同学扫码看吧。
对话过程
完整代码
下面是实际可以运行的完整代码
# 我有一个excle表格，里边包含一组带有编号、X、Y地理坐标（WGS84）的坐标点，另外给定一个含有空间参考坐标系和空间分辨率的栅格数据，想编写一个R语言代码，提取以表格中的坐标点为中心所对应的栅格数据5x5窗口中的像元值，然后存放到表格中相应的坐标点后。
# 加载所需的包
library(readxl)
library(terra)
library(openxlsx)
# 读取Excel表格中的坐标点
coords "./KT1.xlsx") # 替换为实际文件路径
coords
# 将坐标点转换为sf对象（WGS84坐标系）
pts_sf "X", "Y"), crs = "EPSG:4326")
pts_sf
plot(pts_sf)
# 读取栅格数据
raster "./NDVI_KT1.tif") # 替换为实际文件路径
raster
plot(raster)
# 将点转换到栅格数据的坐标系
rast_4326 = project(raster, "EPSG:4326")
rast_4326
pts_coords # 获取坐标矩阵
cells # 获取像元索引
row_col # 转换为行列号
# 初始化存储25个像元值的列表
values_list
# 遍历每个点提取窗口值
for (i in 1:nrow(row_col)) {
# 当前点的中心行列号
center_row
center_col
# 生成5x5窗口的行列范围
rows
cols
# 过滤有效行列（处理边界情况）
valid_rows = 1 & rows <= nrow(rast_4326)]
valid_cols = 1 & cols <= ncol(rast_4326)]
# 跳过完全超出范围的窗口
if (length(valid_rows) == 0 || length(valid_cols) == 0) {
values_list[[i]]
next
}
# 提取值并转换为矩阵
window_vals
# 创建全尺寸5x5矩阵（自动填充NA）
full_matrix
# 计算填充位置（相对窗口的偏移量，索引从1开始）
row_offset # 原始偏移量
col_offset
# 转换为矩阵有效索引（防止越界）
row_index # 调整索引从1开始
col_index
# 有效性校验
valid_pos = 1) & (row_index <= 5) &
(col_index >= 1) & (col_index <= 5)
# 执行矩阵赋值（仅有效位置）
full_matrix[row_index[valid_pos], col_index[valid_pos]]
# 按行展开并保存
values_list[[i]]
}
# 创建最终数据框
coords_rast_df
ID = coords$ID, # 假设原始数据有id列
X = coords[, "X"], # 经度坐标
Y = coords[, "Y"] # 纬度坐标
)
# 添加25个值列
value_cols
colnames(value_cols) "V", 1:25)
coords_rast_df
# 查看结果
head(coords_rast_df)
更多内容欢迎加入地理数据科学培训班：
