在R中使用st_voronoi生成泰森多边形并绘图时遇错，求解决方法

问题背景 #

基于给定点位数据，使用sf包的st_voronoi生成泰森多边形，并用基础plot函数叠加原始点绘制时，出现如下报错：

Error in if (xsize * ysize * n > prod(total_size)) {: missing value where TRUE/FALSE needed
Traceback:

1. plot(voronoi_grid, col = "lightblue", border = "black", lwd = 1.5)
2. plot(voronoi_grid, col = "lightblue", border = "black", lwd = 1.5)
3. plot.sf(voronoi_grid, col = "lightblue", border = "black", lwd = 1.5)
4. .get_layout(st_bbox(x), min(max.plot, length(cols)), par("din"), 
 .     key.pos, key.width)
5. vapply(1:n, function(x) size(x, n, asp), 0)
6. FUN(X[[i]], ...)
7. size(x, n, asp)

用户提供的text.txt数据：

name,long,lat,water_level,elevation,depth
EM_01,18.553392,-34.07027,14.4,20.358,63.0
EM_27,18.574777,-34.068709,16.196,19.966,48.0
EM_29,18.613985,-34.053271,18.766,25.477,39.0
EM_20,18.654089,-34.045177,20.102,36.502,45.0

用户原代码：

library(sf) # 'simple features' representations of spatial objects

file = 'text.txt'
#-- import
cfaq <- read.csv(file, header = 1, sep = ',', dec = '.') 
#- set as a spatial feature with xy coords an existing projection
cfaq.sf <- st_as_sf(cfaq, coords=c("long", "lat"), crs = 4326) #wgs83
#- transform to local crs
cfaq.sf <- st_transform(cfaq.sf, crs = 32734) #utm 34s
# Voronoi tesselation
voronoi_grid <- st_voronoi(cfaq.sf)

# basic plot with points on top of the Voronoi 
plot(voronoi_grid, col = "lightblue", border = "black", lwd = 1.5)
plot(cfaq.sf, col = "red", pch = 16, cex = 2, add = TRUE)
#text(cfaq.sf$long, cfaq.sf$lat, labels = cfaq.sf$name, pos = 3, cex = 0.8)
# Add a legend
#legend("bottomright", legend = "Points", col = "red", pch = 16, bty = "n", pt.cex = 1.5)
# Add a title
title(main = "Voronoi Diagram with Points")

报错原因 #

st_voronoi()默认返回GEOMETRYCOLLECTION类型对象，直接用plot.sf绘制时，会触发布局计算逻辑错误，导致缺失值判断失败。

修正方案 #

核心是将Voronoi几何集合转换为单个多边形的sf对象，同时修正点位标注的坐标提取方式：

完整修正代码 #

library(sf)

# 读取数据
file <- 'text.txt'
cfaq <- read.csv(file, header = TRUE, sep = ',', dec = '.') 

# 转换为sf对象并投影转换
cfaq.sf <- st_as_sf(cfaq, coords = c("long", "lat"), crs = 4326)
cfaq.sf <- st_transform(cfaq.sf, crs = 32734) # UTM 34S

# 生成Voronoi多边形并转换为sf对象
# 1. 生成缓冲范围作为Voronoi裁剪边界（避免无限大多边形）
vor_bbox <- st_union(cfaq.sf) %>% st_buffer(dist = 1000) # 缓冲距离可按需调整
# 2. 生成Voronoi并提取多边形转为sf对象
voronoi_grid <- st_voronoi(st_union(cfaq.sf), envelope = vor_bbox) %>%
  st_collection_extract("POLYGON") %>%
  st_sf()

# 叠加绘制
plot(voronoi_grid$geometry, col = "lightblue", border = "black", lwd = 1.5)
plot(cfaq.sf$geometry, col = "red", pch = 16, cex = 2, add = TRUE)

# 提取坐标添加点位名称
coords <- st_coordinates(cfaq.sf)
text(coords[,1], coords[,2], labels = cfaq.sf$name, pos = 3, cex = 0.8)

# 添加图例和标题
legend("bottomright", legend = "监测点", col = "red", pch = 16, bty = "n", pt.cex = 1.5)
title(main = "监测点泰森多边形分布图")

关键修改说明 #

• Voronoi对象转换：用st_collection_extract("POLYGON")提取单个多边形，再通过st_sf()转为标准sf对象，解决plot布局计算问题。
• 裁剪范围设置：通过st_union+st_buffer生成缓冲范围作为Voronoi的envelope，避免生成无限延伸的多边形，优化绘图效果。
• 点位坐标提取：sf对象坐标存储在geometry字段中，需用st_coordinates()提取后才能用于text()标注。

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