如何删除与 ggplot 中的多边形相交以外的任何网格单元

Question

我根据近海观察到的鸟类数据集创建了一个网格，并使用

mapdata

包创建了海岸的多边形。我想将网格缩小为仅出现在海洋上空或与海岸相交的单元格。

我可以通过在网格上分层多边形来实现这种美感，但这并没有告诉我有多少网格单元满足我的上述标准。我尝试使用

st_difference

剪辑网格，但没有成功，如我的下面的脚本和附图所示。

我使用

dput

为您提供网格数据，但也注释掉了我创建网格所采取的步骤，以防提供对问题的任何见解。

# Load required libraries
library(ggplot2)
library(dplyr)
library(tidyr)
library(maps)
library(mapdata)
library(ggpubr)
library(ggspatial)
library(sf)

# Where I create the grid based on my data. 
# I have pasted the product of this code above ... 
# but am providing it for context.

#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)

## I create the outline of the coastline using the following script.

#load world map data
world_map <- map_data("world2Hires")
# Extract Canada and USA polygons
country_polygons <- world_map[world_map$region %in% c('Canada', 'USA'), ]
country_polygons$long = (360 - country_polygons$long)*-1

#convert to sf object
land_sf <- st_as_sf(country_polygons, coords = c("long", "lat"), 
                    crs= 4326, remove = FALSE)
land_sf <- st_union(land_sf)

#clip grid and show where it does not overlap with land_sf
clipped_grid <- st_difference(grid,land_sf)

lons = c(-75, -50)     # (max longitude, min longitude)
lats = c(40, 55)       # (min latitude, max latitude)

#plot the grid and the land polygon
ggplot() +
  geom_sf(data = clipped_grid, color = "red") +
  theme_classic() + 
  geom_sf(data = land_sf) +
  coord_sf(xlim = lons, ylim = lats) 

数据：

#grid data
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\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    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)))

Answer 1

我们可以使用

sf::st_intersects

来获取哪些多边形（即

grid

）与点（即

land_sf

）相交，然后简单地使用

[

过滤掉不相交的多边形。

clipped_grid <- grid[sf::st_intersects(grid, land_sf, sparse = FALSE),]  

ggplot() +
  geom_sf(data = clipped_grid, color = "red") +
  theme_classic() + 
  geom_sf(data = land_sf) +
  coord_sf(xlim = lons, ylim = lats) 

^{创建于 2024-02-21，使用 reprex v2.0.2}