R语言中如何为多边形提取的栅格值补充对应x、y坐标？

和点要素直接读取矢量点自带坐标的逻辑不同，多边形提取栅格的核心是匹配所有与多边形范围相交的栅格像元，因此补全坐标的本质是获取每个被命中的栅格像元的中心位置信息，有两种成熟的实现路径可适配原有代码逻辑：

基于现有raster包的适配修改 #

• 修正原有代码的命名空间写法小问题：原代码sf:st_read缺少一个冒号，会导致函数找不到，正确写法为sf::st_read
• 调整extract函数的参数，开启像元编号返回开关：把cellnumbers参数设为TRUE，提取结果会新增一列记录每个匹配到的栅格对应的像元编号

  shp <- sf::st_read('example.shp')
  r  <- raster::raster('example.tif')
  # 同时返回多边形ID、栅格值、匹配的栅格像元编号
  extract_res <- raster::extract(r, shp, df=TRUE, cellnumbers=TRUE)
  

• 基于像元编号查询对应坐标：使用xyFromCell()函数输入栅格对象和像元编号，直接返回每个像元的中心x、y坐标，再合并到提取结果数据框即可

  # 获取每个像元的x、y坐标
  cell_coords <- raster::xyFromCell(r, extract_res$cell)
  # 合并坐标到结果表
  final_extract_df <- cbind(extract_res, cell_coords)
  

注意：上述方法返回的结果是「每个多边形范围内的所有相交栅格像元」的明细记录。如果你的需求是先对每个多边形内的栅格做聚合统计（比如求多边形内栅格的平均值、总和），再给每条聚合结果配对应坐标，只需要两步调整：

1. 在提取函数中新增fun参数指定聚合规则，比如fun=mean代表求每个多边形内栅格值的平均
2. 不再读取栅格像元的坐标，而是通过sf::st_centroid()计算每个多边形的质心坐标，再合并到聚合结果中即可

更高效的terra包实现方案 #

raster包目前已停止核心性能迭代，官方推荐使用功能对等、处理效率更高的terra包完成同类操作，代码逻辑更简洁：

library(terra)
# 读取矢量和栅格数据
shp <- terra::vect('example.shp')
r  <- terra::rast('example.tif')
# 提取时直接开启坐标返回参数，结果自动包含每个像元的x、y
extract_res <- terra::extract(r, shp, cells=TRUE, xy=TRUE)

该方案不需要额外做坐标匹配，函数内置的xy=TRUE参数会自动为每一条提取记录附上对应栅格像元的中心坐标，处理大尺寸栅格时的速度通常比raster包快3-5倍。 #

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