Spatial Analysis Exercises in R: 15 Practice Problems

nc_path <- system.file("shape/nc.shp", package = "sf")
ex_1_1 <- st_read(nc_path, quiet = TRUE)
class(ex_1_1)
#> [1] "sf"         "data.frame"
dim(ex_1_1)
#> [1] 100  15

Explanation: st_read() returns an sf object that inherits from data.frame, so every dplyr verb you already know works on it. The geometry column is a list-column of class sfc carrying the polygons; the other 14 columns are the attribute table from the DBF sidecar. Passing quiet = TRUE keeps stdout clean inside scripts and Rmd reports. If st_read() cannot find the file, double-check system.file() returned a non-empty string.

ex_1_2 <- list(
  n         = nrow(ex_1_1),
  geom_type = unique(st_geometry_type(ex_1_1)),
  cols      = head(names(ex_1_1), 6)
)
ex_1_2
#> $n
#> [1] 100
#> $geom_type
#> [1] MULTIPOLYGON
#> Levels: MULTIPOLYGON
#> $cols
#> [1] "AREA"     "PERIMETER" "CNTY_"     "CNTY_ID"  "NAME"      "FIPS"

Explanation: st_geometry_type() returns one value per feature, so unique() gives you the single shared type when the layer is homogeneous. Shapefiles always store polygons as MULTIPOLYGON even when each feature is a single ring, because the format permits multi-part features. Use this audit pattern (count + geometry type + attribute preview) as the first cell of any new spatial notebook so downstream surprises (mixed geometry, missing attributes) surface immediately.

stores <- tibble::tibble(
  store = c("Raleigh", "Charlotte", "Greensboro", "Durham", "Asheville"),
  lon   = c(-78.6382, -80.8431, -79.7920, -78.8986, -82.5515),
  lat   = c( 35.7796,  35.2271,  36.0726,  35.9940,  35.5951)
)
ex_1_3 <- st_as_sf(stores, coords = c("lon", "lat"), crs = 4326)
ex_1_3
#> Simple feature collection with 5 features and 1 field
#> Geometry type: POINT
#> Geodetic CRS:  WGS 84

Explanation: st_as_sf() converts any data frame with coordinate columns to an sf object in one call. Always declare the CRS at construction time: an sf object without a CRS is a silent landmine because st_distance() will return planar units while you assume meters. EPSG 4326 is the WGS84 lon/lat code, the default for GPS-derived data. If the upstream system used EPSG 3857 (web mercator tiles) or a state plane code, pass that instead.

ex_2_1 <- st_crs(ex_1_1)$epsg
ex_2_1
#> [1] 4267

Explanation: st_crs() returns a crs object that prints a WKT block but exposes $epsg, $wkt, $proj4string, and $input accessors for programmatic use. NC's bundled shapefile is in NAD27 (EPSG 4267), a 1927 North American datum that differs from WGS84 by up to 100m, so reprojecting before any cross-dataset analysis is mandatory. When $epsg returns NA, the layer carries a custom WKT without a registered code and you must reproject by passing the full CRS object, not an integer.

ex_2_2 <- st_transform(ex_1_1, 3857)
print(ex_2_2, n = 0)
#> Simple feature collection with 100 features and 14 fields
#> Geometry type: MULTIPOLYGON
#> Bounding box:  xmin: -9531829 ymin: 3858905 xmax: -8420343 ymax: 4324619
#> Projected CRS: WGS 84 / Pseudo-Mercator

Explanation: st_transform() reprojects vector geometries using PROJ under the hood. EPSG 3857 stores coordinates in meters from the equator and the prime meridian, hence the seven-digit values. Web mercator badly distorts area away from the equator (Greenland looks the size of Africa), so never use 3857 for area or distance calculations: it is a display projection only. For numerical work, pick a local projected CRS such as UTM or a state plane (Exercise 2.3).

nc_utm <- st_transform(ex_1_1, 32617)
nc_utm$area_km2 <- units::set_units(st_area(nc_utm), "km^2")

ex_2_3 <- nc_utm |>
  dplyr::arrange(dplyr::desc(area_km2)) |>
  dplyr::select(NAME, area_km2) |>
  head(5)
ex_2_3

Explanation: st_area() returns a units-aware vector in the CRS's native units (m² for UTM). set_units() from the units package converts safely without silent factor mistakes. Using UTM zone 17N for North Carolina keeps area distortion under 0.1% across the state; using web mercator (EPSG 3857) instead would inflate northern counties by ~25%. For an irregular study area that spans multiple UTM zones, prefer an equal-area projection such as Albers (EPSG 5070 for CONUS).

ex_3_1 <- suppressWarnings(
  ex_1_1 |>
    dplyr::select(NAME, BIR74) |>
    st_centroid()
)
head(ex_3_1, 5)

Explanation: st_centroid() replaces each polygon's geometry with its mathematical centroid. The warning ("attribute variables are assumed to be spatially constant") fires because attributes like population are no longer area-aggregated after collapsing to a point, which is a sensible reminder but harmless here. For concave shapes (e.g., a crescent-shaped park) the centroid can fall outside the polygon: use st_point_on_surface() if you need a point guaranteed to lie inside.

ex_3_2 <- ex_1_3 |>
  st_transform(3857) |>
  st_buffer(dist = 30000) |>
  st_union()

length(ex_3_2)
#> [1] 1
print(ex_3_2, n = 0)
#> Geometry set for 1 feature 
#> Geometry type: MULTIPOLYGON
#> Projected CRS: WGS 84 / Pseudo-Mercator

Explanation: st_buffer() requires a projected CRS so the distance argument is in linear meters; calling it on a geodetic (lon/lat) layer would silently treat the distance as decimal degrees and produce wildly wrong shapes. st_union() dissolves overlapping geometries into a single feature, which is the right collapse when you want the total covered area rather than per-store catchments. To keep per-store catchments separate, skip the st_union() step.

sq_a <- st_polygon(list(rbind(c(0,0), c(10,0), c(10,10), c(0,10), c(0,0))))
sq_b <- st_polygon(list(rbind(c(5,5), c(15,5), c(15,15), c(5,15), c(5,5))))
sq_a_sfc <- st_sfc(sq_a, crs = 32617)
sq_b_sfc <- st_sfc(sq_b, crs = 32617)

ex_3_3 <- st_area(st_intersection(sq_a_sfc, sq_b_sfc))
ex_3_3
#> 25 [m^2]

Explanation: Each ring must be closed (first point equals last), or st_polygon() raises a coercion error. st_intersection() returns the geometry of the overlap; if the inputs are disjoint it returns an empty geometry rather than NULL, which is friendlier for piped pipelines. The shared region here is a 5x5 square (area 25), exactly what st_area() reports. For polygon-on-polygon clipping at scale, st_intersection() with prepared = TRUE (default in sf 1.0+) uses spatial indexes for big speedups.

stores_nad27 <- st_transform(ex_1_3, 4267)
ex_4_1 <- st_join(
  stores_nad27,
  ex_1_1 |> dplyr::select(NAME),
  join = st_within
)
ex_4_1

Explanation: st_join() is the spatial analog of a left join: every left-hand feature is preserved, and matching right-hand attributes are attached. The join argument controls the predicate (st_within, st_intersects, st_contains, st_touches, st_nearest_feature for a fallback when nothing intersects). Always reproject both sides to the same CRS first: sf 1.0+ raises an error on mismatched CRS rather than silently treating one as the other. st_within is the right predicate for "point inside polygon" because st_intersects would also match boundary touches.

raleigh <- st_sfc(st_point(c(-78.6382, 35.7796)), crs = 4326) |>
  st_transform(3857)
nc_3857 <- st_transform(ex_1_1, 3857)

within_idx <- st_is_within_distance(nc_3857, raleigh, dist = 80000, sparse = FALSE)[, 1]
ex_4_2 <- sort(ex_1_1$NAME[within_idx])
ex_4_2

Explanation: st_is_within_distance() returns a sparse logical matrix by default; passing sparse = FALSE flips it to a regular logical matrix you can subset with directly. The 80km cutoff is measured edge-to-edge, not centroid-to-centroid, so a large county whose boundary clips the radius is counted even if its centroid lies outside. For a strict centroid-distance test, run st_distance(st_centroid(nc_3857), raleigh) and filter on the numeric result. Always reproject to a meter-based CRS first or the dist argument is silently mistreated.

stores_nad27 <- st_transform(ex_1_3, 4267)

nc_counts <- ex_1_1
nc_counts$n_stores <- lengths(st_intersects(nc_counts, stores_nad27))

ex_4_3 <- nc_counts |>
  sf::st_drop_geometry() |>
  dplyr::filter(n_stores > 0) |>
  dplyr::select(NAME, n_stores) |>
  dplyr::arrange(dplyr::desc(n_stores), NAME) |>
  tibble::as_tibble()
ex_4_3

Explanation: st_intersects(A, B) returns a list where element i holds the indices of B that intersect feature i of A; lengths() collapses that list to a per-feature count without materializing a join table. This pattern scales to millions of features because sf uses a spatial index (STRtree) automatically. st_drop_geometry() strips the sfc column to free memory and lets you treat the result as a pure tibble for reporting.

plot(st_geometry(ex_1_1), col = NA, border = "grey30")
ex_5_1 <- "boundaries_plotted"
ex_5_1
#> [1] "boundaries_plotted"

Explanation: plot() on a full sf object draws one panel per attribute column, which is rarely what you want. st_geometry() extracts the geometry alone so you get a single clean panel. col = NA means no fill, leaving only borders, the standard choice for showing structure before encoding any variable. For multi-layer plots (counties + points + roads), call plot(..., add = TRUE) on subsequent layers, or move to ggplot2 which handles layering more cleanly (next exercise).

ex_5_2 <- ggplot(ex_1_1) +
  geom_sf(aes(fill = BIR74), colour = "white", linewidth = 0.1) +
  scale_fill_viridis_c(option = "plasma", name = "BIR74") +
  theme_minimal() +
  labs(title = "Live births by NC county, 1974 to 1978")

print(ex_5_2)

Explanation: geom_sf() automatically uses each layer's CRS and reprojects on the fly, so you do not need to pre-align layers for plotting. Continuous viridis variants (viridis, plasma, magma, inferno, cividis) are perceptually uniform and colourblind-safe, which matters for any report you will share. colour = "white" plus linewidth = 0.1 produces airy borders that read well at print resolution. For a discrete variable, swap in scale_fill_viridis_d(); for diverging data, use scale_fill_distiller(palette = "RdBu").

ex_5_3 <- ex_5_2 +
  geom_sf(data = ex_1_3, shape = 21, fill = "red", colour = "white", size = 3) +
  geom_sf_text(data = ex_1_3, aes(label = store),
               nudge_y = 0.15, size = 3, colour = "black")

print(ex_5_3)

Explanation: Each geom_sf() layer carries its own data argument, and ggplot2 reprojects every layer to the first layer's CRS automatically: here ex_1_3 is in EPSG 4326 and ex_1_1 in EPSG 4267, but the plot renders cleanly anyway. shape = 21 is the filled circle whose border colour is controlled separately from the fill (unlike shape = 19 which has no border). For non-overlapping label placement on dense maps, swap geom_sf_text() for ggrepel::geom_text_repel() after extracting coordinates with st_coordinates().