按纬度筛选Shapefile：保留30-fm图层及北纬34.4以上区域

你的问题出在使用fortify()将sf空间对象转换为普通数据框后，丢失了空间几何的结构化属性，导致无法正确基于纬度筛选空间区域。直接在sf对象上操作就能解决这个问题，以下是修改后的完整代码和关键步骤说明：

关键修改点 #

• 放弃fortify()，直接保留sf对象进行空间筛选
• 用sf自带的空间函数或坐标提取来筛选纬度高于34.4的区域

修改后的代码 #

packages <- c('rnaturalearth', 'sf', 'tidyverse', 'ggsignif', 
              'data.table', 'patchwork', 'ggOceanMaps',
              "tools", "dplyr", 'sp', 'readr','mapproj',
              'cowplot', 'raster')

installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

invisible(lapply(packages, library, character.only = TRUE))

# 读取Shapefile为sf对象
path.RCA <- ("~/Downloads/Rockfish_Conservation_Area_-_R7_-_CDFW_[ds3144]/")
fnam.RCA <- "Rockfish_Conservation_Area_-_R7_-_CDFW_[ds3144].shp"
RCA <- st_read(dsn = path.RCA, layer = file_path_sans_ext(fnam.RCA))

# 第一步：筛选30-fm且Coastwide的图层（直接在sf对象上操作）
fathom_30 <- RCA %>%
  filter(str_detect(area_name, "30-fm") & str_detect(area_name, "Coastwide"))

# 第二步：筛选纬度高于34.4的区域
# 方法1：提取每个要素的坐标，筛选所有点纬度>34.4的要素（更精准）
fathom_30_filtered <- fathom_30 %>%
  rowwise() %>%
  mutate(min_lat = min(st_coordinates(geometry)[,2]),
         max_lat = max(st_coordinates(geometry)[,2])) %>%
  filter(max_lat > 34.4) %>%
  ungroup() %>%
  select(-min_lat, -max_lat)

# 方法2：用空间裁剪直接截取y>34.4的区域（更快速）
# bbox <- st_bbox(c(xmin = -121.5, ymin = 34.4, xmax = -117, ymax = 36))
# fathom_30_filtered <- st_crop(fathom_30, bbox)

# 绘制地图
world <- ne_countries(scale = "medium", returnclass = "SF")

socal_map <- ggplot() +
  geom_sf(data = world) +
  geom_sf(data = fathom_30_filtered, size = 1, color = 'lightgreen')+
  coord_sf(xlim = c(-117, -121.5), ylim = c(34.4, 36), expand = FALSE)+
  scale_x_continuous(breaks = c(seq(from = -121, to = -117, by = 2))) +
  scale_y_continuous(breaks = c(34.4, 35, 36))+
  theme(legend.position = 'bottom',
        strip.background = element_blank(), 
        legend.title = element_blank(),
        strip.text.x = element_text(
          size = 12, hjust = 0, face = "bold"),
        panel.border = element_rect(color = "red", 
                                    fill = NA, 
                                    linewidth = 1),
        axis.line.y = element_blank(), axis.line.x = element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank())

socal_map

说明 #

• 方法1通过提取每个空间要素的所有坐标点，判断其最大纬度是否高于34.4，确保只保留完全或部分位于目标纬度以上的要素
• 方法2用st_crop直接裁剪出指定纬度范围的区域，操作更高效，适合大范围数据
• 地图的ylim也调整为从34.4开始，聚焦目标区域

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