### Figure3.R
### PCA and histograms of thicket-wide plot soil parameters showing the three
### sampled soils relative to the available data.

### Script for:
# "Restoring South African subtropical succulent thicket using 
# Portulacaria afra: rooting variation across three soil types"
# by AJ Potts, D Liddell, C Clarke, N Galuszynski

### Libraries
library(groundhog)
set.groundhog.folder("C:\\Groundhog_R\\R4.4.1_2024-05-14") # Sets the folder where Groundhog will store library versions
GroundhogDay <- '2024-05-14' # Specifies the Groundhog snapshot date for reproducibility

groundhog.library("tidyverse", GroundhogDay) # Loads the Tidyverse package for data manipulation and visualization
groundhog.library("vegan", GroundhogDay) # Loads the Vegan package for ecological data analysis

# Set your working directory here
# setwd(")
setwd("D:\\Google Drive\\0_SpekboomResearchGroup2.0\\3_Experiments\\NMU22(07)SoilRootingExperiment\\DataForSubmission/") # Sets the working directory for the project

read_csv("twp_soils.csv") -> dat # Reads the soil data from a CSV file into a dataframe

### FIGURE 3 ###
dat %>%
  select(-ID, -LabNumber, -Soil, -Classification, -Comment, -TValue, -Silt_perc, -Sand_perc, -K_perc, -Na_perc, 
         -Ca_perc, -Mg_perc, -N_perc, -Org_C_perc, -StoneVolume_perc) %>% # Removes columns not required for Figure 3
  pivot_longer(cols = 2:ncol(.), names_to = "variable", values_to = "value") %>% # Converts data to long format
  mutate(variable = factor(variable,
                           levels = c("pH_KCl", "C_perc", "Clay_perc", "PBrayII", "EC_mSm", "K_mgkg", "K_plus", "Ca_2plus", "Mg_2plus", "Na_plus"),
                           labels = c("pH (KCl)", "Total C (%)", "Clay (%)", "P Bray II (mg/kg)", "EC (mS.m; 1:5)", "K (mg/kg)", 
                                      "K + (cmol(+)/kg)", "Ca 2+ (cmol(+)/kg)", "Mg 2+ (cmol(+)/kg)", "Na + (cmol(+)/kg)"))) ->
  dat2

dat2 %>%
  filter(Plot %in% c("SiteA (soil1)", "SiteB (soil3)", "SiteC (soil2)")) %>% # Filters data for focal sites
  mutate(Plot = factor(Plot, level = c("SiteA (soil1)", "SiteB (soil3)", "SiteC (soil2)"),
                       label = c("Site A", "Site B", "Site C"))) ->
  dat2_focalsoils

# Creates histogram faceted by variable
dat2 %>%
  filter(!(variable == "EC (mS.m; 1:5)" & value > 100)) %>% # Removes outlier EC values
  ggplot(aes(value)) +
  geom_histogram() +
  geom_segment(data = dat2_focalsoils, aes(x = value, y = 20, xend = value, yend = 2), 
               arrow = arrow(length = unit(0.5, 'cm')), linewidth = 1.9,colour="black") +
  geom_segment(data = dat2_focalsoils, aes(x = value, y = 20, xend = value, yend = 2), 
               arrow = arrow(length = unit(0.5, 'cm')), linewidth = 1.5,colour="white") +
  geom_segment(data = dat2_focalsoils, aes(x = value, y = 20, xend = value, yend = 2, colour = Plot, linetype = Plot), 
               arrow = arrow(length = unit(0.5, 'cm')), linewidth = 1) +
  scale_colour_manual(values = c("Site A" = "grey", "Site B" = "darkgrey", "Site C" = "black")) +
  facet_wrap(~variable, scales = "free", ncol = 2) +
  theme_bw() +
  theme(legend.position = "top", legend.title = element_blank()) +
  xlab("") + ylab("count")
ggsave("Figure3.jpg", units = "cm", height = 21, width = 15) # Saves the plot as a JPG