Create a spatial polygons data frame that preserves overlapping features with another spatial polygons data frame but does not clip polygon extent

Question

I have been using the intersect() function from the raster package in R to clip a spatial polygons data frame (HUC-4 watersheds) to the extent of another spatial polygons data frame (a region consisting of Colorado, Idaho, Montana, Utah, and Wyoming).

I want to preserve the entire extent of the spatial polygons that overlap with the spatial data frame I am clipping to. Using intersect() clips the HUC-4 watersheds so that they do not extend past the extent of the states being clipped to.

The watershed data that I am using can be downloaded from: ftp://rockyftp.cr.usgs.gov/vdelivery/Datasets/Staged/Hydrography/WBD/National/GDB/ (WBD_National_GDB.zip).

The data for the region encompassing Colorado, Utah, Idaho, Wyoming, and Montana was extracted from county data available here: https://catalog.data.gov/dataset/tiger-line-shapefile-2017-nation-u-s-current-county-and-equivalent-national-shapefile.

The code I am using to do the clip with the intersect() function is as follows:

library(raster)
library(dplyr)
library(spdplyr)
library(rgdal)
library(rgeos)

# albers equal area projection
proj <- CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ")

counties <- readOGR(dsn = "./data/tl_2017_us_county/tl_2017_us_county.shp")

# filtering out only counties in our 5 states of interest
counties <- counties %>%
  filter(STATEFP %in% c("08", "16", "30", "49", "56"))

# transforming to albers projection
counties <- spTransform(counties, proj)

# create a region shapefile (to clip watersheds with)
region <- gUnaryUnion(counties)

# Make Region into a SpatialPolygonsDataFrame
row.names(region) <- as.character(1:length(region))
region_data <- c("West")
region_data <- as.data.frame(region_data)
colnames(region_data) <- "Region"

region <- SpatialPolygonsDataFrame(region, region_data)

file <- "./data/WBD_National_GDB/WBD_National_GDB.gdb"

# huc4 watersheds
huc4 <- readOGR(dsn = file, layer = "WBDHU4")

# transforming to albers projection
huc4 <- spTransform(huc4, proj)

# selecting only huc4 watersheds that intersect with our states of interest
huc4_clip <- raster::intersect(huc4, region)

# plot the result
plot(huc4_clip)

I want an output file that does not clip the extent of the spatial polygons that are on the edge of the region of interest, but does not include any spatial polygons that do not directly overlap with the region of interest. Are there any other functions I can use that are similar to intersect() but that do not clip the extent of the spatial polygons on the region border?

Answer 1

If I understand the question properly, you could use function gIntersects to find out which watersheds intersect your region, and then extract only those from the huc4 dataset. In practice, something like this could work:

intersects <- which(gIntersects(huc4, region, byid = TRUE))
huc4_clip  <- huc4[intersects, ]