在R中基于坐标创建10×10km栅格网格并关联数据框

步骤1：创建10×10公里栅格网格数据框 #

由于经纬度是球面坐标，无法直接生成等距的10公里网格，我们需要先将研究区域转换为平面投影（比如UTM），再创建网格并分配编号，最后转回经纬度坐标系。

# 加载所需工具包
library(sf)
library(terra)

# 自定义研究区域的经纬度范围（请根据实际需求修改）
lon_min <- 67
lon_max <- 68
lat_min <- 67
lat_max <- 68

# 创建研究区域边界（采用WGS84经纬度坐标系，EPSG:4326）
region_bbox <- st_bbox(c(xmin = lon_min, xmax = lon_max, ymin = lat_min, ymax = lat_max), 
                       crs = st_crs(4326)) %>%
  st_as_sfc()

# 转换为UTM平面投影（请根据研究区域选择对应EPSG，示例为UTM 36N，EPSG:32636）
utm_crs <- st_crs(32636)
region_utm <- st_transform(region_bbox, utm_crs)

# 生成10×10公里网格（cellsize单位为米，10km=10000m）
grid_utm <- st_make_grid(region_utm, cellsize = c(10000, 10000), what = "polygons") %>%
  st_sf() %>%
  mutate(Grid = row_number())  # 为每个网格分配唯一编号

# 转回WGS84经纬度坐标系，得到最终网格数据框df
df <- st_transform(grid_utm, st_crs(4326))

步骤2：为现有数据框分配网格编号 #

将带经纬度的现有数据转换为空间点对象，通过空间匹配找到每个点所属的网格，提取编号后合并回原数据。

# 示例现有数据框（替换为你的实际数据）
existing_df <- data.frame(
  LONGITUDE = c(67.71150, 67.72150),
  LATITUDE = c(67.71150, 67.72150),
  Year = c(2017, 2017),
  Presence = c("Y", "N"),
  Set = c(1, 2)
)

# 将现有数据转换为空间点对象（WGS84坐标系）
existing_sf <- st_as_sf(existing_df, 
                        coords = c("LONGITUDE", "LATITUDE"), 
                        crs = st_crs(4326))

# 空间连接：匹配每个点所在的网格，获取Grid编号
joined_data <- st_join(existing_sf, df, join = st_within)

# 转换为普通数据框，移除几何列（不需要空间信息时可执行）
final_df <- as.data.frame(joined_data) %>%
  dplyr::select(-geometry)

# 查看结果
print(final_df)

关键注意事项 #

• 投影选择：必须使用平面投影（如UTM）才能保证网格是严格的10×10公里，经纬度坐标系无法实现等距网格。请根据研究区域查询对应的EPSG代码。
• 范围调整：务必将lon_min、lon_max、lat_min、lat_max修改为你的实际研究区域坐标。
• 匹配规则：st_within确保点完全落在网格内，若需要包含边界点的匹配，可调整连接方式为st_intersects。

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