ggplot绘制美国双变量县域地图时边界侵入湖泊问题求助

我尝试使用ggplot制作美国双变量地图，但生成的地图中县域边界超出范围延伸至湖泊内，尤其是五大湖区域。我先后尝试使用fipio::fips_county()和tigris::counties(year = 2020, class = "sf", resolution = "20m")获取shapefile/坐标，两者输出的地图均存在该问题。请问有解决方法吗？

用户提供的示例代码：

library(tidyverse)
library(ggplot2)
library(cowplot)
library(sf)
library(biscale)
library(fipio)

all_counties <- tigris::counties(year = 2020, class = "sf", resolution = "20m")

all_counties <- all_counties %>%
  # mutate(geometry = fips_geometry(GEOID)) %>%
  filter(!grepl("^(02|15)", GEOID))


all_counties <- bi_class(all_counties, x = ALAND, y = AWATER, style = "quantile", dim = 3) 

# create map
map <- ggplot() +
  geom_sf(data = all_counties, mapping = aes(fill = bi_class, geometry=geometry), color = "white", size = 1, show.legend = FALSE) +
  bi_scale_fill(pal = "GrPink", dim = 3) +
  bi_theme()
map

legend <- bi_legend(pal = "GrPink",
                    dim = 3,
                    xlab = "More Land ",
                    ylab = "More Water ",
                    size = 7)
finalPlot <- ggdraw() +
  draw_plot(map, 0, 0, 1, 1) +
  draw_plot(legend, 0.05, .2, 0.2, 0.2)

finalPlot

可行的解决方法 #

方法1：用美国本土陆地边界裁剪县域数据 #

核心是获取不含湖泊的美国陆地边界，将县域数据与陆地做交集运算，剔除侵入湖泊的部分：

library(tidyverse)
library(ggplot2)
library(cowplot)
library(sf)
library(biscale)
library(tigris)

# 获取并过滤县域数据
all_counties <- counties(year = 2020, class = "sf", resolution = "20m") %>%
  filter(!grepl("^(02|15)", GEOID)) # 移除阿拉斯加、夏威夷

# 获取美国本土陆地边界（合并非本土外的州）
us_land <- states(year = 2020, class = "sf", resolution = "20m") %>%
  filter(!grepl("^(02|15|60|66|69|72|78)", GEOID)) %>% # 移除海外领地和非本土州
  st_union() %>% # 合并为单一陆地多边形
  st_sf()

# 裁剪县域数据，仅保留与陆地重叠的部分
all_counties_clipped <- st_intersection(all_counties, us_land)

# 继续双变量地图流程
all_counties_clipped <- bi_class(all_counties_clipped, x = ALAND, y = AWATER, style = "quantile", dim = 3) 

map <- ggplot() +
  geom_sf(data = all_counties_clipped, aes(fill = bi_class), color = "white", size = 0.1, show.legend = FALSE) +
  bi_scale_fill(pal = "GrPink", dim = 3) +
  bi_theme()

legend <- bi_legend(pal = "GrPink",
                    dim = 3,
                    xlab = "更多陆地 ",
                    ylab = "更多水域 ",
                    size = 7)

finalPlot <- ggdraw() +
  draw_plot(map, 0, 0, 1, 1) +
  draw_plot(legend, 0.05, .2, 0.2, 0.2)

finalPlot

方法2：用rmapshaper包精准裁剪 #

如果st_intersection处理效果不佳，可以用rmapshaper的ms_clip工具做更精准的拓扑裁剪：

library(rmapshaper)

# 替换裁剪步骤
all_counties_clipped <- ms_clip(target = all_counties, clip = us_land)

方法3：直接使用预处理好的制图边界数据 #

tidycensus包提供的cartographic boundary数据已经预先处理了湖泊侵入问题，适合直接制图：

library(tidycensus)

# 获取带地理信息的县域数据（已处理湖泊问题）
all_counties <- get_acs(geography = "county", 
                        variables = "B01003_001", # 任意选一个变量，仅为获取地理数据
                        year = 2020,
                        geometry = TRUE,
                        cb = TRUE) %>%
  filter(!grepl("^(02|15)", GEOID))

注：cb = TRUE参数获取的是简化版边界，文件更小，制图速度更快，同时已经剔除了湖泊内的县域边界。

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