r - 根据另一个 sf 特征绘制 sf 对象的子集(对于某些状态)

Question

此问题与 how map certain USDA hardiness zones in R 相关，但我只需要 CA、NM 和 AZ 的 map 。我应该在 filter 命令中输入什么才能仅具有 3 种状态？

library(sf)
library(tidyverse)
library(USAboundaries)

state_zip <- "https://www.weather.gov/source/gis/Shapefiles/County/s_22mr22.zip"
 download.file(state_zip, destfile = ".", junkpaths=TRUE, overwrite=TRUE)
 
 unzip("state_zip.zip", junkpaths = TRUE, exdir = ".")

 state_boundaries <- read_sf(".s_22mr22.shp")

temp_shapefile <- tempfile()
download.file('http://prism.oregonstate.edu/projects/public/phm/phm_us_shp.zip', temp_shapefile)
unzip(temp_shapefile)

# Read full shapefile
shp_hardness <- read_sf('phm_us_shp.shp')

 shp_hardness_subset <- shp_hardness %>%
   filter(str_detect(ZONE, '9b|10a|10b|6a|6b|7a|7b'))

我想要的输出是这张 map ，由耐寒区域着色

 ca.az.nm <- subset(state_boundaries, STATE=="CA" | STATE=="AZ" |  STATE=="NM")
 
 ggplot() +
   geom_sf(data =  ca.az.nm) 

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

Answer 1

当您使用sf时，可用的工具之一是st_filter()，它是一种空间过滤器，允许您选择相交等几何图形。尽管取决于您的需求，但仅此一项就足够了，因为相交的耐寒多边形延伸得很远。在下面的示例中，还包括 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 年 5 月 11 日 reprex v2.0.2}