基于另一个 sf 特征绘制 sf 对象的子集（对于某些状态）

Question

这个问题与how map certain USDA harness zones in R 有关，但我只需要一张 CA、NM 和 AZ 的地图。我要在

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

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

 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_crop

）来做到这一点，因为存在重复的边。我必须将 shapefile 转换为

SpatialPolygonsDataFrame

，剪辑它们，然后再转换回

sf

.

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

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

shp_hardness_subset_sp <- as_Spatial(shp_hardness_subset)
ca.az.nm_sp <- as_Spatial(ca.az.nm)

st_crs(ca.az.nm_sp)==st_crs(shp_hardness_subset_sp)
shp_hardness_subset_sp_crop <- crop(shp_hardness_subset_sp, ca.az.nm_sp)

shp_hardness_subset_sf_crop <- st_as_sf(shp_hardness_subset_sp_crop)

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

Answer 2

当您使用

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-05-11 与 reprex v2.0.2}