This guide walks you through making maps in R — from your first simple map to interactive visualizations. No prior GIS experience needed, just some basic familiarity with R.

1. Key Concepts

Before touching any code, here are a few terms that will come up everywhere.

Spatial data describes things that have a location. A city, a county boundary, a road, a hospital — all of these can be represented as spatial data.

Vector data represents locations as points, lines, or polygons. A city might be a point; a river might be a line; a state boundary is a polygon.

Raster data represents space as a grid of cells (pixels), like satellite imagery or elevation data. This guide focuses on vector data, which is what most social-science mapping involves.

CRS (Coordinate Reference System) is the system that defines how coordinates map onto the Earth’s surface. Two datasets with different CRS values won’t line up on a map. You’ll need to match them using st_transform(). A common one is WGS84, which uses longitude and latitude (EPSG code 4326). CRS can be a surprisingly deep topic — for a more thorough explanation, Reed’s GIS team has put together a helpful overview at https://www.reed.edu/data-at-reed/gis/crs.html.

sf (Simple Features) is the modern standard for working with vector spatial data in R. An sf object looks like a regular data frame, but it has an extra column (usually called geometry) that stores the spatial information.

2. Packages You’ll Need

If you have never used these packages before, you will need to install them first. If you are working on the Reed R server, most packages are already installed for you — however, you will still need to install rnaturalearth, rnaturalearthdata, and ggspatial yourself. If R tells you “there is no package called ‘…’”, that is your signal that something still needs to be installed.

Run the line below once in your console (it is set to eval=FALSE so it does not run every time you knit):

install.packages(c("sf", "ggplot2", "tmap", "leaflet", "tidycensus",
                   "rnaturalearth", "rnaturalearthdata", "ggspatial",
                   "dplyr", "viridis"))

What each package does (best practice is to write package names with curly braces like {sf} to make it clear you are referring to a package):

  • {sf} — reads, writes, and manipulates spatial data
  • {ggplot2} — the main plotting library; handles sf objects with geom_sf()
  • {tmap} — a mapping-focused package; can produce static and interactive maps using the same code
  • {leaflet} — makes interactive, web-style maps (like Google Maps)
  • {tidycensus} — downloads US Census and ACS data with geometry already attached
  • {rnaturalearth} — provides public domain world and country boundary data
  • {ggspatial} — adds scale bars and north arrows to ggplot2 maps
  • {viridis} — colorblind-friendly color palettes

3. Getting Spatial Data into R

Option A: Use built-in or package data

The {rnaturalearth} package is a great starting point. It gives you country and continent boundaries from public domain data. It works alongside {rnaturalearthdata}, which is a companion package that actually stores the boundary data — {rnaturalearth} provides the functions to access it, while {rnaturalearthdata} provides the files themselves. You need both loaded.

library(sf)
library(rnaturalearth)
library(rnaturalearthdata)

# Get world country boundaries as an sf object
world <- ne_countries(scale = "medium", returnclass = "sf")

# Take a look at the data structure
# (use names(world) to see all available columns)
head(world[, 1:5])
## Simple feature collection with 6 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.36621 ymin: -22.40205 xmax: 109.4449 ymax: 41.9062
## Geodetic CRS:  WGS 84
##        featurecla scalerank labelrank sovereignt sov_a3
## 1 Admin-0 country         1         3   Zimbabwe    ZWE
## 2 Admin-0 country         1         3     Zambia    ZMB
## 3 Admin-0 country         1         3      Yemen    YEM
## 4 Admin-0 country         3         2    Vietnam    VNM
## 5 Admin-0 country         5         3  Venezuela    VEN
## 6 Admin-0 country         6         6    Vatican    VAT
##                         geometry
## 1 MULTIPOLYGON (((31.28789 -2...
## 2 MULTIPOLYGON (((30.39609 -1...
## 3 MULTIPOLYGON (((53.08564 16...
## 4 MULTIPOLYGON (((104.064 10....
## 5 MULTIPOLYGON (((-60.82119 9...
## 6 MULTIPOLYGON (((12.43916 41...
class(world)  # should say "sf" and "data.frame"
## [1] "sf"         "data.frame"

Option B: Read a shapefile from disk

If you have a .shp file downloaded from the internet or a government website, read it with st_read(). Shapefiles are one of the most common formats for spatial data — for an overview of spatial file types you might encounter, look out for the upcoming Data at Reed guide on file formats.

my_map <- st_read("path/to/your/file.shp")

Option C: Convert a regular data frame with coordinates

If your data has latitude and longitude columns, you can convert it to an sf object:

library(sf)

df <- data.frame(
  city = c("Portland", "Seattle", "San Francisco"),
  lat  = c(45.52, 47.61, 37.77),
  lon  = c(-122.68, -122.33, -122.42),
  pop  = c(650000, 750000, 870000)
)

# crs = 4326 means WGS84 (standard lat/lon)
cities_sf <- st_as_sf(df, coords = c("lon", "lat"), crs = 4326)
cities_sf
## Simple feature collection with 3 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -122.68 ymin: 37.77 xmax: -122.33 ymax: 47.61
## Geodetic CRS:  WGS 84
##            city    pop              geometry
## 1      Portland 650000 POINT (-122.68 45.52)
## 2       Seattle 750000 POINT (-122.33 47.61)
## 3 San Francisco 870000 POINT (-122.42 37.77)

4. Your First Map with ggplot2 + sf

{ggplot2} maps work just like regular ggplot charts — you add layers with +. The key geometry for maps is geom_sf().

library(ggplot2)

ggplot(data = world) +
  geom_sf()

That’s it. You have a world map. Now let’s make it look better by changing the theme of the background:

ggplot(data = world) +
  geom_sf(fill = "lightgray", color = "white", linewidth = 0.3) +
  theme_minimal() +
  labs(title = "World Map")

To zoom into a specific region, use coord_sf(). You do this by specifying the minimum and maximum longitude (xlim) and latitude (ylim) of the area you want to show:

ggplot(data = world) +
  geom_sf(fill = "lightgray", color = "white") +
  coord_sf(xlim = c(-20, 55), ylim = c(-40, 40)) +
  theme_minimal() +
  labs(title = "Africa")

5. Choropleth Maps

A choropleth map uses color to show how a variable varies across geographic areas — think election maps that have red and blue states, or population density maps with darker colors for denser areas.

The world object we loaded from {rnaturalearth} comes with several built-in variables we can map. Here we will use GDP — the total economic output of each country, stored in a column whose name varies by package version (hence the detection step below). We plot it on a log scale because GDP values span several orders of magnitude, which would otherwise make most countries look identical.

The map below uses the “plasma” color palette from {viridis}: purple indicates lower GDP and yellow indicates higher GDP.

library(viridis)

# Detect GDP column name — the name varies by rnaturalearth version
# names() returns all column names in the data frame
# intersect() finds which of our candidates actually exist, and [1] picks the first match
gdp_col <- intersect(c("gdp_md_est", "GDP_MD", "gdp_md"), names(world))[1]

ggplot(data = world) +
  # .data[[gdp_col]] is how you reference a column by a variable name inside aes()
  geom_sf(aes(fill = .data[[gdp_col]])) +
  scale_fill_viridis_c(
    name   = "GDP (millions USD)",
    trans  = "log",     # log scale: compresses the range so differences are visible
    option = "plasma"   # "plasma" is one of several viridis palettes; try "viridis" or "magma"
  ) +
  theme_minimal() +
  labs(title = "GDP by Country")

The key part is aes(fill = .data[[gdp_col]]). The .data[[...]] syntax is the tidyverse way of saying “look up this column name from a variable” — useful whenever the column name is stored as a string rather than typed directly.

Color palette options:

  • scale_fill_viridis_c() — continuous data, colorblind-safe
  • scale_fill_distiller(palette = "Blues") — sequential ColorBrewer palettes
  • scale_fill_gradient(low = "white", high = "darkblue") — simple two-color gradient
  • scale_fill_brewer(palette = "RdYlGn") — for categorical/factor data (no _c)

6. Adding Map Elements

A map is more useful with a scale bar, north arrow, and proper labels. The {ggspatial} package adds these to ggplot2 maps.

library(ggspatial)
library(dplyr)

africa <- world %>%
  filter(continent == "Africa")

ggplot(data = africa) +
  geom_sf(fill = "lightblue", color = "gray40") +
  # annotation_scale adds a scale bar
  # location = "bl" puts it in the bottom-left; width_hint controls how wide it is (0–1)
  annotation_scale(location = "bl", width_hint = 0.3) +
  # annotation_north_arrow adds a north arrow
  # which_north = "true" points to true north (vs. grid north)
  # north_arrow_fancy_orienteering() is a decorative style; try north_arrow_minimal() for simpler
  annotation_north_arrow(
    location    = "tr",
    which_north = "true",
    style       = north_arrow_fancy_orienteering()
  ) +
  labs(
    title    = "Africa",
    subtitle = "Country boundaries",
    caption  = "Source: Natural Earth"
  ) +
  theme_minimal() +
  theme(
    axis.text  = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank()
  )

To explore all the options available for these functions, run ?annotation_scale or ?annotation_north_arrow in your console, or visit the {ggspatial} documentation.

7. Interactive Maps with leaflet

{leaflet} creates maps you can zoom, pan, and click on.

Note: Run these chunks interactively in RStudio (Viewer pane) or knit to HTML. The chunks below use eval=FALSE so the document knits cleanly — remove it when you want to run them.

Basic interactive map

library(leaflet)

leaflet() %>%
  addTiles() %>%                               # default OpenStreetMap basemap
  setView(lng = -98.5, lat = 39.5, zoom = 4)  # center on the US

Add markers for point locations

library(dplyr)

cities <- data.frame(
  name = c("Portland", "Seattle", "San Francisco"),
  lat  = c(45.52, 47.61, 37.77),
  lon  = c(-122.68, -122.33, -122.42),
  pop  = c(650000, 750000, 870000)
)

leaflet(data = cities) %>%
  addTiles() %>%
  addMarkers(
    lng   = ~lon,
    lat   = ~lat,
    popup = ~paste0("<b>", name, "</b><br>Population: ", format(pop, big.mark = ","))
  )

The ~ before column names tells leaflet to look inside the data frame. popup controls what appears when you click a marker (HTML is supported).

Add polygons

So far we have only added point markers. Polygons are shapes — like country or county boundaries — that cover an area rather than mark a single location. Adding polygon layers lets you display geographic boundaries as interactive shapes that users can hover over or click.

# leaflet requires the CRS to be WGS84 (EPSG 4326) — reproject if needed
africa_ll <- st_transform(africa, 4326)

leaflet(data = africa_ll) %>%
  addTiles() %>%
  addPolygons(
    fillColor   = "steelblue",
    weight      = 1,           # border line thickness
    color       = "white",     # border color
    fillOpacity = 0.7,
    popup       = ~name_long   # show country name on click
  )

Choropleth in leaflet

Now we will color each country’s polygon by its population, using the pop_est column. The popup argument builds a small HTML label that appears when you click a country, showing its name and formatted population.

pal <- colorNumeric(
  palette = "YlOrRd",
  domain  = africa_ll$pop_est
)

leaflet(data = africa_ll) %>%
  addTiles() %>%
  addPolygons(
    fillColor   = ~pal(pop_est),
    weight      = 1,
    color       = "white",
    fillOpacity = 0.8,
    # popup builds an HTML string shown on click: country name + formatted population
    popup       = ~paste0(name_long, "<br>Population: ",
                          format(pop_est, big.mark = ","))
  ) %>%
  addLegend(
    pal      = pal,
    values   = ~pop_est,
    title    = "Population",
    position = "bottomright"
  )

8. Interactive Maps with tmap

{tmap} is built specifically for maps. Its biggest advantage: switch between static and interactive output with one line — same plotting code, different mode.

library(tmap)

tmap_mode("plot") # static mode

tm_shape(world) +
  tm_polygons(
    col     = "pop_est",
    title   = "Population",
    style   = "quantile",
    palette = "Blues"
  ) +
  tm_layout(main.title = "World Population")

Switch to interactive mode:

tmap_mode("view") # interactive mode — same code, different output

tm_shape(world) +
  tm_polygons(
    col     = "pop_est",
    title   = "Population",
    style   = "quantile",
    palette = "Blues"
  )

Classification styles for the style argument. “Breaks” refers to the cut points that divide your data into color groups — different methods produce different groupings:

Style Description
"quantile" Equal-count breaks — each group has the same number of countries (good default)
"equal" Equal-width breaks — each group spans the same numeric range
"jenks" Natural breaks — minimizes differences within groups, maximizes differences between them
"cont" Continuous gradient — no discrete breaks at all
"fixed" Set your own breaks manually

9. Mapping Census Data with tidycensus

{tidycensus} is the easiest way to get US Census or ACS data with geometry already attached. You need a free Census API key: https://api.census.gov/data/key_signup.html

You will receive your key in an email. Register it once in R (run in your console, not in the document):

library(tidycensus)
census_api_key("YOUR_KEY_HERE", install = TRUE)

Download ACS data with geometry

ACS stands for American Community Survey — an ongoing Census Bureau survey that provides detailed demographic and economic estimates for states, counties, and smaller geographies. Variables in the ACS have coded names like B19013_001 (explained below). Setting geometry = TRUE automatically attaches the county shapefile so you can map the results without any extra steps.

library(tidycensus)

# Get median household income by county in Oregon
# geometry = TRUE automatically attaches the county shapefile
or_income <- get_acs(
  geography = "county",
  variables = "B19013_001",  # median household income
  state     = "OR",
  year      = 2022,
  geometry  = TRUE
)

head(or_income)

The returned data frame has these columns:

  • GEOID — a unique numeric identifier for each geographic unit (e.g. county FIPS code)
  • NAME — a human-readable place name (e.g. “Multnomah County, Oregon”)
  • variable — the ACS variable code you requested
  • estimate — the estimated value (e.g. median income in dollars)
  • moe — margin of error for the estimate
  • geometry — the spatial boundary, used for mapping

Map it

ggplot(data = or_income) +
  geom_sf(aes(fill = estimate)) +
  scale_fill_viridis_c(
    name   = "Median income ($)",
    labels = scales::dollar
  ) +
  labs(
    title    = "Median Household Income by County",
    subtitle = "Oregon, 2022 ACS 5-Year Estimates",
    caption  = "Source: US Census Bureau"
  ) +
  theme_minimal() +
  theme(axis.text = element_blank(), axis.ticks = element_blank())

Note: Because the Census chunks require an API key and internet access, they use eval=FALSE and will not produce output when knitting. Once you have a key registered, remove eval=FALSE from the chunk options and they will run normally.

Find variable codes

Use load_variables() to search for what you want:

vars <- load_variables(2022, "acs5", cache = TRUE)

# Search for variables related to education
vars %>%
  filter(grepl("education|school", label, ignore.case = TRUE)) %>%
  head(20)

Common ACS variables

Variable Description
B01003_001 Total population
B19013_001 Median household income
B25077_001 Median home value
B15003_022 Bachelor’s degree (25+)
B23025_005 Unemployed in labor force
B03002_003 White alone, non-Hispanic

National map with shifted Alaska and Hawaii

For US-wide maps, use shift_geo = TRUE to reposition Alaska and Hawaii below the continental US:

state_age <- get_acs(
  geography = "state",
  variables = "B01002_001",  # median age
  year      = 2022,
  geometry  = TRUE,
  shift_geo = TRUE           # shift AK and HI
)

ggplot(data = state_age) +
  geom_sf(aes(fill = estimate)) +
  scale_fill_viridis_c(name = "Median age") +
  labs(title = "Median Age by State, 2022") +
  theme_void()

10. Point Maps

Point maps show individual locations as dots. Dot size and color can encode data.

us <- ne_countries(
  country     = "united states of america",
  scale       = "medium",
  returnclass = "sf"
)

cities <- data.frame(
  city = c("New York", "Los Angeles", "Chicago", "Houston", "Phoenix"),
  lon  = c(-74.01, -118.24, -87.63, -95.37, -112.07),
  lat  = c( 40.71,   34.05,  41.88,  29.76,   33.45),
  pop  = c(8.3, 3.9, 2.7, 2.3, 1.6)  # in millions
)

cities_sf <- st_as_sf(cities, coords = c("lon", "lat"), crs = 4326)

ggplot() +
  geom_sf(data = us, fill = "gray95", color = "gray70") +
  geom_sf(data = cities_sf, aes(size = pop), color = "steelblue", alpha = 0.7) +
  scale_size_continuous(name = "Population (millions)", range = c(3, 12)) +
  geom_sf_text(data = cities_sf, aes(label = city), nudge_y = 1.5, size = 3) +
  labs(title = "Major US Cities by Population") +
  theme_void()

Cropping to a specific area

Sometimes you want to focus on a particular region without filtering your data. You can use coord_sf() to zoom in by specifying longitude (xlim) and latitude (ylim) limits. This keeps all the data intact but only displays the area you care about. The tradeoff is that you are making a deliberate choice to exclude part of the country for visual clarity — always consider whether that is appropriate for your purpose.

ggplot() +
  geom_sf(data = us, fill = "gray95", color = "gray70") +
  geom_sf(data = cities_sf, aes(size = pop), color = "steelblue", alpha = 0.7) +
  scale_size_continuous(name = "Population (millions)", range = c(3, 12)) +
  geom_sf_text(data = cities_sf, aes(label = city), nudge_y = 1.5, size = 3) +
  coord_sf(xlim = c(-125, -66), ylim = c(24, 50)) +  # crop to contiguous US
  labs(title = "Major US Cities by Population (Contiguous US)") +
  theme_void()

11. Saving Your Maps

These chunks use eval=FALSE since you will run them interactively when you want to save output files.

Static maps (ggplot2)

my_map <- ggplot(data = world) + geom_sf()

ggsave("my_map.png", plot = my_map, width = 8, height = 6, dpi = 300)
ggsave("my_map.pdf", plot = my_map, width = 8, height = 6)

Interactive maps (leaflet or tmap)

Note: Interactive maps only work inside HTML documents. If you save to PDF or Word, the map will not appear. Use saveWidget() or tmap_save() to export as a standalone .html file that anyone can open in a browser.

library(htmlwidgets)

# leaflet
my_leaflet_map <- leaflet() %>% addTiles()
saveWidget(my_leaflet_map, "my_map.html")

# tmap interactive
tmap_mode("view")
my_tmap <- tm_shape(world) + tm_polygons()
tmap_save(my_tmap, "my_map.html")

# tmap static
tmap_mode("plot")
tmap_save(my_tmap, "my_map.png", dpi = 300)

12. Common Errors and Fixes

“Error: Column geometry does not exist”
Your data frame has not been converted to an sf object. Use st_as_sf() or make sure you set geometry = TRUE in get_acs().

“although coordinates are longitude/latitude, st_intersects assumes that they are planar”
This is a warning, not an error. It means your CRS uses degrees instead of meters. For most visualizations this is fine to ignore. To silence it, reproject: st_transform(your_sf, 3857).

Map layers don’t line up
Your datasets have different CRS values. Use st_crs(your_sf) to check, then st_transform(your_sf, 4326) to reproject to WGS84.

Leaflet map is blank or doesn’t render
Leaflet requires WGS84 (EPSG 4326). Always run st_transform(your_sf, 4326) before passing data to leaflet.

Map renders very slowly
This usually means your shapefile has very high geometric resolution — it is storing thousands of tiny coordinate points to draw each boundary precisely. For most thematic maps you do not need that level of detail. The {rmapshaper} package can simplify the geometry by removing redundant points, which dramatically speeds up rendering without meaningfully changing the appearance at normal zoom levels:

library(rmapshaper)
# keep = 0.05 retains 5% of vertices; increase toward 1 for more detail
simplified <- ms_simplify(your_sf, keep = 0.05)

get_acs() returns an error about an API key
Make sure you registered your key with census_api_key("YOUR_KEY"). Check with Sys.getenv("CENSUS_API_KEY").

13. Further Resources

Books (free online)

Package documentation

Data sources

  • US Census/ACS data: {tidycensus} package
  • World/country boundaries: {rnaturalearth} package
  • US state and county shapefiles: {tigris} package
  • OpenStreetMap-based data: {osmdata} package

14. Working with Multiple Layers and Filtering Spatial Data

One of the most powerful things you can do in a GIS program — and in R — is stack multiple spatial layers on top of each other and filter them to show only the features you care about. This section translates those ideas from QGIS into R.

Layering multiple shapefiles

In QGIS, you load each dataset as a separate layer and stack them in a panel. In R with {ggplot2}, you do the same thing by adding multiple geom_sf() calls — each one is a new layer drawn on top of the previous one. The order matters: layers listed first are drawn first and can be covered by later ones.

The alpha argument controls opacity, from 0 (fully transparent) to 1 (fully opaque). This is the R equivalent of the opacity slider in QGIS, and it is useful when layers overlap and you want to see through them.

The example below recreates the kind of layered Portland map from the QGIS intro tutorial — a city boundary, roads, and rivers — using data from the {tigris} package, which provides US Census boundary files directly in R (no downloading needed).

# Install {tigris} if you haven't already
install.packages("tigris")
library(tigris)
library(dplyr)
library(ggplot2)
library(sf)

options(tigris_use_cache = TRUE)

# Portland city boundary
portland <- places(state = "OR", cb = TRUE) %>%
  filter(NAME == "Portland")

# Roads in Multnomah County (loading all OR roads would be very slow)
roads <- roads(state = "OR", county = "Multnomah")

# Rivers/water bodies in Multnomah County
rivers <- area_water(state = "OR", county = "Multnomah")

# Make sure all layers share the same CRS
roads  <- st_transform(roads,  st_crs(portland))
rivers <- st_transform(rivers, st_crs(portland))

# Clip roads and rivers to the Portland city boundary
# st_intersection() keeps only the parts of each layer that fall inside Portland
roads_pdx  <- st_intersection(roads,  portland)
rivers_pdx <- st_intersection(rivers, portland)
ggplot() +
  # Layer 1: city boundary — drawn first, sits at the bottom
  geom_sf(data = portland, fill = "lightyellow", color = "black", linewidth = 0.6) +
  # Layer 2: roads — alpha = 0.5 makes them semi-transparent
  geom_sf(data = roads_pdx, color = "gray50", linewidth = 0.2, alpha = 0.5) +
  # Layer 3: rivers — blue fill, semi-transparent so roads show through
  geom_sf(data = rivers_pdx, fill = "steelblue", color = NA, alpha = 0.6) +
  labs(
    title   = "Portland: Roads and Water",
    caption = "Source: US Census Bureau TIGER/Line"
  ) +
  theme_void()

Filtering features with dplyr

In QGIS’s vector data tools, you can select features by their attributes — for example, showing only highways, or only a specific type of road. In R, this is just a filter() call on the data frame before you pass it to ggplot(). Because sf objects are data frames, all standard {dplyr} verbs work on them directly.

# MTFCC is the Census road classification code
# S1100 = primary highways, S1200 = secondary roads
# Filtering to just these removes local streets and reduces clutter
major_roads <- roads_pdx %>%
  filter(MTFCC %in% c("S1100", "S1200"))

ggplot() +
  geom_sf(data = portland, fill = "lightyellow", color = "black", linewidth = 0.6) +
  geom_sf(data = major_roads, color = "gray30", linewidth = 0.5) +
  geom_sf(data = rivers_pdx, fill = "steelblue", color = NA, alpha = 0.7) +
  labs(
    title    = "Portland: Major Roads and Water",
    subtitle = "Primary and secondary roads only",
    caption  = "Source: US Census Bureau TIGER/Line"
  ) +
  theme_void()

Selecting features by location

In QGIS, “select by location” lets you keep only features from one layer that spatially overlap with another — for example, keeping only the roads that fall inside a city boundary. In R this is done with st_intersection() from {sf}, which clips one layer to the shape of another. We already used this above to clip roads and rivers to Portland. Here is a more explicit walkthrough with county boundaries:

multnomah <- counties(state = "OR", cb = TRUE) %>%
  filter(NAME == "Multnomah") %>%
  st_transform(st_crs(roads))

# st_intersection() returns only the parts of `roads` that overlap
# with `multnomah` — anything outside the county is clipped off
roads_multnomah <- st_intersection(roads, multnomah)
ggplot() +
  geom_sf(data = multnomah, fill = "gray95", color = "black") +
  geom_sf(data = roads_multnomah, color = "gray40", linewidth = 0.2, alpha = 0.6) +
  labs(
    title   = "Roads in Multnomah County",
    caption = "Source: US Census Bureau TIGER/Line"
  ) +
  theme_void()

st_intersection() is the workhorse for spatial filtering, but {sf} has a full set of spatial relationship functions worth knowing:

Function What it does
st_intersection(x, y) Keeps parts of x that overlap y, clipped to that shape
st_filter(x, y) Keeps full features from x that touch y, without clipping
st_join(x, y) Joins attributes from y onto features in x that overlap
st_within(x, y) Returns TRUE/FALSE: is each feature in x fully inside y?
st_intersects(x, y) Returns TRUE/FALSE: does each feature in x touch y at all?

For example, to keep whole road segments that intersect Portland without clipping them at the boundary edge:

# st_filter keeps the full geometry of matching features (no clipping)
roads_within <- st_filter(roads, portland)