# MOSAIC_CREATION_TILE_UTILS.R
# ============================
# Tile-based processing utilities for scalable weekly mosaic creation
# 
# STRATEGY: Memory-efficient, scalable approach for large study areas
# Split daily full-extent mosaics into fixed-size tiles (e.g., 5×5 km),
# then process each tile position across all days, and reassemble.
#
# WORKFLOW:
#   Daily full-extent TIFF
#     ↓
#   Split into N×M tiles (based on area size / tile_size)
#   E.g., 50×80 km area with 5 km tiles = 10×16 = 160 tiles
#     ↓
#   For EACH TILE POSITION (1 to 160):
#     - Load that tile from all 7 days
#     - Create MAX composite for that tile
#     - Save tile_NNN_MAX.tif
#     ↓
#   Reassemble all 160 MAX tiles back into full-extent
#     ↓
#   Save weekly full-extent mosaic
#
# KEY BENEFIT: Memory usage ~50 MB per tile (5×5 km), not 1.3 GB (full 50×80 km)
#              Scales automatically: bigger area = more tiles, same memory per tile
#
# TILE_SIZE: Configurable (default 5000 m = 5 km, adjustable via parameter)

#' Simple tile-based processing using terra::makeTiles()
#'
#' Calculates tile grid based on desired tile SIZE (e.g., 5 km), not grid count.
#' This makes it SCALABLE: bigger area = more tiles automatically.
#'
#' @param raster_path Path to raster to split
#' @param output_dir Directory to save tiles
#' @param tile_size_m Tile size in meters (default: 5000 m = 5 km)
#' @return List with tiles (file paths) and metadata
#'
split_raster_into_tiles <- function(raster_path, output_dir, tile_size_m = 5000) {
  # Load raster
  r <- terra::rast(raster_path)
  
  # Calculate how many tiles we need based on raster extent and tile size
  x_range <- terra::ext(r)$xmax - terra::ext(r)$xmin
  y_range <- terra::ext(r)$ymax - terra::ext(r)$ymin
  
  n_tiles_x <- ceiling(x_range / tile_size_m)
  n_tiles_y <- ceiling(y_range / tile_size_m)
  
  safe_log(paste("Splitting into", n_tiles_x, "×", n_tiles_y, "tiles of", 
                 tile_size_m / 1000, "km"))
  
  # Create output directory
  dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  
  # Use terra::makeTiles() - it splits based on n = c(rows, cols)
  tiles <- terra::makeTiles(
    r,
    n = c(n_tiles_y, n_tiles_x),  # rows, cols
    filename = file.path(output_dir, "tile_.tif"),
    overwrite = TRUE
  )
  
  safe_log(paste("Created", length(tiles), "tile files"))
  
  return(list(
    tiles = tiles,
    n_tiles = length(tiles),
    n_tiles_x = n_tiles_x,
    n_tiles_y = n_tiles_y,
    extent = terra::ext(r),
    raster = raster_path
  ))
}


#' Create weekly MAX mosaic using TRUE tile-based processing
#'
#' TILE-BASED WORKFLOW (Memory efficient):
#' 1. Calculate tile grid dynamically from extent and FIXED tile_size (e.g., 5 km)
#' 2. For EACH TILE across all 7 days:
#'    - Load that tile from each daily file (small ~50 MB, not 1.3 GB)
#'    - Create MAX composite for just that tile
#' 3. Reassemble all tiles into final full-extent mosaic
#'
#' SCALABILITY: Fixed 5 km tile size means bigger area = more tiles (automatic scaling)
#'
#' @param dates List from date_list() with days_filter
#' @param merged_final_dir Directory with daily merged full-extent TIFFs
#' @param final_output_dir Directory to store final reassembled mosaic
#' @param file_name_tif Output filename for final mosaic
#' @param tile_size_m Tile size in meters (default: 5000 m = 5 km). Bigger area automatically gets more tiles.
#' @return Path to final reassembled weekly mosaic
#'
create_weekly_mosaic_tiled <- function(dates, merged_final_dir, 
                                       final_output_dir, file_name_tif,
                                       tile_size_m = 5000) {
  
  # Get daily files for this week
  daily_files <- list.files(merged_final_dir, full.names = TRUE, pattern = "\\.tif$")
  daily_files <- daily_files[grepl(paste(dates$days_filter, collapse = "|"), daily_files)]
  
  if (length(daily_files) == 0) {
    safe_log("No daily files found for this week", "ERROR")
    return(NULL)
  }
  
  safe_log(paste("Found", length(daily_files), "daily files for week", dates$week))
  
  # Load first daily file to get extent and calculate tile grid
  safe_log("Step 1: Loading first daily file to establish tile structure")
  first_raster <- terra::rast(daily_files[1])
  
  # Get CRS and convert tile_size_m to appropriate units
  raster_crs <- terra::crs(first_raster)
  
  # If raster is in lat/lon (geographic), convert tile_size_m to degrees
  # 1 degree latitude ≈ 111 km, so 5000 m ≈ 0.045 degrees
  # Check for GEOGCRS (geographic coordinate system) in WKT string
  is_geographic <- grepl("GEOGCRS|longlat|geographic|ANGLEUNIT.*degree", 
                         as.character(raster_crs), ignore.case = TRUE)
  
  if (is_geographic) {
    # Geographic CRS - convert meters to degrees
    # At equator: 1 degree ≈ 111 km
    tile_size_degrees <- tile_size_m / 111000  # 111 km per degree
    safe_log(paste("Raster is in geographic CRS (lat/lon). Converting", tile_size_m / 1000, 
                   "km to", round(tile_size_degrees, 4), "degrees"))
  } else {
    # Projected CRS - use meters directly
    tile_size_degrees <- tile_size_m
    safe_log(paste("Raster is in projected CRS. Using", tile_size_m / 1000, "km directly"))
  }
  
  # Calculate n_tiles dynamically from extent and tile_size
  x_range <- terra::ext(first_raster)$xmax - terra::ext(first_raster)$xmin
  y_range <- terra::ext(first_raster)$ymax - terra::ext(first_raster)$ymin
  
  n_tiles_x <- ceiling(x_range / tile_size_degrees)
  n_tiles_y <- ceiling(y_range / tile_size_degrees)
  n_tiles_total <- n_tiles_x * n_tiles_y
  
  safe_log(paste("Step 2: Creating tile grid:", tile_size_m / 1000, "km tiles"))
  safe_log(paste("  Extent:", round(x_range, 4), "° ×", round(y_range, 4), "°"))
  safe_log(paste("  Grid:", n_tiles_y, "rows ×", n_tiles_x, "columns =", n_tiles_total, "tiles"))
  
  # Calculate tile extents mathematically (no need to save temp files)
  extent <- terra::ext(first_raster)
  x_min <- extent$xmin
  y_min <- extent$ymin
  
  # Create list of tile extents
  tile_extents <- list()
  tile_idx <- 0
  
  for (row in 1:n_tiles_y) {
    for (col in 1:n_tiles_x) {
      tile_idx <- tile_idx + 1
      
      # Calculate this tile's bounds
      tile_xmin <- x_min + (col - 1) * tile_size_degrees
      tile_xmax <- min(tile_xmin + tile_size_degrees, extent$xmax)
      tile_ymin <- y_min + (row - 1) * tile_size_degrees
      tile_ymax <- min(tile_ymin + tile_size_degrees, extent$ymax)
      
      tile_extents[[tile_idx]] <- terra::ext(tile_xmin, tile_xmax, tile_ymin, tile_ymax)
    }
  }
  
  safe_log(paste("Calculated", length(tile_extents), "tile extents"))
  
  # Save tiles to Laravel storage directory
  tile_storage_dir <- file.path("laravel_app", "storage", "app", "angata", "daily_tiles")
  dir.create(tile_storage_dir, recursive = TRUE, showWarnings = FALSE)
  
  # For each tile, load from all daily files and create MAX
  safe_log("Step 3: Processing tiles (load per-tile across all days, create MAX for each)")
  tile_files_list <- list()
  
  for (tile_idx in seq_along(tile_extents)) {
    if (tile_idx %% max(1, ceiling(n_tiles_total / 10)) == 0 || tile_idx == 1) {
      safe_log(paste("  Processing tile", tile_idx, "of", n_tiles_total))
    }
    
    # Get this tile's extent
    tile_extent <- tile_extents[[tile_idx]]
    
    # Load and crop all daily files to this tile extent
    daily_tile_data <- list()
    for (file_idx in seq_along(daily_files)) {
      tryCatch({
        r <- terra::rast(daily_files[file_idx])
        cropped <- terra::crop(r, tile_extent, snap = "in")
        daily_tile_data[[file_idx]] <- cropped
      }, error = function(e) {
        safe_log(paste("Warning: Could not crop tile", tile_idx, "from", 
                       basename(daily_files[file_idx])), "WARNING")
      })
    }
    
    if (length(daily_tile_data) == 0) {
      safe_log(paste("No valid data for tile", tile_idx), "WARNING")
      next
    }
    
    # Create MAX composite for this tile
    if (length(daily_tile_data) == 1) {
      tile_max <- daily_tile_data[[1]]
    } else {
      collection <- terra::sprc(daily_tile_data)
      tile_max <- terra::mosaic(collection, fun = "max")
    }
    
    # Save tile to disk (keeps memory low)
    tile_file <- file.path(tile_storage_dir, sprintf("tile_%04d.tif", tile_idx))
    terra::writeRaster(tile_max, tile_file, overwrite = TRUE)
    tile_files_list[[tile_idx]] <- tile_file
  }
  
  if (length(tile_files_list) == 0) {
    safe_log("No tiles processed successfully", "ERROR")
    return(NULL)
  }
  
  safe_log(paste("Step 4: Reassembling", length(tile_files_list), "tiles from disk into full-extent mosaic"))
  
  # Load all tile files and reassemble
  tile_rasters <- lapply(tile_files_list, terra::rast)
  collection <- terra::sprc(tile_rasters)
  final_mosaic <- terra::mosaic(collection, fun = "first")
  
  safe_log("Step 5: Saving final reassembled mosaic")
  
  # Save
  dir.create(final_output_dir, recursive = TRUE, showWarnings = FALSE)
  final_file_path <- file.path(final_output_dir, file_name_tif)
  
  tryCatch({
    terra::writeRaster(final_mosaic, final_file_path, overwrite = TRUE)
    safe_log(paste("✓ Weekly mosaic saved to:", final_file_path))
  }, error = function(e) {
    safe_log(paste("Error saving mosaic:", e$message), "ERROR")
    return(NULL)
  })
  
  # Cleanup temporary tile files
  unlink(tile_storage_dir, recursive = TRUE)
  safe_log("Cleaned up temporary tile files")
  
  return(final_file_path)
}

#' Load tile MAX rasters for a specific week (for per-tile analysis)
#'
#' @param week Week number
#' @param tile_dir Directory containing tile mosaics (week subdirectories)
#' @return List of tile rasters with tile_id metadata
#'
load_tile_mosaics <- function(week, tile_dir) {
  week_dir <- file.path(tile_dir, paste0("week_", week))
  
  if (!dir.exists(week_dir)) {
    warning(paste("Tile directory not found:", week_dir))
    return(NULL)
  }
  
  # Load all tile files
  tile_files <- list.files(week_dir, pattern = "^tile_.*\\.tif$", full.names = TRUE)
  
  if (length(tile_files) == 0) {
    warning("No tile files found in:", week_dir)
    return(NULL)
  }
  
  # Sort by tile ID
  tile_numbers <- as.numeric(gsub(".*tile_(\\d+).*", "\\1", tile_files))
  tile_files <- tile_files[order(tile_numbers)]
  
  # Load rasters
  tiles_list <- lapply(tile_files, terra::rast)
  names(tiles_list) <- sprintf("tile_%03d", sort(tile_numbers))
  
  safe_log(paste("Loaded", length(tiles_list), "tile mosaics for week", week))
  
  return(tiles_list)
}