From 4cd62ab82e2a59cceeca86698a84cb293ad4e619 Mon Sep 17 00:00:00 2001 From: Timon Date: Wed, 11 Feb 2026 14:32:36 +0100 Subject: [PATCH] Enhance report utility functions and add validation scripts - Updated `create_CI_map` and `create_CI_diff_map` functions to enforce a 1:1 aspect ratio for consistent map sizing. - Modified `ci_plot` function to adjust widths of arranged maps for better layout. - Changed raster merging method in `aggregate_per_field_mosaics_to_farm_level` from `mosaic` to `merge` for improved handling of field data. - Introduced `test_kpi_validation.R` script to validate the structure of KPI RDS files, ensuring expected KPIs are present. - Added `test_overview_maps_aggregation.R` script to test the aggregation pipeline for overview maps, including loading field mosaics, creating a farm-level mosaic, and generating visualizations. --- README.md | 2 + r_app/10_create_per_field_tiffs_utils.R | 37 +- r_app/20_ci_extraction_per_field.R | 2 + r_app/20_ci_extraction_utils.R | 4 + r_app/80_utils_agronomic_support.R | 5 +- ..._CI_report_with_kpis_agronomic_support.Rmd | 489 ++++++++++++++++-- r_app/MANUAL_PIPELINE_RUNNER.R | 4 +- r_app/report_utils.R | 44 +- r_app/test_kpi_validation.R | 155 ++++++ r_app/test_overview_maps_aggregation.R | 371 +++++++++++++ 10 files changed, 1033 insertions(+), 80 deletions(-) create mode 100644 r_app/test_kpi_validation.R create mode 100644 r_app/test_overview_maps_aggregation.R diff --git a/README.md b/README.md index 068ece8..ebada5f 100644 --- a/README.md +++ b/README.md @@ -18,5 +18,7 @@ SmartCane processes 4-band satellite imagery (RGB + NIR) into actionable crop he --- **Last Updated**: February 2026 + **Maintained By**: Resilience BV + **Repository**: https://github.com/TimonWeitkamp/smartcane_experimental_area diff --git a/r_app/10_create_per_field_tiffs_utils.R b/r_app/10_create_per_field_tiffs_utils.R index 38719e5..af490e0 100644 --- a/r_app/10_create_per_field_tiffs_utils.R +++ b/r_app/10_create_per_field_tiffs_utils.R @@ -152,10 +152,10 @@ crop_tiff_to_fields <- function(tif_path, tif_date, fields, output_base_dir) { #' @param field_tiles_dir Character. Target directory for per-field TIFFs #' (e.g., "field_tiles/"). #' @param fields sf object. GeoDataFrame of field boundaries with 'field_name' column. -#' @param field_tiles_ci_dir Character. Optional. Directory where migrated CI-calculated -#' TIFFs are stored. If provided, skips dates -#' already processed and moved to field_tiles_CI/. -#' Default: NULL (process all dates). +#' @param field_tiles_ci_dir Character. Optional. DEPRECATED - no longer used. +#' Kept for backward compatibility but ignored. +#' Script 10 now only checks its own output +#' directory (field_tiles/). #' @param end_date Date. Optional. End date for processing window (YYYY-MM-DD). #' Default: NULL (process all available dates). #' @param offset Integer. Optional. Number of days to look back from end_date. @@ -245,31 +245,10 @@ process_new_merged_tif <- function(merged_tif_dir, field_tiles_dir, fields, fiel for (tif_path in tiff_files) { tif_date <- gsub("\\.tif$", "", basename(tif_path)) - # CHECK 1: Skip if this date was already migrated to field_tiles_CI/ - # (This means Script 20 already processed it and extracted RDS) - if (!is.null(field_tiles_ci_dir) && dir.exists(field_tiles_ci_dir)) { - # Check if ANY field has this date in field_tiles_CI/ - date_migrated <- FALSE - - # Sample check: look for date in field_tiles_CI/*/DATE.tif - sample_field_dirs <- list.dirs(field_tiles_ci_dir, full.names = TRUE, recursive = FALSE) - for (field_dir in sample_field_dirs) { - potential_file <- file.path(field_dir, paste0(tif_date, ".tif")) - if (file.exists(potential_file)) { - date_migrated <- TRUE - break - } - } - - if (date_migrated) { - safe_log(paste("Skipping:", tif_date, "(already migrated and processed by Script 20)"), "INFO") - total_skipped <- total_skipped + 1 - next - } - } - - # CHECK 2: Skip if this date already exists in field_tiles/ - # (means this date has already been processed through Script 10) + # SKIP CHECK: Only check Script 10's own output directory (field_tiles/) + # Do NOT check field_tiles_CI/ — that's Script 20's responsibility. + # Script 10 and Script 20 are independent; Script 10 should not depend on + # Script 20's success/failure. If field_tiles/ is missing, reprocess it. if (dir.exists(field_tiles_dir)) { date_exists_in_field_tiles <- FALSE diff --git a/r_app/20_ci_extraction_per_field.R b/r_app/20_ci_extraction_per_field.R index 701a108..1a85ad0 100644 --- a/r_app/20_ci_extraction_per_field.R +++ b/r_app/20_ci_extraction_per_field.R @@ -150,7 +150,9 @@ main <- function() { ci_raster <- calc_ci_from_raster(raster_4band) # Create 5-band (R, G, B, NIR, CI) + # Explicitly set band names after combining to ensure proper naming five_band <- c(raster_4band, ci_raster) + names(five_band) <- c("Red", "Green", "Blue", "NIR", "CI") # Now process all fields from this single merged TIFF fields_processed_this_date <- 0 diff --git a/r_app/20_ci_extraction_utils.R b/r_app/20_ci_extraction_utils.R index f8a88c5..21fb8dc 100644 --- a/r_app/20_ci_extraction_utils.R +++ b/r_app/20_ci_extraction_utils.R @@ -1013,6 +1013,10 @@ calc_ci_from_raster <- function(raster_obj) { # NDVI ranges -1 to 1 and is different from Chlorophyll Index ci <- nir / green - 1 + # CRITICAL: Explicitly name the CI band before returning + # This ensures proper band naming when combined with other rasters + names(ci) <- "CI" + return(ci) } diff --git a/r_app/80_utils_agronomic_support.R b/r_app/80_utils_agronomic_support.R index 351d0e5..3d92fee 100644 --- a/r_app/80_utils_agronomic_support.R +++ b/r_app/80_utils_agronomic_support.R @@ -543,7 +543,7 @@ create_field_detail_table <- function(field_df, all_kpis, field_boundaries_sf) { by = c("field_idx") ) %>% left_join( - all_kpis$tch_forecasted %>% select(field_idx, tch_forecasted), + all_kpis$tch_forecasted %>% select(field_idx, tch_forecasted, mean_ci), by = c("field_idx") ) %>% left_join( @@ -562,6 +562,7 @@ create_field_detail_table <- function(field_df, all_kpis, field_boundaries_sf) { `Decline Risk` = decline_severity, `Weed Risk` = weed_pressure_risk, `CI Change %` = mean_ci_pct_change, + `Mean CI` = mean_ci, `CV Value` = cv_value ) %>% # Add placeholder columns expected by reporting script (will be populated from other sources) @@ -570,7 +571,7 @@ create_field_detail_table <- function(field_df, all_kpis, field_boundaries_sf) { `Gap Score` = NA_real_ ) %>% select(field_idx, Field, `Field Size (ha)`, `Growth Uniformity`, `Yield Forecast (t/ha)`, - `Gap Score`, `Decline Risk`, `Weed Risk`, `CI Change %`, `CV Value`) + `Gap Score`, `Decline Risk`, `Weed Risk`, `CI Change %`, `Mean CI`, `CV Value`) return(result) } diff --git a/r_app/90_CI_report_with_kpis_agronomic_support.Rmd b/r_app/90_CI_report_with_kpis_agronomic_support.Rmd index 3b6028d..2d59d0d 100644 --- a/r_app/90_CI_report_with_kpis_agronomic_support.Rmd +++ b/r_app/90_CI_report_with_kpis_agronomic_support.Rmd @@ -53,6 +53,7 @@ suppressPackageStartupMessages({ # Visualization library(tmap) # For interactive maps (field boundary visualization) + library(ggspatial) # For basemap tiles and spatial annotations (OSM basemap with ggplot2) # Reporting library(knitr) # For R Markdown document generation (code execution and output) @@ -378,43 +379,101 @@ if (exists("summary_tables") && !is.null(summary_tables) && length(summary_table ## Executive Summary - Key Performance Indicators ```{r combined_kpi_table, echo=FALSE, results='asis'} -# Safely display KPI tables +# Display KPI tables - standardized format with Level, Count, Percent columns if (exists("summary_tables") && !is.null(summary_tables) && length(summary_tables) > 0) { - # Try to combine KPI tables, with fallback if structure is unexpected tryCatch({ - # Build a list of valid dataframes from summary_tables - valid_tables <- list() - - for (kpi_name in names(summary_tables)) { - kpi_df <- summary_tables[[kpi_name]] - - # Skip NULL, empty, or non-dataframe items - if (!is.null(kpi_df) && is.data.frame(kpi_df) && nrow(kpi_df) > 0) { - # Add KPI name as a column if not already present - if (!"KPI" %in% names(kpi_df)) { - display_name <- gsub("_", " ", tools::toTitleCase(gsub("_summary|_data", "", kpi_name))) - kpi_df$KPI <- display_name + # KPI metadata for display + kpi_display_order <- list( + uniformity = list(display = "Field Uniformity", level_col = "Status", count_col = "Field Count"), + area_change = list(display = "Area Change", level_col = "Status", count_col = "Field Count"), + tch_forecast = list(display = "TCH Forecasted", level_col = NULL, count_col = "Fields"), + growth_decline = list(display = "Growth Decline", level_col = "Trend", count_col = "Field Count"), + weed_pressure = list(display = "Weed Presence", level_col = "Risk Level", count_col = "Field Count"), + gap_filling = list(display = "Gap Filling", level_col = NULL, count_col = NULL) + ) + + standardize_kpi <- function(df, level_col, count_col) { + if (is.null(level_col) || !(level_col %in% names(df)) || is.null(count_col) || !(count_col %in% names(df))) { + return(NULL) + } + total <- sum(df[[count_col]], na.rm = TRUE) + total <- ifelse(total == 0, NA_real_, total) + + df %>% + dplyr::transmute( + Level = as.character(.data[[level_col]]), + Count = as.integer(round(as.numeric(.data[[count_col]]))), + Percent = dplyr::if_else( + is.na(total), + NA_real_, + round(Count / total * 100, 1) + ) + ) + } + + combined_df_list <- list() + order_keys <- names(kpi_display_order) + + for (kpi_key in order_keys) { + config <- kpi_display_order[[kpi_key]] + if (!kpi_key %in% names(summary_tables)) next + kpi_df <- summary_tables[[kpi_key]] + if (is.null(kpi_df) || !is.data.frame(kpi_df) || nrow(kpi_df) == 0) next + + kpi_rows <- standardize_kpi(kpi_df, config$level_col, config$count_col) + if (is.null(kpi_rows) && kpi_key %in% c("tch_forecast", "gap_filling")) { + numeric_cols <- names(kpi_df)[vapply(kpi_df, is.numeric, logical(1))] + if (length(numeric_cols) > 0) { + kpi_rows <- tibble::tibble( + Level = numeric_cols, + Count = round(as.numeric(kpi_df[1, numeric_cols]), 2), + Percent = NA_real_ + ) } - valid_tables[[kpi_name]] <- kpi_df + } + + if (!is.null(kpi_rows)) { + kpi_rows$KPI <- config$display + combined_df_list[[length(combined_df_list) + 1]] <- kpi_rows } } - - # Combine all valid tables - if (length(valid_tables) > 0) { - # Use careful bind_rows that handles mismatched columns - combined_df <- dplyr::bind_rows(valid_tables, .id = NULL) - - # Display as flextable - ft <- flextable(combined_df) %>% autofit() - ft + + combined_df <- dplyr::bind_rows(combined_df_list) + + if (nrow(combined_df) > 0) { + combined_df <- combined_df %>% + dplyr::group_by(KPI) %>% + dplyr::mutate( + KPI_display = if_else(dplyr::row_number() == 1, KPI, "") + ) %>% + dplyr::ungroup() %>% + dplyr::select(KPI = KPI_display, Level, Count, Percent) + + ft <- flextable(combined_df) %>% + merge_v(j = "KPI") %>% + autofit() + + kpi_group_sizes <- combined_df %>% + dplyr::group_by(KPI) %>% + dplyr::tally() %>% + dplyr::pull(n) + cum_rows <- cumsum(kpi_group_sizes) + for (i in seq_along(cum_rows)) { + if (i < length(cum_rows)) { + ft <- ft %>% + hline(i = cum_rows[i], border = officer::fp_border(width = 2)) + } + } + + print(ft) } else { - cat("No valid KPI summary tables found.\n") + cat("No valid KPI summary tables found for display.\n") } }, error = function(e) { - safe_log(paste("Error combining KPI tables:", e$message), "WARNING") - cat("KPI summary tables could not be combined for display. Individual KPI sections will be shown below.\n") + safe_log(paste("Error displaying KPI tables:", e$message), "WARNING") + cat("KPI summary tables could not be displayed. Individual KPI sections will be shown below.\n") }) } else { @@ -620,6 +679,372 @@ if (!exists("field_details_table") || is.null(field_details_table)) { } ``` +## Farm-Level Overview Maps + +```{r aggregate_farm_level_rasters, message=FALSE, warning=FALSE, include=FALSE} +# Aggregate per-field weekly mosaics into single farm-level rasters for visualization +# This creates on-the-fly mosaics for current week and historical weeks without saving intermediate files + +tryCatch({ + safe_log("Starting farm-level raster aggregation for overview maps") + + # Helper function to safely aggregate mosaics for a specific week + aggregate_mosaics_safe <- function(week_num, year_num, label) { + tryCatch({ + safe_log(paste("Aggregating mosaics for", label, "(week", week_num, ",", year_num, ")")) + + # Call the utility function from report_utils.R + # This function reads all per-field mosaics and merges them into a single raster + farm_mosaic <- aggregate_per_field_mosaics_to_farm_level( + weekly_mosaic_dir = weekly_CI_mosaic, + target_week = week_num, + target_year = year_num + ) + + if (!is.null(farm_mosaic)) { + safe_log(paste("✓ Successfully aggregated farm mosaic for", label, "")) + return(farm_mosaic) + } else { + safe_log(paste("Warning: Farm mosaic is NULL for", label), "WARNING") + return(NULL) + } + }, error = function(e) { + safe_log(paste("Error aggregating mosaics for", label, ":", e$message), "WARNING") + return(NULL) + }) + } + + # Aggregate mosaics for three weeks: current, week-1, week-3 + farm_mosaic_current <- aggregate_mosaics_safe(current_week, current_iso_year, "current week") + farm_mosaic_minus_1 <- aggregate_mosaics_safe(week_minus_1, week_minus_1_year, "week-1") + farm_mosaic_minus_3 <- aggregate_mosaics_safe(week_minus_3, week_minus_3_year, "week-3") + + # Extract CI band (5th band, or named "CI") from each aggregated mosaic + farm_ci_current <- NULL + farm_ci_minus_1 <- NULL + farm_ci_minus_3 <- NULL + + if (!is.null(farm_mosaic_current)) { + if ("CI" %in% names(farm_mosaic_current)) { + farm_ci_current <- farm_mosaic_current[["CI"]] + } else if (terra::nlyr(farm_mosaic_current) >= 5) { + farm_ci_current <- farm_mosaic_current[[5]] # CI is typically band 5 (R,G,B,NIR,CI) + } + if (!is.null(farm_ci_current)) { + safe_log("✓ Extracted CI band from current week mosaic") + } + } + + if (!is.null(farm_mosaic_minus_1)) { + if ("CI" %in% names(farm_mosaic_minus_1)) { + farm_ci_minus_1 <- farm_mosaic_minus_1[["CI"]] + } else if (terra::nlyr(farm_mosaic_minus_1) >= 5) { + farm_ci_minus_1 <- farm_mosaic_minus_1[[5]] + } + if (!is.null(farm_ci_minus_1)) { + safe_log("✓ Extracted CI band from week-1 mosaic") + } + } + + if (!is.null(farm_mosaic_minus_3)) { + if ("CI" %in% names(farm_mosaic_minus_3)) { + farm_ci_minus_3 <- farm_mosaic_minus_3[["CI"]] + } else if (terra::nlyr(farm_mosaic_minus_3) >= 5) { + farm_ci_minus_3 <- farm_mosaic_minus_3[[5]] + } + if (!is.null(farm_ci_minus_3)) { + safe_log("✓ Extracted CI band from week-3 mosaic") + } + } + + # Calculate difference rasters + farm_ci_diff_week <- NULL + if (!is.null(farm_ci_current) && !is.null(farm_ci_minus_1)) { + farm_ci_diff_week <- farm_ci_current - farm_ci_minus_1 + safe_log("✓ Calculated week-over-week difference raster") + } else { + safe_log("Warning: Cannot calculate week-over-week difference - missing data", "WARNING") + } + + # Reproject rasters and boundaries for ggplot basemap (OSM uses EPSG:4326) + farm_ci_current_ll <- farm_ci_current + farm_ci_diff_week_ll <- farm_ci_diff_week + AllPivots0_ll <- AllPivots0 + target_crs <- "EPSG:4326" + + downsample_raster <- function(r, max_cells = 2000000) { + if (is.null(r)) { + return(NULL) + } + n_cells <- terra::ncell(r) + if (!is.na(n_cells) && n_cells > max_cells) { + fact <- ceiling(sqrt(n_cells / max_cells)) + safe_log(paste("Downsampling raster by factor", fact), "INFO") + return(terra::aggregate(r, fact = fact, fun = mean, na.rm = TRUE)) + } + r + } + + if (!is.null(farm_ci_current) && !terra::is.lonlat(farm_ci_current)) { + farm_ci_current_ll <- terra::project(farm_ci_current, target_crs, method = "bilinear") + } + if (!is.null(farm_ci_diff_week) && !terra::is.lonlat(farm_ci_diff_week)) { + farm_ci_diff_week_ll <- terra::project(farm_ci_diff_week, target_crs, method = "bilinear") + } + if (!is.null(AllPivots0)) { + AllPivots0_ll <- sf::st_transform(AllPivots0, 4326) + } + + farm_ci_current_ll <- downsample_raster(farm_ci_current_ll) + farm_ci_diff_week_ll <- downsample_raster(farm_ci_diff_week_ll) + + # Ensure boundaries align with raster extent to avoid plotting issues + sf::sf_use_s2(FALSE) + if (!is.null(AllPivots0_ll) && !is.null(farm_ci_current_ll)) { + AllPivots0_ll <- sf::st_make_valid(AllPivots0_ll) + crop_bbox_current <- sf::st_as_sfc(sf::st_bbox(terra::ext(farm_ci_current_ll), crs = 4326)) + AllPivots0_ll <- sf::st_intersection(AllPivots0_ll, crop_bbox_current) + AllPivots0_ll <- sf::st_collection_extract(AllPivots0_ll, "POLYGON") + } + + # Prepare boundary and label data frames for ggplot (avoid geom_sf coord resets) + 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) + } + + # Log completion + safe_log("Farm-level raster aggregation complete") + +}, error = function(e) { + safe_log(paste("Error in farm-level raster aggregation:", e$message), "ERROR") +}) +``` + +### Chlorophyll Index (CI) Overview Map - Current Week + +```{r render_farm_ci_map, echo=FALSE, fig.height=5.5, fig.width=6.5, dpi=150, dev='png', message=FALSE, warning=FALSE} +# Create farm-level chlorophyll index map with OpenStreetMap basemap +tryCatch({ + if (!is.null(farm_ci_current_ll)) { + safe_log("Rendering farm-level CI overview map") + + # Convert raster to data.frame for ggplot + ci_df <- as.data.frame(farm_ci_current_ll, xy = TRUE, na.rm = FALSE) + colnames(ci_df) <- c("x", "y", "ci_value") + + # Choose color palette based on colorblind_friendly parameter + if (colorblind_friendly) { + fill_scale <- ggplot2::scale_fill_viridis_c( + name = "Chlorophyll Index (CI)", + limits = c(1, 8), + direction = -1, # Reversed: green=high, yellow/red=low + na.value = "transparent", + oob = scales::squish + ) + } else { + # Use Red-Yellow-Green diverging palette (reversed for intuitive interpretation) + fill_scale <- ggplot2::scale_fill_distiller( + palette = "RdYlGn", + name = "Chlorophyll Index (CI)", + limits = c(1, 8), + direction = 1, # Standard direction for RdYlGn + na.value = "transparent" + ) + } + + # Build the map + ci_ext <- terra::ext(farm_ci_current_ll) + ci_bbox <- sf::st_bbox(c(xmin = ci_ext$xmin, xmax = ci_ext$xmax, ymin = ci_ext$ymin, ymax = ci_ext$ymax), crs = 4326) + + map <- ggplot2::ggplot() + + ggspatial::annotation_map_tile( + type = "osm", + zoom = 14, + alpha = 0.4 + ) + + # Add raster layer with CI values + ggplot2::geom_raster( + data = ci_df, + ggplot2::aes(x = x, y = y, fill = ci_value) + ) + + fill_scale + + 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)) { + map <- 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)) { + map <- map + ggplot2::geom_text( + data = labels_df, + ggplot2::aes(x = X, y = Y, label = sub_field), + size = 3, + color = "black" + ) + } + + map <- map + + # Add scale bar and theme + ggspatial::annotation_scale( + location = "br", + width_hint = 0.25 + ) + + ggplot2::theme_void() + + ggplot2::theme( + legend.position = "bottom", + legend.direction = "horizontal", + legend.title = ggplot2::element_text(size = 10), + legend.text = ggplot2::element_text(size = 9), + plot.title = ggplot2::element_text(hjust = 0.5, size = 12, face = "bold"), + panel.background = ggplot2::element_rect(fill = "white", color = NA) + ) + + ggplot2::labs( + title = paste("Current Week CI Overview - Week", current_week, "of", current_iso_year) + ) + + # Print the map + print(map) + safe_log("✓ Farm-level CI map rendered successfully") + + } else { + safe_log("Farm-level CI raster not available - skipping overview map", "WARNING") + plot(1, type = "n", axes = FALSE, xlab = "", ylab = "") + text(1, 1, "Farm-level CI overview map not available", cex = 1.5, col = "gray") + } +}, error = function(e) { + safe_log(paste("Error rendering farm CI overview map:", e$message), "ERROR") + plot(1, type = "n", axes = FALSE, xlab = "", ylab = "") + text(1, 1, "Error creating CI overview map", cex = 1.5, col = "red") +}) +``` + +### Weekly Chlorophyll Index Difference Map + +```{r render_farm_ci_diff_map, echo=FALSE, fig.height=5.5, fig.width=6.5, dpi=150, dev='png', message=FALSE, warning=FALSE} +# Create farm-level CI difference map (week-over-week change) +tryCatch({ + if (!is.null(farm_ci_diff_week_ll)) { + safe_log("Rendering farm-level CI difference map (week-over-week)") + + # Convert difference raster to data.frame for ggplot + diff_df <- as.data.frame(farm_ci_diff_week_ll, xy = TRUE, na.rm = FALSE) + colnames(diff_df) <- c("x", "y", "ci_diff") + + # Determine color palette based on colorblind_friendly parameter + if (colorblind_friendly) { + # Use plasma for colorblind-friendly diverging visualization + fill_scale <- ggplot2::scale_fill_viridis_c( + name = "CI Change (Week-over-Week)", + option = "plasma", + limits = c(-3, 3), + na.value = "transparent", + oob = scales::squish + ) + } else { + # Use Red-Blue diverging palette (red=decline, blue=increase) + fill_scale <- ggplot2::scale_fill_distiller( + palette = "RdBu", + name = "CI Change (Week-over-Week)", + limits = c(-3, 3), + direction = 1, + na.value = "transparent" + ) + } + + # Build the map + diff_ext <- terra::ext(farm_ci_diff_week_ll) + diff_bbox <- sf::st_bbox(c(xmin = diff_ext$xmin, xmax = diff_ext$xmax, ymin = diff_ext$ymin, ymax = diff_ext$ymax), crs = 4326) + + map <- ggplot2::ggplot() + + ggspatial::annotation_map_tile( + type = "osm", + zoom = 12, + alpha = 0.4 + ) + + # Add raster layer with difference values + ggplot2::geom_raster( + data = diff_df, + ggplot2::aes(x = x, y = y, fill = ci_diff) + ) + + fill_scale + + ggplot2::coord_sf( + crs = 4326, + xlim = c(diff_ext$xmin, diff_ext$xmax), + ylim = c(diff_ext$ymin, diff_ext$ymax), + expand = FALSE + ) + + if (!is.null(bounds_df)) { + map <- 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)) { + map <- map + ggplot2::geom_text( + data = labels_df, + ggplot2::aes(x = X, y = Y, label = sub_field), + size = 3, + color = "black" + ) + } + + map <- map + + # Add scale bar and theme + ggspatial::annotation_scale( + location = "br", + width_hint = 0.25 + ) + + ggplot2::theme_void() + + ggplot2::theme( + legend.position = "bottom", + legend.direction = "horizontal", + legend.title = ggplot2::element_text(size = 10), + legend.text = ggplot2::element_text(size = 9), + plot.title = ggplot2::element_text(hjust = 0.5, size = 12, face = "bold"), + panel.background = ggplot2::element_rect(fill = "white", color = NA) + ) + + ggplot2::labs( + title = paste("Weekly CI Change - Week", current_week, "vs Week", week_minus_1) + ) + + # Print the map + print(map) + safe_log("✓ Farm-level CI difference map rendered successfully") + + } else { + safe_log("Farm-level difference raster not available - skipping difference map", "WARNING") + plot(1, type = "n", axes = FALSE, xlab = "", ylab = "") + text(1, 1, "Week-over-week difference map not available", cex = 1.5, col = "gray") + } +}, error = function(e) { + safe_log(paste("Error rendering farm CI difference map:", e$message), "ERROR") + plot(1, type = "n", axes = FALSE, xlab = "", ylab = "") + text(1, 1, "Error creating CI difference map", cex = 1.5, col = "red") +}) +``` + +\newpage + # Section 2: Field-by-Field Analysis ## Overview of Field-Level Insights @@ -635,7 +1060,7 @@ This section provides detailed, field-specific analyses including chlorophyll in \newpage -```{r generate_field_visualizations, eval=TRUE, fig.height=3.8, fig.width=10, dpi=300, dev='png', message=TRUE, echo=FALSE, warning=TRUE, include=TRUE, results='asis'} +```{r generate_field_visualizations, eval=TRUE, fig.height=3.8, fig.width=6.5, dpi=150, dev='png', message=TRUE, echo=FALSE, warning=TRUE, include=TRUE, results='asis'} # Generate detailed visualizations for each field using purrr::walk tryCatch({ # Prepare merged field list and week/year info @@ -799,7 +1224,7 @@ tryCatch({ }) ``` -```{r generate_subarea_visualizations, echo=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, results='asis', eval=FALSE} +```{r generate_subarea_visualizations, echo=FALSE, fig.height=3.8, fig.width=6.5, dpi=150, message=FALSE, warning=FALSE, results='asis', eval=FALSE} # Alternative visualization grouped by sub-area (disabled by default) tryCatch({ # Group pivots by sub-area @@ -907,12 +1332,12 @@ if (!exists("field_details_table") || is.null(field_details_table)) { # Display the cleaned field table with flextable # Set column widths to fit page (approximately 6.5 inches for 1-inch margins) -col_widths <- c(1.2, 0.9, 1.0, 1.0, 0.8, 0.9, 0.8, 0.7, 0.6) # inches for each column +# Scale widths proportionally: original total = 8.0 inches, scale to 6.2 inches +col_widths <- c(0.97, 0.73, 0.80, 0.80, 0.65, 0.73, 0.65, 0.56, 0.48) # inches for each column ft <- flextable(field_details_clean) %>% set_caption("Detailed Field Performance Summary") %>% - width(width = col_widths) %>% - autofit(add_w = 0, add_h = 0) + width(width = col_widths) ft ``` diff --git a/r_app/MANUAL_PIPELINE_RUNNER.R b/r_app/MANUAL_PIPELINE_RUNNER.R index fa8f69a..b2c20db 100644 --- a/r_app/MANUAL_PIPELINE_RUNNER.R +++ b/r_app/MANUAL_PIPELINE_RUNNER.R @@ -438,8 +438,8 @@ # rmarkdown::render( rmarkdown::render( "r_app/90_CI_report_with_kpis_agronomic_support.Rmd", - params = list(data_dir = "aura", report_date = as.Date("2026-01-01")), - output_file = "SmartCane_Report_agronomic_support_aura_2026-01-01.docx", + params = list(data_dir = "aura", report_date = as.Date("2022-12-08")), + output_file = "SmartCane_Report_agronomic_support_aura_2022-12-08.docx", output_dir = "laravel_app/storage/app/aura/reports" ) # diff --git a/r_app/report_utils.R b/r_app/report_utils.R index da18901..e7e7888 100644 --- a/r_app/report_utils.R +++ b/r_app/report_utils.R @@ -78,7 +78,11 @@ create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = # Add layout configuration to prevent legend rescaling map <- map + tm_layout(legend.position = c("left", "bottom"), legend.outside = FALSE, - inner.margins = 0.05) + inner.margins = 0.05, + asp = 1) # Force 1:1 aspect ratio for consistent sizing + + # Add bounds/view settings for fixed aspect ratio + map <- map + tm_view(asp = 1) # Add borders if requested if (borders) { @@ -146,7 +150,11 @@ create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_lege # Add layout configuration to prevent legend rescaling map <- map + tm_layout(legend.position = c("right", "bottom"), legend.outside = FALSE, - inner.margins = 0.05) + inner.margins = 0.05, + asp = 1) # Force 1:1 aspect ratio for consistent sizing + + # Add bounds/view settings for fixed aspect ratio + map <- map + tm_view(asp = 1) # Add borders if requested if (borders) { @@ -260,17 +268,19 @@ ci_plot <- function(pivotName, week = week, age = age, borders = borders, colorblind = colorblind_friendly) # Create historical maps only if data is available - maps_to_arrange <- list(CImap) - widths_to_use <- c(1) + # Build list with all available maps - order matches original: [m2, m1, current, diff_1w, diff_3w] + # Widths match original hardcoded: c(0.23, 0.18, 0.18, 0.18, 0.23) + maps_to_arrange <- list() + widths_to_use <- c() field_heading_note <- "" - # Try to create 2-week ago map + # Try to create 2-week ago map (legend on left) if (!is.null(singlePivot_m2)) { CImap_m2 <- create_CI_map(singlePivot_m2, AllPivots2, joined_spans2, show_legend = TRUE, legend_is_portrait = TRUE, week = week_minus_2, age = age - 2, borders = borders, colorblind = colorblind_friendly) - maps_to_arrange <- c(list(CImap_m2), maps_to_arrange) - widths_to_use <- c(0.4, widths_to_use) + maps_to_arrange <- c(maps_to_arrange, list(CImap_m2)) + widths_to_use <- c(widths_to_use, 0.24) } # Try to create 1-week ago map @@ -279,25 +289,29 @@ ci_plot <- function(pivotName, show_legend = FALSE, legend_is_portrait = FALSE, week = week_minus_1, age = age - 1, borders = borders, colorblind = colorblind_friendly) maps_to_arrange <- c(maps_to_arrange, list(CImap_m1)) - widths_to_use <- c(widths_to_use, 0.3) + widths_to_use <- c(widths_to_use, 0.17) } + # Always add current week map (center position) + maps_to_arrange <- c(maps_to_arrange, list(CImap)) + widths_to_use <- c(widths_to_use, 0.17) + # Try to create 1-week difference map if (!is.null(abs_CI_last_week)) { CI_max_abs_last_week <- create_CI_diff_map(abs_CI_last_week, AllPivots2, joined_spans2, show_legend = FALSE, legend_is_portrait = FALSE, week_1 = week, week_2 = week_minus_1, age = age, borders = borders, colorblind = colorblind_friendly) maps_to_arrange <- c(maps_to_arrange, list(CI_max_abs_last_week)) - widths_to_use <- c(widths_to_use, 0.3) + widths_to_use <- c(widths_to_use, 0.17) } - # Try to create 3-week difference map + # Try to create 3-week difference map (legend on right) if (!is.null(abs_CI_three_week)) { CI_max_abs_three_week <- create_CI_diff_map(abs_CI_three_week, AllPivots2, joined_spans2, show_legend = TRUE, legend_is_portrait = TRUE, week_1 = week, week_2 = week_minus_3, age = age, borders = borders, colorblind = colorblind_friendly) maps_to_arrange <- c(maps_to_arrange, list(CI_max_abs_three_week)) - widths_to_use <- c(widths_to_use, 0.4) + widths_to_use <- c(widths_to_use, 0.24) } # Normalize widths to sum to 1 @@ -925,19 +939,19 @@ aggregate_per_field_mosaics_to_farm_level <- function( safe_log(paste("Successfully loaded mosaics for", length(raster_list), "fields"), "INFO") - # Create a SpatRasterCollection and mosaic using correct terra syntax + # Create a SpatRasterCollection and merge using correct terra syntax tryCatch({ rsrc <- terra::sprc(raster_list) safe_log(paste("Created SpatRasterCollection with", length(raster_list), "rasters"), "DEBUG") - # Mosaic the rasters - this merges them into a single continuous raster - farm_mosaic <- terra::mosaic(rsrc) + # Merge the rasters into a single continuous raster (no overlap expected between fields) + farm_mosaic <- terra::merge(rsrc) safe_log(paste("Aggregated", length(raster_list), "per-field mosaics into farm-level mosaic"), "INFO") # Verify mosaic was created successfully if (is.null(farm_mosaic)) { - stop("mosaic() returned NULL") + stop("merge() returned NULL") } return(farm_mosaic) diff --git a/r_app/test_kpi_validation.R b/r_app/test_kpi_validation.R new file mode 100644 index 0000000..f59f650 --- /dev/null +++ b/r_app/test_kpi_validation.R @@ -0,0 +1,155 @@ +#!/usr/bin/env Rscript +# Diagnostic script to validate KPI RDS file structure +# Usage: Rscript test_kpi_validation.R [project] [date] + +# Set up arguments +args <- commandArgs(trailingOnly = TRUE) + +if (length(args) < 2) { + cat("Usage: Rscript test_kpi_validation.R [project] [date]\n") + cat("Example: Rscript test_kpi_validation.R aura 2022-11-14\n") + quit(status = 1) +} + +project_dir <- args[1] +report_date <- as.Date(args[2]) + +cat("\n=== KPI RDS Validation ===\n") +cat("Project:", project_dir, "\n") +cat("Date:", report_date, "\n") + +# Load utilities +source("r_app/parameters_project.R") +source("r_app/00_common_utils.R") + +# Set up paths +paths <- setup_project_directories(project_dir) +kpi_data_dir <- paths$kpi_reports_dir + +# Calculate week +current_week <- as.numeric(format(as.Date(report_date), "%V")) +current_year <- as.numeric(format(as.Date(report_date), "%G")) + +kpi_rds_filename <- paste0(project_dir, "_kpi_summary_tables_week", + sprintf("%02d_%d", current_week, current_year), ".rds") +kpi_rds_path <- file.path(kpi_data_dir, kpi_rds_filename) + +cat("\nKPI directory:", kpi_data_dir, "\n") +cat("KPI filename:", kpi_rds_filename, "\n") +cat("Full path:", kpi_rds_path, "\n\n") + +# Check if directory exists +if (!dir.exists(kpi_data_dir)) { + cat("ERROR: KPI directory does not exist!\n") + quit(status = 1) +} + +# List available files +cat("Files in KPI directory:\n") +files <- list.files(kpi_data_dir, pattern = "\\.rds$") +if (length(files) == 0) { + cat(" (none)\n") +} else { + for (f in files) { + cat(" -", f, "\n") + } +} + +# Check if our specific file exists +if (!file.exists(kpi_rds_path)) { + cat("\nWARNING: Expected KPI file not found!\n") + cat("Expected:", kpi_rds_filename, "\n") + quit(status = 1) +} + +cat("\n✓ KPI file found. Loading...\n\n") + +# Load the RDS +loaded_data <- readRDS(kpi_rds_path) + +# Inspect structure +cat("=== RDS Structure ===\n") +cat("Class:", class(loaded_data), "\n") +cat("Length:", length(loaded_data), "\n") +cat("Names:", paste(names(loaded_data), collapse = ", "), "\n\n") + +# Check if new or legacy structure +if (is.list(loaded_data) && "summary_tables" %in% names(loaded_data)) { + cat("✓ New structure detected (has $summary_tables)\n\n") + summary_tables <- loaded_data$summary_tables + + if ("field_details" %in% names(loaded_data)) { + cat("✓ Also has $field_details\n\n") + } +} else { + cat("✓ Legacy structure (direct list of KPI tables)\n\n") + summary_tables <- loaded_data +} + +# Now inspect the summary_tables +cat("=== Available KPI Tables ===\n") +cat("Keys:", paste(names(summary_tables), collapse = ", "), "\n\n") + +# Expected KPIs +expected_kpis <- c( + "uniformity", + "area_change", + "tch_forecasted", + "growth_decline", + "weed_pressure", + "gap_filling" +) + +cat("=== Expected vs Actual ===\n") +for (kpi in expected_kpis) { + # Try both formats + found <- FALSE + actual_key <- NA + + if (kpi %in% names(summary_tables)) { + found <- TRUE + actual_key <- kpi + } else if (paste0(kpi, "_summary") %in% names(summary_tables)) { + found <- TRUE + actual_key <- paste0(kpi, "_summary") + } + + status <- if (found) "✓ FOUND" else "✗ MISSING" + cat(sprintf("%-20s %s", kpi, status)) + if (found) { + cat(" (key: ", actual_key, ")") + } + cat("\n") +} + +cat("\n=== Detailed KPI Contents ===\n") +for (kpi_key in names(summary_tables)) { + kpi_df <- summary_tables[[kpi_key]] + + cat("\n", kpi_key, ":\n", sep="") + cat(" Class:", class(kpi_df), "\n") + cat(" Dimensions:", nrow(kpi_df), "rows ×", ncol(kpi_df), "cols\n") + cat(" Columns:", paste(names(kpi_df), collapse = ", "), "\n") + + if (nrow(kpi_df) > 0) { + cat(" First few rows:\n") + print(head(kpi_df, 3)) + } else { + cat(" (empty dataframe)\n") + } +} + +cat("\n=== Validation Summary ===\n") +missing_count <- sum(!expected_kpis %in% c(names(summary_tables), paste0(expected_kpis, "_summary"))) +if (missing_count == 0) { + cat("✓ All expected KPIs are present!\n") +} else { + cat("✗ Missing", missing_count, "KPI(s):\n") + for (kpi in expected_kpis) { + if (!kpi %in% names(summary_tables) && !paste0(kpi, "_summary") %in% names(summary_tables)) { + cat(" -", kpi, "\n") + } + } +} + +cat("\n") diff --git a/r_app/test_overview_maps_aggregation.R b/r_app/test_overview_maps_aggregation.R new file mode 100644 index 0000000..8236347 --- /dev/null +++ b/r_app/test_overview_maps_aggregation.R @@ -0,0 +1,371 @@ +#!/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")