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_SandyUplands"
user <- "Anna"
#user <- "VPN"
library(dplyr)
library(kableExtra)
library(STMdevelopment)
# 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 deeper sandy soils without significant 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 with mostly gentle slopes on a variety of landforms (n=r length(landforms_unique)). Some areas can have steeper slopes, but these are not typical of the group.
# Slope
slope_ave <- mean(comps@site$slope_r.x)
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 are moderately deep or deeper to bedrock and are composed primarily of alluvium and eolian sediments. Surface and subsurface horizons very high sand content (>75% on average) with textures including loamy sand, sand, and sandy loam. Soils are non-saline and non sodic and can have up to 10% calcium carbonate, but generally have less than 5% carbonates. 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))). Water erosion hazard is moderate and wind erosion hazard is severe.
## 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 sites have high risk for wind erosion and dune mobilization, especially with drought or large disturbances. They also have potential for water erosion issues. Compared to reference states, alternative states may have higher bare ground, annual invasion, woody encroachment, and/or loss of perennial grasses.
State and transition model development
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 <- c("CA_percent_100plus", "CA_percent_200plus")
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 <- quantile(plot_data_first$UCRB_SGUs_ProbMax, probs = 0.1)
plot_data_first <- filter(plot_data_first, UCRB_SGUs_ProbMax > certainty_cutoff)
# 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")))
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 10th 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.7, # 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 <- 4 #2 it really wants to be 2!!
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 <- as.factor(grp_fuzz)
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("Mixed sagebrush, native per.
# grasses, & cheatgrass (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==1)), ")"),
# paste0("Open pinyon-juniper woodland
# (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==2)), ")"),
# paste0("Mixed pinyon-juniper &
# sagebrush (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==3)), ")"),
# paste0("Native perennial grassland
# (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==4)), ")"))
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
colnames(fuzzy_cluster_membership) <- c(paste0("Group_", 1:4, "_membership"), "PlotCode", "Best_group")
group_names <- data.frame(Best_group = 1:5,
Best_group_name = c("C4 & C3 perennial grassland",
"Depauperate grassland",
"Sagebrush & C3 perennial grass shrubland",
"Gambel oak shrubland",
"Non-native annual grassland"))
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)
# sf::st_write(fuzzy_cluster_memb_sf,
# dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps",
# layer = paste0("SW_LoamyUplands_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_LoamyUplands_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",
# layer = paste0("SW_LoamyUplands_fuzzyclust_mike_", Sys.Date(), ".shp"),
# driver = "ESRI Shapefile")
State and transition model {.tabset}
CUSTOM DIAGRAM
# TODO create custom diagram of states
knitr::include_graphics("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Analyses/Semiarid_Warm_SandyUplands_boxandarrow_DRAFT.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"),
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
) %>%
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"
# summarize species occurrance
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
) %>%
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, 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 %>%
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)
# # 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,
ATCA2,
ERNA10,
CHVI8,
EPVI,
PUTR2,
GUSA2,
JUOS,
PIED,
QUGA
) %>%
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,
CHVI8,
EPVI,
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 <- arrange(filter(indicator_descriptions, Indicator_code %in% cover_funcgrps), Indicator_code)$Indicator
names(cover_funcgrps) <- funcgrps_names
cover_spp <- sort(c("CHVI8",
"EPVI",
"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",
"CHVI8",
"EPVI",
"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))
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_discrete(name="Plant community", labels=group_labels_fuzzy) +
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"))
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))
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")
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
State 1
Grassland State
TODO write state description
Community 1.1
Grassland State: native perennial grassland
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)
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 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
Woodland State
Community 2.1
Woodland State: open pinyon-juniper woodland
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)
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
Grass and Shrubland State
TODO write state description
Community 3.1
Grass and Shrubland State: mixed sagebrush, native perennial grasses, and cheatgrass
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)
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 3.1 plant community composition") %>%
collapse_rows(columns = 1, valign = "top")
State 4
Shrubby Woodland State
TODO write state description
Community 4.1
Shrubby Woodland State: mixed pinyon-juniper and sagebrush
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)
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 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")
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 (probably "None")
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
plot_data_descriptive$PlantCommunity_fuzzy_reordered <- factor(plot_data_descriptive$PlantCommunity_fuzzy,
levels = c("1", "3", "4", "2", "5"),
labels = c("Grassland", "Shrubland", "Mesic shrubland",
"Perennial grass loss", "Invaded"),
ordered = T)
pal_veg2_reordered <- c(pal_veg2[1], pal_veg2[3], pal_veg2[4], pal_veg2[2], pal_veg2[5])
# 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"))
highrestiff(plot_obj = ord_plot_all,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Ordination/SW_LoamyUplands_NMDS_ellipses_functypes.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"))
highrestiff(plot_obj = ord_plot_all_nolabels,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Ordination/SW_LoamyUplands_NMDS_ellipses_functypes_nolabels.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){
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_short) +
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 if available
# put them all together
abiotic_plots <- cowplot::plot_grid(aero_plot + xlab(NULL) + scale_x_discrete(labels=NULL), soilstab_plot + xlab(NULL) + scale_x_discrete(labels=NULL),
ppt_plot, tmean_plot,
nrow = 2, rel_heights = c(0.6, 1), align = "v",
labels = "auto")
abiotic_plots
highrestiff(plot_obj = abiotic_plots,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_abiotic_boxplot_2.tif",
width_in = 6,
height_in = 6,
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 = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_climate_boxplot.tif",
width_in = 9,
height_in = 3.5,
resolution_dpi = 300)
soils_plots <- cowplot::plot_grid(aero_plot, soilstab_plot,
nrow = 1,
labels = "auto")
soils_plots
highrestiff(plot_obj = soils_plots,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_soils_boxplot.tif",
width_in = 6,
height_in = 3.5,
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("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/RAP_SWLoamyUplands/NPP",
full.names = T)
npp <- read.csv(NPP_files[1])
npp$Year <- as.numeric(substr(NPP_files[1], 124, 127))
for(file in NPP_files[-1]){
npp_temp <- read.csv(file)
npp_temp$Year <- as.numeric(substr(file, 124, 127))
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
highrestiff(plot_obj = pft_npp_plots,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_pftnpp_boxplot.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
highrestiff(plot_obj = gh_npp_plots,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_ghnpp_boxplot.tif",
width_in = 3,
height_in = 9,
resolution_dpi = 300)
highrestiff(plot_obj = total_npp_plot,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_totalnpp_boxplot_2.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
highrestiff(plot_obj = biodiversity_plots,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_biodiversity_boxplot_2.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("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/RAP_SWLoamyUplands/VegCover",
full.names = T)
VegCover <- read.csv(VegCover_files[1])
VegCover$Year <- as.numeric(substr(VegCover_files[1], 134, 137))
for(file in VegCover_files[-1]){
VegCover_temp <- read.csv(file)
VegCover_temp$Year <- as.numeric(substr(file, 134, 137))
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)
# do we care that most of the linear models are not significant and only about half have an R-squared > 0.1?
# 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
highrestiff(plot_obj = RSvariability_plots,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_RSvariability_boxplot.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
highrestiff(plot_obj = variance_violin,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Violin/SW_LoamyUplands_variance_violin.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
highrestiff(plot_obj = functypes_cover_plot,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_pftcover_boxplot_2.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)),
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_LoamyUplands_burnedplots.csv",
row.names = F)
rm(mtbs)
mtbs_plots$PlantCommunity_fuzzy_reordered <- factor(mtbs_plots$Best_group,
levels = c("1", "3", "4", "2", "5"),
labels = c("Grassland", "Shrubland", "Mesic shrubland",
"Perennial grass loss", "Invaded"),
ordered = T)
mtbs_summary <- mtbs_plots %>%
group_by(PlantCommunity_fuzzy_reordered) %>%
summarise(pct_Burned = (length(which(Times_burned_presampling > 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") +
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
# 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) %>%
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
highrestiff(plot_obj = protection_plot,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/SW_LoamyUplands_PADUSprotection_flipped_barplot.tif",
width_in = 6,
height_in = 3.5,
resolution_dpi = 300)
protection_plot2 <- ggplot(data = GAP_proportions,
aes(fill=ProtectionStatus, y=Percent, x=PlantCommunity_fuzzy)) +
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_short) +
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
highrestiff(plot_obj = protection_plot2,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/SW_LoamyUplands_PADUSprotection_dodge_barplot.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_Grassland = (length(which(PlantCommunity_fuzzy_reordered == "Grassland"))/n())*100,
pct_Shrubland = (length(which(PlantCommunity_fuzzy_reordered == "Shrubland"))/n())*100,
pct_MesicShrubland = (length(which(PlantCommunity_fuzzy_reordered == "Mesic shrubland"))/n())*100,
pct_GrassLoss = (length(which(PlantCommunity_fuzzy_reordered == "Perennial grass loss"))/n())*100,
pct_Invaded = (length(which(PlantCommunity_fuzzy_reordered == "Invaded"))/n())*100) %>%
pivot_longer(data=.,
cols = c(pct_Grassland,
pct_Shrubland,
pct_MesicShrubland,
pct_GrassLoss,
pct_Invaded),
names_to = "PlantCommunity",
values_to = "Percent") %>%
mutate(PlantCommunity_fuzzy_reordered = factor(PlantCommunity,
levels = c("pct_Grassland", "pct_Shrubland", "pct_MesicShrubland",
"pct_GrassLoss", "pct_Invaded"),
labels = c("Grassland", "Shrubland", "Mesic shrubland",
"Perennial grass loss", "Invaded"),
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
highrestiff(plot_obj = disturbance_plots,
file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/SW_LoamyUplands_disturbanceproctection_barplot.tif",
width_in = 7,
height_in = 3.5,
resolution_dpi = 300)
rm(GAP_raster, plot_GAP)
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_SandyUplands" user <- "Anna" #user <- "VPN" library(dplyr) library(kableExtra) library(STMdevelopment)
# 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 deeper sandy soils without significant 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 with mostly gentle slopes on a variety of landforms (n=r length(landforms_unique)). Some areas can have steeper slopes, but these are not typical of the group.
# Slope slope_ave <- mean(comps@site$slope_r.x) 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 are moderately deep or deeper to bedrock and are composed primarily of alluvium and eolian sediments. Surface and subsurface horizons very high sand content (>75% on average) with textures including loamy sand, sand, and sandy loam. Soils are non-saline and non sodic and can have up to 10% calcium carbonate, but generally have less than 5% carbonates. 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))). Water erosion hazard is moderate and wind erosion hazard is severe.
## 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 sites have high risk for wind erosion and dune mobilization, especially with drought or large disturbances. They also have potential for water erosion issues. Compared to reference states, alternative states may have higher bare ground, annual invasion, woody encroachment, and/or loss of perennial grasses.
State and transition model development
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 <- c("CA_percent_100plus", "CA_percent_200plus") 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 <- quantile(plot_data_first$UCRB_SGUs_ProbMax, probs = 0.1) plot_data_first <- filter(plot_data_first, UCRB_SGUs_ProbMax > certainty_cutoff) # 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")))
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 10th 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.7, # 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 <- 4 #2 it really wants to be 2!! 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 <- as.factor(grp_fuzz)
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("Mixed sagebrush, native per. # grasses, & cheatgrass (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==1)), ")"), # paste0("Open pinyon-juniper woodland # (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==2)), ")"), # paste0("Mixed pinyon-juniper & # sagebrush (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==3)), ")"), # paste0("Native perennial grassland # (n = ", nrow(filter(plot_data_first, PlantCommunity_fuzzy==4)), ")")) 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 colnames(fuzzy_cluster_membership) <- c(paste0("Group_", 1:4, "_membership"), "PlotCode", "Best_group") group_names <- data.frame(Best_group = 1:5, Best_group_name = c("C4 & C3 perennial grassland", "Depauperate grassland", "Sagebrush & C3 perennial grass shrubland", "Gambel oak shrubland", "Non-native annual grassland")) 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) # sf::st_write(fuzzy_cluster_memb_sf, # dsn = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/Maps", # layer = paste0("SW_LoamyUplands_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_LoamyUplands_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", # layer = paste0("SW_LoamyUplands_fuzzyclust_mike_", Sys.Date(), ".shp"), # driver = "ESRI Shapefile")
State and transition model {.tabset}
CUSTOM DIAGRAM
# TODO create custom diagram of states knitr::include_graphics("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Analyses/Semiarid_Warm_SandyUplands_boxandarrow_DRAFT.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"), 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 ) %>% 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" # summarize species occurrance 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 ) %>% 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, 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 %>% 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)
# # 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, ATCA2, ERNA10, CHVI8, EPVI, PUTR2, GUSA2, JUOS, PIED, QUGA ) %>% 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, CHVI8, EPVI, 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 <- arrange(filter(indicator_descriptions, Indicator_code %in% cover_funcgrps), Indicator_code)$Indicator names(cover_funcgrps) <- funcgrps_names cover_spp <- sort(c("CHVI8", "EPVI", "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", "CHVI8", "EPVI", "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))
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_discrete(name="Plant community", labels=group_labels_fuzzy) + 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"))
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))
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") 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
State 1
Grassland State
TODO write state description
Community 1.1
Grassland State: native perennial grassland
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) 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 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
Woodland State
Community 2.1
Woodland State: open pinyon-juniper woodland
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) 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
Grass and Shrubland State
TODO write state description
Community 3.1
Grass and Shrubland State: mixed sagebrush, native perennial grasses, and cheatgrass
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) 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 3.1 plant community composition") %>% collapse_rows(columns = 1, valign = "top")
State 4
Shrubby Woodland State
TODO write state description
Community 4.1
Shrubby Woodland State: mixed pinyon-juniper and sagebrush
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) 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 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")
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 (probably "None")
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 plot_data_descriptive$PlantCommunity_fuzzy_reordered <- factor(plot_data_descriptive$PlantCommunity_fuzzy, levels = c("1", "3", "4", "2", "5"), labels = c("Grassland", "Shrubland", "Mesic shrubland", "Perennial grass loss", "Invaded"), ordered = T) pal_veg2_reordered <- c(pal_veg2[1], pal_veg2[3], pal_veg2[4], pal_veg2[2], pal_veg2[5]) # 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")) highrestiff(plot_obj = ord_plot_all, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Ordination/SW_LoamyUplands_NMDS_ellipses_functypes.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")) highrestiff(plot_obj = ord_plot_all_nolabels, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Ordination/SW_LoamyUplands_NMDS_ellipses_functypes_nolabels.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){ 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_short) + 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 if available # put them all together abiotic_plots <- cowplot::plot_grid(aero_plot + xlab(NULL) + scale_x_discrete(labels=NULL), soilstab_plot + xlab(NULL) + scale_x_discrete(labels=NULL), ppt_plot, tmean_plot, nrow = 2, rel_heights = c(0.6, 1), align = "v", labels = "auto") abiotic_plots highrestiff(plot_obj = abiotic_plots, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_abiotic_boxplot_2.tif", width_in = 6, height_in = 6, 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 = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_climate_boxplot.tif", width_in = 9, height_in = 3.5, resolution_dpi = 300) soils_plots <- cowplot::plot_grid(aero_plot, soilstab_plot, nrow = 1, labels = "auto") soils_plots highrestiff(plot_obj = soils_plots, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_soils_boxplot.tif", width_in = 6, height_in = 3.5, 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("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/RAP_SWLoamyUplands/NPP", full.names = T) npp <- read.csv(NPP_files[1]) npp$Year <- as.numeric(substr(NPP_files[1], 124, 127)) for(file in NPP_files[-1]){ npp_temp <- read.csv(file) npp_temp$Year <- as.numeric(substr(file, 124, 127)) 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 highrestiff(plot_obj = pft_npp_plots, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_pftnpp_boxplot.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 highrestiff(plot_obj = gh_npp_plots, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_ghnpp_boxplot.tif", width_in = 3, height_in = 9, resolution_dpi = 300) highrestiff(plot_obj = total_npp_plot, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_totalnpp_boxplot_2.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 highrestiff(plot_obj = biodiversity_plots, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_biodiversity_boxplot_2.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("C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/RAP_SWLoamyUplands/VegCover", full.names = T) VegCover <- read.csv(VegCover_files[1]) VegCover$Year <- as.numeric(substr(VegCover_files[1], 134, 137)) for(file in VegCover_files[-1]){ VegCover_temp <- read.csv(file) VegCover_temp$Year <- as.numeric(substr(file, 134, 137)) 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) # do we care that most of the linear models are not significant and only about half have an R-squared > 0.1? # 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 highrestiff(plot_obj = RSvariability_plots, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_RSvariability_boxplot.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 highrestiff(plot_obj = variance_violin, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Violin/SW_LoamyUplands_variance_violin.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 highrestiff(plot_obj = functypes_cover_plot, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/Boxplot/SW_LoamyUplands_pftcover_boxplot_2.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)), 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_LoamyUplands_burnedplots.csv", row.names = F) rm(mtbs) mtbs_plots$PlantCommunity_fuzzy_reordered <- factor(mtbs_plots$Best_group, levels = c("1", "3", "4", "2", "5"), labels = c("Grassland", "Shrubland", "Mesic shrubland", "Perennial grass loss", "Invaded"), ordered = T) mtbs_summary <- mtbs_plots %>% group_by(PlantCommunity_fuzzy_reordered) %>% summarise(pct_Burned = (length(which(Times_burned_presampling > 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") + 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 # 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) %>% 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 highrestiff(plot_obj = protection_plot, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/SW_LoamyUplands_PADUSprotection_flipped_barplot.tif", width_in = 6, height_in = 3.5, resolution_dpi = 300) protection_plot2 <- ggplot(data = GAP_proportions, aes(fill=ProtectionStatus, y=Percent, x=PlantCommunity_fuzzy)) + 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_short) + 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 highrestiff(plot_obj = protection_plot2, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/SW_LoamyUplands_PADUSprotection_dodge_barplot.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_Grassland = (length(which(PlantCommunity_fuzzy_reordered == "Grassland"))/n())*100, pct_Shrubland = (length(which(PlantCommunity_fuzzy_reordered == "Shrubland"))/n())*100, pct_MesicShrubland = (length(which(PlantCommunity_fuzzy_reordered == "Mesic shrubland"))/n())*100, pct_GrassLoss = (length(which(PlantCommunity_fuzzy_reordered == "Perennial grass loss"))/n())*100, pct_Invaded = (length(which(PlantCommunity_fuzzy_reordered == "Invaded"))/n())*100) %>% pivot_longer(data=., cols = c(pct_Grassland, pct_Shrubland, pct_MesicShrubland, pct_GrassLoss, pct_Invaded), names_to = "PlantCommunity", values_to = "Percent") %>% mutate(PlantCommunity_fuzzy_reordered = factor(PlantCommunity, levels = c("pct_Grassland", "pct_Shrubland", "pct_MesicShrubland", "pct_GrassLoss", "pct_Invaded"), labels = c("Grassland", "Shrubland", "Mesic shrubland", "Perennial grass loss", "Invaded"), 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 highrestiff(plot_obj = disturbance_plots, file = "C:/Users/aknight/Documents/Telework_Backups/V_drive/ANNA_KNIGHT/ESG/STM/Figures/SW_LoamyUplands_disturbanceproctection_barplot.tif", width_in = 7, height_in = 3.5, resolution_dpi = 300) rm(GAP_raster, plot_GAP)
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.