In annack84/STMdevelopment: Develop State-and-Transition Models for the Upper Colorado River Basin
NOTE: This correspondence is intended for communication of project progress among funders and collaborators only. This information is preliminary and is subject to revision. It is being provided to meet the need for timely best science. The information is provided on the condition that neither the U.S. Geological Survey nor the U.S. Government shall be held liable for any damages resulting from the authorized or unauthorized use of the information.
General information
knitr::opts_chunk$set(echo = TRUE)
target_ESG <- "Semiarid_Warm_Shallow_DeepRocky"
user <- "Duniway_et_al_2024" # use this user to get the exact data used for the Feb 2023 JFSP workshops
#user <- "Anna" # use one of these two options to get the most current data
#user <- "VPN"
output_figure_folder <- file.path("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures", target_ESG)
library(dplyr)
library(kableExtra)
library(sf)
library(STMdevelopment)
# library(devtools);load_all()
# TODO insert dot graphic and text about the status of the ESD (e.g. "Provisional"). Might help
# to create a function to automatically insert the description based on the
# status input
EDIT_map(target_ESG = target_ESG, user = user, maxZoomStyle = "public")
Areas shown in blue indicate the maximum mapped extent of this ecological site. Other ecological sites likely occur within the highlighted areas. It is also possible for this ecological site to occur outside of highlighted areas if detailed soil survey has not been completed or recently updated.
Associated ecological site concepts
This ecological site group includes a variety of upland ecological sites that have shallow to deep soils with high rock content. The group was defined by placing ecological sites into groups by separating soils, climate, and geomorphic features to best differentiate reference vegetation production values and documented ecological states as documented by Nauman et al. (2022).
esds <- esd_data_pull(user=user, target_ESG = target_ESG)
comps <- comp_data_pull(user=user, target_ESG = target_ESG)
num_esds <- nrow(esds)
This ecological site group includes r num_esds ecological sites. These sites were correlated to the group based on average soil, climate and geomorphic properties of the soil components linked to each ecological site.
linked_esd_names <- esds[,c("ecoclassid","ecoclassname")] %>%
mutate(ecoclassid = paste0("[", ecoclassid, "](https://edit.jornada.nmsu.edu/catalogs/esd/",
substr(ecoclassid, 2, 5), "/", ecoclassid, ")"))
## Associated ecological sites
format_tables_EDIT_style(data = linked_esd_names,
col.names=c("ESD Code","ESD Name"),
caption = "Ecological Sites associated with this ESG",
row.names = F)
#knitr::kable(x = esds[,c("ecoclassid","ecoclassname")], col.names=c("ESD Code","ESD Name"), caption = NULL,row.names = F)
# TODO create table of dominant plant functional groups and species
esg_prod_list <- esg_production_pull(user=user, target_ESG = target_ESG)
dom_plants_df <- dplyr::tribble(~Item, ~Description,
"Trees", as.character(esg_prod_list[["Tree.df"]]$COMMON_NAME)[1:5],
"Shrubs", as.character(esg_prod_list[["Shrub.df"]]$COMMON_NAME)[1:5],
"Grasses",as.character(esg_prod_list[["Grass.df"]]$COMMON_NAME)[1:5],
"Forbs", as.character(esg_prod_list[["Forb.df"]]$COMMON_NAME)[1:5])
format_tables_EDIT_style(data = dom_plants_df,
caption = "Dominant plant species")
#knitr::kable(x = dom_plants_df, caption = "Table 1. Dominant plant species")
TODO insert "Download full description" link to a PDF
Physiographic features
# TODO create table of landforms, elevations, slopes, water table depth, flooding duration, flooding frequency, ponding frequency
# ideally should be able to tab between US and metric system measurements
landforms <- unlist(strsplit(esds$landfrms[1], split=","))
for(e in 2:nrow(esds)){
esdlfs <- unlist(strsplit(esds[e,c("landfrms")],split=","))
landforms <- append(landforms,esdlfs)
}
landforms <- landforms[!landforms %in% "NA"]
lfsum <- factor(landforms)
lfsum <- reorder(lfsum,lfsum,FUN=function(x) -length(x))
landforms_unique <- unique(landforms)
These sites occur on upland positions on a variety of landforms (n=r length(landforms_unique)). Moderate to steep slopes and rugged topography are common. Less commonly, some sites may occur on high-energy alluvial deposits.
# Slope
slope_ave <- mean(comps@site$slope_r.x, na.rm = T)
slope_75th <- unname(quantile(comps@site$slope_r.x, probs = 0.75, na.rm = T))
slope_max <- max(comps@site$slope_r.x, na.rm = T)
slope_sd <- sd(comps@site$slope_r.x, na.rm = T)
slope_field <- paste("Slopes average",format(round(slope_ave, 1), nsmall = 1),
"% with a standard deviation of",format(round(slope_sd, 1), nsmall = 1),
"and are generally less than", format(round(slope_75th, 1), nsmall = 1),
"but can be up to", format(round(slope_max, 1), nsmall = 1),"%",sep=" ")
phys.df <- dplyr::tribble(~Item, ~Description,
"Landforms (Top 10)", as.character(unique(lfsum))[1:10],
"Slope", slope_field,
"Flooding Frequency","None",
"Aspect", "Aspect is not a significant factor")
format_tables_EDIT_style(data = phys.df, caption = "Representative physiographic features")
#knitr::kable(x = phys.df, caption = "Table 2. Representative physiographic features")
Climatic features
Climate is generally semiarid and warm with aridity index (precipitation / potential evapotranspiration) values averaging r format(round(mean(esds$aimean), 2), nsmall = 2) and ranging from r format(round(min(esds$aimean), 2), nsmall = 2) to r format(round(max(esds$aimean), 2), nsmall = 2). The average maximum temperatures (Celsius) of the hottest month range from r format(round(min(esds$maxtempmean), 2), nsmall = 2) to r format(round(max(esds$maxtempmean), 2), nsmall = 2), and minimum temperatures of the coldest month range from r format(round(min(esds$mintempmean), 2), nsmall = 2) to r format(round(max(esds$mintempmean), 2), nsmall = 2). Warm season (June to September) precipitation makes up r format(round(mean(esds$pptrt*100), 0), nsmall = 0)% of total precipitation on average, but can range from r format(round(min(esds$pptrt*100), 0), nsmall = 0)% to r format(round(max(esds$pptrt*100), 0), nsmall = 0)%.
clim.df <- dplyr::tribble(~Item, ~Description,
"Frost Free Period (days)", paste(format(round(min(comps@site$ffd_r), 0), nsmall = 0),format(round(mean(comps@site$ffd_r), 0), nsmall = 0),format(round(max(comps@site$ffd_r), 0), nsmall = 0),sep=", "),
"Mean Annual Precipitation (in)", paste(format(round(min(comps@site$reannualprecip_r/25.4,na.rm=T), 1), nsmall = 1),format(round(mean(comps@site$reannualprecip_r/25.4,na.rm=T), 1), nsmall = 1),format(round(max(comps@site$reannualprecip_r/25.4,na.rm=T), 1), nsmall = 1),sep=", "))
format_tables_EDIT_style(data = clim.df, caption = "Representative climatic features (min, mean, max)")
#knitr::kable(x = clim.df, caption = "Table 3. Representative climatic features (min, mean, max)")
# TODO create line chart and bar charts of monthly high/low precip that you can tab between
# TODO create line chart and bar charts of monthly high/low temps that you can tab between
Water features
These sites neither benefit significantly from run-in moisture nor experience excessive loss of moisture from runoff.
Soil features
Soils in this group have shallow to deep depth to bedrock and typically have > 30% rock fragments by volume. Surface and subsurface horizons have a variety of textures ranging from sand to clay loam. Soils are non-saline and non sodic and often have high calcium carbonate contents. Soil pH ranges from r min(c(round(min(comps@horizons[comps@horizons$hzdept_r<30,c("ph1to1h2o_r")],na.rm=T), 0), round(min(comps@horizons[comps@horizons$hzdept_r>30,c("ph1to1h2o_r")],na.rm=T), 0))) to r max(c(round(max(comps@horizons[comps@horizons$hzdept_r<30,c("ph1to1h2o_r")],na.rm=T), 1), round(max(comps@horizons[comps@horizons$hzdept_r>30,c("ph1to1h2o_r")],na.rm=T), 1))). Erosion hazard is generally low but bare ground states are possible.
## Prep soil variables for table
surftext <- factor(esds$txtnm_surf)
surftext <- reorder(surftext,surftext,FUN=function(x) -length(x))
surftext_uni <- levels(surftext)
subtext <- factor(esds$txtnm_sub)
subtext <- reorder(subtext,subtext,FUN=function(x) -length(x))
subtext_uni <- levels(subtext)
drainage <- factor(comps@site$drainagecl.x)
drainage <- reorder(drainage,drainage,FUN=function(x) -length(x))
drainage_uni <- levels(drainage)[1:3]
## Make Table of soil descriptors:
# TODO Could update pH to be based on esd averages instead of comps. Also could update PM to grab ssurgo values
soil.df <- dplyr::tribble(~Item, ~Description,
"Parent Material", "alluvium and eolian sediments",
"Surface Texture (0-30cm)", surftext_uni,
"Subsurface Texture (>30cm)",subtext_uni,
"Drainage", drainage_uni,
"Soil Depth", paste(format(round(min(esds$depthmean,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$depthmean,na.rm=T), 0), nsmall = 0), " cm",sep=""),
"Surface Rock Content %vol (0-30cm)", paste(format(round(min(esds$rock_surf,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$rock_surf,na.rm=T), 0), nsmall = 0), "%",sep=""),
"Subsurface Rock Content %vol (>30cm)", paste(format(round(min(esds$rock_sub,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$rock_sub,na.rm=T), 0), nsmall = 0), "%",sep=""),
"Surface Electrical Conductivity (0-30cm)", paste(format(round(min(esds$ec_surf,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$ec_surf,na.rm=T), 1), nsmall = 1), " dS/m",sep=""),
"Subsurface Electrical Conductivity (>30cm)", paste(format(round(min(esds$ec_sub,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$ec_sub,na.rm=T), 1), nsmall = 1), " dS/m",sep=""),
"Surface Sodium Adsorption Ratio (0-30cm)", paste(format(round(min(esds$sar_surf,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$sar_surf,na.rm=T), 1), nsmall = 1),sep=""),
"Subsurface Sodium Adsorption Ratio (>30cm)", paste(format(round(min(esds$sar_sub,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$sar_sub,na.rm=T), 1), nsmall = 1),sep=""),
"Surface 1:1 pH (0-30cm)", paste(format(round(min(comps@horizons[comps@horizons$hzdept_r<30,c("ph1to1h2o_r")],na.rm=T), 0), nsmall = 0),"-",format(round(max(comps@horizons[comps@horizons$hzdept_r<30,c("ph1to1h2o_r")],na.rm=T), 1), nsmall = 1),sep=""),
"Subsurface 1:1 pH (>30cm)", paste(format(round(min(comps@horizons[comps@horizons$hzdept_r>30,c("ph1to1h2o_r")],na.rm=T), 0), nsmall = 0),"-",format(round(max(comps@horizons[comps@horizons$hzdept_r>30,c("ph1to1h2o_r")],na.rm=T), 1), nsmall = 1),sep=""),
"Surface Calcium Carbonate (0-30cm)", paste(format(round(min(esds$caco3_surf,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$caco3_surf,na.rm=T), 0), nsmall = 0), "%",sep=""),
"Subsurface Calcium Carbonate (>30cm)", paste(format(round(min(esds$caco3_sub,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$caco3_sub,na.rm=T), 0), nsmall = 0), "%",sep=""))
format_tables_EDIT_style(data = soil.df, caption = "Representative soil features")
#knitr::kable(x = soil.df, caption = "Table 4. Representative soil features")
Soil Component data
Soil types correlated to this group include r length(unique(comps@site$cokey)) different components mapped in SSURGO. Soil taxonomic units include:
cat(sort(unique(comps@site$compname)), sep = ", ")
Component level soil property depth plots
agg <- aqp::slab(comps, fm= ~ claytotal_r + silttotal_r + sandtotal_r + ec_r + caco3_r)
lattice::xyplot(top ~ p.q50 | variable, data=agg, ylab='Depth (cm)',
xlab='median bounded by 25th and 75th percentiles',
lower=agg$p.q25, upper=agg$p.q75, ylim=c(201,-2),
panel=aqp::panel.depth_function,
alpha=0.25, sync.colors=TRUE,
par.settings=list(superpose.line=list(col='RoyalBlue', lwd=2)),
prepanel=aqp::prepanel.depth_function,
cf=agg$contributing_fraction, cf.col='black', cf.interval=20,
layout=c(5,1), strip=lattice::strip.custom(bg=grey(0.8)),
scales=list(x=list(tick.number=4, alternating=3, relation='free'))
)
Ecological dynamics
These site are susceptible to drought, especially in shallower soils. Erosion issues are possible but rare. Compared to reference states, alternative states may have woody encroachment, high bare ground, and/or loss of biological soil crust and perennial grasses.
State and transition model development
# TODO update with the newly tabulated ESD STMs from Travis, then set to eval = TRUE
apriori_stms <- read.csv(data_file_paths(user)$apriori_stms, stringsAsFactors = F,
na.strings = c("NA", "", " "))
apriori_stms_target <- select(apriori_stms, any_of(c("State", "Plant_community", target_ESG)))
n_ESDs <- apriori_stms_target[apriori_stms_target$Plant_community=="n_ESDs", target_ESG]
colnames(apriori_stms_target) <- c("State", "Plant community", "Proportion of ESDs containing this plant community")
format_tables_EDIT_style(data = apriori_stms_target[-nrow(apriori_stms_target),],
caption = paste("State and plant community summary from Ecological Sites within this ESG.",
n_ESDs, "complete STMs were available to include.")) %>%
collapse_rows(columns = 1, valign = "top")
Reference state production indicators
## Plot of reference production indicators used in paper
resp_vars <- colnames(esds)[1:21]
resp_vars <- resp_vars[!resp_vars %in% c("ecoclassid")]
par(mar = c(4,8,4,2)) ## This works when knitted
boxplot(esds[,resp_vars], horizontal = TRUE,las = 1, main = "Reference Production (lbs/Acre)")
# by functional group cover
indicators <- c("AH_C3NativePerenGrassCover",
"AH_C3IntroducedPerenGrassCover",
"AH_C4NativePerenGrassCover",
"AH_C4IntroducedPerenGrassCover",
"AH_NativePerenForbCover", # Native perennial forbs
"AH_IntroducedPerenForbCover", # Non-native perennial forbs
"AH_NativeAnnGrassCover", # Native annual grasses
"AH_IntroducedAnnGrassCover", # Non-native annual grasses
"AH_NativeAnnForbCover", # Native annual forbs
"AH_IntroducedAnnForbCover", # Non-native annual forbs
"AH_ArtemisiaTridentataCover", # All big sagebrush species lumped together
"BareSoilCover", # Bare soil
"FH_TotalLitterCover", # Litter
#"CA_percent_100plus", # Canopy gaps > 100 cm
#"CA_percent_200plus",
"FH_LichenCover", # Lichen + moss combined cover
"FH_MossCover"#,
#"SoilStab_all" # to represent cyano abundance... I think this is only marginally a community variable. If used in hierarchical clustering and Bray-Curtis NMDS, will need to standardize it to make units more similar to % cover units
)
ann_grass_by_spp <- F
ann_forb_by_spp <- F
per_grass_by_spp <- F
per_forb_by_spp <- F
succulent_by_spp <- F
shrub_by_spp <- T # All shrubs and sub-shrubs by species
subshrub_by_spp <- T
tree_by_spp <- T # All trees by species
opuntia_combined <- T
data_sources <- c(#"BadgerWash",
#"CRCMonitoring", # drop if not needed for spatial representation
"IM_NCPN",
"LMF",
"NRI",
#"Parashant", # drop if not needed for spatial representation
"AIM"#,
#"VanScoyocThesis" # drop if not needed for spatial representation
)
impute_gap_type <- NULL
plot_data <- plot_data_pull(user=user,
target_ESG = target_ESG,
data_sources = data_sources,
indicators = indicators,
ann_grass_by_spp = ann_grass_by_spp,
ann_forb_by_spp = ann_forb_by_spp,
per_grass_by_spp = per_grass_by_spp,
per_forb_by_spp = per_forb_by_spp,
succulent_by_spp = succulent_by_spp,
shrub_by_spp = shrub_by_spp,
subshrub_by_spp = subshrub_by_spp,
tree_by_spp = tree_by_spp,
opuntia_combined = opuntia_combined,
impute_gap_type = impute_gap_type
) %>%
filter(!is.na(Month))
# TODO update I&M species codes to match current USDA Plants codes! probably a step for the plot_data_pull function?
# remove incomplete rows - ordination won't work with NAs
plot_data_clean <- na.omit(plot_data) %>%
# NRI has a separate plot code for each time a plot was sampled, but all other data sets have the same plot
# code for all sampling times. Make a similar code for NRI.
mutate(PlotCode_NoYear = ifelse(test = SourceKey == "NRI_UCRB",
yes = paste0(SourceKey, "_",
SiteName, "_",
substr(PlotID, 5, nchar(PlotID))),
no = PlotCode))
# Some plots were sampled multiple times - I'll include only the first year that a plot has complete data
plot_data_first <- plot_data_clean %>% #plot_data_clean %>%
dplyr::group_by(PlotCode_NoYear) %>%
dplyr::arrange(.data=., Year, .by_group = TRUE) %>%
dplyr::filter(dplyr::row_number()==1) %>%
dplyr::ungroup()
# remove plots with SGU probability (i.e. certainty of prediction) in the 10th percentile
SGU_prob_raster <- raster::raster("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps/UCRB_SGUs_ProbMax/UCRB_SGUs_ProbMax.tif")
file_paths <- data_file_paths(user)
plot_locations <- sf::st_read(dsn = file.path(file_paths$plotnet_processed, "PlotLocations"),
layer = "all_plot-years_2022-11-08",
quiet=TRUE) %>%
distinct()
if("NRI" %in% data_sources){
nri_locations <- sf::st_read(dsn = file_paths$nri,
layer = "NRI_UCRB_plot-years_2022-11-03",
quiet = TRUE)
plot_locations <- dplyr::filter(plot_locations, !grepl(pattern = "^NRI_",
x=PlotCode)) %>%
dplyr::bind_rows(., nri_locations) %>%
distinct()
}
plot_SGU_prob <- sf::st_as_sf(raster::extract(x = SGU_prob_raster,
y = plot_locations,
sp = T))
plot_data_first <- left_join(plot_data_first, select(plot_SGU_prob, PlotCode, UCRB_SGUs_ProbMax)) %>%
distinct()
certainty_cutoff_level <- 0.1
certainty_cutoff <- quantile(plot_data_first$UCRB_SGUs_ProbMax, probs = certainty_cutoff_level, na.rm = T)
plot_data_first <- filter(plot_data_first, UCRB_SGUs_ProbMax > certainty_cutoff) %>%
filter(PlotCode != "AIM_Utah Vernal FO 2019_009") # this one is duplicated in the plot locations with very different coordinates
# make the clustering and ordination data frame - can't include ANY columns except the actual variables!
ord.df <- dplyr::select(plot_data_first, -SourceKey, -PlotID, -SiteName, -PlotName,
-Year, -Longitude_NAD83, -Latitude_NAD83,
-Month, -Day,
-UCRB_SGUs_ProbMax, -geometry, -PlotCode_NoYear
) %>%
tibble::column_to_rownames("PlotCode") # keep an ID code for the plot as the row name to prevent data scrambling problems
# If using soil stability, standardize variables to make units more similar
# Here, dividing by the max possible value for each measurement type (cover values
# divided by 100, soil stability divided by 6)
if("SoilStab_all" %in% indicators){
ord.df.cover <- select(ord.df, -SoilStab_all)/100
ord.df.soilstab <- select(ord.df, SoilStab_all)/6
ord.df <- cbind(ord.df.cover, ord.df.soilstab)
}
# remove rare species
sp_pa <- mutate(rowwise(ord.df), across(everything(), function(x){if(x>0){1}else{0}})) # make a presence/absence data frame
sp_keep <- names(colSums(sp_pa)[which(colSums(sp_pa)>2)])
ord.df <- select(ord.df, all_of(sp_keep))
# remove unknown species if present
ord.df <- select(ord.df, -any_of(c("UNKS", "UNKPG", "UNKAF", "UNKPF", "PG1", "PG01", "AF1")))
# # save clustering data for other analyses - exclude location data to keep kosher
# # with the data sharing restrictions
# ord.df.share <- select(plot_data_first, all_of(c("SourceKey", "PlotCode", "Year", "Month", "Day", "UCRB_SGUs_ProbMax", colnames(ord.df))))
# write.csv(ord.df.share,
# file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Data/Analysis_ready_data/",
# target_ESG, "_FANNYInputData_", Sys.Date(), ".csv"),
# row.names = F)
Field data from were collected by field crews with the Natural Resource Conservation Service's (NRCS) National Resource Inventory (NRI) program, the Bureau of Land Management's (BLM) Assessment, Inventory, and Monitoring (AIM) and Landscape Monitoring Framework (LMF) programs, and the Nationa Park Service's (NPS) Inventory and Monitoring Northern Colorado Plateau Network (NCPN) program between r min(plot_data_clean$Year) and r max(plot_data_clean$Year). Plots within the target ESG with SGU classification probability less than or equal to the r certainty_cutoff_level*100th percentile (r certainty_cutoff%) were excluded, resulting in a total of r nrow(ord.df) plots for this ESG. Plant and soil cover were measured using line-point intercept and canopy gap methods (Herrick et al., 2017). Plot-level indicators (Table \@ref(tab:doc-plot-indicators)) were calculated using the terradactyl package in R (McCord et al., 2022). When plots were sampled in multiple years, only the first sampling year was used in STM development. r if(!is.null(impute_gap_type)){paste("For plots that were missing annual & perennial gap data but had perennial-only gap data (primarily I&M-NCPN), annual & perennial gap values were imputed from perennial-only gaps and annual herbaceous cover using a linear model built from LMF and NRI plots that had both types of gap (r-squared of 0.83 and 0.81 for 100-cm and 200-cm gap predictions, respectively).")}
rm(plot_data_clean, plot_data)
indicator_descriptions <- make_indicator_descriptions(indicators = indicators,
shrub_by_spp = shrub_by_spp,
subshrub_by_spp = subshrub_by_spp,
tree_by_spp = tree_by_spp,
opuntia_combined = opuntia_combined)
format_tables_EDIT_style(data = indicator_descriptions[,-1],
caption = "Indicators used in ordination and clustering analysis to develop plant communities") %>%
collapse_rows(columns = 1, valign = "top")
# are the data clusterable? Use the Hopkin's statistic (Lawson, Richard, and Jurs, 1990; see https://www.datanovia.com/en/lessons/assessing-clustering-tendency/ for explanation)
# H values > 0.5 indicate clustering tendency (i.e. if H is below 0.5, the data set is uniformly distributed and clusters are not meaningful)
res_LoamyUplands <- factoextra::get_clust_tendency(data = ord.df, n=nrow(ord.df)-1, graph = FALSE)
res_LoamyUplands$hopkins_stat
# cluster and indicator species analysis
library(cluster) # for cluster analysis
library(dendextend)
library(labdsv)
library(foreach)
library(vegan)
library(ggplot2)
# initial clustering
dist.df <- vegan::vegdist(ord.df, method = "bray", na.rm = T) # creates Bray-Curtis distance matrix
# This chunk determines the best combination of number of clusters and fuzziness for the data. It is SLOOOOWWW so better to run it before knitting to get the right clustering parameters and then not run it when knitting the descriptive document.
library(parallel)
library(doParallel)
# get index values for all combos of membership exponent and number of clusters
indices_arg <- c("ch", #"duda", won't work with the ranking scheme
"cindex",
#"gamma", # takes a long time
#"beale", # unclear how this one works
"ccc", "ptbiserial",
#"gplus", # takes a long time
"db", "silhouette", "dunn" #,
#"gap" # takes a long time
)
indices_name <- c("CH", #"Duda",
"Cindex",
#"gamma", # takes a long time
#"Beale",
"CCC", "Ptbiserial",
#"gplus", # takes a long time
"DB", "Silhouette", "Dunn" #,
#"gap" # takes a long time
)
memb_exps <- seq(1.1,
1.3, # start with 2.0, can reduce value to rerun if some larger values don't give results
by = 0.1)
cl <- makeCluster(spec = 3, # number of cores to use
type = "PSOCK",
methods = FALSE)
registerDoParallel(cl)
cluster_eval_fanny <- foreach(memb_exp=memb_exps,
.packages = c("STMdevelopment", "dplyr", "cluster"),
.combine = "rbind") %dopar% {
nbmetrics_fanny_temp <- NbClust_fanny(data = ord.df,
diss = dist.df, #chord.df,
distance = NULL,
min.nc = 2,
max.nc = 12,
method = "fanny",
index = "all", # could also try alllong to get a few more metrics - this is SLOW
alphaBeale = 0.1,
memb.exp= memb_exp # 1.2
)
cluster_eval_fanny_temp <- nbmetrics_fanny_temp$All.index %>%
as.data.frame() %>%
mutate(Number_clusters = rownames(.),
Memb_exp = memb_exp)
}
registerDoSEQ()
stopCluster(cl)
cluster_eval_fanny_indices <- select(cluster_eval_fanny, any_of(c("Memb_exp", "Number_clusters", indices_name)))
# calculate percent worst than best value for each memb. exp. X num. clusters combo for each index
pct_off_best_higher <- function(x, all_vals){
best <- max(all_vals)
pct_off <- ((best-x)/best)*100
return(pct_off)
}
pct_off_best_lower <- function(x, all_vals){
best <- min(all_vals)
pct_off <- ((x-best)/best)*100
return(pct_off)
}
cluster_eval_fanny_indices_rank <- cluster_eval_fanny_indices %>%
rowwise() %>%
mutate(CH_rank = pct_off_best_higher(CH, cluster_eval_fanny_indices$CH),
Cindex_rank = pct_off_best_lower(Cindex, cluster_eval_fanny_indices$Cindex),
CCC_rank = pct_off_best_higher(CCC, cluster_eval_fanny_indices$CCC),
Ptbiserial_rank = pct_off_best_higher(Ptbiserial, cluster_eval_fanny_indices$Ptbiserial),
DB_rank = pct_off_best_lower(DB, cluster_eval_fanny_indices$DB),
Silhouette_rank = pct_off_best_higher(Silhouette, cluster_eval_fanny_indices$Silhouette),
Dunn_rank = pct_off_best_higher(Dunn, cluster_eval_fanny_indices$Dunn)
) %>%
mutate(overall_rank = sum(CH_rank, # lowest overall rank wins!
Cindex_rank,
CCC_rank,
Ptbiserial_rank,
DB_rank,
Silhouette_rank,
Dunn_rank,
na.rm = T))
best_n_clust <- as.numeric(cluster_eval_fanny_indices_rank[[which.min(cluster_eval_fanny_indices_rank$overall_rank), "Number_clusters"]])
best_memb_exp <- as.numeric(cluster_eval_fanny_indices_rank[[which.min(cluster_eval_fanny_indices_rank$overall_rank), "Memb_exp"]])
# VISUALLY CHECK for ties and ranking weirdness before proceeding!!
k_groups_fuzzy <- 8 # best is 8 (rank 74), then 9 (rank 87) then 10 (rank 105) then 7 (rank 108) - all with 1.1 memb. exp.
memb_exp <- 1.1 #1.1
set.seed(1)
clust_fuzz <- fanny(x=dist.df,
k = k_groups_fuzzy,
diss = T,
#metric = "euclidean",
memb.exp = # this affects how crisp/fuzzy the clusters are - closer to 1 is crisper
memb_exp
)
grp_fuzz <- clust_fuzz$clustering
plot_data_first$PlantCommunity_fuzzy_original <- grp_fuzz
# double check groups with original grouping - seems like changing the file path
# to the data reset the grouping algorithm somehow so that the same groups are
# made but they are given different numbers by fanny(). If using the Duniway et
# al 2024 data, groups should line up like this:
# Juniper-pinyon woodland, n = 229. Originally cluster #1. Workshop feedback was to lump with 5 Pinyon-juniper woodland
# Open juniper-pinyon shrubland, n = 238. Originally cluster #2.
# C4 grassland, n = 110. Originally cluster #3.
# Annualized juniper savanna, n=113. Originally cluster #4
# Pinyon-juniper woodland, n = 270. Originally cluster #5. Workshop feedback was to lump with 1 Juniper-pinyon woodland
# Gambel oak woodland, n=112. Originally cluster #6.
# Sagebrush shrubland/grassland, n=198. Originally cluster #7.
# Bare juniper savanna, n=186. Originally cluster #8.
# Based on matching up clusters by number of plots, I'll reassign numbers to match
# the original groupings:
updated_cluster_numbers <- data.frame("PlantCommunity_fuzzy_original" = c(1, 2, 3, 4, 5, 6, 7, 8),
"PlantCommunity_fuzzy" = c(1, 4, 5, 2, 6, 7, 3, 8))
plot_data_first <- left_join(plot_data_first, updated_cluster_numbers, by="PlantCommunity_fuzzy_original")
Plots were clustered into ecological communities using the FANNY fuzzy clustering algorithm (Kaufman and Rousseeuw, 1990) in the "cluster" R package (Maechler et al., 2021) with Bray-Curtis distance ("vegan" package; Oksanen et al., 2020). A suite of clustering indices were used to determine the most appropriate number of clusters and degree of fuzziness (i.e. membership exponent), including Calinski and Harabasz (1974), C-index (Hubert and Levin, 1976), cubic clustering criteria (Sarle, 1983), point-biserial correlation (Milligan 1980, 1981), Davies and Bouldin (1979), silhouette width (Rousseeuw 1987), and Dunn's index (1974). Values between 1 and 12 possible clusters and 1.1 to 1.3 membership exponents were tested. The final clustering solution used r k_groups_fuzzy clusters with a membership exponent of r memb_exp.
# another slow chunk - run before knitting and then don't evaluate in the knit step.
NMDS_scree(ord.df)
ord_dims <- 3 # change manually based on the plot results
#run NMS ordination
set.seed(1)
ord <- metaMDS(ord.df,
k=ord_dims, # number of dimensions
trymax = 30)
The clustering results were visualized in ordination spaced, using non-metric multidimensional scaling (NMDS) with Bray-Curtis distance ("vegan" package in R). Based on testing ordination stress for 1 to 5 dimensions, ordination stress is acceptable with r ord_dims dimensions, and stress reduction is minimal with additional dimensions. The stress for the ordination with r ord_dims dimensions is r round(ord$stress, 3).
# 2D plotting
pal_veg2 <- RColorBrewer::brewer.pal(n=k_groups_fuzzy, name = "Dark2")
# # manually update based on indicator species
# group_labels_fuzzy <- c(paste0("Juniper-pinyon
# woodland (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==1)), ")"),
# paste0("Open juniper-pinyon
# shrubland (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==2)), ")"),
# paste0("C4 grassland
# (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==3)), ")"),
# paste0("Annualized juniper
# savanna (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==4)), ")"),
# paste0("Pinyon-juniper
# woodland (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==5)), ")"),
# paste0("Gambel oak
# woodland (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==6)), ")"),
# paste0("Sagebrush shrubland/
# grassland (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==7)), ")"),
# paste0("Bare juniper
# savanna (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==3)), ")")
# )
group_labels_base <- c(paste0("Juniper-pinyon
woodland"),
paste0("Open juniper-pinyon
shrubland"),
paste0("C4 grassland
"),
paste0("Annualized juniper
savanna"),
paste0("Pinyon-juniper
woodland"),
paste0("Gambel oak
woodland"),
paste0("Sagebrush shrubland/
grassland"),
paste0("Bare juniper
savanna)")
)
group_labels_fuzzy <- paste0(group_labels_base,
" (n = ",
c(nrow(filter(plot_data_first, PlantCommunity_fuzzy==1)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==2)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==3)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==4)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==5)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==6)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==7)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==8))
),
")")
#group_labels_fuzzy <- as.character(1:clust_fuzz$k.crisp)
par(mfrow=c(2,2))
# Axes 1x2
plot(ord, choices = c(1,2), type = "n", # plot the axes
xlim = c(-1, 1),
ylim = c(-1, 1))
points(ord, choices = c(1,2), display = "sites", # plot points - can choose "sites" or "species"
col=pal_veg2[plot_data_first$PlantCommunity_fuzzy],
pch = 21, cex = .6, bg=pal_veg2[plot_data_first$PlantCommunity_fuzzy])
ordiellipse(ord, plot_data_first$PlantCommunity_fuzzy, col=pal_veg2, lwd = 2, label = F,
choices = c(1,2)) # plot your groups
# can use oriellipse, orihull, or ordispider to plot groups, depending on your needs
# Axes 3x2
plot(ord, choices = c(3,2), type = "n",
xlim = c(-1, 1),
ylim = c(-1, 1))
points(ord, choices = c(3,2), display = "sites",
col=pal_veg2[plot_data_first$PlantCommunity_fuzzy],
pch = 21, cex = .6, bg=pal_veg2[plot_data_first$PlantCommunity_fuzzy])
ordiellipse(ord, plot_data_first$PlantCommunity_fuzzy, col=pal_veg2, lwd = 2, label = F,
choices = c(3,2))
# Axes 1x3
plot(ord, choices = c(1,3), type = "n",
xlim = c(-1, 1),
ylim = c(-1, 1))
points(ord, choices = c(1,3), display = "sites",
col=pal_veg2[plot_data_first$PlantCommunity_fuzzy],
pch = 21, cex = .6, bg=pal_veg2[plot_data_first$PlantCommunity_fuzzy])
ordiellipse(ord, plot_data_first$PlantCommunity_fuzzy, col=pal_veg2, lwd = 2, label = F,
choices = c(1,3))
# Legend
plot(ord, type = "n", axes=FALSE,
display = "sites",
col=pal_veg2[plot_data_first$PlantCommunity_fuzzy],
xlab = "",
ylab = "")
legend(x="center",
legend = group_labels_fuzzy,
fill = pal_veg2,
title = "Plant communities",
bty = "n",
y.intersp = 1.5)
par(mfrow=c(1,1))
# save clustering results
fuzzy_cluster_membership <- as.data.frame(clust_fuzz$membership)
fuzzy_cluster_membership$PlotCode <- row.names(clust_fuzz$membership)
fuzzy_cluster_membership$best_group <- clust_fuzz$clustering
# update the cluster numbers to match the original groupings (see note above in the fuzzy-cluster chunk)
fuzzy_cluster_membership <- left_join(fuzzy_cluster_membership, updated_cluster_numbers, by = c("best_group" = "PlantCommunity_fuzzy_original"))
# # merge groups 1 and 5 based on workshop feedback
# fuzzy_cluster_membership <- mutate(fuzzy_cluster_membership,
# best_group = ifelse(best_group==1|best_group==5,
# yes = 9,
# no = best_group))
colnames(fuzzy_cluster_membership) <- c(paste0("Group_", 1:8, "_membership"), "PlotCode", "Best_group")
# group_names <- data.frame(Best_group = 1:8,
# Best_group_name = c("Juniper-pinyon woodland",
# "Open juniper-pinyon shrubland",
# "C4 grassland",
# "Annualized juniper savanna",
# "Pinyon-juniper woodland",
# "Gambel oak woodland",
# "Sagebrush shrubland/grassland",
# "Bare juniper savanna")) # update with descriptive names after examining the groups
# reorder group names due to user=Duniway_et_al_2024 mix up
group_names <- data.frame(Best_group = 1:8,
Best_group_name = c("Juniper-pinyon woodland",
"Annualized juniper savanna",
"Pinyon-juniper woodland",
"Open juniper-pinyon shrubland",
"Gambel oak woodland",
"Sagebrush shrubland/grassland",
"C4 grassland",
"Bare juniper savanna")) # update with descriptive names after examining the groups
# group_names <- data.frame(Best_group = 1:9,
# Best_group_name = c("Juniper-pinyon woodland",
# "Open juniper-pinyon shrubland",
# "C4 grassland",
# "Annualized juniper savanna",
# "Pinyon-juniper woodland",
# "Gambel oak woodland",
# "Sagebrush shrubland/grassland",
# "Bare juniper savanna",
# "Pinyon-juniper woodland"))
fuzzy_cluster_membership <- left_join(fuzzy_cluster_membership, group_names)
# get the plot location and date sampled (use first date)
plot_locations_first <- plot_locations %>%
mutate(DateSampled = as.Date(paste(Year, Month, Day, sep = "-"))) %>%
arrange(DateSampled) %>%
group_by(PlotCode) %>%
summarise(DateSampled = first(DateSampled))
fuzzy_cluster_memb_sf <- left_join(fuzzy_cluster_membership, plot_locations_first) %>%
filter(PlotCode != "AIM_Colorado NWDO Little Snake FO 2017_1024" &
PlotCode != "AIM_UT Color Country DO 2019_01") %>% # the geometry of these plot locations is messed up somehow
mutate(post_workshop1_group = ifelse(test = Best_group==1|Best_group==5,
yes = 9,
no = Best_group))
# # save plots - usually only need to run this once
# sf::st_write(fuzzy_cluster_memb_sf,
# dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps/contains_NRI_do_not_sync",
# layer = paste0("SW_ShallowDeepRocky_fuzzyclusters_", Sys.Date(), ".shp"),
# driver = "ESRI Shapefile")
#
# fuzzy_cluster_memb_sf_shareable <- filter(fuzzy_cluster_memb_sf, !grepl(pattern = "NRI", x=PlotCode))
#
# sf::st_write(fuzzy_cluster_memb_sf_shareable,
# dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps",
# layer = paste0("SW_ShallowDeepRocky_fuzzyclust_share_", Sys.Date(), ".shp"),
# driver = "ESRI Shapefile")
#
# states <- sf::st_read("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/GIS/Political boundaries/cb_2018_us_state_5m/cb_2018_us_state_5m.shp")
# fourcorners <- filter(states, NAME %in% c("Utah", "Colorado", "Arizona", "New Mexico"))
#
# fuzzy_cluster_memb_sf_mike <- sf::st_filter(sf::st_as_sf(fuzzy_cluster_memb_sf, crs=sf::st_crs(plot_locations)),
# sf::st_transform(fourcorners, crs=sf::st_crs(plot_locations))) %>%
# mutate(Source = ifelse(test = grepl(x=PlotCode, pattern = "^AIM"),
# yes = "AIM",
# no = NA)) %>%
# mutate(Source = ifelse(test = grepl(x=PlotCode, pattern = "^LMF"),
# yes = "LMF",
# no = Source)) %>%
# mutate(Source = ifelse(test = grepl(x=PlotCode, pattern = "^IM_NCPN"),
# yes = "NPS",
# no = Source)) %>%
# mutate(Source = ifelse(test = grepl(x=PlotCode, pattern = "^NRI"),
# yes = "NRI",
# no = Source))
#
# sf::st_write(fuzzy_cluster_memb_sf_mike,
# dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps/contains_NRI_do_not_sync",
# layer = paste0("SW_ShallowDeepRocky_fuzzyclust_mike_", Sys.Date(), ".shp"),
# driver = "ESRI Shapefile")
State and transition model {.tabset}
CUSTOM DIAGRAM
knitr::include_graphics("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Analyses/Semiarid_Warm_ShallowDeepRocky_boxandarrow_workshop1.JPG")
# TODO legend for custom diagram if needed
# scratch pad:
# All states described here contain some introduced species (e.g. *Bromus tectorum*; see Tables \@ref(tab:comm1-1-dominant-spp), \@ref(tab:comm2-1-dominant-spp), \@ref(tab:comm3-1-dominant-spp), \@ref(tab:comm4-1-dominant-spp), and \@ref(tab:comm5-1-dominant-spp) for details).
The data used to develop this STM were sourced from large monitoring data sets designed to sample landscapes with a variety of current and historic land use. The model does not describe reference communities because these data do not include targeted sampling of relict plant communities with histories of little alteration to disturbance regimes. Several states described here contain some introduced species (e.g. Bromus tectorum; see Tables \@ref(tab:comm3-1-dominant-spp) and \@ref(tab:comm4-1-dominant-spp) for details).
For furthur information on reference states as well as drivers of transitions, see the Ecological Site Descriptions for the sites listed in Table \@ref(tab:assoc-eco-sites).
STANDARD DIAGRAM
Click on state and transition labels to scroll to the respective text.
TODO create boxes for each state that link to the state description below
State descriptive figures
# load descriptive data
plot_data_descriptive <- plot_data_pull(user=user,
target_ESG = target_ESG,
data_sources = data_sources,
indicators = c("TotalFoliarCover",
"AH_WoodyCover",
"AH_ShrubCover",
# TODO what to do about subshrubs??
"AH_TreeCover",
"AH_ForbGrassCover",
"AH_PerenForbGrassCover",
"AH_AnnForbGrassCover",
"SoilStab_all",
"AH_ArtemisiaTridentataCover",
"CA_percent_200plus",
"CA_percent_100plus"
),
ann_grass_by_spp = T,
ann_forb_by_spp = T,
per_grass_by_spp = T,
per_forb_by_spp = T,
succulent_by_spp = T,
shrub_by_spp = T,
subshrub_by_spp = T,
tree_by_spp = T,
opuntia_combined = F,
impute_gap_type = c("CA_percent_100plus", "CA_percent_200plus")
) %>%
filter(!is.na(Month)) %>%
left_join(select(plot_data_first, PlotCode, Year, PlantCommunity_fuzzy) , .)
colnames(plot_data_descriptive)[which(colnames(plot_data_descriptive)=="AH_ArtemisiaTridentataCover")] <- "ARTR2"
# temporary patch - some plot codes are duplicated but have different primary keys
filter(plot_data_first,
PlotCode=="AIM_Colorado NWDO Little Snake FO 2017_1024" |
PlotCode=="AIM_UT Color Country DO 2019_01")$PlotID
plot_data_descriptive <- filter(plot_data_descriptive,
!(PlotCode=="AIM_Colorado NWDO Little Snake FO 2017_1024" & PlotID!="1705311149423160") &
!(PlotCode=="AIM_UT Color Country DO 2019_01" & PlotID!="1906201000185195"))
plot_locations <- cbind(plot_locations, sf::st_coordinates(plot_locations)) %>%
filter(!(PlotCode=="AIM_Colorado NWDO Little Snake FO 2017_1024" &
X!=filter(plot_data_descriptive, PlotCode=="AIM_Colorado NWDO Little Snake FO 2017_1024")$Longitude_NAD83)
& !(PlotCode=="AIM_UT Color Country DO 2019_01" &
X!=filter(plot_data_descriptive, PlotCode=="AIM_UT Color Country DO 2019_01")$Longitude_NAD83))
# summarize species occurance
descriptive.df <- dplyr::select(plot_data_descriptive, -SourceKey, -PlotID, -SiteName, -PlotName,
-Year, -Longitude_NAD83, -Latitude_NAD83,
-Month, -Day, -PlantCommunity_fuzzy, -TotalFoliarCover,
-AH_WoodyCover, -AH_ShrubCover, -AH_TreeCover, -AH_ForbGrassCover,
-AH_PerenForbGrassCover, -AH_AnnForbGrassCover,
-SoilStab_all, -CA_percent_200plus, -CA_percent_100plus
) %>%
tibble::column_to_rownames("PlotCode") # keep an ID code for the plot as the row name to prevent data scrambling problems
sp_pa_desc <- mutate(rowwise(descriptive.df), across(everything(), function(x){if(x>0){1}else{0}})) # make a presence/absence data frame
sp_keep_desc <- names(colSums(sp_pa_desc)[which(colSums(sp_pa_desc)>0)])
descriptive.df <- select(descriptive.df, all_of(sp_keep_desc))
plot_data_descriptive <- select(plot_data_descriptive, SourceKey, PlotID, SiteName, PlotName,
Year, PlotCode, Longitude_NAD83, Latitude_NAD83,
Month, Day, PlantCommunity_fuzzy, TotalFoliarCover,
AH_WoodyCover, AH_ShrubCover, AH_TreeCover, AH_ForbGrassCover,
AH_PerenForbGrassCover, AH_AnnForbGrassCover,
SoilStab_all, CA_percent_100plus, CA_percent_200plus, all_of(sp_keep_desc))
#ord.df <- select(ord.df, -any_of(c("UNKS", "UNKPG", "UNKAF", "UNKPF", "PG1", "PG01", "AF1")))
## load climate data
AI_raster <- raster::raster(file_paths$ai_raster)
plot_AI <- sf::st_as_sf(raster::extract(x = AI_raster,
y = plot_locations,
sp = T)) %>%
mutate(AI = ai_et0_WesternUSEcoregions*0.0001 # AI rasters are actually AI*10,000 so that they
#can be stored as integers - need to put back into standard AI units
) %>%
sf::st_drop_geometry()
rm(AI_raster)
climate_rasters <- raster::stack(file_paths$prism_rasters)
prism_vals <- raster::extract(x=climate_rasters,
y=plot_locations,
sp = T) %>%
as(., "sf")
sf::st_geometry(prism_vals) <- NULL
colnames(prism_vals)[which(colnames(prism_vals) %in% c("PRISM_ppt_30yr_normal_800mM2_annual_asc",
"PRISM_tmean_30yr_normal_800mM2_annual_asc",
"PRISM_tmax_30yr_normal_800mM2_annual_asc",
"PRISM_tmin_30yr_normal_800mM2_annual_asc"
))] <- c("Ann_ppt_mm", "Ann_tmean",
"Ann_tmax", "Ann_tmin")
rm(climate_rasters)
plot_data_descriptive <- left_join(plot_data_descriptive,
select(plot_AI, PlotCode, Year, AI)) %>%
left_join(.,
select(prism_vals, PlotCode, Year, Ann_ppt_mm, Ann_tmean, Ann_tmax, Ann_tmin))
## species richness and diversity
plot_data_descriptive$shannon_diversity <- vegan::diversity(x=descriptive.df, index = "shannon")
plot_data_descriptive$sp_richness <- vegan::specnumber(x=descriptive.df)
# # AERO
# aero <- read.csv(file_paths$aero)
# plot_aero <- aero %>%
# # AIM error patch
# filter(!(ProjectName=="Colorado NWDO Little Snake FO 2017" & PlotID=="1024" &
# Longitude_NAD83!=filter(plot_data_descriptive,
# PlotCode=="AIM_Colorado NWDO Little Snake FO 2017_1024")$Longitude_NAD83)
# & !(ProjectName=="UT Color Country DO 2019" & PlotID=="01" &
# Longitude_NAD83!=filter(plot_data_descriptive,
# PlotCode=="AIM_UT Color Country DO 2019_01")$Longitude_NAD83)) %>%
# select(ProjectName, PlotID, horizontal_flux_total_MN, horizontal_flux_total_LPI, horizontal_flux_total_UPI) %>%
# setNames(c("SiteName", "PlotName", "horizontal_flux_total_MN", "horizontal_flux_total_LPI", "horizontal_flux_total_UPI"))
#
# plot_data_descriptive <- left_join(plot_data_descriptive, plot_aero)
# AERO
aero <- read.csv(file_paths$aero)
aero_aim <- filter(aero, source=="AIM") %>%
mutate(PlotKey = substr(PrimaryKey, 1, nchar(PrimaryKey)-10), # fix the Jornada's data storage errors
Year = as.numeric(substr(PrimaryKey, nchar(PrimaryKey)-9, nchar(PrimaryKey)-6)),) %>%
select(PlotKey, Year, horizontal_flux_total_MN, horizontal_flux_total_LPI, horizontal_flux_total_UPI) %>%
rename("PlotID" = "PlotKey")
aero_lmf <- filter(aero, source=="LMF") %>%
select(PlotKey, horizontal_flux_total_MN, horizontal_flux_total_LPI, horizontal_flux_total_UPI) %>%
rename("PlotID" = "PlotKey")
plot_data_descriptive <- left_join(plot_data_descriptive, aero_aim) %>%
rows_update(x = ., y = aero_lmf, unmatched = "ignore")
# RHEM
rhem_aim <- read.csv(file_paths$rhem$aim) %>%
mutate(PrimaryKey_Jornada = as.character(PrimaryKey),
Year = as.numeric(substr(PrimaryKey, nchar(PrimaryKey)-9, nchar(PrimaryKey)-6)),
PlotID = substr(PrimaryKey, 1, nchar(PrimaryKey)-10)) %>% # make the Jornada's proprietary PrimaryKey system play nice with everyone else's
select(PlotID, Year, Runoff_Long_Term_MEAN, Soil_Loss_Long_Term_MEAN, Sed_Yield_Long_Term_MEAN) %>%
distinct()
rhem_lmf <- read.csv(file_paths$rhem$lmf) %>%
select(PrimaryKey, Runoff_Long_Term_MEAN, Soil_Loss_Long_Term_MEAN, Sed_Yield_Long_Term_MEAN) %>%
distinct() %>%
rename("PlotID" = "PrimaryKey")
rhem_nri <- read.csv(file_paths$rhem$nri) %>%
select(PrimaryKey, Runoff_Long_Term_MEAN, Soil_Loss_Long_Term_MEAN, Sed_Yield_Long_Term_MEAN) %>%
distinct() %>%
rename("PlotID" = "PrimaryKey")
plot_data_descriptive <- left_join(plot_data_descriptive, rhem_aim) %>%
rows_update(x = ., y = rhem_lmf, unmatched = "ignore") %>%
rows_update(x = ., y = rhem_nri, unmatched = "ignore")
# # consider including species present in >10% of plots:
# names(colSums(sp_pa)[which(colSums(sp_pa)>(nrow(plot_data_first)*0.10))])
# # also consider including species with high mean cover across all plots:
# mean_cover <- summarise(plot_data_first, across(.cols = AH_C3IntroducedPerenGrassCover:AH_OpuntiaCover, .fns = mean)) %>%
# select(-any_of(c("Longitude_NAD83", "Latitude_NAD83", "Year", "Month", "Day"))) %>%
# t() %>%
# as.data.frame()
#
# row.names(mean_cover)[which(mean_cover$V1>=0.9)]
# Consider switching to plotly to allow users to hover and see exact quantile and outlier values, e.g.:
# plotly::plot_ly(x=plot_data_descriptive$PlantCommunity_fuzzy,
# y=plot_data_descriptive$horizontal_flux_total_UPI,
# split=plot_data_descriptive$PlantCommunity_fuzzy,
# type='box')
plot_data_first_tall <- select(plot_data_first,
SourceKey, PlotID, SiteName,
PlotName, Year, PlotCode,
#PlantCommunity,
PlantCommunity_fuzzy,
AH_C3IntroducedPerenGrassCover,
AH_C3NativePerenGrassCover,
AH_C4NativePerenGrassCover,
AH_IntroducedAnnForbCover,
AH_IntroducedAnnGrassCover,
AH_IntroducedPerenForbCover,
AH_NativeAnnForbCover,
AH_NativeAnnGrassCover,
AH_NativePerenForbCover,
BareSoilCover,
FH_TotalLitterCover,
AH_ArtemisiaTridentataCover,
#CA_percent_200plus,
AH_C4IntroducedPerenGrassCover,
FH_LichenMossCover,
#CA_percent_100plus,
EPVI,
CEMO2,
CHVI8,
PUTR2,
AMUT,
QUGA,
GUSA2,
JUOS,
PIED
) %>%
left_join(., select(plot_data_descriptive, PlotCode, CA_percent_100plus, CA_percent_200plus)) %>%
tidyr::pivot_longer(cols = c(AH_C3IntroducedPerenGrassCover,
AH_C3NativePerenGrassCover,
AH_C4NativePerenGrassCover,
AH_IntroducedAnnForbCover,
AH_IntroducedAnnGrassCover,
AH_IntroducedPerenForbCover,
AH_NativeAnnForbCover,
AH_NativeAnnGrassCover,
AH_NativePerenForbCover,
BareSoilCover,
FH_TotalLitterCover,
AH_ArtemisiaTridentataCover,
CA_percent_200plus,
AH_C4IntroducedPerenGrassCover,
FH_LichenMossCover,
CA_percent_100plus,
EPVI,
CEMO2,
CHVI8,
PUTR2,
AMUT,
QUGA,
GUSA2,
JUOS,
PIED
),
names_to = "CoverType",
values_to = "PctCover")
spp_list <- read.csv(file_paths$species_list)
cover_funcgrps <- sort(c("AH_C3IntroducedPerenGrassCover",
"AH_C3NativePerenGrassCover",
"AH_C4NativePerenGrassCover",
"AH_IntroducedAnnForbCover",
"AH_IntroducedAnnGrassCover",
"AH_IntroducedPerenForbCover",
"AH_NativeAnnForbCover",
"AH_NativeAnnGrassCover",
"AH_NativePerenForbCover",
"BareSoilCover",
"FH_TotalLitterCover",
"AH_ArtemisiaTridentataCover",
"CA_percent_200plus",
"AH_C4IntroducedPerenGrassCover",
"FH_LichenMossCover",
"CA_percent_100plus"))
funcgrps_names <- indicator_descriptions %>%
bind_rows(., data.frame(Indicator_code = c("CA_percent_100plus", "CA_percent_200plus"),
Indicator = c("Ann. & peren. canopy gaps > 100 cm", "Ann. & peren. canopy gaps > 200 cm"))) %>%
filter(., Indicator_code %in% cover_funcgrps) %>%
arrange(., Indicator_code) %>%
select(Indicator) %>%
unlist()
names(cover_funcgrps) <- funcgrps_names
cover_spp <- sort(c("EPVI",
"CEMO2",
"CHVI8",
"PUTR2",
"AMUT",
"QUGA",
"GUSA2",
"JUOS",
"PIED"
))
spp_names <- paste0( filter(spp_list, SpeciesCode %in% cover_spp & SpeciesState=="AIM")$ScientificName,
" (", filter(spp_list, SpeciesCode %in% cover_spp & SpeciesState=="AIM")$CommonName, ")")
names(cover_spp) <- spp_names
cover_labels <- sort(c(cover_funcgrps, cover_spp))
functypes_grps <- c("TotalFoliarCover", #"AH_WoodyCover", #AH_ForbGrassCover,
"AH_ShrubCover", "AH_TreeCover",
"AH_PerenForbGrassCover", "AH_AnnForbGrassCover"
)
functypes_cover_states <- ggplot(data = select(plot_data_descriptive,
any_of(c("PlantCommunity_fuzzy", functypes_grps))) %>%
tidyr::pivot_longer(cols = all_of(functypes_grps),
names_to = "CoverType",
values_to = "PctCover"),
aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) +
geom_boxplot() +
scale_fill_brewer(type = "qual", palette = "Set3",
name = "Functional type",
labels = c("Annual herbaceous", "Perennial herbaceous",
#"Woody",
"Shrub", "Tree",
"Total foliar cover")
) +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab(NULL) +
ylab("Cover (%)") +
ggtitle("Plant functional groups") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black"),
axis.text.y = element_text(color = "black"))
woody_grps <- c("AH_ArtemisiaTridentataCover",
"EPVI",
"CEMO2",
"CHVI8",
"PUTR2",
"AMUT",
"QUGA",
"GUSA2",
"JUOS",
"PIED"
)
woody_cover_states <- ggplot(data = filter(plot_data_first_tall, CoverType %in% woody_grps),
aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) +
geom_boxplot() +
scale_fill_brewer(type = "qual", palette = "Set3",
name = "Species",
labels = names(cover_labels[which(cover_labels %in% woody_grps)])) +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab(NULL) +
ylab("Cover (%)") +
ggtitle("Woody plants") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black"),
axis.text.y = element_text(color = "black"))
perenherb_grps <- c("AH_C3IntroducedPerenGrassCover", "AH_C3NativePerenGrassCover",
"AH_C4IntroducedPerenGrassCover", "AH_C4NativePerenGrassCover",
"AH_IntroducedPerenForbCover", "AH_NativePerenForbCover")
perenherb_cover_states <- ggplot(data = filter(plot_data_first_tall, CoverType %in% perenherb_grps),
aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) +
geom_boxplot() +
scale_fill_brewer(type = "qual", palette = "Set3",
name = "Functional group",
labels = names(cover_labels[which(cover_labels %in% perenherb_grps)])) +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab(NULL) +
ylab("Cover (%)") +
ggtitle("Perennial herbaceous plants") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black"),
axis.text.y = element_text(color = "black"))
annherb_grps <- c("AH_IntroducedAnnForbCover", "AH_IntroducedAnnGrassCover",
"AH_NativeAnnForbCover", "AH_NativeAnnGrassCover")
annherb_cover_states <- ggplot(data = filter(plot_data_first_tall, CoverType %in% annherb_grps),
aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) +
geom_boxplot() +
scale_fill_brewer(type = "qual", palette = "Set3",
name = "Functional group",
labels = names(cover_labels[which(cover_labels %in% annherb_grps)])) +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab(NULL) +
ylab("Cover (%)") +
ggtitle("Annual herbaceous plants") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black"),
axis.text.y = element_text(color = "black"))
ground_grps <- c("BareSoilCover", "FH_TotalLitterCover", "CA_percent_200plus",
"FH_LichenMossCover", "CA_percent_100plus")
ground_cover_states <- ggplot(data = filter(plot_data_first_tall, CoverType %in% ground_grps),
aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) +
geom_boxplot() +
scale_fill_brewer(type = "qual", palette = "Set3",
name = "Functional group",
labels = names(cover_labels[which(cover_labels %in% ground_grps)])) +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab("Plant community") +
ylab("Cover (%)") +
ggtitle("Ground cover") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black"),
axis.text.y = element_text(color = "black"))
cover_states <- cowplot::plot_grid(functypes_cover_states + theme(legend.position = "none"),
woody_cover_states + theme(legend.position = "none"),
perenherb_cover_states + theme(legend.position = "none"),
annherb_cover_states + theme(legend.position = "none"),
ground_cover_states + theme(legend.position = "none"),
ncol = 1, axis = "lr", labels = "auto")
cover_axes <- cowplot::plot_grid(cowplot::get_legend(functypes_cover_states),
cowplot::get_legend(woody_cover_states),
cowplot::get_legend(perenherb_cover_states),
cowplot::get_legend(annherb_cover_states),
cowplot::get_legend(ground_cover_states),
ncol = 1)
cowplot::plot_grid(cover_states, cover_axes, ncol = 2, rel_widths = c(1, 0.5))
group_labels_fuzzy <- paste0(group_labels_base,
" (n = ",
c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1 & !is.na(AI))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2 & !is.na(AI))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3 & !is.na(AI))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4 & !is.na(AI)))
),
")")
cowplot::plot_grid(
cowplot::plot_grid(
ggplot(data = plot_data_descriptive,
aes(x=PlantCommunity_fuzzy, y=AI)) +
geom_boxplot() +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab(NULL) +
ylab("Aridity index") +
ylim(c(0, 0.45)) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)),
ggplot(data = plot_data_descriptive,
aes(x=PlantCommunity_fuzzy, y=Ann_ppt_mm)) +
geom_boxplot() +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab(NULL) +
ylab("Annual precipitation (mm)") +
ylim(c(0, 650)) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)),
nrow = 1, labels = c("a", "b")),
ggplot(data = select(plot_data_descriptive, PlantCommunity_fuzzy, Ann_tmean, Ann_tmax, Ann_tmin) %>%
tidyr::pivot_longer(c(Ann_tmean, Ann_tmax, Ann_tmin), names_to = "temptype", values_to = "tempdegrees"),
aes(fill=PlantCommunity_fuzzy, y=tempdegrees, x=temptype)) +
geom_boxplot() +
scale_fill_manual(name="Plant community", labels=group_labels_fuzzy, values = pal_veg2) +
scale_x_discrete(labels=c("Maximum", "Mean", "Minimum")) +
xlab(NULL) +
ylab("Temperature (C)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)),
nrow = 2, labels = c(NA, "c"))
group_labels_fuzzy <- paste0(group_labels_base,
" (n = ",
c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1 & !is.na(SoilStab_all))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2 & !is.na(SoilStab_all))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3 & !is.na(SoilStab_all))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4 & !is.na(SoilStab_all)))
),
")")
ggplot(data = plot_data_descriptive,
aes(x=PlantCommunity_fuzzy, y=SoilStab_all)) +
geom_boxplot() +
scale_x_discrete(labels=group_labels_fuzzy) +
ylim(c(0, 6)) +
xlab("Plant community") +
ylab("Mean soil stability rating") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
group_labels_fuzzy <- paste0(group_labels_base,
" (n = ",
c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1 & !is.na(shannon_diversity))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2 & !is.na(shannon_diversity))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3 & !is.na(shannon_diversity))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4 & !is.na(shannon_diversity)))
),
")")
cowplot::plot_grid(ggplot(data = plot_data_descriptive,
aes(x=PlantCommunity_fuzzy, y=shannon_diversity)) +
geom_boxplot() +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab("Plant community") +
ylab("Shannon diversity index") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)),
ggplot(data = plot_data_descriptive,
aes(x=PlantCommunity_fuzzy, y=sp_richness)) +
geom_boxplot() +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab("Plant community") +
ylab("Number of species") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)),
nrow = 1, labels = "auto")
# dust_tall <- plot_data_descriptive %>%
# select(PlotID, PlantCommunity_fuzzy, horizontal_flux_total_MN, horizontal_flux_total_LPI, horizontal_flux_total_UPI) %>%
# tidyr::pivot_longer(cols = c(horizontal_flux_total_MN, horizontal_flux_total_LPI, horizontal_flux_total_UPI),
# names_to = "dust_statistic",
# values_to = "dust_flux")
group_labels_fuzzy <- paste0(group_labels_base,
" (n = ",
c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1 & !is.na(horizontal_flux_total_MN))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2 & !is.na(horizontal_flux_total_MN))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3 & !is.na(horizontal_flux_total_MN))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4 & !is.na(horizontal_flux_total_MN)))
),
")")
ggplot(data = plot_data_descriptive,
aes(x=PlantCommunity_fuzzy, y=horizontal_flux_total_MN)) +
geom_boxplot() +
scale_x_discrete(labels=group_labels_fuzzy) +
#ylim(c(0, 12)) +
xlab("Plant community") +
ylab("Mean horizontal dust flux") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# TODO add log transformation to y axis
group_labels_fuzzy <- paste0(group_labels_base,
" (n = ",
c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1 & !is.na(Runoff_Long_Term_MEAN))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2 & !is.na(Runoff_Long_Term_MEAN))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3 & !is.na(Runoff_Long_Term_MEAN))),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4 & !is.na(Runoff_Long_Term_MEAN)))
),
")")
cowplot::plot_grid(ggplot(data = plot_data_descriptive,
aes(x=PlantCommunity_fuzzy, y=Runoff_Long_Term_MEAN)) +
geom_boxplot() +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab("Plant community") +
ylab("Long-term mean runoff (mm/yr)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)),
ggplot(data = plot_data_descriptive,
aes(x=PlantCommunity_fuzzy, y=Soil_Loss_Long_Term_MEAN)) +
geom_boxplot() +
scale_x_discrete(labels=group_labels_fuzzy) +
xlab("Plant community") +
ylab("Long-term mean soil loss (ton/ha/yr)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)),
nrow = 1, labels = "auto")
State 1
(No name) State
TODO write state description
Community 1.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp1 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==1),
species_cols = colnames(descriptive.df),
user = user)
# save the complete table for field crews
write.csv(x = dom_spp_grp1,
file = file.path(output_figure_folder, "Tables",
paste0(target_ESG, "_1Juniper-pinyonwoodland_species.csv")),
row.names = F)
dom_spp_grp1_pub <- filter(
select(dom_spp_grp1, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots),
pct_plots >=10) %>%
mutate(MeanCover = round(MeanCover, 2),
pct_plots = round(pct_plots, 2)) %>%
setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present"))
format_tables_EDIT_style(data = dom_spp_grp1_pub,
caption = "Community 1.1 plant community composition") %>%
collapse_rows(columns = 1, valign = "top")
# TODO create pie charts of community composition by functional group. May need to change "production" to "percent cover" to fit with our data
# TODO create table of production (or percent cover) by functional type
# knitr::kable(x = dom_plant_spp_tab, caption = "Table 5. Annual production by plant type")
# TODO create table of ground cover by type
# knitr::kable(x = dom_plant_spp_tab, caption = "Table 6. Ground cover")
# TODO create table of func type percent cover at different heights
# knitr::kable(x = dom_plant_spp_tab, caption = "Table 7. Canopy structure (% cover)")
State 2
(No name) State
Community 2.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp2 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==2),
species_cols = colnames(descriptive.df),
user = user)
# save the complete table for field crews
write.csv(x = dom_spp_grp2,
file = file.path(output_figure_folder, "Tables",
paste0(target_ESG, "_2Openjuniper-pinyonshrubland_species.csv")),
row.names = F)
dom_spp_grp2_pub <- filter(
select(dom_spp_grp2, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots),
pct_plots >=10) %>%
mutate(MeanCover = round(MeanCover, 2),
pct_plots = round(pct_plots, 2)) %>%
setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present"))
format_tables_EDIT_style(data = dom_spp_grp2_pub,
caption = "Community 2.1 plant community composition") %>%
collapse_rows(columns = 1, valign = "top")
State 3
(No name) State
TODO write state description
Community 3.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp3 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==3),
species_cols = colnames(descriptive.df),
user = user)
# save the complete table for field crews
write.csv(x = dom_spp_grp3,
file = file.path(output_figure_folder, "Tables",
paste0(target_ESG, "_3C4grassland_species.csv")),
row.names = F)
dom_spp_grp3_pub <- filter(
select(dom_spp_grp3, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots),
pct_plots >=10) %>%
mutate(MeanCover = round(MeanCover, 2),
pct_plots = round(pct_plots, 2)) %>%
setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present"))
format_tables_EDIT_style(data = dom_spp_grp3_pub,
caption = "Community 3.1 plant community composition") %>%
collapse_rows(columns = 1, valign = "top")
State 4
(No name) State
TODO write state description
Community 4.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp4 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==4),
species_cols = colnames(descriptive.df),
user = user)
# save the complete table for field crews
write.csv(x = dom_spp_grp4,
file = file.path(output_figure_folder, "Tables",
paste0(target_ESG, "_4Annualizedjunipersavanna_species.csv")),
row.names = F)
dom_spp_grp4_pub <- filter(
select(dom_spp_grp4, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots),
pct_plots >=10) %>%
mutate(MeanCover = round(MeanCover, 2),
pct_plots = round(pct_plots, 2)) %>%
setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present"))
format_tables_EDIT_style(data = dom_spp_grp4_pub,
caption = "Community 4.1 plant community composition") %>%
collapse_rows(columns = 1, valign = "top")
Additional community tables
# TODO create table of production (or percent cover) by species, with species linking to their USDA Plants profile
# knitr::kable(x = dom_plant_spp_tab, caption = "Table 8. Community 1.1 plant community composition")
State 5
(No name) State
TODO write state description
Community 5.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp5 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==5),
species_cols = colnames(descriptive.df),
user = user)
# save the complete table for field crews
write.csv(x = dom_spp_grp5,
file = file.path(output_figure_folder, "Tables",
paste0(target_ESG, "_5Pinyon-juniperwoodland_species.csv")),
row.names = F)
dom_spp_grp5_pub <- filter(
select(dom_spp_grp5, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots),
pct_plots >=10) %>%
mutate(MeanCover = round(MeanCover, 2),
pct_plots = round(pct_plots, 2)) %>%
setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present"))
format_tables_EDIT_style(data = dom_spp_grp5_pub,
caption = "Community 5.1 plant community composition") %>%
collapse_rows(columns = 1, valign = "top")
State 6
(No name) State
TODO write state description
Community 6.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp6 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==6),
species_cols = colnames(descriptive.df),
user = user)
# save the complete table for field crews
write.csv(x = dom_spp_grp6,
file = file.path(output_figure_folder, "Tables",
paste0(target_ESG, "_6Gambeloakwoodland_species.csv")),
row.names = F)
dom_spp_grp6_pub <- filter(
select(dom_spp_grp6, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots),
pct_plots >=10) %>%
mutate(MeanCover = round(MeanCover, 2),
pct_plots = round(pct_plots, 2)) %>%
setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present"))
format_tables_EDIT_style(data = dom_spp_grp6_pub,
caption = "Community 6.1 plant community composition") %>%
collapse_rows(columns = 1, valign = "top")
Additional community tables
# TODO create table of production (or percent cover) by species, with species linking to their USDA Plants profile
# knitr::kable(x = dom_plant_spp_tab, caption = "Table 8. Community 1.1 plant community composition")
State 7
(No name) State
TODO write state description
Community 7.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp7 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==7),
species_cols = colnames(descriptive.df),
user = user)
# save the complete table for field crews
write.csv(x = dom_spp_grp7,
file = file.path(output_figure_folder, "Tables",
paste0(target_ESG, "_7Sagebrushshrublandgrassland_species.csv")),
row.names = F)
dom_spp_grp7_pub <- filter(
select(dom_spp_grp7, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots),
pct_plots >=10) %>%
mutate(MeanCover = round(MeanCover, 2),
pct_plots = round(pct_plots, 2)) %>%
setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present"))
format_tables_EDIT_style(data = dom_spp_grp7_pub,
caption = "Community 7.1 plant community composition") %>%
collapse_rows(columns = 1, valign = "top")
Additional community tables
# TODO create table of production (or percent cover) by species, with species linking to their USDA Plants profile
# knitr::kable(x = dom_plant_spp_tab, caption = "Table 8. Community 1.1 plant community composition")
State 8
(No name) State
TODO write state description
Community 8.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp8 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==8),
species_cols = colnames(descriptive.df),
user = user)
# save the complete table for field crews
write.csv(x = dom_spp_grp8,
file = file.path(output_figure_folder, "Tables",
paste0(target_ESG, "_8Barejunipersavanna_species.csv")),
row.names = F)
dom_spp_grp8_pub <- filter(
select(dom_spp_grp8, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots),
pct_plots >=10) %>%
mutate(MeanCover = round(MeanCover, 2),
pct_plots = round(pct_plots, 2)) %>%
setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present"))
format_tables_EDIT_style(data = dom_spp_grp8_pub,
caption = "Community 8.1 plant community composition") %>%
collapse_rows(columns = 1, valign = "top")
Additional community tables
# TODO create table of production (or percent cover) by species, with species linking to their USDA Plants profile
# knitr::kable(x = dom_plant_spp_tab, caption = "Table 8. Community 1.1 plant community composition")
Interpretations
Animal community
TODO write wildlife narrative
Hydrological functions
TODO write hydrology narrative
Recreational uses
TODO write recreation narrative
Wood products
TODO write wood production narrative
Other products
TODO grazing narrative could go here
Supporting information
Other references
Herrick, J. E., Van Zee, J. W., McCord, S. E., Courtright, E. M., Karl, J. W., & Burkett, L. M. (2017). Monitoring Manual for Grassland, Shrubland, and Savanna Ecosystems (Second ed. Vol. 1: Core Methods). Las Cruces, New Mexico: USDA-ARS Jornada Experimental Range.
Maechler, M., Rousseeuw, P., Struyf, A., Hubert, M., Hornik, K.(2021). cluster: Cluster Analysis Basics and Extensions. R package version 2.1.2.
McCord, S. E., Brehm, J. R., Burnett, S. H., Dietrich, C., Edwards, B., Metz, L. J., Hernandez Narvaez, M., Pierson, F., Ramirez, K. S., Stauffer, N. G., Webb, N. P., and Tweedie, C. E. (2022). A framework and toolset for standardizing agroecosystem indicators. Ecological Indicators, 144. doi:10.1016/j.ecolind.2022.109511
Nauman, T. W., Burch, S. S., Humphries, J. T., Knight, A. C., & Duniway, M. C. (2022). A Quantitative Soil-Geomorphic Framework for Developing and Mapping Ecological Site Groups. Rangeland Ecology & Management, 81, 9-33. doi:10.1016/j.rama.2021.11.003
Oksanen, J., F. Guillaume Blanchet, Michael Friendly, Roeland Kindt, Pierre Legendre, Dan McGlinn, Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens, Eduard Szoecs and Helene Wagner (2020). vegan: Community Ecology Package. R package version 2.5-7. https://CRAN.R-project.org/package=vegan
Contributors
- Anna Knight
- Travis Nauman
- Samuel Burch
- Mike Duniway
Reference sheet
Interpreting Indicators of Rangeland Health is a qualitative assessment protocol used to determine ecosystem condition based on benchmark characteristics described in the Reference Sheet. A suite of 17 (or more) indicators are typically considered in an assessment. The ecological site(s) representative of an assessment location must be known prior to applying the protocol and must be verified based on soils and climate. Current plant community cannot be used to identify the ecological site.
# TODO create table with author, author contact info, date, approver, approved date, and "Composition (Indicators 10 and 12) based on" info for the IIRH reference sheet
# knitr::kable(x = dom_plant_spp_tab)
Indicators
TODO fill in the 17 IIRH indicators for the reference state if available
TODO link to PDF of the reference sheet
Print options
TODO look into options for printing from the RMarkdown HTML format
# 1. plant functional group join plot with ordination
# 2. lattice of abiotic gradients that distinguish the states
# aero
# rhem if available
# soil stability
# climate?
# 3. NPP
# 4. biodiversity
# 5 remote sensing variance
# 6. Functional type cover grouped by functional type and colored by state
# 7. Disturbance and land use
# MTBS
# Grazing?
# PJ treatments
# SWreGAP protection classes
#group_labels_short <- c("Grassland", "Perennial grass loss", "Shrubland", "Mesic shrubland", "Invaded")
plot_data_descriptive$PlantCommunity_fuzzy_reordered <- factor(plot_data_descriptive$PlantCommunity_fuzzy,
levels = c(as.character(1:8)),
labels = c("Juniper-pinyon woodland",
"Open juniper-pinyon shrubland",
"C4 grassland",
"Annualized juniper savanna",
"Pinyon-juniper woodland",
"Gambel oak woodland",
"Sagebrush shrubland/grassland",
"Bare juniper savanna"),
ordered = T)
pal_veg2_reordered <- pal_veg2
# 1. plant functional group join plot with ordination
species_scores <- as.data.frame(scores(ord, choices = 1:3, display = "species")) # get axis scores for each response type
plot_scores <- as.data.frame(scores(ord, choices = 1:3, display = "sites")) # get the axis scores for each site
plot_scores$PlotCode <- rownames(plot_scores)
plot_scores <- left_join(plot_scores, select(plot_data_descriptive, PlotCode, PlantCommunity_fuzzy, PlantCommunity_fuzzy_reordered))
# this bit gets the arrow length and direction info
ord.fit <- envfit(ord ~ AH_PerenForbGrassCover +
AH_AnnForbGrassCover +
AH_ShrubCover +
AH_TreeCover +
TotalFoliarCover,
data = plot_data_descriptive,
na.rm=T,
choices=1:3)
vect_functypes <- as.data.frame(ord.fit$vectors$arrows*sqrt(ord.fit$vectors$r))
vect_functypes$var <- rownames(vect_functypes)
vect_functypes$r <- ord.fit$vectors$r
vect_functypes$varnames <- c("Perennial herbaceous",
"Annual herbaceous",
"Shrub",
"Tree",
"Total foliar")
# # since I have 3 ordination axes, I'm filtering just to axis 1 and 2 for this plot
vect_functypes_12 <- filter(vect_functypes, abs(NMDS3)<abs(NMDS1)|abs(NMDS3)<abs(NMDS2))
vect_functypes_32 <- filter(vect_functypes, abs(NMDS1)<abs(NMDS3)|abs(NMDS1)<abs(NMDS2))
vect_functypes_13 <- filter(vect_functypes, abs(NMDS2)<abs(NMDS1)|abs(NMDS2)<abs(NMDS3))
# make pretty plot ggplot!
ord_pretty_plot <- function(points = plot_scores,
vectors = vect_functypes,
x_axis = "NMDS1",
y_axis = "NMDS2",
palette = pal_veg2_reordered,
#group_labels = group_labels_short,
labels = T){
ord_plot <- ggplot(data = points, aes_string(x=x_axis, y=y_axis)) +
geom_point(data = points, aes_string(x=x_axis, y=y_axis, color="PlantCommunity_fuzzy_reordered"), size=.5) +
stat_ellipse(data=points, aes_string(x=x_axis, y=y_axis, color="PlantCommunity_fuzzy_reordered"),
level = 0.5, lwd = 1) +
scale_color_manual(values =palette, name="Plant community"#, # set custom colors for ellipses
#labels=group_labels
) +
geom_segment(data = vectors, # make the arrows with environmental vars
aes_string(x="0", xend=x_axis, y="0", yend=y_axis),
arrow = arrow(length = unit(0.25, "cm")), lwd=1, colour="black", inherit.aes = F) +
coord_fixed() +
lims(x=c(-1.25,1.25), y=c(-1.25,1.25)) +
theme_bw() +
theme(axis.text = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
if(labels){
ord_plot <- ord_plot + ggrepel::geom_label_repel(data = vectors, # put labels on the arrows
aes_string(x=x_axis, y=y_axis, label="varnames"), size=3, colour="black")
}
return(ord_plot)
}
ord_plot12 <- ord_pretty_plot(x_axis = "NMDS1", y_axis = "NMDS2", vectors = vect_functypes_12)
ord_plot32 <- ord_pretty_plot(x_axis = "NMDS3", y_axis = "NMDS2", vectors = vect_functypes_32)
ord_plot13 <- ord_pretty_plot(x_axis = "NMDS1", y_axis = "NMDS3", vectors = vect_functypes_13)
ord_plot_legend <- cowplot::get_legend(ord_plot12)
ord_plot_all <- cowplot::plot_grid(ord_plot12 + theme(legend.position = "none"), ord_plot32 + theme(legend.position = "none"),
ord_plot13 + theme(legend.position = "none"), ord_plot_legend,
nrow = 2, labels = c("a", "b", "c"))
today <- Sys.Date()
highrestiff(plot_obj = ord_plot_all,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Ordination/SW_ShallowDeepRocky_NMDS_ellipses_functypes_", today, ".tif"),
width_in = 6,
height_in = 6,
resolution_dpi = 300)
ord_plot12_nolabels <- ord_pretty_plot(x_axis = "NMDS1", y_axis = "NMDS2", vectors = vect_functypes_12, labels = F)
ord_plot32_nolabels <- ord_pretty_plot(x_axis = "NMDS3", y_axis = "NMDS2", vectors = vect_functypes_32, labels = F)
ord_plot13_nolabels <- ord_pretty_plot(x_axis = "NMDS1", y_axis = "NMDS3", vectors = vect_functypes_13, labels = F)
ord_plot_legend_nolabels <- cowplot::get_legend(ord_plot12_nolabels)
ord_plot_all_nolabels <- cowplot::plot_grid(ord_plot12_nolabels + theme(legend.position = "none"), ord_plot32_nolabels + theme(legend.position = "none"),
ord_plot13_nolabels + theme(legend.position = "none"), ord_plot_legend_nolabels,
nrow = 2, labels = c("a", "b", "c"))
today <- Sys.Date()
highrestiff(plot_obj = ord_plot_all_nolabels,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Ordination/SW_ShallowDeepRocky_NMDS_ellipses_functypes_nolabels_", today, ".tif"),
width_in = 6,
height_in = 6,
resolution_dpi = 300)
# 2. lattice of abiotic gradients that distinguish the states
violin_pretty_plot <- function(dat = plot_data_descriptive,
y_var = "horizontal_flux_total_MN",
y_lab = "Total horizontal dust flux",
fill_var = NULL,
fill_name = "Plant community",
scaling = "area"){
if(!is.null(fill_var)){
violin_plot <- ggplot(data = dat,
aes_string(x="PlantCommunity_fuzzy_reordered", y=y_var, fill = fill_var)) +
scale_fill_manual(values =pal_veg2_reordered, name= fill_name, # set custom colors for ellipses
#labels=group_labels_short
)
}else{
violin_plot <- ggplot(data = dat,
aes_string(x="PlantCommunity_fuzzy_reordered", y=y_var))
}
violin_plot <- violin_plot +
geom_violin(scale = scaling, trim = F, na.rm = T) +
geom_boxplot(width=0.1) +
#scale_x_discrete(labels=group_labels_short) +
xlab(NULL) +
ylab(y_lab) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"),
axis.text.y = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none")
return(violin_plot)
}
box_pretty_plot <- function(dat = plot_data_descriptive,
y_var = "horizontal_flux_total_MN",
y_lab = "Total horizontal dust flux",
fill_var = "PlantCommunity_fuzzy_reordered",
palette = pal_veg2_reordered){
group_labels <- paste0(as.character(levels(dat$PlantCommunity_fuzzy_reordered)),
"
(n = ",
c(nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==1 & !is.na(dat[[y_var]]), ]),
nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==2 & !is.na(dat[[y_var]]), ]),
nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==3 & !is.na(dat[[y_var]]), ]),
nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==4 & !is.na(dat[[y_var]]), ]),
nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==5 & !is.na(dat[[y_var]]), ]),
nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==6 & !is.na(dat[[y_var]]), ]),
nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==7 & !is.na(dat[[y_var]]), ]),
nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==8 & !is.na(dat[[y_var]]), ])
),
")")
box_plot <- ggplot(data = dat,
aes_string(x="PlantCommunity_fuzzy_reordered", y=y_var, fill = fill_var)) +
geom_boxplot() +
scale_x_discrete(labels=group_labels) +
scale_fill_manual(values =palette, name="Plant community", # set custom colors for ellipses
#labels=group_labels_short
) +
xlab(NULL) +
ylab(y_lab) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"),
axis.text.y = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none")
return(box_plot)
}
# aero
aero_plot <- box_pretty_plot(dat = filter(plot_data_descriptive, !is.na(horizontal_flux_total_MN)),
y_var = "horizontal_flux_total_MN",
y_lab = "Horizontal dust flux
(g/m*day)"#,
#scaling = "width"
) +
scale_y_log10(breaks = scales::trans_breaks("log10", function(x) 10^x),
labels = scales::trans_format("log10", scales::label_math(10^.x))) #+
#annotation_logticks(sides = "l")
aero_plot
# soil stability
soilstab_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "SoilStab_all",
y_lab = "Soil stability")
soilstab_plot
# climate?
AI_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "AI",
y_lab = "Aridity index")
AI_plot
ppt_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "Ann_ppt_mm",
y_lab = "Mean ann. ppt. (mm)")
ppt_plot
tmean_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "Ann_tmean",
y_lab = "Mean ann. temp. (C)")
tmean_plot
# rhem
rhem_plot_runoff <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "Runoff_Long_Term_MEAN",
y_lab = "Mean runoff
(mm/yr)")
rhem_plot_runoff
rhem_plot_soilloss <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "Soil_Loss_Long_Term_MEAN",
y_lab = "Mean soil loss
(tons/ha/yr)")
rhem_plot_soilloss
# put them all together
abiotic_plots <- cowplot::plot_grid(aero_plot,
soilstab_plot,
rhem_plot_runoff,
rhem_plot_soilloss,
ppt_plot, tmean_plot,
nrow = 3, #rel_heights = c(0.6, 1),
align = "v",
labels = "auto")
abiotic_plots
today <- Sys.Date()
highrestiff(plot_obj = abiotic_plots,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_abiotic_boxplot_2_", today, ".tif"),
width_in = 6,
height_in = 9,
resolution_dpi = 300)
# put them in a different layout because we love drafts
climate_plots <- cowplot::plot_grid(ppt_plot, tmean_plot, AI_plot,
nrow = 1,
labels = "auto")
climate_plots
highrestiff(plot_obj = climate_plots,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_climate_boxplot_", today, ".tif"),
width_in = 9,
height_in = 3.5,
resolution_dpi = 300)
soils_plots <- cowplot::plot_grid(aero_plot, soilstab_plot,
rhem_plot_runoff, rhem_plot_soilloss,
nrow = 2,
labels = "auto")
soils_plots
highrestiff(plot_obj = soils_plots,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_soils_boxplot_", today, ".tif"),
width_in = 6,
height_in = 6,
resolution_dpi = 300)
# 3. NPP
# NPP data is from the RAP 3.0 partitioned NPP data set (Robinson et al., 2019). I downloaded the predicted NPP values within a 35-m radius of the AIM, LMF, NRI, and NCPN plot centers using a Google Earth Engine script. If we continue using NPP for other ESGs, we will probably need to continue downloading point values for RAP products through GEE - the raster file sizes are enormous!
# Robinson, N. P., Allred, B. W., Naugle, D. E., & Jones, M. O. (2019). Patterns of rangeland productivity and land ownership: Implications for conservation and management. Ecologial Applications, 23(3), e01862.
NPP_files <- list.files(paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/RAP_do_not_sync/", target_ESG, "/NPP"),
full.names = T)
npp <- read.csv(NPP_files[1])
npp$Year <- as.numeric(substr(NPP_files[1], 171, 174))
for(file in NPP_files[-1]){
npp_temp <- read.csv(file)
npp_temp$Year <- as.numeric(substr(file, 171, 174))
npp <- bind_rows(npp, npp_temp)
rm(npp_temp)
}
npp <- select(npp, PlotCod, Year, afgNPP, pfgNPP, shrNPP, treNPP) %>%
setNames(c("PlotCode", "Year", "afgNPP", "pfgNPP", "shrNPP", "treNPP"))
plot_data_descriptive <- left_join(plot_data_descriptive, npp)
afg_npp_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "afgNPP",
y_lab = "Ann. herb. NPP (g C/m2)")
afg_npp_plot
pfg_npp_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "pfgNPP",
y_lab = "Per. herb. NPP (g C/m2)")
pfg_npp_plot
shr_npp_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "shrNPP",
y_lab = "Shrub NPP (g C/m2)")
shr_npp_plot
tre_npp_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "treNPP",
y_lab = "Tree NPP (g C/m2)")
tre_npp_plot
pft_npp_plots <- cowplot::plot_grid(afg_npp_plot + xlab(NULL) + scale_x_discrete(labels=NULL), pfg_npp_plot + xlab(NULL) + scale_x_discrete(labels=NULL),
shr_npp_plot, tre_npp_plot,
nrow = 2, rel_heights = c(0.6, 1), align = "v",
labels = "auto")
pft_npp_plots
today <- Sys.Date()
highrestiff(plot_obj = pft_npp_plots,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_pftnpp_boxplot_", today, ".tif"),
width_in = 6,
height_in = 6.5,
resolution_dpi = 300)
# combining npp to herbaceous, woody, and total
plot_data_descriptive$herbNPP <- plot_data_descriptive$afgNPP + plot_data_descriptive$pfgNPP
plot_data_descriptive$woodyNPP <- plot_data_descriptive$shrNPP + plot_data_descriptive$treNPP
plot_data_descriptive$totalNPP <- plot_data_descriptive$herbNPP + plot_data_descriptive$woodyNPP
herb_npp_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "herbNPP",
y_lab = "Herb. NPP (g C/m2)")
herb_npp_plot
woody_npp_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "woodyNPP",
y_lab = "Woody NPP (g C/m2)")
woody_npp_plot
total_npp_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "totalNPP",
y_lab = "Total NPP (g C/m2)")
total_npp_plot
gh_npp_plots <- cowplot::plot_grid(herb_npp_plot + xlab(NULL) + scale_x_discrete(labels=NULL),
woody_npp_plot + xlab(NULL) + scale_x_discrete(labels=NULL),
total_npp_plot,
nrow = 3, rel_heights = c(0.6, 0.6, 1), align = "v",
labels = "auto")
gh_npp_plots
today <- Sys.Date()
highrestiff(plot_obj = gh_npp_plots,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_ghnpp_boxplot_", today, ".tif"),
width_in = 3,
height_in = 9,
resolution_dpi = 300)
highrestiff(plot_obj = total_npp_plot,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_totalnpp_boxplot_2_", today, ".tif"),
width_in = 3,
height_in = 3.5,
resolution_dpi = 300)
rm(npp)
# 4. biodiversity
diversity_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "shannon_diversity",
y_lab = "Shannon diversity index")
diversity_plot
richness_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "sp_richness",
y_lab = "Number of species")
richness_plot
biodiversity_plots <- cowplot::plot_grid(richness_plot, diversity_plot, total_npp_plot,
#nrow = 2, rel_heights = c(0.6, 1),
nrow = 1,
labels = "auto")
biodiversity_plots
today <- Sys.Date()
highrestiff(plot_obj = biodiversity_plots,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_biodiversity_boxplot_", today, ".tif"),
#width_in = 3,
width_in = 9,
#height_in = 6.5,
height_in = 3.5,
resolution_dpi = 300)
# 5. remote sensing variance
# Ideally, we should probably do this from the LandCART total foliar cover predictions, but obtaining those rasters for the whole UCRB is very slow. For now, we'll use the sum of the RAP 3.0 plant functional type cover values for each plot to represent total foliar cover instead. The RAP values were downloaded from within a 35-m radius of each point in this ESG using a Google Earth Engine script. Consider switching to LandCART, since it actually predicts a total foliar cover value, when raster generation and downloading for the UCRB 1995-2019 (April 1 - Nov 1) is complete.
library(broom) # for easy handling of many linear model parts
library(purrr) # for functional programming, eep!
library(tidyr)
VegCover_files <- list.files(paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/RAP_do_not_sync/", target_ESG, "/VegCover"),
full.names = T)
VegCover <- read.csv(VegCover_files[1])
VegCover$Year <- as.numeric(substr(VegCover_files[1], 168, 171))
for(file in VegCover_files[-1]){
VegCover_temp <- read.csv(file)
VegCover_temp$Year <- as.numeric(substr(file, 168, 171))
VegCover <- bind_rows(VegCover, VegCover_temp)
rm(VegCover_temp)
}
VegCover <- select(VegCover, PlotCod, Bst_grp, Bst_gr_, Year, AFG, BGR, LTR, PFG, SHR, TRE) %>%
setNames(c("PlotCode", "PlantCommunity_fuzzy", "PlantCommunity_name", "Year", "AFG", "BGR", "LTR", "PFG", "SHR", "TRE")) %>%
rowwise() %>%
mutate(TotalFoliar = sum(AFG, PFG, SHR, TRE, na.rm = T))
sampleyears <- select(plot_data_descriptive, PlotCode, Year) %>%
setNames(c("PlotCode", "Year_sampled"))
VegCover <- left_join(VegCover, sampleyears)
VegCover_10yrs <- VegCover %>% # select the 10 years up to and including the field sampling year
mutate(Years_to_sampling = Year_sampled - Year) %>%
filter(Years_to_sampling >=0 & Years_to_sampling <10)
VegCover_variance <- VegCover_10yrs %>%
select(PlotCode, PlantCommunity_fuzzy, Year, Years_to_sampling, TotalFoliar) %>%
as_tibble() %>%
nest(data = c(Year, Years_to_sampling, TotalFoliar)) %>% # create an internal data frame for each plot
mutate(model = map(data, ~lm(TotalFoliar ~ Years_to_sampling, data = .))) %>% # create a linear model from each plot's data frame
mutate(data = map(model, augment), # get residuals for each model and add to internal plot data frames
model_glance = map(model, glance)) %>% # get model summary statistics
mutate(resid_variance = map_dbl(data, ~var(.$.resid)), # get the variance of the detrended data (aka residuals)
resid_AC1 = map_dbl(data,
~acf(x = .$.resid, type = "correlation", lag.max = 1, plot = F)$acf[2]) # get the lag-1 autocorrelation of the detrended data
) %>%
select(-data, -model) %>%
unnest(model_glance)
# make figures
plot_data_descriptive <- left_join(plot_data_descriptive, select(VegCover_variance, PlotCode, resid_variance, resid_AC1))
variance_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "resid_variance",
y_lab = "Variance")
variance_plot
AC1_plot <- box_pretty_plot(dat = plot_data_descriptive,
y_var = "resid_AC1",
y_lab = "Lag-1 autocorrelation")
AC1_plot
RSvariability_plots <- cowplot::plot_grid(variance_plot + xlab(NULL) + scale_x_discrete(labels=NULL),
AC1_plot,
nrow = 2, rel_heights = c(0.6, 1),
labels = "auto")
RSvariability_plots
today <- Sys.Date()
highrestiff(plot_obj = RSvariability_plots,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_RSvariability_boxplot_", today, ".tif"),
width_in = 3,
height_in = 6,
resolution_dpi = 300)
variance_violin <- violin_pretty_plot(dat = plot_data_descriptive,
y_var = "resid_variance",
y_lab = "Variance",
fill_var = "PlantCommunity_fuzzy_reordered")
variance_violin
today <- Sys.Date()
highrestiff(plot_obj = variance_violin,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Violin/SW_ShallowDeepRocky_variance_violin_", today, ".tif"),
width_in = 3,
height_in = 3.5,
resolution_dpi = 300)
rm(VegCover, VegCover_10yrs, VegCover_variance)
# 6. Functional type cover grouped by functional type and colored by state
functypes_grps <- c("TotalFoliarCover", #"AH_WoodyCover", #AH_ForbGrassCover,
"AH_ShrubCover", "AH_TreeCover",
"AH_PerenForbGrassCover", "AH_AnnForbGrassCover"
)
functypes_cover_plot <- ggplot(data = select(plot_data_descriptive,
any_of(c("PlantCommunity_fuzzy_reordered", functypes_grps))) %>%
tidyr::pivot_longer(cols = all_of(functypes_grps),
names_to = "CoverType",
values_to = "PctCover"),
aes(fill=PlantCommunity_fuzzy_reordered, y=PctCover, x=CoverType)) +
geom_boxplot() +
scale_fill_manual(values = pal_veg2_reordered,
name = "Plant community"#,
#labels = group_labels_short
) +
scale_x_discrete(labels=c("Annual
herbaceous", "Perennial
herbaceous",
"Shrub", "Tree",
"Total foliar")) +
xlab("Plant functional type") +
ylab("Cover (%)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"),
axis.text.y = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
functypes_cover_plot
today <- Sys.Date()
highrestiff(plot_obj = functypes_cover_plot,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_pftcover_boxplot_", today, ".tif"),
width_in = 6,
height_in = 3.5,
resolution_dpi = 300)
# 7. Disturbance and land use
# MTBS
# Grazing?
# PJ treatments
# SWreGAP protection classes
# MTBS
library(sf)
library(lubridate)
mtbs <- st_read(dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/GIS/Disturbance_LandUse/mtbs_permis_DD_2020.shp")
fuzzy_cluster_memb_sf <- st_as_sf(fuzzy_cluster_memb_sf, crs = st_crs(plot_locations)) %>%
st_transform(., crs=st_crs(mtbs))
mtbs_plots <- st_join(x=fuzzy_cluster_memb_sf, y=mtbs) %>%
select(PlotCode, Best_group, Best_group_name, DateSampled, Ig_Date, geometry) %>%
mutate(Burned_presampling = ifelse(test = Ig_Date < DateSampled,
yes = TRUE,
no = FALSE)) %>%
group_by(PlotCode) %>%
summarise(Best_group = first(Best_group),
Best_group_name = first(Best_group_name),
DateSampled = first(DateSampled),
Times_burned_presampling = length(which(Burned_presampling)),
Times_burned_postsampling = length(which(!Burned_presampling)),
Date_burned_1 = sort(Ig_Date)[1],
Date_burned_2 = sort(Ig_Date)[2]) %>%
st_drop_geometry()
# write.csv(mtbs_plots,
# "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Analyses/SW_ShallowDeepRocky_burnedplots.csv",
# row.names = F)
rm(mtbs)
mtbs_plots$PlantCommunity_fuzzy_reordered <- factor(mtbs_plots$Best_group,
levels = c(as.character(1:8)),
labels = c(c("Juniper-pinyon woodland",
"Open juniper-pinyon shrubland",
"C4 grassland",
"Annualized juniper savanna",
"Pinyon-juniper woodland",
"Gambel oak woodland",
"Sagebrush shrubland/grassland",
"Bare juniper savanna")),
ordered = T)
mtbs_summary <- mtbs_plots %>%
group_by(PlantCommunity_fuzzy_reordered) %>%
summarise(pct_Burned = (length(which(Times_burned_presampling > 0))/n())*100,
pct_Burned_postsampling = (length(which(Times_burned_postsampling > 0))/n())*100)
burn_plot <- ggplot(data = mtbs_summary,
aes(x=PlantCommunity_fuzzy_reordered, y=pct_Burned, fill = PlantCommunity_fuzzy_reordered)) +
geom_col() +
scale_fill_manual(values = pal_veg2_reordered) +
xlab(NULL) +
ylab("% of plots burned pre-sampling") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"),
axis.text.y = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none")
burn_plot
burn_plot_post <- ggplot(data = mtbs_summary,
aes(x=PlantCommunity_fuzzy_reordered, y=pct_Burned_postsampling, fill = PlantCommunity_fuzzy_reordered)) +
geom_col() +
scale_fill_manual(values = pal_veg2_reordered) +
xlab(NULL) +
ylab("% of plots burned post-sampling") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"),
axis.text.y = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none")
burn_plot_post
today <- Sys.Date()
highrestiff(plot_obj = burn_plot,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_MTBSpctburnedpresample_boxplot_", today, ".tif"),
width_in = 3,
height_in = 3.5,
resolution_dpi = 300)
today <- Sys.Date()
highrestiff(plot_obj = burn_plot_post,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_MTBSpctburnedpostsample_boxplot_", today, ".tif"),
width_in = 3,
height_in = 3.5,
resolution_dpi = 300)
# PJ Treatments
WRI_polygons <- st_read(dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps/FickPJ",
layer = "PJ") %>%
st_transform(crs = st_crs(fuzzy_cluster_memb_sf))
sf_use_s2(FALSE)
pj_plots <- st_join(fuzzy_cluster_memb_sf, WRI_polygons, sparse = T) %>%
mutate(YearSampled = year(DateSampled),
TreatedPreSampling = ifelse(test = YearSampled>=yearStart,
yes = TRUE,
no = FALSE)) %>%
mutate(TreatedPreSampling = ifelse(test = is.na(TreatedPreSampling),
yes = FALSE,
no = TreatedPreSampling))
# Look at the types of treatments present, then categorize into broader types for visual summaries
sort(unique(paste(pj_plots$method, pj_plots$treat, sep = " - ")))
# additional documentation available in G:\Base Layers\Disturbance\TREATMENTS.
# For WRI treatments, see also the database at https://wri.utah.gov/wri/project/search.html
pj_plots[pj_plots$name %in% "Monument Ridge Bullhog", "method"] <- "bullhogfull size" # based on manual check of WRI website
pj_plots[pj_plots$name %in% "Monument Ridge Bullhog", "treat"] <- "brush removal"
pj_plots_cats <- pj_plots %>%
filter(!is.na(src)) %>%
mutate(GeneralTreatmenType = ifelse(test = treat %in% c(# my categories: prescribed fire, brush removal, herbaceous removal, seeding, soil improvements. CHECK aerator, mowing, and weeding for each ESG to make sure they still fit as assigned here!!!
"aerator",
"brush removal",
"mowing",
"weeding",
"herbicide"),
yes = "Woody removal",
no = ifelse(test = treat %in% c("seeding"),
yes = "Seeding only",
no = NA)),
TreatmentSubType = ifelse(test = GeneralTreatmenType=="Woody removal" & method %in% c("anchor chain",
"bullhogfull size",
"bullhogskid steer",
"double drum (1-way)",
"lop and scatter",
"mastication",
"mechanical",
"mowing",
"rollerchop",
"thinning"),
yes = "Mechanical",
no = ifelse(test = GeneralTreatmenType=="Woody removal" & method %in% c("herbicide",
"herbicide aerial"),
yes = "Chemical",
no = ifelse(test = GeneralTreatmenType=="Seeding only"
& method %in% c("broadcast"),
yes = "Broadcast seeding",
no = "Unknown"
))))
# consider just printing this as a table
count(st_drop_geometry(pj_plots_cats), GeneralTreatmenType, TreatmentSubType, seeded, burned)
n_plots_treated <- length(unique(filter(pj_plots_cats, TreatedPreSampling==T)$PlotCode))
n_plots_tx_twice_or_more <- pj_plots_cats %>%
filter(TreatedPreSampling==T) %>%
group_by(PlotCode) %>%
summarise(n_diff_treatments = n_distinct(yearStart)) %>%
filter(n_diff_treatments >= 2) %>%
nrow()
plots_per_group <- data.frame(Best_group = c(1:8),
plots_per_group = c(nrow(filter(plot_data_first, PlantCommunity_fuzzy==1)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==2)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==3)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==4)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==5)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==6)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==7)),
nrow(filter(plot_data_first, PlantCommunity_fuzzy==8))
))
pj_plot_pct <- pj_plots_cats %>%
left_join(., plots_per_group) %>%
filter(TreatedPreSampling==T) %>%
st_drop_geometry() %>%
group_by(PlotCode) %>% # next 3 lines select only the treatment closest to the sampling year
arrange(desc(yearStart), .by_group = T) %>%
mutate(Tx_flag = row_number()) %>%
filter(Tx_flag==1) %>%
ungroup() %>%
group_by(Best_group) %>%
summarise(MechanicalUnseeded = 100*length(which(TreatmentSubType=="Mechanical" & seeded==0))/first(plots_per_group),
MechanicalSeeded = 100*length(which(TreatmentSubType=="Mechanical" & seeded==1))/first(plots_per_group),
ChemicalUnseeded = 100*length(which(TreatmentSubType=="Chemical" & seeded==0))/first(plots_per_group),
ChemicalSeeded = 100*length(which(TreatmentSubType=="Chemical" & seeded==1)/first(plots_per_group))
) %>%
tidyr::pivot_longer(cols = c(MechanicalUnseeded, MechanicalSeeded, ChemicalUnseeded, ChemicalSeeded),
names_to = "TreatmentPair",
values_to = "pct") %>%
mutate(Best_group = as.factor(Best_group))
treatments_plot_pcts <- ggplot(data = pj_plot_pct, aes(x = Best_group, y = pct, fill = TreatmentPair)) +
geom_col() +
scale_fill_manual(values = c("#FDBF6F", "#FF7F00", "#CAB2D6", "#6A3D9A"),
name = "Treatment",
labels = c("Chemical woody treatment (seeded)",
"Chemical woody treatment",
"Mechanical woody treatment (seeded)",
"Mechanical woody treatment")) +
scale_x_discrete(name = NULL,
labels = c("Open woodland", "Shrubland", "Invaded", "Grassland")) +
ylab("% of plots treated in each
plant community") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"),
axis.text.y = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right")
treatments_plot_pcts
today <- Sys.Date()
highrestiff(plot_obj = treatments_plot_pcts,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_PJtreatments_pct_barplot_", today, ".tif"),
width_in = 6,
height_in = 3.5,
resolution_dpi = 300)
# GAP protection classes
GAP_raster <- raster::raster("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/GIS/Disturbance_LandUse/PADUS_GAP_Status_Code1.tif")
plot_GAP <- sf::st_as_sf(raster::extract(x = GAP_raster,
y = fuzzy_cluster_memb_sf,
sp = T)) %>%
sf::st_drop_geometry()
plot_data_descriptive <- left_join(plot_data_descriptive, select(plot_GAP, PlotCode, PADUS_GAP_Status_Code1))
GAP_proportions <- plot_data_descriptive %>%
group_by(PlantCommunity_fuzzy_reordered) %>%
# NOTE the PADUS codes are stored in the "GAP_Sts" attribute, but extraction pulls the raster value,
# which is different! Checked attribute table in Arc to sort it out!
summarise(pct_BiodivProt_NatDistPermitted = (length(which(PADUS_GAP_Status_Code1 == 4))/n())*100,
pct_BiodivProt_NatDistSupressed = (length(which(PADUS_GAP_Status_Code1 == 3))/n())*100,
pct_MultiUse_ExtractPermitted = (length(which(PADUS_GAP_Status_Code1 == 1))/n())*100,
pct_NoBiodivProctect = (length(which(PADUS_GAP_Status_Code1 == 2))/n())*100,
pct_PrivateLand = (length(which(is.na(PADUS_GAP_Status_Code1)))/n())*100) %>%
tidyr::pivot_longer(data=.,
cols = c(pct_BiodivProt_NatDistPermitted,
pct_BiodivProt_NatDistSupressed,
pct_MultiUse_ExtractPermitted,
pct_NoBiodivProctect,
pct_PrivateLand),
names_to = "ProtectionStatus",
values_to = "Percent")
protection_plot <- ggplot(data = GAP_proportions,
aes(x=ProtectionStatus, y=Percent, fill=PlantCommunity_fuzzy_reordered)) +
geom_col() +
scale_fill_manual(values = pal_veg2_reordered,
name = "Plant community"
) +
scale_x_discrete(labels=c("1 - Biodiversity protected,
natural disturbance permitted",
"2 - Biodiversity protected,
natural disturbance suppressed",
"3 - Multi-use, extraction
permitted",
"4 - No biodiversity protection",
"Private land")) +
xlab("PADUS protection status") +
ylab("% plots in
plant community") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"),
axis.text.y = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = unit(c(2,2,2,5), "mm")
)
protection_plot
today <- Sys.Date()
highrestiff(plot_obj = protection_plot,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_PADUSprotection_flipped_barplot_", today, ".tif"),
width_in = 6,
height_in = 3.5,
resolution_dpi = 300)
group_labels_fuzzy <- paste0(group_labels_base,
" (n = ",
c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1)),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2)),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3)),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4)),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==5)),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==6)),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==7)),
nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==8))
),
")")
protection_plot2 <- ggplot(data = GAP_proportions,
aes(fill=ProtectionStatus, y=Percent, x=PlantCommunity_fuzzy_reordered)) +
geom_col() +
scale_fill_manual(values = c("seagreen", "yellowgreen", "gray58", "dodgerblue", "palegoldenrod"),
name = "PADUS protection status",
labels = c("1 - Biodiversity protected,
natural disturbance permitted",
"2 - Biodiversity protected,
natural disturbance suppressed",
"3 - Multi-use, extraction
permitted",
"4 - No biodiversity protection",
"Private land")) +
scale_x_discrete(labels= group_labels_fuzzy) +
xlab(NULL) +
ylab("% plots in plant community") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"),
axis.text.y = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.key.size = unit(1.5, 'lines'))
protection_plot2
today <- Sys.Date()
highrestiff(plot_obj = protection_plot2,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_PADUSprotection_dodge_barplot_", today, ".tif"),
width_in = 6,
height_in = 3.5,
resolution_dpi = 300)
GAP_proportions_byclass <- plot_data_descriptive %>%
mutate(PADUS_GAP_Status_lumped = ifelse(test = PADUS_GAP_Status_Code1 %in% c(3,4),
yes = "Protected public",
no = NA)) %>%
mutate(PADUS_GAP_Status_lumped = ifelse(test = PADUS_GAP_Status_Code1 %in% c(1,2),
yes = "Unprotected public",
no = PADUS_GAP_Status_lumped)) %>%
mutate(PADUS_GAP_Status_lumped = ifelse(test = is.na(PADUS_GAP_Status_Code1),
yes = "Private",
no = PADUS_GAP_Status_lumped)) %>%
group_by(PADUS_GAP_Status_lumped) %>%
summarise(pct_JPwoodland = (length(which(PlantCommunity_fuzzy_reordered == "Juniper-pinyon woodland"))/n())*100,
pct_OpenJPshrubland = (length(which(PlantCommunity_fuzzy_reordered == "Open juniper-pinyon shrubland"))/n())*100,
pct_C4grassland = (length(which(PlantCommunity_fuzzy_reordered == "C4 grassland"))/n())*100,
pct_AnnualizedJuniperSavanna = (length(which(PlantCommunity_fuzzy_reordered == "Annualized juniper savanna"))/n())*100,
pct_PJwoodland = (length(which(PlantCommunity_fuzzy_reordered == "Pinyon-juniper woodland"))/n())*100,
pct_GambelOakWoodland = (length(which(PlantCommunity_fuzzy_reordered == "Gambel oak woodland"))/n())*100,
pct_SagebrushShrubland = (length(which(PlantCommunity_fuzzy_reordered == "Sagebrush shrubland/grassland"))/n())*100,
pct_BareJuniperSavanna = (length(which(PlantCommunity_fuzzy_reordered == "Bare juniper savanna"))/n())*100,
) %>%
tidyr::pivot_longer(data=.,
cols = c(pct_JPwoodland,
pct_OpenJPshrubland,
pct_C4grassland,
pct_AnnualizedJuniperSavanna,
pct_PJwoodland,
pct_GambelOakWoodland,
pct_SagebrushShrubland,
pct_BareJuniperSavanna),
names_to = "PlantCommunity",
values_to = "Percent") %>%
mutate(PlantCommunity_fuzzy_reordered = factor(PlantCommunity,
levels = c("pct_JPwoodland",
"pct_OpenJPshrubland",
"pct_C4grassland",
"pct_AnnualizedJuniperSavanna",
"pct_PJwoodland",
"pct_GambelOakWoodland",
"pct_SagebrushShrubland",
"pct_BareJuniperSavanna"),
labels = c("Juniper-pinyon woodland",
"Open juniper-pinyon shrubland",
"C4 grassland",
"Annualized juniper savanna",
"Pinyon-juniper woodland",
"Gambel oak woodland",
"Sagebrush shrubland/grassland",
"Bare juniper savanna"),
ordered = T))
protection_lumped_plot <- ggplot(data = GAP_proportions_byclass,
aes(x=PADUS_GAP_Status_lumped, y=Percent, fill=PlantCommunity_fuzzy_reordered)) +
geom_col() +
scale_fill_manual(values = pal_veg2_reordered,
name = "Plant community"
) +
#scale_x_discrete(labels=c("Protected - public", "Unprotected - public" "Private land")) +
xlab(NULL) +
ylab("% of plots in
protection group") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"),
axis.text.y = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()#,
#plot.margin = unit(c(2,2,2,5), "mm")
)
protection_lumped_plot
disturbance_plots <- cowplot::plot_grid(burn_plot, protection_lumped_plot,
nrow = 1, rel_widths = c(0.7, 1),
labels = "auto")
disturbance_plots
today <- Sys.Date()
highrestiff(plot_obj = disturbance_plots,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_disturbanceproctection_barplot_", today, ".tif"),
width_in = 7,
height_in = 3.5,
resolution_dpi = 300)
rm(GAP_raster, plot_GAP)
# Make materials for spring 2023 workshop handouts
# 1. condensed community tables
# 2. bar charts of functional types and ground cover - separate copy for each state
### community tables
functype_table_palette <- c("#80B1D3", "#80B1D3", "#9ECAE1",
"#FFFFB3", "#FFFFB3", "#FFFFB3",
"#FFFFB3", "#C6DBEF", "#FB8072",
"#BEBADA", "#8DD3C7")
names(functype_table_palette) <- sort(unique(dom_spp_grp1$Indicator))
dom_spp_grp2_workshop <- dom_spp_grp2_pub %>%
filter(`Mean cover (%)` >= 1) %>%
arrange(desc(`Mean cover (%)`))
dom_spp_grp1_workshop <- dom_spp_grp1_pub %>%
filter(`Mean cover (%)` >= 1) %>%
arrange(desc(`Mean cover (%)`))
dom_spp_grp4_workshop <- dom_spp_grp4_pub %>%
filter(`Mean cover (%)` >= 1) %>%
arrange(desc(`Mean cover (%)`))
dom_spp_grp3_workshop <- dom_spp_grp3_pub %>%
filter(`Mean cover (%)` >= 1) %>%
arrange(desc(`Mean cover (%)`))
dom_spp_grp5_workshop <- dom_spp_grp5_pub %>%
filter(`Mean cover (%)` >= 1) %>%
arrange(desc(`Mean cover (%)`))
dom_spp_grp6_workshop <- dom_spp_grp6_pub %>%
filter(`Mean cover (%)` >= 1) %>%
arrange(desc(`Mean cover (%)`))
dom_spp_grp7_workshop <- dom_spp_grp7_pub %>%
filter(`Mean cover (%)` >= 1) %>%
arrange(desc(`Mean cover (%)`))
dom_spp_grp8_workshop <- dom_spp_grp8_pub %>%
filter(`Mean cover (%)` >= 1) %>%
arrange(desc(`Mean cover (%)`))
workingdat <- dom_spp_grp6_workshop
kableExtra::kbl(x = workingdat,
format = "html", row.names = F,
col.names = c("Functional group",
"Scientific name",
"Common name",
"USDA Plants
code",
"Mean
cover (%)",
"% of plots
where present"), align = "c",
caption = #"Open juniper-pinyon woodland" #grp1
#"Mixed shrubs with open juniper-pinyon" #grp2
#"C4 grassland with mixed shrubs and scattered juniper-pinyon" #grp3
#"Annualized invaded mixed shrubland/C3 grassland with scattered PJ" #grp4
#"Open pinyon-juniper woodland" #grp5
"Gambel oak woodland with mixed shrubs and grasses" #grp6
#"Sagebrush with mixed C3 grasses" #grp7
#"Juniper-pinyon with light C3 grass understory" #grp8
) %>%
kable_styling(bootstrap_options = c("bordered", "condensed")) %>%
column_spec(kable_input=., column=4, width = "2cm") %>%
column_spec(kable_input=., column=5, width = "2cm") %>%
column_spec(kable_input=., column=6, width = "2.5cm") %>%
row_spec(kable_input=., row=0, bold=T, background = "white") %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="C3 introduced perennial grasses"), background = functype_table_palette["C3 introduced perennial grasses"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="C3 native perennial grasses"), background = functype_table_palette["C3 native perennial grasses"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="C4 native perennial grasses"), background = functype_table_palette["C4 native perennial grasses"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="Introduced annual forbs"), background = functype_table_palette["Introduced annual forbs"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="Introduced annual grasses"), background = functype_table_palette["Introduced annual grasses"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="Introduced perennial forbs"), background = functype_table_palette["Introduced perennial forbs"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="Native annual forbs"), background = functype_table_palette["Native annual forbs"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="Native annual grasses"), background = functype_table_palette["Native annual grasses"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="Native perennial forbs"), background = functype_table_palette["Native perennial forbs"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="Shrubs"), background = functype_table_palette["Shrubs"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="Succulents"), background = functype_table_palette["Succulents"]) %>%
row_spec(kable_input=., row = which(workingdat$`Functional group`=="Trees"), background = functype_table_palette["Trees"])
### bar charts
functypes_grps <- c("TotalFoliarCover", #"AH_WoodyCover", #AH_ForbGrassCover,
"AH_ShrubCover", "AH_TreeCover",
"AH_PerenForbGrassCover", "AH_AnnForbGrassCover"
)
plot_data_descriptive_grp <- filter(plot_data_descriptive, PlantCommunity_fuzzy==8)
functypes_cover_states <- ggplot(data = select(plot_data_descriptive_grp,
any_of(c("PlantCommunity_fuzzy", functypes_grps))) %>%
tidyr::pivot_longer(cols = all_of(functypes_grps),
names_to = "CoverType",
values_to = "PctCover"),
aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) +
geom_boxplot(outlier.shape = NA) + # make outliers invisible
scale_fill_manual(values = c("#FFFFB3", "#80B1D3", "#FB8072", "#8DD3C7", "#666666"),
name = "Functional type",
labels = c("Annual herbaceous", "Perennial herbaceous",
#"Woody",
"Shrub", "Tree",
"Total foliar cover")
) +
scale_x_discrete(labels=group_labels_fuzzy[8]) +
xlab(NULL) +
ylab("Cover (%)") +
ylim(c(0,100)) +
ggtitle("Plant functional groups") +
theme_bw() +
theme(axis.text.x = element_text( color = "black"),
axis.text.y = element_text(color = "black"))
functypes_cover_states
today <- Sys.Date()
highrestiff(plot_obj = functypes_cover_states,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/JFSP_Workshop_1/PlantFuncTypes_", "8", "_", today, ".tif"),
width_in = 6,
height_in = 3.5,
resolution_dpi = 300)
ground_grps <- c("BareSoilCover", "FH_TotalLitterCover", "CA_percent_200plus",
"FH_LichenMossCover", "CA_percent_100plus"
)
plot_data_first_tall_grp <- filter(plot_data_first_tall, PlantCommunity_fuzzy==8)
ground_cover_states <- ggplot(data = filter(plot_data_first_tall_grp, CoverType %in% ground_grps),
aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) +
geom_boxplot(outlier.shape = NA) +
scale_fill_manual(values = c("#8DA0CB", "#E78AC3", "#A6D854", "#FFD92F", "#E5C494"),
name = "Functional group",
labels = names(cover_labels[which(cover_labels %in% ground_grps)])
) +
scale_x_discrete(labels=group_labels_fuzzy[8]) +
xlab(NULL) +
ylab("Cover (%)") +
ylim(c(0,100)) +
ggtitle("Ground cover") +
theme_bw() +
theme(axis.text.x = element_text(color = "black"),
axis.text.y = element_text(color = "black"))
ground_cover_states
today <- Sys.Date()
highrestiff(plot_obj = ground_cover_states,
file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/JFSP_Workshop_1/GroundCover_", "8", "_", today, ".tif"),
width_in = 6,
height_in = 3.5,
resolution_dpi = 300)
annack84/STMdevelopment documentation built on April 12, 2024, 6:46 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
NOTE: This correspondence is intended for communication of project progress among funders and collaborators only. This information is preliminary and is subject to revision. It is being provided to meet the need for timely best science. The information is provided on the condition that neither the U.S. Geological Survey nor the U.S. Government shall be held liable for any damages resulting from the authorized or unauthorized use of the information.
General information
knitr::opts_chunk$set(echo = TRUE) target_ESG <- "Semiarid_Warm_Shallow_DeepRocky" user <- "Duniway_et_al_2024" # use this user to get the exact data used for the Feb 2023 JFSP workshops #user <- "Anna" # use one of these two options to get the most current data #user <- "VPN" output_figure_folder <- file.path("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures", target_ESG) library(dplyr) library(kableExtra) library(sf) library(STMdevelopment) # library(devtools);load_all()
# TODO insert dot graphic and text about the status of the ESD (e.g. "Provisional"). Might help # to create a function to automatically insert the description based on the # status input
EDIT_map(target_ESG = target_ESG, user = user, maxZoomStyle = "public")
Areas shown in blue indicate the maximum mapped extent of this ecological site. Other ecological sites likely occur within the highlighted areas. It is also possible for this ecological site to occur outside of highlighted areas if detailed soil survey has not been completed or recently updated.
Associated ecological site concepts
This ecological site group includes a variety of upland ecological sites that have shallow to deep soils with high rock content. The group was defined by placing ecological sites into groups by separating soils, climate, and geomorphic features to best differentiate reference vegetation production values and documented ecological states as documented by Nauman et al. (2022).
esds <- esd_data_pull(user=user, target_ESG = target_ESG) comps <- comp_data_pull(user=user, target_ESG = target_ESG) num_esds <- nrow(esds)
This ecological site group includes r num_esds ecological sites. These sites were correlated to the group based on average soil, climate and geomorphic properties of the soil components linked to each ecological site.
linked_esd_names <- esds[,c("ecoclassid","ecoclassname")] %>% mutate(ecoclassid = paste0("[", ecoclassid, "](https://edit.jornada.nmsu.edu/catalogs/esd/", substr(ecoclassid, 2, 5), "/", ecoclassid, ")")) ## Associated ecological sites format_tables_EDIT_style(data = linked_esd_names, col.names=c("ESD Code","ESD Name"), caption = "Ecological Sites associated with this ESG", row.names = F) #knitr::kable(x = esds[,c("ecoclassid","ecoclassname")], col.names=c("ESD Code","ESD Name"), caption = NULL,row.names = F)
# TODO create table of dominant plant functional groups and species esg_prod_list <- esg_production_pull(user=user, target_ESG = target_ESG) dom_plants_df <- dplyr::tribble(~Item, ~Description, "Trees", as.character(esg_prod_list[["Tree.df"]]$COMMON_NAME)[1:5], "Shrubs", as.character(esg_prod_list[["Shrub.df"]]$COMMON_NAME)[1:5], "Grasses",as.character(esg_prod_list[["Grass.df"]]$COMMON_NAME)[1:5], "Forbs", as.character(esg_prod_list[["Forb.df"]]$COMMON_NAME)[1:5]) format_tables_EDIT_style(data = dom_plants_df, caption = "Dominant plant species") #knitr::kable(x = dom_plants_df, caption = "Table 1. Dominant plant species")
TODO insert "Download full description" link to a PDF
Physiographic features
# TODO create table of landforms, elevations, slopes, water table depth, flooding duration, flooding frequency, ponding frequency # ideally should be able to tab between US and metric system measurements landforms <- unlist(strsplit(esds$landfrms[1], split=",")) for(e in 2:nrow(esds)){ esdlfs <- unlist(strsplit(esds[e,c("landfrms")],split=",")) landforms <- append(landforms,esdlfs) } landforms <- landforms[!landforms %in% "NA"] lfsum <- factor(landforms) lfsum <- reorder(lfsum,lfsum,FUN=function(x) -length(x)) landforms_unique <- unique(landforms)
These sites occur on upland positions on a variety of landforms (n=r length(landforms_unique)). Moderate to steep slopes and rugged topography are common. Less commonly, some sites may occur on high-energy alluvial deposits.
# Slope slope_ave <- mean(comps@site$slope_r.x, na.rm = T) slope_75th <- unname(quantile(comps@site$slope_r.x, probs = 0.75, na.rm = T)) slope_max <- max(comps@site$slope_r.x, na.rm = T) slope_sd <- sd(comps@site$slope_r.x, na.rm = T) slope_field <- paste("Slopes average",format(round(slope_ave, 1), nsmall = 1), "% with a standard deviation of",format(round(slope_sd, 1), nsmall = 1), "and are generally less than", format(round(slope_75th, 1), nsmall = 1), "but can be up to", format(round(slope_max, 1), nsmall = 1),"%",sep=" ") phys.df <- dplyr::tribble(~Item, ~Description, "Landforms (Top 10)", as.character(unique(lfsum))[1:10], "Slope", slope_field, "Flooding Frequency","None", "Aspect", "Aspect is not a significant factor") format_tables_EDIT_style(data = phys.df, caption = "Representative physiographic features") #knitr::kable(x = phys.df, caption = "Table 2. Representative physiographic features")
Climatic features
Climate is generally semiarid and warm with aridity index (precipitation / potential evapotranspiration) values averaging r format(round(mean(esds$aimean), 2), nsmall = 2) and ranging from r format(round(min(esds$aimean), 2), nsmall = 2) to r format(round(max(esds$aimean), 2), nsmall = 2). The average maximum temperatures (Celsius) of the hottest month range from r format(round(min(esds$maxtempmean), 2), nsmall = 2) to r format(round(max(esds$maxtempmean), 2), nsmall = 2), and minimum temperatures of the coldest month range from r format(round(min(esds$mintempmean), 2), nsmall = 2) to r format(round(max(esds$mintempmean), 2), nsmall = 2). Warm season (June to September) precipitation makes up r format(round(mean(esds$pptrt*100), 0), nsmall = 0)% of total precipitation on average, but can range from r format(round(min(esds$pptrt*100), 0), nsmall = 0)% to r format(round(max(esds$pptrt*100), 0), nsmall = 0)%.
clim.df <- dplyr::tribble(~Item, ~Description, "Frost Free Period (days)", paste(format(round(min(comps@site$ffd_r), 0), nsmall = 0),format(round(mean(comps@site$ffd_r), 0), nsmall = 0),format(round(max(comps@site$ffd_r), 0), nsmall = 0),sep=", "), "Mean Annual Precipitation (in)", paste(format(round(min(comps@site$reannualprecip_r/25.4,na.rm=T), 1), nsmall = 1),format(round(mean(comps@site$reannualprecip_r/25.4,na.rm=T), 1), nsmall = 1),format(round(max(comps@site$reannualprecip_r/25.4,na.rm=T), 1), nsmall = 1),sep=", ")) format_tables_EDIT_style(data = clim.df, caption = "Representative climatic features (min, mean, max)") #knitr::kable(x = clim.df, caption = "Table 3. Representative climatic features (min, mean, max)")
# TODO create line chart and bar charts of monthly high/low precip that you can tab between
# TODO create line chart and bar charts of monthly high/low temps that you can tab between
Water features
These sites neither benefit significantly from run-in moisture nor experience excessive loss of moisture from runoff.
Soil features
Soils in this group have shallow to deep depth to bedrock and typically have > 30% rock fragments by volume. Surface and subsurface horizons have a variety of textures ranging from sand to clay loam. Soils are non-saline and non sodic and often have high calcium carbonate contents. Soil pH ranges from r min(c(round(min(comps@horizons[comps@horizons$hzdept_r<30,c("ph1to1h2o_r")],na.rm=T), 0), round(min(comps@horizons[comps@horizons$hzdept_r>30,c("ph1to1h2o_r")],na.rm=T), 0))) to r max(c(round(max(comps@horizons[comps@horizons$hzdept_r<30,c("ph1to1h2o_r")],na.rm=T), 1), round(max(comps@horizons[comps@horizons$hzdept_r>30,c("ph1to1h2o_r")],na.rm=T), 1))). Erosion hazard is generally low but bare ground states are possible.
## Prep soil variables for table surftext <- factor(esds$txtnm_surf) surftext <- reorder(surftext,surftext,FUN=function(x) -length(x)) surftext_uni <- levels(surftext) subtext <- factor(esds$txtnm_sub) subtext <- reorder(subtext,subtext,FUN=function(x) -length(x)) subtext_uni <- levels(subtext) drainage <- factor(comps@site$drainagecl.x) drainage <- reorder(drainage,drainage,FUN=function(x) -length(x)) drainage_uni <- levels(drainage)[1:3] ## Make Table of soil descriptors: # TODO Could update pH to be based on esd averages instead of comps. Also could update PM to grab ssurgo values soil.df <- dplyr::tribble(~Item, ~Description, "Parent Material", "alluvium and eolian sediments", "Surface Texture (0-30cm)", surftext_uni, "Subsurface Texture (>30cm)",subtext_uni, "Drainage", drainage_uni, "Soil Depth", paste(format(round(min(esds$depthmean,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$depthmean,na.rm=T), 0), nsmall = 0), " cm",sep=""), "Surface Rock Content %vol (0-30cm)", paste(format(round(min(esds$rock_surf,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$rock_surf,na.rm=T), 0), nsmall = 0), "%",sep=""), "Subsurface Rock Content %vol (>30cm)", paste(format(round(min(esds$rock_sub,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$rock_sub,na.rm=T), 0), nsmall = 0), "%",sep=""), "Surface Electrical Conductivity (0-30cm)", paste(format(round(min(esds$ec_surf,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$ec_surf,na.rm=T), 1), nsmall = 1), " dS/m",sep=""), "Subsurface Electrical Conductivity (>30cm)", paste(format(round(min(esds$ec_sub,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$ec_sub,na.rm=T), 1), nsmall = 1), " dS/m",sep=""), "Surface Sodium Adsorption Ratio (0-30cm)", paste(format(round(min(esds$sar_surf,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$sar_surf,na.rm=T), 1), nsmall = 1),sep=""), "Subsurface Sodium Adsorption Ratio (>30cm)", paste(format(round(min(esds$sar_sub,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$sar_sub,na.rm=T), 1), nsmall = 1),sep=""), "Surface 1:1 pH (0-30cm)", paste(format(round(min(comps@horizons[comps@horizons$hzdept_r<30,c("ph1to1h2o_r")],na.rm=T), 0), nsmall = 0),"-",format(round(max(comps@horizons[comps@horizons$hzdept_r<30,c("ph1to1h2o_r")],na.rm=T), 1), nsmall = 1),sep=""), "Subsurface 1:1 pH (>30cm)", paste(format(round(min(comps@horizons[comps@horizons$hzdept_r>30,c("ph1to1h2o_r")],na.rm=T), 0), nsmall = 0),"-",format(round(max(comps@horizons[comps@horizons$hzdept_r>30,c("ph1to1h2o_r")],na.rm=T), 1), nsmall = 1),sep=""), "Surface Calcium Carbonate (0-30cm)", paste(format(round(min(esds$caco3_surf,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$caco3_surf,na.rm=T), 0), nsmall = 0), "%",sep=""), "Subsurface Calcium Carbonate (>30cm)", paste(format(round(min(esds$caco3_sub,na.rm=T), 0), nsmall = 0),"-",format(round(max(esds$caco3_sub,na.rm=T), 0), nsmall = 0), "%",sep="")) format_tables_EDIT_style(data = soil.df, caption = "Representative soil features") #knitr::kable(x = soil.df, caption = "Table 4. Representative soil features")
Soil Component data
Soil types correlated to this group include r length(unique(comps@site$cokey)) different components mapped in SSURGO. Soil taxonomic units include:
cat(sort(unique(comps@site$compname)), sep = ", ")
Component level soil property depth plots
agg <- aqp::slab(comps, fm= ~ claytotal_r + silttotal_r + sandtotal_r + ec_r + caco3_r) lattice::xyplot(top ~ p.q50 | variable, data=agg, ylab='Depth (cm)', xlab='median bounded by 25th and 75th percentiles', lower=agg$p.q25, upper=agg$p.q75, ylim=c(201,-2), panel=aqp::panel.depth_function, alpha=0.25, sync.colors=TRUE, par.settings=list(superpose.line=list(col='RoyalBlue', lwd=2)), prepanel=aqp::prepanel.depth_function, cf=agg$contributing_fraction, cf.col='black', cf.interval=20, layout=c(5,1), strip=lattice::strip.custom(bg=grey(0.8)), scales=list(x=list(tick.number=4, alternating=3, relation='free')) )
Ecological dynamics
These site are susceptible to drought, especially in shallower soils. Erosion issues are possible but rare. Compared to reference states, alternative states may have woody encroachment, high bare ground, and/or loss of biological soil crust and perennial grasses.
State and transition model development
# TODO update with the newly tabulated ESD STMs from Travis, then set to eval = TRUE apriori_stms <- read.csv(data_file_paths(user)$apriori_stms, stringsAsFactors = F, na.strings = c("NA", "", " ")) apriori_stms_target <- select(apriori_stms, any_of(c("State", "Plant_community", target_ESG))) n_ESDs <- apriori_stms_target[apriori_stms_target$Plant_community=="n_ESDs", target_ESG] colnames(apriori_stms_target) <- c("State", "Plant community", "Proportion of ESDs containing this plant community") format_tables_EDIT_style(data = apriori_stms_target[-nrow(apriori_stms_target),], caption = paste("State and plant community summary from Ecological Sites within this ESG.", n_ESDs, "complete STMs were available to include.")) %>% collapse_rows(columns = 1, valign = "top")
Reference state production indicators
## Plot of reference production indicators used in paper resp_vars <- colnames(esds)[1:21] resp_vars <- resp_vars[!resp_vars %in% c("ecoclassid")] par(mar = c(4,8,4,2)) ## This works when knitted boxplot(esds[,resp_vars], horizontal = TRUE,las = 1, main = "Reference Production (lbs/Acre)")
# by functional group cover indicators <- c("AH_C3NativePerenGrassCover", "AH_C3IntroducedPerenGrassCover", "AH_C4NativePerenGrassCover", "AH_C4IntroducedPerenGrassCover", "AH_NativePerenForbCover", # Native perennial forbs "AH_IntroducedPerenForbCover", # Non-native perennial forbs "AH_NativeAnnGrassCover", # Native annual grasses "AH_IntroducedAnnGrassCover", # Non-native annual grasses "AH_NativeAnnForbCover", # Native annual forbs "AH_IntroducedAnnForbCover", # Non-native annual forbs "AH_ArtemisiaTridentataCover", # All big sagebrush species lumped together "BareSoilCover", # Bare soil "FH_TotalLitterCover", # Litter #"CA_percent_100plus", # Canopy gaps > 100 cm #"CA_percent_200plus", "FH_LichenCover", # Lichen + moss combined cover "FH_MossCover"#, #"SoilStab_all" # to represent cyano abundance... I think this is only marginally a community variable. If used in hierarchical clustering and Bray-Curtis NMDS, will need to standardize it to make units more similar to % cover units ) ann_grass_by_spp <- F ann_forb_by_spp <- F per_grass_by_spp <- F per_forb_by_spp <- F succulent_by_spp <- F shrub_by_spp <- T # All shrubs and sub-shrubs by species subshrub_by_spp <- T tree_by_spp <- T # All trees by species opuntia_combined <- T data_sources <- c(#"BadgerWash", #"CRCMonitoring", # drop if not needed for spatial representation "IM_NCPN", "LMF", "NRI", #"Parashant", # drop if not needed for spatial representation "AIM"#, #"VanScoyocThesis" # drop if not needed for spatial representation ) impute_gap_type <- NULL plot_data <- plot_data_pull(user=user, target_ESG = target_ESG, data_sources = data_sources, indicators = indicators, ann_grass_by_spp = ann_grass_by_spp, ann_forb_by_spp = ann_forb_by_spp, per_grass_by_spp = per_grass_by_spp, per_forb_by_spp = per_forb_by_spp, succulent_by_spp = succulent_by_spp, shrub_by_spp = shrub_by_spp, subshrub_by_spp = subshrub_by_spp, tree_by_spp = tree_by_spp, opuntia_combined = opuntia_combined, impute_gap_type = impute_gap_type ) %>% filter(!is.na(Month))
# TODO update I&M species codes to match current USDA Plants codes! probably a step for the plot_data_pull function? # remove incomplete rows - ordination won't work with NAs plot_data_clean <- na.omit(plot_data) %>% # NRI has a separate plot code for each time a plot was sampled, but all other data sets have the same plot # code for all sampling times. Make a similar code for NRI. mutate(PlotCode_NoYear = ifelse(test = SourceKey == "NRI_UCRB", yes = paste0(SourceKey, "_", SiteName, "_", substr(PlotID, 5, nchar(PlotID))), no = PlotCode)) # Some plots were sampled multiple times - I'll include only the first year that a plot has complete data plot_data_first <- plot_data_clean %>% #plot_data_clean %>% dplyr::group_by(PlotCode_NoYear) %>% dplyr::arrange(.data=., Year, .by_group = TRUE) %>% dplyr::filter(dplyr::row_number()==1) %>% dplyr::ungroup() # remove plots with SGU probability (i.e. certainty of prediction) in the 10th percentile SGU_prob_raster <- raster::raster("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps/UCRB_SGUs_ProbMax/UCRB_SGUs_ProbMax.tif") file_paths <- data_file_paths(user) plot_locations <- sf::st_read(dsn = file.path(file_paths$plotnet_processed, "PlotLocations"), layer = "all_plot-years_2022-11-08", quiet=TRUE) %>% distinct() if("NRI" %in% data_sources){ nri_locations <- sf::st_read(dsn = file_paths$nri, layer = "NRI_UCRB_plot-years_2022-11-03", quiet = TRUE) plot_locations <- dplyr::filter(plot_locations, !grepl(pattern = "^NRI_", x=PlotCode)) %>% dplyr::bind_rows(., nri_locations) %>% distinct() } plot_SGU_prob <- sf::st_as_sf(raster::extract(x = SGU_prob_raster, y = plot_locations, sp = T)) plot_data_first <- left_join(plot_data_first, select(plot_SGU_prob, PlotCode, UCRB_SGUs_ProbMax)) %>% distinct() certainty_cutoff_level <- 0.1 certainty_cutoff <- quantile(plot_data_first$UCRB_SGUs_ProbMax, probs = certainty_cutoff_level, na.rm = T) plot_data_first <- filter(plot_data_first, UCRB_SGUs_ProbMax > certainty_cutoff) %>% filter(PlotCode != "AIM_Utah Vernal FO 2019_009") # this one is duplicated in the plot locations with very different coordinates # make the clustering and ordination data frame - can't include ANY columns except the actual variables! ord.df <- dplyr::select(plot_data_first, -SourceKey, -PlotID, -SiteName, -PlotName, -Year, -Longitude_NAD83, -Latitude_NAD83, -Month, -Day, -UCRB_SGUs_ProbMax, -geometry, -PlotCode_NoYear ) %>% tibble::column_to_rownames("PlotCode") # keep an ID code for the plot as the row name to prevent data scrambling problems # If using soil stability, standardize variables to make units more similar # Here, dividing by the max possible value for each measurement type (cover values # divided by 100, soil stability divided by 6) if("SoilStab_all" %in% indicators){ ord.df.cover <- select(ord.df, -SoilStab_all)/100 ord.df.soilstab <- select(ord.df, SoilStab_all)/6 ord.df <- cbind(ord.df.cover, ord.df.soilstab) } # remove rare species sp_pa <- mutate(rowwise(ord.df), across(everything(), function(x){if(x>0){1}else{0}})) # make a presence/absence data frame sp_keep <- names(colSums(sp_pa)[which(colSums(sp_pa)>2)]) ord.df <- select(ord.df, all_of(sp_keep)) # remove unknown species if present ord.df <- select(ord.df, -any_of(c("UNKS", "UNKPG", "UNKAF", "UNKPF", "PG1", "PG01", "AF1"))) # # save clustering data for other analyses - exclude location data to keep kosher # # with the data sharing restrictions # ord.df.share <- select(plot_data_first, all_of(c("SourceKey", "PlotCode", "Year", "Month", "Day", "UCRB_SGUs_ProbMax", colnames(ord.df)))) # write.csv(ord.df.share, # file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Data/Analysis_ready_data/", # target_ESG, "_FANNYInputData_", Sys.Date(), ".csv"), # row.names = F)
Field data from were collected by field crews with the Natural Resource Conservation Service's (NRCS) National Resource Inventory (NRI) program, the Bureau of Land Management's (BLM) Assessment, Inventory, and Monitoring (AIM) and Landscape Monitoring Framework (LMF) programs, and the Nationa Park Service's (NPS) Inventory and Monitoring Northern Colorado Plateau Network (NCPN) program between r min(plot_data_clean$Year) and r max(plot_data_clean$Year). Plots within the target ESG with SGU classification probability less than or equal to the r certainty_cutoff_level*100th percentile (r certainty_cutoff%) were excluded, resulting in a total of r nrow(ord.df) plots for this ESG. Plant and soil cover were measured using line-point intercept and canopy gap methods (Herrick et al., 2017). Plot-level indicators (Table \@ref(tab:doc-plot-indicators)) were calculated using the terradactyl package in R (McCord et al., 2022). When plots were sampled in multiple years, only the first sampling year was used in STM development. r if(!is.null(impute_gap_type)){paste("For plots that were missing annual & perennial gap data but had perennial-only gap data (primarily I&M-NCPN), annual & perennial gap values were imputed from perennial-only gaps and annual herbaceous cover using a linear model built from LMF and NRI plots that had both types of gap (r-squared of 0.83 and 0.81 for 100-cm and 200-cm gap predictions, respectively).")}
rm(plot_data_clean, plot_data) indicator_descriptions <- make_indicator_descriptions(indicators = indicators, shrub_by_spp = shrub_by_spp, subshrub_by_spp = subshrub_by_spp, tree_by_spp = tree_by_spp, opuntia_combined = opuntia_combined) format_tables_EDIT_style(data = indicator_descriptions[,-1], caption = "Indicators used in ordination and clustering analysis to develop plant communities") %>% collapse_rows(columns = 1, valign = "top")
# are the data clusterable? Use the Hopkin's statistic (Lawson, Richard, and Jurs, 1990; see https://www.datanovia.com/en/lessons/assessing-clustering-tendency/ for explanation) # H values > 0.5 indicate clustering tendency (i.e. if H is below 0.5, the data set is uniformly distributed and clusters are not meaningful) res_LoamyUplands <- factoextra::get_clust_tendency(data = ord.df, n=nrow(ord.df)-1, graph = FALSE) res_LoamyUplands$hopkins_stat
# cluster and indicator species analysis library(cluster) # for cluster analysis library(dendextend) library(labdsv) library(foreach) library(vegan) library(ggplot2) # initial clustering dist.df <- vegan::vegdist(ord.df, method = "bray", na.rm = T) # creates Bray-Curtis distance matrix
# This chunk determines the best combination of number of clusters and fuzziness for the data. It is SLOOOOWWW so better to run it before knitting to get the right clustering parameters and then not run it when knitting the descriptive document. library(parallel) library(doParallel) # get index values for all combos of membership exponent and number of clusters indices_arg <- c("ch", #"duda", won't work with the ranking scheme "cindex", #"gamma", # takes a long time #"beale", # unclear how this one works "ccc", "ptbiserial", #"gplus", # takes a long time "db", "silhouette", "dunn" #, #"gap" # takes a long time ) indices_name <- c("CH", #"Duda", "Cindex", #"gamma", # takes a long time #"Beale", "CCC", "Ptbiserial", #"gplus", # takes a long time "DB", "Silhouette", "Dunn" #, #"gap" # takes a long time ) memb_exps <- seq(1.1, 1.3, # start with 2.0, can reduce value to rerun if some larger values don't give results by = 0.1) cl <- makeCluster(spec = 3, # number of cores to use type = "PSOCK", methods = FALSE) registerDoParallel(cl) cluster_eval_fanny <- foreach(memb_exp=memb_exps, .packages = c("STMdevelopment", "dplyr", "cluster"), .combine = "rbind") %dopar% { nbmetrics_fanny_temp <- NbClust_fanny(data = ord.df, diss = dist.df, #chord.df, distance = NULL, min.nc = 2, max.nc = 12, method = "fanny", index = "all", # could also try alllong to get a few more metrics - this is SLOW alphaBeale = 0.1, memb.exp= memb_exp # 1.2 ) cluster_eval_fanny_temp <- nbmetrics_fanny_temp$All.index %>% as.data.frame() %>% mutate(Number_clusters = rownames(.), Memb_exp = memb_exp) } registerDoSEQ() stopCluster(cl) cluster_eval_fanny_indices <- select(cluster_eval_fanny, any_of(c("Memb_exp", "Number_clusters", indices_name))) # calculate percent worst than best value for each memb. exp. X num. clusters combo for each index pct_off_best_higher <- function(x, all_vals){ best <- max(all_vals) pct_off <- ((best-x)/best)*100 return(pct_off) } pct_off_best_lower <- function(x, all_vals){ best <- min(all_vals) pct_off <- ((x-best)/best)*100 return(pct_off) } cluster_eval_fanny_indices_rank <- cluster_eval_fanny_indices %>% rowwise() %>% mutate(CH_rank = pct_off_best_higher(CH, cluster_eval_fanny_indices$CH), Cindex_rank = pct_off_best_lower(Cindex, cluster_eval_fanny_indices$Cindex), CCC_rank = pct_off_best_higher(CCC, cluster_eval_fanny_indices$CCC), Ptbiserial_rank = pct_off_best_higher(Ptbiserial, cluster_eval_fanny_indices$Ptbiserial), DB_rank = pct_off_best_lower(DB, cluster_eval_fanny_indices$DB), Silhouette_rank = pct_off_best_higher(Silhouette, cluster_eval_fanny_indices$Silhouette), Dunn_rank = pct_off_best_higher(Dunn, cluster_eval_fanny_indices$Dunn) ) %>% mutate(overall_rank = sum(CH_rank, # lowest overall rank wins! Cindex_rank, CCC_rank, Ptbiserial_rank, DB_rank, Silhouette_rank, Dunn_rank, na.rm = T)) best_n_clust <- as.numeric(cluster_eval_fanny_indices_rank[[which.min(cluster_eval_fanny_indices_rank$overall_rank), "Number_clusters"]]) best_memb_exp <- as.numeric(cluster_eval_fanny_indices_rank[[which.min(cluster_eval_fanny_indices_rank$overall_rank), "Memb_exp"]]) # VISUALLY CHECK for ties and ranking weirdness before proceeding!!
k_groups_fuzzy <- 8 # best is 8 (rank 74), then 9 (rank 87) then 10 (rank 105) then 7 (rank 108) - all with 1.1 memb. exp. memb_exp <- 1.1 #1.1 set.seed(1) clust_fuzz <- fanny(x=dist.df, k = k_groups_fuzzy, diss = T, #metric = "euclidean", memb.exp = # this affects how crisp/fuzzy the clusters are - closer to 1 is crisper memb_exp ) grp_fuzz <- clust_fuzz$clustering plot_data_first$PlantCommunity_fuzzy_original <- grp_fuzz # double check groups with original grouping - seems like changing the file path # to the data reset the grouping algorithm somehow so that the same groups are # made but they are given different numbers by fanny(). If using the Duniway et # al 2024 data, groups should line up like this: # Juniper-pinyon woodland, n = 229. Originally cluster #1. Workshop feedback was to lump with 5 Pinyon-juniper woodland # Open juniper-pinyon shrubland, n = 238. Originally cluster #2. # C4 grassland, n = 110. Originally cluster #3. # Annualized juniper savanna, n=113. Originally cluster #4 # Pinyon-juniper woodland, n = 270. Originally cluster #5. Workshop feedback was to lump with 1 Juniper-pinyon woodland # Gambel oak woodland, n=112. Originally cluster #6. # Sagebrush shrubland/grassland, n=198. Originally cluster #7. # Bare juniper savanna, n=186. Originally cluster #8. # Based on matching up clusters by number of plots, I'll reassign numbers to match # the original groupings: updated_cluster_numbers <- data.frame("PlantCommunity_fuzzy_original" = c(1, 2, 3, 4, 5, 6, 7, 8), "PlantCommunity_fuzzy" = c(1, 4, 5, 2, 6, 7, 3, 8)) plot_data_first <- left_join(plot_data_first, updated_cluster_numbers, by="PlantCommunity_fuzzy_original")
Plots were clustered into ecological communities using the FANNY fuzzy clustering algorithm (Kaufman and Rousseeuw, 1990) in the "cluster" R package (Maechler et al., 2021) with Bray-Curtis distance ("vegan" package; Oksanen et al., 2020). A suite of clustering indices were used to determine the most appropriate number of clusters and degree of fuzziness (i.e. membership exponent), including Calinski and Harabasz (1974), C-index (Hubert and Levin, 1976), cubic clustering criteria (Sarle, 1983), point-biserial correlation (Milligan 1980, 1981), Davies and Bouldin (1979), silhouette width (Rousseeuw 1987), and Dunn's index (1974). Values between 1 and 12 possible clusters and 1.1 to 1.3 membership exponents were tested. The final clustering solution used r k_groups_fuzzy clusters with a membership exponent of r memb_exp.
# another slow chunk - run before knitting and then don't evaluate in the knit step. NMDS_scree(ord.df)
ord_dims <- 3 # change manually based on the plot results #run NMS ordination set.seed(1) ord <- metaMDS(ord.df, k=ord_dims, # number of dimensions trymax = 30)
The clustering results were visualized in ordination spaced, using non-metric multidimensional scaling (NMDS) with Bray-Curtis distance ("vegan" package in R). Based on testing ordination stress for 1 to 5 dimensions, ordination stress is acceptable with r ord_dims dimensions, and stress reduction is minimal with additional dimensions. The stress for the ordination with r ord_dims dimensions is r round(ord$stress, 3).
# 2D plotting pal_veg2 <- RColorBrewer::brewer.pal(n=k_groups_fuzzy, name = "Dark2") # # manually update based on indicator species # group_labels_fuzzy <- c(paste0("Juniper-pinyon # woodland (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==1)), ")"), # paste0("Open juniper-pinyon # shrubland (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==2)), ")"), # paste0("C4 grassland # (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==3)), ")"), # paste0("Annualized juniper # savanna (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==4)), ")"), # paste0("Pinyon-juniper # woodland (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==5)), ")"), # paste0("Gambel oak # woodland (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==6)), ")"), # paste0("Sagebrush shrubland/ # grassland (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==7)), ")"), # paste0("Bare juniper # savanna (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==3)), ")") # ) group_labels_base <- c(paste0("Juniper-pinyon woodland"), paste0("Open juniper-pinyon shrubland"), paste0("C4 grassland "), paste0("Annualized juniper savanna"), paste0("Pinyon-juniper woodland"), paste0("Gambel oak woodland"), paste0("Sagebrush shrubland/ grassland"), paste0("Bare juniper savanna)") ) group_labels_fuzzy <- paste0(group_labels_base, " (n = ", c(nrow(filter(plot_data_first, PlantCommunity_fuzzy==1)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==2)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==3)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==4)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==5)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==6)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==7)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==8)) ), ")") #group_labels_fuzzy <- as.character(1:clust_fuzz$k.crisp) par(mfrow=c(2,2)) # Axes 1x2 plot(ord, choices = c(1,2), type = "n", # plot the axes xlim = c(-1, 1), ylim = c(-1, 1)) points(ord, choices = c(1,2), display = "sites", # plot points - can choose "sites" or "species" col=pal_veg2[plot_data_first$PlantCommunity_fuzzy], pch = 21, cex = .6, bg=pal_veg2[plot_data_first$PlantCommunity_fuzzy]) ordiellipse(ord, plot_data_first$PlantCommunity_fuzzy, col=pal_veg2, lwd = 2, label = F, choices = c(1,2)) # plot your groups # can use oriellipse, orihull, or ordispider to plot groups, depending on your needs # Axes 3x2 plot(ord, choices = c(3,2), type = "n", xlim = c(-1, 1), ylim = c(-1, 1)) points(ord, choices = c(3,2), display = "sites", col=pal_veg2[plot_data_first$PlantCommunity_fuzzy], pch = 21, cex = .6, bg=pal_veg2[plot_data_first$PlantCommunity_fuzzy]) ordiellipse(ord, plot_data_first$PlantCommunity_fuzzy, col=pal_veg2, lwd = 2, label = F, choices = c(3,2)) # Axes 1x3 plot(ord, choices = c(1,3), type = "n", xlim = c(-1, 1), ylim = c(-1, 1)) points(ord, choices = c(1,3), display = "sites", col=pal_veg2[plot_data_first$PlantCommunity_fuzzy], pch = 21, cex = .6, bg=pal_veg2[plot_data_first$PlantCommunity_fuzzy]) ordiellipse(ord, plot_data_first$PlantCommunity_fuzzy, col=pal_veg2, lwd = 2, label = F, choices = c(1,3)) # Legend plot(ord, type = "n", axes=FALSE, display = "sites", col=pal_veg2[plot_data_first$PlantCommunity_fuzzy], xlab = "", ylab = "") legend(x="center", legend = group_labels_fuzzy, fill = pal_veg2, title = "Plant communities", bty = "n", y.intersp = 1.5) par(mfrow=c(1,1))
# save clustering results fuzzy_cluster_membership <- as.data.frame(clust_fuzz$membership) fuzzy_cluster_membership$PlotCode <- row.names(clust_fuzz$membership) fuzzy_cluster_membership$best_group <- clust_fuzz$clustering # update the cluster numbers to match the original groupings (see note above in the fuzzy-cluster chunk) fuzzy_cluster_membership <- left_join(fuzzy_cluster_membership, updated_cluster_numbers, by = c("best_group" = "PlantCommunity_fuzzy_original")) # # merge groups 1 and 5 based on workshop feedback # fuzzy_cluster_membership <- mutate(fuzzy_cluster_membership, # best_group = ifelse(best_group==1|best_group==5, # yes = 9, # no = best_group)) colnames(fuzzy_cluster_membership) <- c(paste0("Group_", 1:8, "_membership"), "PlotCode", "Best_group") # group_names <- data.frame(Best_group = 1:8, # Best_group_name = c("Juniper-pinyon woodland", # "Open juniper-pinyon shrubland", # "C4 grassland", # "Annualized juniper savanna", # "Pinyon-juniper woodland", # "Gambel oak woodland", # "Sagebrush shrubland/grassland", # "Bare juniper savanna")) # update with descriptive names after examining the groups # reorder group names due to user=Duniway_et_al_2024 mix up group_names <- data.frame(Best_group = 1:8, Best_group_name = c("Juniper-pinyon woodland", "Annualized juniper savanna", "Pinyon-juniper woodland", "Open juniper-pinyon shrubland", "Gambel oak woodland", "Sagebrush shrubland/grassland", "C4 grassland", "Bare juniper savanna")) # update with descriptive names after examining the groups # group_names <- data.frame(Best_group = 1:9, # Best_group_name = c("Juniper-pinyon woodland", # "Open juniper-pinyon shrubland", # "C4 grassland", # "Annualized juniper savanna", # "Pinyon-juniper woodland", # "Gambel oak woodland", # "Sagebrush shrubland/grassland", # "Bare juniper savanna", # "Pinyon-juniper woodland")) fuzzy_cluster_membership <- left_join(fuzzy_cluster_membership, group_names) # get the plot location and date sampled (use first date) plot_locations_first <- plot_locations %>% mutate(DateSampled = as.Date(paste(Year, Month, Day, sep = "-"))) %>% arrange(DateSampled) %>% group_by(PlotCode) %>% summarise(DateSampled = first(DateSampled)) fuzzy_cluster_memb_sf <- left_join(fuzzy_cluster_membership, plot_locations_first) %>% filter(PlotCode != "AIM_Colorado NWDO Little Snake FO 2017_1024" & PlotCode != "AIM_UT Color Country DO 2019_01") %>% # the geometry of these plot locations is messed up somehow mutate(post_workshop1_group = ifelse(test = Best_group==1|Best_group==5, yes = 9, no = Best_group)) # # save plots - usually only need to run this once # sf::st_write(fuzzy_cluster_memb_sf, # dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps/contains_NRI_do_not_sync", # layer = paste0("SW_ShallowDeepRocky_fuzzyclusters_", Sys.Date(), ".shp"), # driver = "ESRI Shapefile") # # fuzzy_cluster_memb_sf_shareable <- filter(fuzzy_cluster_memb_sf, !grepl(pattern = "NRI", x=PlotCode)) # # sf::st_write(fuzzy_cluster_memb_sf_shareable, # dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps", # layer = paste0("SW_ShallowDeepRocky_fuzzyclust_share_", Sys.Date(), ".shp"), # driver = "ESRI Shapefile") # # states <- sf::st_read("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/GIS/Political boundaries/cb_2018_us_state_5m/cb_2018_us_state_5m.shp") # fourcorners <- filter(states, NAME %in% c("Utah", "Colorado", "Arizona", "New Mexico")) # # fuzzy_cluster_memb_sf_mike <- sf::st_filter(sf::st_as_sf(fuzzy_cluster_memb_sf, crs=sf::st_crs(plot_locations)), # sf::st_transform(fourcorners, crs=sf::st_crs(plot_locations))) %>% # mutate(Source = ifelse(test = grepl(x=PlotCode, pattern = "^AIM"), # yes = "AIM", # no = NA)) %>% # mutate(Source = ifelse(test = grepl(x=PlotCode, pattern = "^LMF"), # yes = "LMF", # no = Source)) %>% # mutate(Source = ifelse(test = grepl(x=PlotCode, pattern = "^IM_NCPN"), # yes = "NPS", # no = Source)) %>% # mutate(Source = ifelse(test = grepl(x=PlotCode, pattern = "^NRI"), # yes = "NRI", # no = Source)) # # sf::st_write(fuzzy_cluster_memb_sf_mike, # dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps/contains_NRI_do_not_sync", # layer = paste0("SW_ShallowDeepRocky_fuzzyclust_mike_", Sys.Date(), ".shp"), # driver = "ESRI Shapefile")
State and transition model {.tabset}
CUSTOM DIAGRAM
knitr::include_graphics("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Analyses/Semiarid_Warm_ShallowDeepRocky_boxandarrow_workshop1.JPG")
# TODO legend for custom diagram if needed # scratch pad: # All states described here contain some introduced species (e.g. *Bromus tectorum*; see Tables \@ref(tab:comm1-1-dominant-spp), \@ref(tab:comm2-1-dominant-spp), \@ref(tab:comm3-1-dominant-spp), \@ref(tab:comm4-1-dominant-spp), and \@ref(tab:comm5-1-dominant-spp) for details).
The data used to develop this STM were sourced from large monitoring data sets designed to sample landscapes with a variety of current and historic land use. The model does not describe reference communities because these data do not include targeted sampling of relict plant communities with histories of little alteration to disturbance regimes. Several states described here contain some introduced species (e.g. Bromus tectorum; see Tables \@ref(tab:comm3-1-dominant-spp) and \@ref(tab:comm4-1-dominant-spp) for details).
For furthur information on reference states as well as drivers of transitions, see the Ecological Site Descriptions for the sites listed in Table \@ref(tab:assoc-eco-sites).
STANDARD DIAGRAM
Click on state and transition labels to scroll to the respective text.
TODO create boxes for each state that link to the state description below
State descriptive figures
# load descriptive data plot_data_descriptive <- plot_data_pull(user=user, target_ESG = target_ESG, data_sources = data_sources, indicators = c("TotalFoliarCover", "AH_WoodyCover", "AH_ShrubCover", # TODO what to do about subshrubs?? "AH_TreeCover", "AH_ForbGrassCover", "AH_PerenForbGrassCover", "AH_AnnForbGrassCover", "SoilStab_all", "AH_ArtemisiaTridentataCover", "CA_percent_200plus", "CA_percent_100plus" ), ann_grass_by_spp = T, ann_forb_by_spp = T, per_grass_by_spp = T, per_forb_by_spp = T, succulent_by_spp = T, shrub_by_spp = T, subshrub_by_spp = T, tree_by_spp = T, opuntia_combined = F, impute_gap_type = c("CA_percent_100plus", "CA_percent_200plus") ) %>% filter(!is.na(Month)) %>% left_join(select(plot_data_first, PlotCode, Year, PlantCommunity_fuzzy) , .) colnames(plot_data_descriptive)[which(colnames(plot_data_descriptive)=="AH_ArtemisiaTridentataCover")] <- "ARTR2" # temporary patch - some plot codes are duplicated but have different primary keys filter(plot_data_first, PlotCode=="AIM_Colorado NWDO Little Snake FO 2017_1024" | PlotCode=="AIM_UT Color Country DO 2019_01")$PlotID plot_data_descriptive <- filter(plot_data_descriptive, !(PlotCode=="AIM_Colorado NWDO Little Snake FO 2017_1024" & PlotID!="1705311149423160") & !(PlotCode=="AIM_UT Color Country DO 2019_01" & PlotID!="1906201000185195")) plot_locations <- cbind(plot_locations, sf::st_coordinates(plot_locations)) %>% filter(!(PlotCode=="AIM_Colorado NWDO Little Snake FO 2017_1024" & X!=filter(plot_data_descriptive, PlotCode=="AIM_Colorado NWDO Little Snake FO 2017_1024")$Longitude_NAD83) & !(PlotCode=="AIM_UT Color Country DO 2019_01" & X!=filter(plot_data_descriptive, PlotCode=="AIM_UT Color Country DO 2019_01")$Longitude_NAD83)) # summarize species occurance descriptive.df <- dplyr::select(plot_data_descriptive, -SourceKey, -PlotID, -SiteName, -PlotName, -Year, -Longitude_NAD83, -Latitude_NAD83, -Month, -Day, -PlantCommunity_fuzzy, -TotalFoliarCover, -AH_WoodyCover, -AH_ShrubCover, -AH_TreeCover, -AH_ForbGrassCover, -AH_PerenForbGrassCover, -AH_AnnForbGrassCover, -SoilStab_all, -CA_percent_200plus, -CA_percent_100plus ) %>% tibble::column_to_rownames("PlotCode") # keep an ID code for the plot as the row name to prevent data scrambling problems sp_pa_desc <- mutate(rowwise(descriptive.df), across(everything(), function(x){if(x>0){1}else{0}})) # make a presence/absence data frame sp_keep_desc <- names(colSums(sp_pa_desc)[which(colSums(sp_pa_desc)>0)]) descriptive.df <- select(descriptive.df, all_of(sp_keep_desc)) plot_data_descriptive <- select(plot_data_descriptive, SourceKey, PlotID, SiteName, PlotName, Year, PlotCode, Longitude_NAD83, Latitude_NAD83, Month, Day, PlantCommunity_fuzzy, TotalFoliarCover, AH_WoodyCover, AH_ShrubCover, AH_TreeCover, AH_ForbGrassCover, AH_PerenForbGrassCover, AH_AnnForbGrassCover, SoilStab_all, CA_percent_100plus, CA_percent_200plus, all_of(sp_keep_desc)) #ord.df <- select(ord.df, -any_of(c("UNKS", "UNKPG", "UNKAF", "UNKPF", "PG1", "PG01", "AF1"))) ## load climate data AI_raster <- raster::raster(file_paths$ai_raster) plot_AI <- sf::st_as_sf(raster::extract(x = AI_raster, y = plot_locations, sp = T)) %>% mutate(AI = ai_et0_WesternUSEcoregions*0.0001 # AI rasters are actually AI*10,000 so that they #can be stored as integers - need to put back into standard AI units ) %>% sf::st_drop_geometry() rm(AI_raster) climate_rasters <- raster::stack(file_paths$prism_rasters) prism_vals <- raster::extract(x=climate_rasters, y=plot_locations, sp = T) %>% as(., "sf") sf::st_geometry(prism_vals) <- NULL colnames(prism_vals)[which(colnames(prism_vals) %in% c("PRISM_ppt_30yr_normal_800mM2_annual_asc", "PRISM_tmean_30yr_normal_800mM2_annual_asc", "PRISM_tmax_30yr_normal_800mM2_annual_asc", "PRISM_tmin_30yr_normal_800mM2_annual_asc" ))] <- c("Ann_ppt_mm", "Ann_tmean", "Ann_tmax", "Ann_tmin") rm(climate_rasters) plot_data_descriptive <- left_join(plot_data_descriptive, select(plot_AI, PlotCode, Year, AI)) %>% left_join(., select(prism_vals, PlotCode, Year, Ann_ppt_mm, Ann_tmean, Ann_tmax, Ann_tmin)) ## species richness and diversity plot_data_descriptive$shannon_diversity <- vegan::diversity(x=descriptive.df, index = "shannon") plot_data_descriptive$sp_richness <- vegan::specnumber(x=descriptive.df) # # AERO # aero <- read.csv(file_paths$aero) # plot_aero <- aero %>% # # AIM error patch # filter(!(ProjectName=="Colorado NWDO Little Snake FO 2017" & PlotID=="1024" & # Longitude_NAD83!=filter(plot_data_descriptive, # PlotCode=="AIM_Colorado NWDO Little Snake FO 2017_1024")$Longitude_NAD83) # & !(ProjectName=="UT Color Country DO 2019" & PlotID=="01" & # Longitude_NAD83!=filter(plot_data_descriptive, # PlotCode=="AIM_UT Color Country DO 2019_01")$Longitude_NAD83)) %>% # select(ProjectName, PlotID, horizontal_flux_total_MN, horizontal_flux_total_LPI, horizontal_flux_total_UPI) %>% # setNames(c("SiteName", "PlotName", "horizontal_flux_total_MN", "horizontal_flux_total_LPI", "horizontal_flux_total_UPI")) # # plot_data_descriptive <- left_join(plot_data_descriptive, plot_aero) # AERO aero <- read.csv(file_paths$aero) aero_aim <- filter(aero, source=="AIM") %>% mutate(PlotKey = substr(PrimaryKey, 1, nchar(PrimaryKey)-10), # fix the Jornada's data storage errors Year = as.numeric(substr(PrimaryKey, nchar(PrimaryKey)-9, nchar(PrimaryKey)-6)),) %>% select(PlotKey, Year, horizontal_flux_total_MN, horizontal_flux_total_LPI, horizontal_flux_total_UPI) %>% rename("PlotID" = "PlotKey") aero_lmf <- filter(aero, source=="LMF") %>% select(PlotKey, horizontal_flux_total_MN, horizontal_flux_total_LPI, horizontal_flux_total_UPI) %>% rename("PlotID" = "PlotKey") plot_data_descriptive <- left_join(plot_data_descriptive, aero_aim) %>% rows_update(x = ., y = aero_lmf, unmatched = "ignore") # RHEM rhem_aim <- read.csv(file_paths$rhem$aim) %>% mutate(PrimaryKey_Jornada = as.character(PrimaryKey), Year = as.numeric(substr(PrimaryKey, nchar(PrimaryKey)-9, nchar(PrimaryKey)-6)), PlotID = substr(PrimaryKey, 1, nchar(PrimaryKey)-10)) %>% # make the Jornada's proprietary PrimaryKey system play nice with everyone else's select(PlotID, Year, Runoff_Long_Term_MEAN, Soil_Loss_Long_Term_MEAN, Sed_Yield_Long_Term_MEAN) %>% distinct() rhem_lmf <- read.csv(file_paths$rhem$lmf) %>% select(PrimaryKey, Runoff_Long_Term_MEAN, Soil_Loss_Long_Term_MEAN, Sed_Yield_Long_Term_MEAN) %>% distinct() %>% rename("PlotID" = "PrimaryKey") rhem_nri <- read.csv(file_paths$rhem$nri) %>% select(PrimaryKey, Runoff_Long_Term_MEAN, Soil_Loss_Long_Term_MEAN, Sed_Yield_Long_Term_MEAN) %>% distinct() %>% rename("PlotID" = "PrimaryKey") plot_data_descriptive <- left_join(plot_data_descriptive, rhem_aim) %>% rows_update(x = ., y = rhem_lmf, unmatched = "ignore") %>% rows_update(x = ., y = rhem_nri, unmatched = "ignore")
# # consider including species present in >10% of plots: # names(colSums(sp_pa)[which(colSums(sp_pa)>(nrow(plot_data_first)*0.10))]) # # also consider including species with high mean cover across all plots: # mean_cover <- summarise(plot_data_first, across(.cols = AH_C3IntroducedPerenGrassCover:AH_OpuntiaCover, .fns = mean)) %>% # select(-any_of(c("Longitude_NAD83", "Latitude_NAD83", "Year", "Month", "Day"))) %>% # t() %>% # as.data.frame() # # row.names(mean_cover)[which(mean_cover$V1>=0.9)] # Consider switching to plotly to allow users to hover and see exact quantile and outlier values, e.g.: # plotly::plot_ly(x=plot_data_descriptive$PlantCommunity_fuzzy, # y=plot_data_descriptive$horizontal_flux_total_UPI, # split=plot_data_descriptive$PlantCommunity_fuzzy, # type='box') plot_data_first_tall <- select(plot_data_first, SourceKey, PlotID, SiteName, PlotName, Year, PlotCode, #PlantCommunity, PlantCommunity_fuzzy, AH_C3IntroducedPerenGrassCover, AH_C3NativePerenGrassCover, AH_C4NativePerenGrassCover, AH_IntroducedAnnForbCover, AH_IntroducedAnnGrassCover, AH_IntroducedPerenForbCover, AH_NativeAnnForbCover, AH_NativeAnnGrassCover, AH_NativePerenForbCover, BareSoilCover, FH_TotalLitterCover, AH_ArtemisiaTridentataCover, #CA_percent_200plus, AH_C4IntroducedPerenGrassCover, FH_LichenMossCover, #CA_percent_100plus, EPVI, CEMO2, CHVI8, PUTR2, AMUT, QUGA, GUSA2, JUOS, PIED ) %>% left_join(., select(plot_data_descriptive, PlotCode, CA_percent_100plus, CA_percent_200plus)) %>% tidyr::pivot_longer(cols = c(AH_C3IntroducedPerenGrassCover, AH_C3NativePerenGrassCover, AH_C4NativePerenGrassCover, AH_IntroducedAnnForbCover, AH_IntroducedAnnGrassCover, AH_IntroducedPerenForbCover, AH_NativeAnnForbCover, AH_NativeAnnGrassCover, AH_NativePerenForbCover, BareSoilCover, FH_TotalLitterCover, AH_ArtemisiaTridentataCover, CA_percent_200plus, AH_C4IntroducedPerenGrassCover, FH_LichenMossCover, CA_percent_100plus, EPVI, CEMO2, CHVI8, PUTR2, AMUT, QUGA, GUSA2, JUOS, PIED ), names_to = "CoverType", values_to = "PctCover") spp_list <- read.csv(file_paths$species_list) cover_funcgrps <- sort(c("AH_C3IntroducedPerenGrassCover", "AH_C3NativePerenGrassCover", "AH_C4NativePerenGrassCover", "AH_IntroducedAnnForbCover", "AH_IntroducedAnnGrassCover", "AH_IntroducedPerenForbCover", "AH_NativeAnnForbCover", "AH_NativeAnnGrassCover", "AH_NativePerenForbCover", "BareSoilCover", "FH_TotalLitterCover", "AH_ArtemisiaTridentataCover", "CA_percent_200plus", "AH_C4IntroducedPerenGrassCover", "FH_LichenMossCover", "CA_percent_100plus")) funcgrps_names <- indicator_descriptions %>% bind_rows(., data.frame(Indicator_code = c("CA_percent_100plus", "CA_percent_200plus"), Indicator = c("Ann. & peren. canopy gaps > 100 cm", "Ann. & peren. canopy gaps > 200 cm"))) %>% filter(., Indicator_code %in% cover_funcgrps) %>% arrange(., Indicator_code) %>% select(Indicator) %>% unlist() names(cover_funcgrps) <- funcgrps_names cover_spp <- sort(c("EPVI", "CEMO2", "CHVI8", "PUTR2", "AMUT", "QUGA", "GUSA2", "JUOS", "PIED" )) spp_names <- paste0( filter(spp_list, SpeciesCode %in% cover_spp & SpeciesState=="AIM")$ScientificName, " (", filter(spp_list, SpeciesCode %in% cover_spp & SpeciesState=="AIM")$CommonName, ")") names(cover_spp) <- spp_names cover_labels <- sort(c(cover_funcgrps, cover_spp)) functypes_grps <- c("TotalFoliarCover", #"AH_WoodyCover", #AH_ForbGrassCover, "AH_ShrubCover", "AH_TreeCover", "AH_PerenForbGrassCover", "AH_AnnForbGrassCover" ) functypes_cover_states <- ggplot(data = select(plot_data_descriptive, any_of(c("PlantCommunity_fuzzy", functypes_grps))) %>% tidyr::pivot_longer(cols = all_of(functypes_grps), names_to = "CoverType", values_to = "PctCover"), aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) + geom_boxplot() + scale_fill_brewer(type = "qual", palette = "Set3", name = "Functional type", labels = c("Annual herbaceous", "Perennial herbaceous", #"Woody", "Shrub", "Tree", "Total foliar cover") ) + scale_x_discrete(labels=group_labels_fuzzy) + xlab(NULL) + ylab("Cover (%)") + ggtitle("Plant functional groups") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black"), axis.text.y = element_text(color = "black")) woody_grps <- c("AH_ArtemisiaTridentataCover", "EPVI", "CEMO2", "CHVI8", "PUTR2", "AMUT", "QUGA", "GUSA2", "JUOS", "PIED" ) woody_cover_states <- ggplot(data = filter(plot_data_first_tall, CoverType %in% woody_grps), aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) + geom_boxplot() + scale_fill_brewer(type = "qual", palette = "Set3", name = "Species", labels = names(cover_labels[which(cover_labels %in% woody_grps)])) + scale_x_discrete(labels=group_labels_fuzzy) + xlab(NULL) + ylab("Cover (%)") + ggtitle("Woody plants") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black"), axis.text.y = element_text(color = "black")) perenherb_grps <- c("AH_C3IntroducedPerenGrassCover", "AH_C3NativePerenGrassCover", "AH_C4IntroducedPerenGrassCover", "AH_C4NativePerenGrassCover", "AH_IntroducedPerenForbCover", "AH_NativePerenForbCover") perenherb_cover_states <- ggplot(data = filter(plot_data_first_tall, CoverType %in% perenherb_grps), aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) + geom_boxplot() + scale_fill_brewer(type = "qual", palette = "Set3", name = "Functional group", labels = names(cover_labels[which(cover_labels %in% perenherb_grps)])) + scale_x_discrete(labels=group_labels_fuzzy) + xlab(NULL) + ylab("Cover (%)") + ggtitle("Perennial herbaceous plants") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black"), axis.text.y = element_text(color = "black")) annherb_grps <- c("AH_IntroducedAnnForbCover", "AH_IntroducedAnnGrassCover", "AH_NativeAnnForbCover", "AH_NativeAnnGrassCover") annherb_cover_states <- ggplot(data = filter(plot_data_first_tall, CoverType %in% annherb_grps), aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) + geom_boxplot() + scale_fill_brewer(type = "qual", palette = "Set3", name = "Functional group", labels = names(cover_labels[which(cover_labels %in% annherb_grps)])) + scale_x_discrete(labels=group_labels_fuzzy) + xlab(NULL) + ylab("Cover (%)") + ggtitle("Annual herbaceous plants") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black"), axis.text.y = element_text(color = "black")) ground_grps <- c("BareSoilCover", "FH_TotalLitterCover", "CA_percent_200plus", "FH_LichenMossCover", "CA_percent_100plus") ground_cover_states <- ggplot(data = filter(plot_data_first_tall, CoverType %in% ground_grps), aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) + geom_boxplot() + scale_fill_brewer(type = "qual", palette = "Set3", name = "Functional group", labels = names(cover_labels[which(cover_labels %in% ground_grps)])) + scale_x_discrete(labels=group_labels_fuzzy) + xlab("Plant community") + ylab("Cover (%)") + ggtitle("Ground cover") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black"), axis.text.y = element_text(color = "black")) cover_states <- cowplot::plot_grid(functypes_cover_states + theme(legend.position = "none"), woody_cover_states + theme(legend.position = "none"), perenherb_cover_states + theme(legend.position = "none"), annherb_cover_states + theme(legend.position = "none"), ground_cover_states + theme(legend.position = "none"), ncol = 1, axis = "lr", labels = "auto") cover_axes <- cowplot::plot_grid(cowplot::get_legend(functypes_cover_states), cowplot::get_legend(woody_cover_states), cowplot::get_legend(perenherb_cover_states), cowplot::get_legend(annherb_cover_states), cowplot::get_legend(ground_cover_states), ncol = 1) cowplot::plot_grid(cover_states, cover_axes, ncol = 2, rel_widths = c(1, 0.5))
group_labels_fuzzy <- paste0(group_labels_base, " (n = ", c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1 & !is.na(AI))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2 & !is.na(AI))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3 & !is.na(AI))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4 & !is.na(AI))) ), ")") cowplot::plot_grid( cowplot::plot_grid( ggplot(data = plot_data_descriptive, aes(x=PlantCommunity_fuzzy, y=AI)) + geom_boxplot() + scale_x_discrete(labels=group_labels_fuzzy) + xlab(NULL) + ylab("Aridity index") + ylim(c(0, 0.45)) + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)), ggplot(data = plot_data_descriptive, aes(x=PlantCommunity_fuzzy, y=Ann_ppt_mm)) + geom_boxplot() + scale_x_discrete(labels=group_labels_fuzzy) + xlab(NULL) + ylab("Annual precipitation (mm)") + ylim(c(0, 650)) + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)), nrow = 1, labels = c("a", "b")), ggplot(data = select(plot_data_descriptive, PlantCommunity_fuzzy, Ann_tmean, Ann_tmax, Ann_tmin) %>% tidyr::pivot_longer(c(Ann_tmean, Ann_tmax, Ann_tmin), names_to = "temptype", values_to = "tempdegrees"), aes(fill=PlantCommunity_fuzzy, y=tempdegrees, x=temptype)) + geom_boxplot() + scale_fill_manual(name="Plant community", labels=group_labels_fuzzy, values = pal_veg2) + scale_x_discrete(labels=c("Maximum", "Mean", "Minimum")) + xlab(NULL) + ylab("Temperature (C)") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)), nrow = 2, labels = c(NA, "c"))
group_labels_fuzzy <- paste0(group_labels_base, " (n = ", c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1 & !is.na(SoilStab_all))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2 & !is.na(SoilStab_all))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3 & !is.na(SoilStab_all))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4 & !is.na(SoilStab_all))) ), ")") ggplot(data = plot_data_descriptive, aes(x=PlantCommunity_fuzzy, y=SoilStab_all)) + geom_boxplot() + scale_x_discrete(labels=group_labels_fuzzy) + ylim(c(0, 6)) + xlab("Plant community") + ylab("Mean soil stability rating") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1))
group_labels_fuzzy <- paste0(group_labels_base, " (n = ", c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1 & !is.na(shannon_diversity))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2 & !is.na(shannon_diversity))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3 & !is.na(shannon_diversity))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4 & !is.na(shannon_diversity))) ), ")") cowplot::plot_grid(ggplot(data = plot_data_descriptive, aes(x=PlantCommunity_fuzzy, y=shannon_diversity)) + geom_boxplot() + scale_x_discrete(labels=group_labels_fuzzy) + xlab("Plant community") + ylab("Shannon diversity index") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)), ggplot(data = plot_data_descriptive, aes(x=PlantCommunity_fuzzy, y=sp_richness)) + geom_boxplot() + scale_x_discrete(labels=group_labels_fuzzy) + xlab("Plant community") + ylab("Number of species") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)), nrow = 1, labels = "auto")
# dust_tall <- plot_data_descriptive %>% # select(PlotID, PlantCommunity_fuzzy, horizontal_flux_total_MN, horizontal_flux_total_LPI, horizontal_flux_total_UPI) %>% # tidyr::pivot_longer(cols = c(horizontal_flux_total_MN, horizontal_flux_total_LPI, horizontal_flux_total_UPI), # names_to = "dust_statistic", # values_to = "dust_flux") group_labels_fuzzy <- paste0(group_labels_base, " (n = ", c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1 & !is.na(horizontal_flux_total_MN))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2 & !is.na(horizontal_flux_total_MN))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3 & !is.na(horizontal_flux_total_MN))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4 & !is.na(horizontal_flux_total_MN))) ), ")") ggplot(data = plot_data_descriptive, aes(x=PlantCommunity_fuzzy, y=horizontal_flux_total_MN)) + geom_boxplot() + scale_x_discrete(labels=group_labels_fuzzy) + #ylim(c(0, 12)) + xlab("Plant community") + ylab("Mean horizontal dust flux") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) # TODO add log transformation to y axis
group_labels_fuzzy <- paste0(group_labels_base, " (n = ", c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1 & !is.na(Runoff_Long_Term_MEAN))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2 & !is.na(Runoff_Long_Term_MEAN))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3 & !is.na(Runoff_Long_Term_MEAN))), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4 & !is.na(Runoff_Long_Term_MEAN))) ), ")") cowplot::plot_grid(ggplot(data = plot_data_descriptive, aes(x=PlantCommunity_fuzzy, y=Runoff_Long_Term_MEAN)) + geom_boxplot() + scale_x_discrete(labels=group_labels_fuzzy) + xlab("Plant community") + ylab("Long-term mean runoff (mm/yr)") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)), ggplot(data = plot_data_descriptive, aes(x=PlantCommunity_fuzzy, y=Soil_Loss_Long_Term_MEAN)) + geom_boxplot() + scale_x_discrete(labels=group_labels_fuzzy) + xlab("Plant community") + ylab("Long-term mean soil loss (ton/ha/yr)") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)), nrow = 1, labels = "auto")
State 1
(No name) State
TODO write state description
Community 1.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp1 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==1), species_cols = colnames(descriptive.df), user = user) # save the complete table for field crews write.csv(x = dom_spp_grp1, file = file.path(output_figure_folder, "Tables", paste0(target_ESG, "_1Juniper-pinyonwoodland_species.csv")), row.names = F) dom_spp_grp1_pub <- filter( select(dom_spp_grp1, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots), pct_plots >=10) %>% mutate(MeanCover = round(MeanCover, 2), pct_plots = round(pct_plots, 2)) %>% setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present")) format_tables_EDIT_style(data = dom_spp_grp1_pub, caption = "Community 1.1 plant community composition") %>% collapse_rows(columns = 1, valign = "top")
# TODO create pie charts of community composition by functional group. May need to change "production" to "percent cover" to fit with our data
# TODO create table of production (or percent cover) by functional type # knitr::kable(x = dom_plant_spp_tab, caption = "Table 5. Annual production by plant type")
# TODO create table of ground cover by type # knitr::kable(x = dom_plant_spp_tab, caption = "Table 6. Ground cover")
# TODO create table of func type percent cover at different heights # knitr::kable(x = dom_plant_spp_tab, caption = "Table 7. Canopy structure (% cover)")
State 2
(No name) State
Community 2.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp2 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==2), species_cols = colnames(descriptive.df), user = user) # save the complete table for field crews write.csv(x = dom_spp_grp2, file = file.path(output_figure_folder, "Tables", paste0(target_ESG, "_2Openjuniper-pinyonshrubland_species.csv")), row.names = F) dom_spp_grp2_pub <- filter( select(dom_spp_grp2, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots), pct_plots >=10) %>% mutate(MeanCover = round(MeanCover, 2), pct_plots = round(pct_plots, 2)) %>% setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present")) format_tables_EDIT_style(data = dom_spp_grp2_pub, caption = "Community 2.1 plant community composition") %>% collapse_rows(columns = 1, valign = "top")
State 3
(No name) State
TODO write state description
Community 3.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp3 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==3), species_cols = colnames(descriptive.df), user = user) # save the complete table for field crews write.csv(x = dom_spp_grp3, file = file.path(output_figure_folder, "Tables", paste0(target_ESG, "_3C4grassland_species.csv")), row.names = F) dom_spp_grp3_pub <- filter( select(dom_spp_grp3, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots), pct_plots >=10) %>% mutate(MeanCover = round(MeanCover, 2), pct_plots = round(pct_plots, 2)) %>% setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present")) format_tables_EDIT_style(data = dom_spp_grp3_pub, caption = "Community 3.1 plant community composition") %>% collapse_rows(columns = 1, valign = "top")
State 4
(No name) State
TODO write state description
Community 4.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp4 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==4), species_cols = colnames(descriptive.df), user = user) # save the complete table for field crews write.csv(x = dom_spp_grp4, file = file.path(output_figure_folder, "Tables", paste0(target_ESG, "_4Annualizedjunipersavanna_species.csv")), row.names = F) dom_spp_grp4_pub <- filter( select(dom_spp_grp4, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots), pct_plots >=10) %>% mutate(MeanCover = round(MeanCover, 2), pct_plots = round(pct_plots, 2)) %>% setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present")) format_tables_EDIT_style(data = dom_spp_grp4_pub, caption = "Community 4.1 plant community composition") %>% collapse_rows(columns = 1, valign = "top")
Additional community tables
# TODO create table of production (or percent cover) by species, with species linking to their USDA Plants profile # knitr::kable(x = dom_plant_spp_tab, caption = "Table 8. Community 1.1 plant community composition")
State 5
(No name) State
TODO write state description
Community 5.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp5 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==5), species_cols = colnames(descriptive.df), user = user) # save the complete table for field crews write.csv(x = dom_spp_grp5, file = file.path(output_figure_folder, "Tables", paste0(target_ESG, "_5Pinyon-juniperwoodland_species.csv")), row.names = F) dom_spp_grp5_pub <- filter( select(dom_spp_grp5, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots), pct_plots >=10) %>% mutate(MeanCover = round(MeanCover, 2), pct_plots = round(pct_plots, 2)) %>% setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present")) format_tables_EDIT_style(data = dom_spp_grp5_pub, caption = "Community 5.1 plant community composition") %>% collapse_rows(columns = 1, valign = "top")
State 6
(No name) State
TODO write state description
Community 6.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp6 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==6), species_cols = colnames(descriptive.df), user = user) # save the complete table for field crews write.csv(x = dom_spp_grp6, file = file.path(output_figure_folder, "Tables", paste0(target_ESG, "_6Gambeloakwoodland_species.csv")), row.names = F) dom_spp_grp6_pub <- filter( select(dom_spp_grp6, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots), pct_plots >=10) %>% mutate(MeanCover = round(MeanCover, 2), pct_plots = round(pct_plots, 2)) %>% setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present")) format_tables_EDIT_style(data = dom_spp_grp6_pub, caption = "Community 6.1 plant community composition") %>% collapse_rows(columns = 1, valign = "top")
Additional community tables
# TODO create table of production (or percent cover) by species, with species linking to their USDA Plants profile # knitr::kable(x = dom_plant_spp_tab, caption = "Table 8. Community 1.1 plant community composition")
State 7
(No name) State
TODO write state description
Community 7.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp7 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==7), species_cols = colnames(descriptive.df), user = user) # save the complete table for field crews write.csv(x = dom_spp_grp7, file = file.path(output_figure_folder, "Tables", paste0(target_ESG, "_7Sagebrushshrublandgrassland_species.csv")), row.names = F) dom_spp_grp7_pub <- filter( select(dom_spp_grp7, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots), pct_plots >=10) %>% mutate(MeanCover = round(MeanCover, 2), pct_plots = round(pct_plots, 2)) %>% setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present")) format_tables_EDIT_style(data = dom_spp_grp7_pub, caption = "Community 7.1 plant community composition") %>% collapse_rows(columns = 1, valign = "top")
Additional community tables
# TODO create table of production (or percent cover) by species, with species linking to their USDA Plants profile # knitr::kable(x = dom_plant_spp_tab, caption = "Table 8. Community 1.1 plant community composition")
State 8
(No name) State
TODO write state description
Community 8.1
(No name) State: (No name community)
TODO write state description
dom_spp_grp8 <- get_dominant_spp(data = filter(plot_data_descriptive, PlantCommunity_fuzzy==8), species_cols = colnames(descriptive.df), user = user) # save the complete table for field crews write.csv(x = dom_spp_grp8, file = file.path(output_figure_folder, "Tables", paste0(target_ESG, "_8Barejunipersavanna_species.csv")), row.names = F) dom_spp_grp8_pub <- filter( select(dom_spp_grp8, Indicator, ScientificName, CommonName, SpeciesCode, MeanCover, pct_plots), pct_plots >=10) %>% mutate(MeanCover = round(MeanCover, 2), pct_plots = round(pct_plots, 2)) %>% setNames(c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "Percent of plots where present")) format_tables_EDIT_style(data = dom_spp_grp8_pub, caption = "Community 8.1 plant community composition") %>% collapse_rows(columns = 1, valign = "top")
Additional community tables
# TODO create table of production (or percent cover) by species, with species linking to their USDA Plants profile # knitr::kable(x = dom_plant_spp_tab, caption = "Table 8. Community 1.1 plant community composition")
Interpretations
Animal community
TODO write wildlife narrative
Hydrological functions
TODO write hydrology narrative
Recreational uses
TODO write recreation narrative
Wood products
TODO write wood production narrative
Other products
TODO grazing narrative could go here
Supporting information
Other references
Herrick, J. E., Van Zee, J. W., McCord, S. E., Courtright, E. M., Karl, J. W., & Burkett, L. M. (2017). Monitoring Manual for Grassland, Shrubland, and Savanna Ecosystems (Second ed. Vol. 1: Core Methods). Las Cruces, New Mexico: USDA-ARS Jornada Experimental Range.
Maechler, M., Rousseeuw, P., Struyf, A., Hubert, M., Hornik, K.(2021). cluster: Cluster Analysis Basics and Extensions. R package version 2.1.2.
McCord, S. E., Brehm, J. R., Burnett, S. H., Dietrich, C., Edwards, B., Metz, L. J., Hernandez Narvaez, M., Pierson, F., Ramirez, K. S., Stauffer, N. G., Webb, N. P., and Tweedie, C. E. (2022). A framework and toolset for standardizing agroecosystem indicators. Ecological Indicators, 144. doi:10.1016/j.ecolind.2022.109511
Nauman, T. W., Burch, S. S., Humphries, J. T., Knight, A. C., & Duniway, M. C. (2022). A Quantitative Soil-Geomorphic Framework for Developing and Mapping Ecological Site Groups. Rangeland Ecology & Management, 81, 9-33. doi:10.1016/j.rama.2021.11.003
Oksanen, J., F. Guillaume Blanchet, Michael Friendly, Roeland Kindt, Pierre Legendre, Dan McGlinn, Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens, Eduard Szoecs and Helene Wagner (2020). vegan: Community Ecology Package. R package version 2.5-7. https://CRAN.R-project.org/package=vegan
Contributors
- Anna Knight
- Travis Nauman
- Samuel Burch
- Mike Duniway
Reference sheet
Interpreting Indicators of Rangeland Health is a qualitative assessment protocol used to determine ecosystem condition based on benchmark characteristics described in the Reference Sheet. A suite of 17 (or more) indicators are typically considered in an assessment. The ecological site(s) representative of an assessment location must be known prior to applying the protocol and must be verified based on soils and climate. Current plant community cannot be used to identify the ecological site.
# TODO create table with author, author contact info, date, approver, approved date, and "Composition (Indicators 10 and 12) based on" info for the IIRH reference sheet # knitr::kable(x = dom_plant_spp_tab)
Indicators
TODO fill in the 17 IIRH indicators for the reference state if available
TODO link to PDF of the reference sheet
Print options
TODO look into options for printing from the RMarkdown HTML format
# 1. plant functional group join plot with ordination # 2. lattice of abiotic gradients that distinguish the states # aero # rhem if available # soil stability # climate? # 3. NPP # 4. biodiversity # 5 remote sensing variance # 6. Functional type cover grouped by functional type and colored by state # 7. Disturbance and land use # MTBS # Grazing? # PJ treatments # SWreGAP protection classes #group_labels_short <- c("Grassland", "Perennial grass loss", "Shrubland", "Mesic shrubland", "Invaded") plot_data_descriptive$PlantCommunity_fuzzy_reordered <- factor(plot_data_descriptive$PlantCommunity_fuzzy, levels = c(as.character(1:8)), labels = c("Juniper-pinyon woodland", "Open juniper-pinyon shrubland", "C4 grassland", "Annualized juniper savanna", "Pinyon-juniper woodland", "Gambel oak woodland", "Sagebrush shrubland/grassland", "Bare juniper savanna"), ordered = T) pal_veg2_reordered <- pal_veg2 # 1. plant functional group join plot with ordination species_scores <- as.data.frame(scores(ord, choices = 1:3, display = "species")) # get axis scores for each response type plot_scores <- as.data.frame(scores(ord, choices = 1:3, display = "sites")) # get the axis scores for each site plot_scores$PlotCode <- rownames(plot_scores) plot_scores <- left_join(plot_scores, select(plot_data_descriptive, PlotCode, PlantCommunity_fuzzy, PlantCommunity_fuzzy_reordered)) # this bit gets the arrow length and direction info ord.fit <- envfit(ord ~ AH_PerenForbGrassCover + AH_AnnForbGrassCover + AH_ShrubCover + AH_TreeCover + TotalFoliarCover, data = plot_data_descriptive, na.rm=T, choices=1:3) vect_functypes <- as.data.frame(ord.fit$vectors$arrows*sqrt(ord.fit$vectors$r)) vect_functypes$var <- rownames(vect_functypes) vect_functypes$r <- ord.fit$vectors$r vect_functypes$varnames <- c("Perennial herbaceous", "Annual herbaceous", "Shrub", "Tree", "Total foliar") # # since I have 3 ordination axes, I'm filtering just to axis 1 and 2 for this plot vect_functypes_12 <- filter(vect_functypes, abs(NMDS3)<abs(NMDS1)|abs(NMDS3)<abs(NMDS2)) vect_functypes_32 <- filter(vect_functypes, abs(NMDS1)<abs(NMDS3)|abs(NMDS1)<abs(NMDS2)) vect_functypes_13 <- filter(vect_functypes, abs(NMDS2)<abs(NMDS1)|abs(NMDS2)<abs(NMDS3)) # make pretty plot ggplot! ord_pretty_plot <- function(points = plot_scores, vectors = vect_functypes, x_axis = "NMDS1", y_axis = "NMDS2", palette = pal_veg2_reordered, #group_labels = group_labels_short, labels = T){ ord_plot <- ggplot(data = points, aes_string(x=x_axis, y=y_axis)) + geom_point(data = points, aes_string(x=x_axis, y=y_axis, color="PlantCommunity_fuzzy_reordered"), size=.5) + stat_ellipse(data=points, aes_string(x=x_axis, y=y_axis, color="PlantCommunity_fuzzy_reordered"), level = 0.5, lwd = 1) + scale_color_manual(values =palette, name="Plant community"#, # set custom colors for ellipses #labels=group_labels ) + geom_segment(data = vectors, # make the arrows with environmental vars aes_string(x="0", xend=x_axis, y="0", yend=y_axis), arrow = arrow(length = unit(0.25, "cm")), lwd=1, colour="black", inherit.aes = F) + coord_fixed() + lims(x=c(-1.25,1.25), y=c(-1.25,1.25)) + theme_bw() + theme(axis.text = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), panel.grid.major = element_blank(), panel.grid.minor = element_blank()) if(labels){ ord_plot <- ord_plot + ggrepel::geom_label_repel(data = vectors, # put labels on the arrows aes_string(x=x_axis, y=y_axis, label="varnames"), size=3, colour="black") } return(ord_plot) } ord_plot12 <- ord_pretty_plot(x_axis = "NMDS1", y_axis = "NMDS2", vectors = vect_functypes_12) ord_plot32 <- ord_pretty_plot(x_axis = "NMDS3", y_axis = "NMDS2", vectors = vect_functypes_32) ord_plot13 <- ord_pretty_plot(x_axis = "NMDS1", y_axis = "NMDS3", vectors = vect_functypes_13) ord_plot_legend <- cowplot::get_legend(ord_plot12) ord_plot_all <- cowplot::plot_grid(ord_plot12 + theme(legend.position = "none"), ord_plot32 + theme(legend.position = "none"), ord_plot13 + theme(legend.position = "none"), ord_plot_legend, nrow = 2, labels = c("a", "b", "c")) today <- Sys.Date() highrestiff(plot_obj = ord_plot_all, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Ordination/SW_ShallowDeepRocky_NMDS_ellipses_functypes_", today, ".tif"), width_in = 6, height_in = 6, resolution_dpi = 300) ord_plot12_nolabels <- ord_pretty_plot(x_axis = "NMDS1", y_axis = "NMDS2", vectors = vect_functypes_12, labels = F) ord_plot32_nolabels <- ord_pretty_plot(x_axis = "NMDS3", y_axis = "NMDS2", vectors = vect_functypes_32, labels = F) ord_plot13_nolabels <- ord_pretty_plot(x_axis = "NMDS1", y_axis = "NMDS3", vectors = vect_functypes_13, labels = F) ord_plot_legend_nolabels <- cowplot::get_legend(ord_plot12_nolabels) ord_plot_all_nolabels <- cowplot::plot_grid(ord_plot12_nolabels + theme(legend.position = "none"), ord_plot32_nolabels + theme(legend.position = "none"), ord_plot13_nolabels + theme(legend.position = "none"), ord_plot_legend_nolabels, nrow = 2, labels = c("a", "b", "c")) today <- Sys.Date() highrestiff(plot_obj = ord_plot_all_nolabels, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Ordination/SW_ShallowDeepRocky_NMDS_ellipses_functypes_nolabels_", today, ".tif"), width_in = 6, height_in = 6, resolution_dpi = 300) # 2. lattice of abiotic gradients that distinguish the states violin_pretty_plot <- function(dat = plot_data_descriptive, y_var = "horizontal_flux_total_MN", y_lab = "Total horizontal dust flux", fill_var = NULL, fill_name = "Plant community", scaling = "area"){ if(!is.null(fill_var)){ violin_plot <- ggplot(data = dat, aes_string(x="PlantCommunity_fuzzy_reordered", y=y_var, fill = fill_var)) + scale_fill_manual(values =pal_veg2_reordered, name= fill_name, # set custom colors for ellipses #labels=group_labels_short ) }else{ violin_plot <- ggplot(data = dat, aes_string(x="PlantCommunity_fuzzy_reordered", y=y_var)) } violin_plot <- violin_plot + geom_violin(scale = scaling, trim = F, na.rm = T) + geom_boxplot(width=0.1) + #scale_x_discrete(labels=group_labels_short) + xlab(NULL) + ylab(y_lab) + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") return(violin_plot) } box_pretty_plot <- function(dat = plot_data_descriptive, y_var = "horizontal_flux_total_MN", y_lab = "Total horizontal dust flux", fill_var = "PlantCommunity_fuzzy_reordered", palette = pal_veg2_reordered){ group_labels <- paste0(as.character(levels(dat$PlantCommunity_fuzzy_reordered)), " (n = ", c(nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==1 & !is.na(dat[[y_var]]), ]), nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==2 & !is.na(dat[[y_var]]), ]), nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==3 & !is.na(dat[[y_var]]), ]), nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==4 & !is.na(dat[[y_var]]), ]), nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==5 & !is.na(dat[[y_var]]), ]), nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==6 & !is.na(dat[[y_var]]), ]), nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==7 & !is.na(dat[[y_var]]), ]), nrow(dat[as.integer(dat$PlantCommunity_fuzzy_reordered)==8 & !is.na(dat[[y_var]]), ]) ), ")") box_plot <- ggplot(data = dat, aes_string(x="PlantCommunity_fuzzy_reordered", y=y_var, fill = fill_var)) + geom_boxplot() + scale_x_discrete(labels=group_labels) + scale_fill_manual(values =palette, name="Plant community", # set custom colors for ellipses #labels=group_labels_short ) + xlab(NULL) + ylab(y_lab) + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") return(box_plot) } # aero aero_plot <- box_pretty_plot(dat = filter(plot_data_descriptive, !is.na(horizontal_flux_total_MN)), y_var = "horizontal_flux_total_MN", y_lab = "Horizontal dust flux (g/m*day)"#, #scaling = "width" ) + scale_y_log10(breaks = scales::trans_breaks("log10", function(x) 10^x), labels = scales::trans_format("log10", scales::label_math(10^.x))) #+ #annotation_logticks(sides = "l") aero_plot # soil stability soilstab_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "SoilStab_all", y_lab = "Soil stability") soilstab_plot # climate? AI_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "AI", y_lab = "Aridity index") AI_plot ppt_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "Ann_ppt_mm", y_lab = "Mean ann. ppt. (mm)") ppt_plot tmean_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "Ann_tmean", y_lab = "Mean ann. temp. (C)") tmean_plot # rhem rhem_plot_runoff <- box_pretty_plot(dat = plot_data_descriptive, y_var = "Runoff_Long_Term_MEAN", y_lab = "Mean runoff (mm/yr)") rhem_plot_runoff rhem_plot_soilloss <- box_pretty_plot(dat = plot_data_descriptive, y_var = "Soil_Loss_Long_Term_MEAN", y_lab = "Mean soil loss (tons/ha/yr)") rhem_plot_soilloss # put them all together abiotic_plots <- cowplot::plot_grid(aero_plot, soilstab_plot, rhem_plot_runoff, rhem_plot_soilloss, ppt_plot, tmean_plot, nrow = 3, #rel_heights = c(0.6, 1), align = "v", labels = "auto") abiotic_plots today <- Sys.Date() highrestiff(plot_obj = abiotic_plots, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_abiotic_boxplot_2_", today, ".tif"), width_in = 6, height_in = 9, resolution_dpi = 300) # put them in a different layout because we love drafts climate_plots <- cowplot::plot_grid(ppt_plot, tmean_plot, AI_plot, nrow = 1, labels = "auto") climate_plots highrestiff(plot_obj = climate_plots, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_climate_boxplot_", today, ".tif"), width_in = 9, height_in = 3.5, resolution_dpi = 300) soils_plots <- cowplot::plot_grid(aero_plot, soilstab_plot, rhem_plot_runoff, rhem_plot_soilloss, nrow = 2, labels = "auto") soils_plots highrestiff(plot_obj = soils_plots, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_soils_boxplot_", today, ".tif"), width_in = 6, height_in = 6, resolution_dpi = 300) # 3. NPP # NPP data is from the RAP 3.0 partitioned NPP data set (Robinson et al., 2019). I downloaded the predicted NPP values within a 35-m radius of the AIM, LMF, NRI, and NCPN plot centers using a Google Earth Engine script. If we continue using NPP for other ESGs, we will probably need to continue downloading point values for RAP products through GEE - the raster file sizes are enormous! # Robinson, N. P., Allred, B. W., Naugle, D. E., & Jones, M. O. (2019). Patterns of rangeland productivity and land ownership: Implications for conservation and management. Ecologial Applications, 23(3), e01862. NPP_files <- list.files(paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/RAP_do_not_sync/", target_ESG, "/NPP"), full.names = T) npp <- read.csv(NPP_files[1]) npp$Year <- as.numeric(substr(NPP_files[1], 171, 174)) for(file in NPP_files[-1]){ npp_temp <- read.csv(file) npp_temp$Year <- as.numeric(substr(file, 171, 174)) npp <- bind_rows(npp, npp_temp) rm(npp_temp) } npp <- select(npp, PlotCod, Year, afgNPP, pfgNPP, shrNPP, treNPP) %>% setNames(c("PlotCode", "Year", "afgNPP", "pfgNPP", "shrNPP", "treNPP")) plot_data_descriptive <- left_join(plot_data_descriptive, npp) afg_npp_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "afgNPP", y_lab = "Ann. herb. NPP (g C/m2)") afg_npp_plot pfg_npp_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "pfgNPP", y_lab = "Per. herb. NPP (g C/m2)") pfg_npp_plot shr_npp_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "shrNPP", y_lab = "Shrub NPP (g C/m2)") shr_npp_plot tre_npp_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "treNPP", y_lab = "Tree NPP (g C/m2)") tre_npp_plot pft_npp_plots <- cowplot::plot_grid(afg_npp_plot + xlab(NULL) + scale_x_discrete(labels=NULL), pfg_npp_plot + xlab(NULL) + scale_x_discrete(labels=NULL), shr_npp_plot, tre_npp_plot, nrow = 2, rel_heights = c(0.6, 1), align = "v", labels = "auto") pft_npp_plots today <- Sys.Date() highrestiff(plot_obj = pft_npp_plots, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_pftnpp_boxplot_", today, ".tif"), width_in = 6, height_in = 6.5, resolution_dpi = 300) # combining npp to herbaceous, woody, and total plot_data_descriptive$herbNPP <- plot_data_descriptive$afgNPP + plot_data_descriptive$pfgNPP plot_data_descriptive$woodyNPP <- plot_data_descriptive$shrNPP + plot_data_descriptive$treNPP plot_data_descriptive$totalNPP <- plot_data_descriptive$herbNPP + plot_data_descriptive$woodyNPP herb_npp_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "herbNPP", y_lab = "Herb. NPP (g C/m2)") herb_npp_plot woody_npp_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "woodyNPP", y_lab = "Woody NPP (g C/m2)") woody_npp_plot total_npp_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "totalNPP", y_lab = "Total NPP (g C/m2)") total_npp_plot gh_npp_plots <- cowplot::plot_grid(herb_npp_plot + xlab(NULL) + scale_x_discrete(labels=NULL), woody_npp_plot + xlab(NULL) + scale_x_discrete(labels=NULL), total_npp_plot, nrow = 3, rel_heights = c(0.6, 0.6, 1), align = "v", labels = "auto") gh_npp_plots today <- Sys.Date() highrestiff(plot_obj = gh_npp_plots, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_ghnpp_boxplot_", today, ".tif"), width_in = 3, height_in = 9, resolution_dpi = 300) highrestiff(plot_obj = total_npp_plot, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_totalnpp_boxplot_2_", today, ".tif"), width_in = 3, height_in = 3.5, resolution_dpi = 300) rm(npp) # 4. biodiversity diversity_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "shannon_diversity", y_lab = "Shannon diversity index") diversity_plot richness_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "sp_richness", y_lab = "Number of species") richness_plot biodiversity_plots <- cowplot::plot_grid(richness_plot, diversity_plot, total_npp_plot, #nrow = 2, rel_heights = c(0.6, 1), nrow = 1, labels = "auto") biodiversity_plots today <- Sys.Date() highrestiff(plot_obj = biodiversity_plots, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_biodiversity_boxplot_", today, ".tif"), #width_in = 3, width_in = 9, #height_in = 6.5, height_in = 3.5, resolution_dpi = 300) # 5. remote sensing variance # Ideally, we should probably do this from the LandCART total foliar cover predictions, but obtaining those rasters for the whole UCRB is very slow. For now, we'll use the sum of the RAP 3.0 plant functional type cover values for each plot to represent total foliar cover instead. The RAP values were downloaded from within a 35-m radius of each point in this ESG using a Google Earth Engine script. Consider switching to LandCART, since it actually predicts a total foliar cover value, when raster generation and downloading for the UCRB 1995-2019 (April 1 - Nov 1) is complete. library(broom) # for easy handling of many linear model parts library(purrr) # for functional programming, eep! library(tidyr) VegCover_files <- list.files(paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/RAP_do_not_sync/", target_ESG, "/VegCover"), full.names = T) VegCover <- read.csv(VegCover_files[1]) VegCover$Year <- as.numeric(substr(VegCover_files[1], 168, 171)) for(file in VegCover_files[-1]){ VegCover_temp <- read.csv(file) VegCover_temp$Year <- as.numeric(substr(file, 168, 171)) VegCover <- bind_rows(VegCover, VegCover_temp) rm(VegCover_temp) } VegCover <- select(VegCover, PlotCod, Bst_grp, Bst_gr_, Year, AFG, BGR, LTR, PFG, SHR, TRE) %>% setNames(c("PlotCode", "PlantCommunity_fuzzy", "PlantCommunity_name", "Year", "AFG", "BGR", "LTR", "PFG", "SHR", "TRE")) %>% rowwise() %>% mutate(TotalFoliar = sum(AFG, PFG, SHR, TRE, na.rm = T)) sampleyears <- select(plot_data_descriptive, PlotCode, Year) %>% setNames(c("PlotCode", "Year_sampled")) VegCover <- left_join(VegCover, sampleyears) VegCover_10yrs <- VegCover %>% # select the 10 years up to and including the field sampling year mutate(Years_to_sampling = Year_sampled - Year) %>% filter(Years_to_sampling >=0 & Years_to_sampling <10) VegCover_variance <- VegCover_10yrs %>% select(PlotCode, PlantCommunity_fuzzy, Year, Years_to_sampling, TotalFoliar) %>% as_tibble() %>% nest(data = c(Year, Years_to_sampling, TotalFoliar)) %>% # create an internal data frame for each plot mutate(model = map(data, ~lm(TotalFoliar ~ Years_to_sampling, data = .))) %>% # create a linear model from each plot's data frame mutate(data = map(model, augment), # get residuals for each model and add to internal plot data frames model_glance = map(model, glance)) %>% # get model summary statistics mutate(resid_variance = map_dbl(data, ~var(.$.resid)), # get the variance of the detrended data (aka residuals) resid_AC1 = map_dbl(data, ~acf(x = .$.resid, type = "correlation", lag.max = 1, plot = F)$acf[2]) # get the lag-1 autocorrelation of the detrended data ) %>% select(-data, -model) %>% unnest(model_glance) # make figures plot_data_descriptive <- left_join(plot_data_descriptive, select(VegCover_variance, PlotCode, resid_variance, resid_AC1)) variance_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "resid_variance", y_lab = "Variance") variance_plot AC1_plot <- box_pretty_plot(dat = plot_data_descriptive, y_var = "resid_AC1", y_lab = "Lag-1 autocorrelation") AC1_plot RSvariability_plots <- cowplot::plot_grid(variance_plot + xlab(NULL) + scale_x_discrete(labels=NULL), AC1_plot, nrow = 2, rel_heights = c(0.6, 1), labels = "auto") RSvariability_plots today <- Sys.Date() highrestiff(plot_obj = RSvariability_plots, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_RSvariability_boxplot_", today, ".tif"), width_in = 3, height_in = 6, resolution_dpi = 300) variance_violin <- violin_pretty_plot(dat = plot_data_descriptive, y_var = "resid_variance", y_lab = "Variance", fill_var = "PlantCommunity_fuzzy_reordered") variance_violin today <- Sys.Date() highrestiff(plot_obj = variance_violin, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Violin/SW_ShallowDeepRocky_variance_violin_", today, ".tif"), width_in = 3, height_in = 3.5, resolution_dpi = 300) rm(VegCover, VegCover_10yrs, VegCover_variance) # 6. Functional type cover grouped by functional type and colored by state functypes_grps <- c("TotalFoliarCover", #"AH_WoodyCover", #AH_ForbGrassCover, "AH_ShrubCover", "AH_TreeCover", "AH_PerenForbGrassCover", "AH_AnnForbGrassCover" ) functypes_cover_plot <- ggplot(data = select(plot_data_descriptive, any_of(c("PlantCommunity_fuzzy_reordered", functypes_grps))) %>% tidyr::pivot_longer(cols = all_of(functypes_grps), names_to = "CoverType", values_to = "PctCover"), aes(fill=PlantCommunity_fuzzy_reordered, y=PctCover, x=CoverType)) + geom_boxplot() + scale_fill_manual(values = pal_veg2_reordered, name = "Plant community"#, #labels = group_labels_short ) + scale_x_discrete(labels=c("Annual herbaceous", "Perennial herbaceous", "Shrub", "Tree", "Total foliar")) + xlab("Plant functional type") + ylab("Cover (%)") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), panel.grid.major = element_blank(), panel.grid.minor = element_blank()) functypes_cover_plot today <- Sys.Date() highrestiff(plot_obj = functypes_cover_plot, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Boxplot/SW_ShallowDeepRocky_pftcover_boxplot_", today, ".tif"), width_in = 6, height_in = 3.5, resolution_dpi = 300) # 7. Disturbance and land use # MTBS # Grazing? # PJ treatments # SWreGAP protection classes # MTBS library(sf) library(lubridate) mtbs <- st_read(dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/GIS/Disturbance_LandUse/mtbs_permis_DD_2020.shp") fuzzy_cluster_memb_sf <- st_as_sf(fuzzy_cluster_memb_sf, crs = st_crs(plot_locations)) %>% st_transform(., crs=st_crs(mtbs)) mtbs_plots <- st_join(x=fuzzy_cluster_memb_sf, y=mtbs) %>% select(PlotCode, Best_group, Best_group_name, DateSampled, Ig_Date, geometry) %>% mutate(Burned_presampling = ifelse(test = Ig_Date < DateSampled, yes = TRUE, no = FALSE)) %>% group_by(PlotCode) %>% summarise(Best_group = first(Best_group), Best_group_name = first(Best_group_name), DateSampled = first(DateSampled), Times_burned_presampling = length(which(Burned_presampling)), Times_burned_postsampling = length(which(!Burned_presampling)), Date_burned_1 = sort(Ig_Date)[1], Date_burned_2 = sort(Ig_Date)[2]) %>% st_drop_geometry() # write.csv(mtbs_plots, # "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Analyses/SW_ShallowDeepRocky_burnedplots.csv", # row.names = F) rm(mtbs) mtbs_plots$PlantCommunity_fuzzy_reordered <- factor(mtbs_plots$Best_group, levels = c(as.character(1:8)), labels = c(c("Juniper-pinyon woodland", "Open juniper-pinyon shrubland", "C4 grassland", "Annualized juniper savanna", "Pinyon-juniper woodland", "Gambel oak woodland", "Sagebrush shrubland/grassland", "Bare juniper savanna")), ordered = T) mtbs_summary <- mtbs_plots %>% group_by(PlantCommunity_fuzzy_reordered) %>% summarise(pct_Burned = (length(which(Times_burned_presampling > 0))/n())*100, pct_Burned_postsampling = (length(which(Times_burned_postsampling > 0))/n())*100) burn_plot <- ggplot(data = mtbs_summary, aes(x=PlantCommunity_fuzzy_reordered, y=pct_Burned, fill = PlantCommunity_fuzzy_reordered)) + geom_col() + scale_fill_manual(values = pal_veg2_reordered) + xlab(NULL) + ylab("% of plots burned pre-sampling") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") burn_plot burn_plot_post <- ggplot(data = mtbs_summary, aes(x=PlantCommunity_fuzzy_reordered, y=pct_Burned_postsampling, fill = PlantCommunity_fuzzy_reordered)) + geom_col() + scale_fill_manual(values = pal_veg2_reordered) + xlab(NULL) + ylab("% of plots burned post-sampling") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") burn_plot_post today <- Sys.Date() highrestiff(plot_obj = burn_plot, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_MTBSpctburnedpresample_boxplot_", today, ".tif"), width_in = 3, height_in = 3.5, resolution_dpi = 300) today <- Sys.Date() highrestiff(plot_obj = burn_plot_post, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_MTBSpctburnedpostsample_boxplot_", today, ".tif"), width_in = 3, height_in = 3.5, resolution_dpi = 300) # PJ Treatments WRI_polygons <- st_read(dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps/FickPJ", layer = "PJ") %>% st_transform(crs = st_crs(fuzzy_cluster_memb_sf)) sf_use_s2(FALSE) pj_plots <- st_join(fuzzy_cluster_memb_sf, WRI_polygons, sparse = T) %>% mutate(YearSampled = year(DateSampled), TreatedPreSampling = ifelse(test = YearSampled>=yearStart, yes = TRUE, no = FALSE)) %>% mutate(TreatedPreSampling = ifelse(test = is.na(TreatedPreSampling), yes = FALSE, no = TreatedPreSampling)) # Look at the types of treatments present, then categorize into broader types for visual summaries sort(unique(paste(pj_plots$method, pj_plots$treat, sep = " - "))) # additional documentation available in G:\Base Layers\Disturbance\TREATMENTS. # For WRI treatments, see also the database at https://wri.utah.gov/wri/project/search.html pj_plots[pj_plots$name %in% "Monument Ridge Bullhog", "method"] <- "bullhogfull size" # based on manual check of WRI website pj_plots[pj_plots$name %in% "Monument Ridge Bullhog", "treat"] <- "brush removal" pj_plots_cats <- pj_plots %>% filter(!is.na(src)) %>% mutate(GeneralTreatmenType = ifelse(test = treat %in% c(# my categories: prescribed fire, brush removal, herbaceous removal, seeding, soil improvements. CHECK aerator, mowing, and weeding for each ESG to make sure they still fit as assigned here!!! "aerator", "brush removal", "mowing", "weeding", "herbicide"), yes = "Woody removal", no = ifelse(test = treat %in% c("seeding"), yes = "Seeding only", no = NA)), TreatmentSubType = ifelse(test = GeneralTreatmenType=="Woody removal" & method %in% c("anchor chain", "bullhogfull size", "bullhogskid steer", "double drum (1-way)", "lop and scatter", "mastication", "mechanical", "mowing", "rollerchop", "thinning"), yes = "Mechanical", no = ifelse(test = GeneralTreatmenType=="Woody removal" & method %in% c("herbicide", "herbicide aerial"), yes = "Chemical", no = ifelse(test = GeneralTreatmenType=="Seeding only" & method %in% c("broadcast"), yes = "Broadcast seeding", no = "Unknown" )))) # consider just printing this as a table count(st_drop_geometry(pj_plots_cats), GeneralTreatmenType, TreatmentSubType, seeded, burned) n_plots_treated <- length(unique(filter(pj_plots_cats, TreatedPreSampling==T)$PlotCode)) n_plots_tx_twice_or_more <- pj_plots_cats %>% filter(TreatedPreSampling==T) %>% group_by(PlotCode) %>% summarise(n_diff_treatments = n_distinct(yearStart)) %>% filter(n_diff_treatments >= 2) %>% nrow() plots_per_group <- data.frame(Best_group = c(1:8), plots_per_group = c(nrow(filter(plot_data_first, PlantCommunity_fuzzy==1)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==2)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==3)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==4)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==5)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==6)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==7)), nrow(filter(plot_data_first, PlantCommunity_fuzzy==8)) )) pj_plot_pct <- pj_plots_cats %>% left_join(., plots_per_group) %>% filter(TreatedPreSampling==T) %>% st_drop_geometry() %>% group_by(PlotCode) %>% # next 3 lines select only the treatment closest to the sampling year arrange(desc(yearStart), .by_group = T) %>% mutate(Tx_flag = row_number()) %>% filter(Tx_flag==1) %>% ungroup() %>% group_by(Best_group) %>% summarise(MechanicalUnseeded = 100*length(which(TreatmentSubType=="Mechanical" & seeded==0))/first(plots_per_group), MechanicalSeeded = 100*length(which(TreatmentSubType=="Mechanical" & seeded==1))/first(plots_per_group), ChemicalUnseeded = 100*length(which(TreatmentSubType=="Chemical" & seeded==0))/first(plots_per_group), ChemicalSeeded = 100*length(which(TreatmentSubType=="Chemical" & seeded==1)/first(plots_per_group)) ) %>% tidyr::pivot_longer(cols = c(MechanicalUnseeded, MechanicalSeeded, ChemicalUnseeded, ChemicalSeeded), names_to = "TreatmentPair", values_to = "pct") %>% mutate(Best_group = as.factor(Best_group)) treatments_plot_pcts <- ggplot(data = pj_plot_pct, aes(x = Best_group, y = pct, fill = TreatmentPair)) + geom_col() + scale_fill_manual(values = c("#FDBF6F", "#FF7F00", "#CAB2D6", "#6A3D9A"), name = "Treatment", labels = c("Chemical woody treatment (seeded)", "Chemical woody treatment", "Mechanical woody treatment (seeded)", "Mechanical woody treatment")) + scale_x_discrete(name = NULL, labels = c("Open woodland", "Shrubland", "Invaded", "Grassland")) + ylab("% of plots treated in each plant community") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "right") treatments_plot_pcts today <- Sys.Date() highrestiff(plot_obj = treatments_plot_pcts, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_PJtreatments_pct_barplot_", today, ".tif"), width_in = 6, height_in = 3.5, resolution_dpi = 300) # GAP protection classes GAP_raster <- raster::raster("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/GIS/Disturbance_LandUse/PADUS_GAP_Status_Code1.tif") plot_GAP <- sf::st_as_sf(raster::extract(x = GAP_raster, y = fuzzy_cluster_memb_sf, sp = T)) %>% sf::st_drop_geometry() plot_data_descriptive <- left_join(plot_data_descriptive, select(plot_GAP, PlotCode, PADUS_GAP_Status_Code1)) GAP_proportions <- plot_data_descriptive %>% group_by(PlantCommunity_fuzzy_reordered) %>% # NOTE the PADUS codes are stored in the "GAP_Sts" attribute, but extraction pulls the raster value, # which is different! Checked attribute table in Arc to sort it out! summarise(pct_BiodivProt_NatDistPermitted = (length(which(PADUS_GAP_Status_Code1 == 4))/n())*100, pct_BiodivProt_NatDistSupressed = (length(which(PADUS_GAP_Status_Code1 == 3))/n())*100, pct_MultiUse_ExtractPermitted = (length(which(PADUS_GAP_Status_Code1 == 1))/n())*100, pct_NoBiodivProctect = (length(which(PADUS_GAP_Status_Code1 == 2))/n())*100, pct_PrivateLand = (length(which(is.na(PADUS_GAP_Status_Code1)))/n())*100) %>% tidyr::pivot_longer(data=., cols = c(pct_BiodivProt_NatDistPermitted, pct_BiodivProt_NatDistSupressed, pct_MultiUse_ExtractPermitted, pct_NoBiodivProctect, pct_PrivateLand), names_to = "ProtectionStatus", values_to = "Percent") protection_plot <- ggplot(data = GAP_proportions, aes(x=ProtectionStatus, y=Percent, fill=PlantCommunity_fuzzy_reordered)) + geom_col() + scale_fill_manual(values = pal_veg2_reordered, name = "Plant community" ) + scale_x_discrete(labels=c("1 - Biodiversity protected, natural disturbance permitted", "2 - Biodiversity protected, natural disturbance suppressed", "3 - Multi-use, extraction permitted", "4 - No biodiversity protection", "Private land")) + xlab("PADUS protection status") + ylab("% plots in plant community") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.margin = unit(c(2,2,2,5), "mm") ) protection_plot today <- Sys.Date() highrestiff(plot_obj = protection_plot, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_PADUSprotection_flipped_barplot_", today, ".tif"), width_in = 6, height_in = 3.5, resolution_dpi = 300) group_labels_fuzzy <- paste0(group_labels_base, " (n = ", c(nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==1)), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==2)), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==3)), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==4)), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==5)), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==6)), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==7)), nrow(filter(plot_data_descriptive, PlantCommunity_fuzzy==8)) ), ")") protection_plot2 <- ggplot(data = GAP_proportions, aes(fill=ProtectionStatus, y=Percent, x=PlantCommunity_fuzzy_reordered)) + geom_col() + scale_fill_manual(values = c("seagreen", "yellowgreen", "gray58", "dodgerblue", "palegoldenrod"), name = "PADUS protection status", labels = c("1 - Biodiversity protected, natural disturbance permitted", "2 - Biodiversity protected, natural disturbance suppressed", "3 - Multi-use, extraction permitted", "4 - No biodiversity protection", "Private land")) + scale_x_discrete(labels= group_labels_fuzzy) + xlab(NULL) + ylab("% plots in plant community") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.key.size = unit(1.5, 'lines')) protection_plot2 today <- Sys.Date() highrestiff(plot_obj = protection_plot2, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_PADUSprotection_dodge_barplot_", today, ".tif"), width_in = 6, height_in = 3.5, resolution_dpi = 300) GAP_proportions_byclass <- plot_data_descriptive %>% mutate(PADUS_GAP_Status_lumped = ifelse(test = PADUS_GAP_Status_Code1 %in% c(3,4), yes = "Protected public", no = NA)) %>% mutate(PADUS_GAP_Status_lumped = ifelse(test = PADUS_GAP_Status_Code1 %in% c(1,2), yes = "Unprotected public", no = PADUS_GAP_Status_lumped)) %>% mutate(PADUS_GAP_Status_lumped = ifelse(test = is.na(PADUS_GAP_Status_Code1), yes = "Private", no = PADUS_GAP_Status_lumped)) %>% group_by(PADUS_GAP_Status_lumped) %>% summarise(pct_JPwoodland = (length(which(PlantCommunity_fuzzy_reordered == "Juniper-pinyon woodland"))/n())*100, pct_OpenJPshrubland = (length(which(PlantCommunity_fuzzy_reordered == "Open juniper-pinyon shrubland"))/n())*100, pct_C4grassland = (length(which(PlantCommunity_fuzzy_reordered == "C4 grassland"))/n())*100, pct_AnnualizedJuniperSavanna = (length(which(PlantCommunity_fuzzy_reordered == "Annualized juniper savanna"))/n())*100, pct_PJwoodland = (length(which(PlantCommunity_fuzzy_reordered == "Pinyon-juniper woodland"))/n())*100, pct_GambelOakWoodland = (length(which(PlantCommunity_fuzzy_reordered == "Gambel oak woodland"))/n())*100, pct_SagebrushShrubland = (length(which(PlantCommunity_fuzzy_reordered == "Sagebrush shrubland/grassland"))/n())*100, pct_BareJuniperSavanna = (length(which(PlantCommunity_fuzzy_reordered == "Bare juniper savanna"))/n())*100, ) %>% tidyr::pivot_longer(data=., cols = c(pct_JPwoodland, pct_OpenJPshrubland, pct_C4grassland, pct_AnnualizedJuniperSavanna, pct_PJwoodland, pct_GambelOakWoodland, pct_SagebrushShrubland, pct_BareJuniperSavanna), names_to = "PlantCommunity", values_to = "Percent") %>% mutate(PlantCommunity_fuzzy_reordered = factor(PlantCommunity, levels = c("pct_JPwoodland", "pct_OpenJPshrubland", "pct_C4grassland", "pct_AnnualizedJuniperSavanna", "pct_PJwoodland", "pct_GambelOakWoodland", "pct_SagebrushShrubland", "pct_BareJuniperSavanna"), labels = c("Juniper-pinyon woodland", "Open juniper-pinyon shrubland", "C4 grassland", "Annualized juniper savanna", "Pinyon-juniper woodland", "Gambel oak woodland", "Sagebrush shrubland/grassland", "Bare juniper savanna"), ordered = T)) protection_lumped_plot <- ggplot(data = GAP_proportions_byclass, aes(x=PADUS_GAP_Status_lumped, y=Percent, fill=PlantCommunity_fuzzy_reordered)) + geom_col() + scale_fill_manual(values = pal_veg2_reordered, name = "Plant community" ) + #scale_x_discrete(labels=c("Protected - public", "Unprotected - public" "Private land")) + xlab(NULL) + ylab("% of plots in protection group") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), panel.grid.major = element_blank(), panel.grid.minor = element_blank()#, #plot.margin = unit(c(2,2,2,5), "mm") ) protection_lumped_plot disturbance_plots <- cowplot::plot_grid(burn_plot, protection_lumped_plot, nrow = 1, rel_widths = c(0.7, 1), labels = "auto") disturbance_plots today <- Sys.Date() highrestiff(plot_obj = disturbance_plots, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/Bar/SW_ShallowDeepRocky_disturbanceproctection_barplot_", today, ".tif"), width_in = 7, height_in = 3.5, resolution_dpi = 300) rm(GAP_raster, plot_GAP)
# Make materials for spring 2023 workshop handouts # 1. condensed community tables # 2. bar charts of functional types and ground cover - separate copy for each state ### community tables functype_table_palette <- c("#80B1D3", "#80B1D3", "#9ECAE1", "#FFFFB3", "#FFFFB3", "#FFFFB3", "#FFFFB3", "#C6DBEF", "#FB8072", "#BEBADA", "#8DD3C7") names(functype_table_palette) <- sort(unique(dom_spp_grp1$Indicator)) dom_spp_grp2_workshop <- dom_spp_grp2_pub %>% filter(`Mean cover (%)` >= 1) %>% arrange(desc(`Mean cover (%)`)) dom_spp_grp1_workshop <- dom_spp_grp1_pub %>% filter(`Mean cover (%)` >= 1) %>% arrange(desc(`Mean cover (%)`)) dom_spp_grp4_workshop <- dom_spp_grp4_pub %>% filter(`Mean cover (%)` >= 1) %>% arrange(desc(`Mean cover (%)`)) dom_spp_grp3_workshop <- dom_spp_grp3_pub %>% filter(`Mean cover (%)` >= 1) %>% arrange(desc(`Mean cover (%)`)) dom_spp_grp5_workshop <- dom_spp_grp5_pub %>% filter(`Mean cover (%)` >= 1) %>% arrange(desc(`Mean cover (%)`)) dom_spp_grp6_workshop <- dom_spp_grp6_pub %>% filter(`Mean cover (%)` >= 1) %>% arrange(desc(`Mean cover (%)`)) dom_spp_grp7_workshop <- dom_spp_grp7_pub %>% filter(`Mean cover (%)` >= 1) %>% arrange(desc(`Mean cover (%)`)) dom_spp_grp8_workshop <- dom_spp_grp8_pub %>% filter(`Mean cover (%)` >= 1) %>% arrange(desc(`Mean cover (%)`)) workingdat <- dom_spp_grp6_workshop kableExtra::kbl(x = workingdat, format = "html", row.names = F, col.names = c("Functional group", "Scientific name", "Common name", "USDA Plants code", "Mean cover (%)", "% of plots where present"), align = "c", caption = #"Open juniper-pinyon woodland" #grp1 #"Mixed shrubs with open juniper-pinyon" #grp2 #"C4 grassland with mixed shrubs and scattered juniper-pinyon" #grp3 #"Annualized invaded mixed shrubland/C3 grassland with scattered PJ" #grp4 #"Open pinyon-juniper woodland" #grp5 "Gambel oak woodland with mixed shrubs and grasses" #grp6 #"Sagebrush with mixed C3 grasses" #grp7 #"Juniper-pinyon with light C3 grass understory" #grp8 ) %>% kable_styling(bootstrap_options = c("bordered", "condensed")) %>% column_spec(kable_input=., column=4, width = "2cm") %>% column_spec(kable_input=., column=5, width = "2cm") %>% column_spec(kable_input=., column=6, width = "2.5cm") %>% row_spec(kable_input=., row=0, bold=T, background = "white") %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="C3 introduced perennial grasses"), background = functype_table_palette["C3 introduced perennial grasses"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="C3 native perennial grasses"), background = functype_table_palette["C3 native perennial grasses"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="C4 native perennial grasses"), background = functype_table_palette["C4 native perennial grasses"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="Introduced annual forbs"), background = functype_table_palette["Introduced annual forbs"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="Introduced annual grasses"), background = functype_table_palette["Introduced annual grasses"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="Introduced perennial forbs"), background = functype_table_palette["Introduced perennial forbs"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="Native annual forbs"), background = functype_table_palette["Native annual forbs"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="Native annual grasses"), background = functype_table_palette["Native annual grasses"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="Native perennial forbs"), background = functype_table_palette["Native perennial forbs"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="Shrubs"), background = functype_table_palette["Shrubs"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="Succulents"), background = functype_table_palette["Succulents"]) %>% row_spec(kable_input=., row = which(workingdat$`Functional group`=="Trees"), background = functype_table_palette["Trees"]) ### bar charts functypes_grps <- c("TotalFoliarCover", #"AH_WoodyCover", #AH_ForbGrassCover, "AH_ShrubCover", "AH_TreeCover", "AH_PerenForbGrassCover", "AH_AnnForbGrassCover" ) plot_data_descriptive_grp <- filter(plot_data_descriptive, PlantCommunity_fuzzy==8) functypes_cover_states <- ggplot(data = select(plot_data_descriptive_grp, any_of(c("PlantCommunity_fuzzy", functypes_grps))) %>% tidyr::pivot_longer(cols = all_of(functypes_grps), names_to = "CoverType", values_to = "PctCover"), aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) + geom_boxplot(outlier.shape = NA) + # make outliers invisible scale_fill_manual(values = c("#FFFFB3", "#80B1D3", "#FB8072", "#8DD3C7", "#666666"), name = "Functional type", labels = c("Annual herbaceous", "Perennial herbaceous", #"Woody", "Shrub", "Tree", "Total foliar cover") ) + scale_x_discrete(labels=group_labels_fuzzy[8]) + xlab(NULL) + ylab("Cover (%)") + ylim(c(0,100)) + ggtitle("Plant functional groups") + theme_bw() + theme(axis.text.x = element_text( color = "black"), axis.text.y = element_text(color = "black")) functypes_cover_states today <- Sys.Date() highrestiff(plot_obj = functypes_cover_states, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/JFSP_Workshop_1/PlantFuncTypes_", "8", "_", today, ".tif"), width_in = 6, height_in = 3.5, resolution_dpi = 300) ground_grps <- c("BareSoilCover", "FH_TotalLitterCover", "CA_percent_200plus", "FH_LichenMossCover", "CA_percent_100plus" ) plot_data_first_tall_grp <- filter(plot_data_first_tall, PlantCommunity_fuzzy==8) ground_cover_states <- ggplot(data = filter(plot_data_first_tall_grp, CoverType %in% ground_grps), aes(x=PlantCommunity_fuzzy, y=PctCover, fill=CoverType)) + geom_boxplot(outlier.shape = NA) + scale_fill_manual(values = c("#8DA0CB", "#E78AC3", "#A6D854", "#FFD92F", "#E5C494"), name = "Functional group", labels = names(cover_labels[which(cover_labels %in% ground_grps)]) ) + scale_x_discrete(labels=group_labels_fuzzy[8]) + xlab(NULL) + ylab("Cover (%)") + ylim(c(0,100)) + ggtitle("Ground cover") + theme_bw() + theme(axis.text.x = element_text(color = "black"), axis.text.y = element_text(color = "black")) ground_cover_states today <- Sys.Date() highrestiff(plot_obj = ground_cover_states, file = paste0("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/", target_ESG, "/JFSP_Workshop_1/GroundCover_", "8", "_", today, ".tif"), width_in = 6, height_in = 3.5, resolution_dpi = 300)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.