library(tidyverse)
library(ggbeeswarm)
library(ggpubr)
library(plotly)
library(htmlwidgets)
library(growthcurver)
library(FactoMineR)
library(ggrepel)
#
theme_AGL <- theme_bw() + 
  theme(strip.background   = element_rect(colour = "black", fill = NA, size = 0.5),
        panel.background   = element_rect(colour = "black", fill = NA, size = 0.5),
        panel.border       = element_rect(colour = "black", size = 0.5),
        panel.grid         = element_line(color  = alpha("black", 0.1), size = 0.5),
        panel.grid.minor.x = element_blank(), 
        panel.grid.minor.y = element_blank())
#
myExpts1 <- c("Metaponto, Italy 2017", "Rosthern, Canada 2017", 
              "Sutherland, Canada 2017", "Sutherland, Canada 2018")
myExpts2 <- c("It17", "Ro17", "Su17", "Su18")
myColors_Expt <- c("darkred", "darkgreen", "darkblue", "steelblue")
myColors_Cluster <- c("red4", "darkorange3", "blue2", "deeppink3", 
                      "steelblue", "darkorchid4", "darkslategray", "chartreuse4")

myLDP <- read.csv("myLDP.csv")
#
it17_f <- read.csv("Data_Manual/LDP_Italy_2017_Phenotypes.csv") 
ro17_f <- read.csv("Data_Manual/LDP_Rosthern_2017_Phenotypes.csv")
su17_f <- read.csv("Data_Manual/LDP_Sutherland_2017_Phenotypes.csv")
su18_f <- read.csv("Data_Manual/LDP_Sutherland_2018_Phenotypes.csv")

# function to read in the drone data
readDrone <- function(folder, colNum, xf) {
  fnames <- list.files(folder)
  xx <- NULL
  for(i in fnames) {
    xi <- read.csv(paste0(folder, i))
    xx <- bind_rows(xx, xi)
  }
  xx <- xx %>% 
    rename(Plot=plot_name, Date=mission_date, Sensor=sensor_type) %>%
    mutate(Plot = as.integer(Plot)) %>%
    gather(Trait, Value, colNum:ncol(.)) %>%
    left_join(xf[,1:10], by = "Plot") %>%
    filter(!is.na(Entry), !is.na(Value)) %>%
    mutate(DAP = as.numeric(as.Date(Date) - as.Date(Planting.Date)),
           Value_f = Value) %>%
    select(Expt, Row, Column, Plot, Rep, Entry, Name, 
           Planting.Date, DAP, Date, Sensor, Trait, Value, Value_f)
  xx
}
# Load drone data
it17_d <- readDrone(folder = "Data_Drone/LDP_Italy_2017/",
                    colNum = 6, xf = it17_f)
ro17_d <- readDrone(folder = "Data_Drone/LDP_Rosthern_2017/",
                    colNum = 7, xf = ro17_f)
su17_d <- readDrone(folder = "Data_Drone/LDP_Sutherland_2017/",
                    colNum = 7, xf = su17_f)
su18_d <- readDrone(folder = "Data_Drone/LDP_Sutherland_2018/",
                    colNum = 7, xf = su18_f)

it17_f <- it17_f %>% gather(Trait, Value, 11:ncol(.))
ro17_f <- ro17_f %>% gather(Trait, Value, 11:ncol(.))
su17_f <- su17_f %>% gather(Trait, Value, 11:ncol(.))
su18_f <- su18_f %>% gather(Trait, Value, 11:ncol(.))

# Function to add data points
traitAdd <- function(xx, myTrait, myDAP, myValue) {
  for(i in unique(xx$Plot)) {
    xi <- xx %>% 
      filter(Trait == myTrait, Plot == i) %>%
      slice(1) %>%
      mutate(DAP = myDAP, Sensor = "NA", 
             Date = as.character(as.Date(Planting.Date) + myDAP), 
             Value = NA, Value_f = myValue)
    xx <- bind_rows(xi,xx)
  }
  xx
}
# Function to add a data point during the "stationary phase"
traitPredict <- function(xx, myDAP, myTrait) {
  for(i in unique(xx$Plot)) {
    xi <- xx %>% filter(Plot == i, Trait == myTrait)
    if(sum(!is.na(xi$Value_f)) > 0) {
      myVal <- xi %>% pull(Value_f) %>% max(na.rm = T)
      xi <- xi %>% slice(1) %>% 
        mutate(DAP = myDAP, Sensor = "NA", 
               Date = as.character(as.Date(Planting.Date) + myDAP), 
               Value = NA,
               Value_f = myVal)
      xx <- bind_rows(xx, xi)
    }
  }
  xx
}
# Function to find outliers after a certain DAP based on decreases from max
outliers_Maturity <- function(xx, myDAP, myThreshold = 1,
                              predict = T,
                              myTrait, myOtherTraits ) {
  yy <- NULL
  for(i in unique(xx$Plot)) {
    myMax <- max(xx %>% filter(Plot == i, Trait == myTrait) %>% pull(Value_f), na.rm = T)
    xi <- xx %>% filter(Plot == i)
    daps_r <- xi %>% 
      filter(DAP >= myDAP, Trait == myTrait, Value_f < myMax * myThreshold ) %>% 
      pull(DAP)
    daps_r <- c(daps_r, unique(xi$DAP)[unique(xi$DAP)>daps_r])
    xi <- xi %>%
      mutate(Value_f = ifelse(DAP %in% daps_r, NA, Value_f))
    if(predict == T) {
      xi <- xi %>%
        mutate(Value_f = ifelse(DAP %in% daps_r & Trait == myTrait, 
                                myMax*0.99, Value_f))
      for(k in myOtherTraits) {
        myMax <- max(xi %>% filter(Trait == k) %>% pull(Value_f), na.rm = T)
        xi <- xi %>%
          mutate(Value_f = ifelse(DAP %in% daps_r & Trait == k, 
                                  myMax*0.99, Value_f))
      }
    }
    yy <- bind_rows(yy, xi)
  }
  yy
}
# Function to identify and fix outliers
outliers_Fix <- function(xx, myTrait) {
  for(i in unique(xx$Plot)) {
    xi <- xx %>% filter(Plot == i, Trait == myTrait) %>%
      mutate(Change = NA, PercChange = NA)
    j<-2
    for(j in 2:nrow(xi)) {
      x1 <- max(xi$Value_f[1:(j-1)], na.rm = T)
      x2 <- xi$Value_f[j]
      xi$Change[j] <- x2 - x1
      xi$PercChange[j] <- 100 * xi$Change[j] / x1
    }
    for(j in 2:(nrow(xi)-1)) {
      xx$Value_f[xx$Plot == i & xx$Trait == myTrait & xx$DAP == xi$DAP[j]] <- 
        ifelse(xi$PercChange[j] > 30 & xi$PercChange[j+1] < -30, NA, xi$Value_f[j])
    }
  }
  xx
}

ggDroneCheck <- function(xd = it17_d, xf = it17_f,
                         myTitle = NULL,
                         myFilename = "DataCheck/RawDrone_Italy_2017.pdf",
                         xt = c("crop_volume_m3_excess.green.based",
                                "crop_area_m2_excess.green.based",
                                "mean_ndyi", 
                                "mean_crop_height_m_excess.green.based",
                                "median_crop_height_m_excess.green.based", 
                                "max_crop_height_m_excess.green.based") ) {
  myColors <- c("darkgreen", "darkorange", "darkred")
  xd <- xd %>% arrange(Entry)
  #
  pdf(myFilename, width = 16, height = 10)
  for(i in unique(xd$Entry)) {
    xdi <- xd %>% filter(Entry == i)
    xfi <- xf %>% filter(Entry == i)
    j <- unique(xdi$Name)
    #
    ggPlotDrone <- function(xti=xt[1]) {
      yMax <- max(xd%>%filter(Trait==xti)%>%pull(Value_f), na.rm = T)
      xMax <- max(c(xd$DAP, xf%>%filter(Trait=="DTM")%>%pull(Value)), na.rm = T)
      xdtf <- xfi %>% filter(Trait == "DTF") 
      xdtm <- xfi %>% filter(Trait == "DTM") 
      xd <- xdi %>% filter(Trait == xti)
      xd2 <- xd %>% filter(!is.na(Value_f))
      ggplot(xd, aes(x = DAP, y = Value)) +
        geom_line(alpha = 0.2, color = "black", size = 1) +
        geom_line(data = xd2, aes(y = Value_f, color = as.factor(Rep)), 
                  size = 1.5, alpha = 0.7) +
        geom_point(aes(y = Value_f), color = "red", size = 0.5) + 
        geom_point(aes(shape = Sensor)) + 
        geom_vline(data = xdtf, aes(xintercept = Value), 
                   alpha = 0.5, color = "darkgreen") +
        geom_vline(data = xdtm, aes(xintercept = Value), 
                   alpha = 0.5, color = "darkorange", lty = 2) +
        facet_grid(. ~ paste(Rep, "-", Plot), scales = "free_y") +
        scale_color_manual(values = myColors) +
        theme_AGL +
        theme(legend.position = "none") +
        coord_cartesian(ylim = c(0,yMax), xlim = c(0,xMax)) + 
        labs(y = xti)
    }
    #
    mp1 <- ggPlotDrone(xti = xt[1])
    mp2 <- ggPlotDrone(xti = xt[2])
    mp3 <- ggPlotDrone(xti = xt[3])
    #
    mp <- ggarrange(mp1, mp2, mp3, ncol = 1, nrow = 3) %>%
      annotate_figure(top = text_grob(paste(myTitle, "Entry", i, j)))
    print(mp)
  }
  dev.off()
}

ggDroneTrait <- function(xx, myTitle, myDAPs = NULL, myTraits, myTraitNames) {
  # Prep data
  if(is.null(myDAPs)) { myDAPs <- unique(xx$DAP)[order(unique(xx$DAP))] }
  xx <- xx %>% filter(Trait %in% myTraits, DAP %in% myDAPs) %>%
    mutate(Trait = plyr::mapvalues(Trait, myTraits, myTraitNames),
           Trait = factor(Trait, levels = myTraitNames),
           Group = factor(paste("Rep", Rep), 
                        levels = c("Rep 1", "Rep 2", "Rep 3", "Outliers")),
           DAP = paste("DAP", DAP),
           DAP = factor(DAP, levels = paste("DAP", myDAPs)))
  xx <- xx %>% filter(!is.na(Value_f))
  # Plot
  ggplot(xx, aes(x = Rep, color = Group, pch = Group)) +
    geom_quasirandom(aes(y = Value_f), alpha = 0.25) +
    facet_grid(Trait ~ DAP, scales = "free_y", switch = "y") +
    scale_color_manual(name = NULL, 
                       values = c("darkslategray","steelblue","darkslategray4","black")) +
    scale_shape_manual(name = NULL, values = c(20,20,20,10)) +
    scale_y_continuous(sec.axis = sec_axis(~ .)) +
    theme_AGL +
    theme(legend.position = "bottom",
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank()) +
    labs(title = myTitle, y = NULL, x = NULL)
}
# Italy 2017
xx <- it17_d %>%
  mutate(DAP = plyr::mapvalues(DAP, c(93, 101, 119, 133, 163), 
                               c(94, 100, 118, 132, 164)))
mp <- ggDroneTrait(xx = xx, myTitle = "Italy 2017", 
                   myDAPs = c(94,100,118,132,164),
                   myTraits = c("crop_volume_m3_excess.green.based",
                                "mean_crop_height_m_excess.green.based",
                                "crop_area_m2_excess.green.based"),
                   myTraitNames = c("Volume", "Height", "Area"))
ggsave("Additional/ggDroneTrait_It17.png", mp, width = 10, height = 4.5)
# Rosthern 2017
mp <- ggDroneTrait(xx = ro17_d, myTitle = "Rosthern 2017", 
                   myDAPs = c(24,32,40,54,61,67,74,80,89,95),
                   myTraits = c("crop_volume_m3_blue.ndvi.based",
                                "mean_crop_height_m_blue.ndvi.based",
                                "crop_area_m2_blue.ndvi.based"),
                   myTraitNames = c("Volume", "Height", "Area"))
ggsave("Additional/ggDroneTrait_Ro17.png", mp, width = 12, height = 3.5)
# Sutherland 2017
mp <- ggDroneTrait(xx = su17_d, myTitle = "Sutherland 2017", 
                   myDAPs = c(32,40,46,51,57,61,68,72,75,77,83),
                   myTraits = c("crop_volume_m3_blue.ndvi.based",
                                "mean_crop_height_m_blue.ndvi.based", 
                                "crop_area_m2_blue.ndvi.based"),
                   myTraitNames = c("Volume", "Height", "Area"))
ggsave("Additional/ggDroneTrait_Su17.png", mp, width = 12, height = 3.5)
# Sutherland 2018
mp <- ggDroneTrait(xx = su18_d, myTitle = "Sutherland 2018", 
                   myDAPs = c(42,55,58,61,65,68,72,76,79,86),
                   myTraits = c("crop_volume_m3_excess.green.based",
                                "mean_crop_height_m_excess.green.based",
                                "crop_area_m2_excess.green.based"),
                   myTraitNames = c("Volume", "Height", "Area"))
ggsave("Additional/ggDroneTrait_Su18.png", mp, width = 12, height = 3.5)

myGrowthCurve <- function(xd, xf, myFolder, myFilename, 
                          myTrait, myOrder = "Entries") {
  #
  xd <- xd %>% filter(Trait == myTrait)
  xm <- xf %>% select(Plot, Rep, Entry, Name, Expt, Planting.Date) %>% 
    filter(!duplicated(Plot))
  #
  xg <- NULL # Growth curve coefficients
  xp <- NULL # Predicted Values
  for(i in unique(xd$Plot)) {
    xgi <- NULL
    xdi <- xd %>% arrange(DAP) %>%
      filter(Plot == i, !is.na(Value_f))
    #
    if(nrow(xdi) > 2) {
      fit <- SummarizeGrowth(data_t = xdi$DAP, 
                             data_n = xdi$Value_f,
                             bg_correct = "none")
      #
      suppressMessages( xgi <- bind_rows(xgi, data.frame(t(c(i, unlist(fit$vals))))) )
      #
      xgi <- xgi %>% rename(y0=n0, y0_se=n0_se, y0_p=n0_p, x_mid=t_mid)
      for(k in 1:(ncol(xgi)-1)) { xgi[,k] <- as.numeric(xgi[,k]) }
      xgi <- xgi %>% mutate(y_mid = k/(1 + ((k - y0)/y0) * exp(-r * x_mid)),
                          b = -((r * x_mid) - y_mid) )
      #
      myDays <- min(xd$DAP):max(xd$DAP)
      xpi <- data.frame(Plot = i, DAP = myDays,
                       Predicted.Values = xgi$k/(1 + ((xgi$k - xgi$y0)/xgi$y0) * exp(-xgi$r * myDays)))
      xp <- bind_rows(xp, xpi)
      #
      gdays <- xpi %>% filter(DAP > (xgi$x_mid - 3), DAP < (xgi$x_mid + 3))
      if(nrow(gdays)>0) {
        myLM <- lm(Predicted.Values ~ DAP, gdays)
        xgi <- xgi %>% mutate(g.b = myLM[[1]][1],
                              g.r = myLM[[1]][2] )
      }
      #
      xpi_max <- xpi %>% filter(Predicted.Values >= xgi$k*0.95) %>% arrange(Predicted.Values) %>% slice(1)
      if(nrow(xpi_max)==0) { 
        xpi_max <- xpi %>% 
          filter(Predicted.Values >= xgi$k*0.94) %>% 
          arrange(Predicted.Values) %>% slice(1) 
      }
      if(nrow(xpi_max)==0) { 
        xpi_max <- xpi %>% 
          filter(Predicted.Values >= xgi$k*0.90) %>% 
          arrange(Predicted.Values) %>% slice(1) 
      }
      xpi_min <- xpi %>% filter(Predicted.Values >= xgi$k*0.05) %>% arrange(Predicted.Values) %>% slice(1)
      xgi$x.max <- ifelse(is.null(nrow(xpi_max)), max(xpi$DAP), xpi_max$DAP)
      xgi$x.min <- xpi_min$DAP
      xgi$x.d <- xgi$x.max - xgi$x_mid 
      xg <- bind_rows(xg, xgi)
    }
  }
  #
  xg[xg=="0"] <- NA
  xg <- xg %>% rename(Plot=1) %>%
    mutate(Plot = as.integer(Plot)) %>%
    left_join(xm, by = "Plot") %>%
    select(Plot, Rep, Entry, Name, Expt, Planting.Date,
           r, b, g.r, g.b, k, x.min, x.mid=x_mid, x.max, 
           x.d, x.gen=t_gen, auc.l=auc_l, auc.e=auc_e, y0, sigma, 
           r.se=r_se, r.p=r_p, k.se=k_se, k.p=k_p, y0.se=y0_se, y0.p=y0_p, df, note)
  xp <- xp %>% 
    left_join(xm, by = "Plot") %>%
    select(Plot, Rep, Entry, Name, Expt, Planting.Date, DAP, Predicted.Values)
  #
  yMax <- max(xd$Value_f, na.rm = T)
  xMax <- max(xd$DAP, na.rm = T)
  xpm <- NULL
  #
  if(myOrder == "Entries") { myOrders <- 1:324 }
  if(myOrder == "DTF") { myOrders <- xg %>% arrange(k) %>% pull(Entries) %>% unique() }
  #
  pdf(paste0(myFolder, "PDF_", myFilename, ".pdf"), width = 8, height = 4)
  for(i in myOrders) {
    xdi <- xd %>% filter(Entry == i, !is.na(Value_f))
    xfi <- xf %>% filter(Entry == i)
    xgi <- xg %>% filter(Entry == i)
    xpi <- xp %>% filter(Entry == i)
    #
    xpmi <- xpi %>% group_by(Entry, Name, DAP) %>%
      summarise(Predicted.Values = mean(Predicted.Values, na.rm = T))
    xpm <- bind_rows(xpm, xpmi)
    #
    xdtf <- xfi %>% filter(Trait == "DTF") 
    xdtm <- xfi %>% filter(Trait == "DTM") 
    xdtf2 <- xdtf %>% group_by(Entry) %>% summarise(Value = mean(Value, na.rm =T))
    xdtm2 <- xdtm %>% group_by(Entry) %>% summarise(Value = mean(Value, na.rm =T))
    myDays <- min(xdi$DAP):max(xdi$DAP)  
    #
    myTitle <- unique(paste(stringr::str_pad(xdi$Entry, 3, pad = "0"), "|", xdi$Name))
    #
    mp1 <- ggplot(xdi, aes(x = DAP)) +
      geom_vline(data = xdtf, aes(xintercept = Value), color = "darkgreen", alpha = 0.5) +
      geom_vline(data = xdtm, aes(xintercept = Value), color = "darkorange", alpha = 0.5) +
      geom_point(aes(y = Value_f), color = "red", alpha = 0.7) +
      geom_point(aes(y = Value), alpha = 0.5) +
      geom_line(data = xpi, aes(y = Predicted.Values), alpha = 0.8) +
      facet_grid(. ~ paste(Rep, Plot, sep = "-")) +
      coord_cartesian(ylim = c(0,yMax), xlim = c(0,xMax)) +
      theme_AGL +
      labs(title = myTitle, y = myTrait)
    mp2 <- ggplot(xpi, aes(x = DAP, y = Predicted.Values)) +
      geom_vline(data = xdtf2, aes(xintercept = Value), color = "darkgreen", alpha = 0.5) +
      geom_vline(data = xdtm2, aes(xintercept = Value), color = "darkorange", alpha = 0.5) +
      geom_line(alpha = 0.3, aes(group = Plot)) +
      geom_line(data = xpmi, color = "red") +
      facet_grid(. ~ paste("Entry",Entry)) +
      coord_cartesian(ylim = c(0,yMax), xlim = c(0,xMax)) +
      theme_AGL +
      labs(title = "")
    mp <- ggarrange(mp1, mp2, widths = c(1,0.5))
    print(mp)
  }
  dev.off()
  #
  list(xd, # input data 
       xg, # Growth curve coefficients
       xp, # Predicted Values
       xpm) #
}
ggGrowthCurve <- function(xpm, myTitle = NULL, 
                          myEntries = c(9, 35, 73, 94, 113),
                          myColors = c("darkgreen", "darkgoldenrod3", "darkred",
                                       "steelblue4", "darkslategray") ) {
  #
  
  if(sum(colnames(xpm)=="Expt")==0) { xpm <- xpm %>% mutate(Expt = myTitle) }
  xpe <- xpm %>% filter(Entry %in% myEntries) %>%
    mutate(Entry = factor(Entry, levels = myEntries))
  mp <- ggplot(xpm, aes(x = DAP, y = Predicted.Values, group = Entry, key1 = Name)) +
    geom_line(alpha = 0.05) +
    geom_line(data = xpe, size = 2, aes(color = paste0(Name," (",Entry,")"))) +
    facet_grid(. ~ Expt, scales = "free_x", space = "free_x") +
    scale_color_manual(name = NULL, values = myColors) +
    theme_AGL + 
    labs(x = "Days After Planting")
  mp
}

# Prep data
xE <- read.csv("https://raw.githubusercontent.com/derekmichaelwright/AGILE_LDP_Phenology/master/data/data_env.csv") %>%
  rename(`Days After Planting`=DaysAfterPlanting) %>%
  filter(Expt %in% myExpts1) %>%
  select(Expt, `Days After Planting`, DayLength, Temp_mean) %>%
  gather(Trait, Value, DayLength, Temp_mean) %>%
  mutate(Unit = ifelse(Trait == "DayLength", "Hours", "Degrees Celcius"),
         Trait = plyr::mapvalues(Trait, c("DayLength", "Temp_mean"), 
                                 c("DayLength", "Mean Temperature")))
yy <- bind_rows(it17_f %>% filter(Trait == "DTF"),
                ro17_f %>% filter(Trait == "DTF"),
                su17_f %>% filter(Trait == "DTF"),
                su18_f %>% filter(Trait == "DTF") ) %>%
  mutate(Expt = plyr::mapvalues(Expt, myExpts2, myExpts1)) %>%
  group_by(Expt, Name) %>% 
  summarise(Value = mean(Value, na.rm = T))
# Plot
mp1 <- ggplot(yy, aes(x = Value, fill = Expt)) +
  geom_density(alpha = 0.7) +
  scale_fill_manual(name = NULL, values = myColors_Expt) +
  scale_x_continuous(breaks = seq(30, 160, by = 10)) +
  theme_AGL +
  theme(legend.position = "none") +
  labs(title = "A) Days to Flower", 
       x = "Days After Planting", y = "Density")
mp2 <- ggplot(xE %>% filter(Trait == "DayLength"), 
              aes(x = `Days After Planting`, y = Value, color = Expt)) +
  geom_line(size = 1.25, alpha = 0.7) +
  scale_color_manual(name = NULL, values = myColors_Expt) +
  theme_AGL +
  theme(legend.position = "bottom") +
  labs(title = "B) Day Length", 
       y = "Hours", x = "Days After Planting")
mp3 <- ggplot(xE %>% filter(Trait == "Mean Temperature"), 
              aes(x = `Days After Planting`, y = Value, color = Expt)) +
  geom_line(size = 0.75, alpha = 0.7) +
  scale_color_manual(name = NULL, values = myColors_Expt) +
  theme_AGL +
  theme(legend.position = "bottom") +
  labs(title = "C) Mean Temperature", 
       y = "Degrees Celcius", x = "Days After Planting")
mp <- ggarrange(mp1, ggarrange(mp2, mp3, ncol = 2, nrow = 1, 
                               common.legend = T, legend = "bottom"), 
                ncol = 1, nrow = 2, heights = c(0.75,1) )
ggsave("Figure_01.png", mp, width = 8, height = 6, dpi = 600, bg = "white")

# Create plotting function
ggGrowthCoefs <- function(xd, xf, myPlot, myYmax = NULL, 
                          minDAP = NULL, maxDAP = NULL, 
                          myTitle = NULL, mySubtitle = F) {
  #
  if(is.null(minDAP)) { minDAP <- min(xd[[3]]$DAP) }
  if(is.null(maxDAP)) { maxDAP <- max(xd[[3]]$DAP) }
  if(is.null(myYmax)) { myYmax <- max(xd[[2]]$k) }
  #
  x1 <- xd[[2]] %>% filter(Plot == myPlot)
  x2 <- xd[[3]] %>% filter(Plot == myPlot)
  #
  if(mySubtitle == T) { 
    mySubtitle <- paste("Plot", x1$Plot, "| Entry", x1$Entry, "|",  x1$Name) 
  } else { mySubtitle <- NULL }
  #
  ggplot(x1) +
    #
    geom_hline(yintercept = 0, lty = 2, alpha = 0.5) +
    geom_line(data = x2, aes(x = DAP, y = Predicted.Values), size = 1.5) +
    #
    geom_hline(aes(yintercept = k, color = "k"), size = 1) +
    geom_abline(aes(slope = g.r, intercept = g.b, color = "g.r"), size = 1) +
    geom_point(aes(x = x.mid, y = k/2)) +
    #
    geom_point(aes(x = x.min), y = 0) +
    geom_text(aes(x = x.min), y = -0.01, label = "x.min") +
    #
    geom_segment(y = 0, aes(x = x.mid, xend = x.mid, yend = k/2), lty = 2) +
    geom_point(y = 0, aes(x = x.mid)) +
    geom_text(y = -0.01, aes(x = x.mid, label = "x.mid")) +
    #
    geom_segment(y = 0, aes(x = x.max, xend = x.max, yend = k), lty = 2) +
    geom_point(aes(x = x.max, y = k)) +
    geom_point(aes(x = x.max, y = 0)) +
    geom_text(aes(x = x.max, y = -0.01), label = "x.max") +
    #
    geom_segment(aes(x = x.mid, xend = x.max, y = k/10, yend = k/10), lty = 2) +
    geom_text(aes(x = x.max-((x.max-x.mid)/2), y = 1.5*k/10), label = "x.d") +
    #
    theme_AGL +
    scale_color_manual(name = NULL, breaks = c("k","r","g.r"),
                       values = c("purple","steelblue","darkgreen")) +
    scale_x_continuous(limits = c(minDAP, maxDAP), expand = c(0,0)) +
    ylim(c(-0.01,myYmax)) +
    labs(title = myTitle, subtitle = mySubtitle,
         y = "Predicted Value (y)", x = "Days After Planting (x)",
         caption = "y = k / (1 + ((k - y0) / y0) * exp(-r * x))")
}
#
mp <- ggGrowthCoefs(xd = ro17_gc.v, xf = ro17_f, myPlot = 6322, 
                     myYmax = 0.175, minDAP = 20, maxDAP = 65)
ggsave("Figure_02.png", mp, width = 6, height = 3, dpi = 600)

ggDroneTrait <- function(xx, myTitle, myDAPs = NULL, myTraits, 
                         myTraitNames = myTraits ) {
  # Prep data
  if(is.null(myDAPs)) { myDAPs <- unique(xx$DAP)[order(unique(xx$DAP))] }
  xx <- xx %>% filter(Trait %in% myTraits, DAP %in% myDAPs) %>%
    mutate(Trait = plyr::mapvalues(Trait, myTraits, myTraitNames),
           Trait = factor(Trait, levels = myTraitNames),
           Rep = factor(paste("Rep", Rep)),
           DAP = paste("DAP", DAP),
           DAP = factor(DAP, levels = paste("DAP", myDAPs)) ) %>%
    filter(!is.na(Value_f), !is.na(Value))
  # Plot
  ggplot(xx, aes(x = Rep, color = Rep, pch = Rep)) +
    geom_quasirandom(aes(y = Value_f), size = 0.3, alpha = 0.7) +
    facet_grid(Trait ~ DAP, scales = "free_y", switch = "y") +
    scale_color_manual(name = NULL, 
                       values = c("darkslategray","slategray4","darkslategray4")) +
    scale_shape_manual(name = NULL, values = c(16,15,18)) +
    scale_y_continuous(sec.axis = sec_axis(~ .)) +
    theme_AGL +
    theme(strip.placement = "outside",
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank()) +
    guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1))) +
    labs(title = myTitle, y = NULL, x = NULL, 
         caption = paste("Number of data points =", nrow(xx)))
}
# Italy 2017
xx <- it17_d %>%
  mutate(DAP = plyr::mapvalues(DAP, c(93, 101, 119, 133, 163), 
                               c(94, 100, 118, 132, 164)))
mp1 <- ggDroneTrait(xx = xx, myTitle = "A) Metaponto, Italy 2017", 
                   myDAPs = c(94,100,118,132,164),
                   myTraits = c("crop_volume_m3_excess.green.based",
                                "mean_crop_height_m_excess.green.based",
                                "crop_area_m2_excess.green.based"),
                   myTraitNames = c("Volume", "Height", "Area"))
# Rosthern 2017
mp2 <- ggDroneTrait(xx = ro17_d, myTitle = "B) Rosthern, Canada", 
                   myDAPs = c(24,32,40,54,61,67,74,80,89,95),
                   myTraits = c("crop_volume_m3_blue.ndvi.based",
                                "mean_crop_height_m_blue.ndvi.based",
                                "crop_area_m2_blue.ndvi.based"),
                   myTraitNames = c("Volume", "Height", "Area"))
# Sutherland 2017
mp3 <- ggDroneTrait(xx = su17_d, myTitle = "C) Sutherland, Canada", 
                   myDAPs = c(32,40,46,51,57,61,68,72,75,77,83),
                   myTraits = c("crop_volume_m3_blue.ndvi.based",
                                "mean_crop_height_m_blue.ndvi.based",
                                "crop_area_m2_blue.ndvi.based"),
                   myTraitNames = c("Volume", "Height", "Area"))
# Sutherland 2018
mp4 <- ggDroneTrait(xx = su18_d, myTitle = "D) Sutherland, Canada 2018", 
                   myDAPs = c(42,55,58,61,65,68,72,76,79,86),
                   myTraits = c("crop_volume_m3_excess.green.based",
                                "mean_crop_height_m_excess.green.based",
                                "crop_area_m2_excess.green.based"),
                   myTraitNames = c("Volume", "Height", "Area"))
#
mp <- ggarrange(mp1, mp2, mp3, mp4, ncol = 1, nrow = 4, 
                common.legend = T, legend = "bottom")
ggsave("Supplemental_Figure_01.png", mp, 
       width = 8, height = 11, dpi = 600, bg = "white")

ggsave("Additional/Supplemental_Figure_01_A.png", mp1, width = 10, height = 4.5)
ggsave("Additional/Supplemental_Figure_01_B.png", mp2, width = 12, height = 3.5)
ggsave("Additional/Supplemental_Figure_01_C.png", mp3, width = 12, height = 3.5)
ggsave("Additional/Supplemental_Figure_01_D.png", mp4, width = 12, height = 3.5)

# Create plotting function
ggDroneCorr <- function(xd, xf, myTitle, myTrait, xLab, yLab, myColor = "darkgreen") {
  # Prep data
  xd <- xd[[2]] %>% select(Plot, k)
  xf <- xf %>% filter(Trait == myTrait) %>% select(Plot, Entry, Name, Value)
  xx <- left_join(xf, xd, by = "Plot") %>%
    group_by(Entry, Name) %>%
    summarise(Value = mean(Value, na.rm = T),
              k = mean(k, na.rm = T)) %>%
    ungroup()
  myCorr <- cor(xx$Value, xx$k, use = "complete.obs")^2
  myLabel1 <- paste("italic(R)^2 ==", round(myCorr,2))
  # Plot
  ggplot(xx, aes(x = Value, y = k)) +
    stat_smooth(geom = "line", method = "lm", se = F, 
                size = 1, color = "black") +
    geom_point(alpha = 0.6, pch = 16, color = myColor) +
    geom_label(label = myLabel1, hjust = 0, parse = T,
              x = -Inf, y = Inf, vjust = 1 ) +
    facet_grid(. ~ as.character(myTitle)) +
    theme_AGL + 
    labs(x = xLab, y = yLab)
}
#
mp1 <- ggDroneCorr(xd = it17_gc.h, xf = it17_f, myColor = myColors_Expt[1],
                   myTitle = "Metaponto, Italy 2017", yLab = "Plot Height (UAV)",
                   myTrait = "Canopy.Height", xLab = "Canopy Height (Manual)")
mp2 <- ggDroneCorr(xd = ro17_gc.h, xf = ro17_f, myColor = myColors_Expt[2],
                   myTitle = "Rosthern, Canada 2017", yLab = "Plot Height (UAV)", 
                   myTrait = "Canopy.Height", xLab = "Canopy Height (Manual)")
mp3 <- ggDroneCorr(xd = su17_gc.h, xf = su17_f, myColor = myColors_Expt[3], 
                   myTitle = "Sutherland, Canada 2017", yLab = "Plot Height (UAV)", 
                   myTrait = "Canopy.Height", xLab = "Canopy Height (Manual)")
mp4 <- ggDroneCorr(xd = su18_gc.h, xf = su18_f, myColor = myColors_Expt[4],
                   myTitle = "Sutherland, Canada 2018", yLab = "Plot Height (UAV)", 
                   myTrait = "CanopyHeight", xLab = "Canopy Height (Manual)")
mp_a <- ggarrange(mp1, mp2, mp3, mp4, ncol = 4, nrow = 1) %>%
  annotate_figure(fig.lab = "A) Plot Height (UAV) x Canopy Height (Manual)",
                  fig.lab.pos = "top.left", fig.lab.size = 13, top = "") + 
  theme(plot.margin=unit(c(0.1,0,0,0), 'cm')) 
#
mp1 <- ggDroneCorr(xd = it17_gc.v, xf = it17_f, myColor = myColors_Expt[1],
                   myTitle = "Metaponto, Italy 2017", yLab = "Plot Volume (UAV)", 
                   myTrait = "Plot.Biomass", xLab = "Total Biomass (Manual)")
mp2 <- ggDroneCorr(xd = ro17_gc.v, xf = ro17_f, myColor = myColors_Expt[2],
                   myTitle = "Rosthern, Canada 2017", yLab = "Plot Volume (UAV)", 
                   myTrait = "Plot.Biomass", xLab = "Total Biomass (Manual)")
mp3 <- ggDroneCorr(xd = su17_gc.v, xf = su17_f, myColor = myColors_Expt[3],
                   myTitle = "Sutherland, Canada 2017", yLab = "Plot Volume (UAV)", 
                   myTrait = "Plot.Biomass", xLab = "Total Biomass (Manual)")
mp4 <- ggDroneCorr(xd = su18_gc.v, xf = su18_f, myColor = myColors_Expt[4],
                   myTitle = "Sutherland, Canada 2018", yLab = "Plot Volume (UAV)", 
                   myTrait = "Plot.Biomass", xLab = "Total Biomass (Manual)")
mp_b <- ggarrange(mp1, mp2, mp3, mp4, ncol = 4, nrow = 1) %>%
  annotate_figure(fig.lab = "B) Plot Volume (UAV) x Biomass (Manual)",
                  fig.lab.pos = "top.left", fig.lab.size = 13, top = "") + 
  theme(plot.margin=unit(c(0.1,0,0,0), 'cm'))
mp <- ggarrange(mp_a, mp_b, nrow = 2)
ggsave("Supplemental_Figure_02.png", mp, width = 10, height = 6, dpi = 600, bg = "white")

# Create plotting function
ggGrowthCurve <- function(xpm, xF, myTitle = NULL, 
                          xlab = "Days After Planting", ylab = "Predicted Values", 
                          myEntries = c(19, 94, 107, 27), #252
                          myColors = c("chartreuse4", "red4", "blue2", "steelblue"), 
                          linesize = 1.75 ) {
  #
  xF <- xF %>% filter(Trait == "DTF", Entry %in% myEntries) %>%
    group_by(Name, Entry) %>% summarise(Value = mean(Value, na.rm = T))
  if(sum(colnames(xpm)=="Expt")==0) { xpm <- xpm %>% mutate(Expt = myTitle) }
  xpe <- xpm %>% filter(Entry %in% myEntries) %>%
    mutate(Entry = factor(Entry, levels = myEntries)) %>%
    arrange(Entry) %>%
    mutate(Name = factor(Name))
  mp <- ggplot(xpm, aes(x = DAP, y = Predicted.Values, group = Entry, key1 = Name)) +
    geom_line(alpha = 0.1) +
    geom_line(data = xpe, size = linesize, alpha = 0.8,
              aes(color = Name)) + 
    geom_point(data = xF, y = 0, size = 1.5, alpha = 0.8, pch = 16,
               aes(x = Value, color = Name)) +
    facet_grid(. ~ Expt, scales = "free_x", space = "free_x") +
    scale_color_manual(name = NULL, values = myColors) +
    theme_AGL + 
    labs(x = xlab, y = ylab)
  mp
}

# Italy 2017
mp1 <- ggGrowthCurve(xpm = it17_gc.v[[4]], xF = it17_f, xlab = "",
                     ylab = "Plot Volume", myTitle = "Metaponto, Italy 2017")
mp2 <- ggGrowthCurve(xpm = it17_gc.h[[4]], xF = it17_f, xlab = "",
                     ylab = "Plot Height", myTitle = "Metaponto, Italy 2017")
mp3 <- ggGrowthCurve(xpm = it17_gc.a[[4]], xF = it17_f, xlab = "",
                     ylab = "Plot Area",  myTitle = "Metaponto, Italy 2017")
# Sutherland 2017
mp4 <- ggGrowthCurve(xpm = su17_gc.v[[4]], xF = su17_f, xlab = "",
                     ylab = "", myTitle = "Sutherland, Canada 2017")
mp5 <- ggGrowthCurve(xpm = su17_gc.h[[4]], xF = su17_f, xlab = "",
                     ylab = "", myTitle = "Sutherland, Canada 2017")
mp6 <- ggGrowthCurve(xpm = su17_gc.a[[4]], xF = su17_f, xlab = "Days After Planting",
                     ylab = "", myTitle = "Sutherland, Canada 2017")
# Rosthern 2017
mp7 <- ggGrowthCurve(xpm = ro17_gc.v[[4]], xF = ro17_f, xlab = "",
                     ylab = "", myTitle = "Rosthern, Canada 2017")
mp8 <- ggGrowthCurve(xpm = ro17_gc.h[[4]], xF = ro17_f, xlab = "",
                     ylab = "", myTitle = "Rosthern, Canada 2017")
mp9 <- ggGrowthCurve(xpm = ro17_gc.a[[4]], xF = ro17_f, xlab = "",
                     ylab = "", myTitle = "Rosthern, Canada 2017")
# 
mp <- ggarrange(mp1, mp4, mp7, mp2, mp5, mp8, mp3, mp6, mp9, 
                ncol = 3, nrow = 3, common.legend = T)
ggsave("Figure_03.png", mp, width = 8, height = 6, dpi = 600, bg = "white")

# Prep data
myColors_Origins <- c("darkgreen","darkorange", "darkred", "black")
myOrigins <- c("India", "Iran", "Canada", "Other")
myExpts3 <- c("Italy 2017", "Rosthern 2017", "Sutherland 2017", "Sutherland 2018")
#
xH <- rbind(it17_gc.h[[2]], ro17_gc.h[[2]], su17_gc.h[[2]]) %>%
  select(Plot, Entry, Name, Expt, Height=k)
xA <- rbind(it17_gc.a[[2]], ro17_gc.a[[2]], su17_gc.a[[2]]) %>%
  select(Plot, Entry, Name, Expt, Area=k)
xV <- rbind(it17_gc.v[[2]], ro17_gc.v[[2]], su17_gc.v[[2]]) %>%
  select(Plot, Entry, Name, Expt, Volume=k)
xS <- rbind(it17_f %>% filter(Trait == "Seed.Mass"),
            ro17_f %>% filter(Trait == "Seed.Mass"),
            su17_f %>% filter(Trait == "Seed.Mass"),
            su18_f %>% filter(Trait == "Seed.Mass")) %>%
  select(Plot, Entry, Name, Expt, Seed.Mass=Value)
xD <- rbind(it17_f %>% filter(Trait == "DTF"),
            ro17_f %>% filter(Trait == "DTF"),
            su17_f %>% filter(Trait == "DTF"),
            su18_f %>% filter(Trait == "DTF")) %>%
  select(Plot, Entry, Name, Expt, DTF=Value) 
xP <- read.csv("myPCA.csv") %>%
  mutate(Cluster = factor(Cluster)) %>%
  select(Entry, Name, Cluster, Origin)
#
xx <- left_join(xH, xA, by = c("Plot","Entry","Name","Expt")) %>%
  left_join(xV, by = c("Plot","Entry","Name","Expt")) %>%
  left_join(xS, by = c("Plot","Entry","Name","Expt")) %>%
  left_join(xD, by = c("Plot","Entry","Name","Expt")) %>%
  group_by(Entry, Name, Expt) %>%
  summarise(Height = mean(Height, na.rm = T),
            Area = mean(Area, na.rm = T),
            Volume = mean(Volume, na.rm = T),
            Seed.Mass = mean(Seed.Mass, na.rm = T),
            DTF = mean(DTF, na.rm = T)) %>%
  ungroup() %>%
  left_join(xP, by = c("Entry", "Name")) %>%
  mutate(Group = ifelse(Origin %in% myOrigins[1:3], Origin, myOrigins[4]),
         Expt = plyr::mapvalues(Expt, myExpts3, myExpts1)) %>%
  arrange(desc(Group))
# Plot
mp1 <- ggplot(xx, aes(x = Area, y = Height, color = Group, alpha = Group)) +
  geom_point(aes(key1 = Origin, key2 = Entry, key3 = Name), pch = 16) +
  facet_grid(. ~ Expt) +
  scale_color_manual(name = "Accession Origin", values = myColors_Origins) +
  scale_alpha_manual(name = "Accession Origin", values = c(0.8,0.8,0.8,0.2)) +
  theme_AGL +
  labs(title = "A) Height x Area", y = "Plot Height", x = "Plot Area")
#
mp2 <- ggplot(xx, aes(x = Volume, y = Seed.Mass, color = Group, alpha = Group)) +
  geom_point(aes(key1 = Origin, key2 = Entry, key3 = Name), pch = 16) +
  facet_grid(. ~ Expt) +
  scale_color_manual(name = "Accession Origin", values = myColors_Origins) +
  scale_alpha_manual(name = "Accession Origin", values = c(0.8,0.8,0.8,0.2)) +
  theme_AGL +
  labs(title = "B) Yield x Volume", x = "Plot Volume", y = "Plot Yield")
#
mp3 <- ggplot(xx, aes(x = DTF, y = Volume, key1 = Origin, color = Group, alpha = Group)) +
  geom_point(aes(key1 = Origin, key2 = Entry, key3 = Name), pch = 16) +
  facet_grid(. ~ Expt, scales = "free") +
  scale_color_manual(name = "Accession Origin", values = myColors_Origins) +
  scale_alpha_manual(name = "Accession Origin", values = c(0.8,0.8,0.8,0.2)) +
  theme_AGL +
  labs(title = "C) Volume x DTF", x = "Days From Sowing To Flower (DTF)", y = "Plot Volume")
#
mp <- ggarrange(mp1, mp2, mp3, nrow = 3, ncol = 1, common.legend = T, legend = "bottom")
ggsave("Figure_04.png", mp, width = 6, height = 7, dpi = 600, bg = "white")
#
ggsave("Additional/Figure_04_A.png", mp1, width = 8, height = 2.5, dpi = 600)
ggsave("Additional/Figure_04_B.png", mp2, width = 8, height = 2.5, dpi = 600)
ggsave("Additional/Figure_04_C.png", mp3, width = 8, height = 2.5, dpi = 600)
#
saveWidget(ggplotly(mp1), file="Additional/Figure_04_A.html")
saveWidget(ggplotly(mp2), file="Additional/Figure_04_B.html")
saveWidget(ggplotly(mp3), file="Additional/Figure_04_C.html")