--- 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.") }) ```