Spatial Subsetting

library(sf)

# Set up paths to data
ard_path <- "data/AnalysisReady/AnalysisReady.gpkg"
ws_path <- "data/source/watersheds/WBDHU12.shp"
segi_path <- "data/source/segi/2023_GS_Groves_OTI_public.gdb"

# Transform watershed data
# open watershed shapefile
ws <- st_read(ws_path)
# View layers in analysis ready gpkg
st_layers(ard_path)
# Open the fire boundary layer
fire <- st_read(ard_path, layer = "fire_boundary")
# Are the CRSs the same?
st_crs(ws) == st_crs(fire)
# Reproject ws to match the fire dataset
ws <- st_transform(ws, crs = st_crs(fire))

# Spatially subset watersheds
# Create a new object that containing the watersheds that overlap with the fire boundary.
ws_fire <- ws[fire, ]
# Save the result in our analysis ready geopackage
st_write(ws_fire, ard_path, layer = "watersheds")

# Transform giant sequoia (segi) points data
st_layers(segi_path)
segi_pts <- st_read(segi_path, layer = "SEGI_OTI_2023_public")
st_crs(segi_pts) == st_crs(fire)

# Get the segi locations that are in the fire perimeter.
# We'll use the piping operator (%>%) to chain together the workflow.
# Ctrl + Shift + m is the shortcut key to %>%
segi_pts_fire <- segi_pts %>% 
  st_transform(crs = st_crs(fire)) %>% 
  st_filter(y = fire, .predicate = st_intersects) %>% 
  st_write(ard_path, layer = "segi_points")

# Transform SEGI polygons in a similar way
st_layers(segi_path)
segi_poly <- st_read(segi_path, layer = "OGS_Groves_2023_97_public")
segi_poly <- st_cast(segi_poly, "MULTIPOLYGON")

segi_poly_fire <- segi_poly %>% 
  st_transform(crs = st_crs(fire)) %>% 
  st_filter(y = fire) %>% 
  st_write(ard_path, layer = "segi_groves")

# View the layers in the geopackage
st_layers(ard_path)

Additional information on geometric operations with sf: https://r-spatial.github.io/sf/reference/geos_binary_ops.html