R中基于网格质心坐标使用st_make_grid自动生成格网

现有5184个规则网格的质心空间点数据，需使用R语言sf包的st_make_grid函数生成恰好包含5184个单元的对应矢量网格，当前手动估算cellsize和offset参数的方式效率低且容易出错，需要实现参数的自动计算。

示例数据结构 #

所用点数据为WGS84坐标系（EPSG:4326），前6行内容如下：

head(df_sf)
Simple feature collection with 6 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -74.29167 ymin: 40.89167 xmax: -73.75 ymax: 41
Geodetic CRS:  WGS 84
  avg.dday30                   geometry
1  30.439710       POINT (-74.29167 41)
2  29.280088       POINT (-74.21667 41)
3 190.002556 POINT (-74.06667 40.99167)
4   9.808251 POINT (-73.76667 40.89167)
5   8.892234 POINT (-73.75833 40.89167)

原有手动参数代码 #

原有代码通过手动估算行列数71、手动调整偏移量生成网格，无法复用：

df_sf <- st_as_sf(df, coords=c("lon","lat"), crs=4326)
grid <- st_make_grid(df_sf,
                     cellsize=c(diff(st_bbox(df_sf)[c(1,3)]),diff(st_bbox(df_sf)[c(2,4)]))/c(71,71),
                     offset = st_bbox(df_sf)[c("xmin", "ymin")]-c(cellsize/2,cellsize/2),
                     what="polygons",
                     square=TRUE
)

规则网格的质心坐标具备等间距分布特征，可直接从质心坐标提取所需参数，无需手动估算：

• 提取所有质心的x、y坐标，去重排序后计算相邻坐标的间距，得到x、y方向的单元格尺寸
• 以最小质心坐标减去半个单元格边长作为网格左下角偏移量，保证质心恰好落在格网中心
• 以去重后的x、y坐标数量作为网格列数、行数，直接传入st_make_grid避免边界扩展导致格网数量不符

完整可运行代码：

library(sf)

# 提取质心坐标矩阵
coords <- st_coordinates(df_sf)
# 坐标保留7位小数消除浮点误差，分别提取排序后的唯一x、y值
x_vals <- sort(unique(round(coords[, "X"], 7)))
y_vals <- sort(unique(round(coords[, "Y"], 7)))

# 自动计算x、y方向单元格大小，取相邻坐标差的中位数过滤异常值
cellsize_x <- median(diff(x_vals))
cellsize_y <- median(diff(y_vals))
cellsize <- c(cellsize_x, cellsize_y)

# 自动计算网格左下角偏移量
offset <- c(
  min(x_vals) - cellsize_x / 2,
  min(y_vals) - cellsize_y / 2
)

# 自动获取网格行列数
grid_n <- c(length(x_vals), length(y_vals))
# 校验质心总数与格网总数匹配
stopifnot(grid_n[1] * grid_n[2] == nrow(df_sf))

# 生成目标网格
grid <- st_make_grid(
  x = df_sf,
  cellsize = cellsize,
  offset = offset,
  n = grid_n,
  what = "polygons",
  square = TRUE
)

# 校验最终格网数量为5184
stopifnot(length(grid) == 5184)

注意事项 #

• 坐标保留小数位数可根据自身数据精度调整，核心是消除浮点计算误差导致的唯一值识别错误
• 若质心为完全规则无误差的网格，可直接取diff(x_vals)的唯一值作为cellsize，无需用中位数
• 传入n参数可强制st_make_grid生成指定行列数的网格，避免函数自动扩展边界导致格网数量超出预期

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