# CROP_MESSAGING_UTILS.R
# ======================
# Utility functions for the SmartCane crop messaging workflow.
# These functions support crop analysis, messaging, and output generation.

#' Convert hectares to acres
#' @param hectares Numeric value in hectares
#' @return Numeric value in acres
hectares_to_acres <- function(hectares) {
  return(hectares * 2.47105)
}

#' Format area with both hectares and acres
#' @param hectares Numeric value in hectares
#' @param precision Number of decimal places (default 1)
#' @return Character string with both measurements
format_area_both <- function(hectares, precision = 1) {
  acres <- hectares_to_acres(hectares)
  return(sprintf("%.1f ha (%.0f acres)", hectares, acres))
}
#' @param message The message to log
#' @param level The log level (default: "INFO")
#' @return NULL (used for side effects)
#'
safe_log <- function(message, level = "INFO") {
  if (exists("log_message")) {
    log_message(message, level)
  } else {
    if (level %in% c("ERROR", "WARNING")) {
      warning(message)
    } else {
      message(message)
    }
  }
}

# 2. Analysis configuration
# -----------------------
# Thresholds for change detection
CI_CHANGE_INCREASE_THRESHOLD <- 0.5
CI_CHANGE_DECREASE_THRESHOLD <- -0.5

# Thresholds for field uniformity (coefficient of variation as decimal)
UNIFORMITY_THRESHOLD <- 0.15  # Below this = good uniformity, above = requires attention
EXCELLENT_UNIFORMITY_THRESHOLD <- 0.08  # Below this = excellent uniformity
POOR_UNIFORMITY_THRESHOLD <- 0.25  # Above this = poor uniformity, urgent attention needed

# Thresholds for spatial clustering (adjusted for agricultural fields)
# Agricultural fields naturally have spatial autocorrelation, so higher thresholds are needed
MORAN_THRESHOLD_HIGH <- 0.95   # Above this = very strong clustering (problematic patterns)
MORAN_THRESHOLD_MODERATE <- 0.85  # Above this = moderate clustering
MORAN_THRESHOLD_LOW <- 0.7     # Above this = normal field continuity

# Threshold for acceptable area percentage
ACCEPTABLE_AREA_THRESHOLD <- 40  # Below this percentage = management issue

#' Calculate uniformity metrics using terra statistics (optimized)
#' @param mean_val Mean CI value from terra
#' @param sd_val Standard deviation from terra
#' @param median_val Median CI value from terra
#' @param min_val Minimum CI value from terra
#' @param max_val Maximum CI value from terra
#' @param values Raw values for quantile calculations only
#' @return List with various uniformity metrics (all scaled to be comparable)
calculate_uniformity_metrics_terra <- function(mean_val, sd_val, median_val, min_val, max_val, values) {

  if (is.na(mean_val) || length(values) < 2) return(list(
    cv = NA, iqr_cv = NA, range_cv = NA,
    mad_cv = NA, percentile_cv = NA, interpretation = "insufficient_data"
  ))

  # 1. Coefficient of variation (from terra) - already normalized
  cv <- sd_val / mean_val

  # 2. IQR-based CV (IQR/median) - using R's built-in IQR function
  iqr_val <- IQR(values, na.rm = TRUE)
  iqr_cv <- iqr_val / median_val

  # 3. Range-based CV (range/mean) - using terra min/max
  range_val <- max_val - min_val
  range_cv <- range_val / mean_val

  # 4. MAD-based CV (MAD/median) - using R's built-in mad function
  mad_val <- mad(values, constant = 1.4826, na.rm = TRUE)  # scaled to match SD for normal distribution
  mad_cv <- mad_val / median_val

  # 5. Percentile-based CV (P90-P10)/mean - using R's built-in quantile
  percentiles <- quantile(values, c(0.1, 0.9), na.rm = TRUE)
  percentile_cv <- (percentiles[2] - percentiles[1]) / mean_val

  # Interpretation based on CV thresholds (all metrics now comparable)
  # CV < 0.15 = Very uniform, 0.15-0.30 = Moderate variation, 0.30-0.50 = High variation, >0.50 = Very high variation
  interpret_uniformity <- function(metric_value) {
    if (is.na(metric_value)) return("unknown")
    if (metric_value < 0.15) return("very uniform")
    if (metric_value < 0.30) return("moderate variation")
    if (metric_value < 0.50) return("high variation")
    return("very high variation")
  }

  return(list(
    cv = cv,
    iqr_cv = iqr_cv,
    range_cv = range_cv,
    mad_cv = mad_cv,
    percentile_cv = percentile_cv,
    cv_interpretation = interpret_uniformity(cv),
    iqr_interpretation = interpret_uniformity(iqr_cv),
    mad_interpretation = interpret_uniformity(mad_cv),
    percentile_interpretation = interpret_uniformity(percentile_cv)
  ))
}

#' Calculate percentage within acceptable range using terra mean
#' Acceptable range = within 25% of the field mean CI value
#' This indicates what percentage of the field has "normal" performance
#' @param mean_val Mean CI value from terra
#' @param values Raw CI values
#' @param threshold_factor Factor to multiply mean by for acceptable range (default 0.25 = 25%)
#' @return Percentage of values within acceptable range
calculate_acceptable_percentage_terra <- function(mean_val, values, threshold_factor = 0.25) {
  values <- values[!is.na(values) & is.finite(values)]
  if (length(values) < 2 || is.na(mean_val)) return(NA)

  threshold <- mean_val * threshold_factor  # 25% of mean as default
  within_range <- abs(values - mean_val) <= threshold
  percentage <- (sum(within_range) / length(values)) * 100
  return(percentage)
}

#' Calculate coefficient of variation for uniformity assessment
#' @param values Numeric vector of CI values
#' @return Coefficient of variation (CV) as decimal
calculate_cv <- function(values) {
  values <- values[!is.na(values) & is.finite(values)]
  if (length(values) < 2) return(NA)
  cv <- sd(values) / mean(values)  # Keep as decimal
  return(cv)
}

#' Calculate Shannon entropy for spatial heterogeneity assessment
#' Higher entropy = more heterogeneous/variable field
#' Lower entropy = more homogeneous/uniform field
#' @param values Numeric vector of CI values
#' @param n_bins Number of bins for histogram (default 10)
#' @return Shannon entropy value
calculate_entropy <- function(values, n_bins = 10) {
  values <- values[!is.na(values) & is.finite(values)]
  if (length(values) < 2) return(NA)

  # Create histogram bins
  value_range <- range(values)
  breaks <- seq(value_range[1], value_range[2], length.out = n_bins + 1)

  # Count values in each bin
  bin_counts <- hist(values, breaks = breaks, plot = FALSE)$counts

  # Calculate probabilities (remove zero counts)
  probabilities <- bin_counts[bin_counts > 0] / sum(bin_counts)

  # Calculate Shannon entropy: H = -sum(p * log(p))
  entropy <- -sum(probabilities * log(probabilities))

  return(entropy)
}

#' Calculate percentage of field with positive vs negative change
#' @param current_values Current week CI values
#' @param previous_values Previous week CI values
#' @return List with percentage of positive and negative change areas
calculate_change_percentages <- function(current_values, previous_values) {
  # Ensure same length (should be from same field boundaries)
  if (length(current_values) != length(previous_values)) {
    return(list(positive_pct = NA, negative_pct = NA, stable_pct = NA))
  }

  # Calculate pixel-wise change
  change_values <- current_values - previous_values
  valid_changes <- change_values[!is.na(change_values) & is.finite(change_values)]

  if (length(valid_changes) < 2) {
    return(list(positive_pct = NA, negative_pct = NA, stable_pct = NA))
  }

  # Count positive, negative, and stable areas
  positive_pct <- sum(valid_changes > 0) / length(valid_changes) * 100
  negative_pct <- sum(valid_changes < 0) / length(valid_changes) * 100
  stable_pct <- sum(valid_changes == 0) / length(valid_changes) * 100

  return(list(
    positive_pct = positive_pct,
    negative_pct = negative_pct,
    stable_pct = stable_pct
  ))
}

#' Calculate spatial autocorrelation (Moran's I) for a field
#' @param ci_raster Terra raster of CI values
#' @param field_boundary Terra vector of field boundary
#' @return List with Moran's I statistic and p-value
calculate_spatial_autocorrelation <- function(ci_raster, field_boundary) {

  tryCatch({
    # Crop and mask raster to field boundary
    field_raster <- terra::crop(ci_raster, field_boundary)
    field_raster <- terra::mask(field_raster, field_boundary)

    # Convert to points for spatial analysis
    raster_points <- terra::as.points(field_raster, na.rm = TRUE)

    # Check if we have enough points
    if (length(raster_points) < 10) {
      return(list(morans_i = NA, p_value = NA, interpretation = "insufficient_data"))
    }

    # Convert to sf for spdep
    points_sf <- sf::st_as_sf(raster_points)

    # Create spatial weights matrix (k-nearest neighbors)
    coords <- sf::st_coordinates(points_sf)

    # Use adaptive number of neighbors based on sample size
    k_neighbors <- min(8, max(4, floor(nrow(coords) / 10)))

    knn_nb <- spdep::knearneigh(coords, k = k_neighbors)
    knn_listw <- spdep::nb2listw(spdep::knn2nb(knn_nb), style = "W", zero.policy = TRUE)

    # Calculate Moran's I
    ci_values <- points_sf[[1]]  # First column contains CI values
    moran_result <- spdep::moran.test(ci_values, knn_listw, zero.policy = TRUE)

    # Interpret results
    morans_i <- moran_result$estimate[1]
    p_value <- moran_result$p.value

    interpretation <- if (is.na(morans_i)) {
      "insufficient_data"
    } else if (p_value > 0.05) {
      "random"  # Not significant spatial pattern
    } else if (morans_i > MORAN_THRESHOLD_HIGH) {
      "very_strong_clustering"  # Very strong clustering - may indicate management issues
    } else if (morans_i > MORAN_THRESHOLD_MODERATE) {
      "strong_clustering"  # Strong clustering - worth monitoring
    } else if (morans_i > MORAN_THRESHOLD_LOW) {
      "normal_continuity"  # Normal field continuity - expected for uniform fields
    } else if (morans_i > 0.3) {
      "weak_clustering"  # Some clustering present
    } else if (morans_i < -0.3) {
      "dispersed"  # Checkerboard pattern
    } else {
      "low_autocorrelation"  # Low spatial autocorrelation
    }

    return(list(
      morans_i = morans_i,
      p_value = p_value,
      interpretation = interpretation
    ))

  }, error = function(e) {
    warning(paste("Error calculating spatial autocorrelation:", e$message))
    return(list(morans_i = NA, p_value = NA, interpretation = "error"))
  })
}

#' Calculate percentage of field in extreme values using simple threshold
#' Hotspots = areas with CI > mean + 1.5*SD (high-performing areas)
#' Coldspots = areas with CI < mean - 1.5*SD (underperforming areas)
#' @param values Numeric vector of CI values
#' @param threshold_multiplier Standard deviation multiplier (default 1.5)
#' @return List with percentage of hotspots and coldspots
calculate_extreme_percentages_simple <- function(values, threshold_multiplier = 1.5) {

  if (length(values) < 10) return(list(hotspot_pct = NA, coldspot_pct = NA, method = "insufficient_data"))

  mean_val <- mean(values, na.rm = TRUE)
  sd_val <- sd(values, na.rm = TRUE)

  # Hotspots: significantly ABOVE average (good performance)
  upper_threshold <- mean_val + (threshold_multiplier * sd_val)
  # Coldspots: significantly BELOW average (poor performance)
  lower_threshold <- mean_val - (threshold_multiplier * sd_val)

  hotspot_pct <- sum(values > upper_threshold, na.rm = TRUE) / length(values) * 100
  coldspot_pct <- sum(values < lower_threshold, na.rm = TRUE) / length(values) * 100

  return(list(
    hotspot_pct = hotspot_pct,
    coldspot_pct = coldspot_pct,
    method = "simple_threshold",
    threshold_used = threshold_multiplier
  ))
}

#' Categorize CI change based on thresholds
#' @param change_value Mean change in CI between weeks
#' @return Character string: "increase", "stable", or "decrease"
categorize_change <- function(change_value) {
  if (is.na(change_value)) return("unknown")
  if (change_value >= CI_CHANGE_INCREASE_THRESHOLD) return("increase")
  if (change_value <= CI_CHANGE_DECREASE_THRESHOLD) return("decrease")
  return("stable")
}

#' Categorize field uniformity based on coefficient of variation and spatial pattern
#' @param cv_value Coefficient of variation (primary uniformity metric)
#' @param spatial_info List with spatial autocorrelation results
#' @param extreme_pct List with hotspot/coldspot percentages
#' @param acceptable_pct Percentage of field within acceptable range
#' @return Character string describing field uniformity pattern
categorize_uniformity_enhanced <- function(cv_value, spatial_info, extreme_pct, acceptable_pct = NA) {

  if (is.na(cv_value)) return("unknown variation")

  # Check for poor uniformity first (urgent issues)
  if (cv_value > POOR_UNIFORMITY_THRESHOLD || (!is.na(acceptable_pct) && acceptable_pct < ACCEPTABLE_AREA_THRESHOLD)) {
    return("poor uniformity - urgent attention needed")
  }

  # Check for excellent uniformity
  if (cv_value <= EXCELLENT_UNIFORMITY_THRESHOLD && (!is.na(acceptable_pct) && acceptable_pct >= 45)) {
    return("excellent uniformity")
  }

  # Check for good uniformity
  if (cv_value <= UNIFORMITY_THRESHOLD) {
    return("good uniformity")
  }

  # Field has moderate variation - determine if localized or distributed
  spatial_pattern <- spatial_info$interpretation
  hotspot_pct <- extreme_pct$hotspot_pct
  coldspot_pct <- extreme_pct$coldspot_pct

  # Determine pattern type based on CV (primary) and spatial pattern (secondary)
  if (spatial_pattern %in% c("very_strong_clustering") && !is.na(hotspot_pct) && (hotspot_pct > 15 || coldspot_pct > 5)) {
    # Very strong clustering with substantial extreme areas - likely problematic
    if (hotspot_pct > coldspot_pct) {
      return("localized high-performing areas")
    } else if (coldspot_pct > hotspot_pct) {
      return("localized problem areas")
    } else {
      return("localized hotspots and coldspots")
    }
  } else if (spatial_pattern %in% c("strong_clustering") && !is.na(hotspot_pct) && (hotspot_pct > 10 || coldspot_pct > 3)) {
    # Strong clustering with moderate extreme areas
    if (hotspot_pct > coldspot_pct) {
      return("localized high-performing areas")
    } else if (coldspot_pct > hotspot_pct) {
      return("localized problem areas")
    } else {
      return("clustered variation")
    }
  } else {
    # Normal field continuity or weak patterns - rely primarily on CV
    return("moderate variation")
  }
}

#' Generate enhanced message based on analysis results including spatial patterns
#' @param uniformity_category Character: enhanced uniformity category with spatial info
#' @param change_category Character: "increase", "stable", or "decrease"
#' @param extreme_pct List with hotspot/coldspot percentages
#' @param acceptable_pct Percentage of field within acceptable range
#' @param morans_i Moran's I value for additional context
#' @param growth_stage Character: growth stage (simplified for now)
#' @return List with message and worth_sending flag
generate_enhanced_message <- function(uniformity_category, change_category, extreme_pct, acceptable_pct = NA, morans_i = NA, growth_stage = "vegetation stage") {

  # Enhanced message matrix based on spatial patterns
  messages <- list()

  # Poor uniformity scenarios (urgent)
  if (uniformity_category == "poor uniformity - urgent attention needed") {
    messages <- list(
      "stable" = list(
        message = "🚨 URGENT: Poor field uniformity detected - immediate management review required",
        worth_sending = TRUE
      ),
      "decrease" = list(
        message = "🚨 CRITICAL: Poor uniformity with declining trend - emergency intervention needed",
        worth_sending = TRUE
      ),
      "increase" = list(
        message = "⚠️ CAUTION: Improving but still poor uniformity - continue intensive monitoring",
        worth_sending = TRUE
      )
    )
  }

  # Excellent uniformity scenarios
  else if (uniformity_category == "excellent uniformity") {
    messages <- list(
      "stable" = list(
        message = "✅ Excellent: Optimal field uniformity and stability",
        worth_sending = FALSE
      ),
      "decrease" = list(
        message = "⚠️ Alert: Excellent uniformity but declining - investigate cause early",
        worth_sending = TRUE
      ),
      "increase" = list(
        message = "🌟 Outstanding: Excellent uniformity with continued improvement",
        worth_sending = FALSE
      )
    )
  }

  # Good uniformity scenarios
  else if (uniformity_category == "good uniformity") {
    # Check for very strong clustering which may indicate management issues
    if (!is.na(morans_i) && morans_i > MORAN_THRESHOLD_HIGH) {
      messages <- list(
        "stable" = list(
          message = "⚠️ Alert: Good uniformity but very strong clustering detected - check management practices",
          worth_sending = TRUE
        ),
        "decrease" = list(
          message = "🚨 Alert: Good uniformity declining with clustering patterns - targeted intervention needed",
          worth_sending = TRUE
        ),
        "increase" = list(
          message = "✅ Good: Improving uniformity but monitor clustering patterns",
          worth_sending = FALSE
        )
      )
    } else {
      messages <- list(
        "stable" = list(
          message = "✅ Good: Stable field with good uniformity",
          worth_sending = FALSE
        ),
        "decrease" = list(
          message = "⚠️ Alert: Good uniformity but declining trend - early intervention recommended",
          worth_sending = TRUE
        ),
        "increase" = list(
          message = "✅ Great: Good uniformity with improvement trend",
          worth_sending = FALSE
        )
      )
    }
  }

  # Moderate variation scenarios
  else if (uniformity_category == "moderate variation") {
    acceptable_msg <- if (!is.na(acceptable_pct) && acceptable_pct < 45) " - low acceptable area" else ""

    messages <- list(
      "stable" = list(
        message = paste0("⚠️ Alert: Moderate field variation detected", acceptable_msg, " - review management uniformity"),
        worth_sending = TRUE
      ),
      "decrease" = list(
        message = paste0("🚨 Alert: Moderate variation with declining trend", acceptable_msg, " - intervention needed"),
        worth_sending = TRUE
      ),
      "increase" = list(
        message = paste0("📈 Monitor: Improving but still moderate variation", acceptable_msg, " - continue optimization"),
        worth_sending = FALSE
      )
    )
  }

  # Localized problem areas
  else if (uniformity_category == "localized problem areas") {
    hotspot_pct <- round(extreme_pct$hotspot_pct, 1)
    coldspot_pct <- round(extreme_pct$coldspot_pct, 1)

    messages <- list(
      "stable" = list(
        message = paste0("🚨 Alert: Problem zones detected (", coldspot_pct, "% underperforming) - targeted intervention needed"),
        worth_sending = TRUE
      ),
      "decrease" = list(
        message = paste0("🚨 URGENT: Problem areas expanding with overall decline (", coldspot_pct, "% affected) - immediate action required"),
        worth_sending = TRUE
      ),
      "increase" = list(
        message = paste0("⚠️ Caution: Overall improvement but ", coldspot_pct, "% problem areas remain - monitor closely"),
        worth_sending = TRUE
      )
    )
  }

  # Localized high-performing areas
  else if (uniformity_category == "localized high-performing areas") {
    hotspot_pct <- round(extreme_pct$hotspot_pct, 1)

    messages <- list(
      "stable" = list(
        message = paste0("💡 Opportunity: ", hotspot_pct, "% of field performing well - replicate conditions in remaining areas"),
        worth_sending = FALSE
      ),
      "decrease" = list(
        message = paste0("⚠️ Alert: High-performing areas (", hotspot_pct, "%) declining - investigate cause to prevent spread"),
        worth_sending = TRUE
      ),
      "increase" = list(
        message = paste0("🌟 Excellent: High-performing areas (", hotspot_pct, "%) expanding - excellent management practices"),
        worth_sending = FALSE
      )
    )
  }

  # Clustered variation (general)
  else if (uniformity_category == "clustered variation") {
    messages <- list(
      "stable" = list(
        message = "⚠️ Alert: Clustered variation detected - investigate spatial management patterns",
        worth_sending = TRUE
      ),
      "decrease" = list(
        message = "🚨 Alert: Clustered decline pattern - targeted investigation needed",
        worth_sending = TRUE
      ),
      "increase" = list(
        message = "📈 Monitor: Clustered improvement - identify and replicate successful practices",
        worth_sending = FALSE
      )
    )
  }

  # Default fallback
  else {
    messages <- list(
      "stable" = list(message = "❓ Field analysis inconclusive - manual review recommended", worth_sending = FALSE),
      "decrease" = list(message = "⚠️ Field showing decline - investigation recommended", worth_sending = TRUE),
      "increase" = list(message = "📈 Field showing improvement", worth_sending = FALSE)
    )
  }

  # Return appropriate message
  if (change_category %in% names(messages)) {
    return(messages[[change_category]])
  } else {
    return(list(
      message = paste("❓ Analysis inconclusive -", uniformity_category, "with", change_category, "trend"),
      worth_sending = FALSE
    ))
  }
}

#' Load and analyze a weekly mosaic for individual fields with spatial analysis
#' @param week_file_path Path to the weekly mosaic file
#' @param field_boundaries_sf SF object with field boundaries
#' @return List with CI statistics per field including spatial metrics
analyze_weekly_mosaic <- function(week_file_path, field_boundaries_sf) {

  if (!file.exists(week_file_path)) {
    warning(paste("Mosaic file not found:", week_file_path))
    return(NULL)
  }

  tryCatch({
    # Load the raster and select only the CI band (5th band)
    mosaic_raster <- terra::rast(week_file_path)
    ci_raster <- mosaic_raster[[5]]  # Select the CI band
    names(ci_raster) <- "CI"

    # Convert field boundaries to terra vect for extraction
    field_boundaries_vect <- terra::vect(field_boundaries_sf)

    # Extract CI values for each field
    field_results <- list()

    for (i in seq_len(nrow(field_boundaries_sf))) {
      field_name <- field_boundaries_sf$field[i]
      sub_field_name <- field_boundaries_sf$sub_field[i]

      # Check and get field area from geojson if available
      field_area_ha <- NA
      if ("area_ha" %in% colnames(field_boundaries_sf)) {
        field_area_ha <- field_boundaries_sf$area_ha[i]
      } else if ("AREA_HA" %in% colnames(field_boundaries_sf)) {
        field_area_ha <- field_boundaries_sf$AREA_HA[i]
      } else if ("area" %in% colnames(field_boundaries_sf)) {
        field_area_ha <- field_boundaries_sf$area[i]
      } else {
        # Calculate area from geometry as fallback
        field_geom <- field_boundaries_sf[i,]
        if (sf::st_is_longlat(field_geom)) {
          # For geographic coordinates, transform to projected for area calculation
          field_geom <- sf::st_transform(field_geom, 3857)  # Web Mercator
        }
        field_area_ha <- as.numeric(sf::st_area(field_geom)) / 10000  # Convert to hectares
      }

      cat("Processing field:", field_name, "-", sub_field_name, "(", round(field_area_ha, 1), "ha)\n")

      # Extract values for this specific field
      field_vect <- field_boundaries_vect[i]

      # Extract with built-in statistics from terra (PRIMARY METHOD)
      terra_stats <- terra::extract(ci_raster, field_vect, fun = c("mean", "sd", "min", "max", "median"), na.rm = TRUE)

      # Extract raw values for additional calculations and validation
      ci_values <- terra::extract(ci_raster, field_vect, fun = NULL)

      # Flatten and clean the values
      field_values <- unlist(ci_values)
      valid_values <- field_values[!is.na(field_values) & is.finite(field_values)]

      if (length(valid_values) > 0) {

        # Use TERRA as primary calculations
        primary_mean <- terra_stats$mean[1]
        primary_sd <- terra_stats$sd[1]
        primary_cv <- primary_sd / primary_mean
        primary_median <- terra_stats$median[1]
        primary_min <- terra_stats$min[1]
        primary_max <- terra_stats$max[1]

        # Manual calculations for validation only
        manual_mean <- mean(valid_values)
        manual_cv <- sd(valid_values) / manual_mean

        basic_stats <- list(
          field = field_name,
          sub_field = sub_field_name,
          # PRIMARY statistics (terra-based)
          mean_ci = primary_mean,
          median_ci = primary_median,
          sd_ci = primary_sd,
          cv = primary_cv,
          min_ci = primary_min,
          max_ci = primary_max,
          # Store raw values for change analysis
          raw_values = valid_values,
          # Other metrics using terra values
          acceptable_pct = calculate_acceptable_percentage_terra(primary_mean, valid_values),
          n_pixels = length(valid_values),
          # Field area from geojson
          field_area_ha = field_area_ha
        )

        # Calculate spatial statistics
        spatial_info <- calculate_spatial_autocorrelation(ci_raster, field_vect)
        extreme_pct <- calculate_extreme_percentages_simple(valid_values)

        # Calculate entropy for additional uniformity measure
        entropy_value <- calculate_entropy(valid_values)

        # Enhanced uniformity categorization
        uniformity_category <- categorize_uniformity_enhanced(
          basic_stats$cv,
          spatial_info,
          extreme_pct,
          basic_stats$acceptable_pct
        )

        # Combine all results
        field_stats <- c(
          basic_stats,
          list(
            spatial_autocorr = spatial_info,
            extreme_percentages = extreme_pct,
            entropy = entropy_value,
            uniformity_category = uniformity_category
          )
        )

        field_results[[paste0(field_name, "_", sub_field_name)]] <- field_stats

      } else {
        warning(paste("No valid CI values found for field:", field_name, sub_field_name))
      }
    }

    return(field_results)

  }, error = function(e) {
    warning(paste("Error analyzing mosaic:", e$message))
    return(NULL)
  })
}

#' Run crop analysis for any estate
#' @param estate_name Character: name of the estate (e.g., "simba", "chemba")
#' @param current_week Numeric: current week number
#' @param previous_week Numeric: previous week number
#' @param year Numeric: year (default 2025)
#' @return List with analysis results
run_estate_analysis <- function(estate_name, current_week, previous_week, year = 2025) {

  cat("=== CROP ANALYSIS MESSAGING SYSTEM ===\n")
  cat("Analyzing:", toupper(estate_name), "estate\n")
  cat("Comparing week", previous_week, "vs week", current_week, "of", year, "\n\n")

  # Set project_dir globally for parameters_project.R
  assign("project_dir", estate_name, envir = .GlobalEnv)

  # Load project configuration
  tryCatch({
    source("parameters_project.R")
    cat("✓ Project configuration loaded\n")
  }, error = function(e) {
    tryCatch({
      source(here::here("r_app", "parameters_project.R"))
      cat("✓ Project configuration loaded from r_app directory\n")
    }, error = function(e) {
      stop("Failed to load project configuration")
    })
  })

  # Verify required variables are available
  if (!exists("weekly_CI_mosaic") || !exists("field_boundaries_sf")) {
    stop("Required project variables not initialized. Check project configuration.")
  }

  # Construct file paths for weekly mosaics
  current_week_file <- sprintf("week_%02d_%d.tif", current_week, year)
  previous_week_file <- sprintf("week_%02d_%d.tif", previous_week, year)

  current_week_path <- file.path(weekly_CI_mosaic, current_week_file)
  previous_week_path <- file.path(weekly_CI_mosaic, previous_week_file)

  cat("Looking for files:\n")
  cat("- Current week:", current_week_path, "\n")
  cat("- Previous week:", previous_week_path, "\n\n")

  # Check if files exist and handle missing data scenarios
  current_exists <- file.exists(current_week_path)
  previous_exists <- file.exists(previous_week_path)

  if (!current_exists) {
    cat("❌ Current week mosaic not found. No analysis possible.\n")
    return(NULL)
  }

  # Analyze both weeks for all fields
  cat("Analyzing weekly mosaics per field...\n")
  current_field_stats <- analyze_weekly_mosaic(current_week_path, field_boundaries_sf)

  if (!previous_exists) {
    cat("⚠️ Previous week mosaic not found (likely due to clouds). Performing spatial-only analysis.\n")
    previous_field_stats <- NULL
  } else {
    previous_field_stats <- analyze_weekly_mosaic(previous_week_path, field_boundaries_sf)
  }

  if (is.null(current_field_stats)) {
    stop("Could not analyze current weekly mosaic")
  }

  # Generate field results
  field_results <- generate_field_results(current_field_stats, previous_field_stats, current_week, previous_week)

  return(list(
    estate_name = estate_name,
    current_week = current_week,
    previous_week = previous_week,
    year = year,
    field_results = field_results,
    current_field_stats = current_field_stats,
    previous_field_stats = previous_field_stats
  ))
}

#' Generate analysis results for all fields
#' @param current_field_stats Analysis results for current week
#' @param previous_field_stats Analysis results for previous week
#' @param current_week Current week number
#' @param previous_week Previous week number
#' @return List with field results
generate_field_results <- function(current_field_stats, previous_field_stats, current_week, previous_week) {

  field_results <- list()

  # Get common field names between both weeks (or all current fields if previous is missing)
  if (!is.null(previous_field_stats)) {
    common_fields <- intersect(names(current_field_stats), names(previous_field_stats))
  } else {
    common_fields <- names(current_field_stats)
  }

  for (field_id in common_fields) {
    current_field <- current_field_stats[[field_id]]
    previous_field <- if (!is.null(previous_field_stats)) previous_field_stats[[field_id]] else NULL

    # Calculate change metrics for this field (only if previous data exists)
    if (!is.null(previous_field)) {
      ci_change <- current_field$mean_ci - previous_field$mean_ci
      change_category <- categorize_change(ci_change)

      # Calculate spatial change percentages
      change_percentages <- calculate_change_percentages(
        current_field$raw_values,
        previous_field$raw_values
      )
    } else {
      # No previous data - spatial analysis only
      ci_change <- NA
      change_category <- "spatial_only"
      change_percentages <- list(positive_pct = NA, negative_pct = NA, stable_pct = NA)
    }

    # Use enhanced uniformity category from current week analysis
    uniformity_category <- current_field$uniformity_category

    # Generate enhanced message for this field
    message_result <- generate_enhanced_message(
      uniformity_category,
      change_category,
      current_field$extreme_percentages,
      current_field$acceptable_pct,
      current_field$spatial_autocorr$morans_i
    )

    # Store results
    field_results[[field_id]] <- list(
      current_stats = current_field,
      previous_stats = previous_field,
      ci_change = ci_change,
      change_category = change_category,
      change_percentages = change_percentages,
      uniformity_category = uniformity_category,
      message_result = message_result
    )
  }

  return(field_results)
}

#' Format analysis results for WhatsApp/Word copy-paste
#' @param analysis_results Results from run_estate_analysis
#' @return Character string with formatted text
format_for_whatsapp <- function(analysis_results) {

  field_results <- analysis_results$field_results
  estate_name <- toupper(analysis_results$estate_name)
  current_week <- analysis_results$current_week
  previous_week <- analysis_results$previous_week

  output <- c()
  output <- c(output, paste("🌾", estate_name, "CROP ANALYSIS"))
  output <- c(output, paste("📅 Week", current_week, "vs Week", previous_week))
  output <- c(output, "")

  # Summary statistics
  alert_count <- sum(sapply(field_results, function(x) x$message_result$worth_sending))
  total_fields <- length(field_results)

  # Calculate total area and area statistics
  total_hectares <- sum(sapply(field_results, function(x) x$current_stats$field_area_ha), na.rm = TRUE)

  output <- c(output, "📊 SUMMARY:")
  output <- c(output, paste("• Estate:", estate_name))
  output <- c(output, paste("• Fields analyzed:", total_fields))
  output <- c(output, paste("• Total area:", format_area_both(total_hectares)))
  output <- c(output, paste("• Alerts needed:", alert_count))
  output <- c(output, "")

  # Field-by-field alerts only
  if (alert_count > 0) {
    output <- c(output, "🚨 PRIORITY FIELDS:")
    for (field_id in names(field_results)) {
      field_info <- field_results[[field_id]]
      if (field_info$message_result$worth_sending) {
        field_name <- paste(field_info$current_stats$field, field_info$current_stats$sub_field, sep="-")
        area <- field_info$current_stats$field_area_ha
        message <- field_info$message_result$message

        output <- c(output, paste("•", field_name, paste0("(", format_area_both(area), "):"), message))
      }
    }
  } else {
    output <- c(output, "✅ No urgent alerts - all fields stable")
  }

  # Quick farm summary
  output <- c(output, "")
  output <- c(output, "📈 QUICK STATS:")

  # Calculate improving vs declining areas (only if temporal data available)
  has_temporal_data <- any(sapply(field_results, function(x) !is.na(x$change_percentages$positive_pct)))

  if (has_temporal_data) {
    total_improving <- sum(sapply(field_results, function(x) {
      if (!is.na(x$change_percentages$positive_pct)) {
        (x$change_percentages$positive_pct / 100) * x$current_stats$field_area_ha
      } else 0
    }), na.rm = TRUE)

    total_declining <- sum(sapply(field_results, function(x) {
      if (!is.na(x$change_percentages$negative_pct)) {
        (x$change_percentages$negative_pct / 100) * x$current_stats$field_area_ha
      } else 0
    }), na.rm = TRUE)

    total_stable <- sum(sapply(field_results, function(x) {
      if (!is.na(x$change_percentages$stable_pct)) {
        (x$change_percentages$stable_pct / 100) * x$current_stats$field_area_ha
      } else 0
    }), na.rm = TRUE)

    improving_pct <- (total_improving / total_hectares) * 100
    declining_pct <- (total_declining / total_hectares) * 100
    stable_pct <- (total_stable / total_hectares) * 100

    output <- c(output, paste("• Improving areas:", format_area_both(total_improving), paste0("(", round(improving_pct, 1), "%)")))
    output <- c(output, paste("• Declining areas:", format_area_both(total_declining), paste0("(", round(declining_pct, 1), "%)")))
    output <- c(output, paste("• Stable areas:", format_area_both(total_stable), paste0("(", round(stable_pct, 1), "%)")))

    # Overall trend
    if (improving_pct > declining_pct) {
      trend_diff <- round(improving_pct - declining_pct, 1)
      output <- c(output, paste("• Trend: ✅ POSITIVE (+", trend_diff, "%)"))
    } else if (declining_pct > improving_pct) {
      trend_diff <- round(declining_pct - improving_pct, 1)
      output <- c(output, paste("• Trend: ⚠️ NEGATIVE (-", trend_diff, "%)"))
    } else {
      output <- c(output, "• Trend: ➖ BALANCED")
    }
  } else {
    output <- c(output, "• Analysis: Spatial patterns only (previous week data unavailable)")
  }

  # Add farm-wide analysis summary
  output <- c(output, "")
  output <- c(output, "=== FARM-WIDE ANALYSIS SUMMARY ===")
  output <- c(output, "")

  # Field uniformity statistics
  excellent_fields <- sum(sapply(field_results, function(x) x$current_stats$cv <= 0.08))
  good_fields <- sum(sapply(field_results, function(x) x$current_stats$cv > 0.08 & x$current_stats$cv <= 0.15))
  moderate_fields <- sum(sapply(field_results, function(x) x$current_stats$cv > 0.15 & x$current_stats$cv <= 0.30))
  poor_fields <- sum(sapply(field_results, function(x) x$current_stats$cv > 0.30))

  output <- c(output, "FIELD UNIFORMITY SUMMARY:")
  output <- c(output, "│ Uniformity Level    │ Count │ Percent │")
  output <- c(output, sprintf("│ Excellent (CV≤0.08) │ %5d │ %6.1f%% │", excellent_fields, (excellent_fields/total_fields)*100))
  output <- c(output, sprintf("│ Good (CV 0.08-0.15) │ %5d │ %6.1f%% │", good_fields, (good_fields/total_fields)*100))
  output <- c(output, sprintf("│ Moderate (CV 0.15-0.30) │ %5d │ %6.1f%% │", moderate_fields, (moderate_fields/total_fields)*100))
  output <- c(output, sprintf("│ Poor (CV>0.30)      │ %5d │ %6.1f%% │", poor_fields, (poor_fields/total_fields)*100))
  output <- c(output, sprintf("│ Total fields        │ %5d │ %6.1f%% │", total_fields, 100.0))
  output <- c(output, "")

  # Farm-wide area change summary
  output <- c(output, "FARM-WIDE AREA CHANGE SUMMARY:")
  output <- c(output, "│ Change Type         │ Area (ha/acres) │ Percent │")

  if (has_temporal_data) {
    output <- c(output, sprintf("│ Improving areas     │ %s │ %6.1f%% │", format_area_both(total_improving), improving_pct))
    output <- c(output, sprintf("│ Stable areas        │ %s │ %6.1f%% │", format_area_both(total_stable), stable_pct))
    output <- c(output, sprintf("│ Declining areas     │ %s │ %6.1f%% │", format_area_both(total_declining), declining_pct))
    output <- c(output, sprintf("│ Total area          │ %s │ %6.1f%% │", format_area_both(total_hectares), 100.0))
  } else {
    output <- c(output, "│ Improving areas     │     N/A         │    N/A │")
    output <- c(output, "│ Stable areas        │     N/A         │    N/A │")
    output <- c(output, "│ Declining areas     │     N/A         │    N/A │")
    output <- c(output, sprintf("│ Total area          │ %s │ %6.1f%% │", format_area_both(total_hectares), 100.0))
  }
  output <- c(output, "")

  # Key insights
  output <- c(output, "KEY INSIGHTS:")
  good_uniformity_pct <- ((excellent_fields + good_fields) / total_fields) * 100
  excellent_uniformity_pct <- (excellent_fields / total_fields) * 100

  output <- c(output, sprintf("• %d%% of fields have good uniformity (CV ≤ 0.15)", round(good_uniformity_pct)))
  output <- c(output, sprintf("• %d%% of fields have excellent uniformity (CV ≤ 0.08)", round(excellent_uniformity_pct)))

  if (has_temporal_data) {
    output <- c(output, sprintf("• %s (%.1f%%) of farm area is improving week-over-week", format_area_both(total_improving), improving_pct))
    output <- c(output, sprintf("• %s (%.1f%%) of farm area is stable week-over-week", format_area_both(total_stable), stable_pct))
    output <- c(output, sprintf("• %s (%.1f%%) of farm area is declining week-over-week", format_area_both(total_declining), declining_pct))
    output <- c(output, sprintf("• Total farm area analyzed: %s", format_area_both(total_hectares)))

    if (improving_pct > declining_pct) {
      trend_diff <- round(improving_pct - declining_pct, 1)
      output <- c(output, sprintf("• Overall trend: POSITIVE (%.1f%% more area improving than declining)", trend_diff))
    } else if (declining_pct > improving_pct) {
      trend_diff <- round(declining_pct - improving_pct, 1)
      output <- c(output, sprintf("• Overall trend: NEGATIVE (%.1f%% more area declining than improving)", trend_diff))
    } else {
      output <- c(output, "• Overall trend: BALANCED (equal improvement and decline)")
    }

    # Add note about 0% decline potentially being due to missing data
    if (declining_pct == 0) {
      output <- c(output, "")
      output <- c(output, "⚠️ IMPORTANT NOTE: 0% decline does NOT necessarily mean all crops are healthy.")
      output <- c(output, "• This may be due to missing satellite data from the previous week (cloud cover)")
      output <- c(output, "• Areas with clouds (CI=0) cannot be analyzed for decline")
      output <- c(output, "• True decline levels may be higher than reported")
    }
  } else {
    output <- c(output, "• Analysis: Spatial patterns only (previous week data unavailable)")
    output <- c(output, "• Total farm area analyzed: %.1f hectares", total_hectares)
    output <- c(output, "• Note: Due to clouds in previous week (CI=0), no decline measurements available")
    output <- c(output, "• This does NOT mean fields didn't decline - only that no comparison data exists")
  }

  # Add KPI Dashboard Tables
  output <- c(output, "")
  output <- c(output, "=== FARM KEY PERFORMANCE INDICATORS ===")
  output <- c(output, "")
  
  # Table 1: Field Performance Distribution & Risk Assessment
  output <- c(output, "FIELD PERFORMANCE INDICATORS")
  output <- c(output, "________________________________________")
  output <- c(output, "UNIFORMITY DISTRIBUTION:                    RISK ASSESSMENT:")
  output <- c(output, "CV Category         Count   Percent         Risk Level        Count   Percent")
  
  # Calculate risk levels based on CV + Moran's I combination
  risk_low <- 0
  risk_moderate <- 0  
  risk_high <- 0
  risk_very_high <- 0
  
  for (field_id in names(field_results)) {
    field_info <- field_results[[field_id]]
    cv <- field_info$current_stats$cv
    morans_i <- field_info$current_stats$spatial_autocorr$morans_i
    
    # Risk logic: Low CV + Low clustering = Low risk, High CV + High clustering = High risk
    if (!is.na(cv) && !is.na(morans_i)) {
      if (cv <= 0.10 && morans_i <= 0.8) {
        risk_low <- risk_low + 1
      } else if (cv <= 0.20 && morans_i <= 0.9) {
        risk_moderate <- risk_moderate + 1
      } else if (cv <= 0.30 || morans_i <= 0.95) {
        risk_high <- risk_high + 1
      } else {
        risk_very_high <- risk_very_high + 1
      }
    } else {
      risk_moderate <- risk_moderate + 1  # Default for missing data
    }
  }
  
  output <- c(output, sprintf("Excellent (CV≤0.08)   %d     %5.1f%%          Low (CV≤0.10)       %d     %5.1f%%", 
                             excellent_fields, (excellent_fields/total_fields)*100, risk_low, (risk_low/total_fields)*100))
  output <- c(output, sprintf("Good (CV 0.08-0.15)   %d     %5.1f%%          Moderate (0.10-0.20) %d     %5.1f%%", 
                             good_fields, (good_fields/total_fields)*100, risk_moderate, (risk_moderate/total_fields)*100))
  output <- c(output, sprintf("Moderate (0.15-0.30)  %d     %5.1f%%          High (0.20-0.30)    %d     %5.1f%%", 
                             moderate_fields, (moderate_fields/total_fields)*100, risk_high, (risk_high/total_fields)*100))
  output <- c(output, sprintf("Poor (CV>0.30)        %d     %5.1f%%          Very High (>0.30)   %d     %5.1f%%", 
                             poor_fields, (poor_fields/total_fields)*100, risk_very_high, (risk_very_high/total_fields)*100))
  output <- c(output, sprintf("Total fields         %d    100.0%%          Total fields       %d    100.0%%", 
                             total_fields, total_fields))
  output <- c(output, "")
  
  # Performance quartiles and CI change patterns
  if (has_temporal_data) {
    # Calculate performance quartiles based on combination of current CI and change
    field_performance <- sapply(field_results, function(x) {
      current_ci <- x$current_stats$mean_ci
      ci_change <- x$ci_change
      # Combine current performance with improvement trend
      performance_score <- current_ci + (ci_change * 0.5)  # Weight change as 50% of current
      return(performance_score)
    })
    
    sorted_performance <- sort(field_performance, decreasing = TRUE)
    q75 <- quantile(sorted_performance, 0.75, na.rm = TRUE)
    q25 <- quantile(sorted_performance, 0.25, na.rm = TRUE)
    
    top_quartile <- sum(field_performance >= q75, na.rm = TRUE)
    bottom_quartile <- sum(field_performance <= q25, na.rm = TRUE)
    middle_quartile <- total_fields - top_quartile - bottom_quartile
    
    avg_ci_top <- mean(sapply(field_results[field_performance >= q75], function(x) x$current_stats$mean_ci), na.rm = TRUE)
    avg_ci_mid <- mean(sapply(field_results[field_performance > q25 & field_performance < q75], function(x) x$current_stats$mean_ci), na.rm = TRUE)
    avg_ci_bot <- mean(sapply(field_results[field_performance <= q25], function(x) x$current_stats$mean_ci), na.rm = TRUE)
    
    output <- c(output, "PERFORMANCE QUARTILES:                      CI CHANGE PATTERNS:")
    output <- c(output, "Quartile            Count   Avg CI          Change Type       Hectares  Percent")
    output <- c(output, sprintf("Top 25%%              %d      %4.1f            Improving areas     %5.1f ha   %5.1f%%", 
                               top_quartile, avg_ci_top, total_improving, improving_pct))
    output <- c(output, sprintf("Average (25-75%%)     %d      %4.1f            Stable areas       %5.1f ha   %5.1f%%", 
                               middle_quartile, avg_ci_mid, total_stable, stable_pct))
    output <- c(output, sprintf("Bottom 25%%           %d      %4.1f            Declining areas     %5.1f ha   %5.1f%%", 
                               bottom_quartile, avg_ci_bot, total_declining, declining_pct))
    output <- c(output, sprintf("Total fields        %d      %4.1f            Total area         %5.1f ha  100.0%%", 
                               total_fields, mean(sapply(field_results, function(x) x$current_stats$mean_ci), na.rm = TRUE), total_hectares))
  }
  output <- c(output, "")
  
  # Table 2: Anomaly Detection & Management Alerts
  output <- c(output, "ANOMALY DETECTION & MANAGEMENT PRIORITIES")
  output <- c(output, "________________________________________")
  
  # Weed detection (CI increase > 1.5)
  weed_fields <- 0
  weed_area <- 0
  harvest_fields <- 0
  harvest_area <- 0
  fallow_fields <- 0
  fallow_area <- 0
  high_hotspot_fields <- 0
  high_hotspot_area <- 0
  
  if (has_temporal_data) {
    for (field_id in names(field_results)) {
      field_info <- field_results[[field_id]]
      ci_change <- field_info$ci_change
      current_ci <- field_info$current_stats$mean_ci
      area <- field_info$current_stats$field_area_ha
      hotspots <- field_info$current_stats$extreme_percentages$hotspot_pct
      
      # Weed detection: CI increase > 1.5
      if (!is.na(ci_change) && ci_change > 1.5) {
        weed_fields <- weed_fields + 1
        weed_area <- weed_area + area
      }
      
      # Harvesting/theft detection: CI decrease > 1.5
      if (!is.na(ci_change) && ci_change < -1.5) {
        harvest_fields <- harvest_fields + 1
        harvest_area <- harvest_area + area
      }
      
      # Fallow detection: CI < 2.0
      if (!is.na(current_ci) && current_ci < 2.0) {
        fallow_fields <- fallow_fields + 1
        fallow_area <- fallow_area + area
      }
      
      # High hotspot detection: > 5%
      if (!is.na(hotspots) && hotspots > 5.0) {
        high_hotspot_fields <- high_hotspot_fields + 1
        high_hotspot_area <- high_hotspot_area + area
      }
    }
  }
  
  output <- c(output, "WEED PRESENCE INDICATORS:                   HARVESTING/THEFT INDICATORS:")
  output <- c(output, "High CI Increase (>1.5):                   High CI Decrease (>1.5):")
  output <- c(output, "Fields to check       Count   Area         Fields to check      Count   Area")
  output <- c(output, sprintf("Potential weed areas    %d     %4.1f ha       Potential harvesting   %d     %4.1f ha", 
                             weed_fields, weed_area, harvest_fields, harvest_area))
  output <- c(output, sprintf("Total monitored fields %d    %5.1f ha       Total monitored fields%d    %5.1f ha", 
                             total_fields, total_hectares, total_fields, total_hectares))
  output <- c(output, "")
  output <- c(output, "FALLOW FIELD DETECTION:                     HOTSPOT ANALYSIS:")
  output <- c(output, "Low CI Fields (<2.0):                       Spatial Clustering:")
  output <- c(output, "Fields to check       Count   Area         High hotspot fields (>5%) Count   Area")
  output <- c(output, sprintf("Potential fallow        %d     %4.1f ha       Spatial issues detected   %d     %4.1f ha", 
                             fallow_fields, fallow_area, high_hotspot_fields, high_hotspot_area))
  output <- c(output, sprintf("Total catchment fields %d    %5.1f ha       Total analyzed fields    %d    %5.1f ha", 
                             total_fields, total_hectares, total_fields, total_hectares))
  output <- c(output, "")
  
  # Table 3: Priority Action Items & Field Rankings
  output <- c(output, "IMMEDIATE ACTION PRIORITIES")
  output <- c(output, "________________________________________")
  
  # Find urgent and monitoring fields
  urgent_fields <- sapply(field_results, function(x) x$message_result$worth_sending && grepl("URGENT", x$message_result$message))
  monitoring_fields <- sapply(field_results, function(x) x$message_result$worth_sending && !grepl("URGENT", x$message_result$message))
  
  output <- c(output, "URGENT INTERVENTIONS:                       MONITORING REQUIRED:")
  output <- c(output, "Field Name          Issue Type    Area      Field Name         Issue Type     Area")
  
  urgent_count <- 0
  monitoring_count <- 0
  
  for (field_id in names(field_results)) {
    if (urgent_fields[field_id]) {
      field_info <- field_results[[field_id]]
      field_name <- paste(field_info$current_stats$field, field_info$current_stats$sub_field, sep="-")
      if (nchar(field_name) > 15) field_name <- substr(field_name, 1, 15)
      area <- field_info$current_stats$field_area_ha
      if (urgent_count == 0) {
        output <- c(output, sprintf("%-15s Poor uniformity %4.1f ha   %-15s %-13s %4.1f ha", 
                                   field_name, area, "", "", 0.0))
      }
      urgent_count <- urgent_count + 1
    }
    
    if (monitoring_fields[field_id]) {
      field_info <- field_results[[field_id]]
      field_name <- paste(field_info$current_stats$field, field_info$current_stats$sub_field, sep="-")
      if (nchar(field_name) > 15) field_name <- substr(field_name, 1, 15)
      area <- field_info$current_stats$field_area_ha
      if (monitoring_count == 0) {
        # Update the previous line to include monitoring field
        last_line <- output[length(output)]
        if (grepl("Poor uniformity", last_line) && grepl("0.0 ha$", last_line)) {
          output[length(output)] <- sprintf("%-15s Poor uniformity %4.1f ha   %-15s %-13s %4.1f ha", 
                                           sub(" .*", "", last_line), 
                                           as.numeric(sub(".*Poor uniformity ([0-9.]+) ha.*", "\\1", last_line)),
                                           field_name, "Moderate var.", area)
        }
      }
      monitoring_count <- monitoring_count + 1
    }
  }
  
  if (urgent_count == 0 && monitoring_count == 0) {
    output <- c(output, "No urgent interventions    -         -      No monitoring required    -         -")
  }
  
  output <- c(output, "")
  
  # Field performance ranking
  if (has_temporal_data) {
    output <- c(output, "FIELD PERFORMANCE RANKING:                  WEEKLY PRIORITIES:")
    output <- c(output, "Rank  Field Name        CI    Status        Priority Level    Fields  Action Required")
    
    # Sort fields by performance score
    field_names <- names(field_performance)
    sorted_indices <- order(field_performance, decreasing = TRUE)
    
    priority_immediate <- sum(urgent_fields)
    priority_weekly <- sum(monitoring_fields)
    priority_routine <- total_fields - priority_immediate - priority_weekly
    
    for (i in 1:min(3, length(sorted_indices))) {
      field_id <- field_names[sorted_indices[i]]
      field_info <- field_results[[field_id]]
      field_name <- paste(field_info$current_stats$field, field_info$current_stats$sub_field, sep="-")
      if (nchar(field_name) > 12) field_name <- substr(field_name, 1, 12)
      
      ci <- field_info$current_stats$mean_ci
      status <- if (field_info$current_stats$cv <= 0.08) "Excellent" else if (field_info$current_stats$cv <= 0.15) "Good" else "Caution"
      
      if (i == 1) {
        output <- c(output, sprintf("%d     %-12s %4.1f   %-9s     Immediate           %d      Field inspection", 
                                   i, field_name, ci, status, priority_immediate))
      } else if (i == 2) {
        output <- c(output, sprintf("%d     %-12s %4.1f   %-9s     This week           %d      Continue monitoring", 
                                   i, field_name, ci, status, priority_weekly))
      } else {
        output <- c(output, sprintf("%d     %-12s %4.1f   %-9s     Monitor             %d      Routine management", 
                                   i, field_name, ci, status, priority_routine))
      }
    }
    
    output <- c(output, sprintf("...                                           Total fields       %d", total_fields))
  }

  return(paste(output, collapse = "\n"))
}

#' Format analysis results as CSV data
#' @param analysis_results Results from run_estate_analysis
#' @return Data frame ready for write.csv
format_as_csv <- function(analysis_results) {

  field_results <- analysis_results$field_results
  estate_name <- analysis_results$estate_name
  current_week <- analysis_results$current_week
  previous_week <- analysis_results$previous_week

  csv_data <- data.frame()

  for (field_id in names(field_results)) {
    field_info <- field_results[[field_id]]

    row_data <- data.frame(
      estate = estate_name,
      field = field_info$current_stats$field,
      sub_field = field_info$current_stats$sub_field,
      area_ha = round(field_info$current_stats$field_area_ha, 2),
      current_week = current_week,
      previous_week = previous_week,
      current_week_ci = round(field_info$current_stats$mean_ci, 3),
      previous_week_ci = if (!is.null(field_info$previous_stats)) round(field_info$previous_stats$mean_ci, 3) else NA,
      ci_change = round(field_info$ci_change, 3),
      change_category = field_info$change_category,
      cv = round(field_info$current_stats$cv, 3),
      uniformity_category = field_info$uniformity_category,
      acceptable_pct = round(field_info$current_stats$acceptable_pct, 1),
      hotspot_pct = round(field_info$current_stats$extreme_percentages$hotspot_pct, 1),
      coldspot_pct = round(field_info$current_stats$extreme_percentages$coldspot_pct, 1),
      morans_i = round(field_info$current_stats$spatial_autocorr$morans_i, 3),
      alert_needed = field_info$message_result$worth_sending,
      message = field_info$message_result$message,
      stringsAsFactors = FALSE
    )

    csv_data <- rbind(csv_data, row_data)
  }

  return(csv_data)
}

#' Format analysis results as markdown table
#' @param analysis_results Results from run_estate_analysis
#' @return Character string with markdown table
format_as_markdown_table <- function(analysis_results) {

  field_results <- analysis_results$field_results
  estate_name <- toupper(analysis_results$estate_name)
  current_week <- analysis_results$current_week
  previous_week <- analysis_results$previous_week

  output <- c()
  output <- c(output, paste("# Crop Analysis Summary -", estate_name, "Estate"))
  output <- c(output, paste("**Analysis Period:** Week", previous_week, "vs Week", current_week))
  output <- c(output, "")
  output <- c(output, "| Field | Area (ha) | Current CI | Change | Uniformity | Alert | Message |")
  output <- c(output, "|-------|-----------|------------|--------|------------|-------|---------|")

  for (field_id in names(field_results)) {
    field_info <- field_results[[field_id]]

    field_name <- paste(field_info$current_stats$field, field_info$current_stats$sub_field, sep="-")
    area <- round(field_info$current_stats$field_area_ha, 1)
    current_ci <- round(field_info$current_stats$mean_ci, 3)
    change <- field_info$change_category
    uniformity <- field_info$uniformity_category
    alert <- if(field_info$message_result$worth_sending) "🚨 YES" else "✅ NO"
    message <- field_info$message_result$message

    row <- paste("|", field_name, "|", area, "|", current_ci, "|", change, "|", uniformity, "|", alert, "|", message, "|")
    output <- c(output, row)
  }

  return(paste(output, collapse = "\n"))
}

#' Create Word document with analysis results
#' @param analysis_results Results from run_estate_analysis
#' @param output_dir Directory to save the Word document
#' @return Path to the created Word document
create_word_document <- function(analysis_results, output_dir) {

  estate_name <- toupper(analysis_results$estate_name)
  current_week <- analysis_results$current_week
  previous_week <- analysis_results$previous_week

  # Create a new Word document
  doc <- officer::read_docx()

  # Add title
  doc <- officer::body_add_par(doc, paste(estate_name, "Crop Analysis Report"), style = "heading 1")

  # Add summary
  field_results <- analysis_results$field_results
  alert_count <- sum(sapply(field_results, function(x) x$message_result$worth_sending))
  total_fields <- length(field_results)
  total_hectares <- sum(sapply(field_results, function(x) x$current_stats$field_area_ha), na.rm = TRUE)

  doc <- officer::body_add_par(doc, "Summary", style = "heading 2")
  doc <- officer::body_add_par(doc, paste("• Fields analyzed:", total_fields))
  doc <- officer::body_add_par(doc, paste("• Total area:", format_area_both(total_hectares)))
  doc <- officer::body_add_par(doc, paste("• Alerts needed:", alert_count))
  doc <- officer::body_add_par(doc, "")

  # Field-by-field alerts only
  if (alert_count > 0) {
    doc <- officer::body_add_par(doc, "Priority Fields", style = "heading 2")
    for (field_id in names(field_results)) {
      field_info <- field_results[[field_id]]
      if (field_info$message_result$worth_sending) {
        field_name <- paste(field_info$current_stats$field, field_info$current_stats$sub_field, sep="-")
        area <- round(field_info$current_stats$field_area_ha, 1)
        message <- field_info$message_result$message

        doc <- officer::body_add_par(doc, paste("•", field_name, paste0("(", format_area_both(area), "):"), message))
      }
    }
    doc <- officer::body_add_par(doc, "")
  } else {
    doc <- officer::body_add_par(doc, "✅ No urgent alerts - all fields stable")
    doc <- officer::body_add_par(doc, "")
  }

  # Quick farm summary
  doc <- officer::body_add_par(doc, "Quick Stats", style = "heading 2")

  # Calculate improving vs declining areas (only if temporal data available)
  has_temporal_data <- any(sapply(field_results, function(x) !is.na(x$change_percentages$positive_pct)))

  if (has_temporal_data) {
    total_improving <- sum(sapply(field_results, function(x) {
      if (!is.na(x$change_percentages$positive_pct)) {
        (x$change_percentages$positive_pct / 100) * x$current_stats$field_area_ha
      } else 0
    }), na.rm = TRUE)

    total_declining <- sum(sapply(field_results, function(x) {
      if (!is.na(x$change_percentages$negative_pct)) {
        (x$change_percentages$negative_pct / 100) * x$current_stats$field_area_ha
      } else 0
    }), na.rm = TRUE)

    total_stable <- sum(sapply(field_results, function(x) {
      if (!is.na(x$change_percentages$stable_pct)) {
        (x$change_percentages$stable_pct / 100) * x$current_stats$field_area_ha
      } else 0
    }), na.rm = TRUE)

    improving_pct <- (total_improving / total_hectares) * 100
    declining_pct <- (total_declining / total_hectares) * 100
    stable_pct <- (total_stable / total_hectares) * 100

    doc <- officer::body_add_par(doc, paste("• Improving areas:", format_area_both(total_improving), paste0("(", round(improving_pct, 1), "%)")))
    doc <- officer::body_add_par(doc, paste("• Stable areas:", format_area_both(total_stable), paste0("(", round(stable_pct, 1), "%)")))
    doc <- officer::body_add_par(doc, paste("• Declining areas:", format_area_both(total_declining), paste0("(", round(declining_pct, 1), "%)")))

    # Overall trend
    if (improving_pct > declining_pct) {
      trend_diff <- round(improving_pct - declining_pct, 1)
      doc <- officer::body_add_par(doc, paste("• Trend: POSITIVE (+", trend_diff, "%)"))
    } else if (declining_pct > improving_pct) {
      trend_diff <- round(declining_pct - improving_pct, 1)
      doc <- officer::body_add_par(doc, paste("• Trend: NEGATIVE (-", trend_diff, "%)"))
    } else {
      doc <- officer::body_add_par(doc, "• Trend: BALANCED")
    }
  } else {
    doc <- officer::body_add_par(doc, "• Analysis: Spatial patterns only (previous week data unavailable)")
  }

  doc <- officer::body_add_par(doc, "")

  # Add farm-wide analysis summary
  doc <- officer::body_add_par(doc, "Farm-Wide Analysis Summary", style = "heading 2")
  doc <- officer::body_add_par(doc, "")

  # Field uniformity statistics
  excellent_fields <- sum(sapply(field_results, function(x) x$current_stats$cv <= 0.08))
  good_fields <- sum(sapply(field_results, function(x) x$current_stats$cv > 0.08 & x$current_stats$cv <= 0.15))
  moderate_fields <- sum(sapply(field_results, function(x) x$current_stats$cv > 0.15 & x$current_stats$cv <= 0.30))
  poor_fields <- sum(sapply(field_results, function(x) x$current_stats$cv > 0.30))

  # Create uniformity table
  uniformity_data <- data.frame(
    "Uniformity Level" = c("Excellent (CV≤0.08)", "Good (CV 0.08-0.15)", "Moderate (CV 0.15-0.30)", "Poor (CV>0.30)", "Total fields"),
    "Fields" = c(excellent_fields, good_fields, moderate_fields, poor_fields, total_fields),
    "Percent" = c(
      round((excellent_fields/total_fields)*100, 1),
      round((good_fields/total_fields)*100, 1),
      round((moderate_fields/total_fields)*100, 1),
      round((poor_fields/total_fields)*100, 1),
      100.0
    ),
    stringsAsFactors = FALSE
  )

  uniformity_ft <- flextable::flextable(uniformity_data)
  uniformity_ft <- flextable::autofit(uniformity_ft)
  uniformity_ft <- flextable::set_header_labels(uniformity_ft,
    "Uniformity.Level" = "Uniformity Level"
  )

  doc <- officer::body_add_par(doc, "Field Uniformity Summary", style = "heading 3")
  doc <- flextable::body_add_flextable(doc, uniformity_ft)
  doc <- officer::body_add_par(doc, "")

  # Farm-wide area change summary
  doc <- officer::body_add_par(doc, "Farm-Wide Area Change Summary", style = "heading 3")

  if (has_temporal_data) {
    change_data <- data.frame(
      "Change Type" = c("Improving areas", "Stable areas", "Declining areas", "Total area"),
      "Area" = c(format_area_both(total_improving), format_area_both(total_stable), format_area_both(total_declining), format_area_both(total_hectares)),
      "Percent" = c(round(improving_pct, 1), round(stable_pct, 1), round(declining_pct, 1), 100.0),
      stringsAsFactors = FALSE
    )
  } else {
    change_data <- data.frame(
      "Change Type" = c("Improving areas", "Stable areas", "Declining areas", "Total area"),
      "Area" = c("N/A", "N/A", "N/A", format_area_both(total_hectares)),
      "Percent" = c("N/A", "N/A", "N/A", 100.0),
      stringsAsFactors = FALSE
    )
  }

  change_ft <- flextable::flextable(change_data)
  change_ft <- flextable::autofit(change_ft)
  change_ft <- flextable::set_header_labels(change_ft,
    "Change.Type" = "Change Type",
    "Area" = "Area (ha/acres)"
  )

  doc <- flextable::body_add_flextable(doc, change_ft)
  doc <- officer::body_add_par(doc, "")

  # Create and add detailed results tables using flextable (split into multiple tables for better formatting)
  csv_data <- format_as_csv(analysis_results)

  # Split data into multiple tables for better readability
  doc <- officer::body_add_par(doc, "Detailed Results", style = "heading 2")

  # Table 2: Current Week Analysis
  current_data <- csv_data[, c("field", "current_week_ci", "cv", "uniformity_category", "acceptable_pct")]
  current_data$current_week_ci <- round(current_data$current_week_ci, 3)
  current_data$cv <- round(current_data$cv, 3)
  current_data$acceptable_pct <- round(current_data$acceptable_pct, 1)

  current_ft <- flextable::flextable(current_data)
  current_ft <- flextable::autofit(current_ft)
  current_ft <- flextable::set_header_labels(current_ft,
    "field" = "Field",
  #  "sub_field" = "Sub-field",
    "current_week_ci" = "Current CI",
    "cv" = "CV",
    "uniformity_category" = "Uniformity",
    "acceptable_pct" = "Acceptable %"
  )
  current_ft <- flextable::theme_vanilla(current_ft)
  current_ft <- flextable::fontsize(current_ft, size = 9)
  current_ft <- flextable::width(current_ft, width = 1.0)  # Set column width

  doc <- officer::body_add_par(doc, "Current Week Analysis", style = "heading 3")
  doc <- flextable::body_add_flextable(doc, current_ft)
  doc <- officer::body_add_par(doc, "")

  # Table 3: Change Analysis (only if temporal data available)
  if (has_temporal_data && any(!is.na(csv_data$ci_change))) {
    change_data <- csv_data[, c("field", "previous_week_ci", "ci_change", "change_category")]
    change_data <- change_data[!is.na(change_data$ci_change), ]  # Remove rows with NA change
    change_data$previous_week_ci <- round(change_data$previous_week_ci, 3)
    change_data$ci_change <- round(change_data$ci_change, 3)

    change_ft <- flextable::flextable(change_data)
    change_ft <- flextable::autofit(change_ft)
    change_ft <- flextable::set_header_labels(change_ft,
      "field" = "Field",
    #  "sub_field" = "Sub-field",
      "previous_week_ci" = "Previous CI",
      "ci_change" = "CI Change",
      "change_category" = "Change Type"
    )
    change_ft <- flextable::theme_vanilla(change_ft)
    change_ft <- flextable::fontsize(change_ft, size = 9)
    change_ft <- flextable::width(change_ft, width = 1.0)  # Set column width

    doc <- officer::body_add_par(doc, "Week-over-Week Change Analysis", style = "heading 3")
    doc <- flextable::body_add_flextable(doc, change_ft)
    doc <- officer::body_add_par(doc, "")
  }

  # Table 4: Spatial Analysis Results (split into two tables for better fit)
  spatial_data <- csv_data[, c("field", "hotspot_pct", "coldspot_pct", "morans_i", "alert_needed")]
  spatial_data$hotspot_pct <- round(spatial_data$hotspot_pct, 1)
  spatial_data$coldspot_pct <- round(spatial_data$coldspot_pct, 1)
  spatial_data$morans_i <- round(spatial_data$morans_i, 3)
  spatial_data$alert_needed <- ifelse(spatial_data$alert_needed, "YES", "NO")

  spatial_ft <- flextable::flextable(spatial_data)
  spatial_ft <- flextable::autofit(spatial_ft)
  spatial_ft <- flextable::set_header_labels(spatial_ft,
    "field" = "Field",
   # "sub_field" = "Sub-field",
    "hotspot_pct" = "Hotspots %",
    "coldspot_pct" = "Coldspots %",
    "morans_i" = "Moran's I",
    "alert_needed" = "Alert"
  )
  spatial_ft <- flextable::theme_vanilla(spatial_ft)
  spatial_ft <- flextable::fontsize(spatial_ft, size = 9)
  spatial_ft <- flextable::width(spatial_ft, width = 0.8)  # Set column width

  doc <- officer::body_add_par(doc, "Spatial Analysis Results", style = "heading 3")
  doc <- flextable::body_add_flextable(doc, spatial_ft)
  doc <- officer::body_add_par(doc, "")

  # Table 5: Alert Messages (separate table for long messages)
  message_data <- csv_data[, c("field","message")]
  message_data$message <- substr(message_data$message, 1, 80)  # Truncate long messages for table fit

  message_ft <- flextable::flextable(message_data)
  message_ft <- flextable::autofit(message_ft)
  message_ft <- flextable::set_header_labels(message_ft,
    "field" = "Field",
  #  "sub_field" = "Sub-field",
    "message" = "Alert Message"
  )
  message_ft <- flextable::theme_vanilla(message_ft)
  message_ft <- flextable::fontsize(message_ft, size = 8)  # Smaller font for messages
  message_ft <- flextable::width(message_ft, width = 2.0)  # Wider column for messages

  doc <- officer::body_add_par(doc, "Alert Messages", style = "heading 3")
  doc <- flextable::body_add_flextable(doc, message_ft)

  # Add interpretation guide for all columns
  doc <- officer::body_add_par(doc, "")
  doc <- officer::body_add_par(doc, "Column Interpretation Guide", style = "heading 3")
  doc <- officer::body_add_par(doc, "")

  # Table 1 interpretation
  doc <- officer::body_add_par(doc, "Field Information Table:", style = "Normal")
  doc <- officer::body_add_par(doc, "• Field/Sub-field: Field identifiers and names")
  doc <- officer::body_add_par(doc, "• Area (ha): Field size in hectares")
  doc <- officer::body_add_par(doc, "• Current/Previous Week: Weeks being compared")
  doc <- officer::body_add_par(doc, "")

  # Table 2 interpretation
  doc <- officer::body_add_par(doc, "Current Week Analysis Table:", style = "Normal")
  doc <- officer::body_add_par(doc, "• Current CI: Crop Index (0-10 scale, higher = healthier crop)")
  doc <- officer::body_add_par(doc, "• CV: Coefficient of Variation (lower = more uniform field)")
  doc <- officer::body_add_par(doc, "• Uniformity: Field uniformity rating (Excellent/Good/Moderate/Poor)")
  doc <- officer::body_add_par(doc, "• Acceptable %: % of field within ±25% of average CI (higher = more uniform)")
  doc <- officer::body_add_par(doc, "")

  # Table 3 interpretation (only if temporal data available)
  if (has_temporal_data && any(!is.na(csv_data$ci_change))) {
    doc <- officer::body_add_par(doc, "Week-over-Week Change Analysis Table:", style = "Normal")
    doc <- officer::body_add_par(doc, "• Previous CI: Crop Index from previous week")
    doc <- officer::body_add_par(doc, "• CI Change: Week-over-week change in CI values")
    doc <- officer::body_add_par(doc, "• Change Type: >+0.5 = Improving, -0.5 to +0.5 = Stable, <-0.5 = Declining")
    doc <- officer::body_add_par(doc, "")
  }

  # Table 4 interpretation
  doc <- officer::body_add_par(doc, "Spatial Analysis Results Table:", style = "Normal")
  doc <- officer::body_add_par(doc, "• Hotspots %: % of field significantly above average (> mean + 1.5×SD)")
  doc <- officer::body_add_par(doc, "• Coldspots %: % of field significantly below average (< mean - 1.5×SD)")
  doc <- officer::body_add_par(doc, "• Moran's I: Spatial autocorrelation (-1 to +1, higher = more clustered)")
  doc <- officer::body_add_par(doc, "• Alert: YES/NO indicating if field needs management attention")
  doc <- officer::body_add_par(doc, "")

  # Table 5 interpretation
  doc <- officer::body_add_par(doc, "Alert Messages Table:", style = "Normal")
  doc <- officer::body_add_par(doc, "• Message: Specific recommendations or warnings for each field")
  doc <- officer::body_add_par(doc, "")

  # Overall interpretation guide
  doc <- officer::body_add_par(doc, "Performance Thresholds:", style = "heading 3")
  doc <- officer::body_add_par(doc, "Acceptable %: >45% = Excellent uniformity, 35-45% = Good, <35% = Needs attention")
  doc <- officer::body_add_par(doc, "CV: <0.08 = Excellent, 0.08-0.15 = Good, 0.15-0.30 = Moderate, >0.30 = Poor")
  doc <- officer::body_add_par(doc, "Moran's I: >0.7 = Strong clustering, 0.3-0.7 = Normal field patterns, <0.3 = Random")
  doc <- officer::body_add_par(doc, "Hotspots/Coldspots: >10% = Significant spatial issues, 3-10% = Monitor, <3% = Normal")

  # Add KPI Dashboard to Word Document
  doc <- officer::body_add_par(doc, "")
  doc <- officer::body_add_par(doc, "Farm Key Performance Indicators", style = "heading 2")
  doc <- officer::body_add_par(doc, "")
  
  # Table 1: Field Performance Distribution & Risk Assessment
  doc <- officer::body_add_par(doc, "Field Performance Indicators", style = "heading 3")
  
  # Calculate risk levels based on CV + Moran's I combination
  risk_low <- 0
  risk_moderate <- 0  
  risk_high <- 0
  risk_very_high <- 0
  
  for (field_id in names(field_results)) {
    field_info <- field_results[[field_id]]
    cv <- field_info$current_stats$cv
    morans_i <- field_info$current_stats$spatial_autocorr$morans_i
    
    # Risk logic: Low CV + Low clustering = Low risk, High CV + High clustering = High risk
    if (!is.na(cv) && !is.na(morans_i)) {
      if (cv <= 0.10 && morans_i <= 0.8) {
        risk_low <- risk_low + 1
      } else if (cv <= 0.20 && morans_i <= 0.9) {
        risk_moderate <- risk_moderate + 1
      } else if (cv <= 0.30 || morans_i <= 0.95) {
        risk_high <- risk_high + 1
      } else {
        risk_very_high <- risk_very_high + 1
      }
    } else {
      risk_moderate <- risk_moderate + 1  # Default for missing data
    }
  }
  
  # Uniformity Distribution Table
  uniformity_kpi_data <- data.frame(
    "CV Category" = c("Excellent (CV≤0.08)", "Good (CV 0.08-0.15)", "Moderate (0.15-0.30)", "Poor (CV>0.30)", "Total fields"),
    "Count" = c(excellent_fields, good_fields, moderate_fields, poor_fields, total_fields),
    "Percent" = c(round((excellent_fields/total_fields)*100, 1), round((good_fields/total_fields)*100, 1), 
                  round((moderate_fields/total_fields)*100, 1), round((poor_fields/total_fields)*100, 1), 100.0),
    stringsAsFactors = FALSE
  )
  
  # Risk Assessment Table
  risk_data <- data.frame(
    "Risk Level" = c("Low (CV≤0.10)", "Moderate (0.10-0.20)", "High (0.20-0.30)", "Very High (>0.30)", "Total fields"),
    "Count" = c(risk_low, risk_moderate, risk_high, risk_very_high, total_fields),
    "Percent" = c(round((risk_low/total_fields)*100, 1), round((risk_moderate/total_fields)*100, 1), 
                  round((risk_high/total_fields)*100, 1), round((risk_very_high/total_fields)*100, 1), 100.0),
    stringsAsFactors = FALSE
  )
  
  uniformity_kpi_ft <- flextable::flextable(uniformity_kpi_data)
  uniformity_kpi_ft <- flextable::autofit(uniformity_kpi_ft)
  risk_ft <- flextable::flextable(risk_data)
  risk_ft <- flextable::autofit(risk_ft)
  
  doc <- officer::body_add_par(doc, "Uniformity Distribution:")
  doc <- flextable::body_add_flextable(doc, uniformity_kpi_ft)
  doc <- officer::body_add_par(doc, "")
  doc <- officer::body_add_par(doc, "Risk Assessment:")
  doc <- flextable::body_add_flextable(doc, risk_ft)
  doc <- officer::body_add_par(doc, "")
  
  # Performance Quartiles (if temporal data available)
  if (has_temporal_data) {
    # Calculate performance quartiles based on combination of current CI and change
    field_performance <- sapply(field_results, function(x) {
      current_ci <- x$current_stats$mean_ci
      ci_change <- x$ci_change
      # Combine current performance with improvement trend
      performance_score <- current_ci + (ci_change * 0.5)  # Weight change as 50% of current
      return(performance_score)
    })
    
    sorted_performance <- sort(field_performance, decreasing = TRUE)
    q75 <- quantile(sorted_performance, 0.75, na.rm = TRUE)
    q25 <- quantile(sorted_performance, 0.25, na.rm = TRUE)
    
    top_quartile <- sum(field_performance >= q75, na.rm = TRUE)
    bottom_quartile <- sum(field_performance <= q25, na.rm = TRUE)
    middle_quartile <- total_fields - top_quartile - bottom_quartile
    
    avg_ci_top <- mean(sapply(field_results[field_performance >= q75], function(x) x$current_stats$mean_ci), na.rm = TRUE)
    avg_ci_mid <- mean(sapply(field_results[field_performance > q25 & field_performance < q75], function(x) x$current_stats$mean_ci), na.rm = TRUE)
    avg_ci_bot <- mean(sapply(field_results[field_performance <= q25], function(x) x$current_stats$mean_ci), na.rm = TRUE)
    
    quartile_data <- data.frame(
      "Quartile" = c("Top 25%", "Average (25-75%)", "Bottom 25%", "Total fields"),
      "Count" = c(top_quartile, middle_quartile, bottom_quartile, total_fields),
      "Avg CI" = c(round(avg_ci_top, 1), round(avg_ci_mid, 1), round(avg_ci_bot, 1), 
                   round(mean(sapply(field_results, function(x) x$current_stats$mean_ci), na.rm = TRUE), 1)),
      stringsAsFactors = FALSE
    )
    
    quartile_ft <- flextable::flextable(quartile_data)
    quartile_ft <- flextable::autofit(quartile_ft)
    
    doc <- officer::body_add_par(doc, "Performance Quartiles:")
    doc <- flextable::body_add_flextable(doc, quartile_ft)
    doc <- officer::body_add_par(doc, "")
  }
  
  # Table 2: Anomaly Detection
  doc <- officer::body_add_par(doc, "Anomaly Detection & Management Priorities", style = "heading 3")
  
  # Calculate anomalies
  weed_fields <- 0
  weed_area <- 0
  harvest_fields <- 0
  harvest_area <- 0
  fallow_fields <- 0
  fallow_area <- 0
  high_hotspot_fields <- 0
  high_hotspot_area <- 0
  
  if (has_temporal_data) {
    for (field_id in names(field_results)) {
      field_info <- field_results[[field_id]]
      ci_change <- field_info$ci_change
      current_ci <- field_info$current_stats$mean_ci
      area <- field_info$current_stats$field_area_ha
      hotspots <- field_info$current_stats$extreme_percentages$hotspot_pct
      
      # Weed detection: CI increase > 1.5
      if (!is.na(ci_change) && ci_change > 1.5) {
        weed_fields <- weed_fields + 1
        weed_area <- weed_area + area
      }
      
      # Harvesting/theft detection: CI decrease > 1.5
      if (!is.na(ci_change) && ci_change < -1.5) {
        harvest_fields <- harvest_fields + 1
        harvest_area <- harvest_area + area
      }
      
      # Fallow detection: CI < 2.0
      if (!is.na(current_ci) && current_ci < 2.0) {
        fallow_fields <- fallow_fields + 1
        fallow_area <- fallow_area + area
      }
      
      # High hotspot detection: > 5%
      if (!is.na(hotspots) && hotspots > 5.0) {
        high_hotspot_fields <- high_hotspot_fields + 1
        high_hotspot_area <- high_hotspot_area + area
      }
    }
  }
  
  anomaly_data <- data.frame(
    "Detection Type" = c("Potential weed areas (CI increase >1.5)", "Potential harvesting (CI decrease >1.5)", 
                         "Potential fallow fields (CI <2.0)", "High hotspot fields (>5%)"),
    "Fields to Check" = c(weed_fields, harvest_fields, fallow_fields, high_hotspot_fields),
    "Area (ha)" = c(round(weed_area, 1), round(harvest_area, 1), round(fallow_area, 1), round(high_hotspot_area, 1)),
    stringsAsFactors = FALSE
  )
  
  anomaly_ft <- flextable::flextable(anomaly_data)
  anomaly_ft <- flextable::autofit(anomaly_ft)
  doc <- flextable::body_add_flextable(doc, anomaly_ft)
  doc <- officer::body_add_par(doc, "")
  
  # Table 3: Priority Actions
  doc <- officer::body_add_par(doc, "Immediate Action Priorities", style = "heading 3")
  
  # Find urgent and monitoring fields
  urgent_fields <- sapply(field_results, function(x) x$message_result$worth_sending && grepl("URGENT", x$message_result$message))
  monitoring_fields <- sapply(field_results, function(x) x$message_result$worth_sending && !grepl("URGENT", x$message_result$message))
  
  urgent_data <- data.frame()
  monitoring_data <- data.frame()
  
  for (field_id in names(field_results)) {
    if (urgent_fields[field_id]) {
      field_info <- field_results[[field_id]]
      field_name <- paste(field_info$current_stats$field, field_info$current_stats$sub_field, sep="-")
      area <- field_info$current_stats$field_area_ha
      urgent_data <- rbind(urgent_data, data.frame(
        "Field Name" = field_name,
        "Issue Type" = "Poor uniformity",
        "Area (ha)" = round(area, 1),
        stringsAsFactors = FALSE
      ))
    }
    
    if (monitoring_fields[field_id]) {
      field_info <- field_results[[field_id]]
      field_name <- paste(field_info$current_stats$field, field_info$current_stats$sub_field, sep="-")
      area <- field_info$current_stats$field_area_ha
      monitoring_data <- rbind(monitoring_data, data.frame(
        "Field Name" = field_name,
        "Issue Type" = "Moderate variation",
        "Area (ha)" = round(area, 1),
        stringsAsFactors = FALSE
      ))
    }
  }
  
  if (nrow(urgent_data) > 0) {
    urgent_ft <- flextable::flextable(urgent_data)
    urgent_ft <- flextable::autofit(urgent_ft)
    doc <- officer::body_add_par(doc, "Urgent Interventions:")
    doc <- flextable::body_add_flextable(doc, urgent_ft)
  } else {
    doc <- officer::body_add_par(doc, "Urgent Interventions: None required")
  }
  
  doc <- officer::body_add_par(doc, "")
  
  if (nrow(monitoring_data) > 0) {
    monitoring_ft <- flextable::flextable(monitoring_data)
    monitoring_ft <- flextable::autofit(monitoring_ft)
    doc <- officer::body_add_par(doc, "Monitoring Required:")
    doc <- flextable::body_add_flextable(doc, monitoring_ft)
  } else {
    doc <- officer::body_add_par(doc, "Monitoring Required: None required")
  }
  
  doc <- officer::body_add_par(doc, "")

  # Add interpretation guide for Table 5
  doc <- officer::body_add_par(doc, "")
  doc <- officer::body_add_par(doc, "Column Guide - Alert Messages:", style = "heading 3")
  doc <- officer::body_add_par(doc, "• Alert Message: Specific recommendations based on field analysis")
  doc <- officer::body_add_par(doc, "• 🚨 URGENT: Immediate management action required")
  doc <- officer::body_add_par(doc, "• ⚠️ ALERT: Early intervention recommended")
  doc <- officer::body_add_par(doc, "• ✅ POSITIVE: Good performance, continue current practices")
  doc <- officer::body_add_par(doc, "• 📈 OPPORTUNITY: Potential for improvement identified")

  # Save the document
  timestamp <- format(Sys.time(), "%Y%m%d_%H%M")
  filename <- paste0("crop_analysis_", estate_name, "_w", current_week, "vs", previous_week, "_", timestamp, ".docx")
  filepath <- file.path(output_dir, filename)

  print(doc, target = filepath)

  return(filepath)
}

#' Save analysis outputs in multiple formats
#' @param analysis_results Results from run_estate_analysis
#' @param output_dir Directory to save files (optional)
#' @return List with file paths created
save_analysis_outputs <- function(analysis_results, output_dir = NULL) {

  estate_name <- analysis_results$estate_name
  current_week <- analysis_results$current_week
  previous_week <- analysis_results$previous_week

  # Create output directory if not specified
  if (is.null(output_dir)) {
    output_dir <- file.path("output", estate_name)
  }
  if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)

  timestamp <- format(Sys.time(), "%Y%m%d_%H%M")
  base_filename <- paste0("crop_analysis_", estate_name, "_w", current_week, "vs", previous_week, "_", timestamp)

  # Generate output formats
  whatsapp_text <- format_for_whatsapp(analysis_results)

  # Save files
  whatsapp_file <- file.path(output_dir, paste0(base_filename, "_whatsapp.txt"))

  writeLines(whatsapp_text, whatsapp_file)

  # Create Word document
  docx_file <- create_word_document(analysis_results, output_dir)

  # Display summary
  cat("\n=== OUTPUT FILES CREATED ===\n")
  cat("📱 WhatsApp format:", whatsapp_file, "\n")
  cat("� Word document:", docx_file, "\n")

  # Display WhatsApp format in console for immediate copy
  cat("\n=== WHATSAPP/WORD READY FORMAT ===\n")
  cat("(Copy text below directly to WhatsApp or Word)\n")
  cat(rep("=", 50), "\n")
  cat(whatsapp_text)
  cat("\n", rep("=", 50), "\n")

  return(list(
    whatsapp_file = whatsapp_file,
    docx_file = docx_file
  ))
}