1033 lines
36 KiB
R
1033 lines
36 KiB
R
# 80_UTILS_AGRONOMIC_SUPPORT.R
|
|
# ============================================================================
|
|
# AURA-SPECIFIC KPI UTILITIES (SCRIPT 80 - CLIENT TYPE: agronomic_support)
|
|
#
|
|
# Contains all 6 AURA 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_field_analysis_agronomic_support()
|
|
# 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)
|
|
library(caret)
|
|
library(CAST)
|
|
|
|
# ============================================================================
|
|
# 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.
|
|
# ============================================================================
|
|
|
|
# ============================================================================
|
|
# AURA 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_band Raster band with CI values
|
|
#'
|
|
#' @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_band = NULL,
|
|
mosaic_dir = NULL, week_file = NULL) {
|
|
result <- data.frame(
|
|
field_idx = integer(),
|
|
cv_value = numeric(),
|
|
morans_i = numeric(),
|
|
uniformity_score = numeric(),
|
|
uniformity_category = character(),
|
|
interpretation = character(),
|
|
stringsAsFactors = FALSE
|
|
)
|
|
|
|
# Determine if we're using per-field structure
|
|
is_per_field <- !is.null(mosaic_dir) && !is.null(week_file)
|
|
|
|
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) {
|
|
result <- rbind(result, data.frame(
|
|
field_idx = field_idx,
|
|
cv_value = NA_real_,
|
|
morans_i = NA_real_,
|
|
uniformity_score = NA_real_,
|
|
uniformity_category = "No data",
|
|
interpretation = "No data",
|
|
stringsAsFactors = FALSE
|
|
))
|
|
next
|
|
}
|
|
|
|
cv_val <- calculate_cv(ci_pixels)
|
|
|
|
# Calculate Moran's I
|
|
morans_i <- NA_real_
|
|
if (is_per_field) {
|
|
# Load individual field raster for per-field structure
|
|
field_name <- field_boundaries_sf$field[field_idx]
|
|
field_mosaic_path <- file.path(mosaic_dir, field_name, week_file)
|
|
|
|
if (file.exists(field_mosaic_path)) {
|
|
tryCatch({
|
|
field_raster <- terra::rast(field_mosaic_path)[["CI"]]
|
|
single_field <- field_boundaries_sf[field_idx, ]
|
|
morans_result <- calculate_spatial_autocorrelation(field_raster, single_field)
|
|
|
|
if (is.list(morans_result)) {
|
|
morans_i <- morans_result$morans_i
|
|
} else {
|
|
morans_i <- morans_result
|
|
}
|
|
}, error = function(e) {
|
|
message(paste(" Warning: Spatial autocorrelation failed for field", field_name, ":", e$message))
|
|
})
|
|
}
|
|
} else if (!is.null(ci_band) && inherits(ci_band, "SpatRaster")) {
|
|
# Use single raster for single-file structure
|
|
tryCatch({
|
|
single_field <- field_boundaries_sf[field_idx, ]
|
|
morans_result <- calculate_spatial_autocorrelation(ci_band, single_field)
|
|
|
|
if (is.list(morans_result)) {
|
|
morans_i <- morans_result$morans_i
|
|
} else {
|
|
morans_i <- morans_result
|
|
}
|
|
}, error = function(e) {
|
|
message(paste(" Warning: Spatial autocorrelation failed for field", field_idx, ":", e$message))
|
|
})
|
|
}
|
|
|
|
# Normalize CV (0-1 scale, invert so lower CV = higher score)
|
|
cv_normalized <- min(cv_val / 0.3, 1)
|
|
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"
|
|
uniformity_category <- "No data"
|
|
} else if (cv_val < 0.08) {
|
|
interpretation <- "Excellent uniformity"
|
|
uniformity_category <- "Excellent"
|
|
} else if (cv_val < 0.15) {
|
|
interpretation <- "Good uniformity"
|
|
uniformity_category <- "Good"
|
|
} else if (cv_val < 0.25) {
|
|
interpretation <- "Acceptable uniformity"
|
|
uniformity_category <- "Acceptable"
|
|
} else if (cv_val < 0.4) {
|
|
interpretation <- "Poor uniformity"
|
|
uniformity_category <- "Poor"
|
|
} else {
|
|
interpretation <- "Very poor uniformity"
|
|
uniformity_category <- "Very poor"
|
|
}
|
|
|
|
result <- rbind(result, data.frame(
|
|
field_idx = field_idx,
|
|
cv_value = cv_val,
|
|
morans_i = morans_i,
|
|
uniformity_score = round(uniformity_score, 3),
|
|
uniformity_category = uniformity_category,
|
|
interpretation = interpretation,
|
|
stringsAsFactors = FALSE
|
|
))
|
|
}
|
|
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) {
|
|
|
|
# Initialize result data frame
|
|
result <- data.frame(
|
|
field_idx = seq_len(nrow(current_stats)),
|
|
mean_ci_pct_change = NA_real_,
|
|
interpretation = NA_character_,
|
|
stringsAsFactors = FALSE
|
|
)
|
|
|
|
# Handle case where previous stats is NULL or empty
|
|
if (is.null(previous_stats) || nrow(previous_stats) == 0) {
|
|
result$interpretation <- "No previous data"
|
|
return(result)
|
|
}
|
|
|
|
# Match fields between current and previous stats
|
|
# Handle both naming conventions (Field_id vs field_idx)
|
|
if ("Field_id" %in% names(current_stats)) {
|
|
current_field_col <- "Field_id"
|
|
prev_field_col <- "Field_id"
|
|
ci_col <- "Mean_CI"
|
|
} else {
|
|
current_field_col <- "field_idx"
|
|
prev_field_col <- "field_idx"
|
|
ci_col <- "mean_ci"
|
|
}
|
|
|
|
# Create lookup for previous stats
|
|
prev_lookup <- setNames(
|
|
previous_stats[[ci_col]],
|
|
previous_stats[[prev_field_col]]
|
|
)
|
|
|
|
# Calculate percentage change for each field
|
|
for (i in seq_len(nrow(current_stats))) {
|
|
current_field_id <- current_stats[[current_field_col]][i]
|
|
current_ci <- current_stats[[ci_col]][i]
|
|
|
|
# Find matching previous CI value
|
|
prev_ci <- prev_lookup[[as.character(current_field_id)]]
|
|
|
|
if (!is.null(prev_ci) && !is.na(prev_ci) && !is.na(current_ci) && prev_ci > 0) {
|
|
# Calculate percentage change
|
|
pct_change <- ((current_ci - prev_ci) / prev_ci) * 100
|
|
result$mean_ci_pct_change[i] <- round(pct_change, 2)
|
|
|
|
# Add interpretation
|
|
if (pct_change > 15) {
|
|
result$interpretation[i] <- "Rapid growth"
|
|
} else if (pct_change > 5) {
|
|
result$interpretation[i] <- "Positive growth"
|
|
} else if (pct_change > -5) {
|
|
result$interpretation[i] <- "Stable"
|
|
} else if (pct_change > -15) {
|
|
result$interpretation[i] <- "Declining"
|
|
} else {
|
|
result$interpretation[i] <- "Rapid decline"
|
|
}
|
|
} else {
|
|
result$interpretation[i] <- "No previous data"
|
|
}
|
|
}
|
|
|
|
return(result)
|
|
}
|
|
|
|
#' KPI 3: Calculate TCH forecasted (tonnes of cane per hectare)
|
|
#'
|
|
#' Projects final harvest tonnage based on historical yield data and CI growth trajectory.
|
|
#' Uses a Random Forest model trained on harvest data to predict yields for mature fields.
|
|
#' Delegates to calculate_yield_prediction_kpi() in 80_utils_common.R.
|
|
#'
|
|
#' @param field_statistics Current field statistics (dataframe with Mean_CI or mean_ci column)
|
|
#' @param harvesting_data Historical harvest data frame (with tonnage_ha column)
|
|
#' @param field_boundaries_sf SF object with field geometries
|
|
#' @param cumulative_CI_vals_dir Directory with combined CI RDS files (optional)
|
|
#' @param data_dir Project data directory (from setup_project_directories or parameters_project.R)
|
|
#' Used to build cumulative_CI_vals_dir path if not provided directly (optional)
|
|
#' @param project_dir Deprecated: only used if data_dir not provided (optional)
|
|
#'
|
|
#' @return Data frame with field-level yield forecasts ready for orchestrator
|
|
#' Columns: field_idx, tch_forecasted (yields in t/ha)
|
|
calculate_tch_forecasted_kpi <- function(field_statistics, harvesting_data = NULL,
|
|
field_boundaries_sf = NULL,
|
|
cumulative_CI_vals_dir = NULL,
|
|
data_dir = NULL,
|
|
project_dir = NULL) {
|
|
|
|
# Use common utils yield prediction function (handles all ML logic)
|
|
# This replaces the previous linear model (TCH = 50 + CI*10) with proper ML prediction
|
|
|
|
# Validate required parameters
|
|
if (is.null(field_boundaries_sf)) {
|
|
safe_log("field_boundaries_sf is NULL in calculate_tch_forecasted_kpi", "WARNING")
|
|
return(data.frame(
|
|
field_idx = integer(),
|
|
tch_forecasted = numeric(),
|
|
stringsAsFactors = FALSE
|
|
))
|
|
}
|
|
|
|
# Determine cumulative CI directory
|
|
if (is.null(cumulative_CI_vals_dir)) {
|
|
# Priority 1: Use provided data_dir parameter
|
|
if (!is.null(data_dir)) {
|
|
cumulative_CI_vals_dir <- file.path(data_dir, "extracted_ci", "cumulative_vals")
|
|
} else if (exists("data_dir", envir = .GlobalEnv)) {
|
|
# Priority 2: Fallback to global data_dir from parameters_project.R
|
|
cumulative_CI_vals_dir <- file.path(get("data_dir", envir = .GlobalEnv), "extracted_ci", "cumulative_vals")
|
|
} else {
|
|
# Priority 3: Last resort - log warning and fail gracefully
|
|
safe_log("Missing project data directory configuration: provide data_dir parameter or ensure parameters_project.R has set data_dir globally", "WARNING")
|
|
safe_log("No training data available for yield prediction", "WARNING")
|
|
return(data.frame(
|
|
field_idx = integer(),
|
|
tch_forecasted = numeric(),
|
|
stringsAsFactors = FALSE
|
|
))
|
|
}
|
|
}
|
|
|
|
# Call the shared yield prediction function from common utils
|
|
yield_result <- calculate_yield_prediction_kpi(field_boundaries_sf, harvesting_data, cumulative_CI_vals_dir)
|
|
|
|
# Extract field-level results from the list
|
|
field_results <- yield_result$field_results
|
|
|
|
# Convert to format expected by orchestrator
|
|
# If no predictions, return empty data frame
|
|
if (is.null(field_results) || nrow(field_results) == 0) {
|
|
return(data.frame(
|
|
field_idx = integer(),
|
|
tch_forecasted = numeric(),
|
|
stringsAsFactors = FALSE
|
|
))
|
|
}
|
|
|
|
# Map field names to field_idx using field_boundaries_sf
|
|
result <- field_results %>%
|
|
mutate(
|
|
field_idx = match(field, field_boundaries_sf$field),
|
|
tch_forecasted = yield_forecast_t_ha
|
|
) %>%
|
|
filter(!is.na(field_idx)) %>%
|
|
select(field_idx, tch_forecasted)
|
|
|
|
# Ensure result has proper structure even if empty
|
|
if (nrow(result) == 0) {
|
|
return(data.frame(
|
|
field_idx = integer(),
|
|
tch_forecasted = numeric(),
|
|
stringsAsFactors = FALSE
|
|
))
|
|
}
|
|
|
|
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)
|
|
}
|
|
|
|
#' KPI 5: Calculate field patchiness (combines fragmentation into Gini-like metric + risk)
|
|
#'
|
|
#' @param weed_kpi_result Data frame from calculate_weed_presence_kpi()
|
|
#' @param mean_ci_values Optional vector of mean CI values per field
|
|
#'
|
|
#' @return Data frame with patchiness indicators (gini_coefficient, patchiness_risk, interpretation)
|
|
calculate_patchiness_kpi <- function(weed_kpi_result, ci_pixels_by_field = NULL, mean_ci_values = NULL) {
|
|
result <- weed_kpi_result %>%
|
|
mutate(
|
|
# Calculate Gini coefficient from CI pixels (proper calculation)
|
|
gini_coefficient = NA_real_,
|
|
mean_ci = if (!is.null(mean_ci_values)) mean_ci_values[field_idx] else NA_real_,
|
|
# Map weed_pressure_risk to patchiness_risk
|
|
patchiness_risk = weed_pressure_risk,
|
|
# Create interpretation based on gini and risk
|
|
patchiness_interpretation = case_when(
|
|
is.na(gini_coefficient) ~ "No data",
|
|
gini_coefficient < 0.2 & patchiness_risk %in% c("Low", "Minimal") ~ "Uniform growth",
|
|
gini_coefficient < 0.4 & patchiness_risk %in% c("Low", "Medium") ~ "Moderate patchiness",
|
|
gini_coefficient >= 0.4 & patchiness_risk %in% c("High", "Medium") ~ "High patchiness",
|
|
TRUE ~ "Mixed heterogeneity"
|
|
)
|
|
) %>%
|
|
select(field_idx, gini_coefficient, mean_ci,
|
|
patchiness_interpretation, patchiness_risk)
|
|
|
|
# Calculate actual Gini coefficients if CI pixels provided
|
|
if (!is.null(ci_pixels_by_field)) {
|
|
for (i in seq_len(nrow(result))) {
|
|
ci_pixels <- ci_pixels_by_field[[i]]
|
|
|
|
if (!is.null(ci_pixels) && length(ci_pixels) > 0) {
|
|
ci_pixels <- ci_pixels[!is.na(ci_pixels)]
|
|
|
|
if (length(ci_pixels) > 1) {
|
|
# Calculate Gini coefficient
|
|
# Formula: Gini = (2 * sum(i * x_i)) / (n * sum(x_i)) - (n+1)/n
|
|
# where x_i are sorted values
|
|
ci_sorted <- sort(ci_pixels)
|
|
n <- length(ci_sorted)
|
|
|
|
numerator <- 2 * sum(seq_len(n) * ci_sorted)
|
|
denominator <- n * sum(ci_sorted)
|
|
|
|
gini <- (numerator / denominator) - (n + 1) / n
|
|
# Clamp to 0-1 range (should be within this by formula but guard against numerical errors)
|
|
gini <- max(0, min(1, gini))
|
|
|
|
result$gini_coefficient[i] <- gini
|
|
|
|
# Update interpretation based on calculated Gini
|
|
result$patchiness_interpretation[i] <- dplyr::case_when(
|
|
gini < 0.2 ~ "Uniform growth",
|
|
gini < 0.4 & result$patchiness_risk[i] %in% c("Low", "Medium", "Minimal") ~ "Moderate patchiness",
|
|
gini >= 0.4 & result$patchiness_risk[i] %in% c("High", "Medium") ~ "High patchiness",
|
|
TRUE ~ "Mixed heterogeneity"
|
|
)
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return(result)
|
|
}
|
|
|
|
|
|
|
|
# ============================================================================
|
|
# KPI ORCHESTRATOR AND REPORTING
|
|
# ============================================================================
|
|
|
|
#' Create summary tables for all 6 KPIs
|
|
#'
|
|
#' @param all_kpis List containing results from all 6 KPI functions
|
|
#'
|
|
#' @return List of summary data frames ready for reporting
|
|
create_summary_tables <- function(all_kpis) {
|
|
kpi_summary <- list(
|
|
uniformity = all_kpis$uniformity %>%
|
|
select(field_idx, cv_value, uniformity_category, interpretation),
|
|
|
|
spatial_clustering = if (!is.null(all_kpis$spatial_clustering) && nrow(all_kpis$spatial_clustering) > 0) {
|
|
all_kpis$spatial_clustering %>%
|
|
select(field_idx, morans_i)
|
|
} else {
|
|
NULL
|
|
},
|
|
|
|
area_change = all_kpis$area_change %>%
|
|
select(field_idx, mean_ci_pct_change, interpretation),
|
|
|
|
tch_forecast = all_kpis$tch_forecasted %>%
|
|
select(field_idx, tch_forecasted),
|
|
|
|
growth_decline = all_kpis$growth_decline %>%
|
|
select(field_idx, four_week_trend, trend_interpretation, decline_severity),
|
|
|
|
patchiness = all_kpis$patchiness %>%
|
|
select(field_idx, gini_coefficient, patchiness_interpretation, patchiness_risk),
|
|
|
|
gap_filling = if (!is.null(all_kpis$gap_filling) && nrow(all_kpis$gap_filling) > 0) {
|
|
all_kpis$gap_filling %>%
|
|
select(field_idx, gap_score, gap_level)
|
|
} else {
|
|
NULL
|
|
}
|
|
)
|
|
return(kpi_summary)
|
|
}
|
|
|
|
#' Create detailed field-by-field KPI report (ALL KPIs in one row)
|
|
#'
|
|
#' @param field_boundaries_sf SF object with field boundaries
|
|
#' @param all_kpis List with all KPI results
|
|
#' @param current_week Current week number
|
|
#' @param current_year Current year
|
|
#'
|
|
#' @return Data frame with one row per field, all KPI columns
|
|
create_field_detail_table <- function(field_boundaries_sf, all_kpis, current_week, current_year) {
|
|
|
|
# Start with field identifiers AND field_idx for joining
|
|
result <- field_boundaries_sf %>%
|
|
sf::st_drop_geometry() %>%
|
|
mutate(
|
|
field_idx = row_number(),
|
|
Field_id = field,
|
|
Field_name = field,
|
|
Week = current_week,
|
|
Year = current_year
|
|
) %>%
|
|
select(field_idx, Field_id, Field_name, Week, Year)
|
|
|
|
# ============================================
|
|
# GROUP 1: FIELD UNIFORMITY (KPI 1)
|
|
# ============================================
|
|
result <- result %>%
|
|
left_join(
|
|
all_kpis$uniformity %>%
|
|
select(field_idx, CV = cv_value,
|
|
Uniformity_Category = uniformity_category),
|
|
by = "field_idx"
|
|
)
|
|
|
|
# ============================================
|
|
# GROUP 2: GROWTH & TREND ANALYSIS (KPI 2 + KPI 4)
|
|
# ============================================
|
|
# KPI 2: Area Change
|
|
result <- result %>%
|
|
left_join(
|
|
all_kpis$area_change %>%
|
|
select(field_idx, Weekly_CI_Change = mean_ci_pct_change,
|
|
Area_Change_Interpretation = interpretation),
|
|
by = "field_idx"
|
|
)
|
|
|
|
# KPI 4: Growth Decline
|
|
result <- result %>%
|
|
left_join(
|
|
all_kpis$growth_decline %>%
|
|
select(field_idx, Four_Week_Trend = four_week_trend,
|
|
Trend_Interpretation = trend_interpretation,
|
|
Decline_Severity = decline_severity),
|
|
by = "field_idx"
|
|
)
|
|
|
|
# ============================================
|
|
# GROUP 3: FIELD HETEROGENEITY/PATCHINESS (KPI 5 + Spatial Clustering)
|
|
# ============================================
|
|
# KPI 5: Field Patchiness
|
|
result <- result %>%
|
|
left_join(
|
|
all_kpis$patchiness %>%
|
|
select(field_idx, Gini_Coefficient = gini_coefficient, Mean_CI = mean_ci,
|
|
Patchiness_Interpretation = patchiness_interpretation,
|
|
Patchiness_Risk = patchiness_risk),
|
|
by = "field_idx"
|
|
)
|
|
|
|
# Moran's I (spatial clustering - used in patchiness calculation)
|
|
if (!is.null(all_kpis$spatial_clustering) && nrow(all_kpis$spatial_clustering) > 0) {
|
|
result <- result %>%
|
|
left_join(
|
|
all_kpis$spatial_clustering %>%
|
|
select(field_idx, Morans_I = morans_i),
|
|
by = "field_idx"
|
|
)
|
|
}
|
|
|
|
# ============================================
|
|
# GROUP 4: YIELD FORECAST (KPI 3)
|
|
# ============================================
|
|
result <- result %>%
|
|
left_join(
|
|
all_kpis$tch_forecasted %>%
|
|
select(field_idx, TCH_Forecasted = tch_forecasted),
|
|
by = "field_idx"
|
|
)
|
|
|
|
# ============================================
|
|
# GROUP 5: DATA QUALITY / GAP FILLING (KPI 6)
|
|
# ============================================
|
|
# Add gap filling if available
|
|
if (!is.null(all_kpis$gap_filling) && nrow(all_kpis$gap_filling) > 0) {
|
|
result <- result %>%
|
|
left_join(
|
|
all_kpis$gap_filling %>%
|
|
select(field_idx, Gap_Score = gap_score, Gap_Level = gap_level),
|
|
by = "field_idx"
|
|
)
|
|
}
|
|
|
|
# Remove field_idx from final output
|
|
result <- result %>%
|
|
select(-field_idx)
|
|
|
|
# Round numeric columns
|
|
result <- result %>%
|
|
mutate(across(where(is.numeric), ~ round(., 2)))
|
|
|
|
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(
|
|
"## AURA 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 field_detail_df Data frame with all KPI columns (one row per field)
|
|
#' @param kpi_summary List with summary tables (optional, for metadata)
|
|
#' @param output_dir Directory for output files
|
|
#' @param week Week number
|
|
#' @param year Year
|
|
#' @param project_dir Project name
|
|
#' @return List of output file paths
|
|
export_kpi_data <- function(field_detail_df, kpi_summary, output_dir, week, year, project_dir) {
|
|
|
|
# Use the common export function from 80_utils_common.R
|
|
export_paths <- export_field_analysis_excel(
|
|
field_df = field_detail_df,
|
|
summary_df = NULL, # No separate summary sheet for agronomic support
|
|
project_dir = project_dir,
|
|
current_week = week,
|
|
year = year,
|
|
reports_dir = output_dir
|
|
)
|
|
|
|
return(export_paths)
|
|
}
|
|
|
|
# ============================================================================
|
|
# ORCHESTRATOR FUNCTION
|
|
# ============================================================================
|
|
|
|
#' Calculate all 6 AURA KPIs
|
|
#'
|
|
#' Main entry point for AURA 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
|
|
#'
|
|
#' @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 AURA KPIs
|
|
#' 4. Creates summary tables
|
|
#' 5. Exports results to Excel/RDS
|
|
#'
|
|
calculate_all_field_analysis_agronomic_support <- 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,
|
|
data_dir = NULL,
|
|
project_dir = NULL
|
|
) {
|
|
|
|
message("\n============ AURA KPI CALCULATION (6 KPIs) ============")
|
|
|
|
# DETECT STRUCTURE FIRST - before any use of is_per_field
|
|
week_file <- sprintf("week_%02d_%d.tif", current_week, current_year)
|
|
field_dirs <- list.dirs(current_mosaic_dir, full.names = FALSE, recursive = FALSE)
|
|
field_dirs <- field_dirs[field_dirs != ""]
|
|
|
|
is_per_field <- length(field_dirs) > 0 &&
|
|
file.exists(file.path(current_mosaic_dir, field_dirs[1], week_file))
|
|
|
|
if (is_per_field) {
|
|
message("Detected per-field mosaic structure...")
|
|
message("Using field-by-field extraction (similar to cane supply workflow)...")
|
|
|
|
# Use the same extraction method as cane supply
|
|
current_stats <- calculate_field_statistics(
|
|
field_boundaries_sf,
|
|
current_week,
|
|
current_year,
|
|
current_mosaic_dir,
|
|
report_date = Sys.Date()
|
|
)
|
|
|
|
# Extract CI pixels for each field from their individual mosaics
|
|
ci_pixels_by_field <- list()
|
|
for (i in seq_len(nrow(field_boundaries_sf))) {
|
|
field_name <- field_boundaries_sf$field[i]
|
|
field_mosaic_path <- file.path(current_mosaic_dir, field_name, week_file)
|
|
|
|
if (file.exists(field_mosaic_path)) {
|
|
tryCatch({
|
|
field_raster <- terra::rast(field_mosaic_path)
|
|
ci_band <- field_raster[["CI"]]
|
|
field_vect <- terra::vect(field_boundaries_sf[i, ])
|
|
ci_pixels_by_field[[i]] <- extract_ci_values(ci_band, field_vect)
|
|
}, error = function(e) {
|
|
message(paste(" Warning: Could not extract CI for field", field_name, ":", e$message))
|
|
ci_pixels_by_field[[i]] <- NULL
|
|
})
|
|
} else {
|
|
ci_pixels_by_field[[i]] <- NULL
|
|
}
|
|
}
|
|
|
|
# For uniformity calculations that need a reference raster, load first available
|
|
current_mosaic <- NULL
|
|
for (field_name in field_dirs) {
|
|
field_mosaic_path <- file.path(current_mosaic_dir, field_name, week_file)
|
|
if (file.exists(field_mosaic_path)) {
|
|
tryCatch({
|
|
current_mosaic <- terra::rast(field_mosaic_path)[["CI"]]
|
|
break
|
|
}, error = function(e) {
|
|
next
|
|
})
|
|
}
|
|
}
|
|
|
|
} else {
|
|
# Single-file mosaic (original behavior)
|
|
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")
|
|
}
|
|
|
|
message("Extracting field statistics from current mosaic...")
|
|
current_stats <- extract_field_statistics_from_ci(current_mosaic, field_boundaries_sf)
|
|
|
|
# Extract CI pixels for each field individually
|
|
ci_pixels_by_field <- list()
|
|
for (i in seq_len(nrow(field_boundaries_sf))) {
|
|
field_vect <- terra::vect(field_boundaries_sf[i, ])
|
|
ci_pixels_by_field[[i]] <- extract_ci_values(current_mosaic, field_vect)
|
|
}
|
|
}
|
|
|
|
# Load previous week mosaic (if available)
|
|
previous_stats <- NULL
|
|
if (!is.null(previous_mosaic_dir) || is_per_field) {
|
|
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, ")..."))
|
|
|
|
if (is_per_field) {
|
|
# Try loading previous week from the same directory structure
|
|
prev_week_file <- sprintf("week_%02d_%d.tif", target_prev$week, target_prev$year)
|
|
prev_field_exists <- any(sapply(field_dirs, function(field) {
|
|
file.exists(file.path(current_mosaic_dir, field, prev_week_file))
|
|
}))
|
|
|
|
if (prev_field_exists) {
|
|
message(" Found previous week per-field mosaics, calculating statistics...")
|
|
previous_stats <- calculate_field_statistics(
|
|
field_boundaries_sf,
|
|
target_prev$week,
|
|
target_prev$year,
|
|
current_mosaic_dir,
|
|
report_date = Sys.Date() - 7
|
|
)
|
|
} else {
|
|
message(" Previous week mosaic not available - skipping area change KPI")
|
|
}
|
|
} else if (!is.null(previous_mosaic_dir)) {
|
|
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...")
|
|
if (is_per_field) {
|
|
uniformity_kpi <- calculate_field_uniformity_kpi(
|
|
ci_pixels_by_field,
|
|
field_boundaries_sf,
|
|
ci_band = NULL,
|
|
mosaic_dir = current_mosaic_dir,
|
|
week_file = week_file
|
|
)
|
|
} else {
|
|
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,
|
|
data_dir = data_dir, project_dir = project_dir)
|
|
|
|
message("Calculating KPI 4: Growth Decline...")
|
|
growth_decline_kpi <- calculate_growth_decline_kpi(
|
|
ci_pixels_by_field
|
|
)
|
|
|
|
message("Calculating KPI 5: Field Patchiness...")
|
|
weed_kpi <- calculate_weed_presence_kpi(ci_pixels_by_field)
|
|
|
|
# Convert weed metrics to patchiness metrics (Gini-like + risk combination)
|
|
mean_ci_values <- current_stats$Mean_CI
|
|
if (is.null(mean_ci_values) || length(mean_ci_values) != nrow(field_boundaries_sf)) {
|
|
mean_ci_values <- rep(NA_real_, nrow(field_boundaries_sf))
|
|
}
|
|
patchiness_kpi <- calculate_patchiness_kpi(weed_kpi, ci_pixels_by_field, mean_ci_values)
|
|
|
|
message("Calculating KPI 6: Gap Filling...")
|
|
# Build list of per-field files for this week
|
|
per_field_files <- c()
|
|
for (field_name in field_dirs) {
|
|
field_mosaic_path <- file.path(current_mosaic_dir, field_name, week_file)
|
|
if (file.exists(field_mosaic_path)) {
|
|
per_field_files <- c(per_field_files, field_mosaic_path)
|
|
}
|
|
}
|
|
|
|
if (length(per_field_files) > 0) {
|
|
# Use the common wrapper function (same as cane supply)
|
|
gap_scores_result <- calculate_gap_scores(per_field_files, field_boundaries_sf)
|
|
|
|
# Convert to the format expected by orchestrator
|
|
gap_filling_kpi <- gap_scores_result %>%
|
|
mutate(field_idx = match(Field_id, field_boundaries_sf$field)) %>%
|
|
select(field_idx, gap_score) %>%
|
|
mutate(
|
|
gap_level = dplyr::case_when(
|
|
gap_score < 10 ~ "Minimal",
|
|
gap_score < 25 ~ "Moderate",
|
|
TRUE ~ "Significant"
|
|
),
|
|
mean_ci = NA_real_,
|
|
outlier_threshold = NA_real_
|
|
)
|
|
} else {
|
|
# Fallback: no per-field files
|
|
gap_filling_kpi <- data.frame(
|
|
field_idx = seq_len(nrow(field_boundaries_sf)),
|
|
gap_score = NA_real_,
|
|
gap_level = NA_character_,
|
|
mean_ci = NA_real_,
|
|
outlier_threshold = NA_real_
|
|
)
|
|
}
|
|
|
|
# Compile results
|
|
all_kpis <- list(
|
|
uniformity = uniformity_kpi,
|
|
area_change = area_change_kpi,
|
|
tch_forecasted = tch_kpi,
|
|
growth_decline = growth_decline_kpi,
|
|
patchiness = patchiness_kpi,
|
|
spatial_clustering = uniformity_kpi %>% select(field_idx, morans_i),
|
|
gap_filling = gap_filling_kpi
|
|
)
|
|
|
|
# Deduplicate KPI dataframes to ensure one row per field_idx
|
|
# (sometimes joins or calculations can create duplicate rows)
|
|
message("Deduplicating KPI results (keeping first occurrence per field)...")
|
|
all_kpis$uniformity <- all_kpis$uniformity %>%
|
|
distinct(field_idx, .keep_all = TRUE)
|
|
all_kpis$area_change <- all_kpis$area_change %>%
|
|
distinct(field_idx, .keep_all = TRUE)
|
|
all_kpis$tch_forecasted <- all_kpis$tch_forecasted %>%
|
|
distinct(field_idx, .keep_all = TRUE)
|
|
all_kpis$growth_decline <- all_kpis$growth_decline %>%
|
|
distinct(field_idx, .keep_all = TRUE)
|
|
all_kpis$patchiness <- all_kpis$patchiness %>%
|
|
distinct(field_idx, .keep_all = TRUE)
|
|
if (!is.null(all_kpis$spatial_clustering)) {
|
|
all_kpis$spatial_clustering <- all_kpis$spatial_clustering %>%
|
|
distinct(field_idx, .keep_all = TRUE)
|
|
}
|
|
all_kpis$gap_filling <- all_kpis$gap_filling %>%
|
|
distinct(field_idx, .keep_all = TRUE)
|
|
|
|
# Built single-sheet field detail table with all KPIs
|
|
message("\nBuilding comprehensive field detail table...")
|
|
field_detail_df <- create_field_detail_table(
|
|
field_boundaries_sf = field_boundaries_sf,
|
|
all_kpis = all_kpis,
|
|
current_week = current_week,
|
|
current_year = current_year
|
|
)
|
|
|
|
# Create summary tables
|
|
message("\nCreating summary tables...")
|
|
kpi_summary <- create_summary_tables(all_kpis)
|
|
|
|
# Export
|
|
message("\nExporting KPI data (single-sheet format)...")
|
|
export_paths <- export_kpi_data(
|
|
field_detail_df = field_detail_df,
|
|
kpi_summary = kpi_summary,
|
|
output_dir = output_dir,
|
|
week = current_week,
|
|
year = current_year,
|
|
project_dir = project_dir
|
|
)
|
|
|
|
message(paste("\n✓ AURA KPI calculation complete. Week", current_week, current_year))
|
|
|
|
return(list(
|
|
field_analysis_df = field_detail_df,
|
|
kpis = all_kpis,
|
|
summary_tables = kpi_summary,
|
|
export_paths = export_paths,
|
|
metadata = list(
|
|
week = current_week,
|
|
year = current_year,
|
|
project = project_dir
|
|
)
|
|
))
|
|
}
|