如何自定义R语言sf包st_union处理XYZ多边形的Z/M值计算逻辑？

方案一：将交叉点Z/M值设为NA #

核心思路是先忽略第三维度完成Union操作，再恢复3D结构并给新增的交叉点赋值NA。具体步骤如下：

• 提取原始多边形的2D几何信息（丢弃Z/M维度）
• 对2D几何执行st_union
• 将Union结果转回3D维度，对比原始坐标，给新增的交叉点的Z/M设为NA

代码示例：

library(sf)

# 生成示例数据（复用你的代码）
m = as.matrix(data.frame(X=c(0,1,1,0,0), Y=c(0,0,1,1,0), Z=c(1,1,1,1,1)))
p = st_polygon(x=list(m), dim="XYZ")
n = 3
l = vector("list", n)
for (i in 1:n)
  l[[i]] = p + 1 * c(runif(2), i)
s = st_sfc(l)

# 步骤1：提取2D几何
s_2d = st_zm(s, drop = TRUE)
# 步骤2：2D Union
union_2d = st_union(s_2d)
# 步骤3：转回3D并处理Z值
union_3d = st_zm(union_2d, dim = "XYZ")

# 获取原始所有坐标的XY-Z映射
original_coords = do.call(rbind, lapply(st_coordinates(s), function(x) x[,c("X","Y","Z")]))
original_coords = unique(original_coords)
original_map = setNames(original_coords$Z, paste(original_coords$X, original_coords$Y, sep = "_"))

# 修改Union后的Z值：原始存在的点保留Z，新增点设为NA
union_coords = st_coordinates(union_3d)
union_coords$Z = sapply(paste(union_coords$X, union_coords$Y, sep = "_"), function(key) {
  ifelse(key %in% names(original_map), original_map[key], NA)
})

# 重新构建3D多边形
union_final = st_as_sfc(st_as_text(st_polygon(list(union_coords[,c("X","Y","Z")]))), dim = "XYZ")

# 查看结果
plot(union_final, col = sf.colors(categorical = TRUE, alpha = .5))
title("st_union with NA for new points' Z")
aa = as.data.frame(st_coordinates(union_final))
text(aa$X, aa$Y, aa$Z)

方案二：自定义交叉点Z/M值的计算逻辑 #

如果需要自定义交叉点的Z/M值（比如取最大值、最小值或其他规则），可以手动提取交叉线段，计算交叉点并赋值自定义值，再重构几何。以下示例以取交叉线段两端Z值的最大值为例：

library(sf)

# 复用示例数据
m = as.matrix(data.frame(X=c(0,1,1,0,0), Y=c(0,0,1,1,0), Z=c(1,1,1,1,1)))
p = st_polygon(x=list(m), dim="XYZ")
n = 3
l = vector("list", n)
for (i in 1:n)
  l[[i]] = p + 1 * c(runif(2), i)
s = st_sfc(l)

# 获取所有多边形的线段（每个多边形拆分为线段）
get_segments = function(poly) {
  coords = st_coordinates(poly)[-nrow(st_coordinates(poly)),] # 去掉闭合点
  lapply(1:(nrow(coords)-1), function(i) {
    st_linestring(rbind(coords[i,c("X","Y","Z")], coords[i+1,c("X","Y","Z")]), dim = "XYZ")
  })
}
all_segments = unlist(lapply(s, get_segments), recursive = FALSE)
all_segments_sfc = st_sfc(all_segments)

# 找到所有相交的线段对，筛选出真正的交叉点（非端点重合）
intersections = st_intersection(all_segments_sfc)
cross_points = intersections[st_dimension(intersections) == 0]

# 自定义计算交叉点的Z值：取两条交叉线段所有端点Z的最大值
custom_z = sapply(cross_points, function(point) {
  # 找到包含该点的线段
  containing_segs = all_segments_sfc[st_contains(all_segments_sfc, point, sparse = FALSE)]
  # 提取这些线段的端点Z值
  seg_coords = lapply(containing_segs, st_coordinates)
  z_values = unlist(lapply(seg_coords, function(x) x[,"Z"]))
  max(z_values) # 此处可替换为自定义逻辑，比如min、自定义函数等
})

# 先获取2D Union结果，再替换交叉点的Z值
union_2d = st_union(st_zm(s, drop = TRUE))
union_coords = st_coordinates(union_2d)
original_coords = do.call(rbind, lapply(st_coordinates(s), function(x) x[,c("X","Y","Z")]))
original_coords = unique(original_coords)

# 替换Union结果中的Z值
union_coords$Z = sapply(1:nrow(union_coords), function(i) {
  xy = st_point(c(union_coords[i,"X"], union_coords[i,"Y"]))
  # 检查是否是交叉点
  match_idx = which(st_distance(xy, cross_points) < 1e-8)
  if(length(match_idx) > 0) {
    custom_z[match_idx]
  } else {
    # 匹配原始坐标的Z值
    original_match = which(abs(original_coords$X - union_coords[i,"X"]) < 1e-8 & abs(original_coords$Y - union_coords[i,"Y"]) < 1e-8)
    original_coords$Z[original_match]
  }
})

# 重构3D多边形
union_custom = st_as_sfc(st_as_text(st_polygon(list(union_coords[,c("X","Y","Z")]))), dim = "XYZ")

# 查看结果
plot(union_custom, col = sf.colors(categorical = TRUE, alpha = .5))
title("st_union with custom Z values for new points")
aa = as.data.frame(st_coordinates(union_custom))
text(aa$X, aa$Y, aa$Z)

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