在R中为奥斯陆海域物种分布地图添加指定网格的技术求助

现有R代码可绘制奥斯陆海域物种分布地图，所有数据CRS为4326。需在地图上叠加Xcm×Xcm的网格，统计有物种分布数据和无数据区域的占比，但尝试多种添加网格的方法均未达到预期效果，怀疑问题源于WKT格式的geometry列，附上相关代码寻求技术指导。

用户尝试的代码 #

# Filter shp2 to include only the geometry for the county of interest
shp2_oslo <- shp2[shp2$FYLKESNAVN == "Oslo", ]
intersect3_oslo <- intersect3 %>%
  filter(FYLKESNAVN == "Oslo",)


# Plot with labels

oslo_map <- ggplot() +
  geom_sf(data = shp2_oslo, color = "gray30", size = 0.2, fill = NA) +
  geom_sf(data = intersect3_oslo, aes(color = occurrences), size = 1, fill = NA) +
  scale_fill_gradient(low = "green", high = "red") +
  scale_color_gradient(low = "green", high = "red") +  # Adding color scale for points
  ggtitle("Oslo") +
  theme_bw() +
  theme(axis.title.x = element_blank(),  # Remove x-axis title
        axis.title.y = element_blank(),  # Remove y-axis title
        plot.title.position = "plot",  # Position the title centered above the plot
        plot.title = element_text(hjust = 0.5))  # Center the title text

# Here are some things I have tried to add the grid, but none of this / none of the variations on this are working the way I want. This one at least separates the WKT to two columns.

df_sf <- intersect3_oslo %>%
  st_coordinates() %>%
  as.data.frame() %>%
  rename(longitude = X, latitude = Y) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

# Here is some other stuff I tried:

intersect3_oslo1 <- intersect3_oslo %>%
  mutate(
    longitude = as.numeric(str_extract(geometry, "[-+]?[0-9]+\\.?[0-9]*")),
    latitude = as.numeric(str_extract(geometry, "\\s([-+]?[0-9]+\\.?[0-9]*)$"))
  )

df_sf <- intersect3_oslo %>%
  st_coordinates() %>%
  as.data.frame() %>%
  rename(longitude = X, latitude = Y)
print(df_sf)
df_sf <- st_as_sf(df_sf, coords = c("longitude", "latitude"), crs = 4326)

grid <- st_make_grid(df_sf, cellsize = 2000)
grid <- st_as_sf(grid)
grid$grid_id <- seq_along(grid)

# Perform spatial join (count points within each grid cell)
grid_with_counts <- st_join(df_sf, grid, intersect3_oslo) %>%
  group_by(grid_id) %>%
  summarise(occurrences = n()) %>%
  ungroup()

ggplot() +
  geom_sf(data = intersect3_oslo, fill = "lightblue") +
  geom_sf(data = grid_with_counts, aes(fill = species_count), color = "grey") +
  scale_fill_gradient(low = "white", high = "red", name = "Species Count") +
  theme_minimal()

df_sf <- intersect3_oslo %>%
  st_coordinates() %>%
  as.data.frame() %>%
  rename(longitude = X, latitude = Y) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

grid <- st_make_grid(intersect3_oslo, cellsize = 2000)
grid <- st_as_sf(grid)
grid$grid_id <- seq_along(grid)
grid_with_counts <- st_join(grid, df_sf, join = st_within) %>%
  group_by(grid_id) %>%
  summarise(occurrences = n()) %>%
  ungroup()

ggplot() +
  geom_sf(data = intersect3_oslo, fill = "lightblue") +
  geom_sf(data = grid_with_counts, aes(fill = occurrences), color = "grey") +
  scale_fill_gradient(low = "white", high = "red", name = "Occurrences Count") +
  theme_minimal()

# Any guidance would be much appreciated.

解决方案 #

核心问题分析 #

1. CRS单位不匹配：数据用的是WGS84（EPSG:4326，经纬度，单位为度），但设置的cellsize=2000是米，球面坐标系下无法直接用米作为网格单位，必须转换到平面投影（如UTM坐标系）。
2. WKT处理冗余：无需手动提取WKT坐标，sf包可直接将WKT字符串转换为sf对象，避免正则提取的误差。
3. 网格范围不合理：基于点数据创建网格会导致网格仅覆盖点的范围，应基于奥斯陆的边界（shp2_oslo）创建网格，确保覆盖整个研究区域。

完整实现代码 #

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

# 1. 确保数据为sf对象（如果geometry是WKT字符串）
intersect3_oslo_sf <- st_as_sf(intersect3_oslo, wkt = "geometry", crs = 4326)
shp2_oslo_sf <- st_as_sf(shp2_oslo, crs = 4326)

# 2. 转换到平面投影（奥斯陆对应UTM33N，EPSG:32633，单位为米）
utm_crs <- 32633
shp2_oslo_utm <- st_transform(shp2_oslo_sf, crs = utm_crs)
intersect3_oslo_utm <- st_transform(intersect3_oslo_sf, crs = utm_crs)

# 3. 设置网格尺寸（X厘米转换为米：例如X=100cm则为1米）
grid_size_cm <- 100  # 替换为你的目标X值
grid_size_m <- grid_size_cm / 100

# 4. 创建覆盖奥斯陆边界的网格
grid <- st_make_grid(shp2_oslo_utm, cellsize = grid_size_m, square = TRUE)
grid_sf <- st_as_sf(grid) %>%
  mutate(grid_id = row_number())

# 5. 统计每个网格内的物种出现次数，标记有无数据
grid_with_data <- st_join(grid_sf, intersect3_oslo_utm, join = st_contains) %>%
  group_by(grid_id) %>%
  summarise(occurrences = n(),
            has_data = ifelse(n() > 0, TRUE, FALSE)) %>%
  ungroup()

# 6. 计算有无数据区域占比
total_grid <- nrow(grid_with_data)
has_data_ratio <- sum(grid_with_data$has_data) / total_grid * 100
no_data_ratio <- 100 - has_data_ratio
cat(sprintf("有数据区域占比：%.2f%%\n无数据区域占比：%.2f%%\n", has_data_ratio, no_data_ratio))

# 7. 绘制带网格的地图（转换回4326用于可视化）
grid_with_data_4326 <- st_transform(grid_with_data, crs = 4326)
shp2_oslo_4326 <- st_transform(shp2_oslo_utm, crs = 4326)
intersect3_oslo_4326 <- st_transform(intersect3_oslo_utm, crs = 4326)

ggplot() +
  geom_sf(data = grid_with_data_4326, aes(fill = has_data), color = "gray50", size = 0.1) +
  geom_sf(data = intersect3_oslo_4326, aes(color = occurrences), size = 0.8) +
  geom_sf(data = shp2_oslo_4326, color = "black", size = 0.3, fill = NA) +
  scale_fill_manual(values = c("FALSE" = "white", "TRUE" = "#ff9999"), 
                    name = "有无数据", labels = c("无", "有")) +
  scale_color_gradient(low = "green", high = "red", name = "出现次数") +
  ggtitle("奥斯陆海域物种分布与网格统计") +
  theme_bw() +
  theme(axis.title = element_blank(),
        plot.title = element_text(hjust = 0.5))

关键说明 #

• 投影转换：UTM坐标系是平面投影，单位为米，能准确创建指定厘米/米尺寸的网格。奥斯陆位于UTM第33N带，对应EPSG:32633。
• 网格创建：基于奥斯陆边界创建网格，确保覆盖整个研究区域，避免遗漏无数据区域。
• 数据统计：通过st_join结合st_contains统计每个网格内的点数，直接标记有无数据，后续可快速计算占比。

内容的提问来源于stack exchange，提问作者gn-gia