# MOSAIC_CREATION_UTILS.R
# ======================
# Utility functions for creating weekly mosaics from daily satellite imagery.
# These functions support cloud cover assessment, date handling, and mosaic creation.

#' Safe logging function that works whether log_message exists or not
#'
#' @param message The message to log
#' @param level The log level (default: "INFO")
#' @return NULL (used for side effects)
#'
safe_log <- function(message, level = "INFO") {
  if (exists("log_message")) {
    log_message(message, level)
  } else {
    if (level %in% c("ERROR", "WARNING")) {
      warning(message)
    } else {
      message(message)
    }
  }
}

#' Generate a sequence of dates for processing
#'
#' @param end_date The end date for the sequence (Date object)
#' @param offset Number of days to look back from end_date
#' @return A list containing week number, year, and a sequence of dates for filtering
#' 
date_list <- function(end_date, offset) {
  # Input validation
  if (!lubridate::is.Date(end_date)) {
    end_date <- as.Date(end_date)
    if (is.na(end_date)) {
      stop("Invalid end_date provided. Expected a Date object or a string convertible to Date.")
    }
  }
  
  offset <- as.numeric(offset)
  if (is.na(offset) || offset < 1) {
    stop("Invalid offset provided. Expected a positive number.")
  }
  
  # Calculate date range
  offset <- offset - 1  # Adjust offset to include end_date
  start_date <- end_date - lubridate::days(offset)
  
  # Extract week and year information
  week <- lubridate::week(start_date)
  year <- lubridate::year(start_date)
  
  # Generate sequence of dates
  days_filter <- seq(from = start_date, to = end_date, by = "day")
  days_filter <- format(days_filter, "%Y-%m-%d")  # Format for consistent filtering
  
  # Log the date range
  safe_log(paste("Date range generated from", start_date, "to", end_date))
  
  return(list(
    "week" = week,
    "year" = year,
    "days_filter" = days_filter,
    "start_date" = start_date,
    "end_date" = end_date
  ))
}

#' Create a weekly mosaic from available VRT files
#'
#' @param dates List from date_list() with date range info
#' @param field_boundaries Field boundaries for image cropping
#' @param daily_vrt_dir Directory containing VRT files
#' @param merged_final_dir Directory with merged final rasters
#' @param output_dir Output directory for weekly mosaics
#' @param file_name_tif Output filename for the mosaic
#' @param create_plots Whether to create visualization plots (default: TRUE)
#' @return The file path of the saved mosaic
#'
create_weekly_mosaic <- function(dates, field_boundaries, daily_vrt_dir, 
                                merged_final_dir, output_dir, file_name_tif,
                                create_plots = TRUE) {
  # Find VRT files for the specified date range
  vrt_list <- find_vrt_files(daily_vrt_dir, dates)
  
  # Find final raster files for fallback
  raster_files_final <- list.files(merged_final_dir, full.names = TRUE, pattern = "\\.tif$")
  
  # Process the mosaic if VRT files are available
  if (length(vrt_list) > 0) {
    safe_log("VRT list created, assessing cloud cover for mosaic creation")
    
    # Calculate cloud cover statistics
    missing_pixels_count <- count_cloud_coverage(vrt_list, field_boundaries)
    
    # Create mosaic based on cloud cover assessment
    mosaic <- create_mosaic(vrt_list, missing_pixels_count, field_boundaries, raster_files_final)
    
  } else {
    safe_log("No VRT files available for the date range, creating empty mosaic", "WARNING")
    
    # Create empty mosaic if no files are available
    if (length(raster_files_final) == 0) {
      stop("No VRT files or final raster files available to create mosaic")
    }
    
    mosaic <- terra::rast(raster_files_final[1]) %>% 
      terra::setValues(0) %>%
      terra::crop(field_boundaries, mask = TRUE)
    
    names(mosaic) <- c("Red", "Green", "Blue", "NIR", "CI")
  }
  
  # Save the mosaic
  file_path <- save_mosaic(mosaic, output_dir, file_name_tif, create_plots)
  
  safe_log(paste("Weekly mosaic processing completed for week", dates$week))
  
  return(file_path)
}

#' Find VRT files within a date range
#'
#' @param vrt_directory Directory containing VRT files
#' @param dates List from date_list() function containing days_filter
#' @return Character vector of VRT file paths
#' 
find_vrt_files <- function(vrt_directory, dates) {
  # Get all VRT files in directory
  vrt_files <- list.files(here::here(vrt_directory), full.names = TRUE)
  
  if (length(vrt_files) == 0) {
    warning("No VRT files found in directory: ", vrt_directory)
    return(character(0))
  }
  
  # Filter files by dates
  vrt_list <- purrr::map(dates$days_filter, ~ vrt_files[grepl(pattern = .x, x = vrt_files)]) %>%
    purrr::compact() %>%
    purrr::flatten_chr()
  
  # Log results
  safe_log(paste("Found", length(vrt_list), "VRT files for the date range"))
  
  return(vrt_list)
}

#' Count missing pixels (clouds) in rasters
#'
#' @param vrt_list List of VRT files to analyze
#' @param field_boundaries Field boundaries vector for masking
#' @return Data frame with cloud coverage statistics
#'
count_cloud_coverage <- function(vrt_list, field_boundaries) {
  if (length(vrt_list) == 0) {
    warning("No VRT files provided for cloud coverage calculation")
    return(NULL)
  }
  
  tryCatch({
    # Calculate total pixel area using the first VRT file
    total_pix_area <- terra::rast(vrt_list[1]) %>% 
      terra::subset(1) %>% 
      terra::setValues(1) %>%
      terra::crop(field_boundaries, mask = TRUE) %>%
      terra::global(fun = "notNA")
    
    # Extract layer 1 from all rasters (for cloud detection)
    layer_5_list <- purrr::map(vrt_list, function(file) {
      terra::rast(file) %>% terra::subset(1)
    }) %>% terra::rast()
    
    # Calculate percentage of missing pixels (clouds)
    missing_pixels_count <- terra::global(layer_5_list, fun = "notNA") %>%
      dplyr::mutate(
        total_pixels = total_pix_area$notNA,
        missing_pixels_percentage = round(100 - ((notNA / total_pix_area$notNA) * 100)),
        thres_5perc = as.integer(missing_pixels_percentage < 5),
        thres_40perc = as.integer(missing_pixels_percentage < 45)
      )
    
    # Log results
    safe_log(paste(
      "Cloud cover assessment completed for", length(vrt_list), "files.",
      sum(missing_pixels_count$thres_5perc), "files with <5% cloud cover,",
      sum(missing_pixels_count$thres_40perc), "files with <45% cloud cover"
    ))
    
    return(missing_pixels_count)
  }, error = function(e) {
    warning("Error in cloud coverage calculation: ", e$message)
    return(NULL)
  })
}

#' Create a mosaic from VRT files based on cloud coverage
#'
#' @param vrt_list List of VRT files to create mosaic from
#' @param missing_pixels_count Cloud coverage statistics from count_cloud_coverage()
#' @param field_boundaries Field boundaries vector for masking (optional)
#' @param raster_files_final List of processed raster files to use as fallback
#' @return A SpatRaster object with the mosaic
#'
create_mosaic <- function(vrt_list, missing_pixels_count, field_boundaries = NULL, raster_files_final = NULL) {
  # If no VRT files, create an empty mosaic
  if (length(vrt_list) == 0) {
    if (length(raster_files_final) == 0 || is.null(field_boundaries)) {
      stop("No VRT files available and no fallback raster files or field boundaries provided")
    }
    
    safe_log("No images available for this period, creating empty mosaic", "WARNING")
    
    x <- terra::rast(raster_files_final[1]) %>% 
      terra::setValues(0) %>%
      terra::crop(field_boundaries, mask = TRUE)
    
    names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
    return(x)
  }
  
  # If missing pixel count was not calculated, use all files
  if (is.null(missing_pixels_count)) {
    safe_log("No cloud coverage data available, using all images", "WARNING")
    
    rsrc <- terra::sprc(vrt_list)
    x <- terra::mosaic(rsrc, fun = "max")
    names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
    return(x)
  }
  
  # Determine best rasters to use based on cloud coverage
  index_5perc <- which(missing_pixels_count$thres_5perc == max(missing_pixels_count$thres_5perc))
  index_40perc <- which(missing_pixels_count$thres_40perc == max(missing_pixels_count$thres_40perc))
  
  # Create mosaic based on available cloud-free images
  if (sum(missing_pixels_count$thres_5perc) > 1) {
    safe_log("Creating max composite from multiple cloud-free images (<5% clouds)")
    
    cloudy_rasters_list <- vrt_list[index_5perc]
    rsrc <- terra::sprc(cloudy_rasters_list)
    x <- terra::mosaic(rsrc, fun = "max")
    
  } else if (sum(missing_pixels_count$thres_5perc) == 1) {
    safe_log("Using single cloud-free image (<5% clouds)")
    
    x <- terra::rast(vrt_list[index_5perc[1]])
    
  } else if (sum(missing_pixels_count$thres_40perc) > 1) {
    safe_log("Creating max composite from partially cloudy images (<40% clouds)", "WARNING")
    
    cloudy_rasters_list <- vrt_list[index_40perc]
    rsrc <- terra::sprc(cloudy_rasters_list)
    x <- terra::mosaic(rsrc, fun = "max")
    
  } else if (sum(missing_pixels_count$thres_40perc) == 1) {
    safe_log("Using single partially cloudy image (<40% clouds)", "WARNING")
    
    x <- terra::rast(vrt_list[index_40perc[1]])
    
  } else {
    safe_log("No cloud-free images available, using all images", "WARNING")
    
    rsrc <- terra::sprc(vrt_list)
    x <- terra::mosaic(rsrc, fun = "max")
  }
  
  # Set consistent layer names
  names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
  return(x)
}

#' Save a mosaic raster to disk
#'
#' @param mosaic_raster A SpatRaster object to save
#' @param output_dir Directory to save the output
#' @param file_name Filename for the output raster
#' @param plot_result Whether to create visualizations (default: FALSE)
#' @return The file path of the saved raster
#'
save_mosaic <- function(mosaic_raster, output_dir, file_name, plot_result = FALSE) {
  # Validate input
  if (is.null(mosaic_raster)) {
    stop("No mosaic raster provided to save")
  }
  
  # Create output directory if it doesn't exist
  dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  
  # Create full file path
  file_path <- here::here(output_dir, file_name)
  
  # Save raster
  terra::writeRaster(mosaic_raster, file_path, overwrite = TRUE)
  
  # Create plots if requested
  if (plot_result) {
    if ("CI" %in% names(mosaic_raster)) {
      terra::plot(mosaic_raster$CI, main = paste("CI map", file_name))
    }
    
    if (all(c("Red", "Green", "Blue") %in% names(mosaic_raster))) {
      terra::plotRGB(mosaic_raster, main = paste("RGB map", file_name))
    }
  }
  
  # Log save completion
  safe_log(paste("Mosaic saved to:", file_path))
  
  return(file_path)
}