Subset a shapefile by lat/long

I have a shapefile from this site that I am working with. I'd like to only plot the a) 30 fathom layer, and b) this layer on latitudes above 34.4. I've been able to subset only the 30 fathom layer, but cannot subset within each region to only include certain latitudes. Below is my current code:

packages <- c('rnaturalearth', 'sf', 'tidyverse', 'ggsignif', 
              'data.table', 'patchwork', 'ggOceanMaps',
               "tools", "dplyr", 'sp', 'readr','mapproj',
              'cowplot', 'raster')

installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

invisible(lapply(packages, library, character.only = TRUE))


path.RCA <- ("~/Downloads/Rockfish_Conservation_Area_-_R7_-_CDFW_[ds3144]/")
fnam.RCA <- "Rockfish_Conservation_Area_-_R7_-_CDFW_[ds3144].shp"
RCA <- st_read(dsn = path.RCA, layer = file_path_sans_ext(fnam.RCA))

dat.RCA <- fortify(RCA)

fathom_30 <- dat.RCA[dat.RCA$area_name %like% "30-fm", ]
fathom_30 <- fathom_30[fathom_30$area_name %like% "Coastwide", ]

world <- ne_countries(scale = "medium", returnclass = "SF")

socal_map <- ggplot() +
  geom_sf(data = world) +
  geom_sf(data = fathom_30, size = 1, color = 'lightgreen')+
  coord_sf(xlim = c(-117, -121.5), ylim = c(31.75, 36), expand = FALSE)+
  scale_x_continuous(breaks = c(seq(from = -121, to = -117, by = 2))) +
  scale_y_continuous(breaks = c(seq(from = 32, to = 36, by = 2)))+
  theme(legend.position = 'bottom',
        strip.background = element_blank(), 
        legend.title = element_blank(),
        strip.text.x = element_text(
          size = 12, hjust = 0, face = "bold"),
        panel.border = element_rect(color = "red", 
                                    fill = NA, 
                                    linewidth = 1),
        axis.line.y = element_blank(), axis.line.x = element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank())

socal_map

Thank you in advance!

As Chris noted in comments, you need to take care of coordinate reference system (CRS). So if you need to crop sf object by a bounding box or some other sf object, both must share the same CRS. And to transform from one to another, we use st_transform(). So you either need to transform your cropping boundaries to EPSG:3857 ((-124.679 34.4) translates to (-13879203 4082640) ) or transform sf object to geographic coordinates, WGS84. We'll go with the latter.

Cropping with geographic coordinates can be bit tricky in sf when s2 is enabled (geometries on a sphere, on by default), so we'll turn it off while manipulating with geometries. One of recent answers that attempts to explain & illustrate this - https://stackoverflow.com/a/78725088/646761

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(dplyr)
library(stringr)
library(ggplot2)

world <-  rnaturalearth::ne_countries(scale = "medium")

# underlying GDAL library supports virtual file systems and can read
# zipped shapefiles by prefixing file name with `/vsizip/`
rca_shp <- "/vsizip/Rockfish_Conservation_Area_-_R7_-_CDFW_[ds3144].zip"

# check layers and projection / CRS
st_layers(rca_shp)
#> Driver: ESRI Shapefile 
#> Available layers:
#>                                        layer_name geometry_type features fields
#> 1 Rockfish_Conservation_Area_-_R7_-_CDFW_[ds3144]   Line String       71      5
#>                   crs_name
#> 1 WGS 84 / Pseudo-Mercator

# load shapefile, transform from pseudo-mercator (EPSG:3857) to 
# geographic coordinates (WGS84)
rca_wgs84 <- 
  read_sf(rca_shp) |> 
  st_transform("wgs84")

# use bounding box of sf object for a crop area template, update latitude limits 
(bbox <- st_bbox(rca_wgs84)) 
#>       xmin       ymin       xmax       ymax 
#> -124.67900   32.35717 -117.25733   42.00000
bbox[["ymin"]] <- 34.4
bbox[["ymax"]] <- 36

# to crop with a bounding box and handle geographic coordinates as 
# planar coordinates, we'll crop with s2 (geometries on a sphere) disabled; 
# filter by area_name, crop by updated bbox
sf_use_s2(use_s2 = FALSE)
#> Spherical geometry (s2) switched off
fathom_30 <- 
  rca_wgs84 |> 
  filter(str_detect(area_name, "30-fm.*Coastwide")) |> 
  st_crop(bbox) 
#> although coordinates are longitude/latitude, st_intersection assumes that they
#> are planar
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
sf_use_s2(use_s2 = TRUE)
#> Spherical geometry (s2) switched on

# after filtering by area_name and cropping into bbox we are left with 
# 2 features: 
fathom_30
#> Simple feature collection with 2 features and 5 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -121.5683 ymin: 34.4 xmax: -119.9203 ymax: 36
#> Geodetic CRS:  WGS 84
#> # A tibble: 2 × 6
#>   OBJECTID area_name      START_Fath Region Shape__Len                  geometry
#> *    <int> <chr>               <int> <chr>       <dbl>          <LINESTRING [°]>
#> 1       52 30-fm (55-m) …         30 South…    548019. (-120.5052 34.45, -120.4…
#> 2       53 30-fm (55-m) …         30 Centr…    266038. (-121.5683 36, -121.5635…

# to visualize
# mapview::mapview(fathom_30, zcol = "Region")

# plot
ggplot() +
  geom_sf(data = world) +
  geom_sf(data = fathom_30, linewidth = 1, aes(color = "fathom_30")) +
  scale_color_manual(name = "Legend title", values = c(fathom_30 = "lightgreen")) +
  coord_sf(xlim = c(-117, -121.5), ylim = c(31.75, 36), expand = FALSE) +
  theme_minimal()

^{Created on 2024-07-23 with reprex v2.1.0}