install.packages(“stoichUtilties”)

install.packages(“tidyverse”)
install.packages(“lubridate”)
install.packages(“sf”)
install.packages(“units”)

# Tidyverse and other support libraries
library(tidyverse)
library(lubridate)
library(sf)
library(units)

# Load stoichUtilities
library(stoichUtilities)

# Load the STOICH data (assumed to be in a directory named "data" inside the home directory)
stoichData <- loadSTOICH(dataPath=file.path(path.expand("~"), "data"))

# Filter the STOICH data
stoichFiltered <- filterSTOICH(dataTables=stoichData, 
                               var="State", 
                               val="North Dakota", 
                               condition="Equal")

# Match organism stoichiometry data with water chemsitry data for samples that weren't taken at the exact same time
stoichPaired <- locateDataPairsSTOICH(stoichData, 
                                      timeDiff=2, 
                                      timeUnits="weeks", 
                                      distance=2, 
                                      pairMethod="Min Time", 
                                      ignoreExisting=TRUE)

# Join all the tables into one large wide table
stoichTable <- joinSTOICH(stoichPaired)

# Count how many sites are in each state.
stateCounts <- stoichData$tbl_Site |>
  select(c("State")) |> 
  drop_na(State) |>
  # Change the data to lowercase to allow the table to merge with the Tidyverse state data
  mutate(State=str_to_lower(State)) |>
  group_by(State) |>
  # Summarise to get the number of sites for each state.
  summarise(nSites=n(), .groups="keep")

# Use state data from Tidyverse and join the site data to it.
us_states <- map_data("state") |>
  left_join(stateCounts, by=c("region"="State"))

# Create a list of states in the Tidyverse dataset so only data from states with boundaries are used.
stateList <- unique(us_states$region)

# Split out the GPS data to map individual sites.
gpsSites <- stoichData$tbl_Site |>
  filter(str_to_lower(State) %in% stateList) |>
  select(c("Latitude", "Longitude")) |>
  drop_na(Latitude, Longitude) |>
  # Round the gps coordinates so points that are close will be merged into a single point.
  mutate(across(everything(), ~round(.x, 2))) |>
  unique()

# Start building the map.
ggplot(data=us_states, mapping=aes(x=long, y=lat, group=group, fill=nSites)) +
  geom_polygon(color="black") +
  theme(legend.position = "bottom") +
  geom_point(data=gpsSites, aes(x=Longitude, y=Latitude, group=NA), size=2, shape=23, fill="orange")

help(loadSTOICH)