R语言：构建多边形覆盖百分比矩阵及优化实现问询

我有两个加拿大统计局的shapefile文件（已导入R，转换CRS至WGS84，修复了无效几何）：

• file_1：包含1507个MULTIPOLYGON，字段含ADAUID、LANDAREA等
• file_2：包含310个MULTIPOLYGON，字段含PCUID、PCNAME等

需求是构建两个百分比矩阵：

1. 每个file_1多边形被file_2多边形覆盖的百分比矩阵（行对应file_1，列对应file_2）
2. 每个file_2多边形被file_1多边形覆盖的百分比矩阵（行对应file_2，列对应file_1）

当前我用双重循环计算两两多边形的交集面积占比，代码正在运行，但不确定实现是否正确，也想知道更高效的优化方法。

一、当前实现的正确性验证 #

你的代码逻辑方向没问题，但存在几个细节问题：

1. 资源浪费：先执行了st_intersection(file_1, file_2)获取所有交集，但后续循环又重复计算两两交集，完全没用到已有的结果。
2. 空交集判断：当两个多边形无交集时，st_intersection返回空sf对象，st_area会返回长度为0的对象，你用length(intersection_area) > 0可以处理，但用!sf::st_is_empty(intersection_area)判断更准确。
3. 需求未完全覆盖：当前代码只计算了file_1被file_2覆盖的百分比，未实现file_2被file_1覆盖的部分。

二、高效优化方案 #

双重循环的时间复杂度为O(n*m)（约467万次计算），改用sf矢量运算结合dplyr/tidyr可以大幅提速，且代码更简洁：

步骤1：预计算多边形面积 #

先计算每个多边形的面积并转为数值类型（避免单位运算警告）：

library(sf)
library(dplyr)
library(tidyr)

# 为file_1添加面积字段
file_1 <- file_1 %>% 
  mutate(area_file1 = as.numeric(st_area(geometry)))

# 为file_2添加面积字段
file_2 <- file_2 %>% 
  mutate(area_file2 = as.numeric(st_area(geometry)))

步骤2：批量计算交集与占比 #

用st_intersection一次性获取所有有交集的多边形对，同时计算两种覆盖百分比：

# 获取所有交集对并计算占比
intersect_pairs <- st_intersection(file_1, file_2) %>%
  mutate(
    intersect_area = as.numeric(st_area(geometry)),
    # file_1被file_2覆盖的百分比
    coverage_file1 = 100 * intersect_area / area_file1,
    # file_2被file_1覆盖的百分比
    coverage_file2 = 100 * intersect_area / area_file2
  ) %>%
  st_drop_geometry() %>%  # 移除几何字段，保留关键数据
  select(ADAUID, PCUID, coverage_file1, coverage_file2)

步骤3：转换为百分比矩阵 #

用pivot_wider将长格式数据转为矩阵，无交集的位置填充0：

# 构建file_1被file_2覆盖的矩阵（行=file_1的ADAUID，列=file_2的PCUID）
matrix_file1 <- intersect_pairs %>%
  pivot_wider(
    id_cols = ADAUID,
    names_from = PCUID,
    values_from = coverage_file1,
    values_fill = 0
  ) %>%
  column_to_rownames("ADAUID") %>%
  as.matrix()

# 构建file_2被file_1覆盖的矩阵（行=file_2的PCUID，列=file_1的ADAUID）
matrix_file2 <- intersect_pairs %>%
  pivot_wider(
    id_cols = PCUID,
    names_from = ADAUID,
    values_from = coverage_file2,
    values_fill = 0
  ) %>%
  column_to_rownames("PCUID") %>%
  as.matrix()

优化优势 #

• 空间索引优化：st_intersection内部会自动使用空间索引，只计算实际有交集的多边形对，避免循环中大量无意义的空交集计算。
• 矢量运算效率：R的矢量运算比显式循环快数倍，尤其适合大规模多边形数据。
• 代码可读性：逻辑清晰，避免嵌套循环的冗余代码，便于维护和修改。

补充说明 #

1. 用ADAUID和PCUID作为矩阵行/列名，比序号更具业务意义，方便后续分析对应区域。
2. 如果需要保留原始行号作为矩阵名，只需在预计算面积时添加row_id = row_number()字段，后续用该字段替代ADAUID/PCUID即可。

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