如何在R中利用经纬度创建7×7平方公里的空间网格？

核心问题分析 #

你的代码存在三个关键问题：

• 输入的DF是普通数据框，不是带空间属性的sf对象，st_make_grid无法识别其空间信息
• 用度作为cellsize单位，经纬度的度在不同纬度对应的实际距离差异极大，无法精准对应7平方公里的网格
• 默认st_make_grid是生成覆盖所有输入点范围的网格，而非以每个点为质心的独立网格

分步解决方案 #

1. 转换数据为sf空间对象 #

首先把普通数据框转换成带WGS84坐标系的sf点对象：

library(sf)

# 转换为sf点对象，注意先传经度(lon)再传纬度(lat)
DF_sf <- st_as_sf(DF, coords = c("lon", "lat"), crs = 4326)

2. 切换到米单位的投影坐标系 #

要创建固定面积的网格，必须使用以米为单位的投影（比如UTM），先自动匹配数据对应的UTM带：

# 计算数据中心经度，确定UTM带
center_lon <- mean(st_coordinates(DF_sf)[,1])
utm_zone <- floor((center_lon + 180)/6) + 1
# 构建北半球UTM的CRS编码（你的数据在北纬43度）
utm_crs <- paste0("EPSG:326", utm_zone)
# 转换投影到UTM
DF_utm <- st_transform(DF_sf, crs = utm_crs)

3. 计算7平方公里网格的边长 #

7平方公里的正方形，边长为sqrt(7)*1000米（约2645.8米）：

cell_side <- sqrt(7) * 1000  # 单位：米

4. 生成以每个点为质心的网格 #

循环每个点，生成对应的正方形网格单元：

grid_list <- list()

for (i in seq_len(nrow(DF_utm))) {
  # 获取当前点的UTM坐标
  point_coords <- st_coordinates(DF_utm[i,])
  # 计算正方形四个顶点的坐标
  x_min <- point_coords[1] - cell_side/2
  x_max <- point_coords[1] + cell_side/2
  y_min <- point_coords[2] - cell_side/2
  y_max <- point_coords[2] + cell_side/2
  # 创建多边形
  square <- st_polygon(list(rbind(
    c(x_min, y_min), c(x_max, y_min),
    c(x_max, y_max), c(x_min, y_max),
    c(x_min, y_min)  # 闭合多边形
  )))
  # 转换为sf对象并添加ID标识
  square_sf <- st_sfc(square, crs = utm_crs)
  grid_list[[i]] <- st_sf(geometry = square_sf, point_id = i)
}

# 合并所有网格单元
MyGrid <- do.call(rbind, grid_list)
# 转换回WGS84经纬度坐标系（方便后续查看）
MyGrid_wgs84 <- st_transform(MyGrid, crs = 4326)

5. 可视化验证 #

# 绘制原始点（红色）
plot(st_geometry(DF_sf), pch = 16, col = "red", main = "7×7km网格（质心为给定经纬度）")
# 叠加网格（蓝色边框）
plot(st_geometry(MyGrid_wgs84), add = TRUE, border = "blue")

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