美国人口普查区域地图几何图形报错问题求助

尝试按人口普查区域绘制美国数据并添加区域轮廓，但在创建区域多边形时遇到拓扑错误。参考相关代码执行到合并区域步骤时抛出错误，具体代码及错误信息如下：

library(tidyverse)
library(tigris)
library(sf)
library(maps)

div_dat <- states(cb = FALSE, resolution = '20m') %>%
  st_drop_geometry() %>%
  select(NAME, DIVISION) %>%
  mutate(ID = tolower(NAME))

# get state data, convert to sf, join with division data
states <- maps::map("state", plot = FALSE, fill = TRUE) %>%
  st_as_sf() %>%
  left_join(div_dat)

# create division polygons
div <- states %>%
  group_by(DIVISION) %>% 
  summarize()

# plot it
ggplot() + 
  theme_void() +
  geom_sf(data = states, 
          aes(fill = as.numeric(DIVISION)), 
          color = 'white') +
  geom_sf(data = div, 
          color = 'black', 
          fill = NA,
          size = 1) +
  scale_fill_viridis_c() +
  coord_sf(crs = 5070) +
  labs(fill = NULL)

执行创建区域多边形的代码时出现以下错误：

Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented,  : 
  Loop 0 is not valid: Edge 236 crosses edge 239

已尝试关闭球面几何和使用st_make_valid()，但问题仍未解决，可通过以下方法修复：

方法1：替换数据源，使用tigris官方Shapefile构建区域 #

maps包的state数据可能存在拓扑瑕疵，直接使用tigris返回的规范sf对象构建区域，从根源避免拓扑错误：

library(tidyverse)
library(tigris)
library(sf)

# 直接从tigris获取带DIVISION信息的states数据
states <- states(cb = FALSE, resolution = '20m') %>%
  select(NAME, DIVISION)

# 用st_union合并同区域的州几何
div <- states %>%
  group_by(DIVISION) %>%
  summarize(geometry = st_union(geometry)) %>%
  ungroup()

# 绘图
ggplot() + 
  theme_void() +
  geom_sf(data = states, 
          aes(fill = as.numeric(DIVISION)), 
          color = 'white') +
  geom_sf(data = div, 
          color = 'black', 
          fill = NA,
          size = 1) +
  scale_fill_viridis_c() +
  coord_sf(crs = 5070) +
  labs(fill = NULL)

方法2：修复maps数据的拓扑问题 #

若必须使用maps包数据，可先转换为平面坐标系再修复拓扑：

library(tidyverse)
library(tigris)
library(sf)
library(maps)

div_dat <- states(cb = FALSE, resolution = '20m') %>%
  st_drop_geometry() %>%
  select(NAME, DIVISION) %>%
  mutate(ID = tolower(NAME))

# 获取state数据后，先转平面CRS再修复拓扑
states <- maps::map("state", plot = FALSE, fill = TRUE) %>%
  st_as_sf() %>%
  st_set_crs(4326) %>%
  st_transform(5070) %>% # 转换为目标平面坐标系
  st_make_valid() %>% # 修复基础拓扑
  st_buffer(0) %>% # 消除自相交边
  left_join(div_dat)

# 创建区域多边形
div <- states %>%
  group_by(DIVISION) %>% 
  summarize(geometry = st_union(geometry)) %>%
  ungroup()

# 绘图
ggplot() + 
  theme_void() +
  geom_sf(data = states, 
          aes(fill = as.numeric(DIVISION)), 
          color = 'white') +
  geom_sf(data = div, 
          color = 'black', 
          fill = NA,
          size = 1) +
  scale_fill_viridis_c() +
  coord_sf(crs = 5070) +
  labs(fill = NULL)

方法3：正确关闭S2球面几何并重新处理 #

若之前关闭S2的方式不正确，可在会话开始时全局禁用：

sf::sf_use_s2(FALSE)

之后再运行原代码，或结合st_union代替默认summarize完成几何合并。

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