# ==============================================================================
# Title: Intra-Crown Phenological Analysis Code
# Description: 
#   This script performs the following steps to reproduce the study's findings:
#   1. Loads high-resolution UAV-derived time-series data.
#   2. Detects phenological events (Flowering/Senescence) using derivative analysis.
#   3. Aggregates grid-level data to individual tree-level statistics.
#   4. Generates key figures (Event detection examples and Statistical validation).
#
# System Requirements: R version 4.0 or higher
# Dependencies: dplyr, ggplot2, readr, tidyr, gridExtra, ggtext
# ==============================================================================

# ------------------------------------------------------------------------------
# 0. Setup and Libraries
# ------------------------------------------------------------------------------
if (!require("pacman")) install.packages("pacman")
pacman::p_load(dplyr, ggplot2, readr, tidyr, gridExtra, ggtext)

# File Paths (Adjust these paths based on your data location)
INPUT_FILE    <- "Input_Data.csv"       # Raw time-series data
OUTPUT_FIGURE_DETECTION <- "Figure_Event_Detection.tiff"
OUTPUT_FIGURE_STATS     <- "Figure_Statistical_Validation.tiff"

# ------------------------------------------------------------------------------
# 1. Data Loading and Preprocessing
# ------------------------------------------------------------------------------
message(">>> [Step 1] Loading and cleaning data...")

if (!file.exists(INPUT_FILE)) stop("Error: Input file not found.")

raw_data <- read_csv(INPUT_FILE, show_col_types = FALSE) %>%
  filter(!is.na(DOY))

# Normalize Species Names and Assign Codes
# RY: Rhododendron yedoense f. poukhanense (Target for Spring)
# AT: Acer triflorum (Target for Autumn)
raw_data <- raw_data %>%
  mutate(
    Species_Code = case_when(
      grepl("Rhododendron|R\\. yedoense", Species) ~ "RY",
      grepl("Acer|A\\. triflorum", Species) ~ "AT",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(Species_Code))

# Check Data Availability
n_ry <- nrow(raw_data[raw_data$Species_Code == "RY" & raw_data$DOY >= 100 & raw_data$DOY <= 150, ])
n_at <- nrow(raw_data[raw_data$Species_Code == "AT" & raw_data$DOY >= 250 & raw_data$DOY <= 300, ])

message(sprintf(" - Loaded Data: %d rows", nrow(raw_data)))
message(sprintf(" - Valid Data Points: RY (Spring) = %d, AT (Autumn) = %d", n_ry, n_at))

if (n_ry == 0 && n_at == 0) stop("Error: No valid data found within the target DOY windows.")

# ------------------------------------------------------------------------------
# 2. Phenological Event Detection (Derivative Analysis)
# ------------------------------------------------------------------------------
message(">>> [Step 2] Detecting phenological events using derivatives...")

# Function: Robust Event Detection
# - Fits a smoothing spline to the time-series vegetation index.
# - Calculates the first derivative.
# - Identifies the local extremum (Max slope for Spring, Min slope for Autumn) within the search window.
get_event_robust <- function(doy, value, mode, window) {
  
  # Filter valid observations
  valid_idx <- !is.na(doy) & !is.na(value)
  doy <- doy[valid_idx]
  value <- value[valid_idx]
  
  # Constraints: Minimum data points
  if(length(doy) < 5) return(NA)
  if(sum(doy >= window[1] & doy <= window[2]) < 3) return(NA)
  
  tryCatch({
    # Smoothing Spline
    fit <- smooth.spline(doy, value, spar = 0.4)
    
    # Derivative Calculation
    pred_doy <- seq(min(doy), max(doy), by = 1)
    preds <- predict(fit, pred_doy, deriv = 1)
    
    df_deriv <- data.frame(doy = preds$x, slope = preds$y)
    
    # Search within the defined window
    df_search <- df_deriv %>% filter(doy >= window[1] & doy <= window[2])
    
    if(nrow(df_search) == 0) return(NA)
    
    # Detect Peak (Spring) or Valley (Autumn)
    target_doy <- if(mode == "increase") {
      df_search$doy[which.max(df_search$slope)]
    } else {
      df_search$doy[which.min(df_search$slope)]
    }
    return(target_doy)
    
  }, error = function(e) return(NA))
}

# 2.1. Process RY (Spring Flowering: Max RCC slope)
results_ry <- raw_data %>%
  filter(Species_Code == "RY") %>%
  group_by(Cell_ID, Tree_ID) %>%
  summarise(
    Species_Code = "RY",
    Event_Date = get_event_robust(DOY, RCC, mode="increase", window=c(100, 150)),
    Solar_Mean = mean(Cum_Duration, na.rm=TRUE),
    .groups = 'drop'
  )

# 2.2. Process AT (Autumn Senescence: Min GCC slope)
results_at <- raw_data %>%
  filter(Species_Code == "AT") %>%
  group_by(Cell_ID, Tree_ID) %>%
  summarise(
    Species_Code = "AT",
    Event_Date = get_event_robust(DOY, GCC, mode="decrease", window=c(250, 300)),
    Solar_Mean = mean(Cum_Duration, na.rm=TRUE),
    .groups = 'drop'
  )

# 2.3. Merge Results (Grid-Level Data)
grid_metrics <- bind_rows(results_ry, results_at) %>%
  filter(!is.na(Event_Date))

message(sprintf(" - Events Detected: %d grid cells", nrow(grid_metrics)))

# ------------------------------------------------------------------------------
# 3. Aggregation: Grid to Individual Tree Level
# ------------------------------------------------------------------------------
message(">>> [Step 3] Aggregating statistics to tree level...")

tree_stats <- grid_metrics %>%
  group_by(Species_Code, Tree_ID) %>%
  summarise(
    N_Cells = n(),
    SD_Solar = sd(Solar_Mean, na.rm = TRUE),    # Intra-crown Microclimatic Heterogeneity
    SD_Pheno = sd(Event_Date, na.rm = TRUE),    # Phenological Buffering (Duration)
    
    # Calculate Phenological Velocity (Sunlit vs Shaded)
    # Sunlit: Top 20% solar duration, Shaded: Bottom 20%
    Date_Sunlit = mean(Event_Date[Solar_Mean >= quantile(Solar_Mean, 0.8, na.rm=T)], na.rm=T),
    Date_Shaded = mean(Event_Date[Solar_Mean <= quantile(Solar_Mean, 0.2, na.rm=T)], na.rm=T),
    .groups = 'drop'
  ) %>%
  filter(N_Cells >= 10, !is.na(SD_Solar), !is.na(SD_Pheno)) %>%
  mutate(
    # Velocity Gap: Positive value means sunlit parts extend the growing season
    Velocity_Gap = case_when(
      Species_Code == "RY" ~ Date_Shaded - Date_Sunlit, # Spring: Sunlit blooms earlier
      Species_Code == "AT" ~ Date_Sunlit - Date_Shaded  # Autumn: Sunlit colors later
    ),
    Label = ifelse(Species_Code == "RY", "R. yedoense f. poukhanense (RY)", "Acer triflorum (AT)")
  )

# Set Factor Levels for Visualization
tree_stats$Label <- factor(tree_stats$Label, levels = c("R. yedoense f. poukhanense (RY)", "Acer triflorum (AT)"))

# ------------------------------------------------------------------------------
# 4. Visualization
# ------------------------------------------------------------------------------
message(">>> [Step 4] Generating Figures...")

# Colors: RY=Orange, AT=Pink
species_colors <- c("R. yedoense f. poukhanense (RY)" = "orange", 
                    "Acer triflorum (AT)" = "deeppink")

# --- Figure A: Phenological Buffering (Regression) ---
plot_buffering <- ggplot(tree_stats, aes(x = SD_Solar, y = SD_Pheno, fill = Label)) +
  geom_smooth(method = "lm", se = TRUE, aes(color = Label), alpha = 0.15, linewidth=0.8) +
  geom_point(shape = 21, size = 4.5, color = "black", stroke = 0.6, alpha = 0.85) +
  facet_wrap(~Label, scales = "free") +
  scale_fill_manual(values = species_colors) +
  scale_color_manual(values = species_colors) +
  labs(
    title = "(A) Phenological Buffering Capacity",
    subtitle = "Linear Regression between Microclimatic Heterogeneity and Phenological Duration",
    x = "Intra-crown Solar Heterogeneity (SD)",
    y = "Phenological Duration (SD of Days)"
  ) +
  theme_bw() +
  theme(legend.position = "none", strip.text = element_text(face = "italic", size = 11),
        plot.title = element_text(face="bold", size=12))

# --- Figure B: Phenological Velocity (Boxplot) ---
plot_velocity <- ggplot(tree_stats, aes(x = Label, y = Velocity_Gap, fill = Label)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray40", linewidth=1) +
  geom_boxplot(alpha = 0.4, outlier.shape = NA, width = 0.5) +
  geom_jitter(width = 0.2, size = 3.5, shape = 21, color = "black", stroke=0.6, alpha = 0.8) +
  scale_fill_manual(values = species_colors) +
  labs(
    title = "(B) Phenological Velocity",
    subtitle = "Gap > 0: Sunlit parts extend growing season",
    x = "",
    y = "Time Lag (Days) [Sunlit Effect]"
  ) +
  theme_bw() +
  theme(legend.position = "none", axis.text.x = element_text(face = "italic", size = 10),
        plot.title = element_text(face="bold", size=12))

# Combine Plots
final_stats_plot <- grid.arrange(plot_buffering, plot_velocity, ncol = 2)

# Save Figures (High Resolution for Publication)
ggsave(OUTPUT_FIGURE_STATS, final_stats_plot, device="tiff", dpi=300, 
       width=12, height=6, units="in", compression="lzw")

message(">>> All processes completed successfully.")
message(sprintf("Output saved: %s", OUTPUT_FIGURE_STATS))