R计算点缓冲区与重叠多边形非重叠面积及intersect空结果解决方案

问题根源定位 #

• 坐标赋值顺序错误：创建空间点对象时，将纬度（lat）设为坐标系x轴、经度（long）设为y轴，与WGS84坐标系默认的（经度=x、纬度=y）规则相悖，经UTM投影转换后，生成的圆形缓冲区与贫民窟shapefile分属完全无关的空间位置，自然无法检测到交集。
• 多边形合并逻辑存在风险：逐次使用gUnion合并多个贫民窟shapefile时，未统一所有输入图层的投影基准，易产生拓扑匹配错误；且逐次合并的处理方式在多边形存在重叠时，易生成无效的冗余几何元素。

实现流程修正 #

• 严格遵循坐标系的坐标顺序规范创建空间点对象，确保WGS84下的坐标值与轴定义匹配。
• 先将所有原始贫民窟shapefile统一转换至目标UTM投影坐标系，再通过一次性合并的方式消除重叠区域，完全匹配“任意数据源记录即判定为贫民窟范围”的判定规则。
• 执行交集计算前，先对比圆形缓冲区与合并后贫民窟图层的空间范围，确认存在重叠区间后再开展面积计算，减少无效运算。

可复用的修正代码 #

# 正确创建WGS84坐标系下的兴趣点对象
lat = c(-22.868879628748203,-22.88511,-22.82166,-22.89692,-22.67945)
long = c(-43.237195000177564,-43.34278,-43.04717,-43.35168,-43.59667)
data_points = cbind.data.frame(long, lat) # 列顺序严格对应x=经度、y=纬度
coordinates(data_points) = c("long", "lat")
proj4string(data_points) = CRS("+init=epsg:4326")

# 定义匹配研究区域的UTM投影（南纬区域需开启south参数）
utm_crs <- CRS("+proj=utm +zone=23 +south +datum=NAD83 +units=m +no_defs +ellps=GRS80")
data_points_utm <- spTransform(data_points, utm_crs)

# 批量生成每个点周边2000米半径的圆形缓冲区
circular_buffers <- list()
for (i in 1:nrow(data_points_utm)) {
  spc <- spCircle(
    radius = 2000,
    centerPoint = c(
      x = data_points_utm@coords[i, "long"],
      y = data_points_utm@coords[i, "lat"]
    ),
    spID = i,
    spUnits = utm_crs
  )
  circular_buffers[[i]] <- spc[["spCircle"]]
}

# 统一所有贫民窟图层的投影后合并，消除重叠冗余
slum_layers <- list(shape1, shape2, shape3, shape4, shape5, shape6)
slum_layers_utm <- lapply(slum_layers, function(shp) spTransform(shp, utm_crs))
all_slums <- do.call(rbind, slum_layers_utm)
all_slums_union <- gUnaryUnion(all_slums)

# 逐点计算缓冲区范围内的贫民窟面积占比
slum_ratio <- numeric(length(circular_buffers))
for (i in 1:length(circular_buffers)) {
  circle <- circular_buffers[[i]]
  # 无空间重叠时占比直接设为0
  if (is.null(intersect(extent(circle), extent(all_slums_union)))) {
    slum_ratio[i] <- 0
    next
  }
  # 计算交集面积并换算占比
  slum_in_circle <- gIntersection(circle, all_slums_union)
  slum_area <- ifelse(is.null(slum_in_circle), 0, area(slum_in_circle))
  circle_area <- area(circle)
  slum_ratio[i] <- slum_area / circle_area
}

运行排查提示 #

• 若交集计算仍返回空值，可调用gIsValid()分别检查圆形缓冲区与合并后贫民窟多边形的拓扑有效性，修复自相交等拓扑错误后再开展计算。
• 所有面积计算需在米制投影坐标系下完成，避免在地理坐标系（角度单位）下直接计算导致结果偏差。

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