Refactor area calculation functions and enhance error handling; add area conversion utility

This commit is contained in:
Timon 2026-02-24 13:07:55 +01:00
parent b487cc983f
commit e0bfbccf0e
6 changed files with 89 additions and 67 deletions

View file

@ -356,8 +356,8 @@ calculate_tch_forecasted_kpi <- function(field_statistics, harvesting_data = NUL
#' Uses FOUR_WEEK_TREND_* thresholds defined in 80_utils_common.R: #' Uses FOUR_WEEK_TREND_* thresholds defined in 80_utils_common.R:
#' - FOUR_WEEK_TREND_STRONG_GROWTH_MIN (0.3) #' - FOUR_WEEK_TREND_STRONG_GROWTH_MIN (0.3)
#' - FOUR_WEEK_TREND_GROWTH_MIN (0.1) #' - FOUR_WEEK_TREND_GROWTH_MIN (0.1)
#' - FOUR_WEEK_TREND_STABLE_THRESHOLD (0.1 and -0.1) #' - FOUR_WEEK_TREND_STABLE_THRESHOLD (0.1; used as lower bound via -FOUR_WEEK_TREND_STABLE_THRESHOLD)
#' - FOUR_WEEK_TREND_WEAK_DECLINE_THRESHOLD (-0.3) #' - FOUR_WEEK_TREND_WEAK_DECLINE_THRESHOLD (-0.1)
#' - FOUR_WEEK_TREND_STRONG_DECLINE_MAX (-0.3) #' - FOUR_WEEK_TREND_STRONG_DECLINE_MAX (-0.3)
#' #'
calculate_growth_decline_kpi <- function(ci_values_list) { calculate_growth_decline_kpi <- function(ci_values_list) {
@ -632,14 +632,25 @@ create_field_detail_table <- function(field_boundaries_sf, all_kpis, current_wee
# ============================================ # ============================================
# GROUP 0b: FIELD AREA (from geometry) # GROUP 0b: FIELD AREA (from geometry)
# ============================================ # ============================================
# Precompute unit label before tryCatch to avoid re-raising errors in error handler
unit_label <- get_area_unit_label(AREA_UNIT_PREFERENCE)
area_col_name <- paste0("Area_", unit_label)
# Calculate field area using unified function and preferred unit # Calculate field area using unified function and preferred unit
tryCatch({ tryCatch({
field_areas <- calculate_area_from_geometry(field_boundaries_sf, unit = AREA_UNIT_PREFERENCE) field_areas <- calculate_area_from_geometry(field_boundaries_sf, unit = AREA_UNIT_PREFERENCE)
unit_label <- get_area_unit_label(AREA_UNIT_PREFERENCE)
area_col_name <- paste0("Area_", unit_label) # Validate that calculated areas match the number of fields
n_fields <- nrow(field_boundaries_sf)
if (length(field_areas) != n_fields) {
stop(
"Mismatch: calculate_area_from_geometry returned ", length(field_areas),
" areas but field_boundaries_sf has ", n_fields, " fields"
)
}
area_df <- data.frame( area_df <- data.frame(
field_idx = seq_along(field_areas), field_idx = seq_len(n_fields),
area_value = field_areas, area_value = field_areas,
stringsAsFactors = FALSE stringsAsFactors = FALSE
) )
@ -649,7 +660,7 @@ create_field_detail_table <- function(field_boundaries_sf, all_kpis, current_wee
left_join(area_df, by = "field_idx") left_join(area_df, by = "field_idx")
}, error = function(e) { }, error = function(e) {
message(paste("Warning: Could not calculate field areas:", e$message)) message(paste("Warning: Could not calculate field areas:", e$message))
result[[paste0("Area_", get_area_unit_label(AREA_UNIT_PREFERENCE))]] <<- NA_real_ result[[area_col_name]] <<- NA_real_
}) })
# ============================================ # ============================================

View file

@ -53,12 +53,22 @@ CI_CHANGE_INCREASE_THRESHOLD <- CI_CHANGE_RAPID_GROWTH_THRESHOLD # Weekly CI
calculate_field_acreages <- function(field_boundaries_sf, unit = AREA_UNIT_PREFERENCE) { calculate_field_acreages <- function(field_boundaries_sf, unit = AREA_UNIT_PREFERENCE) {
tryCatch({ tryCatch({
# Get field identifier (handles pivot.geojson structure) # Get field identifier (handles pivot.geojson structure)
field_names <- field_boundaries_sf %>% # Use intersect() to find the first matching column, avoiding pull() ambiguity
sf::st_drop_geometry() %>% field_names <- {
pull(any_of(c("field", "field_id", "Field_id", "name", "Name"))) available_cols <- intersect(
colnames(field_boundaries_sf),
c("field", "field_id", "Field_id", "name", "Name")
)
if (length(field_names) == 0 || all(is.na(field_names))) { if (length(available_cols) > 0) {
field_names <- seq_len(nrow(field_boundaries_sf)) # Extract the first matching column
field_boundaries_sf %>%
sf::st_drop_geometry() %>%
pull(available_cols[1])
} else {
# Fall back to sequential indices if no field name column exists
seq_len(nrow(field_boundaries_sf))
}
} }
# Use unified area calculation function # Use unified area calculation function
@ -70,15 +80,15 @@ calculate_field_acreages <- function(field_boundaries_sf, unit = AREA_UNIT_PREFE
result_df <- data.frame( result_df <- data.frame(
field = field_names, field = field_names,
area = areas, warning(paste("Could not calculate areas from geometries -", e$message)) stringsAsFactors = FALSE
stringsAsFactors = FALSE
) )
colnames(result_df) <- c("field", col_name) colnames(result_df) <- c("field", col_name)
# Aggregate by field to handle multi-row fields (e.g., sub_fields) # Aggregate by field to handle multi-row fields (e.g., sub_fields)
# Use bare lambda to preserve original column name (not list() which creates suffixes)
result_df %>% result_df %>%
group_by(field) %>% group_by(field) %>%
summarise(across(all_of(col_name), list(~ sum(., na.rm = TRUE))), .groups = "drop") summarise(across(all_of(col_name), ~ sum(., na.rm = TRUE)), .groups = "drop")
}, error = function(e) { }, error = function(e) {
message(paste("Warning: Could not calculate areas from geometries -", e$message)) message(paste("Warning: Could not calculate areas from geometries -", e$message))
unit_label <- get_area_unit_label(unit) unit_label <- get_area_unit_label(unit)
@ -600,21 +610,6 @@ calculate_field_analysis_cane_supply <- function(setup,
reports_dir reports_dir
) )
# cat("\n--- Per-field Results (first 10) ---\n")
# available_cols <- c("Field_id", "Acreage", "Age_week", "Mean_CI", "Four_week_trend", "Status_Alert", "Cloud_category")
# available_cols <- available_cols[available_cols %in% names(field_analysis_df)]
# if (length(available_cols) > 0) {
# print(head(field_analysis_df[, available_cols], 10))
# }
# ========== PHASE 10: CALCULATE FARM-LEVEL KPIS ==========
# farm_kpi_results <- calculate_farm_level_kpis(
# field_analysis_df,
# current_week,
# current_year,
# end_date
# )
# For now, farm-level KPIs are not implemented in CANE_SUPPLY workflow # For now, farm-level KPIs are not implemented in CANE_SUPPLY workflow
farm_kpi_results <- NULL farm_kpi_results <- NULL

View file

@ -86,12 +86,12 @@ PHASE_DEFINITIONS <- data.frame(
#' Calculate field area from geometry in specified unit #' Calculate field area from geometry in specified unit
#' #'
#' Unified function for calculating polygon area from sf or SpatVect geometries. #' Unified function for calculating polygon area from sf or SpatVector geometries.
#' Uses equal-area projection (EPSG:6933) for accurate calculations across all zones. #' Uses equal-area projection (EPSG:6933) for accurate calculations across all zones.
#' #'
#' @param geometry sf or SpatVect object containing field polygons #' @param geometry sf or SpatVector object containing field polygons
#' If sf: must have geometry column (auto-detected) #' If sf: must have geometry column (auto-detected)
#' If SpatVect: terra object with geometry #' If SpatVector: terra object with geometry
#' @param unit Character. Output unit: "hectare" (default) or "acre" #' @param unit Character. Output unit: "hectare" (default) or "acre"
#' #'
#' @return Numeric vector of areas in specified unit #' @return Numeric vector of areas in specified unit
@ -118,7 +118,7 @@ PHASE_DEFINITIONS <- data.frame(
#' areas_ha <- calculate_area_from_geometry(fields_sf, unit = "hectare") #' areas_ha <- calculate_area_from_geometry(fields_sf, unit = "hectare")
#' areas_ac <- calculate_area_from_geometry(fields_sf, unit = "acre") #' areas_ac <- calculate_area_from_geometry(fields_sf, unit = "acre")
#' #'
#' # With SpatVect #' # With SpatVector
#' library(terra) #' library(terra)
#' fields_vect <- vect("pivot.geojson") #' fields_vect <- vect("pivot.geojson")
#' areas_ha <- calculate_area_from_geometry(fields_vect, unit = "hectare") #' areas_ha <- calculate_area_from_geometry(fields_vect, unit = "hectare")
@ -137,12 +137,12 @@ calculate_area_from_geometry <- function(geometry, unit = "hectare") {
# Handle sf object # Handle sf object
geometry_proj <- sf::st_transform(geometry, 6933) geometry_proj <- sf::st_transform(geometry, 6933)
areas_m2 <- as.numeric(sf::st_area(geometry_proj)) areas_m2 <- as.numeric(sf::st_area(geometry_proj))
} else if (inherits(geometry, "SpatVect")) { } else if (inherits(geometry, "SpatVector")) {
# Handle terra SpatVect object # Handle terra SpatVector object
geometry_proj <- terra::project(geometry, "EPSG:6933") geometry_proj <- terra::project(geometry, "EPSG:6933")
areas_m2 <- as.numeric(terra::expanse(geometry_proj)) areas_m2 <- as.numeric(terra::expanse(geometry_proj))
} else { } else {
stop("geometry must be an sf or terra SpatVect object. Got: ", stop("geometry must be an sf or terra SpatVector object. Got: ",
paste(class(geometry), collapse = ", ")) paste(class(geometry), collapse = ", "))
} }
@ -1068,12 +1068,14 @@ calculate_field_statistics <- function(field_boundaries_sf, week_num, year,
} }
# Guard: detect cloud-masked dates (CI == 0 indicates no-data) # Guard: detect cloud-masked dates (CI == 0 indicates no-data)
# When any extracted value is 0, treat the entire date as cloud-masked # Calculate proportion of zero values; only treat as cloud-masked if proportion exceeds threshold
has_zeros <- any(extracted$CI == 0, na.rm = TRUE) prop_zero <- mean(extracted$CI == 0, na.rm = TRUE)
cloud_threshold <- 0.25 # 25% tolerance for zero pixels
if (has_zeros) { if (prop_zero > cloud_threshold) {
# Cloud-masked date: skip temporal analysis, set stats to NA # Cloud-masked date: skip temporal analysis, set stats to NA
message(paste(" [CLOUD] Field", field_name, "- entire date is cloud-masked (CI==0)")) message(paste(" [CLOUD] Field", field_name, "- date is cloud-masked (",
round(prop_zero * 100, 1), "% CI==0)"))
results_list[[length(results_list) + 1]] <- data.frame( results_list[[length(results_list) + 1]] <- data.frame(
Field_id = field_name, Field_id = field_name,
@ -1089,7 +1091,7 @@ calculate_field_statistics <- function(field_boundaries_sf, week_num, year,
next next
} }
ci_vals <- extracted$CI[!is.na(extracted$CI)] ci_vals <- extracted$CI[!is.na(extracted$CI) & extracted$CI > 0]
if (length(ci_vals) == 0) { if (length(ci_vals) == 0) {
next next
@ -1847,21 +1849,8 @@ calculate_yield_prediction_kpi <- function(field_boundaries, harvesting_data, cu
pred_rf_current_season$predicted_Tcha < yield_quartiles[3], na.rm = TRUE) pred_rf_current_season$predicted_Tcha < yield_quartiles[3], na.rm = TRUE)
lowest_25_count <- sum(pred_rf_current_season$predicted_Tcha < yield_quartiles[1], na.rm = TRUE) lowest_25_count <- sum(pred_rf_current_season$predicted_Tcha < yield_quartiles[1], na.rm = TRUE)
# Calculate total area # Calculate total area using centralized function
if (!inherits(field_boundaries, "SpatVector")) { field_areas <- calculate_area_from_geometry(field_boundaries, unit = "hectare")
field_boundaries_vect <- terra::vect(field_boundaries)
} else {
field_boundaries_vect <- field_boundaries
}
# Handle both sf and SpatVector inputs for area calculation
if (inherits(field_boundaries, "sf")) {
field_boundaries_projected <- sf::st_transform(field_boundaries, "EPSG:6933")
field_areas <- sf::st_area(field_boundaries_projected) / 10000 # m² to hectares
} else {
field_boundaries_projected <- terra::project(field_boundaries_vect, "EPSG:6933")
field_areas <- terra::expanse(field_boundaries_projected) / 10000
}
total_area <- sum(as.numeric(field_areas)) total_area <- sum(as.numeric(field_areas))
safe_log(paste("Total area:", round(total_area, 1), "hectares")) safe_log(paste("Total area:", round(total_area, 1), "hectares"))

View file

@ -531,19 +531,20 @@ map_trend_to_arrow <- function(text_vec, include_text = FALSE) {
} }
# Determine category and build output with translated labels # Determine category and build output with translated labels
if (grepl("strong growth", text)) { # Using word-boundary anchored patterns (perl=TRUE) to avoid substring mis-matches
if (grepl("\\bstrong growth\\b", text, perl = TRUE)) {
arrow <- "↑↑" arrow <- "↑↑"
trans_key <- "Strong growth" trans_key <- "Strong growth"
} else if (grepl("slight growth|weak growth|growth|increasing", text)) { } else if (grepl("\\b(?:slight|weak) growth\\b|(?<!no\\s)\\bgrowth\\b|\\bincreasing\\b", text, perl = TRUE)) {
arrow <- "↑" arrow <- "↑"
trans_key <- "Slight growth" trans_key <- "Slight growth"
} else if (grepl("stable|no growth", text)) { } else if (grepl("\\bstable\\b|\\bno growth\\b", text, perl = TRUE)) {
arrow <- "→" arrow <- "→"
trans_key <- "Stable" trans_key <- "Stable"
} else if (grepl("weak decline|slight decline|moderate decline", text)) { } else if (grepl("\\b(?:weak|slight|moderate) decline\\b", text, perl = TRUE)) {
arrow <- "↓" arrow <- "↓"
trans_key <- "Slight decline" trans_key <- "Slight decline"
} else if (grepl("strong decline|severe", text)) { } else if (grepl("\\bstrong decline\\b|\\bsevere\\b", text, perl = TRUE)) {
arrow <- "↓↓" arrow <- "↓↓"
trans_key <- "Strong decline" trans_key <- "Strong decline"
} else { } else {
@ -748,8 +749,7 @@ if (exists("summary_tables") && !is.null(summary_tables) && length(summary_table
# Translate the table for visualization # Translate the table for visualization
names(display_df) <- c(tr_key("KPI"), tr_key("Level"), tr_key("Count"), tr_key("Percent")) names(display_df) <- c(tr_key("KPI"), tr_key("Level"), tr_key("Count"), tr_key("Percent"))
display_df[, 1:2] <- lapply(display_df[, 1:2], function(col) sapply(col, t)) display_df[, 1:2] <- lapply(display_df[, 1:2], function(col) sapply(col, tr_key))
ft <- flextable(display_df) %>% ft <- flextable(display_df) %>%
merge_v(j = tr_key("KPI")) %>% merge_v(j = tr_key("KPI")) %>%
autofit() autofit()
@ -1864,7 +1864,7 @@ metadata_info <- data.frame(
format(Sys.time(), "%Y-%m-%d %H:%M:%S"), format(Sys.time(), "%Y-%m-%d %H:%M:%S"),
paste(tr_key("project"), toupper(project_dir)), paste(tr_key("project"), toupper(project_dir)),
paste(tr_key("week"), current_week, tr_key("of"), year), paste(tr_key("week"), current_week, tr_key("of"), year),
ifelse(exists("AllPivots0"), nrow(AllPivots0 %>% filter(!is.na(field)) %>% group_by(field) %>% summarise()), tr_key("unknown")), if (exists("AllPivots0")) { nrow(AllPivots0 %>% filter(!is.na(field)) %>% group_by(field) %>% summarise()) } else { tr_key("unknown") },
tr_key("next_wed") tr_key("next_wed")
) )
) )

View file

@ -64,8 +64,35 @@ get_area_unit_label <- function(unit = AREA_UNIT_PREFERENCE) {
} else if (unit_lower == "acre") { } else if (unit_lower == "acre") {
return("ac") return("ac")
} else { } else {
warning("Unknown unit: ", unit, ". Defaulting to 'ha'") warning("Unknown unit '", unit, "'. Falling back to AREA_UNIT_PREFERENCE ('", AREA_UNIT_PREFERENCE, "')")
return("ha") return(get_area_unit_label(AREA_UNIT_PREFERENCE))
}
}
#' Get area conversion factor from hectares to requested unit
#'
#' Returns the numeric multiplier to convert 1 hectare into the requested unit.
#' Case-insensitive unit matching. On unknown unit, warns and defaults to 1 (hectare).
#'
#' @param unit Character. Unit preference ("hectare" or "acre")
#' @return Numeric. Conversion factor (1 for hectare, 2.47105 for acre)
#'
#' @details
#' Conversion factors:
#' - "hectare" → 1 (identity)
#' - "acre" → 2.47105 (1 hectare ≈ 2.47105 acres)
#'
#'
#' @export
get_area_conversion_factor <- function(unit = AREA_UNIT_PREFERENCE) {
unit_lower <- tolower(unit)
if (unit_lower == "hectare") {
return(1)
} else if (unit_lower == "acre") {
return(2.47105) # 1 hectare = 2.47105 acres (inverse of 0.404686)
} else {
warning("Unknown unit '", unit, "'. Falling back to AREA_UNIT_PREFERENCE ('", AREA_UNIT_PREFERENCE, "')")
return(get_area_conversion_factor(AREA_UNIT_PREFERENCE))
} }
} }

Binary file not shown.