# REPORT_UTILS.R # ============= # Utility functions for generating SmartCane reports with visualizations. # These functions support the creation of maps, charts and report elements # for the CI_report_dashboard_planet.Rmd document. #' Creates a sub-chunk for use within RMarkdown documents #' #' @param g A ggplot object to render in the sub-chunk #' @param fig_height Height of the figure in inches #' @param fig_width Width of the figure in inches #' @return NULL (writes chunk directly to output) #' subchunkify <- function(g, fig_height=7, fig_width=5) { g_deparsed <- paste0(deparse( function() {g} ), collapse = '') sub_chunk <- paste0(" `","``{r sub_chunk_", floor(runif(1) * 10000), ", fig.height=", fig_height, ", fig.width=", fig_width, ", dpi=300, dev='png', out.width='100%', echo=FALSE}", "\n(", g_deparsed , ")()", "\n`","`` ") cat(knitr::knit(text = knitr::knit_expand(text = sub_chunk), quiet = TRUE)) } #' Creates a Chlorophyll Index map for a pivot #' #' @param pivot_raster The raster data containing CI values #' @param pivot_shape The shape of the pivot field #' @param pivot_spans Additional boundary data for the field #' @param show_legend Whether to show the legend (default: FALSE) #' @param legend_is_portrait Whether to show the legend in portrait orientation (default: FALSE) #' @param legend_position Position for the legend when shown: "left", "right", "top", "bottom" (default: "bottom") #' @param week Week number to display in the title #' @param age Age of the crop in weeks #' @param borders Whether to display field borders (default: FALSE) #' @return A tmap object with the CI map #' create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, legend_position = "bottom", week, age, borders = FALSE, colorblind = FALSE){ # Input validation if (missing(pivot_raster) || is.null(pivot_raster)) { stop("pivot_raster is required") } if (missing(pivot_shape) || is.null(pivot_shape)) { stop("pivot_shape is required") } if (missing(pivot_spans) || is.null(pivot_spans)) { stop("pivot_spans is required") } if (missing(week) || is.null(week)) { stop("week parameter is required") } if (missing(age) || is.null(age)) { stop("age parameter is required") } # Choose palette based on colorblind flag palette <- if (colorblind) "viridis" else "brewer.rd_yl_gn" # Create the base map map <- tm_shape(pivot_raster, unit = "m") # Add raster with continuous spectrum (fixed scale 8-1 for consistent comparison, reversed) map <- map + tm_raster( "CI", col.scale = tm_scale_continuous( values = palette, limits = c(1, 8), ticks = seq(1, 8, by = 1), outliers.trunc = c(TRUE, TRUE) ), col.legend = tm_legend( title = "CI", orientation = if (legend_is_portrait) "portrait" else "landscape", show = show_legend, position = if (show_legend) tm_pos_out(legend_position, "center") else c("left", "bottom"), reverse = TRUE ) ) # Add layout elements map <- map + tm_layout( main.title = paste0("Max CI week ", week,"\n", age, " weeks (", age * 7, " days) old"), main.title.size = 0.7, #legend.height = 0.85, # Constrain vertical legend height to not exceed map asp = 1 # Fixed aspect ratio ) # Add borders if requested if (borders) { map <- map + tm_shape(pivot_shape) + tm_borders(lwd = 3) + tm_text("sub_field", size = 1/2) + tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha = 0.5) } return(map) } #' Creates a Chlorophyll Index difference map between two weeks #' #' @param pivot_raster The raster data containing CI difference values #' @param pivot_shape The shape of the pivot field #' @param pivot_spans Additional boundary data for the field #' @param show_legend Whether to show the legend (default: FALSE) #' @param legend_is_portrait Whether to show the legend in portrait orientation (default: FALSE) #' @param legend_position Position for the legend when shown: "left", "right", "top", "bottom" (default: "bottom") #' @param week_1 First week number for comparison #' @param week_2 Second week number for comparison #' @param age Age of the crop in weeks #' @param borders Whether to display field borders (default: TRUE) #' @return A tmap object with the CI difference map #' create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, legend_position = "bottom", week_1, week_2, age, borders = TRUE, colorblind = FALSE){ # Input validation if (missing(pivot_raster) || is.null(pivot_raster)) { stop("pivot_raster is required") } if (missing(pivot_shape) || is.null(pivot_shape)) { stop("pivot_shape is required") } if (missing(pivot_spans) || is.null(pivot_spans)) { stop("pivot_spans is required") } if (missing(week_1) || is.null(week_1) || missing(week_2) || is.null(week_2)) { stop("week_1 and week_2 parameters are required") } if (missing(age) || is.null(age)) { stop("age parameter is required") } # Choose palette based on colorblind flag palette <- if (colorblind) "plasma" else "brewer.rd_yl_gn" # Create the base map map <- tm_shape(pivot_raster, unit = "m") # Add raster with continuous spectrum (centered at 0 for difference maps, fixed scale, reversed) map <- map + tm_raster( "CI", col.scale = tm_scale_continuous( values = palette, limits = c(-3, 3), ticks = seq(-3, 3, by = 1), midpoint = 0, outliers.trunc = c(TRUE, TRUE) ), col.legend = tm_legend( title = "CI diff.", orientation = if (legend_is_portrait) "portrait" else "landscape", show = show_legend, position = if (show_legend) tm_pos_out(legend_position, "center") else c("left", "bottom"), reverse = TRUE ) ) # Add layout elements map <- map + tm_layout( main.title = paste0("CI change week ", week_1, " - week ", week_2, "\n", age, " weeks (", age * 7, " days) old"), main.title.size = 0.7, #legend.height = 0.85, # Constrain vertical legend height to not exceed map asp = 1 # Fixed aspect ratio ) # Add borders if requested if (borders) { map <- map + tm_shape(pivot_shape) + tm_borders(lwd = 3) + tm_text("sub_field", size = 1/2) + tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha = 0.5) } return(map) } #' Creates a visualization of CI data for a specific pivot field #' #' @param pivotName The name or ID of the pivot field to visualize #' @param field_boundaries Field boundaries spatial data (sf object) #' @param current_ci Current week's Chlorophyll Index raster #' @param ci_minus_1 Previous week's Chlorophyll Index raster #' @param ci_minus_2 Two weeks ago Chlorophyll Index raster #' @param last_week_diff Difference raster between current and last week #' @param three_week_diff Difference raster between current and three weeks ago #' @param harvesting_data Data frame containing field harvesting/planting information #' @param week Current week number #' @param week_minus_1 Previous week number #' @param week_minus_2 Two weeks ago week number #' @param week_minus_3 Three weeks ago week number #' @param borders Whether to display field borders (default: TRUE) #' @param colorblind_friendly Whether to use colorblind-friendly color schemes (default: FALSE) #' @return NULL (adds output directly to R Markdown document) #' ci_plot <- function(pivotName, field_boundaries = AllPivots0, current_ci = CI, ci_minus_1 = CI_m1, ci_minus_2 = CI_m2, last_week_diff = last_week_dif_raster_abs, three_week_diff = three_week_dif_raster_abs, harvesting_data = harvesting_data, week = week, week_minus_1 = week_minus_1, week_minus_2 = week_minus_2, week_minus_3 = week_minus_3, borders = TRUE, colorblind_friendly = FALSE){ # Input validation if (missing(pivotName) || is.null(pivotName) || pivotName == "") { stop("pivotName is required") } if (missing(field_boundaries) || is.null(field_boundaries)) { stop("field_boundaries is required") } if (missing(current_ci) || is.null(current_ci)) { stop("current_ci is required") } # Note: ci_minus_1, ci_minus_2, last_week_diff, three_week_diff are now optional # (may be NULL if historical data is not available for early seasons) if (missing(harvesting_data) || is.null(harvesting_data)) { stop("harvesting_data is required") } # Warn if critical rasters are missing if (is.null(ci_minus_1) || is.null(ci_minus_2)) { safe_log(paste("Warning: Historical CI data missing for field", pivotName, "- will show current week only"), "WARNING") } if (is.null(last_week_diff) || is.null(three_week_diff)) { safe_log(paste("Warning: Difference rasters missing for field", pivotName, "- difference maps skipped"), "WARNING") } # Extract pivot shape and age data tryCatch({ pivotShape <- field_boundaries %>% dplyr::filter(field %in% pivotName) %>% sf::st_transform(terra::crs(current_ci)) age <- harvesting_data %>% dplyr::filter(field %in% pivotName) %>% sort("year") %>% tail(., 1) %>% dplyr::select(age) %>% unique() %>% pull() %>% round() # Filter for the specific pivot AllPivots2 <- field_boundaries %>% dplyr::filter(field %in% pivotName) # Per-field mosaics are already cropped to field boundaries, so use directly without cropping singlePivot <- current_ci singlePivot_m1 <- ci_minus_1 singlePivot_m2 <- ci_minus_2 # Use difference maps directly (already field-specific) abs_CI_last_week <- last_week_diff abs_CI_three_week <- three_week_diff # Get planting date planting_date <- harvesting_data %>% dplyr::filter(field %in% pivotName) %>% ungroup() %>% dplyr::select(season_start) %>% unique() # Create spans for borders joined_spans2 <- field_boundaries %>% sf::st_transform(sf::st_crs(pivotShape)) %>% dplyr::filter(field %in% pivotName) # Create maps conditionally based on data availability # Always create current week map (required) CImap <- create_CI_map(singlePivot, AllPivots2, joined_spans2, show_legend = FALSE, legend_is_portrait = FALSE, week = week, age = age, borders = borders, colorblind = colorblind_friendly) # Create historical maps only if data is available # Build list with all available maps - order matches original: [m2, m1, current, diff_1w, diff_3w] maps_to_arrange <- list() field_heading_note <- "" # 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, legend_position = "left", week = week_minus_2, age = age - 2, borders = borders, colorblind = colorblind_friendly) maps_to_arrange <- c(maps_to_arrange, list(CImap_m2)) } # Try to create 1-week ago map if (!is.null(singlePivot_m1)) { CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2, joined_spans2, 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)) } # Always add current week map (center position) maps_to_arrange <- c(maps_to_arrange, list(CImap)) # 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)) } # 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, legend_position = "right", 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)) } # Add note if historical data is limited if (length(maps_to_arrange) == 1) { field_heading_note <- " (Current week only - historical data not yet available)" } else if (length(maps_to_arrange) < 5) { field_heading_note <- " (Limited historical data)" } # Arrange the maps in a row with more width for first and last (for legends) # Give maps with legends (1st and 5th) more space: 23%, middle maps get 18% each widths <- if (length(maps_to_arrange) == 5) { c(0.23, 0.18, 0.18, 0.18, 0.23) } else if (length(maps_to_arrange) == 4) { c(0.25, 0.25, 0.25, 0.25) # Equal if only 4 maps } else if (length(maps_to_arrange) == 3) { c(0.33, 0.33, 0.34) # Equal if only 3 maps } else if (length(maps_to_arrange) == 2) { c(0.5, 0.5) # Equal if only 2 maps } else { NULL # Single map or other cases } tst <- do.call(tmap_arrange, c(maps_to_arrange, list(nrow = 1, widths = widths))) # Output heading and map to R Markdown age_months <- round(age / 4.348, 1) cat(paste("## Field", pivotName, "-", age, "weeks/", age_months, "months after planting/harvest", field_heading_note, "\n\n")) print(tst) }, error = function(e) { safe_log(paste("Error creating CI plot for pivot", pivotName, ":", e$message), "ERROR") cat(paste("# Field", pivotName, "- Error creating visualization", "\n\n")) cat(paste("Error:", e$message, "\n\n")) }) } #' Creates a plot showing Chlorophyll Index data over time for a pivot field #' #' @param pivotName The name or ID of the pivot field to visualize #' @param ci_quadrant_data Data frame containing CI quadrant data with field, sub_field, Date, DAH, cumulative_CI, value and season columns #' @param plot_type Type of plot to generate: "absolute", "cumulative", or "both" #' @param facet_on Whether to facet the plot by season (TRUE) or overlay all seasons (FALSE) #' @param x_unit Unit for x-axis: "days" for DAH or "weeks" for week number (default: "days") #' @param colorblind_friendly Whether to use colorblind-friendly color schemes (default: FALSE) #' @param show_benchmarks Whether to show historical benchmark lines (default: FALSE) #' @param estate_name Name of the estate for benchmark calculation (required if show_benchmarks = TRUE) #' @param benchmark_percentiles Vector of percentiles for benchmarks (default: c(10, 50, 90)) #' @return NULL (adds output directly to R Markdown document) #' cum_ci_plot <- function(pivotName, ci_quadrant_data = CI_quadrant, plot_type = "absolute", facet_on = FALSE, x_unit = "days", colorblind_friendly = FALSE, show_benchmarks = FALSE, estate_name = NULL, benchmark_percentiles = c(10, 50, 90), benchmark_data = NULL) { # Input validation if (missing(pivotName) || is.null(pivotName) || pivotName == "") { stop("pivotName is required") } if (missing(ci_quadrant_data) || is.null(ci_quadrant_data)) { stop("ci_quadrant_data is required") } if (!plot_type %in% c("absolute", "cumulative", "both")) { stop("plot_type must be one of: 'absolute', 'cumulative', or 'both'") } if (!x_unit %in% c("days", "weeks")) { stop("x_unit must be either 'days' or 'weeks'") } # Filter data for the specified pivot tryCatch({ data_ci <- ci_quadrant_data %>% dplyr::filter(field == pivotName) if (nrow(data_ci) == 0) { safe_log(paste("No CI data found for field", pivotName), "WARNING") return(cum_ci_plot2(pivotName)) # Use fallback function when no data is available } # Process data data_ci2 <- data_ci %>% dplyr::mutate(CI_rate = cumulative_CI / DAH, week = lubridate::week(Date)) %>% dplyr::group_by(field) %>% dplyr::mutate(mean_CIrate_rolling_10_days = zoo::rollapplyr(CI_rate, width = 10, FUN = mean, partial = TRUE), mean_rolling_10_days = zoo::rollapplyr(value, width = 10, FUN = mean, partial = TRUE)) data_ci2 <- data_ci2 %>% dplyr::mutate(season = as.factor(season)) # Compute benchmarks if requested and not provided if (show_benchmarks && is.null(benchmark_data)) { benchmark_data <- compute_ci_benchmarks(ci_quadrant_data, estate_name, benchmark_percentiles) } # Prepare benchmark data for plotting if available if (!is.null(benchmark_data)) { benchmark_data <- benchmark_data %>% dplyr::mutate( ci_type_label = case_when( ci_type == "value" ~ "10-Day Rolling Mean CI", ci_type == "cumulative_CI" ~ "Cumulative CI", TRUE ~ ci_type ), benchmark_label = paste0(percentile, "th Percentile") ) safe_log("Benchmark data prepared for plotting", "INFO") } else if (show_benchmarks) { safe_log("No benchmark data available", "WARNING") } data_ci3 <- tidyr::pivot_longer( data_ci2, cols = c("mean_rolling_10_days", "cumulative_CI"), names_to = "ci_type", # This column will say "mean_rolling_10_days" or "cumulative_CI" values_to = "ci_value" # This column will have the numeric values ) # Prepare date information by season date_preparation_perfect_pivot <- data_ci2 %>% dplyr::group_by(season) %>% dplyr::summarise(min_date = min(Date), max_date = max(Date), days = max_date - min_date) # Get the 3 most recent seasons unique_seasons <- sort(unique(date_preparation_perfect_pivot$season), decreasing = TRUE)[1:3] latest_season <- unique_seasons[1] # Identify the latest season # Create plotting function that uses data_ci3 and filters by ci_type create_plot <- function(ci_type_filter, y_label, title_suffix) { # Filter data based on ci_type plot_data <- data_ci3 %>% dplyr::filter(season %in% unique_seasons, ci_type == ci_type_filter) %>% dplyr::mutate(is_latest = season == latest_season) # Flag for latest season # Determine x-axis variable based on x_unit parameter x_var <- if (x_unit == "days") { if (facet_on) "Date" else "DAH" } else { "week" } x_label <- switch(x_unit, "days" = if (facet_on) "Date" else "Age of Crop (Days)", "weeks" = "Week Number") # Calculate dynamic max values for breaks max_dah <- max(plot_data$DAH, na.rm = TRUE) + 20 max_week <- max(as.numeric(plot_data$week), na.rm = TRUE) + ceiling(20 / 7) # Create plot with either facets by season or overlay by DAH/week if (facet_on) { g <- ggplot2::ggplot(data = plot_data) + ggplot2::facet_wrap(~season, scales = "free_x") + ggplot2::geom_line( ggplot2::aes( x = .data[[x_var]], y = .data[["ci_value"]], col = .data[["sub_field"]], group = .data[["sub_field"]] ) ) + ggplot2::labs(title = paste("Plot of", y_label), color = "Field Name", y = y_label, x = x_label) + ggplot2::scale_x_date(date_breaks = "1 month", date_labels = "%m-%Y", sec.axis = ggplot2::sec_axis(~ ., name = "Age in Months", breaks = scales::breaks_pretty(), labels = function(x) round(as.numeric(x - min(x)) / 30.44, 1))) + ggplot2::theme_minimal() + ggplot2::theme(axis.text.x = ggplot2::element_text(hjust = 0.5), axis.text.x.top = ggplot2::element_text(hjust = 0.5), axis.title.x.top = ggplot2::element_text(size = 8), legend.justification = c(1, 0), legend.position = "inside", legend.position.inside = c(1, 0), legend.title = ggplot2::element_text(size = 8), legend.text = ggplot2::element_text(size = 8)) + ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = TRUE)) } else { # Choose color palette based on colorblind_friendly flag color_scale <- if (colorblind_friendly) { ggplot2::scale_color_brewer(type = "qual", palette = "Set2") } else { ggplot2::scale_color_discrete() } g <- ggplot2::ggplot(data = plot_data) + # Add benchmark lines first (behind season lines) { if (!is.null(benchmark_data) && ci_type_filter %in% benchmark_data$ci_type) { benchmark_subset <- benchmark_data %>% dplyr::filter(ci_type == ci_type_filter) %>% dplyr::mutate( benchmark_x = if (x_var == "DAH") { DAH } else if (x_var == "week") { DAH / 7 # Approximate conversion } else { DAH # For Date, use DAH as is (may not align perfectly) } ) ggplot2::geom_smooth( data = benchmark_subset, ggplot2::aes( x = .data[["benchmark_x"]], y = .data[["benchmark_value"]], group = factor(.data[["percentile"]]) ), color = "gray70", linewidth = 0.5, se = FALSE, inherit.aes = FALSE ) } } + # Plot older seasons with lighter lines ggplot2::geom_line( data = plot_data %>% dplyr::filter(!is_latest), ggplot2::aes( x = .data[[x_var]], y = .data[["ci_value"]], col = .data[["season"]], group = .data[["season"]] ), linewidth = 0.7, alpha = 0.4 ) + # Plot latest season with thicker, more prominent line ggplot2::geom_line( data = plot_data %>% dplyr::filter(is_latest), ggplot2::aes( x = .data[[x_var]], y = .data[["ci_value"]], col = .data[["season"]], group = .data[["season"]] ), linewidth = 1.5, alpha = 1 ) + ggplot2::labs(title = paste("Plot of", y_label, "for Field", pivotName, title_suffix), color = "Season", y = y_label, x = x_label) + color_scale + { if (x_var == "DAH") { ggplot2::scale_x_continuous(breaks = seq(0, 450, by = 50), sec.axis = ggplot2::sec_axis(~ . / 30.44, name = "Age in Months", breaks = seq(0, 14, by = 1))) } else if (x_var == "week") { ggplot2::scale_x_continuous(breaks = seq(0, max_week, by = 4), sec.axis = ggplot2::sec_axis(~ . / 4.348, name = "Age in Months", breaks = seq(0, 14, by = 1))) } } + ggplot2::theme_minimal() + ggplot2::theme(axis.text.x = ggplot2::element_text(hjust = 0.5), axis.text.x.top = ggplot2::element_text(hjust = 0.5), axis.title.x.top = ggplot2::element_text(size = 8), legend.justification = c(1, 0), legend.position = "inside", legend.position.inside = c(1, 0), legend.title = ggplot2::element_text(size = 8), legend.text = ggplot2::element_text(size = 8)) + ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = TRUE)) } # Add y-axis limits for absolute CI (10-day rolling mean) to fix scale at 0-7 if (ci_type_filter == "mean_rolling_10_days") { g <- g + ggplot2::ylim(0, 7) } return(g) } # Generate plots based on plot_type if (plot_type == "absolute") { g <- create_plot("mean_rolling_10_days", "10-Day Rolling Mean CI", "") subchunkify(g, 2.8, 10) } else if (plot_type == "cumulative") { g <- create_plot("cumulative_CI", "Cumulative CI", "") subchunkify(g, 2.8, 10) } else if (plot_type == "both") { # Create faceted plot with both CI types using pivot_longer approach plot_data_both <- data_ci3 %>% dplyr::filter(season %in% unique_seasons) %>% dplyr::mutate( ci_type_label = case_when( ci_type == "mean_rolling_10_days" ~ "10-Day Rolling Mean CI", ci_type == "cumulative_CI" ~ "Cumulative CI", TRUE ~ ci_type ), is_latest = season == latest_season # Flag for latest season ) # Determine x-axis variable based on x_unit parameter x_var <- if (x_unit == "days") { if (facet_on) "Date" else "DAH" } else { "week" } x_label <- switch(x_unit, "days" = if (facet_on) "Date" else "Age of Crop (Days)", "weeks" = "Week Number") # Choose color palette based on colorblind_friendly flag color_scale <- if (colorblind_friendly) { ggplot2::scale_color_brewer(type = "qual", palette = "Set2") } else { ggplot2::scale_color_discrete() } # Calculate dynamic max values for breaks max_dah_both <- max(plot_data_both$DAH, na.rm = TRUE) + 20 max_week_both <- max(as.numeric(plot_data_both$week), na.rm = TRUE) + ceiling(20 / 7) # Create the faceted plot g_both <- ggplot2::ggplot(data = plot_data_both) + # Add benchmark lines first (behind season lines) { if (!is.null(benchmark_data)) { benchmark_subset <- benchmark_data %>% dplyr::mutate( benchmark_x = if (x_var == "DAH") { DAH } else if (x_var == "week") { DAH / 7 } else { DAH }, ci_type_label = case_when( ci_type == "value" ~ "10-Day Rolling Mean CI", ci_type == "cumulative_CI" ~ "Cumulative CI", TRUE ~ ci_type ) ) ggplot2::geom_smooth( data = benchmark_subset, ggplot2::aes( x = .data[["benchmark_x"]], y = .data[["benchmark_value"]], group = factor(.data[["percentile"]]) ), color = "gray70", linewidth = 0.5, se = FALSE, inherit.aes = FALSE ) } } + ggplot2::facet_wrap(~ci_type_label, scales = "free_y") + # Plot older seasons with lighter lines ggplot2::geom_line( data = plot_data_both %>% dplyr::filter(!is_latest), ggplot2::aes( x = .data[[x_var]], y = .data[["ci_value"]], col = .data[["season"]], group = .data[["season"]] ), linewidth = 0.7, alpha = 0.4 ) + # Plot latest season with thicker, more prominent line ggplot2::geom_line( data = plot_data_both %>% dplyr::filter(is_latest), ggplot2::aes( x = .data[[x_var]], y = .data[["ci_value"]], col = .data[["season"]], group = .data[["season"]] ), linewidth = 1.5, alpha = 1 ) + ggplot2::labs(title = paste("CI Analysis for Field", pivotName), color = "Season", y = "CI Value", x = x_label) + color_scale + { if (x_var == "DAH") { ggplot2::scale_x_continuous(breaks = seq(0, 450, by = 50), sec.axis = ggplot2::sec_axis(~ . / 30.44, name = "Age in Months", breaks = seq(0, 14, by = 1))) } else if (x_var == "week") { ggplot2::scale_x_continuous(breaks = seq(0, max_week_both, by = 4), sec.axis = ggplot2::sec_axis(~ . / 4.348, name = "Age in Months", breaks = seq(0, 14, by = 1))) } else if (x_var == "Date") { ggplot2::scale_x_date(breaks = "1 month", date_labels = "%b-%Y", sec.axis = ggplot2::sec_axis(~ ., name = "Age in Months", breaks = scales::breaks_pretty())) } } + ggplot2::theme_minimal() + ggplot2::theme(axis.text.x = ggplot2::element_text(hjust = 0.5), axis.text.x.top = ggplot2::element_text(hjust = 0.5), axis.title.x.top = ggplot2::element_text(size = 8), legend.justification = c(1, 0), legend.position = "inside", legend.position.inside = c(1, 0), legend.title = ggplot2::element_text(size = 8), legend.text = ggplot2::element_text(size = 8)) + ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = TRUE)) # For the rolling mean data, we want to set reasonable y-axis limits # Since we're using free_y scales, each facet will have its own y-axis # The rolling mean will automatically scale to its data range, # but we can ensure it shows the 0-7 context by adding invisible points # Add invisible points to set the y-axis range for rolling mean facet dummy_data <- data.frame( ci_type_label = "10-Day Rolling Mean CI", ci_value = c(0, 7), stringsAsFactors = FALSE ) dummy_data[[x_var]] <- range(plot_data_both[[x_var]], na.rm = TRUE) dummy_data[["season"]] <- factor("dummy", levels = levels(plot_data_both[["season"]])) g_both <- g_both + ggplot2::geom_point( data = dummy_data, ggplot2::aes(x = .data[[x_var]], y = .data[["ci_value"]]), alpha = 0, size = 0 ) # Invisible points to set scale # Display the combined faceted plot subchunkify(g_both, 2.8, 10) } }, error = function(e) { safe_log(paste("Error creating CI trend plot for pivot", pivotName, ":", e$message), "ERROR") cum_ci_plot2(pivotName) # Use fallback function in case of error }) } #' Fallback function for creating CI visualization when data is missing #' #' @param pivotName The name or ID of the pivot field to visualize #' @return NULL (adds output directly to R Markdown document) #' cum_ci_plot2 <- function(pivotName){ # Input validation if (missing(pivotName) || is.null(pivotName) || pivotName == "") { stop("pivotName is required") } # Create a simple plot showing "No data available" tryCatch({ end_date <- Sys.Date() start_date <- end_date %m-% months(11) # 11 months ago from end_date date_seq <- seq.Date(from = start_date, to = end_date, by = "month") midpoint_date <- start_date + (end_date - start_date) / 2 g <- ggplot() + scale_x_date(limits = c(start_date, end_date), date_breaks = "1 month", date_labels = "%m-%Y") + scale_y_continuous(limits = c(0, 4)) + labs(title = paste("14 day rolling MEAN CI rate - Field ", pivotName), x = "Date", y = "CI Rate") + theme_minimal() + theme(axis.text.x = element_text(hjust = 0.5), legend.justification = c(1, 0), legend.position = "inside", legend.position.inside = c(1, 0), legend.title = element_text(size = 8), legend.text = element_text(size = 8)) + annotate("text", x = midpoint_date, y = 2, label = "No data available", size = 6, hjust = 0.5) subchunkify(g, 3.2, 10) }, error = function(e) { safe_log(paste("Error creating fallback CI plot for pivot", pivotName, ":", e$message), "ERROR") cat(paste("No data available for field", pivotName, "\n")) }) } #' Gets the file path for a specific week's mosaic #' #' @param mosaic_path Base directory containing mosaic files #' @param input_date Reference date to calculate from #' @param week_offset Number of weeks to offset from input date (positive or negative) #' @return File path to the requested week's mosaic TIF file #' get_week_path <- function(mosaic_path, input_date, week_offset) { # Input validation if (missing(mosaic_path) || is.null(mosaic_path) || mosaic_path == "") { stop("mosaic_path is required") } if (missing(input_date)) { stop("input_date is required") } tryCatch({ # Convert input_date to Date object (in case it's a string) input_date <- as.Date(input_date) if (is.na(input_date)) { stop("Invalid input_date. Expected a Date object or a string convertible to Date.") } # Validate week_offset week_offset <- as.integer(week_offset) if (is.na(week_offset)) { stop("Invalid week_offset. Expected an integer value.") } # Get the start of the week for the input date (adjust to Monday as the start of the week) start_of_week <- lubridate::floor_date(input_date, unit = "week", week_start = 1) # Calculate the new date after applying the week offset target_date <- start_of_week + lubridate::weeks(week_offset) # Get the week number and year of the target date target_week <- sprintf("%02d", lubridate::isoweek(target_date)) # Left-pad week number with a zero if needed target_year <- lubridate::isoyear(target_date) # Load single-file mosaic for the given week path_to_week <- here::here(mosaic_path, paste0("week_", target_week, "_", target_year, ".tif")) # Log the path calculation safe_log(paste("Calculated path for week", target_week, "of year", target_year, ":", path_to_week), "INFO") # Return the path return(path_to_week) }, error = function(e) { safe_log(paste("Error calculating week path:", e$message), "ERROR") stop(e$message) }) } #' Computes historical percentile benchmarks for CI data per estate #' #' @param ci_quadrant_data Data frame containing CI quadrant data with field, Date, DAH, cumulative_CI, value, season columns #' @param estate_name Name of the estate/client to filter data for #' @param percentiles Vector of percentiles to compute (e.g., c(10, 50, 90)) #' @param min_seasons Minimum number of seasons required for reliable benchmarks (default: 3) #' @return Data frame with DAH, percentile, ci_type, benchmark_value, or NULL if insufficient data #' compute_ci_benchmarks <- function(ci_quadrant_data, estate_name, percentiles = c(10, 50, 90), min_seasons = 3) { # Input validation if (missing(ci_quadrant_data) || is.null(ci_quadrant_data)) { stop("ci_quadrant_data is required") } if (missing(estate_name) || is.null(estate_name) || estate_name == "") { stop("estate_name is required") } if (!all(percentiles >= 0 & percentiles <= 100)) { stop("percentiles must be between 0 and 100") } tryCatch({ # Filter data for the specified estate (assuming estate is not directly in data, but we can infer from context) # Since the data is per field, and fields are unique to estates, we'll use all data but could add estate filtering if available data_filtered <- ci_quadrant_data # Check if we have enough seasons unique_seasons <- unique(data_filtered$season) if (length(unique_seasons) < min_seasons) { safe_log(paste("Insufficient historical seasons for estate", estate_name, ":", length(unique_seasons), "seasons found, need at least", min_seasons), "WARNING") return(NULL) } # Prepare data for both CI types data_prepared <- data_filtered %>% dplyr::ungroup() %>% # Ensure no existing groupings dplyr::select(DAH, value, cumulative_CI, season) %>% tidyr::pivot_longer( cols = c("value", "cumulative_CI"), names_to = "ci_type", values_to = "ci_value" ) %>% dplyr::filter(!is.na(ci_value)) # Remove NA values # Compute percentiles for each DAH and ci_type benchmarks <- data_prepared %>% dplyr::group_by(DAH, ci_type) %>% dplyr::summarise( p10 = tryCatch(quantile(ci_value, 0.1, na.rm = TRUE), error = function(e) NA_real_), p50 = tryCatch(quantile(ci_value, 0.5, na.rm = TRUE), error = function(e) NA_real_), p90 = tryCatch(quantile(ci_value, 0.9, na.rm = TRUE), error = function(e) NA_real_), n_observations = n(), .groups = 'drop' ) %>% dplyr::filter(n_observations >= min_seasons) %>% # Only include DAHs with sufficient data tidyr::pivot_longer( cols = c(p10, p50, p90), names_to = "percentile", values_to = "benchmark_value" ) %>% dplyr::mutate( percentile = case_when( percentile == "p10" ~ 10, percentile == "p50" ~ 50, percentile == "p90" ~ 90 ) ) %>% dplyr::filter(!is.na(benchmark_value)) # Remove any NA benchmarks # Rename columns for clarity benchmarks <- benchmarks %>% dplyr::select(DAH, ci_type, percentile, benchmark_value) safe_log(paste("Computed CI benchmarks for estate", estate_name, "with", length(unique_seasons), "seasons and", nrow(benchmarks), "benchmark points"), "INFO") return(benchmarks) }, error = function(e) { safe_log(paste("Error computing CI benchmarks for estate", estate_name, ":", e$message), "ERROR") print(paste("DEBUG: Error details:", e$message, "Call:", deparse(e$call))) return(NULL) }) } #' Aggregate per-field weekly mosaics into a farm-level mosaic #' #' Reads all per-field mosaic TIFs for a given week and merges them into a single farm-level mosaic. #' This is used for overview maps in the report (Script 90). #' #' Per-field mosaics already have proper geospatial metadata (CRS, geotransform) from Script 10, #' so terra can align them automatically without needing field boundaries or extent information. #' #' @param weekly_mosaic_dir Path to weekly_mosaic directory (e.g., "laravel_app/storage/app/{project}/weekly_mosaic") #' @param target_week ISO week number (e.g., 52) #' @param target_year ISO year (e.g., 2025) #' @return SpatRaster object (5-band: R,G,B,NIR,CI) or NULL if no fields found #' #' @details #' Per-field mosaics are located at: weekly_mosaic/{FIELD}/week_WW_YYYY.tif #' This function: #' 1. Finds all per-field subdirectories #' 2. Loads each field's weekly mosaic #' 3. Merges to a single raster using terra::mosaic() (alignment handled automatically by metadata) #' 4. Returns combined 5-band raster for visualization #' aggregate_per_field_mosaics_to_farm_level <- function( weekly_mosaic_dir, target_week, target_year ) { tryCatch({ # Validate directory exists if (!dir.exists(weekly_mosaic_dir)) { safe_log(paste("Weekly mosaic directory not found:", weekly_mosaic_dir), "WARNING") return(NULL) } # Find all per-field subdirectories (non-TIF files at top level) all_items <- list.files(weekly_mosaic_dir, full.names = FALSE) field_dirs <- all_items[ !grepl("\\.tif$", all_items, ignore.case = TRUE) & dir.exists(file.path(weekly_mosaic_dir, all_items)) ] if (length(field_dirs) == 0) { safe_log(paste("No per-field directories found in", weekly_mosaic_dir), "WARNING") return(NULL) } safe_log(paste("Found", length(field_dirs), "field directories. Aggregating week", sprintf("%02d", target_week), target_year), "INFO") # Collect rasters from each field raster_list <- list() for (field_dir in field_dirs) { field_mosaic_path <- file.path( weekly_mosaic_dir, field_dir, paste0("week_", sprintf("%02d", target_week), "_", target_year, ".tif") ) if (file.exists(field_mosaic_path)) { tryCatch({ r <- terra::rast(field_mosaic_path) raster_list[[field_dir]] <- r safe_log(paste("Loaded mosaic for field:", field_dir), "DEBUG") }, error = function(e) { safe_log(paste("Could not load mosaic for field", field_dir, ":", e$message), "WARNING") }) } } if (length(raster_list) == 0) { safe_log(paste("No field mosaics found for week", sprintf("%02d", target_week), target_year), "WARNING") return(NULL) } safe_log(paste("Successfully loaded mosaics for", length(raster_list), "fields"), "INFO") # 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") # 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("merge() returned NULL") } return(farm_mosaic) }, error = function(e) { safe_log(paste("Error during mosaic creation:", e$message), "ERROR") return(NULL) }) }, error = function(e) { safe_log(paste("Error aggregating per-field mosaics:", e$message), "ERROR") return(NULL) }) } #' Get per-field mosaic path (new per-field architecture) #' #' Returns the path to a per-field weekly mosaic for direct visualization. #' Replaces the old cropping workflow: now we load the field's own mosaic instead of cropping farm-level. #' #' @param weekly_mosaic_dir Path to weekly_mosaic directory #' @param field_name Name of the field #' @param target_week ISO week number #' @param target_year ISO year #' @return Path to field-specific mosaic TIF, or NULL if not found #' get_per_field_mosaic_path <- function( weekly_mosaic_dir, field_name, target_week, target_year ) { path <- file.path( weekly_mosaic_dir, field_name, paste0("week_", sprintf("%02d", target_week), "_", target_year, ".tif") ) if (file.exists(path)) { return(path) } else { safe_log(paste("Per-field mosaic not found for field", field_name, "week", sprintf("%02d", target_week), target_year), "WARNING") return(NULL) } } #' Determine field priority level based on CV and Moran's I #' #' @param cv Coefficient of Variation (uniformity metric) #' @param morans_i Moran's I spatial autocorrelation index #' @return Priority level: 1=Urgent, 2=Monitor, 3=No stress #' get_field_priority_level <- function(cv, morans_i) { # Handle NA values if (is.na(cv) || is.na(morans_i)) return(3) # Default to no stress # Determine priority based on thresholds if (cv < 0.1) { if (morans_i < 0.7) { return(3) # No stress } else if (morans_i <= 0.9) { return(2) # Monitor (young field with some clustering) } else { return(1) # Urgent } } else if (cv <= 0.15) { if (morans_i < 0.7) { return(2) # Monitor } else { return(1) # Urgent } } else { # cv > 0.15 return(1) # Urgent } } #' Generate field-specific KPI summary for display in reports #' #' @param field_name Name of the field to summarize #' @param field_details_table Data frame with field-level KPI details #' @param CI_quadrant Data frame containing CI quadrant data with Date, DAH, season columns #' @param report_date Report date (used for filtering current season data) #' @return Formatted text string with field KPI summary #' generate_field_kpi_summary <- function(field_name, field_details_table, CI_quadrant, report_date = Sys.Date()) { tryCatch({ # Get field age from CI quadrant data for the CURRENT SEASON only # First identify the current season for this field # Get field age from CI quadrant data for the CURRENT SEASON only # First identify the current season for this field current_season_data <- CI_quadrant %>% filter(field == field_name, Date <= as.Date(report_date)) %>% group_by(season) %>% summarise(season_end = max(Date), .groups = 'drop') %>% filter(season == max(season)) if (nrow(current_season_data) == 0) { return(paste("**Field", field_name, "KPIs:** No CI data available for current season")) } current_season <- current_season_data %>% pull(season) # Get the most recent DAH from the current season field_age_data <- CI_quadrant %>% filter(field == field_name, season == current_season) %>% pull(DAH) field_age <- if (length(field_age_data) > 0) max(field_age_data, na.rm = TRUE) else NA_real_ # Filter data for this specific field field_data <- field_details_table %>% filter(Field == field_name) if (nrow(field_data) == 0) { return(paste("**Field", field_name, "KPIs:** Data not available")) } # Aggregate sub-field data for field-level summary # For categorical data, take the most common value or highest risk level field_summary <- field_data %>% summarise( field_size = sum(`Field Size (acres)`, na.rm = TRUE), uniformity_levels = paste(unique(`Growth Uniformity`), collapse = "/"), avg_yield_forecast = ifelse(is.na(`Yield Forecast (t/ha)`[1]), NA, mean(`Yield Forecast (t/ha)`, na.rm = TRUE)), max_gap_score = max(`Gap Score`, na.rm = TRUE), highest_decline_risk = case_when( any(`Decline Risk` == "Very-high") ~ "Very-high", any(`Decline Risk` == "High") ~ "High", any(`Decline Risk` == "Moderate") ~ "Moderate", any(`Decline Risk` == "Low") ~ "Low", TRUE ~ "Unknown" ), highest_weed_risk = case_when( any(`Weed Risk` == "High") ~ "High", any(`Weed Risk` == "Moderate") ~ "Moderate", any(`Weed Risk` == "Low") ~ "Low", TRUE ~ "Unknown" ), avg_mean_ci = mean(`Mean CI`, na.rm = TRUE), avg_cv = mean(`CV Value`, na.rm = TRUE), .groups = 'drop' ) # Apply age-based filtering to yield forecast if (is.na(field_age) || field_age < 240) { field_summary$avg_yield_forecast <- NA_real_ } # Format the summary text yield_text <- if (is.na(field_summary$avg_yield_forecast)) { "Yield Forecast: NA" } else { paste0("Yield Forecast: ", round(field_summary$avg_yield_forecast, 1), " t/ha") } kpi_text <- paste0( "Size: ", round(field_summary$field_size * 0.404686, 1), " ha | Mean CI: ", round(field_summary$avg_mean_ci, 2), " | Growth Uniformity: ", field_summary$uniformity_levels, " | ", yield_text, " | Gap Score: ", round(field_summary$max_gap_score, 1), " | Decline Risk: ", field_summary$highest_decline_risk, " | Weed Risk: ", field_summary$highest_weed_risk ) # Wrap in smaller text HTML tags for Word output kpi_text <- paste0("", kpi_text, "") return(kpi_text) }, error = function(e) { safe_log(paste("Error generating KPI summary for field", field_name, ":", e$message), "ERROR") return(paste("**Field", field_name, "KPIs:** Error generating summary")) }) } #' Normalize field_details_table column structure #' #' Standardizes column names and ensures all expected KPI columns exist. #' Handles Field → Field_id rename and injects missing columns as NA. #' #' @param field_details_table data.frame to normalize #' @return data.frame with standardized column structure normalize_field_details_columns <- function(field_details_table) { if (is.null(field_details_table) || nrow(field_details_table) == 0) { return(field_details_table) } # Rename Field → Field_id if needed if ("Field" %in% names(field_details_table) && !("Field_id" %in% names(field_details_table))) { field_details_table <- field_details_table %>% dplyr::rename(Field_id = Field) } # Ensure all expected KPI columns exist; add as NA if missing expected_cols <- c("Field_id", "Mean_CI", "CV", "TCH_Forecasted", "Gap_Score", "Trend_Interpretation", "Weekly_CI_Change", "Uniformity_Interpretation", "Decline_Severity", "Patchiness_Risk") for (col in expected_cols) { if (!col %in% names(field_details_table)) { field_details_table[[col]] <- NA } } return(field_details_table) }