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 report is provisional and preliminary, and the information is subject to change without notification. The findings have not been peer-reviewed nor Bureau approved.

General information

knitr::opts_chunk$set(echo = TRUE)
target_ESG <- "Semiarid_Warm_SandyUplands_LoamyUplands"
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)

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 ny Nauman et al., (In Review).

esds <- esd_data_pull(user)
comps <- comp_data_pull(user)
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.

## Associated ecological sites
format_tables_EDIT_style(data = esds[,c("ecoclassid","ecoclassname")],
                         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)

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(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 horizons have sand, loamy sand, sandy loam, and loam textures and subsurface horizons are similar but can also include sandy clay loams. Soils are non-saline and non sodic and can have up to 20% calcium carbonate, but generally have less than 6% carbonates. Soil pH ranges from 5.0 to 9.6, but is generally closer to 8.0. 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:

r print(unique(comps@site$compname))

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

TODO narrative about key types of transitions and disturbances

State and transition model development

# TODO pull in table of Ecological Site STMs for this ESG
apriori_stms <- read.csv(data_file_paths(user)$apriori_stms, stringsAsFactors = F,
                         na.strings = c("NA", "", " "), header = F)

apriori_stms_target <- dplyr::filter(apriori_stms, V1 %in% c(NA, "ESG", target_ESG))
apriori_stms_target <-  t(apriori_stms_target)[4:41,]
colnames(apriori_stms_target) <- c("State", "Plant community", "Proportion of ESDs containing this state")
# apriori_stms_target <- apriori_stms_target[c(1,2,41,3:40), ] %>%
#   as.data.frame() %>%
#   mutate(V3 = stringr::str_replace_all(V3, Climate_group_definitions_replace)) %>%
#   mutate(V3 = stringr::str_replace_all(V3, SGU_definitions_replace))

format_tables_EDIT_style(data = apriori_stms_target,
                         caption = "STM summary from Ecological Sites within this ESG") %>%
  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","Ponderosa","Aspen")]
par(mar = c(4,8,4,2)) ## This works when knited
boxplot(esds[,resp_vars], horizontal = TRUE,las = 1, main = "Reference Production (lbs/Acre)")

indicators <- c(#"AH_C4PerenGrassCover", # C4 native perennial grasses TODO calculate C4 NATIVE
                #"AH_C3PerenGrassCover", # C3 native perennial grasses TODO calculate C3 NATIVE
                "AH_C3NativePerenGrassCover", 
                "AH_C3IntroducedPerenGrassCover",
                "AH_C4NativePerenGrassCover", 
                "AH_C4IntroducedPerenGrassCover",
                "AH_IntroducedPerenGrassCover", # Non-native perennial grasses
                "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
                "BareSoilCover", # Bare soil
                "CP_percent_100to200", # Canopy gaps > 100 cm TODO do we want annual or perennial gaps?
                "CP_percent_200plus",
                "FH_LichenCover", # Lichen + moss combined cover
                "FH_MossCover")

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",
                             "IM_NCPN",
                             #"LMF",
                             #"NRI",
                             "Parashant",
                             #"TerrADat", 
                             "VanScoyocThesis"
                           )

plot_data <- plot_data_pull(user=user,
                             target_ESG = target_ESG,
                           data_sources = data_sources,
                           indicators = indicators,
                           shrub_by_spp = shrub_by_spp,
                           subshrub_by_spp = subshrub_by_spp,
                           tree_by_spp = tree_by_spp,
                           opuntia_combined = opuntia_combined
                           )

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)

kbl(x = indicator_descriptions[,-1], format = "html", row.names = F, align = "c",
    caption = "Indicators used in ordination and clustering analysis to develop plant communities") %>%
  kable_styling(bootstrap_options = "bordered") %>%
  row_spec(kable_input=., row = 0:nrow(indicator_descriptions), background = "lightsteelblue") %>%
  collapse_rows(columns = 1, valign = "top")

# remove incomplete rows - Bray Curtis distance can't work with NAs
plot_data_clean <- na.omit(plot_data)

# include only the first year that a plot has complete data
plot_data_first <- plot_data_clean %>%
  dplyr::group_by(PlotCode) %>%
  dplyr::arrange(.data=., Year, .by_group = TRUE) %>%
  dplyr::filter(dplyr::row_number()==1) %>%
  dplyr::ungroup()

# TODO update I&M species codes to match current USDA Plants codes! probably a step for the plot_data_pull function?

# make the ordination data frame
ord.df <- dplyr::select(plot_data_first, -SourceKey, -PlotID, -SiteName, -PlotName,
                        -Year, -Longitude_NAD83, -Latitude_NAD83,
                        -Month, -Day) %>%
  tibble::column_to_rownames("PlotCode") # keep plot code as row name

# Keep only species that occur in at least 2 plots
ord.df <- ord.df[ ,which(colSums(ord.df[,1:ncol(ord.df)] !=0) >= 2)] 

Field data from were collected by field crews with the US Geological Survey and the NPS Inventory and Monitoring Northern Colorado Plateau Network between r min(plot_data_first$Year) and r max(plot_data_first$Year). r nrow(ord.df) plots were included 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 (in development, McCord & Stauffer, 2020). When plots were sampled in multiple years, only the first sampling year was used in STM development.

The R packages “cluster” and “dendextend” were used to perform a hierarchical cluster analysis. The distance matrix was calculated using Bray-Curtis distance and a flexible beta linkage was used for clustering (β = –0.25). To determine the optimal number of groups, we used the “labdsv” R package to calculate indicator species values (IndVal, Dufrene and Legendre 1997) when the cluster analysis results were split into 2 to 12 groups. We then selected the number of groups that minimized mean p-values and maximized the number of significant indicator species (Figure \@ref(fig:clustering)).

The clustering results were visualized in ordination spaced, using non-metric multidimensional scaling (NMDS) with Bray-Curtis distance ("vegan" package in R).

# cluster and indicator species analysis
library(cluster) # for cluster analysis
library(dendextend)
library(labdsv)
library(foreach)
library(vegan)

# initial clustering
dist.df <- vegan::vegdist(ord.df, method = "bray") # creates Bray-Curtis distance matrix
#clust <- hclust(dist.df, method = "average") # use agnes() in cluster package for flexible beta method
clust1 <- agnes(x= dist.df, diss=T, method="flexible", par.method = c(0.625,0.625,-0.25,0)) # flexible beta with beta=-0.25
#plot(clust1)
hclust1 <- as.hclust(clust1) # covert to class hclust for easier plotting
#plot(hclust1, labels = FALSE)

# indicator species
# use IndVal from Dufrene and Legendre (1997)
# groupings from hierarchical clustering (grp)
# ALSO CONSIDER USING MRPP

# create vectors of groupings with different numbers of clusters
grps <- data.frame(site=1:nrow(ord.df))
for (i in 2:12) {
  grps <- cbind(grps, cutree(hclust1, k=i))
}
colnames(grps) <- c("site", paste(2:12, "groups", sep="_"))

# calculate IndVal for each species in each number of clusters
listnames <- paste("indvals", 2:12, sep="")

ivlabs_list <- foreach(i=2:12) %do% indval(x=ord.df, clustering =grps[,i], numitr = 1000)

names(ivlabs_list) <- listnames

# Create a data frame with three columns: number of clusters, average p, number of sig ind spp
ivlabs_df <- data.frame(nclusters=2:12, p_mean=NA, nsignif=NA, sumsignif=NA)
for(i in seq(length(ivlabs_list))){
  ivlabs_df[i,"p_mean"] <- mean(ivlabs_list[[i]]$pval, na.rm = T) 
  ivlabs_df[i, "nsignif"] <- length(which(ivlabs_list[[i]]$pval<=0.05))
  ivlabs_df[i,"sumsignif"] <- sum(ivlabs_list[[i]]$indval, na.rm = T)
}

par(mfrow=c(1, 2))
# plot average p values by number of clusters
plot(x=ivlabs_df$nclusters, y=ivlabs_df$p_mean, xlab="Number of clusters", ylab="Mean p-value")
# plot number of significant indicator species by number of clusters
plot(x=ivlabs_df$nclusters, y=ivlabs_df$nsignif, xlab="Number of clusters", ylab="Number of significant species")
par(mfrow=c(1, 1))

k_groups <- 4 # manually change based on plot results

Based on the results in Figure \@ref(fig:clustering), the optimal number of vegetation groupings is r k_groups (mean p-value = r round(ivlabs_df$p_mean[ivlabs_df$nclusters==k_groups], 3); number of indicator species = r ivlabs_df$nsignif[ivlabs_df$nclusters==k_groups]).

# create groups
grp <- cutree(hclust1, k=k_groups)

# create data frame that identifies indicator species by group:
ivlab <- indval(x=ord.df, clustering = grp, numitr = 1000)

gr <- ivlab$maxcls 
iv <- ivlab$indcls 
pv <- ivlab$pval 

indvalsummary <- data.frame(group=gr, indval=iv, pvalue=pv) #, freq=fr)
indvalsummary <- indvalsummary[order(indvalsummary$group, -indvalsummary$indval),]
indvalsummary$Code <- rownames(indvalsummary)

# add in readable indicator names
indvalsummary <- left_join(x=indvalsummary,
                           y=select(indicator_descriptions, Indicator_code, Indicator),
                           by=c("Code"="Indicator_code")) %>%
  left_join(x=.,
            y=select(read.csv(data_file_paths(user)$species_list), SpeciesCode, ScientificName),
            by=c("Code"="SpeciesCode")) %>%
  mutate(Indicator = ifelse(test = is.na(Indicator), yes = ScientificName, no = Indicator))

if(opuntia_combined){
  indvalsummary$Indicator[which(indvalsummary$Code=="AH_OpuntiaCover")] <- indicator_descriptions$Indicator[which(indicator_descriptions$Indicator_code=="opuntia_combined")]
}

# add in relative abundance
relab <- as.data.frame(ivlab$relabu)
colnames(relab) <- paste0("RelativeAbund_", 1:ncol(relab))
relab$Code <- rownames(relab)
indvalsummary <- left_join(x=indvalsummary, y=relab, by="Code")

# add in relative frequency
relfr <- as.data.frame(ivlab$relfrq)
colnames(relfr) <- paste0("RelativeFreq_", 1:ncol(relfr)) 
relfr$Code <- rownames(relfr)
indvalsummary <- left_join(x=indvalsummary, y=relfr, by="Code")

# create a table for display where relative abundance and relative frequency are only displayed for the group for which the species is an indicator
indvalsummary_disp <- dplyr::select(indvalsummary, group, indval, pvalue, Indicator)
indvalsummary_disp$RelAbund <- NA
for(g in 1:k_groups){
 indvalsummary_disp[indvalsummary_disp$group==g, "RelAbund"] <- 
   indvalsummary[indvalsummary$group==g, paste0("RelativeAbund_", g)]
}

indvalsummary_disp$RelFreq <- NA
for(g in 1:k_groups){
 indvalsummary_disp[indvalsummary_disp$group==g, "RelFreq"] <- 
   indvalsummary[indvalsummary$group==g, paste0("RelativeFreq_", g)]
}

# order data frame by group, then significance; for output file
#indvalsummary_disp <- indvalsummary_disp[order(indvalsummary_disp$group, indvalsummary_disp$indval), ]
indvalsummary_disp <- indvalsummary_disp %>%
  group_by(group) %>%
  arrange(desc(indval), .by_group=T)

indvalsummary_disp_filt <- filter(indvalsummary_disp, pvalue<=0.05) %>%
  mutate(indval=round(indval, 2),
         pvalue=round(pvalue, 3),
         RelAbund=round(RelAbund, 2),
         RelFreq=round(RelFreq, 2)) %>%
  select(group, Indicator, indval, pvalue, RelAbund, RelFreq) %>%
  setNames(c("Group", "Indicator", "Indicator value", "p-value", "Relative abundance", "Relative frequency"))

kbl(x = indvalsummary_disp_filt, format = "html", row.names = F, align = "c",
    caption = "Indicator species for plant community groups. The indicator value of a species is a combined measure of fidelity (always present in the group) and specificity (exclusive to the group), where higher values indicate a stronger association with the given vegetation group. Relative abundance indicates specificity, where 1 means the species only occurs in plots of the given group. Relative frequency indicates fidelity, where 1 means the species always occurs in plots of the given group. Only species with a p-value less than 0.05 are displayed.") %>%
  kable_styling(bootstrap_options = "bordered") %>%
  row_spec(kable_input=., row = 0:nrow(indvalsummary_disp_filt), background = "lightsteelblue") %>%
  collapse_rows(columns = 1, valign = "top")

# Cluster dendrogram plot
hclust1_dend <- as.dendrogram(hclust1)

# veg group color palette
pal_veg <- RColorBrewer::brewer.pal(n=k_groups, name = "Set1")

plot(hclust1_dend, leaflab= "none")
rect.dendrogram(hclust1_dend, 
                k=k_groups, 
                cluster = grp, 
                border = pal_veg[unique(cutree(hclust1_dend, k=k_groups)[order.dendrogram(hclust1_dend)])], 
                text = as.character(unique(cutree(hclust1_dend, k=k_groups)[order.dendrogram(hclust1_dend)])), 
                text_cex = 2, 
                lower_rect = -0.25,
                lwd =2)

#plot(ape::as.phylo(hclust1), tip.color = pal_veg9[grp], direction = "downwards")

NMDS_scree(ord.df)
ord_dims <- 2 # 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) 

Based on Figure \@ref(fig:ordination-scree), 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
plot_data_first$PlantCommunity <- as.factor(grp)

# manually update based on indicator species
group_labels <- c("1 - Perennial grassland with
mixed shrubs and forbs",
                    "2 - Juniper and shrub community
with non-native grasses",
                    "3 - Mixed grassland and shrub community",
                    "4 - Sparse shrub community")

par(mfrow=c(1,2))
# Axes 1x2
plot(ord, choices = c(1,2), type = "n", 
     xlim = c(-1.25, 1.25),
     ylim = c(-1, 1))
points(ord, choices = c(1,2), display = "sites",
       col=pal_veg[plot_data_first$PlantCommunity],
       pch = 21, cex = .6, bg=pal_veg[plot_data_first$PlantCommunity])
ordiellipse(ord, plot_data_first$PlantCommunity, col=pal_veg, lwd = 2, label = T)

# Legend
plot(ord, type = "n", axes=FALSE,
     display = "sites",
     col=pal_veg[plot_data_first$PlantCommunity],
     xlab = "",
     ylab = "")
legend(x="center",
       legend = group_labels,
       fill = pal_veg,
       title = "Plant communities")

par(mfrow=c(1,1))

State and transition model {.tabset}

CUSTOM DIAGRAM

# TODO create custom diagram of states
knitr::include_graphics("C:/Users/aknight/Desktop/Telework_Backups/V_drive/PROJECTS/ESG/STM/SW_SandyLoamyUplands_boxandarrow_DRAFT.JPG")

# TODO legend for custom diagram if needed

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 1

Reference State

TODO write state description

Community 1.1

Reference State

TODO write state description

# 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)")

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

Contributors

• Travis Nauman
• Samuel Burch
• Anna Knight
• 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

annack84/STMdevelopment documentation built on April 12, 2024, 6:46 p.m.