diff --git a/r_app/80_utils_agronomic_support.R b/r_app/80_utils_agronomic_support.R
index 351d0e5..e609d4a 100644
--- a/r_app/80_utils_agronomic_support.R
+++ b/r_app/80_utils_agronomic_support.R
@@ -338,58 +338,58 @@ calculate_weed_presence_kpi <- function(ci_pixels_by_field) {
   return(result)
 }
 
-#' Calculate Gap Filling Score KPI (placeholder)
-#' @param ci_raster Current week CI raster
-#' @param field_boundaries Field boundaries
-#' @return Data frame with field-level gap filling scores
-calculate_gap_filling_kpi <- function(ci_raster, field_boundaries) {
-  # Handle both sf and SpatVector inputs
-  if (!inherits(field_boundaries, "SpatVector")) {
-    field_boundaries_vect <- terra::vect(field_boundaries)
-  } else {
-    field_boundaries_vect <- field_boundaries
-  }
+# #' Calculate Gap Filling Score KPI (placeholder)
+# #' @param ci_raster Current week CI raster
+# #' @param field_boundaries Field boundaries
+# #' @return Data frame with field-level gap filling scores
+# calculate_gap_filling_kpi <- function(ci_raster, field_boundaries) {
+#   # Handle both sf and SpatVector inputs
+#   if (!inherits(field_boundaries, "SpatVector")) {
+#     field_boundaries_vect <- terra::vect(field_boundaries)
+#   } else {
+#     field_boundaries_vect <- field_boundaries
+#   }
 
-  results_list <- list()
+#   results_list <- list()
 
-  # Ensure field_boundaries_vect is valid and matches field_boundaries dimensions
-  n_fields_vect <- length(field_boundaries_vect)
-  n_fields_sf <- nrow(field_boundaries)
+#   # Ensure field_boundaries_vect is valid and matches field_boundaries dimensions
+#   n_fields_vect <- length(field_boundaries_vect)
+#   n_fields_sf <- nrow(field_boundaries)
   
-  if (n_fields_sf != n_fields_vect) {
-    warning(paste("Field boundary mismatch: nrow(field_boundaries)=", n_fields_sf, "vs length(field_boundaries_vect)=", n_fields_vect, ". Using actual SpatVector length."))
-  }
+#   if (n_fields_sf != n_fields_vect) {
+#     warning(paste("Field boundary mismatch: nrow(field_boundaries)=", n_fields_sf, "vs length(field_boundaries_vect)=", n_fields_vect, ". Using actual SpatVector length."))
+#   }
 
-  for (i in seq_len(n_fields_vect)) {
-    field_vect <- field_boundaries_vect[i]
+#   for (i in seq_len(n_fields_vect)) {
+#     field_vect <- field_boundaries_vect[i]
 
-    # Extract CI values using helper function
-    ci_values <- extract_ci_values(ci_raster, field_vect)
-    valid_values <- ci_values[!is.na(ci_values) & is.finite(ci_values)]
+#     # Extract CI values using helper function
+#     ci_values <- extract_ci_values(ci_raster, field_vect)
+#     valid_values <- ci_values[!is.na(ci_values) & is.finite(ci_values)]
 
-    if (length(valid_values) > 1) {
-      # Calculate % of valid (non-NA) values = gap filling success
-      total_pixels <- length(ci_values)
-      valid_pixels <- length(valid_values)
-      gap_filling_success <- (valid_pixels / total_pixels) * 100
-      na_percent <- ((total_pixels - valid_pixels) / total_pixels) * 100
+#     if (length(valid_values) > 1) {
+#       # Calculate % of valid (non-NA) values = gap filling success
+#       total_pixels <- length(ci_values)
+#       valid_pixels <- length(valid_values)
+#       gap_filling_success <- (valid_pixels / total_pixels) * 100
+#       na_percent <- ((total_pixels - valid_pixels) / total_pixels) * 100
 
-      results_list[[length(results_list) + 1]] <- list(
-        field_idx = i,
-        gap_filling_success = round(gap_filling_success, 2),
-        na_percent_pre_interpolation = round(na_percent, 2),
-        mean_ci = round(mean(valid_values), 2)
-      )
-    } else {
-      # Not enough valid data
-      results_list[[length(results_list) + 1]] <- list(
-        field_idx = i,
-        gap_filling_success = NA_real_,
-        na_percent_pre_interpolation = NA_real_,
-        mean_ci = NA_real_
-      )
-    }
-  }
+#       results_list[[length(results_list) + 1]] <- list(
+#         field_idx = i,
+#         gap_filling_success = round(gap_filling_success, 2),
+#         na_percent_pre_interpolation = round(na_percent, 2),
+#         mean_ci = round(mean(valid_values), 2)
+#       )
+#     } else {
+#       # Not enough valid data
+#       results_list[[length(results_list) + 1]] <- list(
+#         field_idx = i,
+#         gap_filling_success = NA_real_,
+#         na_percent_pre_interpolation = NA_real_,
+#         mean_ci = NA_real_
+#       )
+#     }
+#   }
 
   # Convert accumulated list to data frame in a single operation
   field_results <- do.call(rbind, lapply(results_list, as.data.frame))
diff --git a/r_app/80_utils_common.R b/r_app/80_utils_common.R
index 9c9ff91..04b1793 100644
--- a/r_app/80_utils_common.R
+++ b/r_app/80_utils_common.R
@@ -354,6 +354,78 @@ calculate_cv_trend_long_term <- function(cv_values) {
   })
 }
 
+#' Calculate Gap Filling Score KPI (2σ method)
+#' @param ci_raster Current week CI raster
+#' @param field_boundaries Field boundaries
+#' @return Data frame with field-level gap filling scores
+calculate_gap_filling_kpi <- function(ci_raster, field_boundaries) {
+  safe_log("Calculating Gap Filling Score KPI (placeholder)")
+
+  # Handle both sf and SpatVector inputs
+  if (!inherits(field_boundaries, "SpatVector")) {
+    field_boundaries_vect <- terra::vect(field_boundaries)
+  } else {
+    field_boundaries_vect <- field_boundaries
+  }
+
+  # Ensure field_boundaries_vect is valid and matches field_boundaries dimensions
+  n_fields_vect <- length(field_boundaries_vect)
+  n_fields_sf <- nrow(field_boundaries)
+  
+  if (n_fields_sf != n_fields_vect) {
+    warning(paste("Field boundary mismatch: nrow(field_boundaries)=", n_fields_sf, "vs length(field_boundaries_vect)=", n_fields_vect, ". Using actual SpatVector length."))
+  }
+  
+  field_results <- data.frame()
+
+  for (i in seq_len(nrow(field_boundaries))) {
+    field_name <- field_boundaries$field[i]
+    sub_field_name <- field_boundaries$sub_field[i]
+    field_vect <- field_boundaries_vect[i]
+
+    # Extract CI values using helper function
+    ci_values <- extract_ci_values(ci_raster, field_vect)
+    valid_values <- ci_values[!is.na(ci_values) & is.finite(ci_values)]
+
+    if (length(valid_values) > 1) {
+      # Gap score using 2σ below median to detect outliers
+      median_ci <- median(valid_values)
+      sd_ci <- sd(valid_values)
+      outlier_threshold <- median_ci - (2 * sd_ci)
+      low_ci_pixels <- sum(valid_values < outlier_threshold)
+      total_pixels <- length(valid_values)
+      gap_score <- (low_ci_pixels / total_pixels) * 100
+
+      # Classify gap severity
+      gap_level <- dplyr::case_when(
+        gap_score < 10 ~ "Minimal",
+        gap_score < 25 ~ "Moderate", 
+        TRUE ~ "Significant"
+      )
+
+      field_results <- rbind(field_results, data.frame(
+        field = field_name,
+        sub_field = sub_field_name,
+        gap_level = gap_level,
+        gap_score = gap_score,
+        mean_ci = mean(valid_values),
+        outlier_threshold = outlier_threshold
+      ))
+    } else {
+      # Not enough valid data, fill with NA row
+      field_results <- rbind(field_results, data.frame(
+        field = field_name,
+        sub_field = sub_field_name,
+        gap_level = NA_character_,
+        gap_score = NA_real_,
+        mean_ci = NA_real_,
+        outlier_threshold = NA_real_
+      ))
+    }
+  }
+
+
+
 # ============================================================================
 # HELPER FUNCTIONS
 # ============================================================================