I’m having trouble plotting data from the tigris and osmdata packages on the same figure. The goal is to have highways overlaying administrative boundaries.

I have included sample code below. I do not believe the issue is due to the objects having different CRS values. Thanks!

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