Search code examples
rggplot2geometrymapsggspatial

How to remove any grid cells beyond where they intersect with a polygon in ggplot

I have created a grid based on a dataset of birds observed offshore and created a polygon of the coast using the mapdata package. I want to reduce the grid to only cells that occur over ocean or that intersect with the coast.

I can achieve the aesthetic of this by layering the polygon over the grid, but that doesn't show me how many grid cells meet my above criteria. I have tried to clip the grid with no success using st_difference as demonstrated by my below script and the attached plot.

I have used dput to provide you with the grid data but have also commented out the steps I took to create the grid in case that provides any insight into the problem.

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

Data:

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

  • We can use sf::st_intersects to get which polygons (i.e. grid) are intersecting with the points (i.e. land_sf), and then simply use [ to filter out non-intersecting ones.

    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)

    Created on 2024-02-21 with reprex v2.0.2