如何在R语言中创建网格单元并计算单元内点数据的平均值？

以下是实现需求的具体步骤，结合raster、dplyr和ggplot2完成：

1. 加载依赖包 #

首先加载需要的R包：

library(raster)
library(dplyr)
library(ggplot2)
library(sf)

2. 将数据转换为空间点对象 #

raster包基于空间数据结构处理，需先把普通数据框转为带坐标系的空间点：

# 定义空间坐标字段
coordinates(df) <- ~longitude + latitude
# 设置通用的WGS84经纬度坐标系
proj4string(df) <- CRS("+proj=longlat +datum=WGS84")

3. 创建自定义网格 #

根据数据的经纬度范围生成指定分辨率的网格：

# 获取数据的经纬度边界范围
data_extent <- extent(df)
# 设置网格分辨率（示例用0.5度单元格，可根据需求调整大小）
grid_raster <- raster(data_extent, res = 0.5)
# 将栅格转为多边形网格对象
grid_polygons <- as(grid_raster, "SpatialPolygonsDataFrame")
# 给每个网格添加唯一ID，方便后续匹配
grid_polygons$grid_id <- seq_len(nrow(grid_polygons))

4. 统计网格内数据并设置阈值 #

匹配点与网格，计算每个网格内的plot数量和目标变量平均值，再过滤不符合阈值的网格：

# 匹配每个数据点所属的网格ID
df$grid_id <- over(df, grid_polygons)$grid_id

# 转换回普通数据框，先计算每个plot的class平均值（因每个plot包含多个观测）
plot_level_avg <- as.data.frame(df) %>%
  group_by(plot) %>%
  summarise(avg_class = mean(class), .groups = "drop")

# 结合网格ID，统计每个网格的plot数量和整体平均class
grid_stats <- as.data.frame(df) %>%
  select(plot, grid_id) %>%
  distinct() %>% # 去重，确保每个plot仅被计数一次
  left_join(plot_level_avg, by = "plot") %>%
  group_by(grid_id) %>%
  summarise(
    plot_count = n(),
    mean_class = mean(avg_class),
    .groups = "drop"
  )

# 设置最低plot数量阈值，示例设为1，可自行修改
min_plot_threshold <- 1
filtered_grid_stats <- grid_stats %>% filter(plot_count >= min_plot_threshold)

# 将统计结果合并到网格多边形对象中
grid_with_stats <- merge(grid_polygons, filtered_grid_stats, by = "grid_id")

5. 可视化结果 #

转换为sf对象，用ggplot绘制带统计值的网格地图：

# 转为sf格式适配ggplot绘图
grid_sf <- st_as_sf(grid_with_stats)

ggplot() +
  geom_sf(data = grid_sf, aes(fill = mean_class), color = "white", size = 0.2) +
  scale_fill_viridis_c(option = "plasma") +
  labs(title = "网格单元内plot的class平均值",
       fill = "平均class值",
       subtitle = paste("仅显示plot数量≥", min_plot_threshold, "的网格")) +
  theme_minimal()

注意事项 #

• 网格分辨率res参数可根据数据密度调整，数值越小网格越精细
• 如果无需先计算每个plot的平均值，直接对所有观测统计，可跳过plot_level_avg步骤，直接按网格ID汇总所有class值
• 坐标系可根据数据实际情况修改，比如使用UTM投影坐标系

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