By now, you've got the basics of the sf package down. But knowing the functions isn't the same as building a report that actually tells a story. For this lesson, I want us to tackle a real-world scenario: identifying "service gaps" in a city. Specifically, we're going to map out public libraries and see which neighborhoods are falling outside a walkable distance (let's say 1km) from a branch.

Turning a CSV into a Spatial Object

Usually, your data starts as a boring CSV. I've got a dataset of library locations with longitude and latitude columns. The first thing I do is pull that in and tell R that these aren't just numbers—they are coordinates on a sphere.

library(sf)
library(ggplot2)

# Loading the raw data
libs_raw <- read.csv("libraries.csv")

# Converting to an sf object
# I'm assuming WGS84 (EPSG: 4326) because that's the standard for GPS
libs_sf <- st_as_sf(libs_raw, coords = c("longitude", "latitude"), crs = 4326)

It's a common habit to stop here, but we have a problem. We can't accurately calculate "1 kilometer" using degrees of longitude and latitude. If you've ever looked at a map near the poles, you know degrees don't translate to a constant distance.

The CRS Headache

Here is where I usually trip up when I'm rushing. I'll try to create a "walkability buffer" around each library using the raw coordinates. Let's see what happens when I do that.

# My mistake: trying to buffer in degrees
walk_zones <- st_buffer(libs_sf, dist = 1000) 

# Plotting this... and wait, why is the whole world covered?
ggplot() + 
  geom_sf(data = walk_zones, fill = "lightblue", alpha = 0.5)

The map looks like a disaster. Why? Because st_buffer saw dist = 1000 and thought I meant 1,000 degrees, not 1,000 meters. I've essentially buffered my libraries by several times the circumference of the Earth. I've seen juniors spend hours debugging this; the fix is simply transforming the Coordinate Reference System (CRS) to a projected system that uses meters.

Correcting the Projection

I need to move from a geographic system (degrees) to a projected system (meters). Depending on your city, you'd pick a specific UTM zone. For this example, I'll use a generic projection that works for my region.

# Transform to a projected CRS (e.g., EPSG: 3857 for Web Mercator or a local UTM)
libs_projected <- st_transform(libs_sf, crs = 3857)

# Now the buffer actually makes sense
walk_zones_fixed <- st_buffer(libs_projected, dist = 1000)

# Let's check the plot again
ggplot() + 
  geom_sf(data = walk_zones_fixed, fill = "lightblue", alpha = 0.5) + 
  geom_sf(data = libs_projected, color = "red", size = 1) +
  theme_minimal() +
  labs(title = "1km Walkability Catchment for City Libraries")

Finding the Gaps

To make this a "report" and not just a map, we need to find the holes. I'll bring in a polygon of the city boundaries and subtract the walk zones from it. This "difference" operation is where sf really shines.

city_boundary <- st_read("city_limits.shp") |> st_transform(crs = 3857)

# Find areas NOT covered by the library buffers
underserved_areas <- st_difference(city_boundary, st_union(walk_zones_fixed))

# Final visual for the report
ggplot() + 
  geom_sf(data = city_boundary, fill = "grey90") + 
  geom_sf(data = underserved_areas, fill = "salmon", alpha = 0.7) + 
  geom_sf(data = libs_projected, color = "blue", size = 0.5) + 
  theme_void() + 
  labs(subtitle = "Highlighted areas are more than 1km from a library")

That's the workflow: Load, Project, Analyze, and Subtract. Always double-check your CRS before running st_buffer or st_distance, or you'll end up mapping the entire galaxy by accident.

📋 Practical Task

Exercise: Mapping Urban Coffee Shop Clusters in Austin

You have been given a dataset austin_coffee.csv containing the latitude and longitude of all registered coffee shops in Austin, Texas, and a shapefile austin_neighborhoods.shp.

Your goal is to create a visualization that identifies "Coffee Deserts." Complete the following steps in your R script:

• Convert the CSV data into an sf object using the WGS84 coordinate system.
• Transform the coffee shop data into a projected coordinate system (use EPSG: 3081 for Texas Centric Albers) to ensure distance calculations are in meters.
• Create a 500-meter buffer around every coffee shop.
• Using the neighborhood shapefile (also transformed to EPSG: 3081), calculate the total area of each neighborhood that is not covered by any coffee shop buffer.
• Create a final ggplot2 map where neighborhoods are colored by the size of their "coffee desert" area.