Timon
4cd62ab82e
- Updated `create_CI_map` and `create_CI_diff_map` functions to enforce a 1:1 aspect ratio for consistent map sizing. - Modified `ci_plot` function to adjust widths of arranged maps for better layout. - Changed raster merging method in `aggregate_per_field_mosaics_to_farm_level` from `mosaic` to `merge` for improved handling of field data. - Introduced `test_kpi_validation.R` script to validate the structure of KPI RDS files, ensuring expected KPIs are present. - Added `test_overview_maps_aggregation.R` script to test the aggregation pipeline for overview maps, including loading field mosaics, creating a farm-level mosaic, and generating visualizations.
828 lines
29 KiB
R
828 lines
29 KiB
R
# 80_UTILS_AGRONOMIC_SUPPORT.R
|
|
# ============================================================================
|
|
# SPECIFIC KPI UTILITIES (SCRIPT 80 - CLIENT TYPE: agronomic_support)
|
|
#
|
|
# Contains all 6 KPI calculation functions and helpers:
|
|
# - Field uniformity KPI (CV-based, spatial autocorrelation)
|
|
# - Area change KPI (week-over-week CI changes)
|
|
# - TCH forecasted KPI (tonnage projections from harvest data)
|
|
# - Growth decline KPI (trend analysis)
|
|
# - Weed presence KPI (field fragmentation detection)
|
|
# - Gap filling KPI (interpolation quality)
|
|
# - KPI reporting (summary tables, field details, text interpretation)
|
|
# - KPI export (Excel, RDS, data export)
|
|
#
|
|
# Orchestrator: calculate_all_kpis()
|
|
# Dependencies: 00_common_utils.R (safe_log), sourced from common
|
|
# Used by: 80_calculate_kpis.R (when client_type == "agronomic_support")
|
|
# ============================================================================
|
|
|
|
library(terra)
|
|
library(sf)
|
|
library(dplyr)
|
|
library(tidyr)
|
|
library(readxl)
|
|
library(writexl)
|
|
library(spdep)
|
|
|
|
# ============================================================================
|
|
# SHARED HELPER FUNCTIONS (NOW IN 80_UTILS_COMMON.R)
|
|
# ============================================================================
|
|
# The following helper functions have been moved to 80_utils_common.R:
|
|
# - calculate_cv()
|
|
# - calculate_change_percentages()
|
|
# - calculate_spatial_autocorrelation()
|
|
# - extract_ci_values()
|
|
# - calculate_week_numbers()
|
|
# - load_field_ci_raster()
|
|
# - load_weekly_ci_mosaic()
|
|
# - prepare_predictions()
|
|
#
|
|
# These are now sourced from common utils and shared by all client types.
|
|
# ============================================================================
|
|
|
|
#' Prepare harvest predictions and ensure proper alignment with field data
|
|
prepare_predictions <- function(harvest_model, field_data, scenario = "optimistic") {
|
|
if (is.null(harvest_model) || is.null(field_data)) {
|
|
return(NULL)
|
|
}
|
|
|
|
tryCatch({
|
|
scenario_factor <- switch(scenario,
|
|
"pessimistic" = 0.85,
|
|
"realistic" = 1.0,
|
|
"optimistic" = 1.15,
|
|
1.0)
|
|
|
|
predictions <- field_data %>%
|
|
mutate(tch_forecasted = field_data$mean_ci * scenario_factor)
|
|
|
|
return(predictions)
|
|
}, error = function(e) {
|
|
message(paste("Error preparing predictions:", e$message))
|
|
return(NULL)
|
|
})
|
|
}
|
|
|
|
# ============================================================================
|
|
# KPI CALCULATION FUNCTIONS (6 KPIS)
|
|
# ============================================================================
|
|
|
|
#' KPI 1: Calculate field uniformity based on CV and spatial autocorrelation
|
|
#'
|
|
#' Measures how uniform crop development is across the field.
|
|
#' Low CV + high positive Moran's I = excellent uniformity
|
|
#'
|
|
#' @param ci_pixels_by_field List of CI pixel arrays for each field
|
|
#' @param field_boundaries_sf SF object with field geometries
|
|
#' @param ci_raster Raster object with CI values (for spatial autocorrelation)
|
|
#'
|
|
#' @return Data frame with field_idx, cv_value, morans_i, uniformity_score, interpretation
|
|
calculate_field_uniformity_kpi <- function(ci_pixels_by_field, field_boundaries_sf, ci_raster = NULL) {
|
|
results_list <- list()
|
|
|
|
for (field_idx in seq_len(nrow(field_boundaries_sf))) {
|
|
ci_pixels <- ci_pixels_by_field[[field_idx]]
|
|
|
|
if (is.null(ci_pixels) || length(ci_pixels) == 0) {
|
|
results_list[[length(results_list) + 1]] <- list(
|
|
field_idx = field_idx,
|
|
cv_value = NA_real_,
|
|
morans_i = NA_real_,
|
|
uniformity_score = NA_real_,
|
|
interpretation = "No data"
|
|
)
|
|
next
|
|
}
|
|
|
|
cv_val <- calculate_cv(ci_pixels)
|
|
|
|
morans_i <- NA_real_
|
|
if (!is.null(ci_raster)) {
|
|
morans_result <- calculate_spatial_autocorrelation(ci_raster, field_boundaries_sf[field_idx, ])
|
|
if (is.list(morans_result)) {
|
|
morans_i <- morans_result$morans_i
|
|
} else {
|
|
morans_i <- morans_result
|
|
}
|
|
}
|
|
|
|
# Normalize CV (0-1 scale, invert so lower CV = higher score)
|
|
cv_normalized <- min(cv_val / 0.3, 1) # 0.3 = threshold for CV
|
|
cv_score <- 1 - cv_normalized
|
|
|
|
# Normalize Moran's I (-1 to 1 scale, shift to 0-1)
|
|
morans_normalized <- if (!is.na(morans_i)) {
|
|
(morans_i + 1) / 2
|
|
} else {
|
|
0.5
|
|
}
|
|
|
|
uniformity_score <- 0.7 * cv_score + 0.3 * morans_normalized
|
|
|
|
# Interpretation
|
|
if (is.na(cv_val)) {
|
|
interpretation <- "No data"
|
|
} else if (cv_val < 0.08) {
|
|
interpretation <- "Excellent uniformity"
|
|
} else if (cv_val < 0.15) {
|
|
interpretation <- "Good uniformity"
|
|
} else if (cv_val < 0.25) {
|
|
interpretation <- "Acceptable uniformity"
|
|
} else if (cv_val < 0.4) {
|
|
interpretation <- "Poor uniformity"
|
|
} else {
|
|
interpretation <- "Very poor uniformity"
|
|
}
|
|
|
|
results_list[[length(results_list) + 1]] <- list(
|
|
field_idx = field_idx,
|
|
cv_value = cv_val,
|
|
morans_i = morans_i,
|
|
uniformity_score = round(uniformity_score, 3),
|
|
interpretation = interpretation
|
|
)
|
|
}
|
|
|
|
# Convert accumulated list to data frame in a single operation
|
|
result <- do.call(rbind, lapply(results_list, as.data.frame))
|
|
|
|
return(result)
|
|
}
|
|
|
|
#' KPI 2: Calculate area change metric (week-over-week CI changes)
|
|
#'
|
|
#' Tracks the percentage change in CI between current and previous week
|
|
#'
|
|
#' @param current_stats Current week field statistics (from extract_field_statistics_from_ci)
|
|
#' @param previous_stats Previous week field statistics
|
|
#'
|
|
#' @return Data frame with field-level CI changes
|
|
calculate_area_change_kpi <- function(current_stats, previous_stats) {
|
|
result <- calculate_change_percentages(current_stats, previous_stats)
|
|
|
|
# Add interpretation
|
|
result$interpretation <- NA_character_
|
|
|
|
for (i in seq_len(nrow(result))) {
|
|
change <- result$mean_ci_pct_change[i]
|
|
|
|
if (is.na(change)) {
|
|
result$interpretation[i] <- "No previous data"
|
|
} else if (change > 15) {
|
|
result$interpretation[i] <- "Rapid growth"
|
|
} else if (change > 5) {
|
|
result$interpretation[i] <- "Positive growth"
|
|
} else if (change > -5) {
|
|
result$interpretation[i] <- "Stable"
|
|
} else if (change > -15) {
|
|
result$interpretation[i] <- "Declining"
|
|
} else {
|
|
result$interpretation[i] <- "Rapid decline"
|
|
}
|
|
}
|
|
|
|
return(result)
|
|
}
|
|
|
|
#' KPI 3: Calculate TCH forecasted (tonnes of cane per hectare)
|
|
#'
|
|
#' Projects final harvest tonnage based on CI growth trajectory
|
|
#'
|
|
#' @param field_statistics Current field statistics
|
|
#' @param harvesting_data Historical harvest data (with yield observations)
|
|
#' @param field_boundaries_sf Field geometries
|
|
#'
|
|
#' @return Data frame with field-level TCH forecasts
|
|
calculate_tch_forecasted_kpi <- function(field_statistics, harvesting_data = NULL, field_boundaries_sf = NULL) {
|
|
result <- data.frame(
|
|
field_idx = field_statistics$field_idx,
|
|
mean_ci = field_statistics$mean_ci,
|
|
tch_forecasted = NA_real_,
|
|
tch_lower_bound = NA_real_,
|
|
tch_upper_bound = NA_real_,
|
|
confidence = NA_character_,
|
|
stringsAsFactors = FALSE
|
|
)
|
|
|
|
# Base TCH model: TCH = 50 + (CI * 10)
|
|
# This is a simplified model; production use should include more variables
|
|
|
|
for (i in seq_len(nrow(result))) {
|
|
if (is.na(result$mean_ci[i])) {
|
|
result$confidence[i] <- "No data"
|
|
next
|
|
}
|
|
|
|
if (is.na(result$mean_ci[i])) {
|
|
result$tch_forecasted[i] <- NA_real_
|
|
result$tch_lower_bound[i] <- NA_real_
|
|
result$tch_upper_bound[i] <- NA_real_
|
|
result$confidence[i] <- "No data"
|
|
}
|
|
}
|
|
|
|
return(result)
|
|
}
|
|
|
|
#' KPI 4: Calculate growth decline indicator
|
|
#'
|
|
#' Identifies fields with negative growth trajectory
|
|
#'
|
|
#' @param ci_values_list List of CI values for each field (multiple weeks)
|
|
#'
|
|
#' @return Data frame with field-level decline indicators
|
|
calculate_growth_decline_kpi <- function(ci_values_list) {
|
|
result <- data.frame(
|
|
field_idx = seq_len(length(ci_values_list)),
|
|
four_week_trend = NA_real_,
|
|
trend_interpretation = NA_character_,
|
|
decline_severity = NA_character_,
|
|
stringsAsFactors = FALSE
|
|
)
|
|
|
|
for (field_idx in seq_len(length(ci_values_list))) {
|
|
ci_vals <- ci_values_list[[field_idx]]
|
|
|
|
if (is.null(ci_vals) || length(ci_vals) < 2) {
|
|
result$trend_interpretation[field_idx] <- "Insufficient data"
|
|
next
|
|
}
|
|
|
|
ci_vals <- ci_vals[!is.na(ci_vals)]
|
|
if (length(ci_vals) < 2) {
|
|
result$trend_interpretation[field_idx] <- "Insufficient data"
|
|
next
|
|
}
|
|
|
|
# Calculate linear trend
|
|
weeks <- seq_along(ci_vals)
|
|
lm_fit <- lm(ci_vals ~ weeks)
|
|
slope <- coef(lm_fit)["weeks"]
|
|
|
|
result$four_week_trend[field_idx] <- round(as.numeric(slope), 3)
|
|
|
|
if (slope > 0.1) {
|
|
result$trend_interpretation[field_idx] <- "Strong growth"
|
|
result$decline_severity[field_idx] <- "None"
|
|
} else if (slope > 0) {
|
|
result$trend_interpretation[field_idx] <- "Weak growth"
|
|
result$decline_severity[field_idx] <- "None"
|
|
} else if (slope > -0.1) {
|
|
result$trend_interpretation[field_idx] <- "Slight decline"
|
|
result$decline_severity[field_idx] <- "Low"
|
|
} else if (slope > -0.3) {
|
|
result$trend_interpretation[field_idx] <- "Moderate decline"
|
|
result$decline_severity[field_idx] <- "Medium"
|
|
} else {
|
|
result$trend_interpretation[field_idx] <- "Strong decline"
|
|
result$decline_severity[field_idx] <- "High"
|
|
}
|
|
}
|
|
|
|
return(result)
|
|
}
|
|
|
|
#' KPI 5: Calculate weed presence indicator
|
|
#'
|
|
#' Detects field fragmentation/patchiness (potential weed/pest pressure)
|
|
#'
|
|
#' @param ci_pixels_by_field List of CI pixel arrays for each field
|
|
#'
|
|
#' @return Data frame with fragmentation indicators
|
|
calculate_weed_presence_kpi <- function(ci_pixels_by_field) {
|
|
result <- data.frame(
|
|
field_idx = seq_len(length(ci_pixels_by_field)),
|
|
cv_value = NA_real_,
|
|
low_ci_percent = NA_real_,
|
|
fragmentation_index = NA_real_,
|
|
weed_pressure_risk = NA_character_,
|
|
stringsAsFactors = FALSE
|
|
)
|
|
|
|
for (field_idx in seq_len(length(ci_pixels_by_field))) {
|
|
ci_pixels <- ci_pixels_by_field[[field_idx]]
|
|
|
|
if (is.null(ci_pixels) || length(ci_pixels) == 0) {
|
|
result$weed_pressure_risk[field_idx] <- "No data"
|
|
next
|
|
}
|
|
|
|
ci_pixels <- ci_pixels[!is.na(ci_pixels)]
|
|
if (length(ci_pixels) == 0) {
|
|
result$weed_pressure_risk[field_idx] <- "No data"
|
|
next
|
|
}
|
|
|
|
cv_val <- calculate_cv(ci_pixels)
|
|
low_ci_pct <- sum(ci_pixels < 1.5) / length(ci_pixels) * 100
|
|
fragmentation <- cv_val * low_ci_pct / 100
|
|
|
|
result$cv_value[field_idx] <- cv_val
|
|
result$low_ci_percent[field_idx] <- round(low_ci_pct, 2)
|
|
result$fragmentation_index[field_idx] <- round(fragmentation, 3)
|
|
|
|
if (is.na(fragmentation)) {
|
|
result$weed_pressure_risk[field_idx] <- "No data"
|
|
} else if (fragmentation > 0.15) {
|
|
result$weed_pressure_risk[field_idx] <- "High"
|
|
} else if (fragmentation > 0.08) {
|
|
result$weed_pressure_risk[field_idx] <- "Medium"
|
|
} else if (fragmentation > 0.04) {
|
|
result$weed_pressure_risk[field_idx] <- "Low"
|
|
} else {
|
|
result$weed_pressure_risk[field_idx] <- "Minimal"
|
|
}
|
|
}
|
|
|
|
return(result)
|
|
}
|
|
|
|
#' Calculate Gap Filling Score KPI (placeholder)
|
|
#' @param ci_raster Current week CI raster
|
|
#' @param field_boundaries Field boundaries
|
|
#' @return Data frame with field-level gap filling scores
|
|
calculate_gap_filling_kpi <- function(ci_raster, field_boundaries) {
|
|
# Handle both sf and SpatVector inputs
|
|
if (!inherits(field_boundaries, "SpatVector")) {
|
|
field_boundaries_vect <- terra::vect(field_boundaries)
|
|
} else {
|
|
field_boundaries_vect <- field_boundaries
|
|
}
|
|
|
|
results_list <- list()
|
|
|
|
# Ensure field_boundaries_vect is valid and matches field_boundaries dimensions
|
|
n_fields_vect <- length(field_boundaries_vect)
|
|
n_fields_sf <- nrow(field_boundaries)
|
|
|
|
if (n_fields_sf != n_fields_vect) {
|
|
warning(paste("Field boundary mismatch: nrow(field_boundaries)=", n_fields_sf, "vs length(field_boundaries_vect)=", n_fields_vect, ". Using actual SpatVector length."))
|
|
}
|
|
|
|
for (i in seq_len(n_fields_vect)) {
|
|
field_vect <- field_boundaries_vect[i]
|
|
|
|
# Extract CI values using helper function
|
|
ci_values <- extract_ci_values(ci_raster, field_vect)
|
|
valid_values <- ci_values[!is.na(ci_values) & is.finite(ci_values)]
|
|
|
|
if (length(valid_values) > 1) {
|
|
# Calculate % of valid (non-NA) values = gap filling success
|
|
total_pixels <- length(ci_values)
|
|
valid_pixels <- length(valid_values)
|
|
gap_filling_success <- (valid_pixels / total_pixels) * 100
|
|
na_percent <- ((total_pixels - valid_pixels) / total_pixels) * 100
|
|
|
|
results_list[[length(results_list) + 1]] <- list(
|
|
field_idx = i,
|
|
gap_filling_success = round(gap_filling_success, 2),
|
|
na_percent_pre_interpolation = round(na_percent, 2),
|
|
mean_ci = round(mean(valid_values), 2)
|
|
)
|
|
} else {
|
|
# Not enough valid data
|
|
results_list[[length(results_list) + 1]] <- list(
|
|
field_idx = i,
|
|
gap_filling_success = NA_real_,
|
|
na_percent_pre_interpolation = NA_real_,
|
|
mean_ci = NA_real_
|
|
)
|
|
}
|
|
}
|
|
|
|
# Convert accumulated list to data frame in a single operation
|
|
field_results <- do.call(rbind, lapply(results_list, as.data.frame))
|
|
|
|
return(field_results)
|
|
}
|
|
|
|
# ============================================================================
|
|
# KPI ORCHESTRATOR AND REPORTING
|
|
# ============================================================================
|
|
|
|
#' Create summary tables for all 6 KPIs (AGGREGATED farm-level summaries)
|
|
#'
|
|
#' @param all_kpis List containing results from all 6 KPI functions (per-field data)
|
|
#'
|
|
#' @return List of summary data frames ready for reporting (farm-level aggregates)
|
|
create_summary_tables <- function(all_kpis) {
|
|
|
|
# ==========================================
|
|
# 1. UNIFORMITY SUMMARY (count by interpretation)
|
|
# ==========================================
|
|
uniformity_summary <- all_kpis$uniformity %>%
|
|
group_by(interpretation) %>%
|
|
summarise(
|
|
field_count = n(),
|
|
avg_cv = mean(cv_value, na.rm = TRUE),
|
|
avg_morans_i = mean(morans_i, na.rm = TRUE),
|
|
.groups = 'drop'
|
|
) %>%
|
|
rename(
|
|
Status = interpretation,
|
|
`Field Count` = field_count,
|
|
`Avg CV` = avg_cv,
|
|
`Avg Moran's I` = avg_morans_i
|
|
)
|
|
|
|
# ==========================================
|
|
# 2. AREA CHANGE SUMMARY (improving/stable/declining counts)
|
|
# ==========================================
|
|
area_change_summary <- all_kpis$area_change %>%
|
|
group_by(interpretation) %>%
|
|
summarise(
|
|
field_count = n(),
|
|
avg_ci_change = mean(mean_ci_pct_change, na.rm = TRUE),
|
|
.groups = 'drop'
|
|
) %>%
|
|
rename(
|
|
Status = interpretation,
|
|
`Field Count` = field_count,
|
|
`Avg CI Change %` = avg_ci_change
|
|
)
|
|
|
|
# ==========================================
|
|
# 3. TCH FORECAST SUMMARY (yield statistics)
|
|
# ==========================================
|
|
tch_summary <- all_kpis$tch_forecasted %>%
|
|
summarise(
|
|
avg_tch = mean(tch_forecasted, na.rm = TRUE),
|
|
min_tch = min(tch_forecasted, na.rm = TRUE),
|
|
max_tch = max(tch_forecasted, na.rm = TRUE),
|
|
avg_ci = mean(mean_ci, na.rm = TRUE),
|
|
fields_with_data = sum(!is.na(tch_forecasted))
|
|
) %>%
|
|
rename(
|
|
`Avg Forecast (t/ha)` = avg_tch,
|
|
`Min (t/ha)` = min_tch,
|
|
`Max (t/ha)` = max_tch,
|
|
`Avg CI` = avg_ci,
|
|
`Fields` = fields_with_data
|
|
)
|
|
|
|
# ==========================================
|
|
# 4. GROWTH DECLINE SUMMARY (trend interpretation)
|
|
# ==========================================
|
|
growth_summary <- all_kpis$growth_decline %>%
|
|
group_by(trend_interpretation) %>%
|
|
summarise(
|
|
field_count = n(),
|
|
avg_trend = mean(four_week_trend, na.rm = TRUE),
|
|
.groups = 'drop'
|
|
) %>%
|
|
rename(
|
|
Trend = trend_interpretation,
|
|
`Field Count` = field_count,
|
|
`Avg 4-Week Trend` = avg_trend
|
|
)
|
|
|
|
# ==========================================
|
|
# 5. WEED PRESSURE SUMMARY (risk level counts)
|
|
# ==========================================
|
|
weed_summary <- all_kpis$weed_presence %>%
|
|
group_by(weed_pressure_risk) %>%
|
|
summarise(
|
|
field_count = n(),
|
|
avg_fragmentation = mean(fragmentation_index, na.rm = TRUE),
|
|
.groups = 'drop'
|
|
) %>%
|
|
rename(
|
|
`Risk Level` = weed_pressure_risk,
|
|
`Field Count` = field_count,
|
|
`Avg Fragmentation` = avg_fragmentation
|
|
)
|
|
|
|
# ==========================================
|
|
# 6. GAP FILLING SUMMARY
|
|
# ==========================================
|
|
gap_summary <- if (!is.null(all_kpis$gap_filling) && is.data.frame(all_kpis$gap_filling) && nrow(all_kpis$gap_filling) > 0) {
|
|
all_kpis$gap_filling %>%
|
|
summarise(
|
|
avg_gap_filling = mean(gap_filling_success, na.rm = TRUE),
|
|
avg_na_percent = mean(na_percent_pre_interpolation, na.rm = TRUE),
|
|
fields_with_data = n()
|
|
) %>%
|
|
rename(
|
|
`Avg Gap Filling Success %` = avg_gap_filling,
|
|
`Avg NA % Pre-Interpolation` = avg_na_percent,
|
|
`Fields Analyzed` = fields_with_data
|
|
)
|
|
} else {
|
|
data.frame(`Avg Gap Filling Success %` = NA_real_, `Avg NA % Pre-Interpolation` = NA_real_, `Fields Analyzed` = 0, check.names = FALSE)
|
|
}
|
|
|
|
# Return as list (each element is a farm-level summary table)
|
|
kpi_summary <- list(
|
|
uniformity = uniformity_summary,
|
|
area_change = area_change_summary,
|
|
tch_forecast = tch_summary,
|
|
growth_decline = growth_summary,
|
|
weed_pressure = weed_summary,
|
|
gap_filling = gap_summary
|
|
)
|
|
|
|
return(kpi_summary)
|
|
}
|
|
|
|
#' Create detailed field-by-field KPI report
|
|
#'
|
|
#' @param field_df Data frame with field identifiers and acreage
|
|
#' @param all_kpis List with all KPI results
|
|
#' @param field_boundaries_sf SF object with field boundaries
|
|
#'
|
|
#' @return Data frame with one row per field, all KPI columns (renamed for reporting compatibility)
|
|
create_field_detail_table <- function(field_df, all_kpis, field_boundaries_sf) {
|
|
result <- field_df %>%
|
|
left_join(
|
|
all_kpis$uniformity %>% select(field_idx, cv_value, uniformity_interpretation = interpretation),
|
|
by = c("field_idx")
|
|
) %>%
|
|
left_join(
|
|
all_kpis$area_change %>% select(field_idx, mean_ci_pct_change),
|
|
by = c("field_idx")
|
|
) %>%
|
|
left_join(
|
|
all_kpis$tch_forecasted %>% select(field_idx, tch_forecasted, mean_ci),
|
|
by = c("field_idx")
|
|
) %>%
|
|
left_join(
|
|
all_kpis$growth_decline %>% select(field_idx, decline_severity),
|
|
by = c("field_idx")
|
|
) %>%
|
|
left_join(
|
|
all_kpis$weed_presence %>% select(field_idx, weed_pressure_risk),
|
|
by = c("field_idx")
|
|
) %>%
|
|
# Rename columns to match reporting script expectations
|
|
rename(
|
|
Field = field_name,
|
|
`Growth Uniformity` = uniformity_interpretation,
|
|
`Yield Forecast (t/ha)` = tch_forecasted,
|
|
`Decline Risk` = decline_severity,
|
|
`Weed Risk` = weed_pressure_risk,
|
|
`CI Change %` = mean_ci_pct_change,
|
|
`Mean CI` = mean_ci,
|
|
`CV Value` = cv_value
|
|
) %>%
|
|
# Add placeholder columns expected by reporting script (will be populated from other sources)
|
|
mutate(
|
|
`Field Size (ha)` = NA_real_,
|
|
`Gap Score` = NA_real_
|
|
) %>%
|
|
select(field_idx, Field, `Field Size (ha)`, `Growth Uniformity`, `Yield Forecast (t/ha)`,
|
|
`Gap Score`, `Decline Risk`, `Weed Risk`, `CI Change %`, `Mean CI`, `CV Value`)
|
|
|
|
return(result)
|
|
}
|
|
|
|
#' Generate KPI text interpretation for inclusion in Word report
|
|
#'
|
|
#' @param all_kpis List with all KPI results
|
|
#'
|
|
#' @return Character string with formatted KPI summary text
|
|
create_field_kpi_text <- function(all_kpis) {
|
|
text_parts <- c(
|
|
"## KPI ANALYSIS SUMMARY\n",
|
|
"### Field Uniformity\n",
|
|
paste(all_kpis$uniformity$interpretation, collapse = "; "), "\n",
|
|
"### Growth Trends\n",
|
|
paste(all_kpis$growth_decline$trend_interpretation, collapse = "; "), "\n",
|
|
"### Weed/Pest Pressure\n",
|
|
paste(all_kpis$weed_presence$weed_pressure_risk, collapse = "; "), "\n"
|
|
)
|
|
|
|
return(paste(text_parts, collapse = ""))
|
|
}
|
|
|
|
#' Export detailed KPI data to Excel/RDS
|
|
#'
|
|
#' @param all_kpis List with all KPI results (per-field data)
|
|
#' @param kpi_summary List with summary tables (farm-level aggregates)
|
|
#' @param project_dir Project name (for filename)
|
|
#' @param output_dir Directory for output files
|
|
#' @param week Week number
|
|
#' @param year Year
|
|
#' @param field_boundaries_sf SF object with field boundaries (optional, for field_details_table)
|
|
#'
|
|
#' @return List of output file paths
|
|
export_kpi_data <- function(all_kpis, kpi_summary, project_dir, output_dir, week, year, field_boundaries_sf = NULL) {
|
|
# Ensure output directory exists
|
|
if (!dir.exists(output_dir)) {
|
|
dir.create(output_dir, recursive = TRUE)
|
|
}
|
|
|
|
# Create unified field details table if field_boundaries_sf is provided
|
|
field_details_table <- NULL
|
|
if (!is.null(field_boundaries_sf)) {
|
|
tryCatch({
|
|
# Create a basic field_df from the boundaries
|
|
# Robust field name extraction with multiple fallbacks
|
|
field_name <- NA_character_
|
|
|
|
# Check for 'name' column in the data.frame
|
|
if ("name" %in% names(field_boundaries_sf)) {
|
|
field_name <- field_boundaries_sf$name
|
|
} else if ("properties" %in% names(field_boundaries_sf)) {
|
|
# Extract from properties column (may be a list-column)
|
|
props <- field_boundaries_sf$properties
|
|
if (is.list(props) && length(props) > 0 && "name" %in% names(props[[1]])) {
|
|
field_name <- sapply(props, function(x) ifelse(is.null(x$name), NA_character_, x$name))
|
|
} else if (!is.list(props)) {
|
|
# Try direct access if properties is a simple column
|
|
field_name <- props
|
|
}
|
|
}
|
|
|
|
# Ensure field_name is a character vector of appropriate length
|
|
if (length(field_name) != nrow(field_boundaries_sf)) {
|
|
field_name <- rep(NA_character_, nrow(field_boundaries_sf))
|
|
}
|
|
|
|
# Replace only NA elements with fallback names, keeping valid names intact
|
|
na_indices <- which(is.na(field_name))
|
|
if (length(na_indices) > 0) {
|
|
field_name[na_indices] <- paste0("Field_", na_indices)
|
|
}
|
|
|
|
field_df <- data.frame(
|
|
field_idx = 1:nrow(field_boundaries_sf),
|
|
field_name = field_name,
|
|
stringsAsFactors = FALSE
|
|
)
|
|
|
|
field_details_table <- create_field_detail_table(field_df, all_kpis, field_boundaries_sf)
|
|
message(paste("✓ Field details table created with", nrow(field_details_table), "fields"))
|
|
}, error = function(e) {
|
|
message(paste("WARNING: Could not create field_details_table:", e$message))
|
|
})
|
|
}
|
|
|
|
# Export all KPI tables to a single Excel file - use project_dir"
|
|
excel_file <- file.path(output_dir, paste0(project_dir, "_kpi_summary_tables_week", sprintf("%02d_%d", week, year), ".xlsx"))
|
|
|
|
sheets <- list(
|
|
"Uniformity" = as.data.frame(kpi_summary$uniformity),
|
|
"Area_Change" = as.data.frame(kpi_summary$area_change),
|
|
"TCH_Forecast" = as.data.frame(kpi_summary$tch_forecast),
|
|
"Growth_Decline" = as.data.frame(kpi_summary$growth_decline),
|
|
"Weed_Pressure" = as.data.frame(kpi_summary$weed_pressure),
|
|
"Gap_Filling" = as.data.frame(kpi_summary$gap_filling)
|
|
)
|
|
|
|
write_xlsx(sheets, excel_file)
|
|
message(paste("✓ KPI data exported to:", excel_file))
|
|
|
|
# Export to RDS for programmatic access (CRITICAL: Both per-field AND summary tables)
|
|
# The reporting script expects: summary_tables (list of 6 summary tables)
|
|
# We also provide: all_kpis (per-field data) and field_details (unified field view)
|
|
rds_file <- file.path(output_dir, paste0(project_dir, "_kpi_summary_tables_week", sprintf("%02d_%d", week, year), ".rds"))
|
|
|
|
# Create the export structure that reporting scripts expect
|
|
export_data <- list(
|
|
summary_tables = kpi_summary, # Farm-level aggregates (6 KPI summaries)
|
|
all_kpis = all_kpis, # Per-field data (6 KPI per-field tables)
|
|
field_details = field_details_table # Unified field-level detail table
|
|
)
|
|
|
|
saveRDS(export_data, rds_file)
|
|
message(paste("✓ KPI RDS exported to:", rds_file))
|
|
message(" Structure: list($summary_tables, $all_kpis, $field_details)")
|
|
|
|
# Return including field_details for orchestrator to capture
|
|
return(list(
|
|
excel = excel_file,
|
|
rds = rds_file,
|
|
field_details = field_details_table
|
|
))
|
|
}
|
|
|
|
# ============================================================================
|
|
# ORCHESTRATOR FUNCTION
|
|
# ============================================================================
|
|
|
|
#' Calculate all 6 KPIs
|
|
#'
|
|
#' Main entry point for KPI calculation.
|
|
#' This function orchestrates the 6 KPI calculations and returns all results.
|
|
#'
|
|
#' @param field_boundaries_sf SF object with field geometries
|
|
#' @param current_week ISO week number (1-53)
|
|
#' @param current_year ISO week year
|
|
#' @param current_mosaic_dir Directory containing current week's mosaic
|
|
#' @param previous_mosaic_dir Directory containing previous week's mosaic (optional)
|
|
#' @param ci_rds_path Path to combined CI RDS file
|
|
#' @param harvesting_data Data frame with harvest data (optional)
|
|
#' @param output_dir Directory for KPI exports
|
|
#' @param project_dir Project name (for filename in exports)
|
|
#'
|
|
#' @return List with results from all 6 KPI functions
|
|
#'
|
|
#' @details
|
|
#' This function:
|
|
#' 1. Loads current week mosaic and extracts field statistics
|
|
#' 2. (Optionally) loads previous week mosaic for comparison metrics
|
|
#' 3. Calculates all 6 KPIs
|
|
#' 4. Creates summary tables
|
|
#' 5. Exports results to Excel/RDS
|
|
#'
|
|
calculate_all_kpis <- function(
|
|
field_boundaries_sf,
|
|
current_week,
|
|
current_year,
|
|
current_mosaic_dir,
|
|
previous_mosaic_dir = NULL,
|
|
ci_rds_path = NULL,
|
|
harvesting_data = NULL,
|
|
output_dir = NULL,
|
|
project_dir = NULL
|
|
) {
|
|
|
|
message("\n============ KPI CALCULATION (6 KPIs) ============")
|
|
|
|
# Load current week mosaic
|
|
message("Loading current week mosaic...")
|
|
current_mosaic <- load_weekly_ci_mosaic(current_week, current_year, current_mosaic_dir)
|
|
|
|
if (is.null(current_mosaic)) {
|
|
stop("Could not load current week mosaic")
|
|
}
|
|
|
|
# Extract field statistics
|
|
message("Extracting field statistics from current mosaic...")
|
|
current_stats <- extract_field_statistics_from_ci(current_mosaic, field_boundaries_sf)
|
|
ci_pixels_by_field <- extract_ci_values(current_mosaic, field_boundaries_sf)
|
|
|
|
# Load previous week mosaic (if available)
|
|
previous_stats <- NULL
|
|
if (!is.null(previous_mosaic_dir)) {
|
|
target_prev <- calculate_target_week_and_year(current_week, current_year, offset_weeks = 1)
|
|
message(paste("Loading previous week mosaic (week", target_prev$week, target_prev$year, ")..."))
|
|
previous_mosaic <- load_weekly_ci_mosaic(target_prev$week, target_prev$year, previous_mosaic_dir)
|
|
|
|
if (!is.null(previous_mosaic)) {
|
|
previous_stats <- extract_field_statistics_from_ci(previous_mosaic, field_boundaries_sf)
|
|
} else {
|
|
message("Previous week mosaic not available - skipping area change KPI")
|
|
}
|
|
}
|
|
|
|
# Calculate 6 KPIs
|
|
message("\nCalculating KPI 1: Field Uniformity...")
|
|
uniformity_kpi <- calculate_field_uniformity_kpi(ci_pixels_by_field, field_boundaries_sf, current_mosaic)
|
|
|
|
message("Calculating KPI 2: Area Change...")
|
|
if (!is.null(previous_stats)) {
|
|
area_change_kpi <- calculate_area_change_kpi(current_stats, previous_stats)
|
|
} else {
|
|
area_change_kpi <- data.frame(
|
|
field_idx = seq_len(nrow(field_boundaries_sf)),
|
|
mean_ci_pct_change = NA_real_,
|
|
interpretation = rep("No previous data", nrow(field_boundaries_sf))
|
|
)
|
|
}
|
|
|
|
message("Calculating KPI 3: TCH Forecasted...")
|
|
tch_kpi <- calculate_tch_forecasted_kpi(current_stats, harvesting_data, field_boundaries_sf)
|
|
|
|
message("Calculating KPI 4: Growth Decline...")
|
|
growth_decline_kpi <- calculate_growth_decline_kpi(
|
|
ci_pixels_by_field # Would need historical data for real trend
|
|
)
|
|
|
|
message("Calculating KPI 5: Weed Presence...")
|
|
weed_kpi <- calculate_weed_presence_kpi(ci_pixels_by_field)
|
|
|
|
message("Calculating KPI 6: Gap Filling...")
|
|
gap_filling_kpi <- calculate_gap_filling_kpi(current_mosaic, field_boundaries_sf)
|
|
|
|
# Compile results
|
|
all_kpis <- list(
|
|
uniformity = uniformity_kpi,
|
|
area_change = area_change_kpi,
|
|
tch_forecasted = tch_kpi,
|
|
growth_decline = growth_decline_kpi,
|
|
weed_presence = weed_kpi,
|
|
gap_filling = gap_filling_kpi
|
|
)
|
|
|
|
# Create summary tables
|
|
kpi_summary <- create_summary_tables(all_kpis)
|
|
|
|
# Export - pass project_dir for proper filename and field_boundaries_sf for field details table
|
|
if (is.null(project_dir)) {
|
|
project_dir <- "AURA" # Fallback if not provided
|
|
}
|
|
export_result <- export_kpi_data(all_kpis, kpi_summary, project_dir, output_dir, current_week, current_year, field_boundaries_sf)
|
|
|
|
message(paste("\n✓", project_dir, "KPI calculation complete. Week", current_week, current_year, "\n"))
|
|
|
|
# Return combined structure (for integration with 80_calculate_kpis.R)
|
|
# Capture field_details from export_result to propagate it out
|
|
return(list(
|
|
all_kpis = all_kpis,
|
|
summary_tables = kpi_summary,
|
|
field_details = export_result$field_details # Propagate field_details from export_kpi_data
|
|
))
|
|
}
|