如何去除五大湖水生多边形，获取美加陆地县/CGU中心经纬度？

要去除湖泊内的虚假多边形，核心思路是通过陆地边界数据筛选出仅与陆地相交的行政单元，具体步骤如下：

步骤1：加载所需依赖包 #

除了你已使用的包，还需要rnaturalearth获取陆地边界数据（未安装的话先运行install.packages("rnaturalearth")）：

library(ggplot2)
library(dplyr)
library(geodata)
library(sf)
library(rnaturalearth)

步骤2：获取行政单元与陆地边界数据 #

# 获取美加二级行政单元（美国县、加拿大CGU）
aa <- gadm(country = c('CAN', 'USA'), level = 2, path=tempdir()) %>%
  st_as_sf()

# 获取北美陆地边界并转换为与行政单元一致的坐标系
north_america_land <- ne_land(scale = 10, returnclass = "sf") %>%
  filter(continent == "North America") %>%
  st_transform(crs = st_crs(aa))

步骤3：过滤仅与陆地相交的行政单元 #

通过空间相交判断，剔除完全位于水域的虚假单元：

# 标记与陆地有交集的行政单元
land_intersect_flag <- st_intersects(aa, north_america_land, sparse = FALSE)[,1]

# 筛选得到仅包含陆地的行政单元
aa_land_only <- aa[land_intersect_flag, ]

步骤4：计算陆地单元质心并验证可视化 #

# 获取陆地行政单元的中心经纬度
centroids_aa_land <- aa_land_only %>% 
  st_centroid(of_largest_polygon = TRUE)

# 绘制验证地图
ggplot() +
  geom_sf(data = north_america_land, fill = "lightgreen", alpha = 0.3) +
  geom_sf(data = aa_land_only, mapping = aes(geometry = geometry), color = "black") +
  geom_sf(data = centroids_aa_land, mapping = aes(geometry = geometry),
          color = "black", fill = 'red', shape = 21, size = 2) +
  coord_sf(xlim = c(-93, -81.5), ylim = c(41.5, 50), expand = F) +
  theme_bw() +
  theme(plot.margin = grid::unit(c(2,2,0,2), "mm"),
        text = element_text(size = 16),
        axis.text.x = element_text(size = 14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"),
        panel.grid.major = element_blank())

补充优化方案 #

如果需要更严格的筛选（要求单元质心必须位于陆地上），可以替换为以下逻辑：

# 先计算所有单元质心，再筛选质心在陆地上的单元
centroids_all <- aa %>% st_centroid(of_largest_polygon = TRUE)
centroid_land_flag <- st_intersects(centroids_all, north_america_land, sparse = FALSE)[,1]

aa_land_only <- aa[centroid_land_flag, ]
centroids_aa_land <- centroids_all[centroid_land_flag, ]

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