R语言如何实现经纬度数据按0.5x0.5度网格分组

解决方案 #

方法1：数学计算法（最高效，无额外依赖） #

直接通过数值计算得到每个点所属0.5°网格的左下角经纬度，同时生成唯一网格ID，逻辑简单且适配十万级以上大数据量：

library(tidyverse)

# 此处用你提供的样例数据做演示，实际使用时替换为你的dassem数据框即可
dassem <- tibble(
  Latitude = c(31.602, 31.675, 31.747, 32.152, 32.218, 32.285, 32.348, 32.412, 32.475, 32.535, 32.595, 32.677, 32.763),
  Longitude = c(-39.848, -39.467, -39.083, -36.795, -36.408, -36.022, -35.635, -35.247, -34.858, -34.47, -34.082, -33.707, -33.323)
)

pos <- dassem %>%
  mutate(
    # 计算所属0.5°网格的左下角经纬度
    grid_lon = floor(Longitude / 0.5) * 0.5,
    grid_lat = floor(Latitude / 0.5) * 0.5,
    # 生成唯一网格ID，可直接用于分组统计
    grid_id = paste0(grid_lon, "_", grid_lat)
  )

运行后每个点会新增3个字段：grid_lon为网格左边界经度、grid_lat为网格下边界纬度、grid_id为唯一网格标识，三个字段任意一个都可以作为分组依据。如果需要取网格中心点坐标，直接用grid_lon + 0.25和grid_lat + 0.25计算即可。

方法2：cut函数法（适配非对齐自定义网格） #

如果你的网格不需要和全球0.5°标准网格对齐，只想按你的数据经纬度范围生成网格，可以用cut函数，同时修复了你原有代码的边界异常问题：

library(tidyverse)

pos <- dassem %>%
  mutate(
    # 按0.5步长切割经度，include.lowest=T避免边界值被判定为NA
    grid_lon_cut = cut(Longitude, 
                       breaks = seq(floor(min(Longitude)/0.5)*0.5, ceiling(max(Longitude)/0.5)*0.5, by = 0.5),
                       include.lowest = TRUE),
    grid_lat_cut = cut(Latitude,
                       breaks = seq(floor(min(Latitude)/0.5)*0.5, ceiling(max(Latitude)/0.5)*0.5, by = 0.5),
                       include.lowest = TRUE),
    # 生成数字格式的网格ID
    grid_id = as.integer(grid_lon_cut) * 1000 + as.integer(grid_lat_cut)
  )

原有代码问题说明 #

1. 生成grid序列时没有向上扩展边界，刚好等于经纬度最大值的点会找不到符合gridx>x的结果，返回空值导致异常
2. 网格ID的计算逻辑未考虑ycell起始值问题，容易出现ID重复或者跨度过大的问题

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