I see this all the time when developers first move from relational databases to spatial analysis: they treat a spatial join as if it's just a JOIN on a user_id, but with coordinates. They assume that if two things are "related" in the real world, a spatial join will automatically find them.

The "Perfect Overlap" Trap

The biggest misconception is believing that a spatial join finds things that are near each other. In reality, a standard spatial join (like st_join in the sf package) looks for exact geometric intersections. If a point is 0.00001 millimeters outside a polygon boundary, the join considers it non-existent.

Imagine you're analyzing a city's infrastructure. You have a layer of "Subway Station" polygons and a layer of "Coffee Shop" points. You run a join to find all coffee shops "at" the stations. You get zero results. Why? Because the coffee shop's front door is actually on the sidewalk, three feet outside the station's official property line. To the computer, they aren't touching, so they aren't joined.

# This is where most people get stuck
# It only returns points EXACTLY inside the polygon
shops_at_stations <- st_join(coffee_shops, subway_stations)
# Result: An empty data frame or many NAs.

Expanding the Search Area with Buffers

To solve this, we stop looking for intersections and start looking for proximity. This is where buffer analysis comes in. Instead of joining the raw point to the polygon, we create a "buffer zone"—an invisible bubble around our geometry—and join based on that bubble.

I always tell my juniors: buffer first, join second. If you want to find coffee shops within 200 meters of a station, you don't join the shops to the stations; you join the shops to a 200-meter polygon surrounding the stations.

# Create a 200-meter bubble around each station
station_buffers <- st_buffer(subway_stations, dist = 200)

# Now the join actually finds the "nearby" shops
shops_near_stations <- st_join(coffee_shops, station_buffers)

One quick warning: check your Coordinate Reference System (CRS). If your data is in WGS84 (longitude/latitude), dist = 200 means 200 degrees, not 200 meters. You'll end up buffering your coffee shop across three continents. Always transform to a projected coordinate system (like UTM) before buffering.

Managing the "Cartesian Explosion"

Here is a pro tip that will save your RAM: be careful with st_join when using large buffers. If a coffee shop falls within the buffer of five different subway stations, st_join will create five duplicate rows for that one coffee shop.

This is technically correct, but it can inflate your dataset unexpectedly. If you only care if a shop is near any station, I recommend using st_intersects inside a filter or mutate call rather than a full join. This keeps your row count stable while still giving you the boolean answer you need.

📋 Practical Task

Analyzing Pharmacy Accessibility in Downtown Seattle

You have been provided with two spatial data frames: seattle_pharmacies (points) and transit_hubs (points). The data is already projected in EPSG:26910 (NAD83 / UTM zone 10N), so distance units are in meters.

Your Task:

• Create a buffer of 500 meters around every transit hub.
• Perform a spatial join to identify which pharmacies fall within these 500-meter buffers.
• Calculate the percentage of total pharmacies that are "accessible" (within 500m of a hub).
• Return a final data frame containing only the pharmacies that are accessible, including the ID of the transit hub they are closest to.

Hint: Since a pharmacy might be within 500m of multiple hubs, use st_nearest_feature or a grouping method to ensure each pharmacy is only linked to its single closest hub in your final result.