R语言Shapefile相交计算异常：结果不符预期的解决求助

问题根源 #

• 仅返回5个FSA：未正确筛选目标省份（PRUID）的FSA数据，或Shapefile中FSA的PRUID字段映射错误
• 规整重复的百分比：未基于实际相交面积计算占比，而是误用了等分逻辑，或未以ADA的总面积作为占比计算的分母

正确实现代码 #

以下是修正后的calculate_intersection_percentages函数，基于sf和dplyr实现：

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

calculate_intersection_percentages <- function(fsa_shp_path, ada_shp_path, target_pruid) {
  # 1. 读取Shapefile并过滤目标省份数据，统一投影坐标系
  fsa <- st_read(fsa_shp_path, quiet = TRUE) %>%
    filter(PRUID == target_pruid) %>%
    st_transform(3347)  # 转换为加拿大Lambert投影，避免地理坐标系面积计算误差
  
  ada <- st_read(ada_shp_path, quiet = TRUE) %>%
    filter(PRUID == target_pruid) %>%
    st_transform(3347)
  
  # 2. 计算每个ADA的总面积
  ada_areas <- ada %>%
    mutate(ada_total_area = st_area(.)) %>%
    select(ADAUID, ada_total_area)
  
  # 3. 计算所有ADA与FSA的相交区域及面积
  intersections <- st_intersection(ada, fsa) %>%
    mutate(intersect_area = st_area(.)) %>%
    select(ADAUID, FSAUID, intersect_area) %>%
    st_drop_geometry()  # 移除几何属性，仅保留数值数据
  
  # 4. 合并总面积数据，计算相交占比
  intersection_percentages <- intersections %>%
    left_join(ada_areas, by = "ADAUID") %>%
    mutate(percent = as.numeric(intersect_area / ada_total_area) * 100) %>%
    select(ADAUID, FSAUID, percent)
  
  # 5. 转换为矩阵格式（行=ADA，列=FSA），填充无相交的0值
  percentage_matrix <- intersection_percentages %>%
    pivot_wider(names_from = FSAUID, values_from = percent, values_fill = 0) %>%
    column_to_rownames("ADAUID") %>%
    as.matrix()
  
  return(percentage_matrix)
}

关键说明 #

• 统一投影：使用EPSG:3347（加拿大Lambert投影）替代默认地理坐标系，彻底避免经纬度计算面积的系统误差
• 精准筛选：同时对FSA和ADA图层过滤指定target_pruid的记录，解决FSA数量缺失问题
• 真实面积占比：以ADA的实际总面积为分母，相交区域面积为分子计算占比，杜绝规整的错误百分比
• 完整矩阵：用values_fill = 0填充无相交的ADA-FSA组合，保证矩阵维度匹配实际数据量

使用示例 #

# 替换为你的本地Shapefile路径
result_matrix <- calculate_intersection_percentages(
  fsa_shp_path = "./fsa_shapefile.shp",
  ada_shp_path = "./ada_shapefile.shp",
  target_pruid = "10"  # 纽芬兰PRUID
)

# 验证FSA数量是否正确（应返回35）
ncol(result_matrix)

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