diff --git a/.Rbuildignore b/.Rbuildignore
index 0fb64fa..005a054 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,2 +1,3 @@
 ^aspen\.Rproj$
 ^\.Rproj\.user$
+^LICENSE\.md$
diff --git a/.gitignore b/.gitignore
index cd67eac..5b6a065 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,4 @@
 .Rproj.user
+.Rhistory
+.RData
+.Ruserdata
diff --git a/DESCRIPTION b/DESCRIPTION
index 26f2de6..3383d18 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,12 +1,20 @@
 Package: aspen
-Title: What the Package Does (One Line, Title Case)
+Title: Aspen Protocol Data Tools
 Version: 0.0.0.9000
-Authors@R: 
-    person("First", "Last", , "first.last@example.com", role = c("aut", "cre"),
-           comment = c(ORCID = "YOUR-ORCID-ID"))
-Description: What the package does (one paragraph).
-License: `use_mit_license()`, `use_gpl3_license()` or friends to pick a
-    license
+Authors@R: c(person("Sarah", "Wright", email = "sarah_wright@nps.gov", role=c("aut","cre"), comment = c(ORCID = "0009-0004-5060-2189")), 
+           person("Irene", "Foster", email = "irene_foster@partner.nps.gov", role=c("aut"), ,
+           comment = c(ORCID = "0000-0001-9681-4786")))
+Description: This R package is intended to provide automated QA/QC, visualization, and
+    basic analysis and statistical summaries of data collected for the Aspen
+    monitoring protocol.
+License: CC0
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.1
+Imports: 
+    dplyr,
+    fetchagol,
+    magrittr
+Suggests: 
+    testthat (>= 3.0.0)
+Config/testthat/edition: 3
diff --git a/LICENSE.md b/LICENSE.md
new file mode 100644
index 0000000..139c68e
--- /dev/null
+++ b/LICENSE.md
@@ -0,0 +1,43 @@
+## creative commons
+
+# CC0 1.0 Universal
+
+CREATIVE COMMONS CORPORATION IS NOT A LAW FIRM AND DOES NOT PROVIDE LEGAL SERVICES. DISTRIBUTION OF THIS DOCUMENT DOES NOT CREATE AN ATTORNEY-CLIENT RELATIONSHIP. CREATIVE COMMONS PROVIDES THIS INFORMATION ON AN "AS-IS" BASIS. CREATIVE COMMONS MAKES NO WARRANTIES REGARDING THE USE OF THIS DOCUMENT OR THE INFORMATION OR WORKS PROVIDED HEREUNDER, AND DISCLAIMS LIABILITY FOR DAMAGES RESULTING FROM THE USE OF THIS DOCUMENT OR THE INFORMATION OR WORKS PROVIDED HEREUNDER.
+
+### Statement of Purpose
+
+The laws of most jurisdictions throughout the world automatically confer exclusive Copyright and Related Rights (defined below) upon the creator and subsequent owner(s) (each and all, an "owner") of an original work of authorship and/or a database (each, a "Work").
+
+Certain owners wish to permanently relinquish those rights to a Work for the purpose of contributing to a commons of creative, cultural and scientific works ("Commons") that the public can reliably and without fear of later claims of infringement build upon, modify, incorporate in other works, reuse and redistribute as freely as possible in any form whatsoever and for any purposes, including without limitation commercial purposes. These owners may contribute to the Commons to promote the ideal of a free culture and the further production of creative, cultural and scientific works, or to gain reputation or greater distribution for their Work in part through the use and efforts of others.
+
+For these and/or other purposes and motivations, and without any expectation of additional consideration or compensation, the person associating CC0 with a Work (the "Affirmer"), to the extent that he or she is an owner of Copyright and Related Rights in the Work, voluntarily elects to apply CC0 to the Work and publicly distribute the Work under its terms, with knowledge of his or her Copyright and Related Rights in the Work and the meaning and intended legal effect of CC0 on those rights.
+
+1. __Copyright and Related Rights.__ A Work made available under CC0 may be protected by copyright and related or neighboring rights ("Copyright and Related Rights"). Copyright and Related Rights include, but are not limited to, the following:
+
+    i. the right to reproduce, adapt, distribute, perform, display, communicate, and translate a Work;
+
+    ii. moral rights retained by the original author(s) and/or performer(s);
+
+    iii. publicity and privacy rights pertaining to a person's image or likeness depicted in a Work;
+
+    iv. rights protecting against unfair competition in regards to a Work, subject to the limitations in paragraph 4(a), below;
+
+    v. rights protecting the extraction, dissemination, use and reuse of data in a Work;
+
+    vi. database rights (such as those arising under Directive 96/9/EC of the European Parliament and of the Council of 11 March 1996 on the legal protection of databases, and under any national implementation thereof, including any amended or successor version of such directive); and
+
+    vii. other similar, equivalent or corresponding rights throughout the world based on applicable law or treaty, and any national implementations thereof.
+
+2. __Waiver.__ To the greatest extent permitted by, but not in contravention of, applicable law, Affirmer hereby overtly, fully, permanently, irrevocably and unconditionally waives, abandons, and surrenders all of Affirmer's Copyright and Related Rights and associated claims and causes of action, whether now known or unknown (including existing as well as future claims and causes of action), in the Work (i) in all territories worldwide, (ii) for the maximum duration provided by applicable law or treaty (including future time extensions), (iii) in any current or future medium and for any number of copies, and (iv) for any purpose whatsoever, including without limitation commercial, advertising or promotional purposes (the "Waiver"). Affirmer makes the Waiver for the benefit of each member of the public at large and to the detriment of Affirmer's heirs and successors, fully intending that such Waiver shall not be subject to revocation, rescission, cancellation, termination, or any other legal or equitable action to disrupt the quiet enjoyment of the Work by the public as contemplated by Affirmer's express Statement of Purpose.
+
+3. __Public License Fallback.__ Should any part of the Waiver for any reason be judged legally invalid or ineffective under applicable law, then the Waiver shall be preserved to the maximum extent permitted taking into account Affirmer's express Statement of Purpose. In addition, to the extent the Waiver is so judged Affirmer hereby grants to each affected person a royalty-free, non transferable, non sublicensable, non exclusive, irrevocable and unconditional license to exercise Affirmer's Copyright and Related Rights in the Work (i) in all territories worldwide, (ii) for the maximum duration provided by applicable law or treaty (including future time extensions), (iii) in any current or future medium and for any number of copies, and (iv) for any purpose whatsoever, including without limitation commercial, advertising or promotional purposes (the "License"). The License shall be deemed effective as of the date CC0 was applied by Affirmer to the Work. Should any part of the License for any reason be judged legally invalid or ineffective under applicable law, such partial invalidity or ineffectiveness shall not invalidate the remainder of the License, and in such case Affirmer hereby affirms that he or she will not (i) exercise any of his or her remaining Copyright and Related Rights in the Work or (ii) assert any associated claims and causes of action with respect to the Work, in either case contrary to Affirmer's express Statement of Purpose.
+
+4. __Limitations and Disclaimers.__
+
+    a. No trademark or patent rights held by Affirmer are waived, abandoned, surrendered, licensed or otherwise affected by this document.
+
+    b. Affirmer offers the Work as-is and makes no representations or warranties of any kind concerning the Work, express, implied, statutory or otherwise, including without limitation warranties of title, merchantability, fitness for a particular purpose, non infringement, or the absence of latent or other defects, accuracy, or the present or absence of errors, whether or not discoverable, all to the greatest extent permissible under applicable law.
+
+    c. Affirmer disclaims responsibility for clearing rights of other persons that may apply to the Work or any use thereof, including without limitation any person's Copyright and Related Rights in the Work. Further, Affirmer disclaims responsibility for obtaining any necessary consents, permissions or other rights required for any use of the Work.
+
+    d. Affirmer understands and acknowledges that Creative Commons is not a party to this document and has no duty or obligation with respect to this CC0 or use of the Work.
diff --git a/NAMESPACE b/NAMESPACE
index 5dbc7f9..8bfd85f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,8 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
-export(fetchAGOLToken)
-export(fetchAllRecords)
-export(fetchLayerAndTableList)
-export(fetchMetadata)
+export(loadAndWrangleMOJNAspen)
+export(loadAndWrangleUCBNAspen)
+export(writeAspen)
+import(dplyr)
+import(fetchagol)
 importFrom(magrittr,"%<>%")
 importFrom(magrittr,"%>%")
diff --git a/R/aspen_qc.R b/R/aspen_qc.R
new file mode 100644
index 0000000..e00443f
--- /dev/null
+++ b/R/aspen_qc.R
@@ -0,0 +1,80 @@
+
+#' Check that there is at least one tree for every observation entry and flag entries with more than 250 trees
+#' @param data Aspen dataset to run through QC functions
+missingTreeQC <- function(data = fetchAndWrangleAspen()){
+
+  missingTreeQC <- data$data$Observations %>%
+    # Add up all trees in a row
+    dplyr::mutate(totalTreeCount = Class1+Class2+Class3+Class4+Class5+Class6) %>%
+    # Filter for any rows that have less than one trees or more than 250 trees
+    dplyr::filter(totalTreeCount < 1 | totalTreeCount > 250) %>%
+    dplyr::select(Park, Site, VisitDate, SpeciesCode, Class1, Class2, Class3, Class4, Class5, Class6, totalTreeCount)
+
+  return(missingTreeQC)
+}
+
+#' Check that no entries in the observation table have a SpeciesCode that is NA or missing
+#' @param data Aspen dataset to run through QC functions
+treeSpeciesQC <- function(data = fetchAndWrangleAspen()){
+
+  treeSpeciesQC <- data$data$Observations %>%
+    dplyr::filter(SpeciesCode == "" | is.na(SpeciesCode) | SpeciesCode == "NA") %>%
+    dplyr::select(Park, Site, VisitDate, SpeciesCode, Class1, Class2, Class3, Class4, Class5, Class6)
+
+  return(treeSpeciesQC)
+}
+
+#' Check that unknown species entries do not have any live trees
+#' @param data Aspen dataset to run through QC functions
+unknownSpeciesQC <- function(data = fetchAndWrangleAspen()){
+
+  unknownSpeciesQC <- data$data$Observations %>%
+    dplyr::filter(SpeciesCode == "UNK") %>%
+    dplyr::filter(Class1 != 0 | Class2 != 0 | Class3 != 0 | Class4 != 0 | Class5 != 0) %>%
+    dplyr::select(Park, Site, VisitDate, SpeciesCode, Class1, Class2, Class3, Class4, Class5, Class6)
+
+  return(unknownSpeciesQC)
+}
+
+#' Check that all plots have at least one live aspen each time they are visited
+#' @param data Aspen dataset to run through QC functions
+checkAspenQC <- function(data = fetchAndWrangleAspen()){
+
+  checkAspenQC <- data$data$Observations %>%
+    # Filter for aspen trees
+    dplyr::filter(SpeciesCode == "POTR5") %>%
+    dplyr::select(parentglobalid, SpeciesCode, Class1, Class2, Class3, Class4, Class5) %>%
+    # Join observations and visits
+    dplyr::full_join(data$data$SiteVisit, by = c("parentglobalid" = "globalid")) %>%
+    # Filter data for any visits that don't have aspen observations
+    dplyr::filter(is.na(SpeciesCode)) %>%
+    dplyr::select(Park, Site, VisitType, VisitDate, SpeciesCode)
+
+  return(checkAspenQC)
+}
+
+#' Flag large changes in number of tree between years
+#' @param data Aspen dataset to run through QC functions
+treeChangeQC <- function(data = fetchAndWrangleAspen()){
+
+  visits <- data$data$SiteVisit %>%
+    dplyr::group_by(Site) %>%
+    dplyr::arrange(Site, VisitDate) %>%
+    dplyr::mutate(visitNumber = seq_along(Site)) %>%
+    dplyr::select("globalid", "Park", "Site", "VisitType", "VisitDate", "visitNumber")
+
+  treeChangeQC <- data$data$Observations %>%
+    # join with site visit table
+    #dplyr::left_join(visits, by = c("parentglobalid" = "globalid")) %>%
+    dplyr::mutate(totalTrees = Class1+Class2+Class3+Class4+Class5) %>%
+    dplyr::group_by(Site, VisitDate) %>%
+    dplyr::summarise(totalTrees = sum(totalTrees)) %>%
+    dplyr::group_by(Site) %>%
+    dplyr::mutate(previousTrees = dplyr::lag(totalTrees),
+                  treeChange = previousTrees - totalTrees)
+
+    # find percent change (in total number of trees both conifer and aspen or separate??) between year ((year 1 - year 2)/year 1)
+    # flag any with a large change (figure out what counts as a large change)
+}
+
+# Flag if there's more than one visit to a plot in a year
diff --git a/R/aspen_qc_ucbn.R b/R/aspen_qc_ucbn.R
new file mode 100644
index 0000000..f8f8466
--- /dev/null
+++ b/R/aspen_qc_ucbn.R
@@ -0,0 +1,80 @@
+
+#' Check that there is at least one tree for every observation entry and flag entries with more than 250 trees
+#' @param data Aspen dataset to run through QC functions
+missingTreeUCBNQC <- function(data){
+
+  missingTreeQC <- data$data$Observations %>%
+    # Add up all trees in a row
+    dplyr::mutate(totalTreeCount = Class1+Class2+Class3+Class4+Class5+Class6) %>%
+    # Filter for any rows that have less than one trees or more than 250 trees
+    dplyr::filter(totalTreeCount < 1 | totalTreeCount > 250) %>%
+    dplyr::select(Park, Unique_ID, VisitDate, SpeciesCode, Class1, Class2, Class3, Class4, Class5, Class6, totalTreeCount)
+
+  return(missingTreeQC)
+}
+
+#' Check that no entries in the observation table have a SpeciesCode that is NA or missing
+#' @param data Aspen dataset to run through QC functions
+treeSpeciesUCBNQC <- function(data){
+
+  treeSpeciesQC <- data$data$Observations %>%
+    dplyr::filter(SpeciesCode == "" | is.na(SpeciesCode) | SpeciesCode == "NA") %>%
+    dplyr::select(Park, Unique_ID, VisitDate, SpeciesCode, Class1, Class2, Class3, Class4, Class5, Class6)
+
+  return(treeSpeciesQC)
+}
+
+#' Check that unknown species entries do not have any live trees
+#' @param data Aspen dataset to run through QC functions
+unknownSpeciesUCBNQC <- function(data){
+
+  unknownSpeciesQC <- data$data$Observations %>%
+    dplyr::filter(SpeciesCode == "UNK") %>%
+    dplyr::filter(Class1 != 0 | Class2 != 0 | Class3 != 0 | Class4 != 0 | Class5 != 0) %>%
+    dplyr::select(Park, Unique_ID, VisitDate, SpeciesCode, Class1, Class2, Class3, Class4, Class5, Class6)
+
+  return(unknownSpeciesQC)
+}
+
+#' Check that all plots have at least one live aspen for each visit
+#' @param data Aspen dataset to run through QC functions
+checkAspenUCBNQC <- function(data){
+
+  checkAspenQC <- data$data$Observations %>%
+    # Filter for aspen trees
+    dplyr::filter(SpeciesCode == "POPTRE") %>%
+    dplyr::select(parentglobalid, SpeciesCode, Class1, Class2, Class3, Class4, Class5) %>%
+    # Join observations and visits
+    dplyr::full_join(data$data$SiteVisit, by = c("parentglobalid" = "globalid")) %>%
+    # Filter data for any visits that don't have aspen observations
+    dplyr::filter(is.na(SpeciesCode)) %>% # Any visit without aspens should have a NA in SpeciesCode after joining with only aspen observations
+    dplyr::select(Park, Unique_ID, VisitDate, SpeciesCode)
+
+  return(checkAspenQC)
+}
+
+#' Flag large changes in number of tree between years
+#' @param data Aspen dataset to run through QC functions
+treeChangeUCBNQC <- function(data){
+
+  visits <- data$data$SiteVisit %>%
+    dplyr::group_by(Unique_ID) %>%
+    dplyr::arrange(Unique_ID, VisitDate) %>%
+    dplyr::mutate(visitNumber = seq_along(Unique_ID)) %>%
+    dplyr::select("globalid", "Park", "Unique_ID", "VisitDate", "visitNumber")
+
+  treeChangeQC <- data$data$Observations %>%
+    # join with site visit table
+    #dplyr::left_join(visits, by = c("parentglobalid" = "globalid")) %>%
+    dplyr::mutate(totalTrees = Class1+Class2+Class3+Class4+Class5) %>%
+    dplyr::group_by(Unique_ID, VisitDate) %>%
+    dplyr::summarise(totalTrees = sum(totalTrees)) %>%
+    dplyr::group_by(Unique_ID) %>%
+    dplyr::mutate(previousTrees = dplyr::lag(totalTrees),
+                  treeChange = previousTrees - totalTrees)
+
+    # find percent change (in total number of trees both conifer and aspen or separate??) between year ((year 1 - year 2)/year 1)
+    # flag any with a large change (figure out what counts as a large change)
+}
+
+# Flag if there's more than one visit to a plot in a year
diff --git a/R/utils.R b/R/utils.R
new file mode 100644
index 0000000..cbbbee2
--- /dev/null
+++ b/R/utils.R
@@ -0,0 +1,185 @@
+#' @importFrom magrittr %>% %<>%
+#' @import dplyr
+#' @import fetchagol
+
+pkg_globals <- new.env(parent = emptyenv())
+
+
+#' Fetch aspen data from AGOL and do preliminary data wrangling
+#'
+#' @param aspen_url URL to main AGOL aspen database
+#' @param site_url URL to AGOL database for aspen sites (if applicable)
+#' @param agol_username Authentication token (not needed for public layers)
+#'
+#' @return A list of data frames and metadata
+#' @export
+
+# TODO: add variables that were joined to data tables to each metadata table
+loadAndWrangleMOJNAspen <- function(
+    aspen_url = "https://services1.arcgis.com/fBc8EJBxQRMcHlei/arcgis/rest/services/MOJN_Aspen_Test_Visit_NonSpatial_gdb/FeatureServer",
+                                 site_url =  "https://services1.arcgis.com/fBc8EJBxQRMcHlei/arcgis/rest/services/AspenSites2/FeatureServer",
+                                 agol_username = "mojn_data") {
+  flattened_data <- list(data = list(),
+                         metadata = list())
+
+  # Import aspen database
+  raw_data <- fetchagol::fetchRawData(aspen_url, agol_username)
+
+
+  # Imports optional second database and connects it to main database
+  # For MOJN - there's a master list of sites in a different database
+  if(!is.null(site_url)) {
+    # Import aspen site database
+    raw_site_data <- fetchagol::fetchRawData(site_url, agol_username)
+
+    # Add site data to list of other data
+    raw_data$data$AllSites <- raw_site_data$data$`MOJN Aspen Sites Master`
+    # Join site metadata to list of other metadata
+    raw_data$metadata$AllSites <- raw_site_data$metadata$`MOJN Aspen Sites Master`
+  }
+
+  raw_data <- fetchagol::cleanData(raw_data)
+
+
+  flattened_data$metadata <- raw_data$metadata
+
+  # Add optional second database to flattened data
+  if(!is.null(site_url)) {
+  flattened_data$data$AllSites <- raw_data$data$AllSites
+
+  # Join background site information with other tables
+  raw_data$data$SiteVisit <- raw_data$data$SiteVisit %>%
+    dplyr::left_join(dplyr::select(flattened_data$data$AllSites, dplyr::any_of(c("Site", "Status", "Panel", "Stratum", "Zone_", "Community"))),
+                     by = dplyr::join_by("Site"))
+  flattened_data$metadata$SiteVisit$fields <- append(flattened_data$metadata$SiteVisit$fields, flattened_data$metadata$AllSites$fields[c("Status", "Panel", "Stratum", "Zone_", "Community")])
+  flattened_data$metadata$AllSites$table_name <- "AllSites"
+
+  }
+
+
+  # Join background site information with other tables
+  # flattened_data$data$SiteVisit <- raw_data$data$SiteVisit %>%
+  #   dplyr::left_join(dplyr::select(flattened_data$data$AllSites, dplyr::any_of(c("Site", "Status", "Panel", "Stratum", "Zone_", "Community"))),
+  #                    by = dplyr::join_by("Site"))
+  flattened_data$data$SiteVisit <- raw_data$data$SiteVisit
+
+  flattened_data$data$Disturbances <- raw_data$data$SiteVisit %>%
+    dplyr::mutate(parentglobalid = globalid) %>%
+    dplyr::select(parentglobalid, Park, Site, VisitType, VisitDate, FieldSeason) %>%
+    dplyr::right_join(raw_data$data$Disturbances, by = c("parentglobalid" = "parentglobalid"))
+  flattened_data$metadata$Disturbances$fields <- append(flattened_data$metadata$Disturbances$fields, flattened_data$metadata$SiteVisit$fields[c("Park", "VisitType", "VisitDate", "FieldSeason")])
+
+
+  flattened_data$data$Observations <- raw_data$data$SiteVisit %>%
+    dplyr::mutate(parentglobalid = globalid) %>%
+    dplyr::select(parentglobalid, Park, Site, VisitType, VisitDate, FieldSeason) %>%
+    dplyr::right_join(raw_data$data$Observations, by = c("parentglobalid" = "parentglobalid"))
+  flattened_data$metadata$Observations$fields <- append(flattened_data$metadata$Observations$fields, flattened_data$metadata$SiteVisit$fields[c("Park", "Site", "VisitType", "VisitDate", "FieldSeason")])
+
+  flattened_data$data$Pests <- flattened_data$data$Observations %>%
+    dplyr::mutate(parentglobalid = globalid) %>%
+    dplyr::select(parentglobalid, Park, Site, VisitType, VisitDate, FieldSeason, SpeciesCode) %>%
+    dplyr::right_join(raw_data$data$Pests, by = c("parentglobalid" = "parentglobalid"))
+  flattened_data$metadata$Pests$fields <- append(flattened_data$metadata$Pests$fields, flattened_data$metadata$Observations$fields[c("Park", "VisitType", "VisitDate", "FieldSeason", "SpeciesCode")])
+
+  invisible(flattened_data)
+}
+
+#' Fetch aspen data from AGOL and do preliminary data wrangling
+#'
+#' @param aspen_url URL to main AGOL aspen database
+#' @param site_url URL to AGOL database for aspen sites (if applicable)
+#' @param agol_username Authentication token (not needed for public layers)
+#'
+#' @return A list of data frames and metadata
+#' @export
+
+loadAndWrangleUCBNAspen <- function(
+    aspen_url = "https://services1.arcgis.com/fBc8EJBxQRMcHlei/arcgis/rest/services/service_4d6343e9204142928351c52c6f1362c5/FeatureServer",
+                                 site_url =  NULL,
+                                 agol_username = "mojn_data") {
+  flattened_data <- list(data = list(),
+                         metadata = list())
+
+  # Import aspen database
+  raw_data <- fetchagol::fetchRawData(aspen_url, agol_username)
+
+
+  # Imports optional second database and connects it to main database
+  # For MOJN - there's a master list of sites in a different database
+  if(!is.null(site_url)) {
+    # Import aspen site database
+    raw_site_data <- fetchagol::fetchRawData(site_url, agol_username)
+
+    # Add site data to list of other data
+    raw_data$data$AllSites <- raw_site_data$data$`MOJN Aspen Sites Master`
+    # Join site metadata to list of other metadata
+    raw_data$metadata$AllSites <- raw_site_data$metadata$`MOJN Aspen Sites Master`
+  }
+
+  raw_data <- fetchagol::cleanData(raw_data)
+
+  # Add optional second database to flattened data
+  if(!is.null(site_url)) {
+    flattened_data$data$AllSites <- raw_data$data$AllSites
+
+    flattened_data$metadata <- raw_data$metadata
+
+    # Join background site information with other tables and add background information metadata to tables it was added to
+    raw_data$data$SiteVisit <- raw_data$data$SiteVisit %>%
+      dplyr::left_join(dplyr::select(flattened_data$data$AllSites, dplyr::any_of(c("Site", "Status", "Panel", "Stratum", "Zone_", "Community"))),
+                       by = dplyr::join_by("Site"))}
+  flattened_data$metadata$SiteVisit$fields <- append(flattened_data$metadata$SiteVisit$fields, flattened_data$metadata$AllSites$fields[c("Site", "Status", "Panel", "Stratum", "Zone_", "Community")])
+
+
+  # Join background site information with other tables
+  # flattened_data$data$SiteVisit <- raw_data$data$SiteVisit %>%
+  #   dplyr::left_join(dplyr::select(flattened_data$data$AllSites, dplyr::any_of(c("Site", "Status", "Panel", "Stratum", "Zone_", "Community"))),
+  #                    by = dplyr::join_by("Site"))
+  flattened_data$data$SiteVisit <- raw_data$data$SiteVisit
+
+  flattened_data$data$Disturbances <- raw_data$data$SiteVisit %>%
+    dplyr::mutate(parentglobalid = globalid) %>%
+    dplyr::select(parentglobalid, Park, Stand, Transect, VisitDate, PlotNum, Unique_ID) %>%
+    dplyr::right_join(raw_data$data$Disturbances, by = c("parentglobalid" = "parentglobalid"))
+  flattened_data$metadata$Disturbances$fields <- append(flattened_data$metadata$Disturbances$fields, flattened_data$metadata$SiteVisit$fields[c("Park", "Stand", "Transect", "VisitDate", "PlotNum", "Unique_ID")])
+
+  flattened_data$data$Observations <- raw_data$data$SiteVisit %>%
+    dplyr::mutate(parentglobalid = globalid) %>%
+    dplyr::select(parentglobalid, Park, Stand, Transect, VisitDate, PlotNum, Unique_ID) %>%
+    dplyr::right_join(raw_data$data$Observations, by = c("parentglobalid" = "parentglobalid"))
+  flattened_data$metadata$Observations$fields <- append(flattened_data$metadata$Observations$fields, flattened_data$metadata$SiteVisit$fields[c("Park", "Stand", "Transect", "VisitDate", "PlotNum", "Unique_ID")])
+
+  flattened_data$data$Pests <- flattened_data$data$Observations %>%
+    dplyr::mutate(parentglobalid = globalid) %>%
+    dplyr::select(parentglobalid, Park, Stand, Transect, VisitDate, PlotNum, Unique_ID) %>%
+    dplyr::right_join(raw_data$data$Pests, by = c("parentglobalid" = "parentglobalid"))
+  flattened_data$metadata$Pests$fields <- append(flattened_data$metadata$Pests$fields, flattened_data$metadata$SiteVisit$fields[c("Park", "Stand", "Transect", "VisitDate", "PlotNum", "Unique_ID")])
+
+  invisible(flattened_data)
+}
+
+#' Write aspen data to CSV
+#'
+#' @inheritParams fetchagol::writeToFiles
+#'
+#' @export
+#'
+writeAspen <- function(all_data, data_dir = here::here("data", "final"), dictionary_dir = here::here("data", "dictionary"),
+                      dictionary_filenames = c(tables = "data_dictionary_tables.txt",
+                                               attributes = "data_dictionary_attributes.txt",
+                                               categories = "data_dictionary_categories.txt"),
+                      verbose = FALSE, removeColumns = TRUE, cols_to_remove = c("Editor", "Creator"))
+{
+  fetchagol::writeToFiles(all_data = all_data, data_dir = data_dir, dictionary_dir = dictionary_dir, lookup_dir = NA, verbose = verbose, removeColumns = TRUE, cols_to_remove = c("Editor", "Creator"))
+}
+
+
+
+# TODO: should a get_data() function be added like in pine??
+
+
+
+
+
+
diff --git a/README.Rmd b/README.Rmd
new file mode 100644
index 0000000..c885ad2
--- /dev/null
+++ b/README.Rmd
@@ -0,0 +1,42 @@
+---
+title: "Aspen Guide"
+output: github_document
+date: "2024-06-20"
+---
+
+The goal of the aspen package is to provide data QC, visualization, and publication assistance for MOJN aspen data.
+
+### Installation
+
+If you already have a GitHub token for RStudio run the code below to download the aspen package. If you haven’t generated a GitHub token for RStudio, you can follow the guide [here](https://frbcesab.github.io/rsetup/chapters/github-token.html) to set it up and then run the following code to download the aspen package.
+
+```{r, eval=FALSE}
+# install.packages("devtools")
+devtools::install_github(repo = "nationalparkservice/mojn-aspen-rpackage")
+```
+
+
+### Metadata and AGOL Password
+
+Before running any of the functions in the package you should check that the data in you aspen AGOL database is filled out. Instructions for filling out AGOL metadata can be found [here](https://github.com/nationalparkservice/imd-fetchagol-package. 
+
+Before running the functions, it is also recommended that you save the password to your AGOL headless account using the keyring package in R. To add headless account information to the default keyring follow the steps below. The code should only have to be run once.
+
+
+```{r, eval = FALSE}
+# Run the function below. Change the username field to the username of your headless account and input the password when prompted
+keyring::key_set(service = "AGOL", username = "USERNAME", prompt = "Password for headless account: ")
+
+# To check the key was saved run code below, filling in the username to your headless account
+keyring::key_get(service = "AGOL", username = "USERNAME")
+```
+
+### Functions
+
+A short summary of what this package contains. For more information see function documentation.
+
+- loadAndWrangleMOJNAspen(): load MOJN aspen data from AGOL into R and perform some basic data wrangling
+- loadAndWrangleUCBNAspen(): load UCBN aspen data from AGOL into R and perform some basic data wrangling
+- writeAspen(): write aspen data and metadata to CSVs
+
+In addition to these functions the aspen package has some quality control functions, a quality control script, and a visualization script. To learn more about the functions see package documentation. You can find the scripts in the script folder of the [GitHub repo](https://github.com/nationalparkservice/mojn-aspen-rpackage) or if you cloned the GitHub repo to your desktop you can open the files from there. Download, open, and then run the .qmd file to see the results. 
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..9b385b2
--- /dev/null
+++ b/README.md
@@ -0,0 +1,58 @@
+Aspen Guide
+================
+2024-06-20
+
+The goal of the aspen package is to provide data QC, visualization, and
+publication assistance for MOJN aspen data.
+
+### Installation
+
+If you already have a GitHub token for RStudio run the code below to
+download the aspen package. If you haven’t generated a GitHub token for
+RStudio, you can follow the guide
+[here](https://frbcesab.github.io/rsetup/chapters/github-token.html) to
+set it up and then run the following code to download the aspen package.
+
+``` r
+# install.packages("devtools")
+devtools::install_github(repo = "nationalparkservice/mojn-aspen-rpackage")
+```
+
+### Metadata and AGOL Password
+
+Before running any of the functions in the package you should check that
+the data in you aspen AGOL database is filled out. Instructions for
+filling out AGOL metadata can be found
+\[here\](<https://github.com/nationalparkservice/imd-fetchagol-package>.
+
+Before running the functions, it is also recommended that you save the
+password to your AGOL headless account using the keyring package in R.
+To add headless account information to the default keyring follow the
+steps below. The code should only have to be run once.
+
+``` r
+# Run the function below. Change the username field to the username of your headless account and input the password when prompted
+keyring::key_set(service = "AGOL", username = "USERNAME", prompt = "Password for headless account: ")
+
+# To check the key was saved run code below, filling in the username to your headless account
+keyring::key_get(service = "AGOL", username = "USERNAME")
+```
+
+### Functions
+
+A short summary of what this package contains. For more information see
+function documentation.
+
+- loadAndWrangleMOJNAspen(): load MOJN aspen data from AGOL into R and
+  perform some basic data wrangling
+- loadAndWrangleUCBNAspen(): load UCBN aspen data from AGOL into R and
+  perform some basic data wrangling
+- writeAspen(): write aspen data and metadata to CSVs
+
+In addition to these functions the aspen package has some quality
+control functions, a quality control script, and a visualization script.
+To learn more about the functions see package documentation. You can
+find the scripts in the script folder of the [GitHub
+repo](https://github.com/nationalparkservice/mojn-aspen-rpackage) or if
+you cloned the GitHub repo to your desktop you can open the files from
+there. Download, open, and then run the .qmd file to see the results.
diff --git a/man/checkAspenQC.Rd b/man/checkAspenQC.Rd
new file mode 100644
index 0000000..2e5d423
--- /dev/null
+++ b/man/checkAspenQC.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aspen_qc.R
+\name{checkAspenQC}
+\alias{checkAspenQC}
+\title{Check that all plots have at least one live aspen each time they are visited}
+\usage{
+checkAspenQC(data = fetchAndWrangleAspen())
+}
+\arguments{
+\item{data}{Aspen dataset to run through QC functions}
+}
+\description{
+Check that all plots have at least one live aspen each time they are visited
+}
diff --git a/man/checkAspenUCBNQC.Rd b/man/checkAspenUCBNQC.Rd
new file mode 100644
index 0000000..6e133bf
--- /dev/null
+++ b/man/checkAspenUCBNQC.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aspen_qc_ucbn.R
+\name{checkAspenUCBNQC}
+\alias{checkAspenUCBNQC}
+\title{Check that all plots have at least one live aspen for each visit}
+\usage{
+checkAspenUCBNQC(data)
+}
+\arguments{
+\item{data}{Aspen dataset to run through QC functions}
+}
+\description{
+Check that all plots have at least one live aspen for each visit
+}
diff --git a/man/fetchAGOLToken.Rd b/man/fetchAGOLToken.Rd
deleted file mode 100644
index 07e98d2..0000000
--- a/man/fetchAGOLToken.Rd
+++ /dev/null
@@ -1,28 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils.R
-\name{fetchAGOLToken}
-\alias{fetchAGOLToken}
-\title{Get AGOL authentication token}
-\usage{
-fetchAGOLToken(
-  agol_username,
-  agol_password = keyring::key_get(service = "AGOL", username = agol_username),
-  root = "nps.maps.arcgis.com",
-  referer = "https://irma.nps.gov"
-)
-}
-\arguments{
-\item{agol_username}{AGOL headless account username}
-
-\item{agol_password}{AGOL headless account password (do not hard code this into your scripts!)}
-
-\item{root}{NPS users should keep the default. See \url{https://developers.arcgis.com/rest/users-groups-and-items/root.htm} for more information.}
-
-\item{referer}{NPS users should keep the default. See \url{https://developers.arcgis.com/rest/users-groups-and-items/generate-token.htm} for more information.}
-}
-\value{
-An AGOL authentication token
-}
-\description{
-Treat this token as you would a password: don't hard-code it in your scripts or save it to a file. It will expire after 60 minutes.
-}
diff --git a/man/fetchAllRecords.Rd b/man/fetchAllRecords.Rd
deleted file mode 100644
index 56e4a68..0000000
--- a/man/fetchAllRecords.Rd
+++ /dev/null
@@ -1,34 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils.R
-\name{fetchAllRecords}
-\alias{fetchAllRecords}
-\title{Fetch tabular data from AGOL}
-\usage{
-fetchAllRecords(
-  data_path,
-  layer_number,
-  token,
-  geometry = FALSE,
-  where = "1=1",
-  outFields = "*"
-)
-}
-\arguments{
-\item{data_path}{Feature service URL}
-
-\item{layer_number}{Layer number}
-
-\item{token}{Authentication token (not needed for public layers)}
-
-\item{geometry}{Include spatial data columns? Works with points, not tested with other geometry types}
-
-\item{where}{Query clause specifying a subset of rows (optional; defaults to all rows). See AGOL REST API documentation.}
-
-\item{outFields}{String indicating which fields to return (optional; defaults to all fields). See AGOL REST API documentation.}
-}
-\value{
-A tibble
-}
-\description{
-Retrieves tabular data from AGOL layers and tables, even when number of rows exceeds maximum record count.
-}
diff --git a/man/fetchLayerAndTableList.Rd b/man/fetchLayerAndTableList.Rd
deleted file mode 100644
index 15261de..0000000
--- a/man/fetchLayerAndTableList.Rd
+++ /dev/null
@@ -1,19 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils.R
-\name{fetchLayerAndTableList}
-\alias{fetchLayerAndTableList}
-\title{Fetch feature service info from AGOL}
-\usage{
-fetchLayerAndTableList(url, token)
-}
-\arguments{
-\item{url}{Feature service URL}
-
-\item{token}{Authentication token (not needed for public layers)}
-}
-\value{
-A list
-}
-\description{
-Retrieves metadata from AGOL layers and tables.
-}
diff --git a/man/fetchMetadata.Rd b/man/fetchMetadata.Rd
deleted file mode 100644
index 54215e5..0000000
--- a/man/fetchMetadata.Rd
+++ /dev/null
@@ -1,21 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils.R
-\name{fetchMetadata}
-\alias{fetchMetadata}
-\title{Fetch metadata from AGOL}
-\usage{
-fetchMetadata(url, token, layer_number)
-}
-\arguments{
-\item{url}{Feature service URL}
-
-\item{token}{Authentication token (not needed for public layers)}
-
-\item{layer_number}{Optional layer ID}
-}
-\value{
-A list
-}
-\description{
-Retrieves metadata from AGOL layers and tables.
-}
diff --git a/man/loadAndWrangleMOJNAspen.Rd b/man/loadAndWrangleMOJNAspen.Rd
new file mode 100644
index 0000000..04239d1
--- /dev/null
+++ b/man/loadAndWrangleMOJNAspen.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{loadAndWrangleMOJNAspen}
+\alias{loadAndWrangleMOJNAspen}
+\title{Fetch aspen data from AGOL and do preliminary data wrangling}
+\usage{
+loadAndWrangleMOJNAspen(
+  aspen_url =
+    "https://services1.arcgis.com/fBc8EJBxQRMcHlei/arcgis/rest/services/MOJN_Aspen_Test_Visit_NonSpatial_gdb/FeatureServer",
+  site_url =
+    "https://services1.arcgis.com/fBc8EJBxQRMcHlei/arcgis/rest/services/AspenSites2/FeatureServer",
+  agol_username = "mojn_data"
+)
+}
+\arguments{
+\item{aspen_url}{URL to main AGOL aspen database}
+
+\item{site_url}{URL to AGOL database for aspen sites (if applicable)}
+
+\item{agol_username}{Authentication token (not needed for public layers)}
+}
+\value{
+A list of data frames and metadata
+}
+\description{
+Fetch aspen data from AGOL and do preliminary data wrangling
+}
diff --git a/man/loadAndWrangleUCBNAspen.Rd b/man/loadAndWrangleUCBNAspen.Rd
new file mode 100644
index 0000000..f9a494c
--- /dev/null
+++ b/man/loadAndWrangleUCBNAspen.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{loadAndWrangleUCBNAspen}
+\alias{loadAndWrangleUCBNAspen}
+\title{Fetch aspen data from AGOL and do preliminary data wrangling}
+\usage{
+loadAndWrangleUCBNAspen(
+  aspen_url =
+    "https://services1.arcgis.com/fBc8EJBxQRMcHlei/arcgis/rest/services/service_4d6343e9204142928351c52c6f1362c5/FeatureServer",
+  site_url = NULL,
+  agol_username = "mojn_data"
+)
+}
+\arguments{
+\item{aspen_url}{URL to main AGOL aspen database}
+
+\item{site_url}{URL to AGOL database for aspen sites (if applicable)}
+
+\item{agol_username}{Authentication token (not needed for public layers)}
+}
+\value{
+A list of data frames and metadata
+}
+\description{
+Fetch aspen data from AGOL and do preliminary data wrangling
+}
diff --git a/man/missingTreeQC.Rd b/man/missingTreeQC.Rd
new file mode 100644
index 0000000..5916796
--- /dev/null
+++ b/man/missingTreeQC.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aspen_qc.R
+\name{missingTreeQC}
+\alias{missingTreeQC}
+\title{Check that there is at least one tree for every observation entry and flag entries with more than 250 trees}
+\usage{
+missingTreeQC(data = fetchAndWrangleAspen())
+}
+\arguments{
+\item{data}{Aspen dataset to run through QC functions}
+}
+\description{
+Check that there is at least one tree for every observation entry and flag entries with more than 250 trees
+}
diff --git a/man/missingTreeUCBNQC.Rd b/man/missingTreeUCBNQC.Rd
new file mode 100644
index 0000000..cee2dbd
--- /dev/null
+++ b/man/missingTreeUCBNQC.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aspen_qc_ucbn.R
+\name{missingTreeUCBNQC}
+\alias{missingTreeUCBNQC}
+\title{Check that there is at least one tree for every observation entry and flag entries with more than 250 trees}
+\usage{
+missingTreeUCBNQC(data)
+}
+\arguments{
+\item{data}{Aspen dataset to run through QC functions}
+}
+\description{
+Check that there is at least one tree for every observation entry and flag entries with more than 250 trees
+}
diff --git a/man/treeChangeQC.Rd b/man/treeChangeQC.Rd
new file mode 100644
index 0000000..31a386d
--- /dev/null
+++ b/man/treeChangeQC.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aspen_qc.R
+\name{treeChangeQC}
+\alias{treeChangeQC}
+\title{Flag large changes in number of tree between years}
+\usage{
+treeChangeQC(data = fetchAndWrangleAspen())
+}
+\arguments{
+\item{data}{Aspen dataset to run through QC functions}
+}
+\description{
+Flag large changes in number of tree between years
+}
diff --git a/man/treeChangeUCBNQC.Rd b/man/treeChangeUCBNQC.Rd
new file mode 100644
index 0000000..4d2854d
--- /dev/null
+++ b/man/treeChangeUCBNQC.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aspen_qc_ucbn.R
+\name{treeChangeUCBNQC}
+\alias{treeChangeUCBNQC}
+\title{Flag large changes in number of tree between years}
+\usage{
+treeChangeUCBNQC(data)
+}
+\arguments{
+\item{data}{Aspen dataset to run through QC functions}
+}
+\description{
+Flag large changes in number of tree between years
+}
diff --git a/man/treeSpeciesQC.Rd b/man/treeSpeciesQC.Rd
new file mode 100644
index 0000000..aaca94f
--- /dev/null
+++ b/man/treeSpeciesQC.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aspen_qc.R
+\name{treeSpeciesQC}
+\alias{treeSpeciesQC}
+\title{Check that no entries in the observation table have a SpeciesCode that is NA or missing}
+\usage{
+treeSpeciesQC(data = fetchAndWrangleAspen())
+}
+\arguments{
+\item{data}{Aspen dataset to run through QC functions}
+}
+\description{
+Check that no entries in the observation table have a SpeciesCode that is NA or missing
+}
diff --git a/man/treeSpeciesUCBNQC.Rd b/man/treeSpeciesUCBNQC.Rd
new file mode 100644
index 0000000..b3b2bbe
--- /dev/null
+++ b/man/treeSpeciesUCBNQC.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aspen_qc_ucbn.R
+\name{treeSpeciesUCBNQC}
+\alias{treeSpeciesUCBNQC}
+\title{Check that no entries in the observation table have a SpeciesCode that is NA or missing}
+\usage{
+treeSpeciesUCBNQC(data)
+}
+\arguments{
+\item{data}{Aspen dataset to run through QC functions}
+}
+\description{
+Check that no entries in the observation table have a SpeciesCode that is NA or missing
+}
diff --git a/man/unknownSpeciesQC.Rd b/man/unknownSpeciesQC.Rd
new file mode 100644
index 0000000..5f7bc1c
--- /dev/null
+++ b/man/unknownSpeciesQC.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aspen_qc.R
+\name{unknownSpeciesQC}
+\alias{unknownSpeciesQC}
+\title{Check that unknown species entries do not have any live trees}
+\usage{
+unknownSpeciesQC(data = fetchAndWrangleAspen())
+}
+\arguments{
+\item{data}{Aspen dataset to run through QC functions}
+}
+\description{
+Check that unknown species entries do not have any live trees
+}
diff --git a/man/unknownSpeciesUCBNQC.Rd b/man/unknownSpeciesUCBNQC.Rd
new file mode 100644
index 0000000..e04c4ff
--- /dev/null
+++ b/man/unknownSpeciesUCBNQC.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aspen_qc_ucbn.R
+\name{unknownSpeciesUCBNQC}
+\alias{unknownSpeciesUCBNQC}
+\title{Check that unknown species entries do not have any live trees}
+\usage{
+unknownSpeciesUCBNQC(data)
+}
+\arguments{
+\item{data}{Aspen dataset to run through QC functions}
+}
+\description{
+Check that unknown species entries do not have any live trees
+}
diff --git a/man/writeAspen.Rd b/man/writeAspen.Rd
new file mode 100644
index 0000000..d61bb3a
--- /dev/null
+++ b/man/writeAspen.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{writeAspen}
+\alias{writeAspen}
+\title{Write aspen data to CSV}
+\usage{
+writeAspen(
+  all_data,
+  data_dir = here::here("data", "final"),
+  dictionary_dir = here::here("data", "dictionary"),
+  dictionary_filenames = c(tables = "data_dictionary_tables.txt", attributes =
+    "data_dictionary_attributes.txt", categories = "data_dictionary_categories.txt"),
+  verbose = FALSE,
+  removeColumns = TRUE,
+  cols_to_remove = c("Editor", "Creator")
+)
+}
+\arguments{
+\item{all_data}{Output of \code{fetchRawData()}}
+
+\item{data_dir}{Folder to store data csv's in}
+
+\item{dictionary_dir}{Folder to store data dictionaries in}
+
+\item{dictionary_filenames}{Named list with names \code{c("tables", "attributes", "categories")} indicating what to name the tables, attributes, and categories data dictionaries. You are encouraged to keep the default names unless you have a good reason to change them.}
+
+\item{verbose}{Output feedback to console?}
+
+\item{removeColumns}{Should columns be removed?}
+
+\item{cols_to_remove}{Columns that should be removed, \code{c("Editor", "Creator")} is the default because they can contain personally identifiable information}
+}
+\description{
+Write aspen data to CSV
+}
diff --git a/mojn-aspen-rpackage.Rproj b/mojn-aspen-rpackage.Rproj
index 497f8bf..270314b 100644
--- a/mojn-aspen-rpackage.Rproj
+++ b/mojn-aspen-rpackage.Rproj
@@ -18,3 +18,4 @@ StripTrailingWhitespace: Yes
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace
diff --git a/scripts/visualizations/aspen_qc_results.qmd b/scripts/visualizations/aspen_qc_results.qmd
new file mode 100644
index 0000000..b2c9146
--- /dev/null
+++ b/scripts/visualizations/aspen_qc_results.qmd
@@ -0,0 +1,63 @@
+---
+title: "Aspen Quality Control"
+format: 
+  html:
+    embed-resources: true
+execute: 
+  echo: false
+  message: false
+  warning: false
+---
+
+```{r}
+library(tidyverse)
+library(aspen)
+```
+
+
+##### Check Trees
+Check that there is at least one tree for every observation entry and flag entries with more than 250 trees
+
+```{r}
+if(nrow(aspen:::missingTreeQC()) == 0){
+  cat('No errors')
+} else{
+  rmarkdown::paged_table(aspen:::missingTreeQC())
+}
+```
+
+
+##### Check for Missing Species
+Check that no entries in the observation table have a SpeciesCode that is NA or missing
+
+```{r}
+if(nrow(aspen:::treeSpeciesQC()) == 0){
+  cat('No errors')
+} else{
+  rmarkdown::paged_table(aspen:::treeSpeciesQC())
+}
+```
+
+
+##### Check Unknown Trees
+Check that unknown species entries do not have any live trees
+
+```{r}
+if(nrow(aspen:::unknownSpeciesQC()) == 0){
+  cat('No errors')
+} else{
+  rmarkdown::paged_table(aspen:::unknownSpeciesQC())
+}
+```
+
+
+##### Check for Aspens
+Check that all of the plots have at least one aspen each time they are visited
+
+```{r}
+if(nrow(aspen:::checkAspenQC()) == 0){
+  cat('No errors')
+} else{
+  rmarkdown::paged_table(aspen:::checkAspenQC())
+}
+```
diff --git a/scripts/visualizations/aspen_visualizations.qmd b/scripts/visualizations/aspen_visualizations.qmd
new file mode 100644
index 0000000..59f72c9
--- /dev/null
+++ b/scripts/visualizations/aspen_visualizations.qmd
@@ -0,0 +1,1113 @@
+---
+title: "Aspen Preliminary Visualizations"
+format: 
+  html:
+    embed-resources: true
+execute: 
+  echo: false
+  message: false
+  warning: false
+---
+
+```{r, message=FALSE}
+library(tidyverse)
+library(khroma)
+library(viridis)
+library(plotly)
+aspen <- aspen::fetchAndWrangleAspen()
+```
+
+```{r}
+# find number of plots in each watershed
+
+# TODO: update this to Allsites tables once not used sites are removed
+
+aspen$data$AllSites <- aspen$data$AllSites %>%
+  filter(Site %in% aspen$data$SiteVisit$Site) %>%
+      # Find Number of plot in each watershed
+  group_by(Stratum) %>%
+  mutate(numberOfPlotsInWatershed = n()) %>%
+  ungroup() %>%
+  group_by(Community) %>%
+  # Find number of plots in each community
+  mutate(numberOfPlotsInCommunity = n()) %>%
+  ungroup()
+```
+
+
+# Vegetation Community
+
+#### Average Density of Live Aspens
+
+```{r}
+# 4 meter radius of plots
+# area = pi*r^2
+# density = #/ha
+
+# this is finding the density for each plot and then averaging it for the communities
+
+aspenDensity2 <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POTR5" & VisitType == "Primary") %>%
+  # Add up LIVE POTR5 for each row
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/(pi * 4^2/10000),
+         totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Community) %>%
+  summarise(avgDensity = mean(density),
+            sd = sd(density), 
+            numberOfPlots = n()) %>%
+  # Find values for error bars
+  mutate(ymin = avgDensity - sd,
+         ymax = avgDensity + sd)
+
+aspenDensityGraph <- aspenDensity2 %>%
+  ggplot(aes(x= Community, y = avgDensity, label = numberOfPlots)) +
+  geom_col() +
+  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.2) +
+  theme_minimal() +
+  labs(title = "Average Aspen Density",
+       x = "Community", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+plotly::ggplotly(aspenDensityGraph)
+```
+
+
+```{r}
+# # 4 meter radius of plots
+# # area = pi*r^2
+# # density = #/ha
+# 
+# 
+# # this is the density of the total community (add up total number of trees in community and divide by total area of plots)
+# aspenDensity <- aspen$data$Observations %>%
+#   left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+#   filter(VisitType == "Primary") %>%
+#   # filter for POTR5
+#   filter(SpeciesCode == "POTR5") %>%
+#   # Add up LIVE POTR5 for each row
+#   mutate(totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+#   # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+#   group_by(Community) %>%
+#   summarise(density = sum(totalLiveTrees)/((first(numberOfPlotsInCommunity) * pi * 4^2)/10000),
+#             sumLiveTrees = sum(totalLiveTrees),
+#             numberOfPlotsInCommunity = first(numberOfPlotsInCommunity))
+# 
+# 
+# aspenDensityGraph <- aspenDensity %>%
+#   ggplot(aes(x= Community, y = density)) +
+#   geom_col() +
+#   theme_minimal() +
+#   labs(title = "Aspen Density",
+#        x = "Community", y = "Density (#/ha)") + 
+#   theme(plot.title = element_text(hjust = 0.5)) +
+#   theme(axis.text.x = element_text(size = 6, angle = 315))
+# 
+# plotly::ggplotly(aspenDensityGraph)
+```
+
+#### Density of Live Aspens per Plot
+
+```{r}
+aspenDensity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POTR5" & VisitType == "Primary") %>%
+  # Find density for each plot
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/((pi * 4^2)/10000)) %>%
+# Filter out communities with not enough points to graph
+group_by(Community) %>%
+  mutate(numberOfPlots = n()) %>%
+  ungroup() %>%
+  filter(numberOfPlots >2)
+
+
+aspenDensityGraph <- aspenDensity %>%
+  ggplot(aes(x= Community, y = density, fill = Community, label = numberOfPlots)) +
+    geom_violin() +
+    geom_boxplot(width=0.1, color="grey", alpha=0.5) +
+  stat_summary(fun = "mean", colour = "red", size = 2, geom = "point") +
+  scale_fill_viridis(discrete = TRUE) +
+  theme_minimal() +
+  labs(title = "Aspen Density",
+       x = "Community", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0.0)) + 
+  guides(fill="none")
+
+aspenDensityGraph
+```
+
+
+```{r}
+aspenDensityGraph <- aspenDensity %>%
+  ggplot(aes(x= Community, y = density, label = numberOfPlots)) +
+    #geom_violin() +
+  geom_boxplot(width=0.3, outlier.shape = NA) +
+  geom_point(alpha = 0.5, position=position_jitter(height=.5, width=.25)) +
+  theme_minimal() +
+  labs(title = "Aspen Density",
+       x = "Community", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0.0)) + 
+  guides(fill="none")
+
+aspenDensityGraph
+```
+
+
+```{r}
+fig <- aspenDensity %>%
+  plotly::plot_ly(
+    x = ~Community,
+    y = ~density,
+    split = ~Community,
+    color = ~Community,
+    colors = "Dark2",
+    type = 'violin',
+    box = list(visible = T),
+    meanline = list(visible = T)) 
+
+fig <- fig %>%
+  layout(
+    xaxis = list(
+      title = "Community"),
+    yaxis = list(
+      title = "Density",
+      zeroline = F,
+      range = list(0, 40000)),
+    showlegend = FALSE)
+
+fig
+```
+
+
+
+#### Average Density of Conifers
+
+```{r}
+# 4 meter radius of plots
+# area = pi*r^2
+# density = #/ha
+
+# this is finding the density for each plot and then averaging it for the communities
+
+aspenDensity2 <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+  # filter for POTR5
+  filter(SpeciesCode != "POTR5" & SpeciesCode != "BEUC2" & VisitType == "Primary") %>%
+  # Aggregate to get the total number of conifers at each site
+  group_by(Site, Community) %>%
+  summarise(Class1 = sum(Class1),
+            Class2 = sum(Class2),
+            Class3 = sum(Class3),
+            Class4 = sum(Class4),
+            Class5 = sum(Class5)) %>%
+  ungroup() %>%
+    # Add up LIVE POTR5 for each row
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/(pi * 4^2/10000),
+         totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Community) %>%
+  summarise(avgDensity = mean(density),
+            sd = sd(density), 
+            numberOfPlots = n()) %>%
+  # Find values for error bars
+  mutate(ymin = avgDensity - sd,
+         ymax = avgDensity + sd)
+
+aspenDensityGraph <- aspenDensity2 %>%
+  ggplot(aes(x= Community, y = avgDensity, label = numberOfPlots)) +
+  geom_col() +
+  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.2) +
+  theme_minimal() +
+  labs(title = "Average Conifer Density",
+       x = "Community", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+plotly::ggplotly(aspenDensityGraph)
+```
+
+
+```{r}
+# coniferDensity <- aspen$data$Observations %>%
+#   left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+#   # filter for POTR5
+#   filter(SpeciesCode != "POTR5" & SpeciesCode != "BEUC2") %>%
+#   # Add up LIVE POTR5 for each row
+#   mutate(totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+#   # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+#   group_by(Community) %>%
+#   summarise(density = sum(totalLiveTrees)/((first(numberOfPlotsInCommunity) * pi * 4^2)/10000),
+#             sumLiveTrees = sum(totalLiveTrees),
+#             numberOfPlotsInCommunity = first(numberOfPlotsInCommunity))
+# 
+# 
+# coniferDensityGraph <- coniferDensity %>%
+#   ggplot(aes(x= Community, y = density)) +
+#   geom_col() +
+#   theme_minimal() +
+#   labs(title = "Conifer Density",
+#        x = "Community", y = "Density (#/ha)") + 
+#   theme(plot.title = element_text(hjust = 0.5)) +
+#   theme(axis.text.x = element_text(size = 6, angle = 315))
+# 
+# plotly::ggplotly(coniferDensityGraph)
+```
+
+
+#### Density of Conifers per Plot
+
+```{r}
+coniferDensity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+  # filter for POTR5
+  filter(SpeciesCode != "POTR5" & SpeciesCode != "BEUC2" & VisitType == "Primary") %>%
+    # Aggregate to get the total number of conifers at each site
+  group_by(Site, Community) %>%
+  summarise(Class1 = sum(Class1),
+            Class2 = sum(Class2),
+            Class3 = sum(Class3),
+            Class4 = sum(Class4),
+            Class5 = sum(Class5)) %>%
+  ungroup() %>%
+  # Find density for each plot
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/((pi * 4^2)/10000)) %>%
+# Filter out communities with not enough points to graph
+group_by(Community) %>%
+  mutate(numberOfPlots = n()) %>%
+  ungroup() %>%
+  filter(numberOfPlots >2)
+
+
+coniferDensityGraph <- coniferDensity %>%
+  ggplot(aes(x= Community, y = density, fill = Community, label = numberOfPlots)) +
+    geom_violin() +
+    geom_boxplot(width=0.1, color="grey", alpha=0.2) +
+  stat_summary(fun = "mean", colour = "red", size = 2, geom = "point") +
+  scale_fill_viridis(discrete = TRUE) +
+  theme_minimal() +
+  labs(title = "Conifer Density",
+       x = "Community", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0.0)) + 
+  guides(fill="none")
+
+coniferDensityGraph
+```
+
+#### Density of Aspen vs Conifers
+
+```{r}
+compareDensity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+    # filter out one species that is not aspen or conifer
+  filter(SpeciesCode != "BEUC2" & VisitType == "Primary") %>%
+  mutate(isAspen = case_when(
+    SpeciesCode == "POTR5" ~ 'Aspen',
+    .default = 'Conifer'
+  )) %>%
+    # Aggregate to get the total number of conifers at each site
+  group_by(Site, Community, isAspen) %>%
+  summarise(Class1 = sum(Class1),
+            Class2 = sum(Class2),
+            Class3 = sum(Class3),
+            Class4 = sum(Class4),
+            Class5 = sum(Class5)) %>%
+  ungroup() %>%
+  # Add up live trees for each row
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/(pi * 4^2/10000)) %>%
+  group_by(Community, isAspen) %>%
+  summarise(avgDensity = mean(density),
+            sd = sd(density), 
+            numberOfPlots = n()) %>%
+  # Find values for error bars
+  mutate(ymin = avgDensity - sd,
+         ymax = avgDensity + sd)
+
+
+compareDensityGraph <- compareDensity %>%
+  ggplot(aes(x= Community, y = avgDensity, fill = as.factor(isAspen), label = numberOfPlots)) +
+  geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.2, position = position_dodge(0.9)) +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))  +
+  scale_fill_light()
+
+plotly::ggplotly(compareDensityGraph)
+```
+
+
+#### Density of Aspen vs Conifers per Plot
+
+```{r}
+compareDensity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+    # filter out one species that is not aspen or conifer
+  filter(SpeciesCode != "BEUC2" & VisitType == "Primary") %>%
+  mutate(isAspen = case_when(
+    SpeciesCode == "POTR5" ~ 'Aspen',
+    .default = 'Conifer'
+  )) %>% 
+      # Aggregate to get the total number of conifers at each site
+  group_by(Site, Community, isAspen) %>%
+  summarise(Class1 = sum(Class1),
+            Class2 = sum(Class2),
+            Class3 = sum(Class3),
+            Class4 = sum(Class4),
+            Class5 = sum(Class5)) %>%
+  ungroup() %>%
+    # Find density for each plot
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/((pi * 4^2)/10000)) %>%
+  filter(Community != "Curl-leaf Mountain-mahogany Shrubland & Woodland Complex")
+
+compareDensityGraph <- compareDensity %>%
+  ggplot(aes(x= Community, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))  +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))
+
+compareDensityGraph
+```
+
+#### Density of Aspen Age Groups
+
+```{r}
+ageGroupCommunity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POTR5" & VisitType == "Primary") %>%
+  # Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  mutate(density = treeCount/(pi * 4^2/10000)) %>%
+  # summarize POTR5 for each age group divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Community, class) %>%
+    summarise(avgDensity = mean(density),
+            sd = sd(density), 
+            numberOfPlots = n()) %>%
+  # Find values for error bars
+  mutate(ymin = avgDensity - sd,
+         ymax = avgDensity + sd)
+
+ageGroupCommunityGraph <- ageGroupCommunity %>%
+  ggplot(aes(x= Community, y = avgDensity, fill = factor(class))) +
+  geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.3, position = position_dodge(0.9)) +
+  theme_minimal() +
+  labs(title = "Age Class Density",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Age Class") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+plotly::ggplotly(ageGroupCommunityGraph)
+```
+
+
+#### Density of Aspen Age Groups Boxplot
+
+```{r}
+ageGroupCommunity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POTR5" & VisitType == "Primary") %>%
+  # Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  mutate(density = treeCount/(pi * 4^2/10000)) %>%
+  filter(Community != "Curl-leaf Mountain-mahogany Shrubland & Woodland Complex")
+
+
+ageGroupCommunityGraph <- ageGroupCommunity %>%
+  ggplot(aes(x= Community, y = density, fill = factor(class))) +
+  geom_boxplot() +
+  theme_minimal() +
+  labs(title = "Age Class Density",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Age Class") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))
+
+ageGroupCommunityGraph
+```
+
+
+
+```{r}
+# #### Density of Aspen Age Groups Line Graph
+# 
+# ageGroupCommunity <- aspen$data$Observations %>%
+#   left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+#   # filter for POTR5
+#   filter(SpeciesCode == "POTR5") %>%
+#   # Pivot longer to make calculations easier
+#   pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+#   # summarize POTR5 for each age group divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+#   group_by(Community, class) %>%
+#   summarise(density = sum(treeCount)/((first(numberOfPlotsInCommunity) * pi * 4^2)/10000),
+#             sumLiveTrees = sum(treeCount),
+#             numberOfPlotsInCommunity = first(numberOfPlotsInCommunity)) 
+# 
+# 
+# ageGroupCommunityGraph <- ageGroupCommunity %>%
+#   ggplot(aes(x= class, y = density, color = factor(Community))) +
+#   geom_line(aes(group = factor(Community))) + 
+#   #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+#   theme_minimal() +
+#   labs(title = "Age Class Density",
+#        x = "Age Class", y = "Density (#/ha)",
+#        color = "Community") + 
+#   theme(plot.title = element_text(hjust = 0.5)) +
+#   scale_fill_light() + 
+#   guides(color="none")
+# # +
+# #   theme(axis.text.x = element_text(size = 6, angle = 315))
+# 
+# #ageGroupCommunityGraph
+# plotly::ggplotly(ageGroupCommunityGraph)
+```
+
+
+#### Density of Conifer Age Groups
+
+```{r}
+ageGroupCommunity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+  # filter for POTR5
+  filter(SpeciesCode != "POTR5" & SpeciesCode != "BEUC2" & VisitType == "Primary") %>%
+  # Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  mutate(density = treeCount/(pi * 4^2/10000)) 
+# %>%
+#   # summarize POTR5 for each age group divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+#   group_by(Community, class) %>%
+#     summarise(avgDensity = mean(density),
+#             sd = sd(density), 
+#             numberOfPlots = n()) %>%
+#   # Find values for error bars
+#   mutate(ymin = avgDensity - sd,
+#          ymax = avgDensity + sd)
+
+ageGroupCommunityGraph <- ageGroupCommunity %>%
+  ggplot(aes(x= Community, y = density, fill = factor(class))) +
+  geom_boxplot() +
+  # geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  # geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.5, position = position_dodge(0.9)) +
+  theme_minimal() +
+  labs(title = "Age Class Density",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Age Class") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+ageGroupCommunityGraph
+```
+
+
+#### Density of Aspen vs Conifer Age Groups
+
+```{r}
+compareAgeDensity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+    # filter out one species that is not aspen or conifer
+  filter(SpeciesCode != "BEUC2" & VisitType == "Primary") %>%
+  mutate(isAspen = case_when(
+    SpeciesCode == "POTR5" ~ 'Aspen',
+    .default = 'Conifer'
+  )) %>%
+    # Find density for each plot
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/((pi * 4^2)/10000)) %>%  
+# Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  mutate(density = treeCount/(pi * 4^2/10000))
+```
+
+::: panel-tabset
+##### Class 1
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class1") %>%
+  ggplot(aes(x= Community, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  scale_fill_viridis(discrete = TRUE) +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 1",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0)) +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+##### Class 2
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class2") %>%
+  ggplot(aes(x= Community, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 2",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))  +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+
+##### Class 3
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class3") %>%
+  ggplot(aes(x= Community, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 3",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))  +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+
+##### Class 4
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class4") %>%
+  ggplot(aes(x= Community, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 4",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))  +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+
+##### Class 5
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class5") %>%
+  ggplot(aes(x= Community, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 5",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))  +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+
+##### Class 6
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class6") %>%
+  ggplot(aes(x= Community, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 6",
+       x = "Community", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))  +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+
+:::
+
+#### Frequency of Disturbances
+
+```{r}
+# Area of plots in which a disturbance occurs / total area of plots (for each community)
+
+disturbancesFrequency <- aspen$data$Disturbances %>%
+  group_by(Site) %>%
+  # Group data so every site with a disturbance has only one entry
+  summarise(Disturbance = first(Disturbance)) %>%
+  # Join disturbance and all site data
+  full_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+  group_by(Community) %>%
+  # Find out how many plots in a watershed were disturbed
+  summarise(plotsWithDisturbances = sum(!is.na(Disturbance)), 
+            numberOfPlotsInCommunity = first(numberOfPlotsInCommunity),
+            frequency = (first(plotsWithDisturbances) * pi * 4^2)/(first(numberOfPlotsInCommunity) * pi * 4^2)*100)
+
+disturbancesGraph <- disturbancesFrequency %>%
+  ggplot(aes(x = Community, y = frequency)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Disturbance Frequency",
+       x = "Community", y = "Frequency") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+plotly::ggplotly(disturbancesGraph)
+```
+
+
+#### Frequency of Pests
+
+```{r}
+# Area of plots in which a disturbance occurs / total area of plots (for each community)
+
+pestFrequency <- aspen$data$Pests %>%
+  group_by(Site) %>%
+  # Group data so every site with a disturbance has only one entry
+  summarise(Pest = first(Pest)) %>%
+  # Join disturbance and all site data
+  full_join(select(aspen$data$AllSites, any_of(c("Site", "Community", "numberOfPlotsInCommunity"))), by = join_by(Site)) %>%
+  group_by(Community) %>%
+  # Find out how many plots in a watershed were disturbed
+  summarise(plotsWithPests = sum(!is.na(Pest)), 
+            numberOfPlotsInCommunity = first(numberOfPlotsInCommunity),
+            frequency = (first(plotsWithPests) * pi * 4^2)/(first(numberOfPlotsInCommunity) * pi * 4^2)*100)
+
+pestGraph <- pestFrequency %>%
+  ggplot(aes(x = Community, y = frequency)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Pest Frequency",
+       x = "Community", y = "Frequency") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+plotly::ggplotly(pestGraph)
+```
+
+
+# Watershed (Stratum)
+
+#### Density of Live Aspens
+
+```{r}
+# 4 meter radius of plots
+# area = pi*r^2
+# denisty = #/ha
+
+
+aspenDensity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Stratum", "numberOfPlotsInWatershed"))), by = join_by(Site)) %>%
+  # Find Number of plot in each watershed
+  # group_by(Stratum) %>%
+  # mutate(numberOfPlots = n_distinct(Site)) %>%
+  # ungroup() %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POTR5") %>%
+  # Add up LIVE POTR5 for each row
+  mutate(totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Stratum) %>%
+  summarise(density = sum(totalLiveTrees)/((first(numberOfPlotsInWatershed) * pi * 4^2)/10000),
+            sumLiveTrees = sum(totalLiveTrees),
+            numberOfPlotsInWatershed = first(numberOfPlotsInWatershed))
+
+
+aspenDensityGraph <- aspenDensity %>%
+  ggplot(aes(x= Stratum, y = density)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Aspen Density",
+       x = "Watershed", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(aspenDensityGraph)
+```
+
+
+#### Density of Conifers
+
+```{r}
+  
+coniferDensity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Stratum", "numberOfPlotsInWatershed"))), by = join_by(Site)) %>%
+  # Find Number of plot in each watershed
+  # group_by(Stratum) %>%
+  # mutate(numberOfPlotsInWaterShed = n_distinct(Site)) %>%
+  # ungroup() %>%
+  # filter for POTR5
+  filter(SpeciesCode != "POTR5" & SpeciesCode != "BEUC2") %>%
+  # Add up LIVE POTR5 for each row
+  mutate(totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Stratum) %>%
+  summarise(density = sum(totalLiveTrees)/((first(numberOfPlotsInWatershed) * pi * 4^2)/10000),
+            sumLiveTrees = sum(totalLiveTrees),
+            numberOfPlotsInWatershed = first(numberOfPlotsInWatershed))
+
+
+coniferDensityGraph <- coniferDensity %>%
+  ggplot(aes(x= Stratum, y = density)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Conifer Density",
+       x = "Watershed", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(coniferDensityGraph)
+```
+
+
+#### Density of Aspen vs Conifers
+
+```{r}
+compareDensity <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Stratum", "numberOfPlotsInWatershed"))), by = join_by(Site)) %>%
+    # filter out one species that is not aspen or conifer
+  filter(SpeciesCode != "BEUC2") %>%
+  mutate(isAspen = case_when(
+    SpeciesCode == "POTR5" ~ 'Aspen',
+    .default = 'Conifer'
+  )) %>%
+  # Add up live trees for each row
+  mutate(liveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Stratum, isAspen) %>%
+  summarise(density = sum(liveTrees)/((first(numberOfPlotsInWatershed) * pi * 4^2)/10000),
+            totalLiveTrees = sum(liveTrees),
+            numberOfPlotsInWatershed = first(numberOfPlotsInWatershed))
+
+
+compareDensityGraph <- compareDensity %>%
+  ggplot(aes(x= Stratum, y = density, fill = as.factor(isAspen))) +
+  geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density",
+       x = "Watershed", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()
+
+plotly::ggplotly(compareDensityGraph)
+```
+
+
+#### Density of Aspen Age Groups
+
+```{r}
+ageGroupWatershed <- aspen$data$Observations %>%
+  left_join(select(aspen$data$AllSites, any_of(c("Site", "Stratum", "numberOfPlotsInWatershed"))), by = join_by(Site)) %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POTR5") %>%
+  # Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  # summarize POTR5 for each age group divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Stratum, class) %>%
+  summarise(density = sum(treeCount)/((first(numberOfPlotsInWatershed) * pi * 4^2)/10000),
+            sumLiveTrees = sum(treeCount),
+            numberOfPlotsInWatershed = first(numberOfPlotsInWatershed)) 
+
+
+ageGroupWatershedGraph <- ageGroupWatershed %>%
+  ggplot(aes(x= Stratum, y = density, fill = factor(class))) +
+  geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Age Class Density",
+       x = "Watershed", y = "Density (#/ha)",
+       fill = "Age Class") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()
+
+plotly::ggplotly(ageGroupWatershedGraph)
+```
+
+
+
+
+#### Frequency of Disturbances
+
+```{r}
+# Area of plots in which a disturbance occurs / total area of plots (for each watershed)
+
+disturbancesFrequency <- aspen$data$Disturbances %>%
+  group_by(Site) %>%
+  # Group data so every site with a disturbance has only one entry
+  summarise(Disturbance = first(Disturbance)) %>%
+  # Join disturbance and all site data
+  full_join(select(aspen$data$AllSites, any_of(c("Site", "Stratum", "numberOfPlotsInWatershed"))), by = join_by(Site)) %>%
+  group_by(Stratum) %>%
+  # Find out how many plots in a watershed were disturbed
+  summarise(plotsWithDisturbances = sum(!is.na(Disturbance)), 
+            numberOfPlotsInWatershed = first(numberOfPlotsInWatershed),
+            frequency = (first(plotsWithDisturbances) * pi * 4^2)/(first(numberOfPlotsInWatershed) * pi * 4^2)*100)
+
+disturbancesGraph <- disturbancesFrequency %>%
+  ggplot(aes(x = Stratum, y = frequency)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Disturbance Frequency",
+       x = "Watershed", y = "Frequency") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(disturbancesGraph)
+```
+
+
+#### Frequency of Pests
+
+```{r}
+# Area of plots in which a disturbance occurs / total area of plots (for each watershed)
+
+pestFrequency <- aspen$data$Pests %>%
+  group_by(Site) %>%
+  # Group data so every site with a disturbance has only one entry
+  summarise(Pest = first(Pest)) %>%
+  # Join disturbance and all site data
+  full_join(select(aspen$data$AllSites, any_of(c("Site", "Stratum", "numberOfPlotsInWatershed"))), by = join_by(Site)) %>%
+  group_by(Stratum) %>%
+  # Find out how many plots in a watershed were disturbed
+  summarise(plotsWithPests = sum(!is.na(Pest)), 
+            numberOfPlotsInWatershed = first(numberOfPlotsInWatershed),
+            frequency = (first(plotsWithPests) * pi * 4^2)/(first(numberOfPlotsInWatershed) * pi * 4^2)*100)
+
+pestGraph <- pestFrequency %>%
+  ggplot(aes(x = Stratum, y = frequency)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Pest Frequency",
+       x = "Watershed", y = "Frequency") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(pestGraph)
+```
+
+
+# Park
+
+#### Density of Live Aspens
+
+```{r}
+# 4 meter radius of plots
+# area = pi*r^2
+# denisty = #/ha
+
+
+  
+aspenParkDensity <- aspen$data$Observations %>%
+  # Find number of plot in each park
+  group_by(Park) %>%
+  mutate(numberOfPlots = n_distinct(Site)) %>%
+  ungroup() %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POTR5") %>%
+  # Add up LIVE POTR5 for each row
+  mutate(totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Park) %>%
+  summarise(density = sum(totalLiveTrees)/((first(numberOfPlots) * pi * 4^2)/10000),
+            sumLiveTrees = sum(totalLiveTrees),
+            numberOfPlots = first(numberOfPlots))
+
+
+aspenParkDensityGraph <- aspenParkDensity %>%
+  ggplot(aes(x= Park, y = density)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Aspen Density",
+       x = "Park", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(aspenParkDensityGraph)
+```
+
+#### Density of Conifers
+
+```{r}
+# 4 meter radius of plots
+# area = pi*r^2
+# denisty = #/ha
+
+
+  
+coniferParkDensity <- aspen$data$Observations %>%
+  # Find number of plot in each park
+  group_by(Park) %>%
+  mutate(numberOfPlots = n_distinct(Site)) %>%
+  ungroup() %>%
+  # filter for POTR5
+  filter(SpeciesCode != "POTR5" & SpeciesCode != "BEUC2") %>%
+  # Add up LIVE POTR5 for each row
+  mutate(totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Park) %>%
+  summarise(density = sum(totalLiveTrees)/((first(numberOfPlots) * pi * 4^2)/10000),
+            sumLiveTrees = sum(totalLiveTrees),
+            numberOfPlots = first(numberOfPlots))
+
+
+coniferParkDensityGraph <- coniferParkDensity %>%
+  ggplot(aes(x= Park, y = density)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Conifer Density",
+       x = "Park", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(coniferParkDensityGraph)
+```
+
+
+#### Density of Aspen vs Conifers
+
+```{r}
+compareDensity <- aspen$data$Observations %>%
+    # Find number of plot in each park
+  group_by(Park) %>%
+  mutate(numberOfPlots = n_distinct(Site)) %>%
+  ungroup() %>%
+  # filter out one species that is not aspen or conifer
+  filter(SpeciesCode != "BEUC2") %>%
+  mutate(isAspen = case_when(
+    SpeciesCode == "POTR5" ~ 'Aspen',
+    .default = 'Conifer'
+  )) %>%
+  # Add up live trees for each row
+  mutate(liveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Park, isAspen) %>%
+  summarise(density = sum(liveTrees)/((first(numberOfPlots) * pi * 4^2)/10000),
+            totalLiveTrees = sum(liveTrees),
+            numberOfPlots = first(numberOfPlots))
+
+
+compareDensityGraph <- compareDensity %>%
+  ggplot(aes(x= Park, y = density, fill = as.factor(isAspen))) +
+  geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density",
+       x = "Park", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()
+
+plotly::ggplotly(compareDensityGraph)
+```
+
+#### Density of Aspen Age Groups
+
+```{r}
+ageGroupPark <- aspen$data$Observations %>%
+  # Find number of plot in each park
+  group_by(Park) %>%
+  mutate(numberOfPlots = n_distinct(Site)) %>%
+  ungroup() %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POTR5") %>%
+  # Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  # summarize POTR5 for each age group divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Park, class) %>%
+  summarise(density = sum(treeCount)/((first(numberOfPlots) * pi * 4^2)/10000),
+            sumLiveTrees = sum(treeCount),
+            numberOfPlots = first(numberOfPlots)) 
+
+
+ageGroupParkGraph <- ageGroupPark %>%
+  ggplot(aes(x= Park, y = density, fill = factor(class))) +
+  geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Age Class Density",
+       x = "Park", y = "Density (#/ha)",
+       fill = "Age Class") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()
+
+plotly::ggplotly(ageGroupParkGraph)
+```
+
+
+#### Frequency of Disturbances
+
+```{r}
+# Area of plots in which a disturbance occurs / total area of plots (for each park)
+
+disturbancesFrequency <- aspen$data$Disturbances %>%
+  group_by(Site) %>%
+  # Group data so every site with a disturbance has only one entry
+  summarise(Disturbance = first(Disturbance)) %>%
+  # Join disturbance and all site data
+  full_join(select(aspen$data$AllSites, any_of(c("Site", "Park"))), by = join_by(Site)) %>%
+  group_by(Park) %>%
+  # Find out how many plots in a watershed were disturbed
+  summarise(plotsWithDisturbances = sum(!is.na(Disturbance)), 
+            numberOfPlotsInPark = n(),
+            frequency = (first(plotsWithDisturbances) * pi * 4^2)/(first(numberOfPlotsInPark) * pi * 4^2)*100)
+
+disturbancesGraph <- disturbancesFrequency %>%
+  ggplot(aes(x = Park, y = frequency)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Disturbance Frequency",
+       x = "Park", y = "Frequency") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(disturbancesGraph)
+```
+
+#### Frequency of Pests
+
+```{r}
+# Area of plots in which a disturbance occurs / total area of plots (for each park)
+
+pestFrequency <- aspen$data$Pests %>%
+  group_by(Site) %>%
+  # Group data so every site with a disturbance has only one entry
+  summarise(Pest = first(Pest)) %>%
+  # Join disturbance and all site data
+  full_join(select(aspen$data$AllSites, any_of(c("Site", "Park"))), by = join_by(Site)) %>%
+  group_by(Park) %>%
+  # Find out how many plots in a watershed were disturbed
+  summarise(plotsWithPests = sum(!is.na(Pest)), 
+            numberOfPlotsInPark = n(),
+            frequency = (first(plotsWithPests) * pi * 4^2)/(first(numberOfPlotsInPark) * pi * 4^2)*100)
+
+pestGraph <- pestFrequency %>%
+  ggplot(aes(x = Park, y = frequency)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Pest Frequency",
+       x = "Park", y = "Frequency") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(pestGraph)
+```
+
diff --git a/scripts/visualizations/aspen_visualizations_ucbn.qmd b/scripts/visualizations/aspen_visualizations_ucbn.qmd
new file mode 100644
index 0000000..19a8ada
--- /dev/null
+++ b/scripts/visualizations/aspen_visualizations_ucbn.qmd
@@ -0,0 +1,823 @@
+---
+title: "Aspen Preliminary Visualizations"
+format: 
+  html:
+    embed-resources: true
+execute: 
+  echo: false
+  message: false
+  warning: false
+---
+
+```{r, message=FALSE}
+library(tidyverse)
+library(khroma)
+library(viridis)
+library(plotly)
+aspen <- aspen::fetchAndWrangleUCBNAspen()
+```
+
+```{r}
+# find number of plots in each watershed
+
+numberOfPlots <- aspen$data$SiteVisit %>%
+  group_by(Stand) %>%
+  summarise(numberOfPlotsInStand = n_distinct(Unique_ID))
+```
+
+
+# Stand
+
+#### Average Density of Live Aspens
+
+```{r}
+# 4 meter radius of plots
+# area = pi*r^2
+# density = #/ha
+
+# this is finding the density for each plot and then averaging it for the communities
+
+aspenDensity2 <- aspen$data$Observations %>%
+  #left_join(select(aspen$data$SiteVisit, any_of(c("numberOfPlotsInStand", "globalid"))), by = c("parentglobalid" = "globalid")) %>%
+  left_join(select(numberOfPlots, any_of(c("Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POPTRE") %>%
+  # Add up LIVE POTR5 for each row
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/(pi * 4^2/10000),
+         totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Stand) %>%
+  summarise(avgDensity = mean(density),
+            sd = sd(density), 
+            numberOfPlots = n()) %>%
+  # Find values for error bars
+  mutate(ymin = avgDensity - sd,
+         ymax = avgDensity + sd)
+
+aspenDensityGraph <- aspenDensity2 %>%
+  ggplot(aes(x= Stand, y = avgDensity, label = numberOfPlots)) +
+  geom_col() +
+  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.2) +
+  theme_minimal() +
+  labs(title = "Average Aspen Density",
+       x = "Stand", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+plotly::ggplotly(aspenDensityGraph)
+```
+
+#### Density of Live Aspens Boxplot
+
+```{r}
+aspenDensity <- aspen$data$Observations %>%
+  #left_join(select(aspen$data$SiteVisit, any_of(c("numberOfPlotsInStand", "globalid"))), by = c("parentglobalid" = "globalid")) %>%
+  left_join(select(numberOfPlots, any_of(c("Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POPTRE") %>%
+  # Find density for each plot
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/((pi * 4^2)/10000)) %>%
+# Filter out communities with not enough points to graph
+group_by(Stand) %>%
+  mutate(numberOfPlots = n()) %>%
+  ungroup() %>%
+  filter(numberOfPlots >2)
+
+
+# aspenDensityGraph <- aspenDensity %>%
+#   ggplot(aes(x= Stand, y = density, fill = Stand, label = numberOfPlots)) +
+#     geom_violin() +
+#     geom_boxplot(width=0.1, color="grey", alpha=0.5) +
+#   stat_summary(fun = "mean", colour = "red", size = 2, geom = "point") +
+#   scale_fill_viridis(discrete = TRUE) +
+#   theme_minimal() +
+#   labs(title = "Aspen Density",
+#        x = "Stand", y = "Density (#/ha)") + 
+#   theme(plot.title = element_text(hjust = 0.5)) +
+#   theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0.0)) + 
+#   guides(fill="none")
+# 
+# aspenDensityGraph
+```
+
+
+```{r}
+aspenDensityGraph <- aspenDensity %>%
+  ggplot(aes(x= Stand, y = density, label = numberOfPlots)) +
+    #geom_violin() +
+  geom_boxplot(width=0.3, outlier.shape = NA) +
+  geom_point(alpha = 0.5, position=position_jitter(height=.5, width=.25)) +
+  theme_minimal() +
+  labs(title = "Aspen Density",
+       x = "Stand", y = "Density (#/ha)") +
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0.0)) + 
+  guides(fill="none")
+
+aspenDensityGraph
+```
+
+
+```{r}
+fig <- aspenDensity %>%
+  plotly::plot_ly(
+    x = ~Stand,
+    y = ~density,
+    split = ~Stand,
+    color = ~Stand,
+    colors = "Dark2",
+    type = 'violin',
+    box = list(visible = T),
+    meanline = list(visible = T)) 
+
+fig <- fig %>%
+  layout(
+    xaxis = list(
+      title = "Stand"),
+    yaxis = list(
+      title = "Density",
+      zeroline = F,
+      range = list(0, 40000)),
+    showlegend = FALSE)
+
+fig
+```
+
+
+
+#### Average Density of Conifers
+
+```{r}
+# 4 meter radius of plots
+# area = pi*r^2
+# density = #/ha
+
+# this is finding the density for each plot and then averaging it for the communities
+
+aspenDensity2 <- aspen$data$Observations %>%
+  # left_join(select(aspen$data$SiteVisit, any_of(c("numberOfPlotsInStand", "globalid"))), by = c("parentglobalid" = "globalid")) %>%
+  left_join(select(numberOfPlots, any_of(c("Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+  # filter for POTR5
+  filter(SpeciesCode != "POPTRE" & !is.na(SpeciesCode)) %>%
+  # Aggregate to get the total number of conifers at each site
+  group_by(Unique_ID, Stand) %>%
+  summarise(Class1 = sum(Class1),
+            Class2 = sum(Class2),
+            Class3 = sum(Class3),
+            Class4 = sum(Class4),
+            Class5 = sum(Class5)) %>%
+  ungroup() %>%
+    # Add up LIVE POTR5 for each row
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/(pi * 4^2/10000),
+         totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Stand) %>%
+  summarise(avgDensity = mean(density),
+            sd = sd(density), 
+            numberOfPlots = n()) %>%
+  # Find values for error bars
+  mutate(ymin = avgDensity - sd,
+         ymax = avgDensity + sd)
+
+aspenDensityGraph <- aspenDensity2 %>%
+  ggplot(aes(x= Stand, y = avgDensity, label = numberOfPlots)) +
+  geom_col() +
+  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.2) +
+  theme_minimal() +
+  labs(title = "Average Conifer Density",
+       x = "Stand", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+plotly::ggplotly(aspenDensityGraph)
+```
+
+
+
+```{r}
+#### Density of Conifers per Plot
+
+
+# coniferDensity <- aspen$data$Observations %>%
+#   #left_join(select(aspen$data$SiteVisit, any_of(c("numberOfPlotsInStand", "globalid"))), by = c("parentglobalid" = "globalid")) %>%
+#   left_join(select(numberOfPlots, any_of(c("Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+#   # filter for POTR5
+#   filter(SpeciesCode != "POPTRE" & !is.na(SpeciesCode)) %>%
+#     # Aggregate to get the total number of conifers at each site
+#   group_by(Unique_ID, Stand) %>%
+#   summarise(Class1 = sum(Class1),
+#             Class2 = sum(Class2),
+#             Class3 = sum(Class3),
+#             Class4 = sum(Class4),
+#             Class5 = sum(Class5)) %>%
+#   ungroup() %>%
+#   # Find density for each plot
+#   mutate(density = (Class1+Class2+Class3+Class4+Class5)/((pi * 4^2)/10000)) %>%
+# # Filter out communities with not enough points to graph
+# group_by(Stand) %>%
+#   mutate(numberOfPlots = n()) %>%
+#   ungroup() %>%
+#   filter(numberOfPlots >2)
+# 
+# 
+# coniferDensityGraph <- coniferDensity %>%
+#   ggplot(aes(x= Stand, y = density, fill = Stand, label = numberOfPlots)) +
+#     geom_violin() +
+#     geom_boxplot(width=0.1, color="grey", alpha=0.2) +
+#   stat_summary(fun = "mean", colour = "red", size = 2, geom = "point") +
+#   scale_fill_viridis(discrete = TRUE) +
+#   theme_minimal() +
+#   labs(title = "Conifer Density",
+#        x = "Stand", y = "Density (#/ha)") + 
+#   theme(plot.title = element_text(hjust = 0.5)) +
+#   theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0.0)) + 
+#   guides(fill="none")
+# 
+# coniferDensityGraph
+```
+
+#### Density of Aspen vs Conifers
+
+```{r}
+compareDensity <- aspen$data$Observations %>%
+  #left_join(select(aspen$data$SiteVisit, any_of(c("numberOfPlotsInStand", "globalid"))), by = c("parentglobalid" = "globalid")) %>%
+  left_join(select(numberOfPlots, any_of(c("Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+    # filter out one species that is not aspen or conifer
+  filter(!is.na(SpeciesCode)) %>%
+  mutate(isAspen = case_when(
+    SpeciesCode == "POPTRE" ~ 'Aspen',
+    .default = 'Conifer'
+  )) %>%
+    # Aggregate to get the total number of conifers at each site
+  group_by(Unique_ID, Stand, isAspen) %>%
+  summarise(Class1 = sum(Class1),
+            Class2 = sum(Class2),
+            Class3 = sum(Class3),
+            Class4 = sum(Class4),
+            Class5 = sum(Class5)) %>%
+  ungroup() %>%
+  # Add up live trees for each row
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/(pi * 4^2/10000)) %>%
+  group_by(Stand, isAspen) %>%
+  summarise(avgDensity = mean(density),
+            sd = sd(density), 
+            numberOfPlots = n()) %>%
+  # Find values for error bars
+  mutate(ymin = avgDensity - sd,
+         ymax = avgDensity + sd)
+
+
+compareDensityGraph <- compareDensity %>%
+  ggplot(aes(x= Stand, y = avgDensity, fill = as.factor(isAspen), label = numberOfPlots)) +
+  geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.2, position = position_dodge(0.9)) +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))  +
+  scale_fill_light()
+
+plotly::ggplotly(compareDensityGraph)
+```
+
+
+#### Density of Aspen vs Conifers per Plot
+
+```{r}
+compareDensity <- aspen$data$Observations %>%
+  #left_join(select(aspen$data$SiteVisit, any_of(c("numberOfPlotsInStand", "globalid"))), by = c("parentglobalid" = "globalid")) %>%
+  left_join(select(numberOfPlots, any_of(c("Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+    # filter out one species that is not aspen or conifer
+  filter(!is.na(SpeciesCode)) %>%
+  mutate(isAspen = case_when(
+    SpeciesCode == "POPTRE" ~ 'Aspen',
+    .default = 'Conifer'
+  )) %>% 
+      # Aggregate to get the total number of conifers at each site
+  group_by(Unique_ID, Stand, isAspen) %>%
+  summarise(Class1 = sum(Class1),
+            Class2 = sum(Class2),
+            Class3 = sum(Class3),
+            Class4 = sum(Class4),
+            Class5 = sum(Class5)) %>%
+  ungroup() %>%
+    # Find density for each plot
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/((pi * 4^2)/10000)) 
+# %>%
+#   filter(Community != "Curl-leaf Mountain-mahogany Shrubland & Woodland Complex")
+
+compareDensityGraph <- compareDensity %>%
+  ggplot(aes(x= Stand, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))  +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))
+
+compareDensityGraph
+```
+
+#### Density of Aspen Age Groups
+
+```{r}
+ageGroupCommunity <- aspen$data$Observations %>%
+  #left_join(select(aspen$data$SiteVisit, any_of(c("numberOfPlotsInStand", "globalid"))), by = c("parentglobalid" = "globalid")) %>%
+  left_join(select(numberOfPlots, any_of(c("Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POPTRE") %>%
+  # Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  mutate(density = treeCount/(pi * 4^2/10000)) %>%
+  # summarize POTR5 for each age group divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Stand, class) %>%
+    summarise(avgDensity = mean(density),
+            sd = sd(density), 
+            numberOfPlots = n()) %>%
+  # Find values for error bars
+  mutate(ymin = avgDensity - sd,
+         ymax = avgDensity + sd)
+
+ageGroupCommunityGraph <- ageGroupCommunity %>%
+  ggplot(aes(x= Stand, y = avgDensity, fill = factor(class))) +
+  geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.3, position = position_dodge(0.9)) +
+  theme_minimal() +
+  labs(title = "Age Class Density",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Age Class") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+plotly::ggplotly(ageGroupCommunityGraph)
+```
+
+
+#### Density of Aspen Age Groups Boxplot
+
+```{r}
+ageGroupCommunity <- aspen$data$Observations %>%
+  #left_join(select(aspen$data$SiteVisit, any_of(c("numberOfPlotsInStand", "globalid"))), by = c("parentglobalid" = "globalid")) %>%
+  left_join(select(numberOfPlots, any_of(c("Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POPTRE") %>%
+  # Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  mutate(density = treeCount/(pi * 4^2/10000)) 
+# %>%
+#   filter(Community != "Curl-leaf Mountain-mahogany Shrubland & Woodland Complex")
+
+
+ageGroupCommunityGraph <- ageGroupCommunity %>%
+  ggplot(aes(x= Stand, y = density, fill = factor(class))) +
+  geom_boxplot() +
+  theme_minimal() +
+  labs(title = "Age Class Density",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Age Class") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))
+
+ageGroupCommunityGraph
+```
+
+#### Density of Conifer Age Groups
+
+```{r}
+ageGroupCommunity <- aspen$data$Observations %>%
+  #left_join(select(aspen$data$SiteVisit, any_of(c("numberOfPlotsInStand", "globalid"))), by = c("parentglobalid" = "globalid")) %>%
+  left_join(select(numberOfPlots, any_of(c("Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+  # filter for POTR5
+  filter(SpeciesCode != "POPTRE" & !is.na(SpeciesCode)) %>%
+  # Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  mutate(density = treeCount/(pi * 4^2/10000)) 
+# %>%
+#   # summarize POTR5 for each age group divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+#   group_by(Community, class) %>%
+#     summarise(avgDensity = mean(density),
+#             sd = sd(density), 
+#             numberOfPlots = n()) %>%
+#   # Find values for error bars
+#   mutate(ymin = avgDensity - sd,
+#          ymax = avgDensity + sd)
+
+ageGroupCommunityGraph <- ageGroupCommunity %>%
+  ggplot(aes(x= Stand, y = density, fill = factor(class))) +
+  geom_boxplot() +
+  # geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  # geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.5, position = position_dodge(0.9)) +
+  theme_minimal() +
+  labs(title = "Age Class Density",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Age Class") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+ageGroupCommunityGraph
+```
+
+
+#### Density of Aspen vs Conifer Age Groups
+
+```{r}
+compareAgeDensity <- aspen$data$Observations %>%
+  #left_join(select(aspen$data$SiteVisit, any_of(c("numberOfPlotsInStand", "globalid"))), by = c("parentglobalid" = "globalid")) %>%
+  left_join(select(numberOfPlots, any_of(c("Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+    # filter out one species that is not aspen or conifer
+  filter(!is.na(SpeciesCode)) %>%
+  mutate(isAspen = case_when(
+    SpeciesCode == "POPTRE" ~ 'Aspen',
+    .default = 'Conifer'
+  )) %>%
+    # Find density for each plot
+  mutate(density = (Class1+Class2+Class3+Class4+Class5)/((pi * 4^2)/10000)) %>%  
+# Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  mutate(density = treeCount/(pi * 4^2/10000))
+```
+
+::: panel-tabset
+##### Class 1
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class1") %>%
+  ggplot(aes(x= Stand, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  scale_fill_viridis(discrete = TRUE) +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 1",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0)) +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+##### Class 2
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class2") %>%
+  ggplot(aes(x= Stand, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 2",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))  +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+
+##### Class 3
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class3") %>%
+  ggplot(aes(x= Stand, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 3",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))  +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+
+##### Class 4
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class4") %>%
+  ggplot(aes(x= Stand, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 4",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))  +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+
+##### Class 5
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class5") %>%
+  ggplot(aes(x= Stand, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 5",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))  +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+
+##### Class 6
+```{r}
+compareAgeDensityGraph <- compareAgeDensity %>%
+  filter(class == "Class6") %>%
+  ggplot(aes(x= Stand, y = density, fill = as.factor(isAspen))) +
+  geom_boxplot() +
+  #geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density Class 6",
+       x = "Stand", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()  +
+  theme(axis.text.x = element_text(size = 6, angle = 315, hjust = 0))  +
+  ylim(c(0, 20000))
+
+compareAgeDensityGraph
+```
+
+:::
+
+#### Frequency of Disturbances
+
+```{r}
+# Area of plots in which a disturbance occurs / total area of plots (for each community)
+
+disturbancesFrequency <- aspen$data$Disturbances %>%
+  group_by(Unique_ID, Stand) %>%
+  # Group data so every site with a disturbance has only one entry
+  summarise(Disturbance = first(Disturbance)) %>%
+  # Join disturbance and all site data
+  left_join(select(numberOfPlots, any_of(c("Unique_ID", "Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+  group_by(Stand) %>%
+  # Find out how many plots in a watershed were disturbed
+  summarise(plotsWithDisturbances = sum(!is.na(Disturbance)), 
+            numberOfPlotsInStand = first(numberOfPlotsInStand),
+            frequency = (first(plotsWithDisturbances) * pi * 4^2)/(first(numberOfPlotsInStand) * pi * 4^2)*100)
+
+disturbancesGraph <- disturbancesFrequency %>%
+  ggplot(aes(x = Stand, y = frequency)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Disturbance Frequency",
+       x = "Stand", y = "Frequency") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+plotly::ggplotly(disturbancesGraph)
+```
+
+
+#### Frequency of Pests
+
+```{r}
+# Area of plots in which a disturbance occurs / total area of plots (for each community)
+
+pestFrequency <- aspen$data$Pests %>%
+  group_by(Unique_ID, Stand) %>%
+  # Group data so every site with a disturbance has only one entry
+  summarise(Pest = first(Pest)) %>%
+  # Join disturbance and all site data
+  left_join(select(numberOfPlots, any_of(c("Unique_ID", "Stand", "numberOfPlotsInStand"))), by = join_by(Stand)) %>%
+  group_by(Stand) %>%
+  # Find out how many plots in a watershed were disturbed
+  summarise(plotsWithPests = sum(!is.na(Pest)), 
+            numberOfPlotsInStand = first(numberOfPlotsInStand),
+            frequency = (first(plotsWithPests) * pi * 4^2)/(first(numberOfPlotsInStand) * pi * 4^2)*100)
+
+pestGraph <- pestFrequency %>%
+  ggplot(aes(x = Stand, y = frequency)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Pest Frequency",
+       x = "Stand", y = "Frequency") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  theme(axis.text.x = element_text(size = 6, angle = 315))
+
+plotly::ggplotly(pestGraph)
+```
+
+# Park
+
+#### Density of Live Aspens
+
+```{r}
+# 4 meter radius of plots
+# area = pi*r^2
+# denisty = #/ha
+
+
+  
+aspenParkDensity <- aspen$data$Observations %>%
+  # Find number of plot in each park
+  group_by(Park) %>%
+  mutate(numberOfPlots = n_distinct(Unique_ID)) %>%
+  ungroup() %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POPTRE") %>%
+  # Add up LIVE POTR5 for each row
+  mutate(totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Park) %>%
+  summarise(density = sum(totalLiveTrees)/((first(numberOfPlots) * pi * 4^2)/10000),
+            sumLiveTrees = sum(totalLiveTrees),
+            numberOfPlots = first(numberOfPlots))
+
+
+aspenParkDensityGraph <- aspenParkDensity %>%
+  ggplot(aes(x= Park, y = density)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Aspen Density",
+       x = "Park", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(aspenParkDensityGraph)
+```
+
+#### Density of Conifers
+
+```{r}
+# 4 meter radius of plots
+# area = pi*r^2
+# denisty = #/ha
+
+
+  
+coniferParkDensity <- aspen$data$Observations %>%
+  # Find number of plot in each park
+  group_by(Park) %>%
+  mutate(numberOfPlots = n_distinct(Unique_ID)) %>%
+  ungroup() %>%
+  # filter for POTR5
+  filter(SpeciesCode != "POPTRE" & SpeciesCode != "BEUC2") %>%
+  # Add up LIVE POTR5 for each row
+  mutate(totalLiveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Park) %>%
+  summarise(density = sum(totalLiveTrees)/((first(numberOfPlots) * pi * 4^2)/10000),
+            sumLiveTrees = sum(totalLiveTrees),
+            numberOfPlots = first(numberOfPlots))
+
+
+coniferParkDensityGraph <- coniferParkDensity %>%
+  ggplot(aes(x= Park, y = density)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Conifer Density",
+       x = "Park", y = "Density (#/ha)") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(coniferParkDensityGraph)
+```
+
+
+#### Density of Aspen vs Conifers
+
+```{r}
+compareDensity <- aspen$data$Observations %>%
+    # Find number of plot in each park
+  group_by(Park) %>%
+  mutate(numberOfPlots = n_distinct(Unique_ID)) %>%
+  ungroup() %>%
+  # filter out one species that is not aspen or conifer
+  filter(SpeciesCode != "BEUC2") %>%
+  mutate(isAspen = case_when(
+    SpeciesCode == "POPTRE" ~ 'Aspen',
+    .default = 'Conifer'
+  )) %>%
+  # Add up live trees for each row
+  mutate(liveTrees = Class1+Class2+Class3+Class4+Class5) %>%
+  # summarize total POTR5 divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Park, isAspen) %>%
+  summarise(density = sum(liveTrees)/((first(numberOfPlots) * pi * 4^2)/10000),
+            totalLiveTrees = sum(liveTrees),
+            numberOfPlots = first(numberOfPlots))
+
+
+compareDensityGraph <- compareDensity %>%
+  ggplot(aes(x= Park, y = density, fill = as.factor(isAspen))) +
+  geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Aspen and Conifer Density",
+       x = "Park", y = "Density (#/ha)",
+       fill = "Species") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()
+
+plotly::ggplotly(compareDensityGraph)
+```
+
+#### Density of Aspen Age Groups
+
+```{r}
+ageGroupPark <- aspen$data$Observations %>%
+  # Find number of plot in each park
+  group_by(Park) %>%
+  mutate(numberOfPlots = n_distinct(Unique_ID)) %>%
+  ungroup() %>%
+  # filter for POTR5
+  filter(SpeciesCode == "POPTRE") %>%
+  # Pivot longer to make calculations easier
+  pivot_longer(cols = Class1:Class6, names_to = "class", values_to = "treeCount") %>%
+  # summarize POTR5 for each age group divided by area converted to ha (area = (number of plots * pi * r^2)/10000 )
+  group_by(Park, class) %>%
+  summarise(density = sum(treeCount)/((first(numberOfPlots) * pi * 4^2)/10000),
+            sumLiveTrees = sum(treeCount),
+            numberOfPlots = first(numberOfPlots)) 
+
+
+ageGroupParkGraph <- ageGroupPark %>%
+  ggplot(aes(x= Park, y = density, fill = factor(class))) +
+  geom_col(position = position_dodge(preserve = "single"), stat = "identity") +
+  theme_minimal() +
+  labs(title = "Age Class Density",
+       x = "Park", y = "Density (#/ha)",
+       fill = "Age Class") + 
+  theme(plot.title = element_text(hjust = 0.5)) +
+  scale_fill_light()
+
+plotly::ggplotly(ageGroupParkGraph)
+```
+
+
+#### Frequency of Disturbances
+
+```{r}
+# Area of plots in which a disturbance occurs / total area of plots (for each park)
+
+disturbancesFrequency <- aspen$data$Disturbances %>%
+  group_by(Park, Unique_ID) %>%
+  # Group data so every site with a disturbance has only one entry
+  summarise(Disturbance = first(Disturbance)) %>%
+  # # Join disturbance and all site data
+  # full_join(select(aspen$data$AllSites, any_of(c("Site", "Park"))), by = join_by(Site)) %>%
+  group_by(Park) %>%
+  # Find out how many plots in a watershed were disturbed
+  summarise(plotsWithDisturbances = sum(!is.na(Disturbance)), 
+            numberOfPlotsInPark = n(),
+            frequency = (first(plotsWithDisturbances) * pi * 4^2)/(first(numberOfPlotsInPark) * pi * 4^2)*100)
+
+disturbancesGraph <- disturbancesFrequency %>%
+  ggplot(aes(x = Park, y = frequency)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Disturbance Frequency",
+       x = "Park", y = "Frequency") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(disturbancesGraph)
+```
+
+#### Frequency of Pests
+
+```{r}
+# Area of plots in which a disturbance occurs / total area of plots (for each park)
+
+pestFrequency <- aspen$data$Pests %>%
+  group_by(Park, Unique_ID) %>%
+  # Group data so every site with a disturbance has only one entry
+  summarise(Pest = first(Pest)) %>%
+  # Join disturbance and all site data
+  # full_join(select(aspen$data$AllSites, any_of(c("Site", "Park"))), by = join_by(Site)) %>%
+  group_by(Park) %>%
+  # Find out how many plots in a watershed were disturbed
+  summarise(plotsWithPests = sum(!is.na(Pest)), 
+            numberOfPlotsInPark = n(),
+            frequency = (first(plotsWithPests) * pi * 4^2)/(first(numberOfPlotsInPark) * pi * 4^2)*100)
+
+pestGraph <- pestFrequency %>%
+  ggplot(aes(x = Park, y = frequency)) +
+  geom_col() +
+  theme_minimal() +
+  labs(title = "Pest Frequency",
+       x = "Park", y = "Frequency") + 
+  theme(plot.title = element_text(hjust = 0.5))
+
+plotly::ggplotly(pestGraph)
+```
+
diff --git a/scripts/visualizations/qcResults.qmd b/scripts/visualizations/qcResults.qmd
new file mode 100644
index 0000000..5932777
--- /dev/null
+++ b/scripts/visualizations/qcResults.qmd
@@ -0,0 +1,78 @@
+---
+title: "QC Results"
+format: 
+  html:
+    embed-resources: true
+execute: 
+  echo: false
+---
+
+
+```{r, message=FALSE, echo=FALSE, include=FALSE}
+library(tidyverse)
+library(aspen)
+
+# Import Data
+# 
+# data_dir = "data"
+# 
+# # Extract the path of the first Access database in the data folder folder of the project directory
+# # Or the database path can be input manually
+# #database_dir <- common::file.find(here::here("data"), "*.accdb", up = 0)[1]
+# 
+# fiveneedlepine::loadPine(database_dir)
+aspen <- fetchAndWrangleUCBNAspen()
+```
+
+# Aspen QC
+
+##### Check that there is at least one tree for every observation entry and flag entries with more than 250 trees
+```{r}
+if(nrow(aspen:::missingTreeUCBNQC(aspen)) == 0){
+  cat('No errors')
+} else{
+  rmarkdown::paged_table(aspen:::missingTreeUCBNQC(aspen))
+}
+```
+
+
+
+##### Check that no entries in the observation table have a SpeciesCode that is NA or missing
+```{r}
+if(nrow(aspen:::treeSpeciesUCBNQC(aspen)) == 0){
+  cat('No errors')
+} else{
+  rmarkdown::paged_table(aspen:::treeSpeciesUCBNQC(aspen))
+}
+```
+
+
+##### Check that unknown species entries do not have any live trees
+```{r}
+if(nrow(aspen:::unknownSpeciesUCBNQC(aspen)) == 0){
+  cat('No errors')
+} else{
+  rmarkdown::paged_table(aspen:::unknownSpeciesUCBNQC(aspen))
+}
+```
+
+
+##### Check that all plots have at least one live aspen for each visit
+```{r}
+if(nrow(aspen:::checkAspenUCBNQC(aspen)) == 0){
+  cat('No errors')
+} else{
+  rmarkdown::paged_table(aspen:::checkAspenUCBNQC(aspen))
+}
+```
+
+
+##### Flag large changes in number of tree between years
+```{r}
+# if(nrow(fiveneedlepine:::treeChangeUCBNQC(aspen)) == 0){
+#   cat('No errors')
+# } else{
+#   rmarkdown::paged_table(fiveneedlepine:::treeChangeUCBNQC(aspen))
+# }
+```
+
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..e90bb1c
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,12 @@
+# This file is part of the standard setup for testthat.
+# It is recommended that you do not modify it.
+#
+# Where should you do additional test configuration?
+# Learn more about the roles of various files in:
+# * https://r-pkgs.org/testing-design.html#sec-tests-files-overview
+# * https://testthat.r-lib.org/articles/special-files.html
+
+library(testthat)
+library(aspen)
+
+test_check("aspen")
diff --git a/tests/testthat/test-aspen_qc.R b/tests/testthat/test-aspen_qc.R
new file mode 100644
index 0000000..5cac498
--- /dev/null
+++ b/tests/testthat/test-aspen_qc.R
@@ -0,0 +1,61 @@
+
+# things to test in qc functions
+
+  # check that each data table has the correct number of columns and that the columns have the right names
+  # check that the return type from function is correct
+  # make fake data and then test number of rows?? - nah probs not
+
+
+data <- loadAndWrangleMOJNAspen()
+
+test_that("Test that missingTreeQC() works", {
+  # Run function
+  returnData <- missingTreeQC(data)
+
+  # Compare expected and actual column names
+  actual_cols <- names(returnData)
+  expected_cols <- c("Park", "Site", "VisitDate", "SpeciesCode", "Class1", "Class2", "Class3", "Class4", "Class5", "Class6", "totalTreeCount")
+  expect_equal(actual_cols, expected_cols)
+
+  # Check that function returns a database
+  expect_s3_class(returnData, "data.frame")
+})
+
+test_that("Test that treeSpeciesQC() works", {
+  # Run function
+  returnData <- treeSpeciesQC(data)
+
+  # Compare expected and actual column names
+  actual_cols <- names(returnData)
+  expected_cols <- c("Park", "Site", "VisitDate", "SpeciesCode", "Class1", "Class2", "Class3", "Class4", "Class5", "Class6")
+  expect_equal(actual_cols, expected_cols)
+
+  # Check that function returns a database
+  expect_s3_class(returnData, "data.frame")
+})
+
+test_that("Test that unknownSpeciesQC() works", {
+  # Run function
+  returnData <- unknownSpeciesQC(data)
+
+  # Compare expected and actual column names
+  actual_cols <- names(returnData)
+  expected_cols <- c("Park", "Site", "VisitDate", "SpeciesCode", "Class1", "Class2", "Class3", "Class4", "Class5", "Class6")
+  expect_equal(actual_cols, expected_cols)
+
+  # Check that function returns a database
+  expect_s3_class(returnData, "data.frame")
+})
+
+test_that("Test that checkAspenQC() works", {
+  # Run function
+  returnData <- checkAspenQC(data)
+
+  # Compare expected and actual column names
+  actual_cols <- names(returnData)
+  expected_cols <- c("Park", "Site","VisitType", "VisitDate", "SpeciesCode")
+  expect_equal(actual_cols, expected_cols)
+
+  # Check that function returns a database
+  expect_s3_class(returnData, "data.frame")
+})
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
new file mode 100644
index 0000000..fe8cbc2
--- /dev/null
+++ b/tests/testthat/test-utils.R
@@ -0,0 +1,25 @@
+
+
+# Things to test in fetchAndWrangleAspen - only function in utils rn
+  # check that $data and $meta data exists
+  # check that there are the correct number of tables in data and metadata
+  # check that each data table has the correct number of columns and that the columns have the right names
+  # check that the return type from function is correct
+
+data <- loadAndWrangleMOJNAspen()
+
+test_that("Test that loadAndWrangleMOJNAspen() works", {
+  # Compare expected and actual names in return object
+  actual_names <- names(data)
+  expected_names <- c("data", "metadata")
+  expect_equal(actual_names, expected_names)
+
+  # Compare expected and actual names of the aspen data frames
+  actual_cols <- names(data$data)
+  expected_cols <- c("AllSites", "SiteVisit", "Disturbances", "Observations", "Pests")
+  expect_equal(actual_cols, expected_cols)
+
+  # Check that the first object in the data list is a database
+  returnType <- data$data[[1]]
+  expect_s3_class(returnType, "data.frame")
+})