Spatial Data 2 (reference card)

R set-up

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

Mapping polygon data

Load and plot polygon data with geom_sf

harv_soils <- read_sf("data/harv/harv_soils.shp")

ggplot() +
  geom_sf(data = harv_soils, mapping = aes(fill = TYPE_)) +
  scale_fill_viridis_d()

Use facet_wrap to show each category in its own subplot

ggplot() +
  geom_sf(data = harv_boundary) +
  geom_sf(data = harv_soils, fill = "blue") +
  facet_wrap(vars(TYPE_))

Aggregating raster data inside polygons

Use extract with fun = mean to get average raster values per polygon, then convert back to sf

harv_soils_elevs <- extract(harv_dtm, harv_soils, fun = mean, bind = TRUE) |>
  st_as_sf()

ggplot() +
  geom_sf(data = harv_soils_elevs, mapping = aes(fill = harv_dtmfull)) +
  scale_fill_viridis_c()

Controlling plot projection

By default geom_sf displays axes in lat/long. Use coord_sf to display in the CRS of your data instead

ggplot() +
  geom_spatraster(data = harv_dtm) +
  geom_sf(data = harv_soils, alpha = 0) +
  coord_sf(datum = st_crs(harv_dtm))

Rotate x-axis labels to prevent overlap using theme

ggplot() +
  geom_spatraster(data = harv_dtm) +
  geom_sf(data = harv_soils, alpha = 0) +
  coord_sf(datum = st_crs(harv_dtm)) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

Creating vector data from a CSV

Use read_sf with options to load lat/long columns as spatial data. No spaces in options arguments. Use crs = 4326 for lat/long data.

harv_plots <- read_sf("data/harv/harv_plots.csv",
                      options = c("X_POSSIBLE_NAMES=longitude",
                                  "Y_POSSIBLE_NAMES=latitude"),
                      crs = 4326)

Then reproject and use like any other vector data

harv_plots_utm <- st_transform(harv_plots, st_crs(harv_dtm))

Cropping spatial data

Crop raster to a vector polygon

harv_dtm_cropped <- crop(harv_dtm, harv_boundary, mask = TRUE)

mask = TRUE sets values outside the polygon to NA. Use na.value = "transparent" to hide them in plots

ggplot() +
  geom_spatraster(data = harv_dtm_cropped) +
  scale_fill_viridis_c(na.value = "transparent") +
  geom_sf(data = harv_boundary, fill = "transparent")

Crop to a bounding box

Create a bounding box with st_bbox, then crop

bbox <- st_bbox(c(xmin = 731000, ymin = 4713000,
                  xmax = 732000, ymax = 4714000),
                crs = st_crs(harv_dtm))

harv_dtm_small <- crop(harv_dtm, bbox)
harv_soils_small <- st_crop(harv_soils, bbox)

Saving spatial data

Save a raster as GeoTIFF

writeRaster(harv_dtm_small, "harv_dtm_small.tif")

Save vector data as shapefile

write_sf(harv_soils_small, "harv_soils_small.shp")

Save a map the same way as any other ggplot

ggsave("harv_map.png")