#Original code written by Jenna Messick, revised by John Unterschuetz and Abby Moore

library(ggplot2) #for all graphics
library(gridExtra) #for arranging plots on a page
library(maps) #for making maps
library(mapdata) #for the state and county boundaries
library(ggpubr) #for writing equations in plots
library(ggspatial) #for north arrow and scale bar on maps
library(sf) #for making ggplot actually recognize that we are dealing with maps and not random shapes
library(dplyr) #to filter data
library(tidyverse) #to add temperature data from a different spreadsheet into the sp_data matrix--not sure if I actually need this

rm(list=ls())

#Change the species name here.
#sp_name <- "Liatris_punctata"
#sp_name <- "Ratibida_columnifera"
#sp_name <- "Vernonia_baldwinii"
sp_name <- "Grindelia_ciliata"
#Change the folder for the input and output files here.
IOFolder <- "C:/myfiles/botany/lab/John/phenology_paper/analyses/"
#time_period <- "annual"
#time_period <- "winter"
#time_period <- "spring"
time_period <- "summer"
#time_period <- "fall"

#It is expecting a tab-delimitted text file in the IOFolder with a name of the form [species name]_data.txt, so Liatris_punctata_data.txt
#It expects that there will be columns labeled long for longitude, lat for latitude, sname for species name, Dayofyear for Julian day of year,
#and PhenoPhase for phenology score (0 to 4).  If you want to check the counties (highly recommended!), then the county needs to be in
#a column labeled County.  All additional columns will be ignored.
#The latitiude and longitude must not have the degree symbol, and west must be negative, if plotting within the US.
#(If not plotting within the US, coordinates must be changed according to the base map used.)
#Phenophases are as follows:
#0: not yet flowering
#1: First Flowers
#2: Peak Flowering
#3: Last Flowers
#4: Flowering Ended

#reading in the data
InFileName <- paste(IOFolder, sp_name, "_data.txt", sep = '')
#data for the species
sp_data <- read.table(paste(IOFolder, sp_name, "_data.txt", sep = ''), sep="\t", header = TRUE)
#the temperature data file
all_temps <- read.table(paste(IOFolder, "OK_", time_period, "_temp_all.txt", sep = ''), sep="\t", header = TRUE)
#the file showing which climate division each county is in
div_data <- read.table(paste(IOFolder, "OK_climate_divisions.txt", sep = ''), sep="\t", header = TRUE)
#the file showing the annual temperatures for each climate division
div_temps <- read.table(paste(IOFolder, "OK_", time_period, "_temp_div.txt", sep = ''), sep="\t", header = TRUE)


#A couple of checks
head(sp_data)
#There should not be any lines marked NA at the end of the file.
#If there are, you have blank lines and you need to fix your text file.
sp_data$PhenoPhase

#removing all lines with PhenoPhase as NA, after this check
sp_data <- filter(sp_data, PhenoPhase != "NA")
#checking to make sure that worked
sp_data$PhenoPhase
#making PhenoPhase into a factor
sp_data$PhenoPhase <- as.factor(sp_data$PhenoPhase)
glimpse(sp_data)

#also checking the temperature data
head(all_temps)
glimpse(all_temps)

#and the climate division table
head(div_data)
glimpse(div_data)

#and the temperature data per climate division
head(div_temps)
glimpse(div_temps)

#adding the climate division to the sp_data
sp_data <- sp_data %>% left_join(div_data, by = "County")
glimpse(sp_data)

#now adding temperature from all_temps to the sp_data data matrix
#This joins temp_data to sp_data with the row of temp_data for tha correct year placed with sp_data
sp_data <- sp_data %>% left_join(all_temps, by = "Year")
glimpse(sp_data)

#now adding the climate-division specific temperature from div_temps to the sp_data data matrix
#This time we need to join not just by Year, but also by Div_Num.
sp_data <- sp_data %>% left_join(div_temps, by = c("Year", "Div_Num"))
glimpse(sp_data)

#subsetting data:
Pre <- subset(sp_data, PhenoPhase == "0")
First <- subset(sp_data, PhenoPhase == "1")
Peak <- subset(sp_data, PhenoPhase == "2")
Last <- subset(sp_data, PhenoPhase == "3")
Ended <- subset(sp_data, PhenoPhase == "4")



#####*****Making the scatter plots with regression lines*****

#figuring out the overall first and last days so that the axes can be the same in all plots
DayMin <- min(sp_data$Dayofyear)
DayMin
DayMax <- max(sp_data$Dayofyear)
DayMax

YearMin <- min(sp_data$Year)
YearMin
YearMax <- max(sp_data$Year)
YearMax

TempMin <- min(sp_data$Temp_all)
TempMin
TempMax <- max(sp_data$Temp_all)
TempMax

#now for minimum and maximum temperatures across all climate divisions
TempMinDiv <- min(sp_data$Div_Temp, na.rm = TRUE)
TempMinDiv
TempMaxDiv <- max(sp_data$Div_Temp, na.rm = TRUE)
TempMaxDiv


####Scatter plots of Dayofyear ~Year###
#graph for first flowers
model_First <- lm(Dayofyear ~ Year, data = First)
summary(model_First)

#various ways to add equation: https://stackoverflow.com/questions/7549694/add-regression-line-equation-and-r2-on-graph, https://groups.google.com/forum/#!topic/ggplot2/1TgH-kG5XMA
#However, I can't change the number of significant digits for the equation for the line, just for R^2 and p.
#I am using this option for r^2 and p, though.
#Another option, the one I used for the equation.
#https://www.roelpeters.be/how-to-add-a-regression-equation-and-r-squared-in-ggplot2/

#version of equation with r^2 in equation
#eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
#                   list(a = format(unname(coef(model_First)[1]), digits = 3),
#                        b = format(unname(coef(model_First)[2]), digits = 3),
#                        r2 = format(summary(model_First)$r.squared, digits = 3)))

#making the equation to add to the plot:
#correct equation if the slope is positive:
if (coef(model_First)[2] > 0) {
  eq_First <- substitute(italic(y) == a + b %.% italic(x), 
                         list(a = format(unname(coef(model_First)[1]), digits = 3), b = format(unname(coef(model_First)[2]), digits = 3)))
}
#correct equation if the slope is negative:
if (coef(model_First)[2] < 0) {
  eq_First <- substitute(italic(y) == a - b %.% italic(x), 
                         list(a = format(unname(coef(model_First)[1]), digits = 3), b = format(unname(-coef(model_First)[2]), digits = 3)))
}

Scatter_First <- ggplot(First, aes(x=Year, y = Dayofyear)) + geom_point() + ylim(DayMin,DayMax) +
  labs(title = "First Flowers", y = "Day of Year") + geom_smooth(method = 'lm', se = FALSE, fullrange = TRUE) +
  scale_x_continuous(limits = c(1900, 2020), breaks = seq(1900, 2020, 20)) +
  #scale_y_continuous can be specified for the individual plants, if desired
  #scale_y_continuous(limits = c(150, 325), breaks = seq(150, 325, 25)) +
  #if scale_y_continuous is specified, the positions of the r2 and p value as well as the equation for the regression line also need to be changed.
  #stat_cor(label.y = 325, digits = 3, aes(label = paste(..rr.label..))) +
  #stat_regline_equation(label.y = 315) + aes(label = as.character(as.expression(eq_First))) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_First))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_First

#Now repeating all of this for the other PhenoPhases, but without the comments.

#graph for pre-flowering
model_Pre <- lm(Dayofyear ~ Year, data = Pre)
if (coef(model_Pre)[2] > 0) {
  eq_Pre <- substitute(italic(y) == a + b %.% italic(x), 
                       list(a = format(unname(coef(model_Pre)[1]), digits = 3), b = format(unname(coef(model_Pre)[2]), digits = 3)))
}
if (coef(model_Pre)[2] < 0) {
  eq_Pre <- substitute(italic(y) == a - b %.% italic(x), 
                       list(a = format(unname(coef(model_Pre)[1]), digits = 3), b = format(unname(-coef(model_Pre)[2]), digits = 3)))
}
Scatter_Pre <- ggplot(Pre, aes(x=Year, y = Dayofyear)) + geom_point() +  xlim(YearMin,YearMax) + ylim(DayMin,DayMax) +
  labs(title = "Pre-Flowering", y = "Day of Year") + geom_smooth(method = 'lm', se = FALSE, fullrange = TRUE) +
  scale_x_continuous(limits = c(1900, 2020), breaks = seq(1900, 2020, 20)) +
  stat_cor(label.y = (DayMax), digits = 3, aes(label = paste(..rr.label..,))) +
  stat_regline_equation(label.y = (DayMax - 10)) + aes(label = as.character(as.expression(eq_Pre))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Pre

#graph for peak flowering
model_Peak <- lm(Dayofyear ~ Year, data = Peak)
if (coef(model_Peak)[2] > 0) {
  eq_Peak <- substitute(italic(y) == a + b %.% italic(x), 
                        list(a = format(unname(coef(model_Peak)[1]), digits = 3), b = format(unname(coef(model_Peak)[2]), digits = 3)))
}
if (coef(model_Peak)[2] < 0) {
  eq_Peak <- substitute(italic(y) == a - b %.% italic(x), 
                        list(a = format(unname(coef(model_Peak)[1]), digits = 3), b = format(unname(-coef(model_Peak)[2]), digits = 3)))
}
Scatter_Peak <- ggplot(Peak, aes(x=Year, y = Dayofyear)) + geom_point() +  xlim(YearMin,YearMax) + ylim(DayMin,DayMax) +
  labs(title = "Peak Flowering", y = "Day of Year") + geom_smooth(method = 'lm', se = FALSE, fullrange = TRUE) +
  scale_x_continuous(limits = c(1900, 2020), breaks = seq(1900, 2020, 20)) +
  #scale_y_continuous(limits = c(150, 325), breaks = seq(150, 325, 25)) +
  stat_cor(label.y = (DayMax), digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  #stat_cor(label.y = (325), digits = 3, aes(label = paste(..rr.label..))) +
  stat_regline_equation(label.y = (DayMax - 10)) + aes(label = as.character(as.expression(eq_Peak))) +
  #stat_regline_equation(label.y = (315)) + aes(label = as.character(as.expression(eq_Peak))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Peak

#graph for last flowers
model_Last <- lm(Dayofyear ~ Year, data = Last)
if (coef(model_Last)[2] > 0) {
  eq_Last <- substitute(italic(y) == a + b %.% italic(x), 
                        list(a = format(unname(coef(model_Last)[1]), digits = 3), b = format(unname(coef(model_Last)[2]), digits = 3)))
}
if (coef(model_Last)[2] < 0) {
  eq_Last <- substitute(italic(y) == a - b %.% italic(x), 
                        list(a = format(unname(coef(model_Last)[1]), digits = 3), b = format(unname(-coef(model_Last)[2]), digits = 3)))
}
Scatter_Last <- ggplot(Last, aes(x=Year, y = Dayofyear)) + geom_point() +  xlim(YearMin,YearMax) + ylim(DayMin,DayMax) +
  labs(title = "Last Flowers", y = "Day of Year") + geom_smooth(method = 'lm', se = FALSE, fullrange = TRUE) +
  scale_x_continuous(limits = c(1900, 2020), breaks = seq(1900, 2020, 20)) +
  #scale_y_continuous(limits = c(150, 325), breaks = seq(150, 325, 25)) +
  stat_cor(label.y = (DayMax), digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  #stat_cor(label.y = (325), digits = 3, aes(label = paste(..rr.label..))) +
  stat_regline_equation(label.y = (DayMax - 10)) + aes(label = as.character(as.expression(eq_Last))) +
  #stat_regline_equation(label.y = (315)) + aes(label = as.character(as.expression(eq_Last))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Last

#graph for flowering ended
model_Ended <- lm(Dayofyear ~ Year, data = Ended)
if (coef(model_Ended)[2] > 0) {
  eq_Ended <- substitute(italic(y) == a + b %.% italic(x), 
                         list(a = format(unname(coef(model_Ended)[1]), digits = 3), b = format(unname(coef(model_Ended)[2]), digits = 3)))
}
if (coef(model_Ended)[2] < 0) {
  eq_Ended <- substitute(italic(y) == a - b %.% italic(x), 
                         list(a = format(unname(coef(model_Ended)[1]), digits = 3), b = format(unname(-coef(model_Ended)[2]), digits = 3)))
}
Scatter_Ended <- ggplot(Ended, aes(x=Year, y = Dayofyear)) + geom_point() +  xlim(YearMin,YearMax) + ylim(DayMin,DayMax) +
  labs(title = "Flowering Ended", y = "Day of Year") + geom_smooth(method = 'lm', se = FALSE, fullrange = TRUE) +
  scale_x_continuous(limits = c(1900, 2020), breaks = seq(1900, 2020, 20)) +
  #scale_y_continuous(limits = c(125, 300), breaks = seq(125, 300, 25)) +
  stat_cor(label.y = (DayMax), digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  #stat_cor(label.y = (300), digits = 3, aes(label = paste(..rr.label..))) +
  stat_regline_equation(label.y = (DayMax - 10)) + aes(label = as.character(as.expression(eq_Ended))) +
  #stat_regline_equation(label.y = (290)) + aes(label = as.character(as.expression(eq_Ended))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Ended

#Print the plot with the selected graphs

ScatterPlotFN <- paste(IOFolder, sp_name, "_ScatterPlot_Year.pdf", sep = '')

pdf(file=ScatterPlotFN, width = 10, height = 8)
if (dim(Ended)[1] <= 2) {
  grid.arrange(Scatter_First, Scatter_Peak, Scatter_Last, ncol = 2, nrow = 2)
}
if (dim(Ended)[1] > 2) {
  grid.arrange(Scatter_First, Scatter_Peak, Scatter_Last, Scatter_Ended, ncol = 2, nrow = 2)
}

dev.off()

#####Scatter plots of Dayofyear ~ Temp_all###--Exact same as for year
#first flowers, no comments, see comments for the previous scatter plots (versus year, instead of versus temp)
model_First <- lm(Dayofyear ~ Temp_all, data = First)
summary(model_First)
if (coef(model_First)[2] > 0) {
  eq_First <- substitute(italic(y) == a + b %.% italic(x), 
                         list(a = format(unname(coef(model_First)[1]), digits = 3), b = format(unname(coef(model_First)[2]), digits = 3)))
}
if (coef(model_First)[2] < 0) {
  eq_First <- substitute(italic(y) == a - b %.% italic(x), 
                         list(a = format(unname(coef(model_First)[1]), digits = 3), b = format(unname(-coef(model_First)[2]), digits = 3)))
}
Scatter_First <- ggplot(First, aes(x=Temp_all, y = Dayofyear)) + geom_point() + xlim(TempMin,TempMax) + ylim(DayMin,DayMax) +
  labs(title = "First Flowers", y = "Day of Year", x = "Yearly Mean Temperature, Statewide") + 
  geom_smooth(method = 'lm', se = FALSE) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_First))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_First

#graph for pre-flowering
model_Pre <- lm(Dayofyear ~ Temp_all, data = Pre)
summary(model_Pre)
if (coef(model_Pre)[2] > 0) {
  eq_Pre <- substitute(italic(y) == a + b %.% italic(x), 
                       list(a = format(unname(coef(model_Pre)[1]), digits = 3), b = format(unname(coef(model_Pre)[2]), digits = 3)))
}
if (coef(model_Pre)[2] < 0) {
  eq_Pre <- substitute(italic(y) == a - b %.% italic(x), 
                       list(a = format(unname(coef(model_Pre)[1]), digits = 3), b = format(unname(-coef(model_Pre)[2]), digits = 3)))
}
Scatter_Pre <- ggplot(Pre, aes(x=Temp_all, y = Dayofyear)) + geom_point() + xlim(TempMin,TempMax) + ylim(DayMin,DayMax) +
  labs(title = "Pre-Flowering", y = "Day of Year", x = "Yearly Mean Temperature, Statewide") + 
  geom_smooth(method = 'lm', se = FALSE) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_Pre))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Pre

#graph for peak flowering
model_Peak <- lm(Dayofyear ~ Temp_all, data = Peak)
summary(model_Peak)
if (coef(model_Peak)[2] > 0) {
  eq_Peak <- substitute(italic(y) == a + b %.% italic(x), 
                        list(a = format(unname(coef(model_Peak)[1]), digits = 3), b = format(unname(coef(model_Peak)[2]), digits = 3)))
}
if (coef(model_Peak)[2] < 0) {
  eq_Peak <- substitute(italic(y) == a - b %.% italic(x), 
                        list(a = format(unname(coef(model_Peak)[1]), digits = 3), b = format(unname(-coef(model_Peak)[2]), digits = 3)))
}
Scatter_Peak <- ggplot(Peak, aes(x=Temp_all, y = Dayofyear)) + geom_point() + xlim(TempMin,TempMax) + ylim(DayMin,DayMax) +
  labs(title = "Peak Flowering", y = "Day of Year", x = "Yearly Mean Temperature, Statewide") + 
  geom_smooth(method = 'lm', se = FALSE) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_Peak))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Peak

#graph for last flowers
model_Last <- lm(Dayofyear ~ Temp_all, data = Last)
summary(model_Last)
if (coef(model_Last)[2] > 0) {
  eq_Last <- substitute(italic(y) == a + b %.% italic(x), 
                        list(a = format(unname(coef(model_Last)[1]), digits = 3), b = format(unname(coef(model_Last)[2]), digits = 3)))
}
if (coef(model_Last)[2] < 0) {
  eq_Last <- substitute(italic(y) == a - b %.% italic(x), 
                        list(a = format(unname(coef(model_Last)[1]), digits = 3), b = format(unname(-coef(model_Last)[2]), digits = 3)))
}
Scatter_Last <- ggplot(Last, aes(x=Temp_all, y = Dayofyear)) + geom_point() + xlim(TempMin,TempMax) + ylim(DayMin,DayMax) +
  labs(title = "Last Flowers", y = "Day of Year", x = "Yearly Mean Temperature, Statewide") + 
  geom_smooth(method = 'lm', se = FALSE) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_Last))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Last

#graph for flowering ended
model_Ended <- lm(Dayofyear ~ Temp_all, data = Ended)
summary(model_Ended)
if (coef(model_Ended)[2] > 0) {
  eq_Ended <- substitute(italic(y) == a + b %.% italic(x), 
                         list(a = format(unname(coef(model_Ended)[1]), digits = 3), b = format(unname(coef(model_Ended)[2]), digits = 3)))
}
if (coef(model_Ended)[2] < 0) {
  eq_Ended <- substitute(italic(y) == a - b %.% italic(x), 
                         list(a = format(unname(coef(model_Ended)[1]), digits = 3), b = format(unname(-coef(model_Ended)[2]), digits = 3)))
}
Scatter_Ended <- ggplot(Ended, aes(x=Temp_all, y = Dayofyear)) + geom_point() + xlim(TempMin,TempMax) + ylim(DayMin,DayMax) +
  labs(title = "Flowering Ended", y = "Day of Year", x = "Yearly Mean Temperature, Statewide") + 
  geom_smooth(method = 'lm', se = FALSE) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_Ended))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Ended

#Print the plot with the selected graphs

ScatterPlotFN <- paste(IOFolder, sp_name, "_", time_period, "_ScatterPlot_Temp_All.pdf", sep = '')

pdf(file=ScatterPlotFN, width = 10, height = 8)
if (dim(Ended)[1] <= 2) {
  grid.arrange(Scatter_First, Scatter_Peak, Scatter_Last, ncol = 2, nrow = 2)
}
if (dim(Ended)[1] > 2) {
  grid.arrange(Scatter_First, Scatter_Peak, Scatter_Last, Scatter_Ended, ncol = 2, nrow = 2)
}

dev.off()

#####Scatter plots of Dayofyear ~ Div_Temp###--Exact same as for year
#first flowers, no comments, see comments for the previous scatter plots (versus year, instead of versus temp)
model_First <- lm(Dayofyear ~ Div_Temp, data = First)
summary(model_First)
if (coef(model_First)[2] > 0) {
  eq_First <- substitute(italic(y) == a + b %.% italic(x), 
                         list(a = format(unname(coef(model_First)[1]), digits = 3), b = format(unname(coef(model_First)[2]), digits = 3)))
}
if (coef(model_First)[2] < 0) {
  eq_First <- substitute(italic(y) == a - b %.% italic(x), 
                         list(a = format(unname(coef(model_First)[1]), digits = 3), b = format(unname(-coef(model_First)[2]), digits = 3)))
}
Scatter_First <- ggplot(First, aes(x=Div_Temp, y = Dayofyear)) + geom_point() + xlim(TempMinDiv,TempMaxDiv) + ylim(DayMin,DayMax) +
  labs(title = "First Flowers", y = "Day of Year", x = "Yearly Mean Temperature, by Climate Division") + 
  geom_smooth(method = 'lm', se = FALSE) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_First))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_First

#graph for pre-flowering
model_Pre <- lm(Dayofyear ~ Div_Temp, data = Pre)
summary(model_Pre)
if (coef(model_Pre)[2] > 0) {
  eq_Pre <- substitute(italic(y) == a + b %.% italic(x), 
                       list(a = format(unname(coef(model_Pre)[1]), digits = 3), b = format(unname(coef(model_Pre)[2]), digits = 3)))
}
if (coef(model_Pre)[2] < 0) {
  eq_Pre <- substitute(italic(y) == a - b %.% italic(x), 
                       list(a = format(unname(coef(model_Pre)[1]), digits = 3), b = format(unname(-coef(model_Pre)[2]), digits = 3)))
}
Scatter_Pre <- ggplot(Pre, aes(x=Div_Temp, y = Dayofyear)) + geom_point() + xlim(TempMinDiv,TempMaxDiv) + ylim(DayMin,DayMax) +
  labs(title = "Pre-Flowering", y = "Day of Year", x = "Yearly Mean Temperature, by Climate Division") + 
  geom_smooth(method = 'lm', se = FALSE) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_Pre))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Pre

#graph for peak flowering
model_Peak <- lm(Dayofyear ~ Div_Temp, data = Peak)
summary(model_Peak)
if (coef(model_Peak)[2] > 0) {
  eq_Peak <- substitute(italic(y) == a + b %.% italic(x), 
                        list(a = format(unname(coef(model_Peak)[1]), digits = 3), b = format(unname(coef(model_Peak)[2]), digits = 3)))
}
if (coef(model_Peak)[2] < 0) {
  eq_Peak <- substitute(italic(y) == a - b %.% italic(x), 
                        list(a = format(unname(coef(model_Peak)[1]), digits = 3), b = format(unname(-coef(model_Peak)[2]), digits = 3)))
}
Scatter_Peak <- ggplot(Peak, aes(x=Div_Temp, y = Dayofyear)) + geom_point() + xlim(TempMinDiv,TempMaxDiv) + ylim(DayMin,DayMax) +
  labs(title = "Peak Flowering", y = "Day of Year", x = "Yearly Mean Temperature, by Climate Division") + 
  geom_smooth(method = 'lm', se = FALSE) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_Peak))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Peak

#graph for last flowers
model_Last <- lm(Dayofyear ~ Div_Temp, data = Last)
summary(model_Last)
if (coef(model_Last)[2] > 0) {
  eq_Last <- substitute(italic(y) == a + b %.% italic(x), 
                        list(a = format(unname(coef(model_Last)[1]), digits = 3), b = format(unname(coef(model_Last)[2]), digits = 3)))
}
if (coef(model_Last)[2] < 0) {
  eq_Last <- substitute(italic(y) == a - b %.% italic(x), 
                        list(a = format(unname(coef(model_Last)[1]), digits = 3), b = format(unname(-coef(model_Last)[2]), digits = 3)))
}
Scatter_Last <- ggplot(Last, aes(x=Div_Temp, y = Dayofyear)) + geom_point() + xlim(TempMinDiv,TempMaxDiv) + ylim(DayMin,DayMax) +
  labs(title = "Last Flowers", y = "Day of Year", x = "Yearly Mean Temperature, by Climate Division") + 
  geom_smooth(method = 'lm', se = FALSE) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_Last))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Last

#graph for flowering ended
model_Ended <- lm(Dayofyear ~ Div_Temp, data = Ended)
summary(model_Ended)
if (coef(model_Ended)[2] > 0) {
  eq_Ended <- substitute(italic(y) == a + b %.% italic(x), 
                         list(a = format(unname(coef(model_Ended)[1]), digits = 3), b = format(unname(coef(model_Ended)[2]), digits = 3)))
}
if (coef(model_Ended)[2] < 0) {
  eq_Ended <- substitute(italic(y) == a - b %.% italic(x), 
                         list(a = format(unname(coef(model_Ended)[1]), digits = 3), b = format(unname(-coef(model_Ended)[2]), digits = 3)))
}
Scatter_Ended <- ggplot(Ended, aes(x=Div_Temp, y = Dayofyear)) + geom_point() + xlim(TempMinDiv,TempMaxDiv) + ylim(DayMin,DayMax) +
  labs(title = "Flowering Ended", y = "Day of Year", x = "Yearly Mean Temperature, by Climate Division") + 
  geom_smooth(method = 'lm', se = FALSE) +
  stat_cor(label.y = DayMax, digits = 3, aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = DayMax-10) + aes(label = as.character(as.expression(eq_Ended))) +
  theme(panel.background = element_rect(fill = NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank())
Scatter_Ended

#Print the plot with the selected graphs

ScatterPlotFN <- paste(IOFolder, sp_name, "_", time_period, "_ScatterPlot_Div_Temp.pdf", sep = '')

pdf(file=ScatterPlotFN, width = 10, height = 8)
if (dim(Ended)[1] <= 2) {
  grid.arrange(Scatter_First, Scatter_Peak, Scatter_Last, ncol = 2, nrow = 2)
}
if (dim(Ended)[1] > 2) {
  grid.arrange(Scatter_First, Scatter_Peak, Scatter_Last, Scatter_Ended, ncol = 2, nrow = 2)
}

dev.off()

#####*****Making the box plots*****

#calculating median peak-flowering value for paper
median(Peak$Dayofyear)

#Y-axis labels for if there are no pre-flowering specimens
if (dim(Pre)[1] == 0) {
  MyLabs <- c("First","Peak","Last","Ended")
  MyLims <- c("1", "2", "3", "4")
}
#Y-axis labels for if there are pre-flowering specimens
if (dim(Pre)[1] != 0) {
  MyLabs <- c("Pre-flowering", "First","Peak","Last","Ended")
  MyLims <- c("0", "1", "2", "3", "4")
}

#drawing the box plot
spBoxPlot <- ggplot(sp_data, aes(y = PhenoPhase, x = Dayofyear, group=PhenoPhase)) + geom_boxplot() + 
  labs(x = "Day of Year", y = "Phenophase") + 
  scale_x_continuous(limits=c(0,366)) + 
  scale_y_discrete(limits=MyLims, labels=MyLabs) +
  theme(
    panel.background = element_rect(fill = NA, colour = "black"),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank()
  )
spBoxPlot

BoxPlotFN <- paste(IOFolder, sp_name, "_BoxPlot.pdf", sep = '')

pdf(file=BoxPlotFN, width = 4, height = 3)
spBoxPlot
dev.off()

#####*****Making the map*****
#####*maps::map() has a conflict with purrr:map(), which is automatically loaded when tidyverse is loaded
#####*So I need to specify maps::map here.

#getting the state boundary, and making it an sf object
states <- st_as_sf(maps::map("state", plot = FALSE, fill = TRUE))
OK <- subset(states, grepl("oklahoma", states$ID))
#getting the county boundaries for within Oklahoma as sf objects
counties <- st_as_sf(maps::map("county", plot = FALSE, fill = TRUE))
OK_counties <- subset(counties, grepl("oklahoma", counties$ID))

#making the base map, without any points on it
#this follows https://www.r-spatial.org/r/2018/10/25/ggplot2-sf.html
OK_base <- ggplot(data=OK) + geom_sf() + theme_bw() + 
  geom_sf(color = "black", fill = "gray") +
  geom_sf(data = OK_counties, fill = NA, color = "white") + #making county boundaries white
  geom_sf(data = OK, fill = NA, color = "black") + #writing over the state boundary over them in black
  labs(x = "Longitude", y = "Latitude") +
  scale_y_continuous(breaks = c(37, 36, 35, 34)) +
  annotation_north_arrow(pad_x = unit(0.5, "cm"), pad_y = unit(1.5, "cm"), height = unit(2, "cm")) + 
  annotation_scale(location = "bl", pad_x = unit(0.5, "cm"), pad_y = unit(0.5, "cm"))
OK_base

#Specifying colors, shapes, and labels
#first, if there are no pre-flowering specimens
if (dim(Pre)[1] == 0) {
  MyCols <- c("blue", "orange", "magenta", "green")
  MyShapes <- c(25, 22, 23, 24)
  MyLabs <- c("First","Peak","Last","Ended")
}
#then if there are pre-flowering specimens
if (dim(Pre)[1] != 0) {
  MyCols <- c("red", "blue", "orange", "magenta", "green")
  MyShapes <- c(21, 25, 22, 23, 24)
  MyLabs <- c("Pre-flowering", "First","Peak","Last","Ended")
}
#making the map with colored points
OK_points <- OK_base + 
  geom_point(data = sp_data, aes(x = long, y = lat, group = sname, fill = PhenoPhase, shape = PhenoPhase),
             size = 2) + 
  scale_fill_manual(values=MyCols, name = "Phenophase", labels = MyLabs) + 
  scale_shape_manual(values=MyShapes, name = "Phenophase", labels = MyLabs) +
  theme(legend.position = c(0.2,0.5))
OK_points

#print map

Map1FN <- paste(IOFolder, sp_name, "_Map_PhenoPhases.pdf", sep = '')
pdf(file=Map1FN, width = 8, height = 6)
OK_points
dev.off()


#Georeferencing check, assuming that there is a column with the county the plants are supposed to be in.
#This loop just makes one plot per county, showing the individuals that are supposed to be in that county.
#Each plot has its own page, so it is easy to check them.
Sp_counties <- unique(sp_data$County)
Map2FN <- paste(IOFolder, sp_name, "_Map_counties.pdf", sep = '')
pdf(file=Map2FN, width = 8, height = 6)
for (County_Name in Sp_counties){
  County_Map <- OK_base + geom_point(data = subset(sp_data, County == County_Name), aes(x = long, y = lat, group = sname), fill = "blue", size = 2) +
    labs(title = County_Name)
  print(County_Map)
}

dev.off()