From bb2a599075b48e6c28ec1ea8a2b1ea3a9d3132f4 Mon Sep 17 00:00:00 2001 From: Timon Date: Wed, 23 Apr 2025 09:47:19 +0200 Subject: [PATCH] Enhanced SmartCane executive summary report with explanatory text and fixed priority map coloring Added explanatory text for all visualizations Fixed priority map color scheme (red=high priority, green=low priority) Improved error handling in farm health data calculations Added fallback mechanisms for missing data --- r_app/CI_report_dashboard_planet.Rmd | 2 +- r_app/CI_report_dashboard_planet_enhanced.Rmd | 1145 +++++++++++++ r_app/CI_report_executive_summary.Rmd | 1463 +++++++++++++++++ r_app/ci_extraction.R | 15 +- r_app/interpolate_growth_model.R | 14 +- r_app/mosaic_creation.R | 14 +- r_app/packages.R | 9 +- r_app/tests/test_report_utils.R | 280 ++++ 8 files changed, 2934 insertions(+), 8 deletions(-) create mode 100644 r_app/CI_report_dashboard_planet_enhanced.Rmd create mode 100644 r_app/CI_report_executive_summary.Rmd create mode 100644 r_app/tests/test_report_utils.R diff --git a/r_app/CI_report_dashboard_planet.Rmd b/r_app/CI_report_dashboard_planet.Rmd index b064334..71de512 100644 --- a/r_app/CI_report_dashboard_planet.Rmd +++ b/r_app/CI_report_dashboard_planet.Rmd @@ -437,7 +437,7 @@ tryCatch({ dplyr::ungroup() # Check if tonnage_ha is empty - if (all(is.na(CI_quadrant$tonnage_ha))) { + if (all(is.na(harvesting_data$tonnage_ha))) { safe_log("Lacking historic harvest data, please provide for yield prediction calculation", "WARNING") knitr::knit_exit() # Exit the chunk if tonnage_ha is empty } diff --git a/r_app/CI_report_dashboard_planet_enhanced.Rmd b/r_app/CI_report_dashboard_planet_enhanced.Rmd new file mode 100644 index 0000000..489b4be --- /dev/null +++ b/r_app/CI_report_dashboard_planet_enhanced.Rmd @@ -0,0 +1,1145 @@ +--- +params: + ref: "word-styles-reference-var1.docx" + output_file: CI_report.docx + report_date: "2024-08-28" + data_dir: "Chemba" + mail_day: "Wednesday" + borders: TRUE + use_breaks: FALSE +output: + # html_document: + # toc: yes + # df_print: paged + word_document: + reference_docx: !expr file.path("word-styles-reference-var1.docx") + toc: yes +editor_options: + chunk_output_type: console +--- + +```{r setup_parameters, include=FALSE} +# Set up basic report parameters from input values +report_date <- params$report_date +mail_day <- params$mail_day +borders <- params$borders +use_breaks <- params$use_breaks # Whether to use breaks or continuous spectrum in visualizations + +# Environment setup notes (commented out) +# # Activeer de renv omgeving +# renv::activate() +# renv::deactivate() +# # Optioneel: Herstel de omgeving als dat nodig is +# # Je kunt dit commentaar geven als je het normaal niet wilt uitvoeren +# renv::restore() +``` + +```{r load_libraries, message=FALSE, warning=FALSE, include=FALSE} +# Configure knitr options +knitr::opts_chunk$set(warning = FALSE, message = FALSE) + +# Path management +library(here) + +# Spatial data libraries +library(sf) +library(terra) +library(exactextractr) +# library(raster) - Removed as it's no longer maintained + +# Data manipulation and visualization +library(tidyverse) # Includes dplyr, ggplot2, etc. +library(tmap) +library(lubridate) +library(zoo) + +# Machine learning +library(rsample) +library(caret) +library(randomForest) +library(CAST) + +# Parallel processing +library(future) +library(furrr) +library(progressr) + +# Load custom utility functions +tryCatch({ + source("report_utils.R") +}, error = function(e) { + message(paste("Error loading report_utils.R:", e$message)) + # Try alternative path if the first one fails + tryCatch({ + source(here::here("r_app", "report_utils.R")) + }, error = function(e) { + stop("Could not load report_utils.R from either location: ", e$message) + }) +}) +``` + +```{r setup_parallel_processing, message=FALSE, warning=FALSE, include=FALSE} +# Set up parallel processing to speed up report generation +# Use half of available cores to avoid overloading the system +num_workers <- parallel::detectCores() / 2 +num_workers <- floor(max(1, num_workers)) # At least 1, but no more than half of cores + +# Set up future plan for parallel processing +future::plan(future::multisession, workers = num_workers) +safe_log(paste("Set up parallel processing with", num_workers, "workers")) + +# Configure progressr reporting +progressr::handlers(progressr::handler_progress( + format = "[:bar] :percent :eta :message", + width = 60 +)) +``` + +```{r initialize_project_config, message=FALSE, warning=FALSE, include=FALSE} +# Set the project directory from parameters +project_dir <- params$data_dir + +# Source project parameters with error handling +tryCatch({ + source(here::here("r_app", "parameters_project.R")) +}, error = function(e) { + stop("Error loading parameters_project.R: ", e$message) +}) + +# Log initial configuration +safe_log("Starting the R Markdown script") +safe_log(paste("mail_day params:", params$mail_day)) +safe_log(paste("report_date params:", params$report_date)) +safe_log(paste("mail_day variable:", mail_day)) +``` + +```{r calculate_dates_and_weeks, message=FALSE, warning=FALSE, include=FALSE} +# Set locale for consistent date formatting +Sys.setlocale("LC_TIME", "C") + +# Initialize date variables from parameters +today <- as.character(report_date) +mail_day_as_character <- as.character(mail_day) + +# Calculate week days +report_date_as_week_day <- weekdays(lubridate::ymd(today)) +days_of_week <- c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday") + +# Calculate initial week number +week <- lubridate::week(today) +safe_log(paste("Initial week calculation:", week, "today:", today)) + +# Calculate previous dates for comparisons +today_minus_1 <- as.character(lubridate::ymd(today) - 7) +today_minus_2 <- as.character(lubridate::ymd(today) - 14) +today_minus_3 <- as.character(lubridate::ymd(today) - 21) + +# Log the weekday calculations for debugging +safe_log(paste("Report date weekday:", report_date_as_week_day)) +safe_log(paste("Weekday index:", which(days_of_week == report_date_as_week_day))) +safe_log(paste("Mail day:", mail_day_as_character)) +safe_log(paste("Mail day index:", which(days_of_week == mail_day_as_character))) + +# Adjust week calculation based on mail day +if (which(days_of_week == report_date_as_week_day) > which(days_of_week == mail_day_as_character)) { + safe_log("Adjusting weeks because of mail day") + week <- lubridate::week(today) + 1 + today_minus_1 <- as.character(lubridate::ymd(today)) + today_minus_2 <- as.character(lubridate::ymd(today) - 7) + today_minus_3 <- as.character(lubridate::ymd(today) - 14) +} + +# Generate subtitle for report +subtitle_var <- paste("Report generated on", Sys.Date()) + +# Calculate week numbers for previous weeks +week_minus_1 <- week - 1 +week_minus_2 <- week - 2 +week_minus_3 <- week - 3 + +# Format current week with leading zeros +week <- sprintf("%02d", week) + +# Get years for each date +year <- lubridate::year(today) +year_1 <- lubridate::year(today_minus_1) +year_2 <- lubridate::year(today_minus_2) +year_3 <- lubridate::year(today_minus_3) +``` + +```{r data, message=TRUE, warning=TRUE, include=FALSE} +# Load CI index data with error handling +tryCatch({ + CI_quadrant <- readRDS(here::here(cumulative_CI_vals_dir, "All_pivots_Cumulative_CI_quadrant_year_v2.rds")) + safe_log("Successfully loaded CI quadrant data") +}, error = function(e) { + stop("Error loading CI quadrant data: ", e$message) +}) + +# Get file paths for different weeks using the utility function +tryCatch({ + path_to_week_current = get_week_path(weekly_CI_mosaic, today, 0) + path_to_week_minus_1 = get_week_path(weekly_CI_mosaic, today, -1) + path_to_week_minus_2 = get_week_path(weekly_CI_mosaic, today, -2) + path_to_week_minus_3 = get_week_path(weekly_CI_mosaic, today, -3) + + # Log the calculated paths + safe_log("Required mosaic paths:") + safe_log(paste("Path to current week:", path_to_week_current)) + safe_log(paste("Path to week minus 1:", path_to_week_minus_1)) + safe_log(paste("Path to week minus 2:", path_to_week_minus_2)) + safe_log(paste("Path to week minus 3:", path_to_week_minus_3)) + + # Validate that files exist + if (!file.exists(path_to_week_current)) warning("Current week mosaic file does not exist: ", path_to_week_current) + if (!file.exists(path_to_week_minus_1)) warning("Week minus 1 mosaic file does not exist: ", path_to_week_minus_1) + if (!file.exists(path_to_week_minus_2)) warning("Week minus 2 mosaic file does not exist: ", path_to_week_minus_2) + if (!file.exists(path_to_week_minus_3)) warning("Week minus 3 mosaic file does not exist: ", path_to_week_minus_3) + + # Load raster data with terra functions + CI <- terra::rast(path_to_week_current)$CI + CI_m1 <- terra::rast(path_to_week_minus_1)$CI + CI_m2 <- terra::rast(path_to_week_minus_2)$CI + CI_m3 <- terra::rast(path_to_week_minus_3)$CI + +}, error = function(e) { + stop("Error loading raster data: ", e$message) +}) +``` + +```{r calculate_difference_rasters, message=TRUE, warning=TRUE, include=FALSE} +# Calculate difference rasters for comparisons +tryCatch({ + # Calculate weekly difference + last_week_dif_raster_abs <- (CI - CI_m1) + safe_log("Calculated weekly difference raster") + + # Calculate three-week difference + three_week_dif_raster_abs <- (CI - CI_m3) + safe_log("Calculated three-week difference raster") +}, error = function(e) { + safe_log(paste("Error calculating difference rasters:", e$message), "ERROR") + # Create placeholder rasters if calculations fail + if (!exists("last_week_dif_raster_abs")) { + last_week_dif_raster_abs <- CI * 0 + } + if (!exists("three_week_dif_raster_abs")) { + three_week_dif_raster_abs <- CI * 0 + } +}) +``` + +```{r load_field_boundaries, message=TRUE, warning=TRUE, include=FALSE} +# Load field boundaries from parameters +tryCatch({ + AllPivots0 <- field_boundaries_sf + safe_log("Successfully loaded field boundaries") +}, error = function(e) { + stop("Error loading field boundaries: ", e$message) +}) +``` + +```{r calculate_field_health_scores, message=FALSE, warning=FALSE, include=FALSE} +# Calculate health scores for all fields +tryCatch({ + # Get list of all fields + all_fields <- unique(AllPivots0$field) + + # Process field health scores + safe_log("Calculating field health scores...") + + # Use future_map instead of map for parallel processing + field_health_scores <- furrr::future_map_dfr(all_fields, function(field_name) { + tryCatch({ + # Get field data + field_shape <- AllPivots0 %>% dplyr::filter(field == field_name) + + # Get field age from harvesting data + field_age_data <- harvesting_data %>% + dplyr::filter(field == field_name) %>% + dplyr::arrange(desc(season_start)) %>% + dplyr::slice(1) + + # Set default age if not available + field_age_weeks <- if(nrow(field_age_data) > 0 && !is.na(field_age_data$age)) { + field_age_data$age + } else { + 10 # Default age if not available + } + + # Crop and mask rasters for this field + ci_current <- terra::crop(CI, field_shape) %>% terra::mask(., field_shape) + ci_change <- terra::crop(last_week_dif_raster_abs, field_shape) %>% terra::mask(., field_shape) + + # Generate health score + health_data <- generate_field_health_score(ci_current, ci_change, field_age_weeks) + + # Return as a data frame row + data.frame( + field = field_name, + health_score = health_data$score, + health_status = health_data$status, + ci_score = health_data$components$ci, + change_score = health_data$components$change, + uniformity_score = health_data$components$uniformity, + age_weeks = field_age_weeks + ) + }, error = function(e) { + safe_log(paste("Error calculating health score for field", field_name, ":", e$message), "ERROR") + # Return NA values if error occurs + data.frame( + field = field_name, + health_score = NA, + health_status = "Error", + ci_score = NA, + change_score = NA, + uniformity_score = NA, + age_weeks = NA + ) + }) + }, .options = furrr::furrr_options(seed = TRUE)) + + # Add recommendations based on health status + field_health_scores <- field_health_scores %>% + dplyr::mutate(recommendation = dplyr::case_when( + health_status == "Critical" ~ "Immediate inspection needed", + health_status == "Needs Attention" ~ "Schedule inspection this week", + health_status == "Fair" ~ "Monitor closely", + health_status == "Good" ~ "Regular monitoring", + health_status == "Excellent" ~ "Maintain current practices", + TRUE ~ "Status unknown - inspect field" + )) + + safe_log("Health scores calculation completed") +}, error = function(e) { + safe_log(paste("Error in health score calculation:", e$message), "ERROR") + # Create empty dataframe if calculation failed + field_health_scores <- data.frame( + field = character(), + health_score = numeric(), + health_status = character(), + recommendation = character() + ) +}) +``` + +```{r helper_functions, message=FALSE, warning=FALSE, include=FALSE} +#' Generate a field health score based on CI values and trends +#' +#' @param ci_current Current CI raster +#' @param ci_change CI change raster +#' @param field_age_weeks Field age in weeks +#' @return List containing score, status, and component scores +#' +generate_field_health_score <- function(ci_current, ci_change, field_age_weeks) { + # Get mean CI value for the field + mean_ci <- terra::global(ci_current, "mean", na.rm=TRUE)[[1]] + + # Get mean CI change + mean_change <- terra::global(ci_change, "mean", na.rm=TRUE)[[1]] + + # Get CI uniformity (coefficient of variation) + ci_sd <- terra::global(ci_current, "sd", na.rm=TRUE)[[1]] + ci_uniformity <- ifelse(mean_ci > 0, ci_sd / mean_ci, 1) + + # Calculate base score from current CI (scale 0-5) + # Adjusted for crop age - expectations increase with age + expected_ci <- min(5, field_age_weeks / 10) # Simple linear model + ci_score <- max(0, min(5, 5 - 2 * abs(mean_ci - expected_ci))) + + # Add points for positive change (scale 0-3) + change_score <- max(0, min(3, 1 + mean_change)) + + # Add points for uniformity (scale 0-2) + uniformity_score <- max(0, min(2, 2 * (1 - ci_uniformity))) + + # Calculate total score (0-10) + total_score <- ci_score + change_score + uniformity_score + + # Create status label + status <- dplyr::case_when( + total_score >= 8 ~ "Excellent", + total_score >= 6 ~ "Good", + total_score >= 4 ~ "Fair", + total_score >= 2 ~ "Needs Attention", + TRUE ~ "Critical" + ) + + # Return results + return(list( + score = round(total_score, 1), + status = status, + components = list( + ci = round(ci_score, 1), + change = round(change_score, 1), + uniformity = round(uniformity_score, 1) + ) + )) +} + +#' Create an irrigation recommendation map +#' +#' @param ci_current Current CI raster +#' @param ci_change CI change raster +#' @param field_shape Field boundary shape +#' @param title Map title +#' @return A tmap object with irrigation recommendations +#' +create_irrigation_map <- function(ci_current, ci_change, field_shape, title = "Irrigation Priority Zones") { + # Create a new raster for irrigation recommendations + irrigation_priority <- ci_current * 0 + + # Extract values for processing + ci_values <- terra::values(ci_current) + change_values <- terra::values(ci_change) + + # Create priority zones: + # 3 = High priority (low CI, negative trend) + # 2 = Medium priority (low CI but stable, or good CI with negative trend) + # 1 = Low priority (watch, good CI with slight decline) + # 0 = No action needed (good CI, stable/positive trend) + priority_values <- rep(NA, length(ci_values)) + + # High priority: Low CI (< 2) and negative change (< 0) + high_priority <- which(ci_values < 2 & change_values < 0 & !is.na(ci_values) & !is.na(change_values)) + priority_values[high_priority] <- 3 + + # Medium priority: Low CI (< 2) with stable/positive change, or moderate CI (2-4) with significant negative change (< -1) + medium_priority <- which( + (ci_values < 2 & change_values >= 0 & !is.na(ci_values) & !is.na(change_values)) | + (ci_values >= 2 & ci_values < 4 & change_values < -1 & !is.na(ci_values) & !is.na(change_values)) + ) + priority_values[medium_priority] <- 2 + + # Low priority (watch): Moderate/good CI (>= 2) with mild negative change (-1 to 0) + low_priority <- which( + ci_values >= 2 & change_values < 0 & change_values >= -1 & !is.na(ci_values) & !is.na(change_values) + ) + priority_values[low_priority] <- 1 + + # No action needed: Good CI (>= 2) with stable/positive change (>= 0) + no_action <- which(ci_values >= 2 & change_values >= 0 & !is.na(ci_values) & !is.na(change_values)) + priority_values[no_action] <- 0 + + # Set values in the irrigation priority raster + terra::values(irrigation_priority) <- priority_values + + # Create the map + tm_shape(irrigation_priority) + + tm_raster( + style = "cat", + palette = c("#1a9850", "#91cf60", "#fc8d59", "#d73027"), + labels = c("No Action", "Watch", "Medium Priority", "High Priority"), + title = "Irrigation Need" + ) + + tm_shape(field_shape) + + tm_borders(lwd = 2) + + tm_layout( + main.title = title, + legend.outside = FALSE, + legend.position = c("left", "bottom") + ) +} + +#' Simple mock function to get weather data for a field +#' In a real implementation, this would fetch data from a weather API +#' +#' @param start_date Start date for weather data +#' @param end_date End date for weather data +#' @param lat Latitude of the field center +#' @param lon Longitude of the field center +#' @return A data frame of weather data +#' +get_weather_data <- function(start_date, end_date, lat = -16.1, lon = 34.7) { + # This is a mock implementation - in production, you'd replace with actual API call + # to a service like OpenWeatherMap, NOAA, or other weather data provider + + # Create date sequence + dates <- seq.Date(from = as.Date(start_date), to = as.Date(end_date), by = "day") + n_days <- length(dates) + + # Generate some random but realistic weather data with seasonal patterns + # More rain in summer, less in winter (Southern hemisphere) + month_nums <- as.numeric(format(dates, "%m")) + + # Simplified seasonal patterns - adjust for your local climate + is_rainy_season <- month_nums %in% c(11, 12, 1, 2, 3, 4) + + # Generate rainfall - more in rainy season, occasional heavy rainfall + rainfall <- numeric(n_days) + rainfall[is_rainy_season] <- pmax(0, rnorm(sum(is_rainy_season), mean = 4, sd = 8)) + rainfall[!is_rainy_season] <- pmax(0, rnorm(sum(!is_rainy_season), mean = 0.5, sd = 2)) + + # Add some rare heavy rainfall events + heavy_rain_days <- sample(which(is_rainy_season), size = max(1, round(sum(is_rainy_season) * 0.1))) + rainfall[heavy_rain_days] <- rainfall[heavy_rain_days] + runif(length(heavy_rain_days), 20, 50) + + # Generate temperatures - seasonal variation + temp_mean <- 18 + 8 * sin((month_nums - 1) * pi/6) # Peak in January (month 1) + temp_max <- temp_mean + rnorm(n_days, mean = 5, sd = 1) + temp_min <- temp_mean - rnorm(n_days, mean = 5, sd = 1) + + # Create weather data frame + weather_data <- data.frame( + date = dates, + rainfall_mm = round(rainfall, 1), + temp_max_c = round(temp_max, 1), + temp_min_c = round(temp_min, 1), + temp_mean_c = round((temp_max + temp_min) / 2, 1) + ) + + return(weather_data) +} + +#' Creates a weather summary visualization integrated with CI data +#' +#' @param pivotName Name of the pivot field +#' @param ci_data CI quadrant data +#' @param days_to_show Number of days of weather to show +#' @return ggplot object +#' +create_weather_ci_plot <- function(pivotName, ci_data = CI_quadrant, days_to_show = 30) { + # Get field data + field_data <- ci_data %>% + dplyr::filter(field == pivotName) %>% + dplyr::arrange(Date) %>% + dplyr::filter(!is.na(value)) + + if (nrow(field_data) == 0) { + return(ggplot() + + annotate("text", x = 0, y = 0, label = "No data available") + + theme_void()) + } + + # Get the latest date and 30 days before + latest_date <- max(field_data$Date, na.rm = TRUE) + start_date <- latest_date - days_to_show + + # Filter for recent data only + recent_field_data <- field_data %>% + dplyr::filter(Date >= start_date) + + # Get center point coordinates for the field (would be calculated from geometry in production) + # This is mocked for simplicity + lat <- -16.1 # Mock latitude + lon <- 34.7 # Mock longitude + + # Get weather data + weather_data <- get_weather_data(start_date, latest_date, lat, lon) + + # Aggregate CI data to daily mean across subfields if needed + daily_ci <- recent_field_data %>% + dplyr::group_by(Date) %>% + dplyr::summarize(mean_ci = mean(value, na.rm = TRUE)) + + # Create combined plot with dual y-axis + g <- ggplot() + + # Rainfall as bars + geom_col(data = weather_data, aes(x = date, y = rainfall_mm), + fill = "#1565C0", alpha = 0.7, width = 0.7) + + # CI as a line + geom_line(data = daily_ci, aes(x = Date, y = mean_ci * 10), + color = "#2E7D32", size = 1) + + geom_point(data = daily_ci, aes(x = Date, y = mean_ci * 10), + color = "#2E7D32", size = 2) + + # Temperature range as ribbon + geom_ribbon(data = weather_data, + aes(x = date, ymin = temp_min_c, ymax = temp_max_c), + fill = "#FF9800", alpha = 0.2) + + # Primary y-axis (rainfall) + scale_y_continuous( + name = "Rainfall (mm)", + sec.axis = sec_axis(~./10, name = "Chlorophyll Index & Temperature (°C)") + ) + + labs( + title = paste("Field", pivotName, "- Weather and CI Relationship"), + subtitle = paste("Last", days_to_show, "days"), + x = "Date" + ) + + theme_minimal() + + theme( + axis.title.y.left = element_text(color = "#1565C0"), + axis.title.y.right = element_text(color = "#2E7D32"), + legend.position = "bottom" + ) + + return(g) +} +``` + +`r subtitle_var` + +\pagebreak +# Explanation of the Report + +This report provides a detailed analysis of your sugarcane fields based on satellite imagery, helping you monitor crop health and development throughout the growing season. The data is processed weekly to give you timely insights for optimal farm management decisions. + +## What is the Chlorophyll Index (CI)? + +The **Chlorophyll Index (CI)** is a vegetation index that measures the relative amount of chlorophyll in plant leaves. Chlorophyll is the green pigment responsible for photosynthesis in plants. Higher CI values indicate: + +* Greater photosynthetic activity +* Healthier plant tissue +* Better nitrogen uptake +* More vigorous crop growth + +CI values typically range from 0 (bare soil or severely stressed vegetation) to 7+ (very healthy, dense vegetation). For sugarcane, values between 3-7 generally indicate good crop health, depending on the growth stage. + +## What You'll Find in This Report: + +1. **Chlorophyll Index Overview Map**: A comprehensive view of all your fields showing current CI values. This helps identify which fields are performing well and which might need attention. + +2. **Weekly Difference Map**: Shows changes in CI values over the past week. Positive values (green) indicate improving crop health, while negative values (red) may signal stress or decline. + +3. **Field-by-Field Analysis**: Detailed maps for each field showing: + * CI values for the current week and two previous weeks + * Week-to-week changes in CI values + * Three-week change in CI values to track longer-term trends + +4. **Growth Trend Graphs**: Time-series visualizations showing how CI values have changed throughout the growing season for each section of your fields. + +5. **Yield Prediction**: For mature crops (over 300 days), we provide estimated yield predictions based on historical data and current CI measurements. + +Use these insights to identify areas that may need irrigation, fertilization, or other interventions, and to track the effectiveness of your management practices over time. + +\pagebreak +# Chlorophyll Index (CI) Overview Map - Current Week +```{r render_ci_overview_map, echo=FALSE, fig.height=6.8, fig.width=9, message=FALSE, warning=FALSE} +# Create overview chlorophyll index map +tryCatch({ + # Base shape + map <- tmap::tm_shape(CI, unit = "m") + + # Add raster layer with either breaks or continuous spectrum based on parameter + if (use_breaks) { + map <- map + tmap::tm_raster(breaks = c(0,0.5,1,2,3,4,5,6,7,Inf), + palette = "RdYlGn", + midpoint = NA, + legend.is.portrait = FALSE, + title = "Chlorophyll Index (CI)") + } else { + map <- map + tmap::tm_raster(palette = "RdYlGn", + style = "cont", + midpoint = NA, + legend.is.portrait = FALSE, + title = "Chlorophyll Index (CI)") + } + + # Complete the map with layout and other elements + map <- map + tmap::tm_layout(legend.outside = TRUE, + legend.outside.position = "bottom", + legend.show = TRUE) + + tmap::tm_scale_bar(position = tm_pos_out("right", "bottom"), text.color = "black") + + tmap::tm_compass(position = tm_pos_out("right", "bottom"), text.color = "black") + + tmap::tm_shape(AllPivots0) + + tmap::tm_borders(col = "black") + + tmap::tm_text("sub_field", size = 0.6, col = "black") + + # Print the map + print(map) +}, error = function(e) { + safe_log(paste("Error creating 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) +}) +``` +\newpage + +# Weekly Chlorophyll Index Difference Map +```{r render_ci_difference_map, echo=FALSE, fig.height=6.8, fig.width=9, message=FALSE, warning=FALSE} +# Create chlorophyll index difference map +tryCatch({ + # Base shape + map <- tmap::tm_shape(last_week_dif_raster_abs, unit = "m") + + # Add raster layer with either breaks or continuous spectrum based on parameter + if (use_breaks) { + map <- map + tmap::tm_raster(breaks = c(-3,-2,-1,0,1,2,3), + palette = "RdYlGn", + midpoint = 0, + legend.is.portrait = FALSE, + title = "Chlorophyll Index (CI) Change") + } else { + map <- map + tmap::tm_raster(palette = "RdYlGn", + style = "cont", + midpoint = 0, + legend.is.portrait = FALSE, + title = "Chlorophyll Index (CI) Change") + } + + # Complete the map with layout and other elements + map <- map + tmap::tm_layout(legend.outside = TRUE, + legend.outside.position = "bottom", + legend.show = TRUE) + + tmap::tm_scale_bar(position = tm_pos_out("right", "bottom"), text.color = "black") + + tmap::tm_compass(position = tm_pos_out("right", "bottom"), text.color = "black") + + tmap::tm_shape(AllPivots0) + + tmap::tm_borders(col = "black") + + tmap::tm_text("sub_field", size = 0.6, col = "black") + + # Print the map + print(map) +}, error = function(e) { + safe_log(paste("Error creating 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) +}) +``` +\newpage + +\pagebreak +# Field Health Overview + +The Field Health Scorecard provides an at-a-glance view of all your fields' current health status. Each field is scored on a scale of 0-10 based on: + +- **Current CI Value** (0-5 points): How well the field's chlorophyll levels match expectations for its growth stage +- **Recent CI Change** (0-3 points): Whether the field is improving or declining over the last week +- **Field Uniformity** (0-2 points): How consistent the chlorophyll levels are across the field + +This helps you quickly identify which fields need attention and which are performing well. + +```{r render_health_scorecard, echo=FALSE, fig.height=6, fig.width=10, message=FALSE, warning=FALSE} +# Create field health scorecard visualization +tryCatch({ + # Sort fields by health score + sorted_health_scores <- field_health_scores %>% + dplyr::arrange(desc(health_score)) + + # Create color mapping for status categories + status_colors <- c( + "Excellent" = "#1a9850", + "Good" = "#91cf60", + "Fair" = "#fee08b", + "Needs Attention" = "#fc8d59", + "Critical" = "#d73027", + "Error" = "#999999" + ) + + # Create the bar chart + g <- ggplot2::ggplot(sorted_health_scores, + ggplot2::aes(x = reorder(field, health_score), + y = health_score, + fill = health_status)) + + ggplot2::geom_bar(stat = "identity") + + ggplot2::geom_text(ggplot2::aes(label = health_score), + hjust = -0.2, + size = 3) + + ggplot2::coord_flip() + + ggplot2::scale_fill_manual(values = status_colors) + + ggplot2::scale_y_continuous(limits = c(0, 11)) + # Add space for labels + ggplot2::labs(title = "Field Health Scores", + x = "", + y = "Health Score (0-10)", + fill = "Status") + + ggplot2::theme_minimal() + + ggplot2::theme( + plot.title = ggplot2::element_text(face = "bold", size = 14), + axis.text.y = ggplot2::element_text(size = 10), + legend.position = "bottom" + ) + + # Print the chart + print(g) + + # Create and print the table with recommendations + health_table <- sorted_health_scores %>% + dplyr::select(field, health_score, health_status, recommendation, age_weeks) %>% + dplyr::rename( + "Field" = field, + "Score" = health_score, + "Status" = health_status, + "Recommendation" = recommendation, + "Age (Weeks)" = age_weeks + ) + + knitr::kable(health_table, + caption = "Field Health Status and Recommendations", + digits = 1) + +}, error = function(e) { + safe_log(paste("Error creating health scorecard:", e$message), "ERROR") + plot(1, type="n", axes=FALSE, xlab="", ylab="") + text(1, 1, "Error creating health scorecard visualization", cex=1.5) +}) +``` + +\pagebreak +# Irrigation Priority Map + +This map highlights areas that may need irrigation based on current CI values and recent changes. The irrigation priority is determined by combining current crop health with recent trends: + +- **High Priority (Red)**: Low CI values with declining trends - these areas need immediate attention +- **Medium Priority (Orange)**: Either low CI with stable/improving trends or moderate CI with significant decline +- **Watch (Yellow)**: Areas with acceptable CI but showing slight negative trends +- **No Action Needed (Green)**: Areas with good CI values and stable or improving trends + +```{r render_irrigation_map, echo=FALSE, fig.height=6, fig.width=9, message=FALSE, warning=FALSE} +# Create overall irrigation priority map +tryCatch({ + # Create the map + irrigation_map <- create_irrigation_map( + ci_current = CI, + ci_change = last_week_dif_raster_abs, + field_shape = AllPivots0, + title = "Farm-Wide Irrigation Priority Zones" + ) + + # Add field labels and borders + irrigation_map <- irrigation_map + + tm_shape(AllPivots0) + + tm_borders(col = "black") + + tm_text("field", size = 0.6) + + tm_layout(legend.outside = TRUE, + legend.outside.position = "bottom") + + tm_scale_bar(position = tm_pos_out("right", "bottom")) + + # Print the map + print(irrigation_map) +}, error = function(e) { + safe_log(paste("Error creating irrigation priority map:", e$message), "ERROR") + plot(1, type="n", axes=FALSE, xlab="", ylab="") + text(1, 1, "Error creating irrigation priority map", cex=1.5) +}) +``` + +\pagebreak +# Weather and CI Relationship + +This section shows the relationship between recent weather patterns and crop health. Understanding this relationship can help identify whether changes in CI are due to weather factors or other issues that may require management intervention. + +```{r render_weather_integration, echo=FALSE, fig.height=7, fig.width=10, message=FALSE, warning=FALSE} +# Create weather-CI relationship visualization for a few sample fields +tryCatch({ + # Get top fields in different health categories to show as examples + sample_fields <- field_health_scores %>% + dplyr::group_by(health_status) %>% + dplyr::slice_head(n = 1) %>% + dplyr::ungroup() %>% + dplyr::pull(field) + + # If we have more than 3 fields, just show 3 for brevity + if(length(sample_fields) > 3) { + sample_fields <- sample_fields[1:3] + } + + # If no sample fields are available, use the first field in the data + if(length(sample_fields) == 0 && nrow(AllPivots_merged) > 0) { + sample_fields <- AllPivots_merged$field[1] + } + + # Create weather plots for each sample field + for(field_name in sample_fields) { + # Create the weather-CI plot + weather_plot <- create_weather_ci_plot( + pivotName = field_name, + ci_data = CI_quadrant, + days_to_show = 60 # Show last 60 days + ) + + # Print the plot + print(weather_plot) + } + + # Add explanation if using mock weather data + cat("*Note: Weather data shown is representative and may vary from actual field conditions.*\n") + cat("*For production use, this would be connected to local weather station data or weather APIs.*\n") + +}, error = function(e) { + safe_log(paste("Error creating weather integration:", e$message), "ERROR") + cat("Error generating weather relationship visualization. See log for details.") +}) +``` + +\newpage + +```{r generate_field_visualizations, eval=TRUE, fig.height=3.8, fig.width=10, message=FALSE,echo=FALSE, warning=FALSE, include=TRUE, results='asis'} +# Generate detailed visualizations for each field +tryCatch({ + # Merge field polygons for processing + AllPivots_merged <- AllPivots0 %>% + dplyr::group_by(field) %>% + dplyr::summarise(.groups = 'drop') + + # Log start time for performance measurement + start_time <- Sys.time() + safe_log(paste("Starting field visualization generation at", start_time)) + + # Setup progress tracking + p <- progressr::progressor(steps = nrow(AllPivots_merged)) + + # Generate field-specific visualizations + field_results <- furrr::future_map(AllPivots_merged$field, function(field_name) { + tryCatch({ + # Update progress + p(sprintf("Processing field %s", field_name)) + + # Temporary list to store outputs + outputs <- list() + outputs$field_name <- field_name + + # Get field data + field_shape <- AllPivots0 %>% dplyr::filter(field == field_name) + + # Add irrigation priority map for this field + field_ci <- terra::crop(CI, field_shape) %>% terra::mask(., field_shape) + field_change <- terra::crop(last_week_dif_raster_abs, field_shape) %>% terra::mask(., field_shape) + + # Store plot objects + outputs$irrigation_map <- create_irrigation_map( + field_ci, + field_change, + field_shape, + title = paste("Field", field_name, "- Irrigation Priority") + ) + + return(outputs) + }, error = function(e) { + safe_log(paste("Error processing field visualization for", field_name, ":", e$message), "ERROR") + return(list( + field_name = field_name, + error = e$message + )) + }) + }, .options = furrr::furrr_options(seed = TRUE)) + + # Log performance metrics + end_time <- Sys.time() + processing_time <- as.numeric(difftime(end_time, start_time, units="secs")) + safe_log(paste("Field visualization processing completed in", round(processing_time, 2), "seconds")) + safe_log(paste("Average time per field:", round(processing_time / nrow(AllPivots_merged), 2), "seconds")) + + # Generate detailed plots for each field using standard sequential processing + # This part still uses sequential processing because the ci_plot function renders directly to the document + safe_log("Starting sequential rendering of field plots") + purrr::walk(AllPivots_merged$field, function(field_name) { + tryCatch({ + cat("\n") # Add an empty line for better spacing + + # First, add field header and retrieve the field-specific irrigation map + cat(paste("## Field", field_name, "\n\n")) + + # Find the irrigation map for this field + field_result <- field_results[[which(AllPivots_merged$field == field_name)]] + + # If we have irrigation data for this field, show it + if (!is.null(field_result$irrigation_map)) { + cat("\n### Irrigation Priority Map\n\n") + print(field_result$irrigation_map) + cat("\n") + } + + # Call ci_plot with explicit parameters + ci_plot( + pivotName = field_name, + 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, + use_breaks = use_breaks, + borders = borders + ) + + cat("\n") + + # Call cum_ci_plot with explicit parameters + cum_ci_plot( + pivotName = field_name, + ci_quadrant_data = CI_quadrant, + plot_type = "value", + facet_on = FALSE + ) + + }, error = function(e) { + safe_log(paste("Error generating plots for field", field_name, ":", e$message), "ERROR") + cat(paste("## Error generating plots for field", field_name, "\n")) + cat(paste(e$message, "\n")) + }) + }) + +}, error = function(e) { + safe_log(paste("Error in field visualization section:", e$message), "ERROR") + cat("Error generating field plots. See log for details.\n") +}) +``` + +```{r generate_subarea_visualizations, echo=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, results='asis', eval=FALSE} +# Alternative visualization grouped by sub-area (disabled by default) +tryCatch({ + # Group pivots by sub-area + pivots_grouped <- AllPivots0 + + # Iterate over each subgroup + for (subgroup in unique(pivots_grouped$sub_area)) { + # Add subgroup heading + cat("\n") + cat("## Subgroup: ", subgroup, "\n") + + # Filter data for current subgroup + subset_data <- dplyr::filter(pivots_grouped, sub_area == subgroup) + + # Generate visualizations for each field in the subgroup + purrr::walk(subset_data$field, function(field_name) { + cat("\n") + ci_plot(field_name) + cat("\n") + cum_ci_plot(field_name) + cat("\n") + }) + + # Add page break after each subgroup + cat("\\pagebreak\n") + } +}, error = function(e) { + safe_log(paste("Error in subarea visualization section:", e$message), "ERROR") + cat("Error generating subarea plots. See log for details.\n") +}) +``` + +# Yield prediction +The below table shows estimates of the biomass if you would harvest them now. + +```{r yield_data_training, message=FALSE, warning=FALSE, include=FALSE} +# Load and prepare yield prediction data with error handling +tryCatch({ + # Load CI quadrant data and fill missing values + CI_quadrant <- readRDS(here::here(cumulative_CI_vals_dir, "All_pivots_Cumulative_CI_quadrant_year_v2.rds")) %>% + dplyr::group_by(model) %>% + tidyr::fill(field, sub_field, .direction = "downup") %>% + dplyr::ungroup() + + # Check if tonnage_ha is empty + if (all(is.na(CI_quadrant$tonnage_ha))) { + safe_log("Lacking historic harvest data, please provide for yield prediction calculation", "WARNING") + knitr::knit_exit() # Exit the chunk if tonnage_ha is empty + } + + # Rename year column to season for consistency + harvesting_data <- harvesting_data %>% dplyr::rename(season = year) + + # Join CI and yield data + CI_and_yield <- dplyr::left_join(CI_quadrant, harvesting_data, by = c("field", "sub_field", "season")) %>% + dplyr::group_by(sub_field, season) %>% + dplyr::slice(which.max(DOY)) %>% + dplyr::select(field, sub_field, tonnage_ha, cumulative_CI, DOY, season, sub_area) %>% + dplyr::mutate(CI_per_day = cumulative_CI / DOY) + + # Define predictors and response variables + predictors <- c("cumulative_CI", "DOY", "CI_per_day") + response <- "tonnage_ha" + + # Prepare test and validation datasets + CI_and_yield_test <- CI_and_yield %>% + as.data.frame() %>% + dplyr::filter(!is.na(tonnage_ha)) + + CI_and_yield_validation <- CI_and_yield_test + + # Prepare prediction dataset (fields without harvest data) + prediction_yields <- CI_and_yield %>% + as.data.frame() %>% + dplyr::filter(is.na(tonnage_ha)) + + # Configure model training parameters + ctrl <- caret::trainControl( + method = "cv", + savePredictions = TRUE, + allowParallel = TRUE, + number = 5, + verboseIter = TRUE + ) + + # Train the model with feature selection + set.seed(202) # For reproducibility + model_ffs_rf <- CAST::ffs( + CI_and_yield_test[, predictors], + CI_and_yield_test[, response], + method = "rf", + trControl = ctrl, + importance = TRUE, + withinSE = TRUE, + tuneLength = 5, + na.rm = TRUE + ) + + # Function to prepare predictions with consistent naming and formatting + prepare_predictions <- function(predictions, newdata) { + return(predictions %>% + as.data.frame() %>% + dplyr::rename(predicted_Tcha = ".") %>% + dplyr::mutate( + sub_field = newdata$sub_field, + field = newdata$field, + Age_days = newdata$DOY, + total_CI = round(newdata$cumulative_CI, 0), + predicted_Tcha = round(predicted_Tcha, 0), + season = newdata$season + ) %>% + dplyr::select(field, sub_field, Age_days, total_CI, predicted_Tcha, season) %>% + dplyr::left_join(., newdata, by = c("field", "sub_field", "season")) + ) + } + + # Predict yields for the validation dataset + pred_ffs_rf <- prepare_predictions(stats::predict(model_ffs_rf, newdata = CI_and_yield_validation), CI_and_yield_validation) + + # Predict yields for the current season (focus on mature fields over 300 days) + pred_rf_current_season <- prepare_predictions(stats::predict(model_ffs_rf, newdata = prediction_yields), prediction_yields) %>% + dplyr::filter(Age_days > 300) %>% + dplyr::mutate(CI_per_day = round(total_CI / Age_days, 1)) + + safe_log("Successfully completed yield prediction calculations") + +}, error = function(e) { + safe_log(paste("Error in yield prediction:", e$message), "ERROR") + # Create empty dataframes to prevent errors in subsequent chunks + pred_ffs_rf <- data.frame() + pred_rf_current_season <- data.frame() +}) +``` + +```{r plotting_yield_data, echo=FALSE, fig.height=5, fig.width=8, message=FALSE, warning=FALSE} +# Display yield prediction visualizations with error handling +tryCatch({ + if (nrow(pred_ffs_rf) > 0) { + # Plot model performance (predicted vs actual) + ggplot2::ggplot(pred_ffs_rf, ggplot2::aes(y = predicted_Tcha, x = tonnage_ha)) + + ggplot2::geom_point(size = 2, alpha = 0.6) + + ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "red") + + ggplot2::scale_x_continuous(limits = c(0, 200)) + + ggplot2::scale_y_continuous(limits = c(0, 200)) + + ggplot2::labs(title = "Model Performance: \nPredicted vs Actual Tonnage/ha", + x = "Actual tonnage/ha (Tcha)", + y = "Predicted tonnage/ha (Tcha)") + + ggplot2::theme_minimal() + } + + if (nrow(pred_rf_current_season) > 0) { + # Plot predicted yields by age + ggplot2::ggplot(pred_rf_current_season, ggplot2::aes(x = Age_days, y = predicted_Tcha)) + + ggplot2::geom_point(size = 2, alpha = 0.6) + + ggplot2::labs(title = "Predicted Yields for Fields Over 300 Days \nOld Yet to Be Harvested", + x = "Age (days)", + y = "Predicted tonnage/ha (Tcha)") + + ggplot2::scale_y_continuous(limits = c(0, 200)) + + ggplot2::theme_minimal() + + # Display prediction table + knitr::kable(pred_rf_current_season, + digits = 0, + caption = "Predicted Tonnage/ha for Fields Over 300 Days Old") + } else { + cat("No fields over 300 days old without harvest data available for yield prediction.") + } +}, error = function(e) { + safe_log(paste("Error in yield prediction visualization:", e$message), "ERROR") + cat("Error generating yield prediction visualizations. See log for details.") +}) +``` + diff --git a/r_app/CI_report_executive_summary.Rmd b/r_app/CI_report_executive_summary.Rmd new file mode 100644 index 0000000..1e21152 --- /dev/null +++ b/r_app/CI_report_executive_summary.Rmd @@ -0,0 +1,1463 @@ +--- +params: + ref: "word-styles-reference-var1.docx" + output_file: CI_report.docx + report_date: "2024-08-28" + data_dir: "Chemba" + mail_day: "Wednesday" + borders: TRUE + use_breaks: FALSE +output: + # html_document: + # toc: yes + # df_print: paged + word_document: + reference_docx: !expr file.path("word-styles-reference-var1.docx") + toc: yes +editor_options: + chunk_output_type: console +--- + +```{r setup_parameters, include=FALSE} +# Set up basic report parameters from input values +report_date <- params$report_date +mail_day <- params$mail_day +borders <- params$borders +use_breaks <- params$use_breaks # Whether to use breaks or continuous spectrum in visualizations + +# Environment setup notes (commented out) +# # Activeer de renv omgeving +# renv::activate() +# renv::deactivate() +# # Optioneel: Herstel de omgeving als dat nodig is +# # Je kunt dit commentaar geven als je het normaal niet wilt uitvoeren +# renv::restore() +``` + +```{r load_libraries, message=FALSE, warning=FALSE, include=FALSE} +# Configure knitr options +knitr::opts_chunk$set(warning = FALSE, message = FALSE) + +# Path management +library(here) + +# Spatial data libraries +library(sf) +library(terra) +library(exactextractr) +# library(raster) - Removed as it's no longer maintained + +# Data manipulation and visualization +library(tidyverse) # Includes dplyr, ggplot2, etc. +library(tmap) +library(lubridate) +library(zoo) + +# Machine learning +library(rsample) +library(caret) +library(randomForest) +library(CAST) + +# Load custom utility functions +tryCatch({ + source("report_utils.R") +}, error = function(e) { + message(paste("Error loading report_utils.R:", e$message)) + # Try alternative path if the first one fails + tryCatch({ + source(here::here("r_app", "report_utils.R")) + }, error = function(e) { + stop("Could not load report_utils.R from either location: ", e$message) + }) +}) +``` + +```{r initialize_project_config, message=FALSE, warning=FALSE, include=FALSE} +# Set the project directory from parameters +project_dir <- params$data_dir + +# Source project parameters with error handling +tryCatch({ + source(here::here("r_app", "parameters_project.R")) +}, error = function(e) { + stop("Error loading parameters_project.R: ", e$message) +}) + +# Log initial configuration +safe_log("Starting the R Markdown script") +safe_log(paste("mail_day params:", params$mail_day)) +safe_log(paste("report_date params:", params$report_date)) +safe_log(paste("mail_day variable:", mail_day)) +``` + +```{r calculate_dates_and_weeks, message=FALSE, warning=FALSE, include=FALSE} +# Set locale for consistent date formatting +Sys.setlocale("LC_TIME", "C") + +# Initialize date variables from parameters +today <- as.character(report_date) +mail_day_as_character <- as.character(mail_day) + +# Calculate week days +report_date_as_week_day <- weekdays(lubridate::ymd(today)) +days_of_week <- c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday") + +# Calculate initial week number +week <- lubridate::week(today) +safe_log(paste("Initial week calculation:", week, "today:", today)) + +# Calculate previous dates for comparisons +today_minus_1 <- as.character(lubridate::ymd(today) - 7) +today_minus_2 <- as.character(lubridate::ymd(today) - 14) +today_minus_3 <- as.character(lubridate::ymd(today) - 21) + +# Log the weekday calculations for debugging +safe_log(paste("Report date weekday:", report_date_as_week_day)) +safe_log(paste("Weekday index:", which(days_of_week == report_date_as_week_day))) +safe_log(paste("Mail day:", mail_day_as_character)) +safe_log(paste("Mail day index:", which(days_of_week == mail_day_as_character))) + +# Adjust week calculation based on mail day +if (which(days_of_week == report_date_as_week_day) > which(days_of_week == mail_day_as_character)) { + safe_log("Adjusting weeks because of mail day") + week <- lubridate::week(today) + 1 + today_minus_1 <- as.character(lubridate::ymd(today)) + today_minus_2 <- as.character(lubridate::ymd(today) - 7) + today_minus_3 <- as.character(lubridate::ymd(today) - 14) +} + +# Generate subtitle for report +subtitle_var <- paste("Report generated on", Sys.Date()) + +# Calculate week numbers for previous weeks +week_minus_1 <- week - 1 +week_minus_2 <- week - 2 +week_minus_3 <- week - 3 + +# Format current week with leading zeros +week <- sprintf("%02d", week) + +# Get years for each date +year <- lubridate::year(today) +year_1 <- lubridate::year(today_minus_1) +year_2 <- lubridate::year(today_minus_2) +year_3 <- lubridate::year(today_minus_3) +``` + +```{r data, message=TRUE, warning=TRUE, include=FALSE} +# Load CI index data with error handling +tryCatch({ + CI_quadrant <- readRDS(here::here(cumulative_CI_vals_dir, "All_pivots_Cumulative_CI_quadrant_year_v2.rds")) + safe_log("Successfully loaded CI quadrant data") +}, error = function(e) { + stop("Error loading CI quadrant data: ", e$message) +}) + +# Get file paths for different weeks using the utility function +tryCatch({ + path_to_week_current = get_week_path(weekly_CI_mosaic, today, 0) + path_to_week_minus_1 = get_week_path(weekly_CI_mosaic, today, -1) + path_to_week_minus_2 = get_week_path(weekly_CI_mosaic, today, -2) + path_to_week_minus_3 = get_week_path(weekly_CI_mosaic, today, -3) + + # Log the calculated paths + safe_log("Required mosaic paths:") + safe_log(paste("Path to current week:", path_to_week_current)) + safe_log(paste("Path to week minus 1:", path_to_week_minus_1)) + safe_log(paste("Path to week minus 2:", path_to_week_minus_2)) + safe_log(paste("Path to week minus 3:", path_to_week_minus_3)) + + # Validate that files exist + if (!file.exists(path_to_week_current)) warning("Current week mosaic file does not exist: ", path_to_week_current) + if (!file.exists(path_to_week_minus_1)) warning("Week minus 1 mosaic file does not exist: ", path_to_week_minus_1) + if (!file.exists(path_to_week_minus_2)) warning("Week minus 2 mosaic file does not exist: ", path_to_week_minus_2) + if (!file.exists(path_to_week_minus_3)) warning("Week minus 3 mosaic file does not exist: ", path_to_week_minus_3) + + # Load raster data with terra functions + CI <- terra::rast(path_to_week_current)$CI + CI_m1 <- terra::rast(path_to_week_minus_1)$CI + CI_m2 <- terra::rast(path_to_week_minus_2)$CI + CI_m3 <- terra::rast(path_to_week_minus_3)$CI + +}, error = function(e) { + stop("Error loading raster data: ", e$message) +}) +``` + +```{r calculate_difference_rasters, message=TRUE, warning=TRUE, include=FALSE} +# Calculate difference rasters for comparisons +tryCatch({ + # Calculate weekly difference + last_week_dif_raster_abs <- (CI - CI_m1) + safe_log("Calculated weekly difference raster") + + # Calculate three-week difference + three_week_dif_raster_abs <- (CI - CI_m3) + safe_log("Calculated three-week difference raster") +}, error = function(e) { + safe_log(paste("Error calculating difference rasters:", e$message), "ERROR") + # Create placeholder rasters if calculations fail + if (!exists("last_week_dif_raster_abs")) { + last_week_dif_raster_abs <- CI * 0 + } + if (!exists("three_week_dif_raster_abs")) { + three_week_dif_raster_abs <- CI * 0 + } +}) +``` + +```{r load_field_boundaries, message=TRUE, warning=TRUE, include=FALSE} +# Load field boundaries from parameters +tryCatch({ + AllPivots0 <- field_boundaries_sf + safe_log("Successfully loaded field boundaries") +}, error = function(e) { + stop("Error loading field boundaries: ", e$message) +}) +``` + +```{r create_farm_health_data, message=FALSE, warning=FALSE, include=FALSE} +# Create farm health summary data from scratch +tryCatch({ + # Ensure we have the required data + if (!exists("AllPivots0") || !exists("CI") || !exists("CI_m1") || !exists("harvesting_data")) { + stop("Required input data (field boundaries, CI data, or harvesting data) not available") + } + + safe_log("Starting to calculate farm health data") + + # Get unique field names + fields <- unique(AllPivots0$field) + safe_log(paste("Found", length(fields), "unique fields")) + + # Initialize result dataframe + farm_health_data <- data.frame( + field = character(), + mean_ci = numeric(), + ci_change = numeric(), + ci_uniformity = numeric(), + status = character(), + anomaly_type = character(), + priority_level = numeric(), + age_weeks = numeric(), + harvest_readiness = character(), + stringsAsFactors = FALSE + ) + + # Process each field with robust error handling + for (field_name in fields) { + tryCatch({ + safe_log(paste("Processing field:", field_name)) + + # Get field boundary + field_shape <- AllPivots0 %>% dplyr::filter(field == field_name) + + # Skip if field shape is empty + if (nrow(field_shape) == 0) { + safe_log(paste("Empty field shape for", field_name), "WARNING") + next + } + + # Get field age from harvesting data - use direct filtering to avoid dplyr errors + field_age_data <- NULL + if (exists("harvesting_data") && !is.null(harvesting_data) && nrow(harvesting_data) > 0) { + field_age_data <- harvesting_data[harvesting_data$field == field_name, ] + if (nrow(field_age_data) > 0) { + field_age_data <- field_age_data[order(field_age_data$season_start, decreasing = TRUE), ][1, ] + } + } + + # Default age if not available + field_age_weeks <- if (!is.null(field_age_data) && nrow(field_age_data) > 0 && !is.na(field_age_data$age)) { + field_age_data$age + } else { + 10 # Default age + } + + # Extract CI values using terra's extract function which is more robust + ci_values <- terra::extract(CI, field_shape) + ci_prev_values <- terra::extract(CI_m1, field_shape) + + # Check if we got valid data + if (nrow(ci_values) == 0 || nrow(ci_prev_values) == 0) { + safe_log(paste("No CI data extracted for field", field_name), "WARNING") + # Add a placeholder row with Unknown status + farm_health_data <- rbind(farm_health_data, data.frame( + field = field_name, + mean_ci = NA, + ci_change = NA, + ci_uniformity = NA, + status = "Unknown", + anomaly_type = "Unknown", + priority_level = 5, # Low priority + age_weeks = field_age_weeks, + harvest_readiness = "Unknown", + stringsAsFactors = FALSE + )) + next + } + + # Calculate metrics - Handle NA values properly + ci_column <- if ("CI" %in% names(ci_values)) "CI" else colnames(ci_values)[1] + ci_prev_column <- if ("CI" %in% names(ci_prev_values)) "CI" else colnames(ci_prev_values)[1] + + mean_ci <- mean(ci_values[[ci_column]], na.rm=TRUE) + mean_ci_prev <- mean(ci_prev_values[[ci_prev_column]], na.rm=TRUE) + ci_change <- mean_ci - mean_ci_prev + ci_sd <- sd(ci_values[[ci_column]], na.rm=TRUE) + ci_uniformity <- ci_sd / max(0.1, mean_ci) # Avoid division by zero + + # Handle NaN or Inf results + if (is.na(mean_ci) || is.na(ci_change) || is.na(ci_uniformity) || + is.nan(mean_ci) || is.nan(ci_change) || is.nan(ci_uniformity) || + is.infinite(mean_ci) || is.infinite(ci_change) || is.infinite(ci_uniformity)) { + safe_log(paste("Invalid calculation results for field", field_name), "WARNING") + # Add a placeholder row with Unknown status + farm_health_data <- rbind(farm_health_data, data.frame( + field = field_name, + mean_ci = NA, + ci_change = NA, + ci_uniformity = NA, + status = "Unknown", + anomaly_type = "Unknown", + priority_level = 5, # Low priority + age_weeks = field_age_weeks, + harvest_readiness = "Unknown", + stringsAsFactors = FALSE + )) + next + } + + # Determine field status + status <- dplyr::case_when( + mean_ci >= 5 ~ "Excellent", + mean_ci >= 3.5 ~ "Good", + mean_ci >= 2 ~ "Fair", + mean_ci >= 1 ~ "Poor", + TRUE ~ "Critical" + ) + + # Determine anomaly type + anomaly_type <- dplyr::case_when( + ci_change > 2 ~ "Potential Weed Growth", + ci_change < -2 ~ "Potential Weeding/Harvesting", + ci_uniformity > 0.5 ~ "High Variability", + mean_ci < 1 ~ "Low Vigor", + TRUE ~ "None" + ) + + # Calculate priority level (1-5, with 1 being highest priority) + priority_score <- dplyr::case_when( + mean_ci < 1 ~ 1, # Critical - highest priority + anomaly_type == "Potential Weed Growth" ~ 2, + anomaly_type == "High Variability" ~ 3, + ci_change < -1 ~ 4, + TRUE ~ 5 # No urgent issues + ) + + # Determine harvest readiness + harvest_readiness <- dplyr::case_when( + field_age_weeks >= 52 & mean_ci >= 4 ~ "Ready for harvest", + field_age_weeks >= 48 & mean_ci >= 3.5 ~ "Approaching harvest", + field_age_weeks >= 40 & mean_ci >= 3 ~ "Mid-maturity", + field_age_weeks >= 12 ~ "Growing", + TRUE ~ "Early stage" + ) + + # Add to summary data + farm_health_data <- rbind(farm_health_data, data.frame( + field = field_name, + mean_ci = round(mean_ci, 2), + ci_change = round(ci_change, 2), + ci_uniformity = round(ci_uniformity, 2), + status = status, + anomaly_type = anomaly_type, + priority_level = priority_score, + age_weeks = field_age_weeks, + harvest_readiness = harvest_readiness, + stringsAsFactors = FALSE + )) + + }, error = function(e) { + safe_log(paste("Error processing field", field_name, ":", e$message), "ERROR") + # Add a placeholder row with Error status + farm_health_data <<- rbind(farm_health_data, data.frame( + field = field_name, + mean_ci = NA, + ci_change = NA, + ci_uniformity = NA, + status = "Unknown", + anomaly_type = "Unknown", + priority_level = 5, # Low priority since we don't know the status + age_weeks = NA, + harvest_readiness = "Unknown", + stringsAsFactors = FALSE + )) + }) + } + + # Make sure we have data for all fields + if (nrow(farm_health_data) == 0) { + safe_log("No farm health data was created", "ERROR") + stop("Failed to create farm health data") + } + + # Sort by priority level + farm_health_data <- farm_health_data %>% dplyr::arrange(priority_level, field) + + safe_log(paste("Successfully created farm health data for", nrow(farm_health_data), "fields")) + +}, error = function(e) { + safe_log(paste("Error creating farm health data:", e$message), "ERROR") + # Create an empty dataframe that can be filled by the verification chunk +}) +``` + +```{r verify_farm_health_data, message=FALSE, warning=FALSE, include=FALSE} +# Verify farm_health_data exists and has content +if (!exists("farm_health_data") || nrow(farm_health_data) == 0) { + safe_log("farm_health_data not found or empty, generating default data", "WARNING") + + # Create minimal fallback data + tryCatch({ + # Get fields from boundaries + fields <- unique(AllPivots0$field) + + # Create basic data frame with just field names + farm_health_data <- data.frame( + field = fields, + mean_ci = rep(NA, length(fields)), + ci_change = rep(NA, length(fields)), + ci_uniformity = rep(NA, length(fields)), + status = rep("Unknown", length(fields)), + anomaly_type = rep("Unknown", length(fields)), + priority_level = rep(5, length(fields)), # Low priority + age_weeks = rep(NA, length(fields)), + harvest_readiness = rep("Unknown", length(fields)), + stringsAsFactors = FALSE + ) + + safe_log("Created fallback farm_health_data with basic field information") + }, error = function(e) { + safe_log(paste("Error creating fallback farm_health_data:", e$message), "ERROR") + farm_health_data <<- data.frame( + field = character(), + mean_ci = numeric(), + ci_change = numeric(), + ci_uniformity = numeric(), + status = character(), + anomaly_type = character(), + priority_level = numeric(), + age_weeks = numeric(), + harvest_readiness = character(), + stringsAsFactors = FALSE + ) + }) +} +``` + +```{r calculate_farm_health, message=FALSE, warning=FALSE, include=FALSE} +# Calculate farm health summary metrics +tryCatch({ + # Generate farm health summary data + farm_health_data <- generate_farm_health_summary( + field_boundaries = AllPivots0, + ci_current = CI, + ci_previous = CI_m1, + harvesting_data = harvesting_data + ) + + # Log the summary data + safe_log(paste("Generated farm health summary with", nrow(farm_health_data), "fields")) + +}, error = function(e) { + safe_log(paste("Error in farm health calculation:", e$message), "ERROR") + # Create empty dataframe if calculation failed + farm_health_data <- data.frame( + field = character(), + mean_ci = numeric(), + ci_change = numeric(), + ci_uniformity = numeric(), + status = character(), + anomaly_type = character(), + priority_level = numeric(), + age_weeks = numeric(), + harvest_readiness = character(), + stringsAsFactors = FALSE + ) +}) +``` + +```{r executive_summary_functions, message=FALSE, warning=FALSE, include=FALSE} +# EXECUTIVE SUMMARY HELPER FUNCTIONS + +#' Generate a summary of farm health status +#' +#' @param field_boundaries Field boundaries spatial data (sf object) +#' @param ci_current Current CI raster +#' @param ci_previous Previous week's CI raster +#' @param harvesting_data Data frame with harvesting information +#' @return A data frame with farm status summary metrics +#' +generate_farm_health_summary <- function(field_boundaries, ci_current, ci_previous, harvesting_data) { + # Generate a summary data frame of farm health by field + tryCatch({ + # Get unique field names + fields <- unique(field_boundaries$field) + + # Initialize result dataframe + summary_data <- data.frame( + field = character(), + mean_ci = numeric(), + ci_change = numeric(), + ci_uniformity = numeric(), + status = character(), + anomaly_type = character(), + priority_level = numeric(), + age_weeks = numeric(), + harvest_readiness = character(), + stringsAsFactors = FALSE + ) + + # Process each field with better error handling + for (field_name in fields) { + tryCatch({ + # Get field boundary + field_shape <- field_boundaries %>% dplyr::filter(field == field_name) + + # Skip if field shape is empty + if (nrow(field_shape) == 0) { + safe_log(paste("Empty field shape for", field_name), "WARNING") + next + } + + # Get field age from harvesting data + field_age_data <- harvesting_data %>% + dplyr::filter(field == field_name) %>% + dplyr::arrange(desc(season_start)) %>% + dplyr::slice(1) + + # Default age if not available + field_age_weeks <- if (nrow(field_age_data) > 0 && !is.na(field_age_data$age)) { + field_age_data$age + } else { + 10 # Default age + } + + # Extract CI values for this field using extract instead of crop/mask to avoid pointer issues + # This is more robust than the crop+mask approach + field_bbox <- sf::st_bbox(field_shape) + extent_vec <- c(field_bbox$xmin, field_bbox$xmax, field_bbox$ymin, field_bbox$ymax) + + # Use terra extract function instead of crop+mask + ci_values <- terra::extract(ci_current, field_shape) + ci_prev_values <- terra::extract(ci_previous, field_shape) + + # Calculate metrics + mean_ci <- mean(ci_values$CI, na.rm=TRUE) + mean_ci_prev <- mean(ci_prev_values$CI, na.rm=TRUE) + ci_change <- mean_ci - mean_ci_prev + ci_sd <- sd(ci_values$CI, na.rm=TRUE) + ci_uniformity <- ci_sd / max(0.1, mean_ci) # Avoid division by zero + + # Determine field status + status <- dplyr::case_when( + mean_ci >= 5 ~ "Excellent", + mean_ci >= 3.5 ~ "Good", + mean_ci >= 2 ~ "Fair", + mean_ci >= 1 ~ "Poor", + TRUE ~ "Critical" + ) + + # Determine anomaly type + anomaly_type <- dplyr::case_when( + ci_change > 2 ~ "Potential Weed Growth", + ci_change < -2 ~ "Potential Weeding/Harvesting", + ci_uniformity > 0.5 ~ "High Variability", + mean_ci < 1 ~ "Low Vigor", + TRUE ~ "None" + ) + + # Calculate priority level (1-5, with 1 being highest priority) + priority_score <- dplyr::case_when( + mean_ci < 1 ~ 1, # Critical - highest priority + anomaly_type == "Potential Weed Growth" ~ 2, + anomaly_type == "High Variability" ~ 3, + ci_change < -1 ~ 4, + TRUE ~ 5 # No urgent issues + ) + + # Determine harvest readiness + harvest_readiness <- dplyr::case_when( + field_age_weeks >= 52 & mean_ci >= 4 ~ "Ready for harvest", + field_age_weeks >= 48 & mean_ci >= 3.5 ~ "Approaching harvest", + field_age_weeks >= 40 & mean_ci >= 3 ~ "Mid-maturity", + field_age_weeks >= 12 ~ "Growing", + TRUE ~ "Early stage" + ) + + # Add to summary data + summary_data <- rbind(summary_data, data.frame( + field = field_name, + mean_ci = round(mean_ci, 2), + ci_change = round(ci_change, 2), + ci_uniformity = round(ci_uniformity, 2), + status = status, + anomaly_type = anomaly_type, + priority_level = priority_score, + age_weeks = field_age_weeks, + harvest_readiness = harvest_readiness, + stringsAsFactors = FALSE + )) + + }, error = function(e) { + safe_log(paste("Error calculating health score for field", field_name, ":", e$message), "ERROR") + # Add a row with NA values for this field to ensure it still appears in outputs + summary_data <- rbind(summary_data, data.frame( + field = field_name, + mean_ci = NA, + ci_change = NA, + ci_uniformity = NA, + status = "Error", + anomaly_type = "Error", + priority_level = 1, # High priority because it needs investigation + age_weeks = NA, + harvest_readiness = "Unknown", + stringsAsFactors = FALSE + )) + }) + } + + # Sort by priority level + summary_data <- summary_data %>% dplyr::arrange(priority_level, field) + + return(summary_data) + + }, error = function(e) { + safe_log(paste("Error in generate_farm_health_summary:", e$message), "ERROR") + return(data.frame( + field = character(), + mean_ci = numeric(), + ci_change = numeric(), + ci_uniformity = numeric(), + status = character(), + anomaly_type = character(), + priority_level = numeric(), + age_weeks = numeric(), + harvest_readiness = character(), + stringsAsFactors = FALSE + )) + }) +} + +#' Create a farm-wide anomaly detection map +#' +#' @param ci_current Current CI raster +#' @param ci_previous Previous week's CI raster +#' @param field_boundaries Field boundaries spatial data (sf object) +#' @return A tmap object with anomaly visualization +#' +create_anomaly_map <- function(ci_current, ci_previous, field_boundaries) { + tryCatch({ + # Calculate difference raster + ci_diff <- ci_current - ci_previous + + # Create a categorical raster for anomalies + anomaly_raster <- ci_current * 0 # Initialize with same extent/resolution + + # Extract values to manipulate + diff_values <- terra::values(ci_diff) + curr_values <- terra::values(ci_current) + + # Define anomaly categories: + # 4: Significant growth (potential weeds) - CI increase > 2 + # 3: Moderate growth - CI increase 1-2 + # 2: Stable - CI change between -1 and 1 + # 1: Moderate decline - CI decrease 1-2 + # 0: Significant decline (potential weeding/harvesting) - CI decrease > 2 + + # Apply classification + anomaly_values <- rep(NA, length(diff_values)) + + # Significant growth (potential weeds) + sig_growth <- which(diff_values > 2 & !is.na(diff_values)) + anomaly_values[sig_growth] <- 4 + + # Moderate growth + mod_growth <- which(diff_values > 1 & diff_values <= 2 & !is.na(diff_values)) + anomaly_values[mod_growth] <- 3 + + # Stable + stable <- which(diff_values >= -1 & diff_values <= 1 & !is.na(diff_values)) + anomaly_values[stable] <- 2 + + # Moderate decline + mod_decline <- which(diff_values < -1 & diff_values >= -2 & !is.na(diff_values)) + anomaly_values[mod_decline] <- 1 + + # Significant decline (potential weeding) + sig_decline <- which(diff_values < -2 & !is.na(diff_values)) + anomaly_values[sig_decline] <- 0 + + # Set values in raster + terra::values(anomaly_raster) <- anomaly_values + + # Create anomaly map + map <- tm_shape(anomaly_raster) + + tm_raster( + style = "cat", + palette = c("#d73027", "#fc8d59", "#ffffbf", "#91cf60", "#1a9850"), + labels = c("Significant Decline", "Moderate Decline", "Stable", "Moderate Growth", "Significant Growth"), + title = "Weekly CI Change" + ) + + tm_shape(field_boundaries) + + tm_borders(col = "black", lwd = 1.5) + + tm_text("field", size = 0.6) + + tm_layout( + main.title = "Farm-Wide Anomaly Detection", + legend.outside = TRUE, + legend.outside.position = "bottom" + ) + + tm_scale_bar(position = tm_pos_out("right", "bottom")) + + return(map) + + }, error = function(e) { + safe_log(paste("Error in create_anomaly_map:", e$message), "ERROR") + return(NULL) + }) +} + +#' Create a choropleth map of field health status +#' +#' @param field_boundaries Field boundaries with health data +#' @param attribute Field to visualize (e.g., "priority_level", "mean_ci") +#' @param title Map title +#' @param palette Color palette to use +#' @param legend_title Legend title +#' @return A tmap object +#' +create_field_status_map <- function(field_boundaries, health_data, attribute, + title = "Field Status Overview", + palette = "RdYlGn", + legend_title = "Status") { + tryCatch({ + # Join health data to field boundaries + field_data <- field_boundaries %>% + dplyr::left_join(health_data, by = "field") + + # Create style based on attribute type + if (attribute == "status") { + # Categorical styling for status + map <- tm_shape(field_data) + + tm_fill( + col = attribute, + palette = c("Critical" = "#d73027", "Poor" = "#fc8d59", + "Fair" = "#ffffbf", "Good" = "#91cf60", "Excellent" = "#1a9850", + "Error" = "#999999"), # Added Error category + title = legend_title + ) + } else if (attribute == "priority_level") { + # Numeric with custom breaks for priority (5 to 1, with 1 being highest priority) + map <- tm_shape(field_data) + + tm_fill( + col = attribute, + palette = "-RdYlGn", # Reversed so red is high priority + breaks = c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5), + labels = c("Critical", "High", "Medium", "Low", "Minimal"), + title = legend_title + ) + } else if (attribute == "anomaly_type") { + # Categorical styling for anomalies + map <- tm_shape(field_data) + + tm_fill( + col = attribute, + palette = c("Potential Weed Growth" = "#d73027", + "Potential Weeding/Harvesting" = "#4575b4", + "High Variability" = "#f46d43", + "Low Vigor" = "#fee090", + "None" = "#91cf60", + "Error" = "#999999"), # Added Error category + title = legend_title + ) + } else if (attribute == "harvest_readiness") { + # Categorical styling for harvest readiness + map <- tm_shape(field_data) + + tm_fill( + col = attribute, + palette = c("Ready for harvest" = "#1a9850", + "Approaching harvest" = "#91cf60", + "Mid-maturity" = "#ffffbf", + "Growing" = "#fc8d59", + "Early stage" = "#d73027", + "Unknown" = "#999999"), # Added Unknown category + title = legend_title + ) + } else { + # Default numerical styling + map <- tm_shape(field_data) + + tm_fill( + col = attribute, + palette = palette, + title = legend_title, + style = "cont", + na.color = "#999999" # Color for NA values + ) + } + + # Complete the map with borders and labels + map <- map + + tm_borders(col = "black", lwd = 1) + + tm_text("field", size = 0.7) + + tm_layout( + main.title = title, + legend.outside = TRUE, + legend.outside.position = "bottom" + ) + + tm_scale_bar(position = tm_pos_out("right", "bottom")) + + return(map) + + }, error = function(e) { + safe_log(paste("Error in create_field_status_map:", e$message), "ERROR") + return(NULL) + }) +} + +#' Create a summary statistics visualization +#' +#' @param health_data Farm health summary data +#' @return A ggplot2 object +#' +create_summary_stats <- function(health_data) { + tryCatch({ + # Handle empty dataframe case + if (nrow(health_data) == 0) { + return(ggplot2::ggplot() + + ggplot2::annotate("text", x = 0, y = 0, label = "No field data available") + + ggplot2::theme_void()) + } + + # Count fields by status + status_counts <- health_data %>% + dplyr::group_by(status) %>% + dplyr::summarise(count = n()) %>% + dplyr::mutate(status = factor(status, levels = c("Excellent", "Good", "Fair", "Poor", "Critical", "Error"))) + + # Create colors for status categories + status_colors <- c( + "Excellent" = "#1a9850", + "Good" = "#91cf60", + "Fair" = "#ffffbf", + "Poor" = "#fc8d59", + "Critical" = "#d73027", + "Error" = "#999999" + ) + + # Create bar chart + p <- ggplot2::ggplot(status_counts, ggplot2::aes(x = status, y = count, fill = status)) + + ggplot2::geom_bar(stat = "identity") + + ggplot2::scale_fill_manual(values = status_colors) + + ggplot2::labs( + title = "Field Status Summary", + x = "Status", + y = "Number of Fields", + fill = "Field Status" + ) + + ggplot2::theme_minimal() + + ggplot2::theme( + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + legend.position = "bottom" + ) + + return(p) + + }, error = function(e) { + safe_log(paste("Error in create_summary_stats:", e$message), "ERROR") + return(ggplot2::ggplot() + + ggplot2::annotate("text", x = 0, y = 0, label = paste("Error:", e$message)) + + ggplot2::theme_void()) + }) +} + +#' Create a bar chart of fields requiring attention +#' +#' @param health_data Farm health summary data +#' @param max_fields Maximum number of fields to display +#' @return A ggplot2 object +#' +create_priority_fields_chart <- function(health_data, max_fields = 10) { + tryCatch({ + # Handle empty dataframe case + if (nrow(health_data) == 0) { + return(ggplot2::ggplot() + + ggplot2::annotate("text", x = 0, y = 0, label = "No field data available") + + ggplot2::theme_void()) + } + + # Filter for fields that need attention (priority 1-3) + priority_fields <- health_data %>% + dplyr::filter(priority_level <= 3) %>% + dplyr::arrange(priority_level) %>% + dplyr::slice_head(n = max_fields) + + # If no priority fields, return message + if (nrow(priority_fields) == 0) { + return(ggplot2::ggplot() + + ggplot2::annotate("text", x = 0, y = 0, label = "No priority fields requiring attention") + + ggplot2::theme_void()) + } + + # Create priority labels + priority_fields$priority_label <- factor( + dplyr::case_when( + priority_fields$priority_level == 1 ~ "Critical", + priority_fields$priority_level == 2 ~ "High", + priority_fields$priority_level == 3 ~ "Medium", + TRUE ~ "Low" + ), + levels = c("Critical", "High", "Medium", "Low") + ) + + # Priority colors + priority_colors <- c( + "Critical" = "#d73027", + "High" = "#fc8d59", + "Medium" = "#fee090", + "Low" = "#91cf60" + ) + + # Create chart + p <- ggplot2::ggplot(priority_fields, + ggplot2::aes(x = reorder(field, -priority_level), + y = mean_ci, + fill = priority_label)) + + ggplot2::geom_bar(stat = "identity") + + ggplot2::geom_text(ggplot2::aes(label = anomaly_type), + position = ggplot2::position_stack(vjust = 0.5), + size = 3, angle = 90, hjust = 0) + + ggplot2::scale_fill_manual(values = priority_colors) + + ggplot2::labs( + title = "Priority Fields Requiring Attention", + subtitle = "With anomaly types and CI values", + x = "Field", + y = "Chlorophyll Index (CI)", + fill = "Priority" + ) + + ggplot2::theme_minimal() + + ggplot2::theme( + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + legend.position = "bottom" + ) + + return(p) + + }, error = function(e) { + safe_log(paste("Error in create_priority_fields_chart:", e$message), "ERROR") + return(ggplot2::ggplot() + + ggplot2::annotate("text", x = 0, y = 0, label = paste("Error:", e$message)) + + ggplot2::theme_void()) + }) +} + +#' Creates a harvest readiness visualization +#' +#' @param health_data Farm health summary data +#' @return A ggplot2 object +create_harvest_readiness_chart <- function(health_data) { + tryCatch({ + # Handle empty dataframe case + if (nrow(health_data) == 0) { + return(ggplot2::ggplot() + + ggplot2::annotate("text", x = 0, y = 0, label = "No field data available") + + ggplot2::theme_void()) + } + + # Count fields by harvest readiness + harvest_counts <- health_data %>% + dplyr::group_by(harvest_readiness) %>% + dplyr::summarise(count = n()) + + # Order factor levels + harvest_order <- c("Ready for harvest", "Approaching harvest", "Mid-maturity", "Growing", "Early stage", "Unknown") + harvest_counts$harvest_readiness <- factor(harvest_counts$harvest_readiness, levels = harvest_order) + + # Create colors for harvest readiness categories + harvest_colors <- c( + "Ready for harvest" = "#1a9850", + "Approaching harvest" = "#91cf60", + "Mid-maturity" = "#ffffbf", + "Growing" = "#fc8d59", + "Early stage" = "#d73027", + "Unknown" = "#999999" + ) + + # Create pie chart + p <- ggplot2::ggplot(harvest_counts, ggplot2::aes(x="", y=count, fill=harvest_readiness)) + + ggplot2::geom_bar(stat="identity", width=1) + + ggplot2::coord_polar("y", start=0) + + ggplot2::scale_fill_manual(values = harvest_colors) + + ggplot2::labs( + title = "Harvest Readiness Overview", + fill = "Harvest Stage" + ) + + ggplot2::theme_minimal() + + ggplot2::theme( + axis.title.x = ggplot2::element_blank(), + axis.title.y = ggplot2::element_blank(), + panel.border = ggplot2::element_blank(), + panel.grid = ggplot2::element_blank(), + axis.ticks = ggplot2::element_blank(), + axis.text = ggplot2::element_blank(), + plot.title = ggplot2::element_text(size=14, face="bold") + ) + + return(p) + + }, error = function(e) { + safe_log(paste("Error in create_harvest_readiness_chart:", e$message), "ERROR") + return(ggplot2::ggplot() + + ggplot2::annotate("text", x = 0, y = 0, label = paste("Error:", e$message)) + + ggplot2::theme_void()) + }) +} + +#' Generate recommendations based on farm health +#' +#' @param health_data Farm health summary data +#' @return HTML formatted recommendations +generate_executive_recommendations <- function(health_data) { + tryCatch({ + # Handle empty dataframe case + if (nrow(health_data) == 0) { + return("

Executive Recommendations

No field data available to generate recommendations.

") + } + + # Count fields by priority level + priority_counts <- health_data %>% + dplyr::group_by(priority_level) %>% + dplyr::summarise(count = n()) + + # Get critical and high priority fields + critical_fields <- health_data %>% + dplyr::filter(priority_level == 1) %>% + dplyr::pull(field) + + high_priority_fields <- health_data %>% + dplyr::filter(priority_level == 2) %>% + dplyr::pull(field) + + # Count harvest-ready fields + harvest_ready <- health_data %>% + dplyr::filter(harvest_readiness == "Ready for harvest") %>% + dplyr::pull(field) + + approaching_harvest <- health_data %>% + dplyr::filter(harvest_readiness == "Approaching harvest") %>% + dplyr::pull(field) + + # Count anomalies by type + anomaly_counts <- health_data %>% + dplyr::filter(anomaly_type != "None" & anomaly_type != "Error") %>% + dplyr::group_by(anomaly_type) %>% + dplyr::summarise(count = n()) + + # Generate HTML recommendations + html_output <- "
" + html_output <- paste0(html_output, "

Executive Recommendations

") + + # Priority recommendations + html_output <- paste0(html_output, "

Priority Actions:

") + + # Anomaly notifications + if (nrow(anomaly_counts) > 0) { + html_output <- paste0(html_output, "

Anomaly Notifications:

") + } + + # Farm status summary + html_output <- paste0(html_output, "

Farm Status Overview:

") + + return(html_output) + + }, error = function(e) { + safe_log(paste("Error in generate_executive_recommendations:", e$message), "ERROR") + return("

Error generating recommendations.

") + }) +} +``` + +`r subtitle_var` + +\pagebreak +# Explanation of the Report + +This report provides a detailed analysis of your sugarcane fields based on satellite imagery, helping you monitor crop health and development throughout the growing season. The data is processed weekly to give you timely insights for optimal farm management decisions. + +## What is the Chlorophyll Index (CI)? + +The **Chlorophyll Index (CI)** is a vegetation index that measures the relative amount of chlorophyll in plant leaves. Chlorophyll is the green pigment responsible for photosynthesis in plants. Higher CI values indicate: + +* Greater photosynthetic activity +* Healthier plant tissue +* Better nitrogen uptake +* More vigorous crop growth + +CI values typically range from 0 (bare soil or severely stressed vegetation) to 7+ (very healthy, dense vegetation). For sugarcane, values between 3-7 generally indicate good crop health, depending on the growth stage. + +# Executive Dashboard + +## Farm Health Status + +The map below shows the overall health status of all fields based on current Chlorophyll Index values. This provides a quick overview of which areas of your farm are performing well and which might need intervention. + +**How it works:** Field health status is determined by the average Chlorophyll Index (CI) value across each field: +- **Excellent** (dark green): CI ≥ 5.0 +- **Good** (light green): CI 3.5-4.99 +- **Fair** (yellow): CI 2.0-3.49 +- **Poor** (orange): CI 1.0-1.99 +- **Critical** (red): CI < 1.0 + +Fields with higher CI values indicate better crop vigor and photosynthetic activity, which typically correlate with healthier plants. + +```{r render_field_status_map, echo=FALSE, fig.height=6, fig.width=9, message=FALSE, warning=FALSE} +# Create field status map +tryCatch({ + # Create and display the field status map + field_status_map <- create_field_status_map( + field_boundaries = AllPivots0, + health_data = farm_health_data, + attribute = "status", + title = "Field Health Status Overview", + palette = "RdYlGn", + legend_title = "Health Status" + ) + + # Print the map + print(field_status_map) +}, error = function(e) { + safe_log(paste("Error creating field status map:", e$message), "ERROR") + plot(1, type="n", axes=FALSE, xlab="", ylab="") + text(1, 1, "Error creating field status map", cex=1.5) +}) +``` + +## Management Priorities + +This map highlights which fields require priority management attention based on current health indicators and trends. Fields in red require immediate attention, while green fields are performing well with minimal intervention needed. + +**How it works:** Priority levels are calculated based on a combination of factors: +- **Critical Priority** (dark red): Fields with CI < 1.0 or critical health issues +- **High Priority** (red): Fields with potential weed growth (CI increase > 2) +- **Medium Priority** (orange): Fields with high internal variability +- **Low Priority** (light green): Fields with moderate decline in CI +- **Minimal Priority** (dark green): Stable, healthy fields + +The priority algorithm considers both absolute CI values and week-to-week changes to identify fields that need immediate management attention. + +```{r render_priority_map, echo=FALSE, fig.height=6, fig.width=9, message=FALSE, warning=FALSE} +# Create priority management map +tryCatch({ + # Fix the priority mapping so red = high priority, green = low priority + # Reverse the priority levels before mapping (1=critical becomes 5, 5=minimal becomes 1) + farm_health_data$display_priority <- 6 - farm_health_data$priority_level + + # Create and display the priority map with corrected priority levels + priority_map <- tm_shape(AllPivots0 %>% dplyr::left_join(farm_health_data, by = "field")) + + tm_fill( + col = "display_priority", + palette = "RdYlGn", # Now properly oriented: red = high priority, green = low priority + breaks = c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5), + labels = c("Minimal", "Low", "Medium", "High", "Critical"), + title = "Priority Level" + ) + + tm_borders(col = "black", lwd = 1) + + tm_text("field", size = 0.7) + + tm_layout( + main.title = "Field Management Priority", + legend.outside = TRUE, + legend.outside.position = "bottom" + ) + + tm_scale_bar(position = tm_pos_out("right", "bottom")) + + # Print the map + print(priority_map) +}, error = function(e) { + safe_log(paste("Error creating priority map:", e$message), "ERROR") + plot(1, type="n", axes=FALSE, xlab="", ylab="") + text(1, 1, "Error creating priority map", cex=1.5) +}) +``` + +\pagebreak +## Crop Anomaly Detection + +The map below highlights potential anomalies in your fields that may require investigation. Areas with sudden changes in CI values could indicate weeding activities, rapid weed growth, or other management interventions. + +**How it works:** This map compares current week's CI values with those from the previous week: +- **Significant Growth** (dark green): CI increase > 2 units (potential weed growth) +- **Moderate Growth** (light green): CI increase of 1-2 units +- **Stable** (yellow): CI change between -1 and +1 units +- **Moderate Decline** (orange): CI decrease of 1-2 units +- **Significant Decline** (red): CI decrease > 2 units (potential weeding/harvesting activities) + +Areas with significant growth (dark green) may indicate rapid weed growth that requires monitoring, while significant declines (red) often indicate recent management activities like weeding or harvesting. + +```{r render_anomaly_map, echo=FALSE, fig.height=6, fig.width=9, message=FALSE, warning=FALSE} +# Create anomaly detection map +tryCatch({ + # Create and display the anomaly map + anomaly_map <- create_anomaly_map( + ci_current = CI, + ci_previous = CI_m1, + field_boundaries = AllPivots0 + ) + + # Print the map + print(anomaly_map) +}, error = function(e) { + safe_log(paste("Error creating anomaly map:", e$message), "ERROR") + plot(1, type="n", axes=FALSE, xlab="", ylab="") + text(1, 1, "Error creating anomaly map", cex=1.5) +}) +``` + +\pagebreak +## Harvest Planning + +This map shows the harvest readiness status of all fields, helping you plan harvest operations and logistics. Fields in dark green are ready for harvest, while those in yellow through red are at earlier growth stages. + +**How it works:** Harvest readiness is determined by combining field age and CI values: +- **Ready for harvest** (dark green): Fields ≥52 weeks old with CI ≥4.0 +- **Approaching harvest** (light green): Fields ≥48 weeks old with CI ≥3.5 +- **Mid-maturity** (yellow): Fields ≥40 weeks old with CI ≥3.0 +- **Growing** (orange): Fields ≥12 weeks old +- **Early stage** (red): Fields <12 weeks old + +This classification helps prioritize harvesting operations and logistical planning by identifying fields that are at optimal maturity for maximum sugar content. + +```{r render_harvest_map, echo=FALSE, fig.height=6, fig.width=9, message=FALSE, warning=FALSE} +# Create harvest planning map +tryCatch({ + # Create and display the harvest readiness map + harvest_map <- create_field_status_map( + field_boundaries = AllPivots0, + health_data = farm_health_data, + attribute = "harvest_readiness", + title = "Harvest Readiness Status", + palette = "RdYlGn", + legend_title = "Harvest Status" + ) + + # Print the map + print(harvest_map) +}, error = function(e) { + safe_log(paste("Error creating harvest map:", e$message), "ERROR") + plot(1, type="n", axes=FALSE, xlab="", ylab="") + text(1, 1, "Error creating harvest map", cex=1.5) +}) +``` + +\pagebreak +## Field Status Summary + +The charts below provide an overview of your farm's health and harvest readiness status, showing the distribution of fields across different health categories and maturity stages. + +**How the Field Status Chart works:** This bar chart displays the count of fields in each health status category, based on the same CI thresholds described in the Farm Health Status section: +- **Excellent** (dark green): CI ≥ 5.0 +- **Good** (light green): CI 3.5-4.99 +- **Fair** (yellow): CI 2.0-3.49 +- **Poor** (orange): CI 1.0-1.99 +- **Critical** (red): CI < 1.0 + +**How the Harvest Readiness Chart works:** This pie chart shows the distribution of fields by harvest readiness, allowing you to see at a glance how many fields are in each stage of development. Fields are categorized based on both age and CI values as described in the Harvest Planning section above. + +```{r render_status_charts, echo=FALSE, fig.height=5, fig.width=10, message=FALSE, warning=FALSE} +# Create field status summary visualization +tryCatch({ + # Create field status charts + status_chart <- create_summary_stats(farm_health_data) + + # Print the chart + print(status_chart) + + # Create a second row with harvest readiness chart + harvest_chart <- create_harvest_readiness_chart(farm_health_data) + + # Print the chart + print(harvest_chart) +}, error = function(e) { + safe_log(paste("Error creating status summary charts:", e$message), "ERROR") + plot(1, type="n", axes=FALSE, xlab="", ylab="") + text(1, 1, "Error creating status summary charts", cex=1.5) +}) +``` + +## Priority Fields Requiring Attention + +The chart below highlights fields that require immediate management attention based on their health scores and anomaly detection. These should be prioritized for field inspections. + +**How it works:** This chart shows fields with priority levels 1-3 (critical, high, and medium): +- Fields are ordered by priority level, with the most critical fields on the left +- Bar height represents the Chlorophyll Index (CI) value +- Bar colors indicate priority level: red (critical), orange (high), yellow (medium) +- Text labels show the detected anomaly type for each field + +The table below the chart provides detailed metrics for these priority fields, including CI values, weekly changes, anomaly types, and harvest status. Only fields requiring management attention (priority levels 1-3) are included. + +```{r render_priority_fields_chart, echo=FALSE, fig.height=5, fig.width=10, message=FALSE, warning=FALSE} +# Create priority fields chart +tryCatch({ + # Create and display priority fields chart + priority_chart <- create_priority_fields_chart(farm_health_data) + + # Print the chart + print(priority_chart) + + # Create a table of priority fields + priority_table <- farm_health_data %>% + dplyr::filter(priority_level <= 3) %>% + dplyr::arrange(priority_level, field) %>% + dplyr::select( + Field = field, + Status = status, + `CI Value` = mean_ci, + `Weekly Change` = ci_change, + `Anomaly Type` = anomaly_type, + `Age (Weeks)` = age_weeks, + `Harvest Status` = harvest_readiness + ) + + # Display the table if there are priority fields + if (nrow(priority_table) > 0) { + knitr::kable(priority_table, caption = "Priority Fields Requiring Management Attention") + } else { + cat("No priority fields requiring immediate attention this week.") + } + +}, error = function(e) { + safe_log(paste("Error creating priority fields chart:", e$message), "ERROR") + cat("Error generating priority fields visualization. See log for details.") +}) +``` + +\pagebreak +## Management Recommendations + +```{r render_recommendations, echo=FALSE, results='asis', message=FALSE, warning=FALSE} +# Generate executive recommendations +tryCatch({ + # Create and display recommendations + recommendations_html <- generate_executive_recommendations(farm_health_data) + + # Print the HTML recommendations + cat(recommendations_html) +}, error = function(e) { + safe_log(paste("Error creating recommendations:", e$message), "ERROR") + cat("

Error generating recommendations. Please see system administrator.

") +}) +``` + +## Yield Prediction Overview + +This section provides yield predictions for mature fields (over 300 days old) based on their Chlorophyll Index values and growth patterns. These predictions can help with harvest planning and yield forecasting. + +```{r render_yield_summary, echo=FALSE, fig.height=5, fig.width=10, message=FALSE, warning=FALSE} +# Create yield summary +tryCatch({ + if (exists("pred_rf_current_season") && nrow(pred_rf_current_season) > 0) { + # Calculate total estimated production + total_yield <- sum(pred_rf_current_season$predicted_Tcha, na.rm = TRUE) + + # Create summary box + cat("
") + cat("

Yield Summary

") + cat("") + cat("
") + + # Display yield prediction table + harvest_ready_fields <- pred_rf_current_season %>% + dplyr::arrange(desc(predicted_Tcha)) %>% + dplyr::select( + Field = field, + `Sub Field` = sub_field, + `Age (Days)` = Age_days, + `Cumulative CI` = total_CI, + `Predicted Yield (Tonnes/ha)` = predicted_Tcha + ) + + knitr::kable(harvest_ready_fields, + caption = "Predicted Yields for Harvest-Ready Fields", + digits = 1) + } else { + cat("
") + cat("

Yield Summary

") + cat("

No fields currently meet harvest readiness criteria (>300 days) for yield prediction.

") + cat("
") + } +}, error = function(e) { + safe_log(paste("Error creating yield summary:", e$message), "ERROR") + cat("

Error generating yield summary. Please see system administrator.

") +}) +``` + diff --git a/r_app/ci_extraction.R b/r_app/ci_extraction.R index dfb8889..ed2f447 100644 --- a/r_app/ci_extraction.R +++ b/r_app/ci_extraction.R @@ -61,8 +61,19 @@ main <- function() { # 3. Initialize project configuration # -------------------------------- new_project_question <- FALSE - source("r_app/parameters_project.R") - source("r_app/ci_extraction_utils.R") + + tryCatch({ + source("parameters_project.R") + source("ci_extraction_utils.R") + }, error = function(e) { + warning("Default source files not found. Attempting to source from 'r_app' directory.") + tryCatch({ + source("r_app/parameters_project.R") + source("r_app/ci_extraction_utils.R") + }, error = function(e) { + stop("Failed to source required files from both default and 'r_app' directories.") + }) + }) # 4. Generate date list for processing # --------------------------------- diff --git a/r_app/interpolate_growth_model.R b/r_app/interpolate_growth_model.R index c34db80..d55dd1e 100644 --- a/r_app/interpolate_growth_model.R +++ b/r_app/interpolate_growth_model.R @@ -33,8 +33,18 @@ main <- function() { } # Initialize project configuration and load utility functions - source("r_app/parameters_project.R") - source("r_app/growth_model_utils.R") + tryCatch({ + source("parameters_project.R") + source("ci_extraction_utils.R") + }, error = function(e) { + warning("Default source files not found. Attempting to source from 'r_app' directory.") + tryCatch({ + source("r_app/parameters_project.R") + source("r_app/ci_extraction_utils.R") + }, error = function(e) { + stop("Failed to source required files from both default and 'r_app' directories.") + }) + }) log_message("Starting CI growth model interpolation") diff --git a/r_app/mosaic_creation.R b/r_app/mosaic_creation.R index cc98263..ed20ba2 100644 --- a/r_app/mosaic_creation.R +++ b/r_app/mosaic_creation.R @@ -67,8 +67,18 @@ main <- function() { # 3. Initialize project configuration # -------------------------------- - source("r_app/parameters_project.R") - source("r_app/mosaic_creation_utils.R") + tryCatch({ + source("parameters_project.R") + source("ci_extraction_utils.R") + }, error = function(e) { + warning("Default source files not found. Attempting to source from 'r_app' directory.") + tryCatch({ + source("r_app/parameters_project.R") + source("r_app/ci_extraction_utils.R") + }, error = function(e) { + stop("Failed to source required files from both default and 'r_app' directories.") + }) + }) # 4. Generate date range for processing # --------------------------------- diff --git a/r_app/packages.R b/r_app/packages.R index 51bc3ec..d6534fa 100644 --- a/r_app/packages.R +++ b/r_app/packages.R @@ -56,7 +56,7 @@ load_smartcane_packages <- function(verbose = FALSE) { "sf", # Simple Features for spatial vector data "terra", # Raster data processing "exactextractr", # Fast extraction from rasters - "raster", # Legacy raster package (for compatibility) + "tmap", # Thematic mapping for spatial visualization # Data manipulation "tidyverse", # Collection of data manipulation packages @@ -64,12 +64,19 @@ load_smartcane_packages <- function(verbose = FALSE) { "readxl", # Excel file reading "stringr", # String manipulation "purrr", # Functional programming tools + "zoo", # Time series processing with rolling functions # Visualization "ggplot2", # Advanced plotting "leaflet", # Interactive maps "plotly", # Interactive plots + # Machine learning and statistics + "caret", # Classification and regression training + "rsample", # Data sampling for modeling + "randomForest", # Random forest implementation + "CAST", # Feature selection for spatial data + # Project management "here", # Path handling diff --git a/r_app/tests/test_report_utils.R b/r_app/tests/test_report_utils.R new file mode 100644 index 0000000..611f7db --- /dev/null +++ b/r_app/tests/test_report_utils.R @@ -0,0 +1,280 @@ +# test_report_utils.R +# +# Tests for visualization functions in report_utils.R +# + +# Load the test framework +source("tests/test_framework.R") + +# Set up test environment +env <- setup_test_env() + +# Required libraries for testing +library(testthat) +library(terra) +library(sf) +library(dplyr) +library(ggplot2) + +# Load the functions to test +source("../report_utils.R") + +# Create mock data for testing +create_mock_data <- function() { + # Create a simple raster for testing + r <- terra::rast(nrows=10, ncols=10, xmin=0, xmax=10, ymin=0, ymax=10, vals=1:100) + names(r) <- "CI" + + # Create a simple field boundary + field_boundaries <- sf::st_sf( + field = c("Field1", "Field2"), + sub_field = c("A", "B"), + geometry = sf::st_sfc( + sf::st_polygon(list(rbind(c(1,1), c(5,1), c(5,5), c(1,5), c(1,1)))), + sf::st_polygon(list(rbind(c(6,6), c(9,6), c(9,9), c(6,9), c(6,6)))) + ), + crs = sf::st_crs(r) + ) + + # Create mock harvest data + harvesting_data <- data.frame( + field = c("Field1", "Field2"), + sub_field = c("A", "B"), + age = c(100, 150), + season_start = as.Date(c("2023-01-01", "2023-02-01")), + year = c(2023, 2023) + ) + + # Create mock CI quadrant data + ci_quadrant <- data.frame( + field = rep(c("Field1", "Field2"), each=10), + sub_field = rep(c("A", "B"), each=10), + Date = rep(seq(as.Date("2023-01-01"), by="week", length.out=10), 2), + DOY = rep(1:10, 2), + cumulative_CI = rep(cumsum(1:10), 2), + value = rep(1:10, 2), + season = rep(2023, 20), + model = rep(c("northwest", "northeast", "southwest", "southeast"), 5) + ) + + return(list( + raster = r, + field_boundaries = field_boundaries, + harvesting_data = harvesting_data, + ci_quadrant = ci_quadrant + )) +} + +# Test the create_CI_map function +test_that("create_CI_map creates a valid tmap object", { + mock_data <- create_mock_data() + + # Test with all required parameters + map <- create_CI_map( + pivot_raster = mock_data$raster, + pivot_shape = mock_data$field_boundaries[1,], + pivot_spans = mock_data$field_boundaries[1,], + week = "01", + age = 10, + borders = TRUE, + use_breaks = TRUE + ) + + # Check if it returned a tmap object + expect_true("tmap" %in% class(map)) + + # Test with missing parameters + expect_error(create_CI_map(pivot_shape = mock_data$field_boundaries[1,], + pivot_spans = mock_data$field_boundaries[1,], + week = "01", age = 10), + "pivot_raster is required") + + expect_error(create_CI_map(pivot_raster = mock_data$raster, + pivot_spans = mock_data$field_boundaries[1,], + week = "01", age = 10), + "pivot_shape is required") +}) + +# Test the create_CI_diff_map function +test_that("create_CI_diff_map creates a valid tmap object", { + mock_data <- create_mock_data() + + # Test with all required parameters + map <- create_CI_diff_map( + pivot_raster = mock_data$raster, + pivot_shape = mock_data$field_boundaries[1,], + pivot_spans = mock_data$field_boundaries[1,], + week_1 = "01", + week_2 = "02", + age = 10, + borders = TRUE, + use_breaks = TRUE + ) + + # Check if it returned a tmap object + expect_true("tmap" %in% class(map)) + + # Test with missing parameters + expect_error(create_CI_diff_map(pivot_shape = mock_data$field_boundaries[1,], + pivot_spans = mock_data$field_boundaries[1,], + week_1 = "01", week_2 = "02", age = 10), + "pivot_raster is required") + + expect_error(create_CI_diff_map(pivot_raster = mock_data$raster, + pivot_spans = mock_data$field_boundaries[1,], + week_1 = "01", age = 10), + "week_1 and week_2 parameters are required") +}) + +# Test the ci_plot function +test_that("ci_plot handles input parameters correctly", { + mock_data <- create_mock_data() + + # Capture output to avoid cluttering the test output + temp_file <- tempfile() + sink(temp_file) + + # Test with all required parameters - should not throw an error + expect_error( + ci_plot( + pivotName = "Field1", + field_boundaries = mock_data$field_boundaries, + current_ci = mock_data$raster, + ci_minus_1 = mock_data$raster, + ci_minus_2 = mock_data$raster, + last_week_diff = mock_data$raster, + three_week_diff = mock_data$raster, + harvesting_data = mock_data$harvesting_data, + week = "01", + week_minus_1 = "52", + week_minus_2 = "51", + week_minus_3 = "50", + use_breaks = TRUE, + borders = TRUE + ), + NA # Expect no error + ) + + # Test with missing parameters + expect_error( + ci_plot(), + "pivotName is required" + ) + + # Test with invalid field name + expect_error( + ci_plot( + pivotName = "NonExistentField", + field_boundaries = mock_data$field_boundaries, + current_ci = mock_data$raster, + ci_minus_1 = mock_data$raster, + ci_minus_2 = mock_data$raster, + last_week_diff = mock_data$raster, + three_week_diff = mock_data$raster, + harvesting_data = mock_data$harvesting_data + ), + regexp = NULL # We expect some error related to the field not being found + ) + + # Reset output + sink() + unlink(temp_file) +}) + +# Test the cum_ci_plot function +test_that("cum_ci_plot handles input parameters correctly", { + mock_data <- create_mock_data() + + # Capture output to avoid cluttering the test output + temp_file <- tempfile() + sink(temp_file) + + # Test with all required parameters - should not throw an error + expect_error( + cum_ci_plot( + pivotName = "Field1", + ci_quadrant_data = mock_data$ci_quadrant, + plot_type = "value", + facet_on = FALSE, + x_unit = "days" + ), + NA # Expect no error + ) + + # Test with different plot types + expect_error( + cum_ci_plot( + pivotName = "Field1", + ci_quadrant_data = mock_data$ci_quadrant, + plot_type = "CI_rate" + ), + NA # Expect no error + ) + + expect_error( + cum_ci_plot( + pivotName = "Field1", + ci_quadrant_data = mock_data$ci_quadrant, + plot_type = "cumulative_CI" + ), + NA # Expect no error + ) + + # Test with invalid plot type + expect_error( + cum_ci_plot( + pivotName = "Field1", + ci_quadrant_data = mock_data$ci_quadrant, + plot_type = "invalid_type" + ), + "plot_type must be one of: 'value', 'CI_rate', or 'cumulative_CI'" + ) + + # Test with missing parameters + expect_error( + cum_ci_plot(), + "pivotName is required" + ) + + # Reset output + sink() + unlink(temp_file) +}) + +# Test the get_week_path function +test_that("get_week_path returns correct path", { + # Test with valid inputs + path <- get_week_path( + mosaic_path = "ci_max_mosaics", + input_date = "2023-01-15", + week_offset = 0 + ) + + # Extract the week number and year from the path + expect_match(path, "week_02_2023\\.tif$", all = FALSE) # Week 2 of 2023 + + # Test with offset + path_minus_1 <- get_week_path( + mosaic_path = "ci_max_mosaics", + input_date = "2023-01-15", + week_offset = -1 + ) + expect_match(path_minus_1, "week_01_2023\\.tif$", all = FALSE) + + # Test with missing parameters + expect_error( + get_week_path(input_date = "2023-01-15", week_offset = 0), + "mosaic_path is required" + ) + + expect_error( + get_week_path(mosaic_path = "ci_max_mosaics", week_offset = 0), + "input_date is required" + ) +}) + +# Clean up +teardown_test_env() + +# Print success message +cat("Report utility function tests completed successfully\n") \ No newline at end of file