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