如何在R中叠加绘制tigris与osmdata的地图？

我在R中尝试把osmdata包获取的高速公路图层叠加在tigris包的行政边界图层上，但两者无法正确显示在同一图表里。我已经确认两个图层的CRS都是NAD83，排除了坐标系不匹配的问题，以下是示例代码及效果图：

library(tidyverse)
# library(censusapi)
library(tigris)
library(sf)
library(ggmap)
library(osmdata)

# Define bounding box

vec_bbox <- c(
  xmin = -71.4, 
  ymin = 42.3,
  xmax = -71.0, # Longitude range
  ymax = 42.55  # Latitude range
)

res_st_bbox <- st_bbox(
  vec_bbox, crs = "NAD83"
) 

res_sfc_bbox <- st_as_sfc(res_st_bbox)

# Highways
osm_data_highways <- opq(bbox = vec_bbox) |> # "Boston, MA") %>%
  add_osm_feature(key = "highway",
                  value = c("motorway") # , "trunk", "primary"
                  ) %>%
  osmdata_sf()

highways_sf <- osm_data_highways$osm_lines

highways_sf <- st_transform(highways_sf, crs = "NAD83")

# Get Town Data

ma_towns <- county_subdivisions(state = "MA",
                                county = "Middlesex",
                                year = 2019)

ma_towns_zoomed <- ma_towns[st_intersects(ma_towns, res_sfc_bbox, sparse = FALSE), ]


ggplot() +
  geom_sf(data = ma_towns_zoomed,color = "black",
          fill = NA) +
  geom_sf(data = highways_sf, color = "blue") +
  theme_bw()

内容的提问来源于stack exchange，提问作者Ryan Baxter-King