构建瑞典数据10×10km网格单元及生成Grid cell ID

问题需求 #

我有一份包含纬度、经度及计数信息的Dataframe，需要为瑞典区域数据创建10×10km网格单元，生成Grid cell ID并统计每个网格内的采样点数量，但编写的R代码无法实现需求，寻求正确方案。

原始数据（Sweden_MM4） #

Taxon_id,validScientificName,individualCount,longitude,latitude,coordinatePrecision,County,Municipality,County2_Parrish
200816,Volucella pellucens,1,17.4133208,62.6939581,100,Västernorrland Timrå,Medelpad,Ljustorp
100395,Arctophila superbiens,6,13.9191707,55.6894566,50,Skåne Sjöbo,Skåne,Lövestad
100395,Arctophila superbiens,5,13.9476204,55.6882818,50,Skåne Tomelilla,Skåne,Andrarum
100395,Arctophila superbiens,8,13.9478629,55.6885198,50,Skåne Tomelilla,Skåne,Andrarum
100395,Arctophila superbiens,1,13.9480606,55.6887375,50,Skåne Tomelilla,Skåne,Andrarum
200816,Volucella pellucens,2,16.2606655,58.2811607,241,Östergötland Åtvidaberg,Östergötland,Yxnerum
200814,Volucella bombylans,1,18.6925487,60.1198292,10,Stockholm Norrtälje,Uppland,Häverö

原错误代码 #

library(sf)
# Convert Sweden_MM4  to an sf object
my_data_sf <- st_as_sf(Sweden_MM4, coords = c("longitude", "latitude"), crs = 4326)
my_data_sf <- st_transform(my_data_sf, crs = "+proj=utm +zone=33 +datum=WGS84")  # Set the CRS to an appropriate projection that uses meters 
# Create a grid of 10x10 km using st_make_grid
grid_size_km <- 10
grid_size_m <- grid_size_km * 1000
grid_polygons <- st_make_grid(my_data_sf, cellsize = grid_size_m, what = "polygons")
# Set the CRS of the grid to match your data
grid_polygons <- st_set_crs(grid_polygons, st_crs(my_data_sf))
# Perform a spatial overlay to associate each data point with the corresponding grid cell ID
grid_overlay <- st_intersection(grid_polygons, my_data_sf)
str(my_data_sf)

# Plot the data points
data_plot <- ggplot(Sweden_MM4, aes(x = longitude, y = latitude)) +
  geom_point()
# Plot the grid overlay
grid_plot <- ggplot(grid_overlay) +
  geom_sf(fill = "transparent", color = "blue", size = 0.5)
# Combine the data plot and the grid plot
# Combine the data plot and the grid plot side by side
combined_plot <- plot_grid(data_plot, grid_plot, labels = c("Data Points", "Grid Overlay"), ncol = 2)

# Display the combined plot
print(combined_plot)

预期输出 #

Taxon_id,validScientificName,individualCount,longitude,latitude,coordinatePrecision,Gridcell_REFID,County,Municipality,County2_Parrish
200816,Volucella pellucens,1,17.4133208,62.6939581,100,,Västernorrland Timrå,Medelpad,Ljustorp
100395,Arctophila superbiens,6,13.9191707,55.6894566,50,,Skåne Sjöbo,Skåne,Lövestad
100395,Arctophila superbiens,5,13.9476204,55.6882818,50,,Skåne Tomelilla,Skåne,Andrarum
100395,Arctophila superbiens,8,13.9478629,55.6885198,50,,Skåne Tomelilla,Skåne,Andrarum
100395,Arctophila superbiens,1,13.9480606,55.6887375,50,,Skåne Tomelilla,Skåne,Andrarum
200816,Volucella pellucens,2,16.2606655,58.2811607,241,,Östergötland Åtvidaberg,Östergötland,Yxnerum
200814,Volucella bombylans,1,18.6925487,60.1198292,10,,Stockholm Norrtälje,Uppland,Häverö

正确解决方案 #

原代码的问题在于：st_intersection仅实现了网格与点的空间切割，未生成网格ID并关联回原始数据，也未完成网格内采样点的统计。以下是修正后的完整代码：

library(sf)
library(dplyr)
library(ggplot2)
library(cowplot)

# 加载原始数据（如果未提前加载）
Sweden_MM4 <- read.csv(text = "Taxon_id,validScientificName,individualCount,longitude,latitude,coordinatePrecision,County,Municipality,County2_Parrish
200816,Volucella pellucens,1,17.4133208,62.6939581,100,Västernorrland Timrå,Medelpad,Ljustorp
100395,Arctophila superbiens,6,13.9191707,55.6894566,50,Skåne Sjöbo,Skåne,Lövestad
100395,Arctophila superbiens,5,13.9476204,55.6882818,50,Skåne Tomelilla,Skåne,Andrarum
100395,Arctophila superbiens,8,13.9478629,55.6885198,50,Skåne Tomelilla,Skåne,Andrarum
100395,Arctophila superbiens,1,13.9480606,55.6887375,50,Skåne Tomelilla,Skåne,Andrarum
200816,Volucella pellucens,2,16.2606655,58.2811607,241,Östergötland Åtvidaberg,Östergötland,Yxnerum
200814,Volucella bombylans,1,18.6925487,60.1198292,10,Stockholm Norrtälje,Uppland,Häverö")

# 转换为空间对象并投影到米制坐标系
my_data_sf <- st_as_sf(Sweden_MM4, coords = c("longitude", "latitude"), crs = 4326)
# 瑞典跨UTM 32N和33N，这里用33N覆盖大部分区域，全区域可改用SWEREF99投影
my_data_sf <- st_transform(my_data_sf, crs = "+proj=utm +zone=33 +datum=WGS84")

# 创建10×10km网格并添加唯一ID
grid_size_km <- 10
grid_size_m <- grid_size_km * 1000
grid_polygons <- st_make_grid(my_data_sf, cellsize = grid_size_m, what = "polygons") %>%
  st_sf() %>%
  mutate(Gridcell_REFID = row_number())

# 将采样点与所属网格关联
points_with_grid <- st_join(my_data_sf, grid_polygons, join = st_within)

# 统计每个网格的采样点数量和个体总数
grid_stats <- points_with_grid %>%
  as.data.frame() %>%
  group_by(Gridcell_REFID) %>%
  summarize(
    sample_count = n(),
    total_individuals = sum(individualCount)
  ) %>%
  ungroup()

# 将网格ID合并回原始数据框，生成符合预期的输出
final_data <- Sweden_MM4 %>%
  bind_cols(Gridcell_REFID = points_with_grid$Gridcell_REFID)

# 查看结果
print(final_data)
print(grid_stats)

# 可视化展示
# 采样点及所属网格ID标注
data_plot <- ggplot(final_data, aes(x = longitude, y = latitude, label = Gridcell_REFID)) +
  geom_point(size = 3) +
  geom_text(hjust = 1.2, vjust = 0, color = "red") +
  labs(title = "采样点及所属网格ID")

# 网格覆盖与采样点分布
grid_plot <- ggplot() +
  geom_sf(data = grid_polygons, fill = "transparent", color = "blue", size = 0.5) +
  geom_sf_text(data = grid_polygons, aes(label = Gridcell_REFID), color = "blue") +
  geom_sf(data = my_data_sf, color = "black", size = 2) +
  labs(title = "10×10km网格覆盖")

# 合并展示两个图
combined_plot <- plot_grid(data_plot, grid_plot, ncol = 2)
print(combined_plot)

关键修正说明 #

1. 网格ID生成：将生成的网格转为sf数据框，添加Gridcell_REFID字段作为网格唯一标识。
2. 空间关联优化：使用st_join结合st_within，高效将每个采样点匹配到所属网格，避免冗余数据。
3. 网格统计：通过分组统计，得到每个网格的采样点数量和个体总数。
4. 原始数据保留：将网格ID合并回原始Dataframe，完全匹配预期输出格式。

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