Я использовал функцию intersect() из растрового пакета в R, чтобы обрезать фрейм данных пространственных полигонов (водоразделы HUC-4) до экстента другого фрейма данных пространственных полигонов (регион, состоящий из Колорадо, Айдахо, Монтаны, Юта и Вайоминг).
Я хочу сохранить весь экстент пространственных многоугольников, которые перекрываются с пространственным фреймом данных, к которому я привязываюсь.Использование intersect() обрезает водоразделы HUC-4 таким образом, чтобы они не выходили за пределы ограничиваемых состояний.
Используемые мной данные водораздела можно загрузить с: ftp: //rockyftp.cr.usgs.gov/vdelivery/Datasets/Staged/Hydrography/WBD/National/GDB/ (WBD_National_GDB.zip).
Данные для региона, охватывающего Колорадо, штат Юта, Айдахо,Вайоминг и Монтана были извлечены из данных графства, доступных здесь: https://catalog.data.gov/dataset/tiger-line-shapefile-2017-nation-u-s-current-county-and-equivalent-national-shapefile.
Код, который я использую для клипа с функцией intersect(), выглядит следующим образом:
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Вам нужен выходной файл, который не обрезает экстент пространственных многоугольников, которые находятся на краю интересующей области, но не включает в себя пространственные многоугольники, которые непосредственно не перекрываются с интересующей областью.Могу ли я использовать какие-либо другие функции, которые похожи на intersect(), но которые не ограничивают экстент пространственных многоугольников на границе региона?