SmartCane/r_app/test_overview_maps_aggregation.R

#!/usr/bin/env Rscript

# ==============================================================================
# TEST SCRIPT: Farm-Level Mosaic Aggregation for Overview Maps
# ==============================================================================
# Purpose: Test each step of the aggregation pipeline independently
# ==============================================================================

# Parse arguments
args <- commandArgs(trailingOnly = TRUE)
project_dir <- if (length(args) > 0) args[1] else "aura"
report_date_str <- if (length(args) > 1) args[2] else "2022-12-08"

cat("\n========== Testing Overview Maps Aggregation ==========\n")
cat(paste("Project:", project_dir, "\n"))
cat(paste("Report Date:", report_date_str, "\n\n"))
cat(paste("Project:", project_dir, "\n"))
cat(paste("Report Date:", report_date_str, "\n"))
cat(paste(strrep("═", 80), "\n\n"))

# Load libraries
suppressPackageStartupMessages({
  library(here)
  library(sf)
  library(terra)
  library(tidyverse)
  library(lubridate)
  library(ggspatial)
})

# Load project config
tryCatch({
  source(here::here("r_app", "parameters_project.R"))
  source(here::here("r_app", "00_common_utils.R"))
}, error = function(e) {
  stop("Error loading project utilities: ", e$message)
})

# Set up paths
paths <- setup_project_directories(project_dir)
weekly_CI_mosaic <- paths$weekly_mosaic_dir

# Calculate week/year from report_date
report_date_obj <- as.Date(report_date_str)
current_week <- lubridate::isoweek(report_date_obj)
current_iso_year <- lubridate::isoyear(report_date_obj)

cat(paste(strrep("=", 80), "\n"))
cat(paste("STEP 1: Check Directory Structure\n"))
cat(paste(strrep("=", 80), "\n"))

cat(paste("\nweekly_CI_mosaic path:", weekly_CI_mosaic, "\n"))
cat(paste("Directory exists:", dir.exists(weekly_CI_mosaic), "\n"))

if (!dir.exists(weekly_CI_mosaic)) {
  cat("ERROR: weekly_mosaic directory not found!\n")
  quit(status = 1)
}

# List contents
all_items <- list.files(weekly_CI_mosaic, full.names = FALSE)
cat(paste("\nTotal items in weekly_mosaic/:", length(all_items), "\n"))
cat("First 10 items:\n")
for (i in 1:min(10, length(all_items))) {
  cat(paste("  ", all_items[i], "\n"))
}

# Find field directories
field_dirs <- all_items[
  !grepl("\\.tif$", all_items, ignore.case = TRUE) & 
  dir.exists(file.path(weekly_CI_mosaic, all_items))
]

cat(paste("\nField directories found:", length(field_dirs), "\n"))
if (length(field_dirs) > 0) {
  cat("First 10 field directories:\n")
  for (i in 1:min(10, length(field_dirs))) {
    cat(paste("  ", field_dirs[i], "\n"))
  }
}

cat(paste(strrep("=", 80), "\n"))
cat(paste("STEP 2: Check Weekly Mosaic Files for Target Week\n"))
cat(paste(strrep("=", 80), "\n"))

cat(paste("\nTarget week:", sprintf("%02d", current_week), "\n"))
cat(paste("Target year:", current_iso_year, "\n\n"))

# Check which fields have mosaic files for this week
found_files <- 0
missing_files <- 0

for (field_dir in field_dirs[1:min(10, length(field_dirs))]) {
  expected_file <- paste0("week_", sprintf("%02d", current_week), "_", current_iso_year, ".tif")
  full_path <- file.path(weekly_CI_mosaic, field_dir, expected_file)
  
  if (file.exists(full_path)) {
    cat(paste("  ✓ FOUND:", field_dir, "/", expected_file, "\n"))
    found_files <- found_files + 1
  } else {
    cat(paste("  ✗ MISSING:", field_dir, "/", expected_file, "\n"))
    missing_files <- missing_files + 1
    
    # List what actually exists in this field's directory
    field_path <- file.path(weekly_CI_mosaic, field_dir)
    field_contents <- list.files(field_path, full.names = FALSE)
    if (length(field_contents) > 0) {
      cat(paste("     Available:", paste(field_contents[1:min(3, length(field_contents))], collapse = ", "), "\n"))
    }
  }
}

cat(paste("\nFound: ", found_files, " files | Missing: ", missing_files, "\n"))

if (found_files == 0) {
  cat("\nERROR: No weekly mosaic files found for this week/year combination!\n")
  cat("Check if Script 40 (mosaic_creation) has been run for this week.\n")
  quit(status = 1)
}

cat("\n================================================================================\n")
cat("STEP 3: Load Individual Field Mosaics\n")
cat("================================================================================\n")

# Load all available mosaics
raster_list <- list()
loaded_count <- 0

for (field_dir in field_dirs) {
  full_path <- file.path(weekly_CI_mosaic, field_dir, 
                         paste0("week_", sprintf("%02d", current_week), "_", current_iso_year, ".tif"))
  
  if (file.exists(full_path)) {
    tryCatch({
      r <- terra::rast(full_path)
      raster_list[[field_dir]] <- r
      loaded_count <- loaded_count + 1
      
      if (loaded_count <= 5) {
        cat(paste("  ✓", field_dir, "- Raster loaded\n"))
        cat(paste("    Dimensions:", dim(r)[1], "×", dim(r)[2], "\n"))
        cat(paste("    Bands:", terra::nlyr(r), "\n"))
        cat(paste("    Band names:", paste(names(r), collapse = ", "), "\n"))
        cat(paste("    CRS:", terra::crs(r), "\n\n"))
      }
    }, error = function(e) {
      cat(paste("  ✗", field_dir, "- ERROR loading:", e$message, "\n"))
    })
  }
}

cat(paste("\nSuccessfully loaded:", loaded_count, "field mosaics\n"))

if (loaded_count == 0) {
  cat("\nERROR: Could not load any field mosaics!\n")
  quit(status = 1)
}

cat("\n================================================================================\n")
cat("STEP 4: Test Mosaic Aggregation\n")
cat("================================================================================\n")

cat(paste("\nAttempting to mosaic", length(raster_list), "rasters...\n"))

tryCatch({
  # Create SpatRasterCollection
  cat("  Creating SpatRasterCollection...\n")
  rsrc <- terra::sprc(raster_list)
  cat(paste("  ✓ SpatRasterCollection created with", length(raster_list), "rasters\n\n"))
  
  # Mosaic
  cat("  Mosaicing rasters...\n")
  farm_mosaic <- terra::merge(rsrc)
  cat("  ✓ Mosaic successful!\n\n")
  
  cat(paste("Farm mosaic dimensions:", dim(farm_mosaic)[1], "×", dim(farm_mosaic)[2], "\n"))
  cat(paste("Bands:", terra::nlyr(farm_mosaic), "\n"))
  cat(paste("Band names:", paste(names(farm_mosaic), collapse = ", "), "\n"))
  cat(paste("CRS:", terra::crs(farm_mosaic), "\n"))
  
}, error = function(e) {
  cat(paste("✗ ERROR during mosaicing:", e$message, "\n"))
  quit(status = 1)
})

cat("\n================================================================================\n")
cat("STEP 5: Extract CI Band\n")
cat("================================================================================\n")

tryCatch({
  if ("CI" %in% names(farm_mosaic)) {
    cat("  CI band found by name\n")
    farm_ci <- farm_mosaic[["CI"]]
  } else if (terra::nlyr(farm_mosaic) >= 5) {
    cat("  CI band not named, using band 5\n")
    farm_ci <- farm_mosaic[[5]]
  } else {
    stop("Could not find CI band (expected band 5 or named 'CI')")
  }
  
  cat(paste("  ✓ CI band extracted\n"))
  cat(paste("    Dimensions:", dim(farm_ci)[1], "×", dim(farm_ci)[2], "\n"))
  cat(paste("    Data range:", round(terra::minmax(farm_ci)[1], 2), "to", round(terra::minmax(farm_ci)[2], 2), "\n"))
  cat(paste("    NA values:", sum(is.na(terra::values(farm_ci))), "\n\n"))
  
}, error = function(e) {
  cat(paste("✗ ERROR extracting CI band:", e$message, "\n"))
  quit(status = 1)
})

cat(paste(strrep("=", 80), "\n"))
cat(paste("STEP 6: Load Field Boundaries for Visualization\n"))
cat(paste(strrep("=", 80), "\n"))

tryCatch({
  boundaries_result <- load_field_boundaries(paths$data_dir)
  
  if (is.list(boundaries_result) && "field_boundaries_sf" %in% names(boundaries_result)) {
    field_boundaries_sf <- boundaries_result$field_boundaries_sf
  } else {
    field_boundaries_sf <- boundaries_result
  }
  
  if (nrow(field_boundaries_sf) == 0) {
    stop("No field boundaries loaded")
  }
  
  AllPivots0 <- field_boundaries_sf %>%
    dplyr::filter(!is.na(field), !is.na(sub_field))
  
  cat(paste("  ✓ Field boundaries loaded\n"))
  cat(paste("    Fields:", nrow(AllPivots0), "\n"))
  cat(paste("    CRS:", sf::st_crs(AllPivots0)$epsg, "\n\n"))
  
}, error = function(e) {
  cat(paste("✗ ERROR loading field boundaries:", e$message, "\n"))
  AllPivots0 <- NULL
})

cat("\n================================================================================\n")
cat("STEP 7: Test ggplot Visualization\n")
cat("================================================================================\n")

tryCatch({
  cat("  Reprojecting raster and boundaries to EPSG:4326 for OSM basemap...\n")
  target_crs <- "EPSG:4326"
  farm_ci_ll <- farm_ci
  AllPivots0_ll <- AllPivots0

  if (!terra::is.lonlat(farm_ci)) {
    farm_ci_ll <- terra::project(farm_ci, target_crs, method = "bilinear")
  }
  if (!is.null(AllPivots0)) {
    AllPivots0_ll <- sf::st_transform(AllPivots0, 4326)
  }

  # Ensure boundaries align with raster extent to avoid plotting issues
  sf::sf_use_s2(FALSE)
  if (!is.null(AllPivots0_ll)) {
    AllPivots0_ll <- sf::st_make_valid(AllPivots0_ll)
    crop_bbox_current <- sf::st_as_sfc(sf::st_bbox(terra::ext(farm_ci_ll), crs = 4326))
    AllPivots0_ll <- sf::st_intersection(AllPivots0_ll, crop_bbox_current)
    AllPivots0_ll <- sf::st_collection_extract(AllPivots0_ll, "POLYGON")
  }

  bounds_df <- NULL
  labels_df <- NULL
  if (!is.null(AllPivots0_ll)) {
    bounds_coords <- sf::st_coordinates(AllPivots0_ll)
    bounds_df <- as.data.frame(bounds_coords)
    bounds_df$group <- interaction(bounds_df$L1, bounds_df$L2, drop = TRUE)
    label_pts <- sf::st_point_on_surface(AllPivots0_ll)
    labels_df <- cbind(as.data.frame(sf::st_coordinates(label_pts)), sub_field = label_pts$sub_field)
  }

  cat("  Converting raster to data.frame...\n")
  ci_df <- as.data.frame(farm_ci_ll, xy = TRUE, na.rm = FALSE)
  colnames(ci_df) <- c("x", "y", "ci_value")
  
  cat(paste("  Data.frame dimensions:", nrow(ci_df), "rows ×", ncol(ci_df), "columns\n"))
  cat(paste("  Non-NA pixels:", sum(!is.na(ci_df$ci_value)), "\n\n"))
  
  cat("  Building ggplot map with OSM basemap...\n")
  
  ci_ext <- terra::ext(farm_ci_ll)
  map <- ggplot2::ggplot() +
    ggspatial::annotation_map_tile(
      type = "osm",
      zoom = 14,
      alpha = 0.4
    ) +
    ggplot2::geom_raster(
      data = ci_df,
      ggplot2::aes(x = x, y = y, fill = ci_value)
    ) +
    ggplot2::scale_fill_viridis_c(
      name = "Chlorophyll Index (CI)",
      limits = c(1, 8),
      direction = -1,
      na.value = "transparent",
      oob = scales::squish
    ) +
    ggplot2::coord_sf(
      crs = 4326,
      xlim = c(ci_ext$xmin, ci_ext$xmax),
      ylim = c(ci_ext$ymin, ci_ext$ymax),
      expand = FALSE
    )
  
  if (!is.null(bounds_df)) {
    map4 <- map + ggplot2::geom_path(
      data = bounds_df,
      ggplot2::aes(x = X, y = Y, group = group),
      color = "black",
      linewidth = 0.4
    )
  }

  if (!is.null(labels_df)) {
    map5 <- map4 + ggplot2::geom_text(
      data = labels_df,
      ggplot2::aes(x = X, y = Y, label = sub_field),
      size = 3,
      color = "black"
    )
  }
  
  map6 <- map5 +
    ggspatial::annotation_scale(
      location = "br",
      width_hint = 0.25
    ) +
    ggplot2::theme_void() +
    ggplot2::theme(
      legend.position = "bottom",
      legend.direction = "horizontal",
      plot.title = ggplot2::element_text(hjust = 0.5, size = 12, face = "bold")
    ) +
    ggplot2::labs(
      title = paste("Test: Farm-Level CI Overview - Week", sprintf("%02d", current_week), "of", current_iso_year)
    )
  
  cat("  ✓ Map object created successfully!\n\n")
  
  # Try to save the map
  output_path <- paste0("test_overview_map_", project_dir, "_w", sprintf("%02d", current_week), "_", current_iso_year, ".png")
  cat(paste("  Saving test map to:", output_path, "\n"))
  
  tryCatch({
    ggplot2::ggsave(output_path, map, width = 12, height = 10, dpi = 150)
    cat(paste("  ✓ Map saved successfully!\n"))
  }, error = function(e) {
    cat(paste("  ✗ Could not save map:", e$message, "\n"))
  })
  
}, error = function(e) {
  cat(paste("✗ ERROR in ggplot visualization:", e$message, "\n"))
  cat(paste("   Full error:", deparse(e), "\n"))
  quit(status = 1)
})

cat("\n================================================================================\n")
cat("SUCCESS: All steps completed!\n")
cat("================================================================================\n")
cat("Summary:\n")
cat(paste("  - Loaded", loaded_count, "field mosaics\n"))
cat(paste("  - Created farm-level mosaic\n"))
cat(paste("  - Extracted CI band\n"))
cat(paste("  - Created ggplot visualization with OSM basemap\n"))
cat("\nThe aggregation pipeline is working correctly.\n")
cat("If the report still shows no maps, check the report date/week combination.\n")