使用R基于质心坐标列表创建空间方格网格的实现方法咨询

你提供的点数据为WGS84坐标系下的规则格网点，经纬度间隔为0.5度，可按照以下步骤生成符合要求的网格： #

完整实现代码 #

# 加载依赖包
library(sf)
library(dplyr)

# --------------------------
# 步骤1：准备点数据（如果你已有sf格式的点对象可跳过本段）
# 示例数据构造，替换为你自己的全量数据读取逻辑即可
points_sf <- data.frame(
  cell = 1:6,
  lon = c(19.75, 20.25, 18.25, 18.75, 19.25, 19.75),
  lat = c(-34.75, -34.75, -34.25, -34.25, -34.25, -34.25)
) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) # 4326为WGS84对应的EPSG编码

# --------------------------
# 步骤2：计算网格参数
# 提取点坐标
coords <- st_coordinates(points_sf)
# 自动计算x、y方向的单元格大小，无需手动输入
x_step <- diff(sort(unique(coords[,1])))[1]
y_step <- diff(sort(unique(coords[,2])))[1]
cellsize <- c(x_step, y_step)
# 计算网格偏移量：默认st_make_grid从bbox左下角起生成网格，偏移半个单元格即可让点落在网格中心
offset <- c(min(coords[,1]) - x_step/2, min(coords[,2]) - y_step/2)

# --------------------------
# 步骤3：生成并匹配网格
# 生成全量正方形网格
grid_sf <- st_make_grid(
  x = points_sf,
  cellsize = cellsize,
  offset = offset,
  what = "polygons"
) %>% st_sf()
# 将原始点的属性（如cell编号）匹配到对应网格
grid_sf <- st_join(grid_sf, points_sf)
# 过滤无对应点的空网格（如果你的点没有覆盖bbox内所有格网位置可保留本步）
grid_sf <- filter(grid_sf, !is.na(cell))

# --------------------------
# 验证：检查网格中心是否与原始点完全匹配
grid_centers <- st_centroid(grid_sf)
all.equal(st_geometry(grid_centers), st_geometry(points_sf[match(grid_sf$cell, points_sf$cell), ]))
# 返回TRUE即为匹配正确

效果参考 #

#

注意事项 #

• 本方案默认输入点为x、y方向间隔固定的规则格网点，与你提供的示例数据特征一致
• 若全量非洲数据运算时出现经纬度面积形变问题，可先将点投影至非洲专用等面积投影（如EPSG:102022）生成网格后，再转回WGS84坐标系
• 生成的网格为无重叠、无缝隙的相邻正方形，符合st_make_grid的输出样式

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