R 基于另一个 sf 特征绘制 sf 对象的子集(对于某些状态)

当您使用sf时，可用的工具之一是st_filter()，这是一种空间过滤器，允许您 Select 相交的几何图形.尽管根据您的需要，这一点可能就足够了，因为相交的坚硬多边形会延伸很远.在下面的示例中，还包括ca.az.nm个多边形的相交和ca.az.nm个边界框的裁剪.

library(sf)
library(dplyr)


url_lst = list(state_zip = "https://www.weather.gov/source/gis/Shapefiles/County/s_22mr22.zip",
               phm_zip = 'http://prism.oregonstate.edu/projects/public/phm/phm_us_shp.zip')

tempd <- tempdir()

zip_lst <- lapply(url_lst, \(url){
  destfile <- file.path(tempdir(),basename(url))
  if(!file.exists(destfile)) download.file(url, destfile = destfile, mode = "wb", overwrite=TRUE)
  destfile
})

state_boundaries <- read_sf(paste0("/vsizip/",zip_lst$state_zip))
shp_hardness <- read_sf(paste0("/vsizip/",zip_lst$phm_zip))

# check crs
st_crs(state_boundaries) == st_crs(shp_hardness)
#> [1] TRUE

ca.az.nm <- state_boundaries %>% filter(STATE %in% c("CA", "AZ", "NM"))

shp_hardness_subset <- shp_hardness %>%
  filter(str_detect(ZONE, '9b|10a|10b|6a|6b|7a|7b')) %>% 
  # sf complains about few geometries in phm_us_shp.zip 
  st_make_valid() %>% 
  # spatial filter by ca.az.nm, default predicate is "intersects"
  st_filter(ca.az.nm)

p1 <- ggplot() +
  geom_sf(data =  ca.az.nm) + 
  geom_sf(data =  shp_hardness_subset, aes(fill = ZONE)) +
  theme_bw()

# intersection by ca.az.nm:
shp_hardness_subset_clip <- st_intersection(shp_hardness_subset, ca.az.nm)
p2 <- ggplot() +
  geom_sf(data =  ca.az.nm) + 
  geom_sf(data =  shp_hardness_subset_clip, aes(fill = ZONE)) +
  theme_bw()

# crop by bounding box of ca.az.nm:
shp_hardness_subset_crop <- st_crop(shp_hardness_subset, ca.az.nm) 
p3 <- ggplot() +
  geom_sf(data =  ca.az.nm) + 
  geom_sf(data =  shp_hardness_subset_crop, aes(fill = ZONE)) +
  theme_bw() 

patchwork::wrap_plots(p1, p2, p3, ncol = 1, guides = "collect")

^{创建于2023-05-11年第reprex v2.0.2页}