如何在R中为两个相似但非一致的Shapefile计算Jaccard指数？

问题描述 #

拥有两个高度重叠但结构不同的Shapefile（记为x和y）：y的区域数量更多（例如包含城市辖区，而x仅包含更大的行政区）。核心需求是：为y中的每个区域，找到x里与其交集面积最大的对应区域，然后计算这一对区域的Jaccard指数。

此前尝试的步骤在计算并集时遇到问题：使用st_union(y,x)生成的并集是所有区域的整体合并，无法通过st_join关联到单个y区域对应的配对并集，导致后续无法正确计算分母。

解决方案：利用Jaccard公式绕开并集计算 #

Jaccard指数的公式可转换为：

J(A,B) = |A∩B| / (|A| + |B| - |A∩B|)
其中|A|是y区域的面积，|B|是对应x区域的面积，|A∩B|是两者的交集面积。无需单独计算并集，直接用该公式即可，步骤如下：

1. 计算原始区域的面积 #

先给x和y各自添加面积列：

library(sf)
library(dplyr)

# 替换<y_id_col>和<x_id_col>为数据集的唯一ID列名
x <- x %>% 
  mutate(area_x = st_area(.),
         x_id = <x_id_col>)
y <- y %>% 
  mutate(area_y = st_area(.),
         y_id = <y_id_col>)

2. 计算所有区域对的交集并筛选最大交集 #

生成y与x的所有交集，筛选出每个y区域对应的最大交集记录：

# 计算所有交集，保留双方ID、面积及交集面积
all_intersects <- st_intersection(y, x) %>%
  mutate(intersect_area = st_area(.)) %>%
  select(y_id, x_id, area_y, area_x, intersect_area)

# 筛选每个y区域的最大交集记录
max_intersect_pairs <- all_intersects %>%
  group_by(y_id) %>%
  filter(intersect_area == max(intersect_area)) %>%
  ungroup()

3. 关联回原数据集并计算Jaccard指数 #

将筛选后的结果关联到y，计算最终的Jaccard指数：

y_with_jaccard <- y %>%
  left_join(max_intersect_pairs, by = "y_id") %>%
  # 计算Jaccard指数，处理可能的NA（无交集的区域）
  mutate(jaccard = ifelse(is.na(intersect_area), 0, intersect_area / (area_y + area_x - intersect_area)))

关键说明 #

• 避免使用st_union计算并集：st_union(y,x)会将所有区域合并为单一几何对象，无法对应到单个区域对的并集，这是此前步骤出错的核心原因。
• 单位一致性：st_area返回的面积带有单位（如平方米），R会自动处理单位运算，无需额外转换。
• 无交集处理：若y中某个区域与x无交集，Jaccard指数设为0，可根据需求调整逻辑。

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