Spatial Data (reference card)

R set-up

library(terra)
library(tidyterra)
library(sf)
library(ggplot2)

Importing spatial data

dsm_harv <- rast("data/neon-airborne/harv_dsmCrop.tif")

Metadata is important to describe the context of spatial data. - bands - projection - units - min, max, mean

Plotting spatial data with ggplot

Special geom for plotting terra rasters

ggplot() +
  geom_spatraster(data = dtm_harv)

Change color ramp using scale function. fill colors the inside of each raster pixel. viridis is a good color ramp; _c indicates continuous scale.

ggplot() +
  geom_spatraster(data = dtm_harv) +
  scale_fill_viridis_c()

When converted to a data frame we can use it as a normal table

dsm_harv_df = as.data.frame(dsm_harv, xy = TRUE)

ggplot() +
  geom_histogram(data = dsm_harv_df,
                 aes(x = HARV_dsmCrop))

DSM is a Digital Surface Model: elevation of top physical point
DTM is Digital Terrain Model: elevation of the ground
CHM is Canopy Height Model: difference between DSM and DTM

dsm_harv <- rast("data/neon-airborne/harv_dsmCrop.tif")
dtm_harv <- rast("data/neon-airborne/harv_dtmcrop.tif")
chm_harv <- dsm_harv - dtm_harv

ggplot() +
  geom_spatraster(data = chm_harv)

Vector data

Vector data includes points, lines, and polygons
shapefiles are one common format — a set of multiple files with the same name but different extensions

Read vector data using sf

plots_harv <- read_sf("data/harv/harv_plots.shp")

The result is a data frame with one row per feature, a geometry column storing spatial info, and regular columns (fields) like plot_id and plot_type.

Plot using geom_sf

ggplot() +
  geom_sf(data = plots_harv)

Color by a field

ggplot() +
  geom_sf(data = plots_harv, mapping = aes(color = plot_type))

Combining spatial layers

boundary_harv <- read_sf("data/harv/harv_boundary.shp")

ggplot() +
  geom_sf(data = boundary_harv, fill = "transparent") +
  geom_sf(data = plots_harv)

Layer order matters — later layers plot on top. Use alpha for transparency or fill = "transparent" for outlines only.

Projections

Raster and vector data may use different coordinate reference systems (CRS). Check with:

st_crs(dtm_harv)
st_crs(plots_harv)

Reproject vector data to match the raster using st_transform

plots_harv_utm <- st_transform(plots_harv, st_crs(dtm_harv))

Now they can be plotted together

ggplot() +
  geom_spatraster(data = dtm_harv) +
  geom_sf(data = plots_harv_utm)

UTM is the most common CRS in ecological research — it preserves distance and works well at local/regional scales
Lat/long (WGS84) is common for data collection but doesn’t preserve distance or area

Extract raster data

Extract raster values at point locations

plot_elevations <- extract(dtm_harv, plots_harv_utm)

Use bind = TRUE and st_as_sf() to return a spatial data frame with extracted values attached

plot_harv_elevations <- extract(dtm_harv, plots_harv_utm, bind = TRUE) |>
  st_as_sf() |>
  rename(elevations = harv_dtmcrop)

Extract with buffer

Buffer points to extract values within a surrounding area

plots_harv_buffered <- st_buffer(plots_harv_utm, dist = 25)

Units for dist match the CRS units (meters for UTM). Then extract the mean value within each buffer

plot_elevations_buffered <- extract(dtm_harv, plots_harv_buffered, fun = mean, bind = TRUE) |>
  st_as_sf()

Map of point data

library(maps)
us_map = map_data("usa")

ggplot() +
  geom_polygon(data = us_map,
               aes(x = long, y = lat, group = group),
               fill = "grey") +
  coord_quickmap()

Add occurrence data on top

library(spocc)
do_gbif = occ(query = "Dipodomys ordii",
              from = "gbif",
              limit = 1000,
              has_coords = TRUE)
do_data = data.frame(dipo_df$gbif$data)

ggplot() +
  geom_polygon(data = us_map,
               aes(x = long, y = lat, group = group),
               fill = "grey") +
  geom_point(data = do_data,
             aes(x = Dipodomys_ordii.longitude,
                 y = Dipodomys_ordii.latitude)) +
  coord_quickmap()

Namespacing

library(dplyr)
library(terra)

terra can conflict with dplyr’s select. Use dplyr::select() to be explicit.

Making your own vector data

do_data_spat <- vect(
  do_data,
  geom = c('Dipodomys_ordii.longitude', 'Dipodomys_ordii.latitude'),
  crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

Once created, you can reproject with project() and extract values from rasters.