Search code examples
rggplot2subsetgeospatialr-sf

Plotting a subset of a sf object (for certain states) based on another sf feature

This question is related to how map certain USDA hardiness zones in R, but I need a map of CA, NM, and AZ only. What do I put into the filter command to have the 3 states only?

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

My desired output is this map, colored by the hardiness zones

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

enter image description here

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

enter image description here

Solution

  • As you are using sf, one of the available tools is st_filter(), a spatial filter that allows you to select geometries that are, for example, intersecting. Though depending on your needs, this alone might be enough as intersecting hardiness polygons extend quite far. In the example below, intersection by ca.az.nm polygons and crop by ca.az.nm bounding box are also included.

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

    Created on 2023-05-11 with reprex v2.0.2