========

SETUP

========

FILEPATH:

filepath = "/Users/emvanmetre/Desktop/Arctic/X-Y_Plots/data/processed/Combined/"

PICK DATE RANGES:

years = c("2022", "2023")
seasons = c("Spring", "Summer","Fall")
sites = c("TNHA", "BEO")

# put season values for the season that has the start of the data
date_start = "2022-06-01" # data starts in June 2022 (YEAR-MO-DY) where day is always 01
# put season values for the season that has the last of the data
date_end = "2023-11-01" # data ends after November of 2023 (will get data up UNTIL date_end not after)

PICK OUTPUT:

scree = F # scree plot
eigen = T # eigenvectors and eigenvalues

Dates (for selecting date ranges; no need to edit)

spring_months = c("March", "April", "May")
summer_months = c("June","July","August")
fall_months = c("September", "October", "November")
winter_months = c("December", "January", "February")

spring_dates = data.frame(months=spring_months, start=c("-03-01","-04-01","-05-01"), end=c("-04-01","-05-01","-06-01"))

summer_dates = data.frame(months=summer_months, start=c("-06-01","-07-01","-08-01"), end=c("-07-01","-08-01","-09-01"))

fall_dates = data.frame(months=fall_months, start=c("-09-01","-10-01","-11-01"), end=c("-10-01","-11-01","-12-01"))

winter_dates = data.frame(months=winter_months, start=c("-12-01","-01-01","-02-01"), end=c("-01-01","-02-01","-03-01"))

all_dates = data.frame(matrix(nrow = 0, ncol = 4))
for(yur in years){
  # spring
  if("Spring" %in% seasons){
    for(i in c(1:3)){
      curdate = as.Date(paste0(yur, spring_dates[i, 2]))
      if(curdate >= date_start && curdate < date_end){
        all_dates <- rbind(all_dates, c(spring_dates[i,1], paste0(yur, spring_dates[i, 2]), paste0(yur, spring_dates[i, 3]),"Spring"))
      }
    }
  }
  if("Summer" %in% seasons){
    for(i in c(1:3)){
      curdate = as.Date(paste0(yur, summer_dates[i, 2]))
      if(curdate >= date_start && curdate < date_end){
        all_dates <- rbind(all_dates, c(summer_dates[i,1], paste0(yur, summer_dates[i, 2]), paste0(yur, summer_dates[i, 3]),"Summer"))
      }
    }
  }
  if("Fall" %in% seasons){
    for(i in c(1:3)){
      curdate = as.Date(paste0(yur, fall_dates[i, 2]))
      if(curdate >= date_start && curdate < date_end){
        all_dates <- rbind(all_dates, c(fall_dates[i,1], paste0(yur, fall_dates[i, 2]), paste0(yur, fall_dates[i, 3]),"Fall"))
      }
    }
  }
  if("Winter" %in% seasons){
    for(i in c(1:3)){
      curdate = as.Date(paste0(yur, winter_dates[i, 2]))
      if(curdate >= date_start && curdate < date_end){
        all_dates <- rbind(all_dates, c(winter_dates[i,1], paste0(yur, winter_dates[i, 2]), paste0(yur, winter_dates[i, 3]),"Winter"))
      }
    }
  }
  
}
colnames(all_dates) = c("months","start","end","szn")

Packages

library(dplyr)
library(lubridate)
library(tidyverse)

=========================

DEFINE FUNCTIONS

=========================

Temporal Range: Season Vertical Spatial Range: 30-45 cm Horizontal Spatial Range: stations across site (TNHA, SSMH, BEO) –> Average Total Site –> North vs South (except for BEO)

Filter by Site and Join Tables

pick_site <- function(cursite){
  gtfile = paste0("GroundTemperature_",szn,yr,"_DAILY.csv")
  airfile = paste0("AirTemperature_",szn,yr,"_DAILY.csv")
  vwcfile = paste0("VWC_",szn,yr,"_DAILY.csv")
  solfile = paste0("Solar_",szn,yr,"_DAILY.csv")
  windspeed = paste0("WindSpeed_",szn,yr,"_DAILY.csv")
  winddir = paste0("WindDirection_",szn,yr,"_DAILY.csv")
  
  grndtmp <<- read.csv(paste0(filepath, gtfile))
  airtmp <<- read.csv(paste0(filepath, airfile))
  vwc <<- read.csv(paste0(filepath, vwcfile))
  solar <<- read.csv(paste0(filepath, solfile))
  wind <<- read.csv(paste0(filepath, windspeed))
  
  pca_ground = grndtmp %>% filter(site == cursite) %>% filter(depth == "30cm")
  pca_air = airtmp %>% filter(site == cursite)
  pca_wind = wind %>% filter(site == cursite)
  pca_solar = solar %>% filter(site == cursite)
  pca_vwc = vwc %>% filter(site==cursite) %>% filter(depth == "30-45cm")
  
  big_df <<- full_join(pca_ground, pca_air, by=c("Date", "site", "station", "fullname")) %>% 
    select(Date, site, station, fullname, avg.x, avg.y) %>% rename("groundtemp"="avg.x", "airtemp"="avg.y")
  big_df <<- full_join(big_df, pca_wind, by=c("Date", "site", "station", "fullname")) %>% 
    select(Date, site, station, fullname, groundtemp, airtemp, avg) %>% rename("windspeed" = avg)
  big_df <<- full_join(big_df, pca_solar, by=c("Date", "site", "station", "fullname")) %>% 
    select(Date, site, station, fullname, groundtemp, airtemp, windspeed, avg)  %>% rename("solar" = avg)
  big_df <<- full_join(big_df, pca_vwc, by=c("Date", "site", "station", "fullname")) %>% 
    select(Date, site, station, fullname, groundtemp, airtemp, windspeed, solar, avg)  %>% rename("vwc" = avg)
  
  if(szn == "Winter"){
    big_df <<- big_df %>% select(-solar)
  }
  return(big_df)
}

Filter by Date Range

pick_dates <- function(datemin, datemax, big_df){
  pca_df <<- big_df %>% filter(Date >= datemin) %>% filter(Date < datemax)
  
  # get rid of NAs
  pca_df <<- na.omit(pca_df)
  pca_df <<- unique(pca_df)
  return(pca_df)
}

Calculate PCA

calc_pca <- function(pca_df){
  pca <<- prcomp(pca_df[,5:ncol(pca_df)], center=TRUE, scale.=TRUE)

  #take out variables
  sd <- pca$sdev
  loads <<- pca$rotation
  rownames(loads) <<- colnames(pca_df[5:ncol(pca_df)])
  scores <<- pca$x
  
  var <- sd^2
  varPercent <- var/sum(var) * 100
  
  return(list("pca"=pca, "loads"=loads))
}

Make Scree Plot

make_scree <- function(pca){
  sd <- pca$sdev
  
  var <- sd^2
  varPercent <- var/sum(var) * 100
  
  barplot(varPercent, xlab="PC", ylab="Percent Variance", names.arg=1:length(varPercent), 
          las=1, ylim=c(0, max(varPercent)), col="gray")
  abline(h=1/ncol(pca_df[5:ncol(pca_df)])*100, col="red")
}

Display Eigenvectors and Eigenvalues

make_eigen <- function(pca){
  eigenvectors <- pca$rotation
  print("Eigenvectors (Loadings):")
  print(eigenvectors)
  
  print("Loadings Cutoff:")
  sqrt(1/ncol(pca_df[5:ncol(pca_df)])) # cutoff for "important" loadings
  
  # Access the eigenvalues (variances of the principal components)
  eigenvalues <- (pca$sdev)^2
  print("Eigenvalues:")
  print(eigenvalues)
}

===============

PCA PLOTS