Timon/SmartCane
1
1
Fork 0
Code Issues Pull requests Actions Packages Projects Releases Wiki Activity

SC-112 Phase 2 Complete: Restructure Script 80 utilities by client type

Browse source 
- Consolidate 80_kpi_utils.R, 80_weekly_stats_utils.R, 80_report_building_utils.R into three client-aware files:
  - 80_utils_common.R (50+ functions): Shared utilities, constants, and helpers
  - 80_utils_agronomic_support.R: AURA-specific KPI calculations (6 KPIs)
  - 80_utils_cane_supply.R: ANGATA placeholder for future expansion
- Move all internal helpers (calculate_cv, calculate_change_percentages, calculate_spatial_autocorrelation, extract_ci_values, etc.) to common
- Add MORAN_THRESHOLD_* constants to common
- Fix parameters_project.R field boundaries path (removed extra 'Data' directory)
- Update 80_calculate_kpis.R with conditional client-type sourcing logic
- Validate both ANGATA (cane_supply) and AURA (agronomic_support) workflows with comprehensive testing
- All 96+ functions accounted for; ready for production use
This commit is contained in:
Timon 2026-02-03 17:21:35 +01:00
parent 20e3d5b4a3
commit 8d560ff155
10 changed files with 1651 additions and 3786 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

252
r_app/40_mosaic_creation.R
View file

@ -1,252 +0,0 @@
# filepath: c:\Users\timon\Resilience BV\4020 SCane ESA DEMO - Documenten\General\4020 SCDEMO Team\4020 TechnicalData\WP3\smartcane\r_app\mosaic_creation.R
#
# MOSAIC_CREATION.R
# ===============
# This script creates weekly mosaics from daily satellite imagery.
# It handles command-line arguments and initiates the mosaic creation process.
#
# Usage: Rscript mosaic_creation.R [end_date] [offset] [project_dir] [file_name] [data_source]
# - end_date: End date for processing (YYYY-MM-DD format)
# - offset: Number of days to look back from end_date (typically 7 for one week)
# - project_dir: Project directory name (e.g., "aura", "angata", "chemba", "esa")
# - file_name: Optional custom output file name (leave empty "" to use default: week_WW_YYYY.tif)
# - data_source: Optional data source folder (e.g., "merged_tif" or "merged_tif_8b")
# If not provided, auto-detects which folder contains actual data
#
# Examples:
# & "C:\Program Files\R\R-4.4.3\bin\x64\Rscript.exe" r_app/40_mosaic_creation.R 2026-01-12 7 aura
# & "C:\Program Files\R\R-4.4.3\bin\x64\Rscript.exe" r_app/40_mosaic_creation.R 2025-12-24 7 aura "" "merged_tif"
#
# 1. Load required packages
# -----------------------
suppressPackageStartupMessages({
library(sf)
library(terra)
library(tidyverse)
library(lubridate)
library(here)
})
# 2. Process command line arguments and run mosaic creation
# ------------------------------------------------------
main <- function() {
# Capture command line arguments
args <- commandArgs(trailingOnly = TRUE)
# Process project_dir argument with default
if (length(args) >= 3 && !is.na(args[3])) {
project_dir <- as.character(args[3])
} else if (exists("project_dir", envir = .GlobalEnv)) {
project_dir <- get("project_dir", envir = .GlobalEnv)
} else {
# Default project directory
project_dir <- "angata"
message("No project_dir provided. Using default:", project_dir)
}
# Make project_dir available globally so parameters_project.R can use it
assign("project_dir", project_dir, envir = .GlobalEnv)
# Process end_date argument with default
if (length(args) >= 1 && !is.na(args[1])) {
# Parse date explicitly in YYYY-MM-DD format from command line
end_date <- as.Date(args[1], format = "%Y-%m-%d")
if (is.na(end_date)) {
message("Invalid end_date provided. Using current date.")
end_date <- Sys.Date()
}
} else if (exists("end_date_str", envir = .GlobalEnv)) {
end_date <- as.Date(get("end_date_str", envir = .GlobalEnv))
} else {
# Default to current date if no argument is provided
end_date <- Sys.Date()
message("No end_date provided. Using current date: ", format(end_date))
}
# Process offset argument with default
if (length(args) >= 2 && !is.na(args[2])) {
offset <- as.numeric(args[2])
if (is.na(offset) || offset <= 0) {
message("Invalid offset provided. Using default (7 days).")
offset <- 7
}
} else {
# Default to 7 days if no argument is provided
offset <- 7
message("No offset provided. Using default:", offset, "days")
}
# Process data_source argument (optional, passed from pipeline)
# If provided, use it; otherwise auto-detect
data_source_from_args <- NULL
if (length(args) >= 5 && !is.na(args[5]) && nchar(args[5]) > 0) {
data_source_from_args <- as.character(args[5])
message("Data source explicitly provided via arguments: ", data_source_from_args)
}
# 3. Initialize project configuration
# --------------------------------
# Detect which data source directory exists (merged_tif or merged_tif_8b)
# IMPORTANT: Only consider a folder as valid if it contains actual files
laravel_storage <- here::here("laravel_app/storage/app", project_dir)
# Load centralized path structure
tryCatch({
source("r_app/parameters_project.R")
paths <- setup_project_directories(project_dir)
}, error = function(e) {
message("Note: Could not open files from r_app directory")
message("Attempting to source from default directory instead...")
tryCatch({
source("parameters_project.R")
paths <- setup_project_directories(project_dir)
message("✓ Successfully sourced files from default directory")
}, error = function(e) {
stop("Failed to source required files from both 'r_app' and default directories.")
})
})
data_source <- if (has_8b_data) {
message("Auto-detected data source: merged_tif_8b (8-band optimized) - contains files")
"merged_tif_8b"
} else if (has_legacy_data) {
message("Auto-detected data source: merged_tif (legacy 4-band) - contains files")
"merged_tif"
} else {
message("Warning: No valid data source found (both folders empty or missing). Using default: merged_tif")
"merged_tif"
}
}
# Set global data_source for parameters_project.R
assign("data_source", data_source, envir = .GlobalEnv)
tryCatch({
source("r_app/parameters_project.R")
source("r_app/00_common_utils.R")
source("r_app/40_mosaic_creation_utils.R")
safe_log(paste("Successfully sourced files from 'r_app' directory."))
}, error = function(e) {
message("Note: Could not open files from r_app directory")
message("Attempting to source from default directory instead...")
tryCatch({
source("parameters_project.R")
paths <- setup_project_directories(project_dir)
message("✓ Successfully sourced files from default directory")
}, error = function(e) {
stop("Failed to source required files from both 'r_app' and default directories.")
})
})
# Use centralized paths (no need to manually construct or create dirs)
merged_final <- paths$growth_model_interpolated_dir # or merged_final_tif if needed
daily_vrt <- paths$vrt_dir
safe_log(paste("Using growth model/mosaic directory:", merged_final))
safe_log(paste("Using daily VRT directory:", daily_vrt))
# 4. Generate date range for processing
# ---------------------------------
dates <- date_list(end_date, offset)
safe_log(paste("Processing data for week", dates$week, "of", dates$year))
# Create output filename
# Only use custom filename if explicitly provided (not empty string)
file_name_tif <- if (length(args) >= 4 && !is.na(args[4]) && nchar(args[4]) > 0) {
as.character(args[4])
} else {
paste0("week_", sprintf("%02d", dates$week), "_", dates$year, ".tif")
}
safe_log(paste("Output will be saved as:", file_name_tif))
# 5. Create weekly mosaics - route based on project tile detection
# ---------------------------------------------------------------
# The use_tile_mosaic flag is auto-detected by parameters_project.R
# based on whether tiles exist in merged_final_tif/
if (!exists("use_tile_mosaic")) {
# Fallback detection if flag not set (shouldn't happen)
merged_final_dir <- file.path(laravel_storage, "merged_final_tif")
tile_detection <- detect_tile_structure_from_merged_final(merged_final_dir)
use_tile_mosaic <- tile_detection$has_tiles
}
if (use_tile_mosaic) {
# TILE-BASED APPROACH: Create per-tile weekly MAX mosaics
# This is used for projects like Angata with large ROIs requiring spatial partitioning
# Input data comes from merged_final_tif/{grid_size}/{DATE}/{DATE}_XX.tif (5-band tiles from script 20)
tryCatch({
safe_log("Starting per-tile mosaic creation (tile-based approach)...")
# Detect grid size from merged_final_tif folder structure
# Expected: merged_final_tif/5x5/ or merged_final_tif/10x10/ etc.
merged_final_base <- file.path(laravel_storage, "merged_final_tif")
grid_subfolders <- list.dirs(merged_final_base, full.names = FALSE, recursive = FALSE)
# Look for grid size patterns like "5x5", "10x10", "20x20"
grid_patterns <- grep("^\\d+x\\d+$", grid_subfolders, value = TRUE)
if (length(grid_patterns) == 0) {
stop("No grid size subfolder found in merged_final_tif/ (expected: 5x5, 10x10, etc.)")
}
grid_size <- grid_patterns[1] # Use first grid size found
safe_log(paste("Detected grid size:", grid_size))
# Point to the grid-specific merged_final_tif directory
merged_final_with_grid <- file.path(merged_final_base, grid_size)
# Set output directory for per-tile mosaics, organized by grid size (from centralized paths)
# Output: weekly_tile_max/{grid_size}/week_WW_YYYY_TT.tif
tile_output_base <- file.path(paths$weekly_tile_max_dir, grid_size)
# Note: no dir.create needed - paths$weekly_tile_max_dir already created by setup_project_directories()
dir.create(tile_output_base, recursive = TRUE, showWarnings = FALSE) # Create grid-size subfolder
created_tile_files <- create_weekly_mosaic_from_tiles(
dates = dates,
merged_final_dir = merged_final_with_grid,
tile_output_dir = tile_output_base,
field_boundaries = field_boundaries
)
safe_log(paste("✓ Per-tile mosaic creation completed - created",
length(created_tile_files), "tile files"))
}, error = function(e) {
safe_log(paste("ERROR in tile-based mosaic creation:", e$message), "ERROR")
traceback()
stop("Mosaic creation failed")
})
} else {
# SINGLE-FILE APPROACH: Create single weekly mosaic file
# This is used for legacy projects (ESA, Chemba, Aura) expecting single-file output
tryCatch({
safe_log("Starting single-file mosaic creation (backward-compatible approach)...")
# Set output directory for single-file mosaics (from centralized paths)
single_file_output_dir <- paths$weekly_mosaic_dir
created_file <- create_weekly_mosaic(
dates = dates,
field_boundaries = field_boundaries,
daily_vrt_dir = daily_vrt,
merged_final_dir = merged_final,
output_dir = single_file_output_dir,
file_name_tif = file_name_tif,
create_plots = FALSE
)
safe_log(paste("✓ Single-file mosaic creation completed:", created_file))
}, error = function(e) {
safe_log(paste("ERROR in single-file mosaic creation:", e$message), "ERROR")
traceback()
stop("Mosaic creation failed")
})
}
}
if (sys.nframe() == 0) {
main()
}

286
r_app/40_mosaic_creation_tile_utils.R
View file

@ -1,286 +0,0 @@
# MOSAIC_CREATION_TILE_UTILS.R
# ============================
# Tile-based processing utilities for scalable weekly mosaic creation
#
# STRATEGY: Memory-efficient, scalable approach for large study areas
# Split daily full-extent mosaics into fixed-size tiles (e.g., 5×5 km),
# then process each tile position across all days, and reassemble.
#
# WORKFLOW:
# Daily full-extent TIFF
# ↓
# Split into N×M tiles (based on area size / tile_size)
# E.g., 50×80 km area with 5 km tiles = 10×16 = 160 tiles
# ↓
# For EACH TILE POSITION (1 to 160):
# - Load that tile from all 7 days
# - Create MAX composite for that tile
# - Save tile_NNN_MAX.tif
# ↓
# Reassemble all 160 MAX tiles back into full-extent
# ↓
# Save weekly full-extent mosaic
#
# KEY BENEFIT: Memory usage ~50 MB per tile (5×5 km), not 1.3 GB (full 50×80 km)
# Scales automatically: bigger area = more tiles, same memory per tile
#
# TILE_SIZE: Configurable (default 5000 m = 5 km, adjustable via parameter)
#' Simple tile-based processing using terra::makeTiles()
#'
#' Calculates tile grid based on desired tile SIZE (e.g., 5 km), not grid count.
#' This makes it SCALABLE: bigger area = more tiles automatically.
#'
#' @param raster_path Path to raster to split
#' @param output_dir Directory to save tiles
#' @param tile_size_m Tile size in meters (default: 5000 m = 5 km)
#' @return List with tiles (file paths) and metadata
#'
split_raster_into_tiles <- function(raster_path, output_dir, tile_size_m = 5000) {
# Load raster
r <- terra::rast(raster_path)
# Calculate how many tiles we need based on raster extent and tile size
x_range <- terra::ext(r)$xmax - terra::ext(r)$xmin
y_range <- terra::ext(r)$ymax - terra::ext(r)$ymin
n_tiles_x <- ceiling(x_range / tile_size_m)
n_tiles_y <- ceiling(y_range / tile_size_m)
safe_log(paste("Splitting into", n_tiles_x, "×", n_tiles_y, "tiles of",
tile_size_m / 1000, "km"))
# Create output directory
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
# Use terra::makeTiles() - it splits based on n = c(rows, cols)
tiles <- terra::makeTiles(
r,
n = c(n_tiles_y, n_tiles_x), # rows, cols
filename = file.path(output_dir, "tile_.tif"),
overwrite = TRUE
)
safe_log(paste("Created", length(tiles), "tile files"))
return(list(
tiles = tiles,
n_tiles = length(tiles),
n_tiles_x = n_tiles_x,
n_tiles_y = n_tiles_y,
extent = terra::ext(r),
raster = raster_path
))
}
#' Create weekly MAX mosaic using TRUE tile-based processing
#'
#' TILE-BASED WORKFLOW (Memory efficient):
#' 1. Calculate tile grid dynamically from extent and FIXED tile_size (e.g., 5 km)
#' 2. For EACH TILE across all 7 days:
#' - Load that tile from each daily file (small ~50 MB, not 1.3 GB)
#' - Create MAX composite for just that tile
#' 3. Reassemble all tiles into final full-extent mosaic
#'
#' SCALABILITY: Fixed 5 km tile size means bigger area = more tiles (automatic scaling)
#'
#' @param dates List from date_list() with days_filter
#' @param merged_final_dir Directory with daily merged full-extent TIFFs
#' @param final_output_dir Directory to store final reassembled mosaic
#' @param file_name_tif Output filename for final mosaic
#' @param tile_size_m Tile size in meters (default: 5000 m = 5 km). Bigger area automatically gets more tiles.
#' @return Path to final reassembled weekly mosaic
#'
create_weekly_mosaic_tiled <- function(dates, merged_final_dir,
final_output_dir, file_name_tif,
tile_size_m = 5000) {
# Get daily files for this week
daily_files <- list.files(merged_final_dir, full.names = TRUE, pattern = "\\.tif$")
daily_files <- daily_files[grepl(paste(dates$days_filter, collapse = "|"), daily_files)]
if (length(daily_files) == 0) {
safe_log("No daily files found for this week", "ERROR")
return(NULL)
}
safe_log(paste("Found", length(daily_files), "daily files for week", dates$week))
# Load first daily file to get extent and calculate tile grid
safe_log("Step 1: Loading first daily file to establish tile structure")
first_raster <- terra::rast(daily_files[1])
# Get CRS and convert tile_size_m to appropriate units
raster_crs <- terra::crs(first_raster)
# If raster is in lat/lon (geographic), convert tile_size_m to degrees
# 1 degree latitude ≈ 111 km, so 5000 m ≈ 0.045 degrees
# Check for GEOGCRS (geographic coordinate system) in WKT string
is_geographic <- grepl("GEOGCRS|longlat|geographic|ANGLEUNIT.*degree",
as.character(raster_crs), ignore.case = TRUE)
if (is_geographic) {
# Geographic CRS - convert meters to degrees
# At equator: 1 degree ≈ 111 km
tile_size_degrees <- tile_size_m / 111000 # 111 km per degree
safe_log(paste("Raster is in geographic CRS (lat/lon). Converting", tile_size_m / 1000,
"km to", round(tile_size_degrees, 4), "degrees"))
} else {
# Projected CRS - use meters directly
tile_size_degrees <- tile_size_m
safe_log(paste("Raster is in projected CRS. Using", tile_size_m / 1000, "km directly"))
}
# Calculate n_tiles dynamically from extent and tile_size
x_range <- terra::ext(first_raster)$xmax - terra::ext(first_raster)$xmin
y_range <- terra::ext(first_raster)$ymax - terra::ext(first_raster)$ymin
n_tiles_x <- ceiling(x_range / tile_size_degrees)
n_tiles_y <- ceiling(y_range / tile_size_degrees)
n_tiles_total <- n_tiles_x * n_tiles_y
safe_log(paste("Step 2: Creating tile grid:", tile_size_m / 1000, "km tiles"))
safe_log(paste(" Extent:", round(x_range, 4), "° ×", round(y_range, 4), "°"))
safe_log(paste(" Grid:", n_tiles_y, "rows ×", n_tiles_x, "columns =", n_tiles_total, "tiles"))
# Calculate tile extents mathematically (no need to save temp files)
extent <- terra::ext(first_raster)
x_min <- extent$xmin
y_min <- extent$ymin
# Create list of tile extents
tile_extents <- list()
tile_idx <- 0
for (row in 1:n_tiles_y) {
for (col in 1:n_tiles_x) {
tile_idx <- tile_idx + 1
# Calculate this tile's bounds
tile_xmin <- x_min + (col - 1) * tile_size_degrees
tile_xmax <- min(tile_xmin + tile_size_degrees, extent$xmax)
tile_ymin <- y_min + (row - 1) * tile_size_degrees
tile_ymax <- min(tile_ymin + tile_size_degrees, extent$ymax)
tile_extents[[tile_idx]] <- terra::ext(tile_xmin, tile_xmax, tile_ymin, tile_ymax)
}
}
safe_log(paste("Calculated", length(tile_extents), "tile extents"))
# Save tiles to Laravel storage directory
tile_storage_dir <- file.path("laravel_app", "storage", "app", "angata", "daily_tiles")
dir.create(tile_storage_dir, recursive = TRUE, showWarnings = FALSE)
# For each tile, load from all daily files and create MAX
safe_log("Step 3: Processing tiles (load per-tile across all days, create MAX for each)")
tile_files_list <- list()
for (tile_idx in seq_along(tile_extents)) {
if (tile_idx %% max(1, ceiling(n_tiles_total / 10)) == 0 || tile_idx == 1) {
safe_log(paste(" Processing tile", tile_idx, "of", n_tiles_total))
}
# Get this tile's extent
tile_extent <- tile_extents[[tile_idx]]
# Load and crop all daily files to this tile extent
daily_tile_data <- list()
for (file_idx in seq_along(daily_files)) {
tryCatch({
r <- terra::rast(daily_files[file_idx])
cropped <- terra::crop(r, tile_extent, snap = "in")
daily_tile_data[[file_idx]] <- cropped
}, error = function(e) {
safe_log(paste("Warning: Could not crop tile", tile_idx, "from",
basename(daily_files[file_idx])), "WARNING")
})
}
if (length(daily_tile_data) == 0) {
safe_log(paste("No valid data for tile", tile_idx), "WARNING")
next
}
# Create MAX composite for this tile
if (length(daily_tile_data) == 1) {
tile_max <- daily_tile_data[[1]]
} else {
collection <- terra::sprc(daily_tile_data)
tile_max <- terra::mosaic(collection, fun = "max")
}
# Save tile to disk (keeps memory low)
tile_file <- file.path(tile_storage_dir, sprintf("tile_%04d.tif", tile_idx))
terra::writeRaster(tile_max, tile_file, overwrite = TRUE)
tile_files_list[[tile_idx]] <- tile_file
}
if (length(tile_files_list) == 0) {
safe_log("No tiles processed successfully", "ERROR")
return(NULL)
}
safe_log(paste("Step 4: Reassembling", length(tile_files_list), "tiles from disk into full-extent mosaic"))
# Load all tile files and reassemble
tile_rasters <- lapply(tile_files_list, terra::rast)
collection <- terra::sprc(tile_rasters)
final_mosaic <- terra::mosaic(collection, fun = "first")
safe_log("Step 5: Saving final reassembled mosaic")
# Save
dir.create(final_output_dir, recursive = TRUE, showWarnings = FALSE)
final_file_path <- file.path(final_output_dir, file_name_tif)
tryCatch({
terra::writeRaster(final_mosaic, final_file_path, overwrite = TRUE)
safe_log(paste("✓ Weekly mosaic saved to:", final_file_path))
}, error = function(e) {
safe_log(paste("Error saving mosaic:", e$message), "ERROR")
return(NULL)
})
# Cleanup temporary tile files
unlink(tile_storage_dir, recursive = TRUE)
safe_log("Cleaned up temporary tile files")
return(final_file_path)
}
#' Load tile MAX rasters for a specific week (for per-tile analysis)
#'
#' @param week Week number
#' @param tile_dir Directory containing tile mosaics (week subdirectories)
#' @return List of tile rasters with tile_id metadata
#'
load_tile_mosaics <- function(week, tile_dir) {
week_dir <- file.path(tile_dir, paste0("week_", week))
if (!dir.exists(week_dir)) {
warning(paste("Tile directory not found:", week_dir))
return(NULL)
}
# Load all tile files
tile_files <- list.files(week_dir, pattern = "^tile_.*\\.tif$", full.names = TRUE)
if (length(tile_files) == 0) {
warning("No tile files found in:", week_dir)
return(NULL)
}
# Sort by tile ID
tile_numbers <- as.numeric(gsub(".*tile_(\\d+).*", "\\1", tile_files))
tile_files <- tile_files[order(tile_numbers)]
# Load rasters
tiles_list <- lapply(tile_files, terra::rast)
names(tiles_list) <- sprintf("tile_%03d", sort(tile_numbers))
safe_log(paste("Loaded", length(tiles_list), "tile mosaics for week", week))
return(tiles_list)
}

779
r_app/40_mosaic_creation_utils.R
View file

@ -1,779 +0,0 @@
# MOSAIC_CREATION_UTILS.R
# ======================
# Utility functions for creating weekly mosaics from daily satellite imagery.
# These functions support cloud cover assessment, date handling, and mosaic creation.
#' Detect whether a project uses tile-based or single-file mosaic approach (utility version)
#'
#' @param merged_final_tif_dir Directory containing merged_final_tif files
#' @return List with has_tiles (logical), detected_tiles (vector), total_files (count)
#'
detect_tile_structure_from_files <- function(merged_final_tif_dir) {
# Check if directory exists
if (!dir.exists(merged_final_tif_dir)) {
return(list(has_tiles = FALSE, detected_tiles = character(), total_files = 0))
}
# List all .tif files in merged_final_tif
tif_files <- list.files(merged_final_tif_dir, pattern = "\\.tif$", full.names = FALSE)
if (length(tif_files) == 0) {
return(list(has_tiles = FALSE, detected_tiles = character(), total_files = 0))
}
# Check if ANY file matches tile naming pattern: *_XX.tif (where XX is 2 digits)
# Tile pattern examples: 2025-11-27_00.tif, 2025-11-27_01.tif, week_50_2024_00.tif
tile_pattern <- "_(\\d{2})\\.tif$"
tile_files <- tif_files[grepl(tile_pattern, tif_files)]
has_tiles <- length(tile_files) > 0
return(list(
has_tiles = has_tiles,
detected_tiles = tile_files,
total_files = length(tif_files)
))
}
#' Generate a sequence of dates for processing
#'
#' @param end_date The end date for the sequence (Date object)
#' @param offset Number of days to look back from end_date
#' @return A list containing week number, year, and a sequence of dates for filtering
#'
# NOTE: date_list() is now in 00_common_utils.R - import from there
# This function was duplicated and has been consolidated
#' Create a weekly mosaic from available VRT files
#'
#' @param dates List from date_list() with date range info
#' @param field_boundaries Field boundaries for image cropping
#' @param daily_vrt_dir Directory containing VRT files
#' @param merged_final_dir Directory with merged final rasters
#' @param output_dir Output directory for weekly mosaics
#' @param file_name_tif Output filename for the mosaic
#' @param create_plots Whether to create visualization plots (default: TRUE)
#' @return The file path of the saved mosaic
#'
create_weekly_mosaic <- function(dates, field_boundaries, daily_vrt_dir,
merged_final_dir, output_dir, file_name_tif,
create_plots = FALSE) {
# NOTE: VRT files are legacy code - we no longer create or use them
# Get dates directly from the dates parameter instead
dates_to_check <- dates$days_filter
# Find final raster files for fallback
raster_files_final <- list.files(merged_final_dir, full.names = TRUE, pattern = "\\.tif$")
# Process the mosaic if we have dates to check
if (length(dates_to_check) > 0) {
safe_log("Processing dates, assessing cloud cover for mosaic creation")
# Calculate aggregated cloud cover statistics (returns data frame for image selection)
cloud_coverage_stats <- count_cloud_coverage(dates_to_check, merged_final_dir, field_boundaries)
# Create mosaic based on cloud cover assessment
mosaic <- create_mosaic(raster_files_final, cloud_coverage_stats, field_boundaries)
} else {
safe_log("No dates available for the date range, creating empty mosaic with NA values", "WARNING")
# Create empty mosaic if no files are available
if (length(raster_files_final) == 0) {
stop("No VRT files or final raster files available to create mosaic")
}
mosaic <- terra::rast(raster_files_final[1])
mosaic <- terra::setValues(mosaic, NA)
mosaic <- terra::crop(mosaic, field_boundaries, mask = TRUE)
names(mosaic) <- c("Red", "Green", "Blue", "NIR", "CI")
}
# Save the mosaic (without mask files to avoid breaking other scripts)
file_path <- save_mosaic(mosaic, output_dir, file_name_tif, create_plots, save_mask = FALSE)
safe_log(paste("Weekly mosaic processing completed for week", dates$week))
return(file_path)
}
#' Find VRT files within a date range
#'
#' @param vrt_directory Directory containing VRT files
#' @param dates List from date_list() function containing days_filter
#' @return Character vector of VRT file paths
#'
find_vrt_files <- function(vrt_directory, dates) {
# Get all VRT files in directory
# Note: vrt_directory is already a full/relative path from parameters_project.R
# Don't wrap it in here::here() again - that would create an incorrect path
vrt_files <- list.files(vrt_directory, full.names = TRUE)
if (length(vrt_files) == 0) {
warning("No VRT files found in directory: ", vrt_directory)
return(character(0))
}
# Filter files by dates
vrt_list <- purrr::map(dates$days_filter, ~ vrt_files[grepl(pattern = .x, x = vrt_files)]) %>%
purrr::compact() %>%
purrr::flatten_chr()
# Log results
safe_log(paste("Found", length(vrt_list), "VRT files for the date range"))
return(vrt_list)
}
#' Count missing pixels (clouds) in rasters - per field analysis using actual TIF files
#'
#' @param dates_to_check Character vector of dates in YYYY-MM-DD format to check for cloud coverage
#' @param merged_final_dir Directory containing the actual TIF files (e.g., merged_final_tif)
#' @param field_boundaries Field boundaries (sf object) for calculating cloud cover over field areas only
#' @return Data frame with aggregated cloud statistics for each TIF file (used for mosaic selection)
#'
count_cloud_coverage <- function(dates_to_check, merged_final_dir = NULL, field_boundaries = NULL) {
if (length(dates_to_check) == 0) {
warning("No dates provided for cloud coverage calculation")
return(NULL)
}
tryCatch({
# Build list of actual TIF files from dates
# TIF filenames are like "2025-12-18.tif"
tif_files <- paste0(here::here(merged_final_dir), "/", dates_to_check, ".tif")
# Check which TIF files exist
tif_exist <- file.exists(tif_files)
if (!any(tif_exist)) {
warning("No TIF files found in directory: ", merged_final_dir)
return(NULL)
}
tif_files <- tif_files[tif_exist]
safe_log(paste("Found", length(tif_files), "TIF files for cloud coverage assessment"))
# Initialize list to store aggregated results
aggregated_results <- list()
# Process each TIF file
for (tif_idx in seq_along(tif_files)) {
tif_file <- tif_files[tif_idx]
tryCatch({
# Load the TIF file (typically has 5 bands: R, G, B, NIR, CI)
current_raster <- terra::rast(tif_file)
# Extract the CI band (last band)
ci_band <- current_raster[[terra::nlyr(current_raster)]]
# Create a unique field mask for THIS raster's extent
# This handles cases where rasters have different extents due to missing data
total_notna <- NA_real_
total_pixels <- NA_real_
if (!is.null(field_boundaries)) {
tryCatch({
# Create mask specific to this raster's grid
field_mask <- terra::rasterize(field_boundaries, ci_band, field = 1)
# Count pixels within field boundaries (for this specific raster)
total_pixels <- terra::global(field_mask, fun = "notNA")$notNA
# Cloud coverage calculated only over field areas
ci_field_masked <- terra::mask(ci_band, field_mask, maskvalue = NA)
total_notna <- terra::global(ci_field_masked, fun = "notNA")$notNA
}, error = function(e) {
# If field mask creation fails, fall back to entire raster
safe_log(paste("Could not create field mask for", basename(tif_file), ":", e$message), "WARNING")
})
}
# If field mask failed, use entire raster
if (is.na(total_notna)) {
total_notna <- terra::global(ci_band, fun = "notNA")$notNA
total_pixels <- terra::ncell(ci_band)
}
# Calculate cloud coverage percentage (missing = clouds)
missing_pct <- round(100 - ((total_notna / total_pixels) * 100))
aggregated_results[[tif_idx]] <- data.frame(
filename = basename(tif_file),
notNA = total_notna,
total_pixels = total_pixels,
missing_pixels_percentage = missing_pct,
thres_5perc = as.integer(missing_pct < 5),
thres_40perc = as.integer(missing_pct < 45),
stringsAsFactors = FALSE
)
}, error = function(e) {
safe_log(paste("Error processing TIF", basename(tif_file), ":", e$message), "WARNING")
aggregated_results[[tif_idx]] <<- data.frame(
filename = basename(tif_file),
notNA = NA_real_,
total_pixels = NA_real_,
missing_pixels_percentage = 100,
thres_5perc = 0,
thres_40perc = 0,
stringsAsFactors = FALSE
)
})
}
# Combine all aggregated results
aggregated_df <- if (length(aggregated_results) > 0) {
do.call(rbind, aggregated_results)
} else {
data.frame()
}
# Log results
safe_log(paste("Cloud coverage assessment completed for", length(dates_to_check), "dates"))
# Return aggregated data only
return(aggregated_df)
}, error = function(e) {
warning("Error in cloud coverage calculation: ", e$message)
return(NULL)
})
}
#' Create a mosaic from merged_final_tif files based on cloud coverage
#'
#' @param tif_files List of processed TIF files (5 bands: R, G, B, NIR, CI)
#' @param cloud_coverage_stats Cloud coverage statistics from count_cloud_coverage()
#' @param field_boundaries Field boundaries for masking (optional)
#' @return A SpatRaster object with 5 bands (Red, Green, Blue, NIR, CI)
#'
create_mosaic <- function(tif_files, cloud_coverage_stats, field_boundaries = NULL) {
# If no TIF files, return NULL
if (length(tif_files) == 0) {
safe_log("No TIF files available for mosaic creation", "ERROR")
return(NULL)
}
# Validate cloud coverage stats
mosaic_type <- "Unknown" # Track what type of mosaic is being created
if (is.null(cloud_coverage_stats) || nrow(cloud_coverage_stats) == 0) {
safe_log("No cloud coverage statistics available, using all files", "WARNING")
rasters_to_use <- tif_files
mosaic_type <- paste("all", length(tif_files), "available images")
} else {
# Determine best rasters to use based on cloud coverage thresholds
# Count how many images meet each threshold
num_5perc <- sum(cloud_coverage_stats$thres_5perc, na.rm = TRUE)
num_40perc <- sum(cloud_coverage_stats$thres_40perc, na.rm = TRUE)
if (num_5perc > 1) {
# Multiple images with <5% cloud coverage
safe_log(paste("Creating max composite from", num_5perc, "cloud-free images (<5% clouds)"))
mosaic_type <- paste(num_5perc, "cloud-free images (<5% clouds)")
best_coverage <- which(cloud_coverage_stats$thres_5perc > 0)
} else if (num_5perc == 1) {
# Single image with <5% cloud coverage
safe_log("Using single cloud-free image (<5% clouds)")
mosaic_type <- "single cloud-free image (<5% clouds)"
best_coverage <- which(cloud_coverage_stats$thres_5perc > 0)
} else if (num_40perc > 1) {
# Multiple images with <40% cloud coverage
safe_log(paste("Creating max composite from", num_40perc, "partially cloudy images (<40% clouds)"), "WARNING")
mosaic_type <- paste(num_40perc, "partially cloudy images (<40% clouds)")
best_coverage <- which(cloud_coverage_stats$thres_40perc > 0)
} else if (num_40perc == 1) {
# Single image with <40% cloud coverage
safe_log("Using single partially cloudy image (<40% clouds)", "WARNING")
mosaic_type <- "single partially cloudy image (<40% clouds)"
best_coverage <- which(cloud_coverage_stats$thres_40perc > 0)
} else {
# No cloud-free images available
safe_log("No cloud-free images available, using all images", "WARNING")
mosaic_type <- paste("all", nrow(cloud_coverage_stats), "available images")
best_coverage <- seq_len(nrow(cloud_coverage_stats))
}
# Get filenames of best-coverage images
# Match by full filename from cloud stats to TIF files
rasters_to_use <- character()
for (idx in best_coverage) {
# Get the full filename from cloud coverage stats
cc_filename <- cloud_coverage_stats$filename[idx]
# Find matching TIF file by full filename
matching_tif <- NULL
for (tif_file in tif_files) {
tif_basename <- basename(tif_file)
if (tif_basename == cc_filename) {
matching_tif <- tif_file
break
}
}
if (!is.null(matching_tif)) {
rasters_to_use <- c(rasters_to_use, matching_tif)
} else {
safe_log(paste("Warning: Could not find TIF file matching cloud stats entry:", cc_filename), "WARNING")
}
}
if (length(rasters_to_use) == 0) {
safe_log("Could not match cloud coverage stats to TIF files, using all files", "WARNING")
rasters_to_use <- tif_files
mosaic_type <- paste("all", length(tif_files), "available images")
}
}
# Load and mosaic the selected rasters
if (length(rasters_to_use) == 1) {
# Single file - just load it
safe_log(paste("Using single image for mosaic:", basename(rasters_to_use)))
mosaic <- terra::rast(rasters_to_use[1])
} else {
# Multiple files - merge handles different extents/grids automatically
safe_log(paste("Creating mosaic from", length(rasters_to_use), "images"))
# Load all rasters with error handling - only keep successful loads
all_rasters <- Filter(Negate(is.null), lapply(rasters_to_use, function(f) {
tryCatch({
terra::rast(f)
}, error = function(e) {
safe_log(paste("Warning: Could not load", basename(f), ":", e$message), "WARNING")
NULL # Return NULL on error, will be filtered out
})
}))
# Check what we loaded
safe_log(paste("Loaded", length(all_rasters), "valid rasters from", length(rasters_to_use), "files"))
if (length(all_rasters) == 0) {
safe_log("No valid rasters to merge", "WARNING")
return(NULL)
}
# Merge all rasters (handles different extents and grids automatically)
if (length(all_rasters) == 1) {
mosaic <- all_rasters[[1]]
safe_log("Using single raster after filtering")
} else {
# Create max composite: take maximum value at each pixel across all dates
# This skips clouds (low/zero CI values) in favor of clear pixels from other dates
mosaic <- tryCatch({
safe_log(paste("Creating max composite from", length(all_rasters), "images to fill clouds"))
# Check if all rasters have identical grids (extent and resolution)
# This is likely for per-tile mosaics from the same tiling scheme
reference_raster <- all_rasters[[1]]
ref_ext <- terra::ext(reference_raster)
ref_res <- terra::res(reference_raster)
grids_match <- all(sapply(all_rasters[-1], function(r) {
isTRUE(all.equal(terra::ext(r), ref_ext, tolerance = 1e-6)) &&
isTRUE(all.equal(terra::res(r), ref_res, tolerance = 1e-6))
}))
if (grids_match) {
# All rasters have matching grids - no cropping/resampling needed!
safe_log("All rasters have identical grids - stacking directly for max composite")
raster_collection <- terra::sprc(all_rasters)
max_mosaic <- terra::mosaic(raster_collection, fun = "max")
} else {
# Grids don't match - need to crop and resample
safe_log("Rasters have different grids - cropping and resampling to common extent")
# Get extent from field boundaries if available, otherwise use raster union
if (!is.null(field_boundaries)) {
crop_extent <- terra::ext(field_boundaries)
safe_log("Using field boundaries extent for consistent area across all dates")
} else {
# Use union of all extents (covers all data)
crop_extent <- terra::ext(all_rasters[[1]])
for (i in 2:length(all_rasters)) {
crop_extent <- terra::union(crop_extent, terra::ext(all_rasters[[i]]))
}
safe_log("Using raster union extent")
}
# Crop all rasters to common extent
cropped_rasters <- lapply(all_rasters, function(r) {
terra::crop(r, crop_extent)
})
# Resample all cropped rasters to match the first one's grid
reference_grid <- cropped_rasters[[1]]
aligned_rasters <- lapply(cropped_rasters, function(r) {
if (isTRUE(all.equal(terra::ext(r), terra::ext(reference_grid), tolerance = 1e-6)) &&
isTRUE(all.equal(terra::res(r), terra::res(reference_grid), tolerance = 1e-6))) {
return(r) # Already aligned
}
terra::resample(r, reference_grid, method = "near")
})
# Create max composite using mosaic on aligned rasters
raster_collection <- terra::sprc(aligned_rasters)
max_mosaic <- terra::mosaic(raster_collection, fun = "max")
}
max_mosaic
}, error = function(e) {
safe_log(paste("Max composite creation failed:", e$message), "WARNING")
safe_log("Using first raster only as fallback")
all_rasters[[1]]
})
safe_log(paste("Max composite created - taking clearest pixel at each location"))
}
# Ensure we have exactly the required bands: Red, Green, Blue, NIR, CI
required_bands <- c("Red", "Green", "Blue", "NIR", "CI")
available_bands <- names(mosaic)
# Check if all required bands are present
if (!all(required_bands %in% available_bands)) {
safe_log(paste("Warning: Not all required bands found. Available:", paste(available_bands, collapse = ", ")), "WARNING")
}
# Select only the required bands in the correct order
if (all(required_bands %in% available_bands)) {
mosaic <- mosaic[[required_bands]]
safe_log("Selected Red, Green, Blue, NIR, CI bands")
} else {
safe_log(paste("Warning: mosaic has", terra::nlyr(mosaic), "bands, expected 5 (R, G, B, NIR, CI)"), "WARNING")
if (terra::nlyr(mosaic) > 5) {
# Keep only first 5 bands as fallback
mosaic <- terra::subset(mosaic, 1:5)
safe_log("Keeping only first 5 bands as fallback")
}
}
}
# Crop/mask to field boundaries if provided
if (!is.null(field_boundaries)) {
tryCatch({
mosaic <- terra::crop(mosaic, field_boundaries, mask = TRUE)
safe_log("Mosaic cropped to field boundaries")
}, error = function(e) {
safe_log(paste("Could not crop to field boundaries:", e$message), "WARNING")
# Return uncropped mosaic
})
}
# Log final mosaic summary
safe_log(paste("✓ Mosaic created from", mosaic_type, "-", terra::nlyr(mosaic),
"bands,", nrow(mosaic), "x", ncol(mosaic), "pixels"))
return(mosaic)
}
#' Save a mosaic raster to disk
#'
#' @param mosaic_raster A SpatRaster object to save
#' @param output_dir Directory to save the output
#' @param file_name Filename for the output raster
#' @param plot_result Whether to create visualizations (default: FALSE)
#' @param save_mask Whether to save cloud masks separately (default: FALSE)
#' @return The file path of the saved raster
#'
save_mosaic <- function(mosaic_raster, output_dir, file_name, plot_result = FALSE, save_mask = FALSE) {
# Validate input
if (is.null(mosaic_raster)) {
stop("No mosaic raster provided to save")
}
# Create output directory if it doesn't exist
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
# Create full file path - use file.path() since output_dir may be absolute path
# Ensure file_name has .tif extension
if (!grepl("\\.tif$|\\.TIF$", file_name)) {
file_name <- paste0(file_name, ".tif")
}
file_path <- file.path(output_dir, file_name)
# Get cloud mask if it exists
cloud_mask <- attr(mosaic_raster, "cloud_mask")
# Save raster
terra::writeRaster(mosaic_raster, file_path, overwrite = TRUE)
# Save cloud mask if available and requested
if (!is.null(cloud_mask) && save_mask) {
# Create mask filename by adding _mask before extension
mask_file_name <- gsub("\\.(tif|TIF)$", "_mask.\\1", file_name)
mask_file_path <- here::here(output_dir, mask_file_name)
# Save the mask
terra::writeRaster(cloud_mask, mask_file_path, overwrite = TRUE)
safe_log(paste("Cloud/shadow mask saved to:", mask_file_path))
} else if (!is.null(cloud_mask)) {
safe_log("Cloud mask available but not saved (save_mask = FALSE)")
}
# Create plots if requested
if (plot_result) {
# Plot the CI band
if ("CI" %in% names(mosaic_raster)) {
terra::plot(mosaic_raster$CI, main = paste("CI map", file_name))
}
# Plot RGB image
if (all(c("Red", "Green", "Blue") %in% names(mosaic_raster))) {
terra::plotRGB(mosaic_raster, main = paste("RGB map", file_name))
}
# Plot cloud mask if available
if (!is.null(cloud_mask)) {
terra::plot(cloud_mask, main = paste("Cloud/shadow mask", file_name),
col = c("red", "green"))
}
# If we have both RGB and cloud mask, create a side-by-side comparison
if (all(c("Red", "Green", "Blue") %in% names(mosaic_raster)) && !is.null(cloud_mask)) {
old_par <- par(mfrow = c(1, 2))
terra::plotRGB(mosaic_raster, main = "RGB Image")
# Create a colored mask for visualization (red = cloud/shadow, green = clear)
mask_plot <- cloud_mask
terra::plot(mask_plot, main = "Cloud/Shadow Mask", col = c("red", "green"))
par(old_par)
}
}
# Log save completion
safe_log(paste("Mosaic saved to:", file_path))
return(file_path)
}
#' Create weekly mosaic from pre-split tiles with MAX aggregation
#'
#' This function processes tiles created by Script 01 and processed by Script 02.
#' For each of the 25 tiles independently:
#' 1. Collects that tile from all dates in the range
#' 2. Calculates cloud coverage per date
#' 3. Uses create_mosaic logic to select best dates (cloud-clean preferred)
#' 4. Creates MAX composite for that tile
#' 5. Saves to weekly_tile_max/tile_XX.tif
#'
#' Input: merged_final_tif/[DATE]/[TILE_01.tif, TILE_02.tif, ..., TILE_25.tif]
#' Output: weekly_tile_max/tile_01.tif through tile_25.tif (25 weekly MAX tiles)
#'
#' @param dates List from date_list() containing days_filter vector
#' @param merged_final_dir Directory containing processed tiles (merged_final_tif)
#' @param tile_output_dir Directory to save weekly MAX tiles (weekly_tile_max)
#' @param field_boundaries Field boundaries for cloud coverage calculation (optional)
#' @return List of paths to created tile files
#'
create_weekly_mosaic_from_tiles <- function(dates, merged_final_dir, tile_output_dir, field_boundaries = NULL) {
safe_log("Starting per-tile mosaic creation with cloud-based date selection...")
# Create output directory if needed
dir.create(tile_output_dir, recursive = TRUE, showWarnings = FALSE)
# Step 1: Discover all tiles from all dates and group by tile ID
tile_groups <- list() # Structure: tile_groups$tile_01 = list of files for that tile across dates
for (date in dates$days_filter) {
date_dir <- file.path(merged_final_dir, date)
if (!dir.exists(date_dir)) {
safe_log(paste(" Skipping date:", date, "- directory not found"), "WARNING")
next
}
tile_files <- list.files(date_dir, pattern = "\\.tif$", full.names = TRUE)
if (length(tile_files) == 0) {
safe_log(paste(" No tiles found for date:", date), "WARNING")
next
}
# Extract tile ID from each filename (e.g., "2026-01-02_01.tif" → tile ID is "01")
for (tile_file in tile_files) {
# Extract tile number from filename
tile_basename <- basename(tile_file)
tile_id <- gsub(".*_([0-9]+)\\.tif", "\\1", tile_basename)
if (!tile_id %in% names(tile_groups)) {
tile_groups[[tile_id]] <- list()
}
tile_groups[[tile_id]][[length(tile_groups[[tile_id]]) + 1]] <- tile_file
}
}
if (length(tile_groups) == 0) {
stop("No tiles found in date range")
}
safe_log(paste("Found", length(tile_groups), "unique tiles across all dates"))
# Step 2: Process each tile independently
created_tiles <- character()
for (tile_id in names(tile_groups)) {
tile_files_for_this_id <- unlist(tile_groups[[tile_id]])
safe_log(paste("Processing tile", tile_id, "with", length(tile_files_for_this_id), "dates"))
# Step 2a: Calculate cloud coverage for this tile across all dates
cloud_stats_this_tile <- tryCatch({
count_cloud_coverage_for_tile(
tile_files = tile_files_for_this_id,
field_boundaries = field_boundaries
)
}, error = function(e) {
safe_log(paste(" Error calculating cloud coverage for tile", tile_id, "-", e$message), "WARNING")
NULL
})
if (is.null(cloud_stats_this_tile) || nrow(cloud_stats_this_tile) == 0) {
safe_log(paste(" No valid cloud stats for tile", tile_id, "- using all available dates"), "WARNING")
cloud_stats_this_tile <- NULL
}
# Step 2b: Create MAX mosaic for this tile using create_mosaic logic
tile_mosaic <- tryCatch({
create_mosaic(
tif_files = tile_files_for_this_id,
cloud_coverage_stats = cloud_stats_this_tile,
field_boundaries = NULL # Don't crop individual tiles
)
}, error = function(e) {
safe_log(paste(" Error creating mosaic for tile", tile_id, "-", e$message), "WARNING")
NULL
})
if (is.null(tile_mosaic)) {
safe_log(paste(" Failed to create mosaic for tile", tile_id, "- skipping"), "WARNING")
next
}
# DEBUG: Check mosaic content before saving
safe_log(paste(" DEBUG: Mosaic tile", tile_id, "dimensions:", nrow(tile_mosaic), "x", ncol(tile_mosaic)))
safe_log(paste(" DEBUG: Mosaic tile", tile_id, "bands:", terra::nlyr(tile_mosaic)))
# Check first band values
band1 <- tile_mosaic[[1]]
band1_min <- terra::global(band1, fun = "min", na.rm = TRUE)$min
band1_max <- terra::global(band1, fun = "max", na.rm = TRUE)$max
safe_log(paste(" DEBUG: Band 1 MIN=", round(band1_min, 2), "MAX=", round(band1_max, 2)))
# Step 2c: Save this tile's weekly MAX mosaic
# Filename format: week_WW_YYYY_TT.tif (e.g., week_02_2026_01.tif for week 2, 2026, tile 1)
tile_filename <- paste0("week_", sprintf("%02d", dates$week), "_", dates$year, "_",
sprintf("%02d", as.integer(tile_id)), ".tif")
tile_output_path <- file.path(tile_output_dir, tile_filename)
tryCatch({
terra::writeRaster(tile_mosaic, tile_output_path, overwrite = TRUE)
safe_log(paste(" ✓ Saved tile", tile_id, "to:", tile_filename))
created_tiles <- c(created_tiles, tile_output_path)
}, error = function(e) {
safe_log(paste(" Error saving tile", tile_id, "-", e$message), "WARNING")
})
}
safe_log(paste("✓ Created", length(created_tiles), "weekly MAX tiles in", tile_output_dir))
return(created_tiles)
}
#' Calculate cloud coverage for a single tile across multiple dates
#'
#' Helper function for per-tile cloud assessment.
#' Takes tile files from different dates and calculates cloud coverage for each.
#' Cloud coverage is calculated ONLY within field boundaries, so total_pixels
#' varies per tile based on which fields are present in that tile area.
#'
#' @param tile_files Character vector of tile file paths from different dates
#' @param field_boundaries Field boundaries for analysis (required for per-field counting)
#' @return Data frame with cloud stats for each date/tile
#'
count_cloud_coverage_for_tile <- function(tile_files, field_boundaries = NULL) {
if (length(tile_files) == 0) {
warning("No tile files provided for cloud coverage calculation")
return(NULL)
}
aggregated_results <- list()
for (idx in seq_along(tile_files)) {
tile_file <- tile_files[idx]
tryCatch({
# Load the tile
current_raster <- terra::rast(tile_file)
# Extract the CI band (last band in 5-band output)
ci_band <- current_raster[[terra::nlyr(current_raster)]]
# Calculate cloud coverage within field boundaries
if (!is.null(field_boundaries)) {
# Create a reference raster template (same extent/resolution as ci_band, but independent of its data)
# This ensures field_mask shows the potential field area even if ci_band is entirely cloudy (all NA)
ref_template <- terra::rast(terra::ext(ci_band), resolution = terra::res(ci_band),
crs = terra::crs(ci_band))
terra::values(ref_template) <- 1 # Fill entire raster with 1s
# Crop and mask to field boundaries: keeps 1 inside fields, NA outside
# This is independent of ci_band's data - represents the potential field area
field_mask <- terra::crop(ref_template, field_boundaries, mask = TRUE)
# Count total potential field pixels from the mask (independent of clouds)
total_pixels <- terra::global(field_mask, fun = "notNA")$notNA
# Now crop and mask CI band to field boundaries to count actual valid (non-cloud) pixels
ci_masked <- terra::crop(ci_band, field_boundaries, mask = TRUE)
total_notna <- terra::global(ci_masked, fun = "notNA")$notNA
} else {
# If no field boundaries provided, use entire tile
total_notna <- terra::global(ci_band, fun = "notNA")$notNA
total_pixels <- terra::ncell(ci_band)
}
# Cloud coverage percentage (missing = clouds)
missing_pct <- round(100 - ((total_notna / total_pixels) * 100))
aggregated_results[[idx]] <- data.frame(
filename = basename(tile_file), # Keep full filename: 2026-01-07_03.tif
notNA = total_notna,
total_pixels = total_pixels,
missing_pixels_percentage = missing_pct,
thres_5perc = as.integer(missing_pct < 5),
thres_40perc = as.integer(missing_pct < 45),
stringsAsFactors = FALSE
)
}, error = function(e) {
safe_log(paste("Error processing tile:", basename(tile_file), "-", e$message), "WARNING")
aggregated_results[[idx]] <<- data.frame(
filename = tile_file,
notNA = NA_real_,
total_pixels = NA_real_,
missing_pixels_percentage = 100,
thres_5perc = 0,
thres_40perc = 0,
stringsAsFactors = FALSE
)
})
}
# Combine results
aggregated_df <- if (length(aggregated_results) > 0) {
do.call(rbind, aggregated_results)
} else {
data.frame()
}
return(aggregated_df)
}

32
r_app/80_calculate_kpis.R
View file

@ -134,17 +134,35 @@ tryCatch({
stop("Error loading 00_common_utils.R: ", e$message)
})
# ============================================================================
# LOAD CLIENT-AWARE UTILITIES
# ============================================================================
# All clients use the common utilities (shared statistical functions, reporting)
tryCatch({
source(here("r_app", "80_weekly_stats_utils.R"))
source(here("r_app", "80_utils_common.R"))
}, error = function(e) {
stop("Error loading 80_weekly_stats_utils.R: ", e$message)
stop("Error loading 80_utils_common.R: ", e$message)
})
tryCatch({
source(here("r_app", "80_report_building_utils.R"))
}, error = function(e) {
stop("Error loading 80_report_building_utils.R: ", e$message)
})
# Client-specific utilities based on client_config$script_90_compatible
# script_90_compatible = TRUE -> AURA workflow (6 KPIs)
# script_90_compatible = FALSE -> CANE_SUPPLY workflow (weekly stats + basic reporting)
if (client_config$script_90_compatible) {
message("Loading AURA client utilities (80_utils_agronomic_support.R)...")
tryCatch({
source(here("r_app", "80_utils_agronomic_support.R"))
}, error = function(e) {
stop("Error loading 80_utils_agronomic_support.R: ", e$message)
})
} else {
message("Loading CANE_SUPPLY client utilities (80_utils_cane_supply.R)...")
tryCatch({
source(here("r_app", "80_utils_cane_supply.R"))
}, error = function(e) {
stop("Error loading 80_utils_cane_supply.R: ", e$message)
})
}
# ============================================================================
# PHASE AND STATUS TRIGGER DEFINITIONS

1508
r_app/80_kpi_utils.R
View file

File diff suppressed because it is too large Load diff

258
r_app/80_report_building_utils.R
View file

@ -1,258 +0,0 @@
# 80_REPORT_BUILDING_UTILS.R
# ============================================================================
# UTILITY FUNCTIONS FOR REPORT GENERATION AND EXCEL/CSV EXPORT
#
# This file contains reusable functions for:
# - Field analysis summary generation
# - Excel/CSV/RDS export functionality
# - Farm-level KPI aggregation and summary
# - Tile-based KPI extraction (alternative calculation method)
#
# Used by: 80_calculate_kpis.R, run_full_pipeline.R, other reporting scripts
# ============================================================================
# ============================================================================
# SUMMARY GENERATION
# ============================================================================
generate_field_analysis_summary <- function(field_df) {
message("Generating summary statistics...")
total_acreage <- sum(field_df$Acreage, na.rm = TRUE)
germination_acreage <- sum(field_df$Acreage[field_df$Phase == "Germination"], na.rm = TRUE)
tillering_acreage <- sum(field_df$Acreage[field_df$Phase == "Tillering"], na.rm = TRUE)
grand_growth_acreage <- sum(field_df$Acreage[field_df$Phase == "Grand Growth"], na.rm = TRUE)
maturation_acreage <- sum(field_df$Acreage[field_df$Phase == "Maturation"], na.rm = TRUE)
unknown_phase_acreage <- sum(field_df$Acreage[field_df$Phase == "Unknown"], na.rm = TRUE)
harvest_ready_acreage <- sum(field_df$Acreage[field_df$Status_trigger == "harvest_ready"], na.rm = TRUE)
stress_acreage <- sum(field_df$Acreage[field_df$Status_trigger == "stress_detected_whole_field"], na.rm = TRUE)
recovery_acreage <- sum(field_df$Acreage[field_df$Status_trigger == "strong_recovery"], na.rm = TRUE)
growth_on_track_acreage <- sum(field_df$Acreage[field_df$Status_trigger == "growth_on_track"], na.rm = TRUE)
germination_complete_acreage <- sum(field_df$Acreage[field_df$Status_trigger == "germination_complete"], na.rm = TRUE)
germination_started_acreage <- sum(field_df$Acreage[field_df$Status_trigger == "germination_started"], na.rm = TRUE)
no_trigger_acreage <- sum(field_df$Acreage[is.na(field_df$Status_trigger)], na.rm = TRUE)
clear_fields <- sum(field_df$Cloud_category == "Clear view", na.rm = TRUE)
partial_fields <- sum(field_df$Cloud_category == "Partial coverage", na.rm = TRUE)
no_image_fields <- sum(field_df$Cloud_category == "No image available", na.rm = TRUE)
total_fields <- nrow(field_df)
clear_acreage <- sum(field_df$Acreage[field_df$Cloud_category == "Clear view"], na.rm = TRUE)
partial_acreage <- sum(field_df$Acreage[field_df$Cloud_category == "Partial coverage"], na.rm = TRUE)
no_image_acreage <- sum(field_df$Acreage[field_df$Cloud_category == "No image available"], na.rm = TRUE)
summary_df <- data.frame(
Category = c(
"--- PHASE DISTRIBUTION ---",
"Germination",
"Tillering",
"Grand Growth",
"Maturation",
"Unknown phase",
"--- STATUS TRIGGERS ---",
"Harvest ready",
"Stress detected",
"Strong recovery",
"Growth on track",
"Germination complete",
"Germination started",
"No trigger",
"--- CLOUD COVERAGE (FIELDS) ---",
"Clear view",
"Partial coverage",
"No image available",
"--- CLOUD COVERAGE (ACREAGE) ---",
"Clear view",
"Partial coverage",
"No image available",
"--- TOTAL ---",
"Total Acreage"
),
Acreage = c(
NA,
round(germination_acreage, 2),
round(tillering_acreage, 2),
round(grand_growth_acreage, 2),
round(maturation_acreage, 2),
round(unknown_phase_acreage, 2),
NA,
round(harvest_ready_acreage, 2),
round(stress_acreage, 2),
round(recovery_acreage, 2),
round(growth_on_track_acreage, 2),
round(germination_complete_acreage, 2),
round(germination_started_acreage, 2),
round(no_trigger_acreage, 2),
NA,
paste0(clear_fields, " fields"),
paste0(partial_fields, " fields"),
paste0(no_image_fields, " fields"),
NA,
round(clear_acreage, 2),
round(partial_acreage, 2),
round(no_image_acreage, 2),
NA,
round(total_acreage, 2)
),
stringsAsFactors = FALSE
)
return(summary_df)
}
# ============================================================================
# EXPORT FUNCTIONS
# ============================================================================
export_field_analysis_excel <- function(field_df, summary_df, project_dir, current_week, year, reports_dir) {
message("Exporting per-field analysis to Excel, CSV, and RDS...")
field_df_rounded <- field_df %>%
mutate(across(where(is.numeric), ~ round(., 2)))
# Handle NULL summary_df
summary_df_rounded <- if (!is.null(summary_df)) {
summary_df %>%
mutate(across(where(is.numeric), ~ round(., 2)))
} else {
NULL
}
output_subdir <- file.path(reports_dir, "kpis", "field_analysis")
if (!dir.exists(output_subdir)) {
dir.create(output_subdir, recursive = TRUE)
}
excel_filename <- paste0(project_dir, "_field_analysis_week", sprintf("%02d_%d", current_week, year), ".xlsx")
excel_path <- file.path(output_subdir, excel_filename)
excel_path <- normalizePath(excel_path, winslash = "\\", mustWork = FALSE)
# Build sheets list dynamically
sheets <- list(
"Field Data" = field_df_rounded
)
if (!is.null(summary_df_rounded)) {
sheets[["Summary"]] <- summary_df_rounded
}
write_xlsx(sheets, excel_path)
message(paste("✓ Field analysis Excel exported to:", excel_path))
kpi_data <- list(
field_analysis = field_df_rounded,
field_analysis_summary = summary_df_rounded,
metadata = list(
current_week = current_week,
year = year,
project = project_dir,
created_at = Sys.time()
)
)
rds_filename <- paste0(project_dir, "_kpi_summary_tables_week", sprintf("%02d_%d", current_week, year), ".rds")
rds_path <- file.path(reports_dir, "kpis", rds_filename)
saveRDS(kpi_data, rds_path)
message(paste("✓ Field analysis RDS exported to:", rds_path))
csv_filename <- paste0(project_dir, "_field_analysis_week", sprintf("%02d_%d", current_week, year), ".csv")
csv_path <- file.path(output_subdir, csv_filename)
write_csv(field_df_rounded, csv_path)
message(paste("✓ Field analysis CSV exported to:", csv_path))
return(list(excel = excel_path, rds = rds_path, csv = csv_path))
}
# ============================================================================
# TILE-BASED KPI EXTRACTION (Alternative calculation method)
# ============================================================================
# [COMMENTED OUT / UNUSED - kept for reference]
# These functions provide tile-based extraction as an alternative to field_statistics approach
# Currently replaced by calculate_field_statistics() in 80_weekly_stats_utils.R
# Uncomment if parallel processing of tiles is needed in future
# calculate_field_kpis_from_tiles <- function(tile_dir, week_num, year, field_boundaries_sf, tile_grid) {
# message("Calculating field-level KPI statistics from tiles...")
#
# tile_pattern <- sprintf("week_%02d_%d_([0-9]{2})\\.tif", week_num, year)
# tile_files <- list.files(tile_dir, pattern = tile_pattern, full.names = TRUE)
#
# if (length(tile_files) == 0) {
# message("No tiles found for week", week_num, year)
# return(NULL)
# }
#
# message(paste("Processing", length(tile_files), "tiles in parallel..."))
#
# field_kpi_list <- furrr::future_map(
# tile_files,
# ~ process_single_kpi_tile(
# tile_file = .,
# field_boundaries_sf = field_boundaries_sf,
# tile_grid = tile_grid
# ),
# .progress = TRUE,
# .options = furrr::furrr_options(seed = TRUE)
# )
#
# field_kpi_stats <- dplyr::bind_rows(field_kpi_list)
#
# if (nrow(field_kpi_stats) == 0) {
# message(" No KPI data extracted from tiles")
# return(NULL)
# }
#
# message(paste(" Extracted KPI stats for", length(unique(field_kpi_stats$field)), "unique fields"))
# return(field_kpi_stats)
# }
# process_single_kpi_tile <- function(tile_file, field_boundaries_sf, tile_grid) {
# # Helper function for calculate_field_kpis_from_tiles
# tryCatch({
# tile_basename <- basename(tile_file)
# tile_raster <- terra::rast(tile_file)
# ci_band <- tile_raster[[1]]
#
# field_bbox <- sf::st_bbox(field_boundaries_sf)
# ci_cropped <- terra::crop(ci_band, terra::ext(field_bbox), snap = "out")
#
# extracted_vals <- terra::extract(ci_cropped, field_boundaries_sf, fun = "mean", na.rm = TRUE)
#
# tile_results <- data.frame()
# tile_id_match <- as.numeric(sub(".*_(\\d{2})\\.tif$", "\\1", tile_basename))
#
# for (field_idx in seq_len(nrow(field_boundaries_sf))) {
# field_id <- field_boundaries_sf$field[field_idx]
# mean_ci <- extracted_vals[field_idx, 2]
#
# if (is.na(mean_ci)) {
# next
# }
#
# tile_results <- rbind(tile_results, data.frame(
# field = field_id,
# tile_id = tile_id_match,
# tile_file = tile_basename,
# mean_ci = round(mean_ci, 4),
# stringsAsFactors = FALSE
# ))
# }
#
# return(tile_results)
#
# }, error = function(e) {
# message(paste(" Warning: Error processing tile", basename(tile_file), ":", e$message))
# return(data.frame())
# })
# }
# calculate_and_export_farm_kpis <- function(report_date, project_dir, field_boundaries_sf,
# harvesting_data, cumulative_CI_vals_dir,
# weekly_CI_mosaic, reports_dir, current_week, year,
# tile_grid, use_tile_mosaic = FALSE, tile_grid_size = "5x5") {
# # Farm-level KPI calculation using tile-based extraction (alternative approach)
# # [Implementation kept as reference for alternative calculation method]
# }

641
r_app/80_utils_agronomic_support.R Normal file
View file

@ -0,0 +1,641 @@
# 80_UTILS_AGRONOMIC_SUPPORT.R
# ============================================================================
# AURA-SPECIFIC KPI UTILITIES (SCRIPT 80 - CLIENT TYPE: agronomic_support)
#
# Contains all 6 AURA KPI calculation functions and helpers:
# - Field uniformity KPI (CV-based, spatial autocorrelation)
# - Area change KPI (week-over-week CI changes)
# - TCH forecasted KPI (tonnage projections from harvest data)
# - Growth decline KPI (trend analysis)
# - Weed presence KPI (field fragmentation detection)
# - Gap filling KPI (interpolation quality)
# - KPI reporting (summary tables, field details, text interpretation)
# - KPI export (Excel, RDS, data export)
#
# Orchestrator: calculate_all_kpis()
# Dependencies: 00_common_utils.R (safe_log), sourced from common
# Used by: 80_calculate_kpis.R (when client_type == "agronomic_support")
# ============================================================================
library(terra)
library(sf)
library(dplyr)
library(tidyr)
library(readxl)
library(writexl)
library(spdep)
library(caret)
library(CAST)
# ============================================================================
# SHARED HELPER FUNCTIONS (NOW IN 80_UTILS_COMMON.R)
# ============================================================================
# The following helper functions have been moved to 80_utils_common.R:
# - calculate_cv()
# - calculate_change_percentages()
# - calculate_spatial_autocorrelation()
# - extract_ci_values()
# - calculate_week_numbers()
# - load_field_ci_raster()
# - load_weekly_ci_mosaic()
# - prepare_predictions()
#
# These are now sourced from common utils and shared by all client types.
# ============================================================================
#' Prepare harvest predictions and ensure proper alignment with field data
prepare_predictions <- function(harvest_model, field_data, scenario = "optimistic") {
if (is.null(harvest_model) || is.null(field_data)) {
return(NULL)
}
tryCatch({
scenario_factor <- switch(scenario,
"pessimistic" = 0.85,
"realistic" = 1.0,
"optimistic" = 1.15,
1.0)
predictions <- field_data %>%
mutate(tch_forecasted = field_data$mean_ci * scenario_factor)
return(predictions)
}, error = function(e) {
message(paste("Error preparing predictions:", e$message))
return(NULL)
})
}
# ============================================================================
# AURA KPI CALCULATION FUNCTIONS (6 KPIS)
# ============================================================================
#' KPI 1: Calculate field uniformity based on CV and spatial autocorrelation
#'
#' Measures how uniform crop development is across the field.
#' Low CV + high positive Moran's I = excellent uniformity
#'
#' @param ci_pixels_by_field List of CI pixel arrays for each field
#' @param field_boundaries_sf SF object with field geometries
#' @param ci_band Raster band with CI values
#'
#' @return Data frame with field_idx, cv_value, morans_i, uniformity_score, interpretation
calculate_field_uniformity_kpi <- function(ci_pixels_by_field, field_boundaries_sf, ci_band = NULL) {
result <- data.frame(
field_idx = integer(),
cv_value = numeric(),
morans_i = numeric(),
uniformity_score = numeric(),
interpretation = character(),
stringsAsFactors = FALSE
)
for (field_idx in seq_len(nrow(field_boundaries_sf))) {
ci_pixels <- ci_pixels_by_field[[field_idx]]
if (is.null(ci_pixels) || length(ci_pixels) == 0) {
result <- rbind(result, data.frame(
field_idx = field_idx,
cv_value = NA_real_,
morans_i = NA_real_,
uniformity_score = NA_real_,
interpretation = "No data",
stringsAsFactors = FALSE
))
next
}
cv_val <- calculate_cv(ci_pixels)
morans_i <- NA_real_
if (!is.null(ci_band)) {
morans_i <- calculate_spatial_autocorrelation(ci_pixels, field_boundaries_sf[field_idx, ])
}
# Normalize CV (0-1 scale, invert so lower CV = higher score)
cv_normalized <- min(cv_val / 0.3, 1) # 0.3 = threshold for CV
cv_score <- 1 - cv_normalized
# Normalize Moran's I (-1 to 1 scale, shift to 0-1)
morans_normalized <- if (!is.na(morans_i)) {
(morans_i + 1) / 2
} else {
0.5
}
uniformity_score <- 0.7 * cv_score + 0.3 * morans_normalized
# Interpretation
if (is.na(cv_val)) {
interpretation <- "No data"
} else if (cv_val < 0.08) {
interpretation <- "Excellent uniformity"
} else if (cv_val < 0.15) {
interpretation <- "Good uniformity"
} else if (cv_val < 0.25) {
interpretation <- "Acceptable uniformity"
} else if (cv_val < 0.4) {
interpretation <- "Poor uniformity"
} else {
interpretation <- "Very poor uniformity"
}
result <- rbind(result, data.frame(
field_idx = field_idx,
cv_value = cv_val,
morans_i = morans_i,
uniformity_score = round(uniformity_score, 3),
interpretation = interpretation,
stringsAsFactors = FALSE
))
}
return(result)
}
#' KPI 2: Calculate area change metric (week-over-week CI changes)
#'
#' Tracks the percentage change in CI between current and previous week
#'
#' @param current_stats Current week field statistics (from extract_field_statistics_from_ci)
#' @param previous_stats Previous week field statistics
#'
#' @return Data frame with field-level CI changes
calculate_area_change_kpi <- function(current_stats, previous_stats) {
result <- calculate_change_percentages(current_stats, previous_stats)
# Add interpretation
result$interpretation <- NA_character_
for (i in seq_len(nrow(result))) {
change <- result$mean_ci_pct_change[i]
if (is.na(change)) {
result$interpretation[i] <- "No previous data"
} else if (change > 15) {
result$interpretation[i] <- "Rapid growth"
} else if (change > 5) {
result$interpretation[i] <- "Positive growth"
} else if (change > -5) {
result$interpretation[i] <- "Stable"
} else if (change > -15) {
result$interpretation[i] <- "Declining"
} else {
result$interpretation[i] <- "Rapid decline"
}
}
return(result)
}
#' KPI 3: Calculate TCH forecasted (tonnes of cane per hectare)
#'
#' Projects final harvest tonnage based on CI growth trajectory
#'
#' @param field_statistics Current field statistics
#' @param harvesting_data Historical harvest data (with yield observations)
#' @param field_boundaries_sf Field geometries
#'
#' @return Data frame with field-level TCH forecasts
calculate_tch_forecasted_kpi <- function(field_statistics, harvesting_data = NULL, field_boundaries_sf = NULL) {
result <- data.frame(
field_idx = field_statistics$field_idx,
mean_ci = field_statistics$mean_ci,
tch_forecasted = NA_real_,
tch_lower_bound = NA_real_,
tch_upper_bound = NA_real_,
confidence = NA_character_,
stringsAsFactors = FALSE
)
# Base TCH model: TCH = 50 + (CI * 10)
# This is a simplified model; production use should include more variables
for (i in seq_len(nrow(result))) {
if (is.na(result$mean_ci[i])) {
result$confidence[i] <- "No data"
next
}
ci_val <- result$mean_ci[i]
# Simple linear model
tch_est <- 50 + (ci_val * 10)
# Confidence interval based on CI range
tch_lower <- tch_est * 0.85
tch_upper <- tch_est * 1.15
result$tch_forecasted[i] <- round(tch_est, 2)
result$tch_lower_bound[i] <- round(tch_lower, 2)
result$tch_upper_bound[i] <- round(tch_upper, 2)
result$confidence[i] <- "Medium"
}
return(result)
}
#' KPI 4: Calculate growth decline indicator
#'
#' Identifies fields with negative growth trajectory
#'
#' @param ci_values_list List of CI values for each field (multiple weeks)
#'
#' @return Data frame with field-level decline indicators
calculate_growth_decline_kpi <- function(ci_values_list) {
result <- data.frame(
field_idx = seq_len(length(ci_values_list)),
four_week_trend = NA_real_,
trend_interpretation = NA_character_,
decline_severity = NA_character_,
stringsAsFactors = FALSE
)
for (field_idx in seq_len(length(ci_values_list))) {
ci_vals <- ci_values_list[[field_idx]]
if (is.null(ci_vals) || length(ci_vals) < 2) {
result$trend_interpretation[field_idx] <- "Insufficient data"
next
}
ci_vals <- ci_vals[!is.na(ci_vals)]
if (length(ci_vals) < 2) {
result$trend_interpretation[field_idx] <- "Insufficient data"
next
}
# Calculate linear trend
weeks <- seq_along(ci_vals)
lm_fit <- lm(ci_vals ~ weeks)
slope <- coef(lm_fit)["weeks"]
result$four_week_trend[field_idx] <- round(as.numeric(slope), 3)
if (slope > 0.1) {
result$trend_interpretation[field_idx] <- "Strong growth"
result$decline_severity[field_idx] <- "None"
} else if (slope > 0) {
result$trend_interpretation[field_idx] <- "Weak growth"
result$decline_severity[field_idx] <- "None"
} else if (slope > -0.1) {
result$trend_interpretation[field_idx] <- "Slight decline"
result$decline_severity[field_idx] <- "Low"
} else if (slope > -0.3) {
result$trend_interpretation[field_idx] <- "Moderate decline"
result$decline_severity[field_idx] <- "Medium"
} else {
result$trend_interpretation[field_idx] <- "Strong decline"
result$decline_severity[field_idx] <- "High"
}
}
return(result)
}
#' KPI 5: Calculate weed presence indicator
#'
#' Detects field fragmentation/patchiness (potential weed/pest pressure)
#'
#' @param ci_pixels_by_field List of CI pixel arrays for each field
#'
#' @return Data frame with fragmentation indicators
calculate_weed_presence_kpi <- function(ci_pixels_by_field) {
result <- data.frame(
field_idx = seq_len(length(ci_pixels_by_field)),
cv_value = NA_real_,
low_ci_percent = NA_real_,
fragmentation_index = NA_real_,
weed_pressure_risk = NA_character_,
stringsAsFactors = FALSE
)
for (field_idx in seq_len(length(ci_pixels_by_field))) {
ci_pixels <- ci_pixels_by_field[[field_idx]]
if (is.null(ci_pixels) || length(ci_pixels) == 0) {
result$weed_pressure_risk[field_idx] <- "No data"
next
}
ci_pixels <- ci_pixels[!is.na(ci_pixels)]
if (length(ci_pixels) == 0) {
result$weed_pressure_risk[field_idx] <- "No data"
next
}
cv_val <- calculate_cv(ci_pixels)
low_ci_pct <- sum(ci_pixels < 1.5) / length(ci_pixels) * 100
fragmentation <- cv_val * low_ci_pct / 100
result$cv_value[field_idx] <- cv_val
result$low_ci_percent[field_idx] <- round(low_ci_pct, 2)
result$fragmentation_index[field_idx] <- round(fragmentation, 3)
if (fragmentation > 0.15) {
result$weed_pressure_risk[field_idx] <- "High"
} else if (fragmentation > 0.08) {
result$weed_pressure_risk[field_idx] <- "Medium"
} else if (fragmentation > 0.04) {
result$weed_pressure_risk[field_idx] <- "Low"
} else {
result$weed_pressure_risk[field_idx] <- "Minimal"
}
}
return(result)
}
#' KPI 6: Calculate gap filling quality (data interpolation success)
#'
#' Measures how well cloud/missing data was interpolated during growth model
#'
#' @param ci_rds_path Path to combined CI RDS file (before/after interpolation)
#'
#' @return Data frame with gap-filling quality metrics
calculate_gap_filling_kpi <- function(ci_rds_path) {
if (!file.exists(ci_rds_path)) {
return(NULL)
}
tryCatch({
ci_data <- readRDS(ci_rds_path)
# ci_data should be a wide matrix: fields × weeks
# NA values = missing data before interpolation
# (Gap filling is done during growth model stage)
result <- data.frame(
field_idx = seq_len(nrow(ci_data)),
na_percent_pre_interpolation = NA_real_,
na_percent_post_interpolation = NA_real_,
gap_filling_success = NA_character_,
stringsAsFactors = FALSE
)
for (field_idx in seq_len(nrow(ci_data))) {
na_count <- sum(is.na(ci_data[field_idx, ]))
na_pct <- na_count / ncol(ci_data) * 100
if (na_pct == 0) {
result$gap_filling_success[field_idx] <- "No gaps (100% data)"
} else if (na_pct < 10) {
result$gap_filling_success[field_idx] <- "Excellent"
} else if (na_pct < 25) {
result$gap_filling_success[field_idx] <- "Good"
} else if (na_pct < 40) {
result$gap_filling_success[field_idx] <- "Fair"
} else {
result$gap_filling_success[field_idx] <- "Poor"
}
result$na_percent_pre_interpolation[field_idx] <- round(na_pct, 2)
}
return(result)
}, error = function(e) {
message(paste("Error calculating gap filling KPI:", e$message))
return(NULL)
})
}
# ============================================================================
# KPI ORCHESTRATOR AND REPORTING
# ============================================================================
#' Create summary tables for all 6 KPIs
#'
#' @param all_kpis List containing results from all 6 KPI functions
#'
#' @return List of summary data frames ready for reporting
create_summary_tables <- function(all_kpis) {
kpi_summary <- list(
uniformity = all_kpis$uniformity %>%
select(field_idx, cv_value, morans_i, uniformity_score, interpretation),
area_change = all_kpis$area_change %>%
select(field_idx, mean_ci_pct_change, interpretation),
tch_forecast = all_kpis$tch_forecasted %>%
select(field_idx, mean_ci, tch_forecasted, tch_lower_bound, tch_upper_bound, confidence),
growth_decline = all_kpis$growth_decline %>%
select(field_idx, four_week_trend, trend_interpretation, decline_severity),
weed_pressure = all_kpis$weed_presence %>%
select(field_idx, fragmentation_index, weed_pressure_risk),
gap_filling = all_kpis$gap_filling %>%
select(field_idx, na_percent_pre_interpolation, gap_filling_success)
)
return(kpi_summary)
}
#' Create detailed field-by-field KPI report
#'
#' @param field_df Data frame with field identifiers and acreage
#' @param all_kpis List with all KPI results
#' @param field_boundaries_sf SF object with field boundaries
#'
#' @return Data frame with one row per field, all KPI columns
create_field_detail_table <- function(field_df, all_kpis, field_boundaries_sf) {
result <- field_df %>%
left_join(
all_kpis$uniformity %>% select(field_idx, cv_value, uniformity_interpretation = interpretation),
by = c("field_idx")
) %>%
left_join(
all_kpis$area_change %>% select(field_idx, mean_ci_pct_change),
by = c("field_idx")
) %>%
left_join(
all_kpis$tch_forecasted %>% select(field_idx, tch_forecasted),
by = c("field_idx")
) %>%
left_join(
all_kpis$growth_decline %>% select(field_idx, decline_severity),
by = c("field_idx")
) %>%
left_join(
all_kpis$weed_presence %>% select(field_idx, weed_pressure_risk),
by = c("field_idx")
)
return(result)
}
#' Generate KPI text interpretation for inclusion in Word report
#'
#' @param all_kpis List with all KPI results
#'
#' @return Character string with formatted KPI summary text
create_field_kpi_text <- function(all_kpis) {
text_parts <- c(
"## AURA KPI ANALYSIS SUMMARY\n",
"### Field Uniformity\n",
paste(all_kpis$uniformity$interpretation, collapse = "; "), "\n",
"### Growth Trends\n",
paste(all_kpis$growth_decline$trend_interpretation, collapse = "; "), "\n",
"### Weed/Pest Pressure\n",
paste(all_kpis$weed_presence$weed_pressure_risk, collapse = "; "), "\n"
)
return(paste(text_parts, collapse = ""))
}
#' Export detailed KPI data to Excel/RDS
#'
#' @param all_kpis List with all KPI results
#' @param kpi_summary List with summary tables
#' @param output_dir Directory for output files
#' @param week Week number
#' @param year Year
#'
#' @return List of output file paths
export_kpi_data <- function(all_kpis, kpi_summary, output_dir, week, year) {
kpi_subdir <- file.path(output_dir, "kpis")
if (!dir.exists(kpi_subdir)) {
dir.create(kpi_subdir, recursive = TRUE)
}
# Export all KPI tables to a single Excel file
excel_file <- paste0(kpi_subdir, "/AURA_KPI_week_", sprintf("%02d_%d", week, year), ".xlsx")
sheets <- list(
"Uniformity" = as.data.frame(kpi_summary$uniformity),
"Area_Change" = as.data.frame(kpi_summary$area_change),
"TCH_Forecast" = as.data.frame(kpi_summary$tch_forecast),
"Growth_Decline" = as.data.frame(kpi_summary$growth_decline),
"Weed_Pressure" = as.data.frame(kpi_summary$weed_pressure),
"Gap_Filling" = as.data.frame(kpi_summary$gap_filling)
)
write_xlsx(sheets, excel_file)
message(paste("✓ AURA KPI data exported to:", excel_file))
# Also export to RDS for programmatic access
rds_file <- paste0(kpi_subdir, "/AURA_KPI_week_", sprintf("%02d_%d", week, year), ".rds")
saveRDS(all_kpis, rds_file)
message(paste("✓ AURA KPI RDS exported to:", rds_file))
return(list(excel = excel_file, rds = rds_file))
}
# ============================================================================
# ORCHESTRATOR FUNCTION
# ============================================================================
#' Calculate all 6 AURA KPIs
#'
#' Main entry point for AURA KPI calculation.
#' This function orchestrates the 6 KPI calculations and returns all results.
#'
#' @param field_boundaries_sf SF object with field geometries
#' @param current_week ISO week number (1-53)
#' @param current_year ISO week year
#' @param current_mosaic_dir Directory containing current week's mosaic
#' @param previous_mosaic_dir Directory containing previous week's mosaic (optional)
#' @param ci_rds_path Path to combined CI RDS file
#' @param harvesting_data Data frame with harvest data (optional)
#' @param output_dir Directory for KPI exports
#'
#' @return List with results from all 6 KPI functions
#'
#' @details
#' This function:
#' 1. Loads current week mosaic and extracts field statistics
#' 2. (Optionally) loads previous week mosaic for comparison metrics
#' 3. Calculates all 6 AURA KPIs
#' 4. Creates summary tables
#' 5. Exports results to Excel/RDS
#'
calculate_all_kpis <- function(
field_boundaries_sf,
current_week,
current_year,
current_mosaic_dir,
previous_mosaic_dir = NULL,
ci_rds_path = NULL,
harvesting_data = NULL,
output_dir = file.path(PROJECT_DIR, "output")
) {
message("\n============ AURA KPI CALCULATION (6 KPIs) ============")
# Load current week mosaic
message("Loading current week mosaic...")
current_mosaic <- load_weekly_ci_mosaic(current_mosaic_dir, current_week, current_year)
if (is.null(current_mosaic)) {
stop("Could not load current week mosaic")
}
# Extract field statistics
message("Extracting field statistics from current mosaic...")
current_stats <- extract_field_statistics_from_ci(current_mosaic, field_boundaries_sf)
ci_pixels_by_field <- extract_ci_values(current_mosaic, field_boundaries_sf)
# Load previous week mosaic (if available)
previous_stats <- NULL
if (!is.null(previous_mosaic_dir)) {
target_prev <- calculate_target_week_and_year(current_week, current_year, offset_weeks = 1)
message(paste("Loading previous week mosaic (week", target_prev$week, target_prev$year, ")..."))
previous_mosaic <- load_weekly_ci_mosaic(previous_mosaic_dir, target_prev$week, target_prev$year)
if (!is.null(previous_mosaic)) {
previous_stats <- extract_field_statistics_from_ci(previous_mosaic, field_boundaries_sf)
} else {
message("Previous week mosaic not available - skipping area change KPI")
}
}
# Calculate 6 KPIs
message("\nCalculating KPI 1: Field Uniformity...")
uniformity_kpi <- calculate_field_uniformity_kpi(ci_pixels_by_field, field_boundaries_sf, current_mosaic)
message("Calculating KPI 2: Area Change...")
if (!is.null(previous_stats)) {
area_change_kpi <- calculate_area_change_kpi(current_stats, previous_stats)
} else {
area_change_kpi <- data.frame(
field_idx = seq_len(nrow(field_boundaries_sf)),
mean_ci_pct_change = NA_real_,
interpretation = rep("No previous data", nrow(field_boundaries_sf))
)
}
message("Calculating KPI 3: TCH Forecasted...")
tch_kpi <- calculate_tch_forecasted_kpi(current_stats, harvesting_data, field_boundaries_sf)
message("Calculating KPI 4: Growth Decline...")
growth_decline_kpi <- calculate_growth_decline_kpi(
list(ci_pixels_by_field) # Would need historical data for real trend
)
message("Calculating KPI 5: Weed Presence...")
weed_kpi <- calculate_weed_presence_kpi(ci_pixels_by_field)
message("Calculating KPI 6: Gap Filling...")
gap_filling_kpi <- calculate_gap_filling_kpi(ci_rds_path)
# Compile results
all_kpis <- list(
uniformity = uniformity_kpi,
area_change = area_change_kpi,
tch_forecasted = tch_kpi,
growth_decline = growth_decline_kpi,
weed_presence = weed_kpi,
gap_filling = gap_filling_kpi
)
# Create summary tables
kpi_summary <- create_summary_tables(all_kpis)
# Export
export_paths <- export_kpi_data(all_kpis, kpi_summary, output_dir, current_week, current_year)
message(paste("\n✓ AURA KPI calculation complete. Week", current_week, current_year, "\n"))
return(all_kpis)
}

210
r_app/80_utils_cane_supply.R Normal file
View file

@ -0,0 +1,210 @@
# 80_UTILS_CANE_SUPPLY.R
# ============================================================================
# CANE SUPPLY CLIENT-SPECIFIC UTILITIES (SCRIPT 80 - CLIENT TYPE: cane_supply)
#
# Contains ANGATA and other cane supply-specific KPI and reporting functions.
#
# Currently, CANE_SUPPLY clients use the common utilities from 80_utils_common.R:
# - Weekly statistics (calculate_field_statistics, calculate_kpi_trends)
# - Field analysis reporting (generate_field_analysis_summary)
# - Excel export (export_field_analysis_excel)
#
# This file is structured to accommodate future ANGATA-specific functionality such as:
# - Custom yield models
# - Harvest readiness criteria
# - Supply chain integration hooks
# - ANGATA-specific alerting and messaging
#
# Orchestrator: (Placeholder - uses common functions)
# Dependencies: 00_common_utils.R, 80_utils_common.R
# Used by: 80_calculate_kpis.R (when client_type == "cane_supply")
# ============================================================================
library(terra)
library(sf)
library(dplyr)
library(tidyr)
library(readxl)
library(writexl)
# ============================================================================
# ANGATA-SPECIFIC HELPER FUNCTIONS (Placeholder Section)
# ============================================================================
#' Placeholder: ANGATA harvest readiness assessment
#'
#' Future implementation will integrate ANGATA-specific harvest readiness criteria:
#' - Maturation phase detection (CI threshold-based)
#' - Field age tracking (days since planting)
#' - Weather-based ripeness indicators (if available)
#' - Historical yield correlations
#'
#' @param field_ci CI values for the field
#' @param field_age_days Days since planting
#'
#' @return Character string with harvest readiness assessment
assess_harvest_readiness <- function(field_ci, field_age_days = NULL) {
# Placeholder implementation
# Real version would check:
# - Mean CI > 3.5 (maturation threshold)
# - Age > 350 days
# - Weekly growth rate < threshold (mature plateau)
if (is.null(field_ci) || all(is.na(field_ci))) {
return("No data available")
}
mean_ci <- mean(field_ci, na.rm = TRUE)
if (mean_ci > 3.5) {
return("Ready for harvest")
} else if (mean_ci > 2.5) {
return("Approaching harvest readiness")
} else {
return("Not ready - continue monitoring")
}
}
#' Placeholder: ANGATA supply chain status flags
#'
#' Future implementation will add supply chain-specific status indicators:
#' - Harvest scheduling readiness
#' - Equipment availability impact
#' - Transportation/logistics flags
#' - Quality parameter flags
#'
#' @param field_analysis Data frame with field analysis results
#'
#' @return Data frame with supply chain status columns
assess_supply_chain_status <- function(field_analysis) {
# Placeholder: return field analysis as-is
# Real version would add columns for:
# - schedule_ready (bool)
# - harvest_window_days (numeric)
# - transportation_priority (char)
# - quality_flags (char)
return(field_analysis)
}
# ============================================================================
# ORCHESTRATOR FOR CANE_SUPPLY WORKFLOWS
# ============================================================================
#' Orchestrate ANGATA weekly field analysis and reporting
#'
#' Main entry point for CANE_SUPPLY (ANGATA, etc.) workflows.
#' Currently uses common utilities; future versions will add client-specific logic.
#'
#' @param field_boundaries_sf SF object with field geometries
#' @param current_week ISO week number (1-53)
#' @param current_year ISO week year
#' @param mosaic_dir Directory containing weekly mosaics
#' @param field_boundaries_path Path to field GeoJSON
#' @param harvesting_data Data frame with harvest data (optional)
#' @param output_dir Directory for exports
#' @param data_dir Base data directory
#'
#' @return List with field analysis results
#'
#' @details
#' This function:
#' 1. Loads weekly mosaic and extracts field statistics
#' 2. Calculates field statistics (using common utilities)
#' 3. Prepares field analysis summary
#' 4. Exports to Excel/CSV/RDS
#' 5. (Future) Applies ANGATA-specific assessments
#'
calculate_field_analysis_cane_supply <- function(
field_boundaries_sf,
current_week,
current_year,
mosaic_dir,
field_boundaries_path = NULL,
harvesting_data = NULL,
output_dir = file.path(PROJECT_DIR, "output"),
data_dir = NULL
) {
message("\n============ CANE SUPPLY FIELD ANALYSIS (ANGATA, etc.) ============")
# Load current week mosaic
message("Loading current week mosaic...")
current_mosaic <- load_weekly_ci_mosaic(mosaic_dir, current_week, current_year)
if (is.null(current_mosaic)) {
warning(paste("Could not load current week mosaic for week", current_week, current_year))
return(NULL)
}
# Extract field statistics
message("Extracting field statistics from current mosaic...")
field_stats <- extract_field_statistics_from_ci(current_mosaic, field_boundaries_sf)
# Load previous week stats for comparison
message("Loading historical data for trends...")
target_prev <- calculate_target_week_and_year(current_week, current_year, offset_weeks = 1)
previous_stats <- NULL
previous_mosaic <- load_weekly_ci_mosaic(mosaic_dir, target_prev$week, target_prev$year)
if (!is.null(previous_mosaic)) {
previous_stats <- extract_field_statistics_from_ci(previous_mosaic, field_boundaries_sf)
}
# Calculate 4-week historical trend
message("Calculating field trends...")
ci_rds_path <- file.path(data_dir, "combined_CI", "combined_CI_data.rds")
field_analysis <- calculate_field_statistics(
field_stats = field_stats,
previous_stats = previous_stats,
week_num = current_week,
year = current_year,
ci_rds_path = ci_rds_path,
field_boundaries_sf = field_boundaries_sf,
harvesting_data = harvesting_data
)
if (is.null(field_analysis)) {
message("Could not generate field analysis")
return(NULL)
}
# Generate summary
message("Generating summary statistics...")
summary_df <- generate_field_analysis_summary(field_analysis)
# Export
message("Exporting field analysis...")
export_paths <- export_field_analysis_excel(
field_analysis,
summary_df,
PROJECT_DIR,
current_week,
current_year,
output_dir
)
message(paste("\n✓ CANE_SUPPLY field analysis complete. Week", current_week, current_year, "\n"))
result <- list(
field_analysis = field_analysis,
summary = summary_df,
exports = export_paths
)
return(result)
}
# ============================================================================
# FUTURE EXTENSION POINTS
# ============================================================================
# Placeholder for ANGATA-specific utilities that may be added in future:
# - Custom yield models based on ANGATA historical data
# - Field condition thresholds specific to ANGATA growing practices
# - Integration with ANGATA harvest scheduling system
# - WhatsApp messaging templates for ANGATA supply chain stakeholders
# - Cost/benefit analysis for ANGATA operational decisions
# These functions can be added here as ANGATA requirements evolve.

1469
r_app/80_weekly_stats_utils.R → r_app/80_utils_common.R
View file

File diff suppressed because it is too large Load diff

2
r_app/parameters_project.R
View file

@ -445,7 +445,7 @@ load_field_boundaries <- function(data_dir) {
if (use_pivot_2) {
field_boundaries_path <- here(data_dir, "pivot_2.geojson")
} else {
field_boundaries_path <- here(data_dir, "Data", "pivot.geojson")
field_boundaries_path <- here(data_dir, "pivot.geojson")
}
if (!file.exists(field_boundaries_path)) {