如何在ggplot中移除与多边形不相交的网格单元

问题描述 #

我基于近海观测的鸟类数据集创建了网格，并使用mapdata包生成了海岸多边形。我希望将网格筛选为仅包含海洋区域或与海岸相交的单元。
我可以通过将多边形叠加在网格上来实现视觉效果，但无法统计符合条件的网格单元数量。我尝试使用st_difference裁剪网格但未成功，以下是我的脚本及生成的图表，同时提供了网格数据和创建网格的步骤注释。

创建网格的代码（上下文参考） #

# 加载所需库
library(ggplot2)
library(dplyr)
library(tidyr)
library(maps)
library(mapdata)
library(ggpubr)
library(ggspatial)
library(sf)

# 基于鸟类数据集创建网格的代码
# 下方已提供这段代码生成的网格数据，此处仅作上下文参考
# data_sf = st_as_sf(bird_data, 
#                    coords = c("decimalLongitude", "decimalLatitude"),
#                    crs = 4326)
# bbox <- st_bbox(data_sf)
# tmp_grid <- st_make_grid(bbox, cellsize=5)
# grid <- st_as_sf(tmp_grid)

创建海岸多边形与裁剪网格的代码 #

## 创建海岸线轮廓的代码

# 加载高精度世界地图数据
world_map <- map_data("world2Hires")
# 提取加拿大和美国的区域多边形
country_polygons <- world_map[world_map$region %in% c('Canada', 'USA'), ]
# 转换经度格式（适配西经范围）
country_polygons$long = (360 - country_polygons$long)*-1

# 将地图数据转换为sf空间对象
land_sf <- st_as_sf(country_polygons, coords = c("long", "lat"), 
                    crs= 4326, remove = FALSE)
# 合并陆地多边形
land_sf <- st_union(land_sf)

# 裁剪网格，保留不与陆地重叠的部分
clipped_grid <- st_difference(grid,land_sf)

# 设置绘图范围
lons = c(-75, -50)     # 经度范围（最大，最小）
lats = c(40, 55)       # 纬度范围（最小，最大）

# 绘制网格与陆地多边形的叠加图
ggplot() +
  geom_sf(data = clipped_grid, color = "red") +
  theme_classic() + 
  geom_sf(data = land_sf) +
  coord_sf(xlim = lons, ylim = lats) 

生成的叠加图表 #

网格数据 #

# 网格数据
grid = structure(list(x = structure(list(structure(list(structure(c(-73.9167, 
-68.9167, -68.9167, -73.9167, -73.9167, 40.0167, 40.0167, 45.0167, 
45.0167, 40.0167), dim = c(5L, 2L))), class = c("XY", "POLYGON", 
"sfg")), structure(list(structure(c(-68.9167, -63.9167, -63.9167, 
-68.9167, -68.9167, 40.0167, 40.0167, 45.0167, 45.0167, 40.0167
), dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-63.9167, -58.9167, -58.9167, -63.9167, -63.9167, 
    40.0167, 40.0167, 45.0167, 45.0167, 40.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-58.9167, -53.9167, -53.9167, -58.9167, -58.9167, 
    40.0167, 40.0167, 45.0167, 45.0167, 40.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-53.9167, -48.9167, -48.9167, -53.9167, -53.9167, 
    40.0167, 40.0167, 45.0167, 45.0167, 40.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-73.9167, -68.9167, -68.9167, -73.9167, -73.9167, 
    45.0167, 45.0167, 50.0167, 50.0167, 45.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-68.9167, -63.9167, -63.9167, -68.9167, -68.9167, 
    45.0167, 45.0167, 50.0167, 50.0167, 45.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-63.9167, -58.9167, -58.9167, -63.9167, -63.9167, 
    45.0167, 45.0167, 50.0167, 50.0167, 45.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-58.9167, -53.9167, -53.9167, -58.9167, -58.9167, 
    45.0167, 45.0167, 50.0167, 50.0167, 45.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-53.9167, -48.9167, -48.9167, -53.9167, -53.9167, 
    45.0167, 45.0167, 50.0167, 50.0167, 45.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-73.9167, -68.9167, -68.9167, -73.9167, -73.9167, 
    50.0167, 50.0167, 55.0167, 55.0167, 50.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-68.9167, -63.9167, -63.9167, -68.9167, -68.9167, 
    50.0167, 50.0167, 55.0167, 55.0167, 50.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-63.9167, -58.9167, -58.9167, -63.9167, -63.9167, 
    50.0167, 50.0167, 55.0167, 55.0167, 50.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-58.9167, -53.9167, -53.9167, -58.9167, -58.9167, 
    50.0167, 50.0167, 55.0167, 55.0167, 50.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(-53.9167, -48.9167, -48.9167, -53.9167, -53.9167, 
    50.0167, 50.0167, 55.0167, 55.0167, 50.0167), dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg"))), class = c("sfc_POLYGON", 
"sfc"), precision = 0, bbox = structure(c(xmin = -73.9167, ymin = 40.0167, 
xmax = -48.9167, ymax = 55.0167), class = "bbox"), crs = structure(list(
    input = "EPSG:4326", wkt = "GEOGCRS[\"WGS 84\",
    ENSEMBLE[\"World Geodetic System 1984 ensemble\",
        MEMBER[\"World Geodetic System 1984 (Transit)\"],
        MEMBER[\"World Geodetic System 1984 (G730)\"],
        MEMBER[\"World Geodetic System 1984 (G873)\"],
        MEMBER[\"World Geodetic System 1984 (G1150)\"],
        MEMBER[\"World Geodetic System 1984 (G1674)\"],
        MEMBER[\"World Geodetic System 1984 (G1762)\"],
        MEMBER[\"World Geodetic System 1984 (G2139)\"],
        ELLIPSOID[\"WGS 84\",6378137,298.257223563,
            LENGTHUNIT[\"metre\",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM[\"Greenwich\",0,
        ANGLEUNIT[\"degree\",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS[\"geodetic latitude (Lat)\",north,
            ORDER[1],
            ANGLEUNIT[\"degree\",0.0174532925199433]],
        AXIS[\"geodetic longitude (Lon)\",east,
            ORDER[2],
            ANGLEUNIT[\"degree\",0.0174532925199433]],
    USAGE[
        SCOPE[\"Horizontal component of 3D system.\"],
        AREA[\"World.\"],
        BBOX[-90,-180,90,180]],
    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, 
15L), class = c("sf", "data.frame"), sf_column = "x", agr = structure(integer(0), class = "factor", levels = c("constant", 
"aggregate", "identity"), names = character(0)))

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