# --- Step 1: Load necessary libraries ---
library(readr)   # For reading CSV files
library(dplyr)   # For data manipulation
library(zoo)     # For rolling sums

# --- Step 2: Import your dataset ---
# Replace "climatedataset.csv" with your actual file path if needed
rain_data <- read_csv("climatedataset.csv")

# --- Step 3: Convert 'Date' column to Date class ---
# Assuming Date is in YYYYMMDD format like 19830101
rain_data$Date <- as.Date(as.character(rain_data$Date), format = "%Y%m%d")

# --- Step 4: Check data loaded correctly ---
head(rain_data)

# --- Step 5: Filter data for the year of interest (e.g., 1983) ---
year_of_interest <- 1983
rain_year <- rain_data %>% filter(format(Date, "%Y") == as.character(year_of_interest))

# --- Step 6: Define the precise onset detection function ---
detect_onset_precise <- function(rainfall, 
                                 onset_cum_threshold = 25, 
                                 onset_window = 5, 
                                 min_rainy_days = 3,  # including first day
                                 min_first_day_rain = 1,
                                 dry_threshold = 1,
                                 max_dry_spell_length = 6,  # less than 7 days allowed
                                 dry_spell_check_days = 30) {
  n <- length(rainfall)
  rolling_sum <- zoo::rollapply(rainfall, width = onset_window, FUN = sum, align = "left", fill = NA)
  
  for (i in 1:(n - onset_window - dry_spell_check_days + 1)) {
    # Check cumulative rainfall in 5-day window
    if (!is.na(rolling_sum[i]) && rolling_sum[i] >= onset_cum_threshold) {
      window_rain <- rainfall[i:(i + onset_window - 1)]
      
      # Check first day rainfall >= 1 mm
      if (window_rain[1] < min_first_day_rain) {
        next  # skip this window
      }
      
      # Check at least two other days (besides first) with rainfall >= 1 mm
      rainy_days_other <- sum(window_rain[-1] >= min_first_day_rain)
      if (rainy_days_other < (min_rainy_days - 1)) {
        next  # skip this window
      }
      
      # Check dry spell in the next 30 days after the 5-day window
      dry_period <- rainfall[(i + onset_window):(i + onset_window + dry_spell_check_days - 1)]
      
      # Find lengths of consecutive dry days (< 1 mm)
      rle_dry <- rle(dry_period < dry_threshold)
      
      # Check if any dry spell length >= 7 days
      if (any(rle_dry$values == TRUE & rle_dry$lengths >= 7)) {
        next  # reject this window due to long dry spell
      }
      
      # If all criteria met, return index of first day of onset
      return(i)
    }
  }
  return(NA)  # no onset detected
}

# --- Step 7: Define function to detect onset within a date range ---
detect_onset_period <- function(data, start_date, end_date) {
  period_data <- data %>% filter(Date >= start_date & Date <= end_date)
  onset_idx <- detect_onset_precise(period_data$Rainfall)
  if (!is.na(onset_idx)) {
    return(period_data$Date[onset_idx])
  } else {
    return(NA)
  }
}

# --- Step 8: Define rainy season periods for bimodal climate ---
long_rains_start <- as.Date(paste0(year_of_interest, "-03-01"))
long_rains_end <- as.Date(paste0(year_of_interest, "-06-30"))

short_rains_start <- as.Date(paste0(year_of_interest, "-09-01"))
short_rains_end <- as.Date(paste0(year_of_interest, "-12-31"))

# --- Step 9: Detect long rains onset ---
long_rains_onset <- detect_onset_period(rain_year, long_rains_start, long_rains_end)

# --- Step 10: Detect short rains onset ---
short_rains_onset <- detect_onset_period(rain_year, short_rains_start, short_rains_end)

# --- Step 11: Print results ---
cat("Year:", year_of_interest, "\n")
cat("Long rains onset:", ifelse(is.na(long_rains_onset), "No onset detected", as.character(long_rains_onset)), "\n")
cat("Short rains onset:", ifelse(is.na(short_rains_onset), "No onset detected", as.character(short_rains_onset)), "\n")

## Inspect rainfall data for March–June 1983

library(ggplot2)

ggplot(rain_year %>% filter(Date >= long_rains_start & Date <= long_rains_end), aes(x = Date, y = Rainfall)) +
  geom_col(fill = "blue") +
  labs(title = "Rainfall March-June 1983", y = "Daily Rainfall (mm)")

#Inspect rainfall data for September-December 1983

library(ggplot2)

ggplot(rain_year %>% filter(Date >= short_rains_start & Date <= short_rains_end), aes(x = Date, y = Rainfall)) +
  geom_col(fill = "blue") +
  labs(title = "Rainfall September-December 1983", y = "Daily Rainfall (mm)")

#Inspect rainfall data for the whole year 1983

library(ggplot2)
library(dplyr)

# Assuming rain_year is your filtered data for the year 1983
# If not filtered yet, filter as:
# rain_year <- rain_data %>% filter(format(Date, "%Y") == "1983")

# Plot daily rainfall for the whole year
ggplot(rain_year, aes(x = Date, y = Rainfall)) +
  geom_col(fill = "blue") +
  labs(title = paste("Daily Rainfall for Year", unique(format(rain_year$Date, "%Y"))),
       x = "Date",
       y = "Rainfall (mm)") +
  theme_minimal() +
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# --- Load libraries ---
library(readr)
library(dplyr)
library(zoo)

# --- Import your dataset ---
rain_data <- read_csv("climatedataset.csv")
rain_data$Date <- as.Date(as.character(rain_data$Date), format = "%Y%m%d")

# --- Define onset detection function as per your precise criteria ---
detect_onset_precise <- function(rainfall, 
                                 onset_cum_threshold = 25, 
                                 onset_window = 5, 
                                 min_rainy_days = 3,  # including first day
                                 min_first_day_rain = 1,
                                 dry_threshold = 1,
                                 max_dry_spell_length = 6,  # less than 7 days allowed
                                 dry_spell_check_days = 30) {
  n <- length(rainfall)
  rolling_sum <- zoo::rollapply(rainfall, width = onset_window, FUN = sum, align = "left", fill = NA)
  
  for (i in 1:(n - onset_window - dry_spell_check_days + 1)) {
    if (!is.na(rolling_sum[i]) && rolling_sum[i] >= onset_cum_threshold) {
      window_rain <- rainfall[i:(i + onset_window - 1)]
      if (window_rain[1] < min_first_day_rain) next
      rainy_days_other <- sum(window_rain[-1] >= min_first_day_rain)
      if (rainy_days_other < (min_rainy_days - 1)) next
      
      dry_period <- rainfall[(i + onset_window):(i + onset_window + dry_spell_check_days - 1)]
      rle_dry <- rle(dry_period < dry_threshold)
      if (any(rle_dry$values == TRUE & rle_dry$lengths >= 7)) next
      
      return(i)
    }
  }
  return(NA)
}

# --- Function to detect onset within a period for one year ---
detect_onset_period <- function(data, start_date, end_date) {
  period_data <- data %>% filter(Date >= start_date & Date <= end_date)
  onset_idx <- detect_onset_precise(period_data$Rainfall)
  if (!is.na(onset_idx)) {
    return(period_data$Date[onset_idx])
  } else {
    return(NA)
  }
}

# --- Years to analyze ---
years <- 1983:2021

# --- Initialize results dataframe ---
results <- data.frame(
  Year = integer(),
  LongRainsOnset = as.Date(character()),
  ShortRainsOnset = as.Date(character()),
  stringsAsFactors = FALSE
)

# --- Loop over years ---
for (yr in years) {
  rain_year <- rain_data %>% filter(format(Date, "%Y") == as.character(yr))
  
  long_rains_start <- as.Date(paste0(yr, "-03-01"))
  long_rains_end <- as.Date(paste0(yr, "-06-30"))
  short_rains_start <- as.Date(paste0(yr, "-09-01"))
  short_rains_end <- as.Date(paste0(yr, "-12-31"))
  
  long_onset <- detect_onset_period(rain_year, long_rains_start, long_rains_end)
  short_onset <- detect_onset_period(rain_year, short_rains_start, short_rains_end)
  
  results <- rbind(results, data.frame(
    Year = yr,
    LongRainsOnset = long_onset,
    ShortRainsOnset = short_onset
  ))
}

# --- View results ---
print(results)

# Optional: save results to CSV
# write.csv(results, "rainy_season_onsets_1983_2021.csv", row.names = FALSE)

# --- Load necessary libraries ---
library(readr)
library(dplyr)
library(zoo)

# --- Import your dataset ---
# Replace "climatedataset.csv" with the correct path to your file
rain_data <- read_csv("climatedataset.csv")

# --- Convert 'Date' column to Date class ---
rain_data$Date <- as.Date(as.character(rain_data$Date), format = "%Y%m%d")

# --- Add Month column for Ra assignment ---
rain_data <- rain_data %>%
  mutate(Month = as.integer(format(Date, "%m")))

# --- Create Ra lookup table for Kigali (approximate FAO56 values at 0° latitude) ---
ra_lookup <- data.frame(
  Month = 1:12,
  Ra = c(36.2, 37.5, 37.9, 36.8, 34.8, 33.4, 33.9, 35.7, 37.2, 37.4, 36.3, 35.6)  # MJ m-2 day-1
)

# --- Join Ra values by Month ---
rain_data <- rain_data %>%
  left_join(ra_lookup, by = "Month")

# --- Calculate mean temperature and PET ---
rain_data <- rain_data %>%
  mutate(
    Tmean = (Tmax + Tmin) / 2,
    temp_diff = pmax(Tmax - Tmin - 0.0123 * Rainfall, 0),
    PET = 0.00102 * Ra * (temp_diff)^0.76 * (Tmean + 17)
  )

# --- Calculate rolling 10-day rainfall sum (including current day) ---
rain_data <- rain_data %>%
  arrange(Date) %>%
  mutate(
    Rainfall_10d = rollapply(Rainfall, width = 10, FUN = sum, align = "right", fill = NA, na.rm = TRUE)
  )

# --- Mark cessation candidate days where 10-day rainfall < 0.5 * PET ---
rain_data <- rain_data %>%
  mutate(
    CessationCandidate = ifelse(Rainfall_10d < 0.5 * PET, TRUE, FALSE)
  )

# --- Define cessation detection function for a given period ---
detect_cessation_period <- function(data, start_date, end_date) {
  period_data <- data %>% filter(Date >= start_date & Date <= end_date)
  cessation_days <- period_data %>% filter(CessationCandidate == TRUE)
  
  if (nrow(cessation_days) == 0) {
    return(NA)  # No cessation detected in this period
  }
  
  # Return earliest cessation candidate date (customize logic if needed)
  return(min(cessation_days$Date))
}

# --- Apply cessation detection for year 1983 ---

year_of_interest <- 1983

# Define bimodal rainy season periods for 1983
long_rains_start <- as.Date(paste0(year_of_interest, "-03-01"))
long_rains_end   <- as.Date(paste0(year_of_interest, "-06-30"))

short_rains_start <- as.Date(paste0(year_of_interest, "-09-01"))
short_rains_end   <- as.Date(paste0(year_of_interest, "-12-31"))

# Detect cessation dates
cessation_long <- detect_cessation_period(rain_data, long_rains_start, long_rains_end)
cessation_short <- detect_cessation_period(rain_data, short_rains_start, short_rains_end)

# Print results
cat("Year:", year_of_interest, "\n")
cat("Long rains cessation:", ifelse(is.na(cessation_long), "No cessation detected", as.character(cessation_long)), "\n")
cat("Short rains cessation:", ifelse(is.na(cessation_short), "No cessation detected", as.character(cessation_short)), "\n")

##for visualisation

library(ggplot2)

# Filter data for long rains period in 1983
long_rains_data <- rain_data %>%
  filter(Date >= long_rains_start & Date <= long_rains_end)

ggplot(long_rains_data, aes(x = Date)) +
  geom_col(aes(y = Rainfall), fill = "blue", alpha = 0.6) +
  geom_line(aes(y = PET), color = "red", size = 1) +
  geom_vline(xintercept = as.numeric(cessation_long), linetype = "dashed", color = "darkgreen") +
  labs(title = "Long Rains 1983: Rainfall, PET and Cessation",
       y = "Rainfall (mm) / PET (MJ m-2 day-1)") +
  theme_minimal()

library(ggplot2)

# Filter data for short rains period in 1983
short_rains_data <- rain_data %>%
  filter(Date >= short_rains_start & Date <= short_rains_end)

ggplot(short_rains_data, aes(x = Date)) +
  geom_col(aes(y = Rainfall), fill = "blue", alpha = 0.6) +
  geom_line(aes(y = PET), color = "red", size = 1) +
  geom_vline(xintercept = as.numeric(cessation_short), linetype = "dashed", color = "darkgreen") +
  labs(title = "Short Rains 1983: Rainfall, PET and Cessation",
       y = "Rainfall (mm) / PET (MJ m-2 day-1)",
       x = "Date") +
  theme_minimal()

library(ggplot2)
library(dplyr)

# Filter data for the year of interest, e.g., 1983
year_of_interest <- 1983
rain_year <- rain_data %>% filter(format(Date, "%Y") == as.character(year_of_interest))

# Plot daily Rainfall as bars and PET as a red line over the whole year
ggplot(rain_year, aes(x = Date)) +
  geom_col(aes(y = Rainfall), fill = "blue", alpha = 0.6) +
  geom_line(aes(y = PET), color = "red", size = 1) +
  labs(title = paste("Daily Rainfall and PET for Year", year_of_interest),
       y = "Rainfall (mm) / PET (MJ m⁻² day⁻¹)",
       x = "Date") +
  theme_minimal() +
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# --- Load necessary libraries ---
library(readr)
library(dplyr)
library(zoo)
library(ggplot2)
library(tidyr)

# --- Step 1: Import your dataset ---
rain_data <- read_csv("climatedataset.csv")
rain_data$Date <- as.Date(as.character(rain_data$Date), format = "%Y%m%d")

# --- Step 2: Prepare data ---
rain_data <- rain_data %>%
  mutate(
    Year = as.integer(format(Date, "%Y")),
    Month = as.integer(format(Date, "%m"))
  )

# --- Step 3: Assign extraterrestrial radiation (Ra) for Kigali ---
ra_lookup <- data.frame(
  Month = 1:12,
  Ra = c(36.2, 37.5, 37.9, 36.8, 34.8, 33.4, 33.9, 35.7, 37.2, 37.4, 36.3, 35.6)  # MJ m-2 day-1
)

rain_data <- rain_data %>%
  left_join(ra_lookup, by = "Month")

# --- Step 4: Calculate PET ---
rain_data <- rain_data %>%
  mutate(
    Tmean = (Tmax + Tmin) / 2,
    temp_diff = pmax(Tmax - Tmin - 0.0123 * Rainfall, 0),
    PET = 0.00102 * Ra * (temp_diff)^0.76 * (Tmean + 17)
  )

# --- Step 5: Calculate rolling 10-day rainfall sum ---
rain_data <- rain_data %>%
  arrange(Date) %>%
  group_by(Year) %>%
  mutate(
    Rainfall_10d = rollapply(Rainfall, width = 10, FUN = sum, align = "right", fill = NA, na.rm = TRUE)
  ) %>%
  ungroup()

# --- Step 6: Mark cessation candidate days ---
rain_data <- rain_data %>%
  mutate(
    CessationCandidate = ifelse(Rainfall_10d < 0.5 * PET, TRUE, FALSE)
  )

# --- Step 7: Define cessation detection function ---
detect_cessation_period <- function(data, start_date, end_date) {
  period_data <- data %>% filter(Date >= start_date & Date <= end_date)
  cessation_days <- period_data %>% filter(CessationCandidate == TRUE)
  if (nrow(cessation_days) == 0) return(NA)
  return(min(cessation_days$Date))
}

# --- Step 8: Loop over years to detect cessation for bimodal seasons ---
years <- unique(rain_data$Year)
results <- data.frame(
  Year = integer(),
  LongRainsCessation = as.Date(character()),
  ShortRainsCessation = as.Date(character()),
  stringsAsFactors = FALSE
)

for (yr in years) {
  long_start <- as.Date(paste0(yr, "-03-01"))
  long_end <- as.Date(paste0(yr, "-06-30"))
  short_start <- as.Date(paste0(yr, "-09-01"))
  short_end <- as.Date(paste0(yr, "-12-31"))
  
  data_year <- rain_data %>% filter(Year == yr)
  
  long_cessation <- detect_cessation_period(data_year, long_start, long_end)
  short_cessation <- detect_cessation_period(data_year, short_start, short_end)
  
  results <- rbind(results, data.frame(
    Year = yr,
    LongRainsCessation = long_cessation,
    ShortRainsCessation = short_cessation
  ))
}

# --- Step 9: Print summary of cessation dates ---
print(results)