Week 3: Discussion Section

Practice vector operations

Author

Alessandra Vidal Meza

Published

November 21, 2025

Jack and Laura Dangermond Preserve near Lompoc, CA

In this section, we’ll explore key vector operations inside the sf package using data from California. California is considered a key global biodiversity hotspot. It is home to an exceptional variety of species and ecosystems as well as high endemism. Its mix of elevation, from the highest to the lowest points in the continental U.S., and its location, near both the ocean and mountains, make it geodiverse and so biologically special. As a result, California has one of the most extensive and advanced systems of protected areas in the U.S. It currently conserves more than 25% of its land. At large, area-based conservation is complex and implementation of conservation sites has been a way for governments to displace Indigenous peoples and local communities.

Learning objectives

  • Explore topological relationships with sf functions
  • Explore distance relationships with sf functions
  • Learn about spatial and distance-based joins
  • Practice writing error/warning messages and unit tests to diagnose outputs

1. Get started

  • Create a version-controlled R Project
  • Add (at least) a subfolder to your R project: data
  • Create a Quarto document

Let’s load all necessary packages:

library(here)
library(tidyverse)
library(sf)
library(tmap)

You will be working with the following datasets:

  • Santa Barbara County’s City Boundaries (Santa Barbara County)
  • California Protected Areas Database (CPAD)
  • iNaturalist Research-grade Observations, 2020-2024 (via rinat)

Next, let’s download our data. Unzip and move this to your version-controlled R Project’s data folder.

sb_protected_areas <- st_read(here("data", "cpad_super_units_sb.shp")) |>  
  st_transform("ESRI:102009")

sb_city_boundaries <- st_read(here("data", "sb_city_boundaries_2003.shp")) |> 
  st_transform("ESRI:102009")

sb_county_boundary <- st_read(here("data", "sb_county_boundary_2020.shp")) |> 
  st_transform("ESRI:102009")

aves <- st_read(here("data", "aves_observations_2020_2024.shp")) |> 
  st_transform("ESRI:102009")

2. Find bird observations within PAs in Santa Barbara County

  • Find how many bird observations are within protected areas in Santa Barbara County and show how your output differs with a spatial subset versus a spatial join
aves_PA_subset <- sb_protected_areas[aves, ] # Subset

nrow(aves_PA_subset) # Check number of rows
[1] 35
tm_shape(sb_county_boundary) +
  tm_fill() +
  tm_shape(sb_protected_areas) +
    tm_borders(lwd = 1, col = "#fb8500") +
    tm_fill(col = "#fb8500", alpha = 0.2) +
  tm_shape(aves_PA_subset) +
    tm_dots(col = "#023047")

aves_PA_join <- st_join(aves, sb_protected_areas) # default: `st_intersects` 

nrow(aves_PA_join) # Check number of rows
[1] 500
tm_shape(sb_county_boundary) +
  tm_fill() +
  tm_shape(sb_protected_areas) +
    tm_borders(lwd = 1, col = "#fb8500") +
    tm_fill(col = "#fb8500", alpha = 0.2) +
  tm_shape(aves_PA_join) +
    tm_dots(col = "#023047")

  • Make a 5 km buffer around the protected areas and subset buffered geometries to only those with bird observations
st_crs(sb_protected_areas)$units # Check if units are in meters
[1] "m"
PA_buffer_5km <- st_buffer(sb_protected_areas, dist = 5000) # Create 5000 m buffer around PAs
aves_buffer_subset <- PA_buffer_5km[aves, ] # Subset

nrow(aves_buffer_subset) # Check number of rows
[1] 327
tm_shape(sb_county_boundary) +
  tm_fill() +
  tm_shape(sb_protected_areas) +
    tm_borders(lwd = 1, col = "#fb8500") +
    tm_fill(col = "#fb8500", alpha = 0.2) +
  tm_shape(aves_buffer_subset) +
  tm_dots(col = "#023047")

3. Find PAs within 15km of Goleta

  • Find the protected areas within 15 km of a city in Santa Barbara County
    • Hint: Use dplyr::filter() to select a city from sb_city_boundaries
# Subset SB county to Goleta
goleta <- sb_city_boundaries |> 
  dplyr::filter(NAME == "Goleta")

st_crs(goleta)$units  # Check if units are in meters
[1] "m"
goleta_buffer_15km <- st_buffer(goleta, dist = 15000) # Create 15 km buffer around Goleta
  • Explore different outputs for #3 with st_intersects(), st_intersection(), and st_within()
goleta_PAs_within <- st_within(sb_protected_areas, goleta_buffer_15km)
goleta_PAs_intersect <- st_intersects(sb_protected_areas, goleta_buffer_15km)
goleta_PAs_intersection <- st_intersection(sb_protected_areas, goleta_buffer_15km)
  • Practice a distance-based join with st_is_within_distance()
goleta_PAs_join <- st_join(sb_protected_areas, goleta, st_is_within_distance, dist = 15000)

3. Find distance between Goleta and Dangermond Preserve

  • Find the distance between your city of choice and a protected area of your choice, using the geometries’ edges and the centroid
# Subset PA to Dangermond Preserve
dangermond <- sb_protected_areas |>
  dplyr::filter(UNIT_NAME == "Jack and Laura Dangermond Preserve")
danger_dist <- st_distance(goleta, dangermond) # Compute the distance between geometries edges
# Calculate the geometric center
dangermond_centroid <- st_centroid(dangermond)
goleta_centroid <- st_centroid(goleta)

danger_dist_centroid <- st_distance(goleta_centroid, dangermond_centroid) # Compute the distance between geometries edges
# Check if the distance matrices are equal
danger_dist == danger_dist_centroid
      [,1]
[1,] FALSE