基于不规则四边形四角在ggmap中生成网格并提取质心绘图

核心思路 #

行列均匀但单元格大小不规则的网格，本质是在网格四条边之间按比例线性插值：先沿南北方向生成11条东西向样线的端点，再沿每条样线的东西方向插值得到11个站点；之后通过相邻站点坐标的平均值计算每个10×10单元格的中心。

代码实现示例 #

假设你有单个网格的四个角GPS坐标（替换为你实际的经纬度），以下是完整的R代码流程：

library(ggmap)
library(dplyr)

# 1. 定义单个网格的四个角坐标（示例数据，替换为你的实际坐标）
grid_corners <- list(
  nw = c(lon = -122.4194, lat = 37.7749), # 西北角
  ne = c(lon = -122.3900, lat = 37.7749), # 东北角
  sw = c(lon = -122.4194, lat = 37.7600), # 西南角
  se = c(lon = -122.3900, lat = 37.7600)  # 东南角
)

# 2. 生成11×11的站点坐标（网格的所有角点）
# 生成行/列的比例序列（0到1，共11个点）
p_seq <- seq(0, 1, length.out = 11)
q_seq <- seq(0, 1, length.out = 11)

# 生成所有站点
grid_stations <- expand.grid(p = p_seq, q = q_seq) %>%
  rowwise() %>%
  mutate(
    # 计算当前行的左右端点坐标
    lon_left = grid_corners$nw["lon"]*(1-p) + grid_corners$sw["lon"]*p,
    lat_left = grid_corners$nw["lat"]*(1-p) + grid_corners$sw["lat"]*p,
    lon_right = grid_corners$ne["lon"]*(1-p) + grid_corners$se["lon"]*p,
    lat_right = grid_corners$ne["lat"]*(1-p) + grid_corners$se["lat"]*p,
    # 沿当前样线插值得到站点坐标
    lon = lon_left*(1-q) + lon_right*q,
    lat = lat_left*(1-q) + lat_right*q,
    # 标记样线和站点编号（方便定位特定点）
    line_num = which(p_seq == p),
    station_num = which(q_seq == q)
  ) %>%
  ungroup() %>%
  select(line_num, station_num, lon, lat)

# 3. 提取10×10单元格的中心坐标
grid_centers <- grid_stations %>%
  filter(line_num < 11, station_num < 11) %>%
  rowwise() %>%
  mutate(
    # 取当前站点与右下相邻站点的坐标平均值作为中心
    center_lon = (lon + grid_stations$lon[grid_stations$line_num == line_num+1 & grid_stations$station_num == station_num+1])/2,
    center_lat = (lat + grid_stations$lat[grid_stations$line_num == line_num+1 & grid_stations$station_num == station_num+1])/2
  ) %>%
  ungroup() %>%
  select(line_num, station_num, center_lon, center_lat)

# 4. 绘制地图并添加站点/中心
# 获取底图（替换为你需要的区域）
map <- get_map(location = c(lon = mean(c(grid_corners$nw["lon"], grid_corners$se["lon"])), 
                            lat = mean(c(grid_corners$nw["lat"], grid_corners$se["lat"]))),
               zoom = 14, maptype = "roadmap")

# 绘制地图 + 所有站点 + 特定中心/站点
ggmap(map) +
  geom_point(data = grid_stations, aes(x = lon, y = lat), color = "blue", size = 1) +
  # 示例：绘制第5条样线的第7个站点
  geom_point(data = filter(grid_stations, line_num == 5, station_num ==7), 
             aes(x = lon, y = lat), color = "red", size = 3) +
  # 示例：绘制第3-7行、第2-6列的单元格中心
  geom_point(data = filter(grid_centers, line_num %in% 3:7, station_num %in%2:6), 
             aes(x = center_lon, y = center_lat), color = "green", size = 2, shape = 17)

多网格处理 #

如果你有三个独立网格，只需将每个网格的四角坐标放入列表，循环执行上述流程即可：

# 示例：三个网格的四角坐标列表
multi_grids <- list(
  grid1 = list(nw = c(lon = -122.4194, lat = 37.7749), ne = c(lon = -122.3900, lat = 37.7749),
               sw = c(lon = -122.4194, lat = 37.7600), se = c(lon = -122.3900, lat = 37.7600)),
  grid2 = list(nw = c(lon = -122.3800, lat = 37.7749), ne = c(lon = -122.3500, lat = 37.7749),
               sw = c(lon = -122.3800, lat = 37.7600), se = c(lon = -122.3500, lat = 37.7600)),
  grid3 = list(nw = c(lon = -122.4194, lat = 37.7500), ne = c(lon = -122.3900, lat = 37.7500),
               sw = c(lon = -122.4194, lat = 37.7350), se = c(lon = -122.3900, lat = 37.7350))
)

# 循环处理每个网格
multi_stations <- lapply(names(multi_grids), function(grid_name) {
  corners <- multi_grids[[grid_name]]
  expand.grid(p = p_seq, q = q_seq) %>%
    rowwise() %>%
    mutate(
      lon_left = corners$nw["lon"]*(1-p) + corners$sw["lon"]*p,
      lat_left = corners$nw["lat"]*(1-p) + corners$sw["lat"]*p,
      lon_right = corners$ne["lon"]*(1-p) + corners$se["lon"]*p,
      lat_right = corners$ne["lat"]*(1-p) + corners$se["lat"]*p,
      lon = lon_left*(1-q) + lon_right*q,
      lat = lat_left*(1-q) + lat_right*q,
      line_num = which(p_seq == p),
      station_num = which(q_seq == q),
      grid_id = grid_name
    ) %>%
    ungroup() %>%
    select(grid_id, line_num, station_num, lon, lat)
}) %>% bind_rows()

# 后续绘制同理，只需按grid_id区分即可

关键说明 #

• 定位特定样线站点：直接通过line_num和station_num筛选grid_stations数据框即可，比如filter(grid_stations, line_num == 3, station_num == 8)。
• 单元格中心计算：取当前站点与右下角相邻站点的坐标平均值，确保对应10×10的单元格。
• 坐标适配：如果你的GPS坐标是（纬度，经度），注意调整代码中的lon/lat顺序，确保与ggmap的要求一致（ggmap默认先经度后纬度）。

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