EDS 223: Discussion Week 10 - Workflows and Maps

Prerequsites

library(tidyverse)
library(patchwork)
library(sf)
library(terra)
library(stars)
library(ggspatial)

Load in the CPAD_2023a_SuperUnits.shp and the ghm.tif files. ghm.tif Transform both to EPSG:4326.

cpad_super <- st_read(here::here("discussion", "data", "CPAD_2023a_SuperUnits.shp"), quiet = TRUE) %>%
  sf::st_transform("EPSG:4326") %>%
  janitor::clean_names() %>%
  mutate(ID = row_number())

ghm <- rast(here::here("discussion", "data", "gHM_masked.tif")) %>%
  project("EPSG:4326")

Exercises

Let’s make nice plots of the California Protected areas by access level

p1 = ggplot(data = cpad_super) +
  geom_sf(aes(color = access_typ, fill = access_typ)) +
  theme_bw() +
  labs(
    color = "Access Type",
    fill = "Access Type"
  ) +
  annotation_scale(plot_unit = "km") +
  annotation_north_arrow(
    location = "tr",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = ggspatial::north_arrow_nautical(
      fill =
        c("grey40", "white"),
      line_col = "grey20"
    )
  ) +
  coord_sf() +
  scale_color_viridis_d() +
  scale_fill_viridis_d()

p1 +
  facet_wrap(~access_typ) +
  theme(strip.background = element_rect(fill = "transparent"))

Let’s try plotting the ghm layers nicely too!

ggplot() +
    geom_stars(data = st_as_stars(ghm)) +
  coord_equal() +
  theme_bw() +
  labs(
    x = "",
    y = "",
    fill = "Global Human Modification"
  ) +
  scale_fill_viridis_c() +
  annotation_scale(plot_unit = "km") +
  annotation_north_arrow(
    location = "tr",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = ggspatial::north_arrow_nautical(
      fill =
        c("grey40", "white"),
      line_col = "grey20"
    )
  )

Create a function to take 2 data sets (1 polygon and 1 raster) and create a boxplot of the values based on a specific layer

summary_boxplot = function(polygon, raster, my_layer, my_label) {
  
  # rasterize polygon by layer
  id_rast = rasterize(polygon, raster, field = "suid_nma")
  
  #do mean zonal statistics
  zonal_layer <<- zonal(raster, id_rast, fun = "mean", na.rm = TRUE)
  
  #join with polygon database
  poly_join <<- full_join(polygon, zonal_layer) %>% 
    select(suid_nma, gHM, my_layer)
  
  #create boxplot based on your layer
  p1 = ggplot(poly_join) +
    geom_boxplot(aes(gHM, .data[[my_layer]])) +
    theme_bw() +
    labs(x = "Human Modification Index", 
         y = my_label)
  
  return(p1)
}

Let’s select some layers and use our new function!

names(cpad_super)

##  [1] "suid_nma"   "access_typ" "park_name"  "park_url"   "mng_ag_id" 
##  [6] "mng_agncy"  "mng_ag_lev" "mng_ag_typ" "agncy_web"  "layer"     
## [11] "acres"      "label_name" "yr_est"     "gap1_acres" "gap2_acres"
## [16] "gap3_acres" "gap4_acres" "gap_tot_ac" "shape_leng" "shape_area"
## [21] "geometry"   "ID"

access = summary_boxplot(cpad_super, ghm, "access_typ", "Access Type")

access

layer = summary_boxplot(cpad_super, ghm, "layer", "Management Agency Type")

layer

EDS 223: Discussion Week 10 - Workflows and Maps

Allie Caughman

2023-11-30

Introduction

Prerequsites

Exercises