knitr::opts_chunk$set(warning = FALSE, message = FALSE)

# Load necessary libraries
rm(list = ls())
library(sf)
library(ggplot2)
library(dplyr)
library(ggspatial)
library(ggrepel)

outfig <- "H:/My Drive/Manuscript/BeavePond_westernUS/Figs"

1. load data

boundaries <- st_read("H:/My Drive/Data/GIS_Data/Downloaded/Tiger/tl_2017_us_state/tl_2017_us_state.shp")

## Reading layer `tl_2017_us_state' from data source 
##   `H:\My Drive\Data\GIS_Data\Downloaded\Tiger\tl_2017_us_state\tl_2017_us_state.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 56 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.2311 ymin: -14.60181 xmax: 179.8597 ymax: 71.43979
## Geodetic CRS:  NAD83

names(boundaries)

##  [1] "REGION"   "DIVISION" "STATEFP"  "STATENS"  "GEOID"    "STUSPS"  
##  [7] "NAME"     "LSAD"     "MTFCC"    "FUNCSTAT" "ALAND"    "AWATER"  
## [13] "INTPTLAT" "INTPTLON" "geometry"

# Load US boundaries and filter for the lower 48 states only
us_boundaries <- st_read("H:/My Drive/Data/GIS_Data/Downloaded/Tiger/tl_2017_us_state/tl_2017_us_state.shp") %>%
  filter(!NAME %in% c("Alaska", 
                      "Hawaii", 
                      "Puerto Rico",
                      "Commonwealth of the Northern Mariana Islands",
                      "District of Columbia",
                      "United States Virgin Islands",
                      "American Samoa",
                      "West Virginia",
                      "Guam"))

## Reading layer `tl_2017_us_state' from data source 
##   `H:\My Drive\Data\GIS_Data\Downloaded\Tiger\tl_2017_us_state\tl_2017_us_state.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 56 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.2311 ymin: -14.60181 xmax: 179.8597 ymax: 71.43979
## Geodetic CRS:  NAD83

# load ecoregion and eeager polygon
ecoregions <- st_read("H:/My Drive/Data/GIS_Data/Downloaded/Ecoregions/NA_CEC_Eco_Level3/NA_CEC_Eco_Level3.shp")

## Reading layer `NA_CEC_Eco_Level3' from data source 
##   `H:\My Drive\Data\GIS_Data\Downloaded\Ecoregions\NA_CEC_Eco_Level3\NA_CEC_Eco_Level3.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2548 features and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -4334052 ymin: -3313739 xmax: 3324076 ymax: 4267265
## Projected CRS: Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Area

# Load EEAGER Polygon 
eeager_area <- st_read("H:/My Drive/Projects/Beavers/Data/EEAGER_Fairfax_2023/Processed/EEAGER_GEEapp_Polygons.shp")

## Reading layer `EEAGER_GEEapp_Polygons' from data source 
##   `H:\My Drive\Projects\Beavers\Data\EEAGER_Fairfax_2023\Processed\EEAGER_GEEapp_Polygons.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 12 features and 11 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -121.7745 ymin: 38.863 xmax: -105.6599 ymax: 48.92571
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

# load laramie polygon 
laramie_area <- st_read("H:/My Drive/Projects/Beavers/Data/BeaverDamCatalog/DamTraces_20240617_processed/Laramie_study_area/Laramie_study_area.shp")

## Reading layer `Laramie_study_area' from data source 
##   `H:\My Drive\Projects\Beavers\Data\BeaverDamCatalog\DamTraces_20240617_processed\Laramie_study_area\Laramie_study_area.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 10 fields
## Geometry type: LINESTRING
## Dimension:     XYZ
## Bounding box:  xmin: -105.4625 ymin: 41.14582 xmax: -105.2909 ymax: 41.306
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

# crs 
albers_crs <- st_crs("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

us_boundaries <- st_transform(us_boundaries, crs = albers_crs)
ecoregions <- st_transform(ecoregions, crs = albers_crs)
eeager_area <- st_transform(eeager_area, crs = albers_crs)
laramie_area <- st_transform(laramie_area, crs = albers_crs)

2. study area

# some EEAGER polygons were excluded 
names(eeager_area)

##  [1] "Name"       "FolderPath" "SymbolID"   "AltMode"    "Base"      
##  [6] "Clamped"    "Extruded"   "Snippet"    "PopupInfo"  "Shape_Leng"
## [11] "Shape_Area" "geometry"

study_area <- eeager_area %>% 
  select(-any_of("Shape_Area")) %>% 
  filter(Name != "I") %>%  # channel dams 
  filter(Name != "L") # 2017 - mapped but not ideal

# Increase the size of the study area polygons by applying a buffer (e.g., 0.05 degrees or meters depending on CRS)
# Adjust the buffer size based on your CRS and desired visibility

3. buffer study area

names(study_area)

##  [1] "Name"       "FolderPath" "SymbolID"   "AltMode"    "Base"      
##  [6] "Clamped"    "Extruded"   "Snippet"    "PopupInfo"  "Shape_Leng"
## [11] "geometry"

study_area_simplified <- study_area %>% 
  dplyr::select(Name,Shape_Leng,geometry) %>%
  dplyr::mutate(Name = if_else(Name == "Outline of Study Area", "M", Name))

st_crs(study_area_simplified)  #Albers Equal Area, unit is meters, so don't need to reproject before applying the buffer

## Coordinate Reference System:
##   User input: +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
##   wkt:
## PROJCRS["unknown",
##     BASEGEOGCRS["unknown",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8901]]],
##     CONVERSION["unknown",
##         METHOD["Albers Equal Area",
##             ID["EPSG",9822]],
##         PARAMETER["Latitude of false origin",37.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-96,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",29.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",45.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

st_geometry_type(study_area_simplified) # polygon and linestring

##  [1] POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON
## [10] POLYGON
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

study_area_buffered <- st_buffer(study_area_simplified, dist = 20000) # 30,000 meters buffer

4. ecoregion III

# Filter and prepare ecoregions by NA_L3CODE
ecoregions <- ecoregions %>% 
  filter(!is.na(NA_L3NAME)) %>% 
  mutate(NA_L3NAME = as.factor(NA_L3NAME))

# Filter ecoregions to only those containing the study area
ecoregions_filtered <- ecoregions %>%
  st_filter(study_area, .predicate = st_intersects)

# Create a new column combining NA_L3NAME and NA_L3CODE
ecoregions_filtered <- ecoregions_filtered %>%
  mutate(Name_Code = paste0(NA_L3NAME, " (", NA_L3CODE, ")"))

5. state boundary

# Define the states to highlight
highlight_states <- c("Oregon", "Wyoming", "Montana", "Colorado")
# Create a new layer with only the highlighted states
highlighted_states <- us_boundaries %>%
  filter(NAME %in% highlight_states)

highlighted_states_centroids <- st_centroid(highlighted_states)

6. deal with the order of legend

# Add a new fill category to ecoregions_filtered
ecoregions_filtered <- ecoregions_filtered %>%
  mutate(Fill_Category = as.character(Name_Code))

# Add a new fill category to study_area_buffered
study_area_buffered <- study_area_buffered %>%
  mutate(Fill_Category = "Study Areas")

# Combine into a single dataframe (optional)
combined_data <- bind_rows(
  ecoregions_filtered %>% select(Fill_Category, geometry),
  study_area_buffered %>% select(Fill_Category, geometry)
)

# Define the desired order of fill categories
desired_order <- c("Study Areas", "Wyoming Basin (10.1.4)","Northwestern Glaciated Plains (9.3.1)",
                   "Canadian Rockies (6.2.4)","Eastern Cascades Slopes and Foothills (6.2.8)",
                   "Southern Rockies (6.2.14)")
                   
# Apply the factor with the desired order
combined_data <- combined_data %>%
  mutate(Fill_Category = factor(Fill_Category, levels = desired_order))

ecoregion_colors <- c(
  "Wyoming Basin (10.1.4)" = "#e78ac3",    
  "Northwestern Glaciated Plains (9.3.1)" = "#a6d854",  
  "Southern Rockies (6.2.14)" = "#fc8d62",  
  "Canadian Rockies (6.2.4)" = "#8da0cb",   
  "Eastern Cascades Slopes and Foothills (6.2.8)" = "#66c2a5" 
)

# Combine colors, placing "Study Area" first
combined_colors <- c("Study Areas" = "black", ecoregion_colors)

7. study area map

# Split layers
ecoregions_layer <- combined_data %>% 
  filter(Fill_Category != "Study Areas") %>%
  mutate(Legend_Category = Fill_Category)

study_area_layer <- combined_data %>% 
  filter(Fill_Category == "Study Areas") %>%
  mutate(Legend_Category = "Study Areas")

# Define colors
combined_colors <- c(
  "Study Areas" = "black",
  "Canadian Rockies (6.2.4)" = "#8da0cb",
  "Wyoming Basin (10.1.4)" = "#e78ac3",
  "Eastern Cascades Slopes and Foothills (6.2.8)" = "#66c2a5",
  "Northwestern Glaciated Plains (9.3.1)" = "#a6d854",
  "Southern Rockies (6.2.14)" = "#fc8d62"
)

# Plot
ggplot() +
  # Base layers
  geom_sf(data = us_boundaries, color = "gray40", fill = "gray95") +
  geom_sf(data = highlighted_states, fill = "gray80") + 
  
  # Ecoregions as filled polygons
  geom_sf(data = ecoregions_layer, aes(fill = Legend_Category), color = NA) +
  
  # Study areas as points (centroids)
  geom_sf(
    data = study_area_layer %>% st_centroid(),
    aes(fill = Legend_Category),
    shape = 21,  # filled circle
    size = 2,
    stroke = 0.2,
    color = "black"
  ) +
  
  # State labels
  geom_sf_text(data = highlighted_states_centroids, aes(label = STUSPS),
               size = 5, fontface = "bold", hjust = 0, nudge_x = 0.2) +
  
  # Unified fill legend
  scale_fill_manual(
    name = NULL,
    values = combined_colors,
    breaks = names(combined_colors),  # explicit order
    guide = guide_legend(
      ncol = 2,
      byrow = TRUE,
      override.aes = list(
        shape = c(21, rep(22, 5)),                 # 21: circle for Study Areas, 22: square for others
        size = c(3, rep(5, 5)),                    # smaller for Study Area
        color = c("black", rep(NA, 5)),            # outline only for Study Area
        fill = c("black", combined_colors[-1])     # correct fill colors
    ))) +
  
  # Layout
  annotation_scale(location = "bl") +
  annotation_north_arrow(location = "tr", which_north = "true",
                         style = north_arrow_fancy_orienteering) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.spacing.y = unit(0.2, "cm"),
    legend.text = element_text(size = 10),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    legend.key.height = unit(0.1, "in"),  
    legend.key.width = unit(0.1, "in")
  )

# ggsave(paste0(outfig,"/study_area_map_color2_update1.png"), plot = last_plot(), width = 6.5, height = 6, units = "in", dpi = 600)

Study area

Luwen Wan

2024-10-27

1. load data

2. study area

3. buffer study area

4. ecoregion III

5. state boundary

6. deal with the order of legend

7. study area map