From 1c080709655430952c9057e3cb839adb5deea6c5 Mon Sep 17 00:00:00 2001 From: Timon Date: Thu, 29 Jan 2026 21:49:15 +0100 Subject: [PATCH] script 40 per-field weekly mosaics - tested with aura --- r_app/40_mosaic_creation_per_field.R | 173 +++++++++++++ r_app/40_mosaic_creation_per_field_utils.R | 275 +++++++++++++++++++++ 2 files changed, 448 insertions(+) create mode 100644 r_app/40_mosaic_creation_per_field.R create mode 100644 r_app/40_mosaic_creation_per_field_utils.R diff --git a/r_app/40_mosaic_creation_per_field.R b/r_app/40_mosaic_creation_per_field.R new file mode 100644 index 0000000..e7bb27d --- /dev/null +++ b/r_app/40_mosaic_creation_per_field.R @@ -0,0 +1,173 @@ +# 40_MOSAIC_CREATION_PER_FIELD.R +# =============================== +# Per-Field Weekly Mosaic Creation +# +# Creates weekly mosaics FROM per-field daily CI TIFFs (output from Script 20) +# TO per-field weekly CI TIFFs (input for Scripts 90/91 reporting). +# +# ARCHITECTURE: +# Input: field_tiles_CI/{FIELD}/{DATE}.tif (5-band daily, per-field from Script 20) +# Output: weekly_mosaic/{FIELD}/week_WW_YYYY.tif (5-band weekly, per-field) +# +# USAGE: +# & "C:\Program Files\R\R-4.4.3\bin\x64\Rscript.exe" r_app/40_mosaic_creation_per_field.R [end_date] [offset] [project_dir] +# +# ARGUMENTS: +# end_date: End date for processing (YYYY-MM-DD format, default: today) +# offset: Days to look back from end_date (typically 7 for one week, default: 7) +# project_dir: Project directory (e.g., "aura", "angata", "chemba", "esa", default: "angata") +# +# EXAMPLES: +# & "C:\Program Files\R\R-4.4.3\bin\x64\Rscript.exe" r_app/40_mosaic_creation_per_field.R 2026-01-12 7 aura +# & "C:\Program Files\R\R-4.4.3\bin\x64\Rscript.exe" r_app/40_mosaic_creation_per_field.R 2025-12-31 7 angata + +# 1. Load required packages +# ----------------------- +suppressPackageStartupMessages({ + library(sf) + library(terra) + library(tidyverse) + library(lubridate) + library(here) +}) + +# 2. Main execution function +# ------------------------- +main <- function() { + + cat("\n") + cat("========================================================\n") + cat(" Script 40: Per-Field Weekly Mosaic Creation\n") + cat("========================================================\n\n") + + # Capture command line arguments + args <- commandArgs(trailingOnly = TRUE) + + # ==== Process Arguments ==== + + # Project directory + 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 { + project_dir <- "angata" + message("[INFO] No project_dir provided. Using default: angata") + } + + assign("project_dir", project_dir, envir = .GlobalEnv) + message(paste("[INFO] Project:", project_dir)) + + # End date + if (length(args) >= 1 && !is.na(args[1])) { + end_date <- as.Date(args[1], format = "%Y-%m-%d") + if (is.na(end_date)) { + message("[WARNING] Invalid end_date. Using current date.") + end_date <- Sys.Date() + } + } else { + end_date <- Sys.Date() + message(paste("[INFO] No end_date provided. Using current date:", format(end_date))) + } + + # Offset (days back) + if (length(args) >= 2 && !is.na(args[2])) { + offset <- as.numeric(args[2]) + if (is.na(offset) || offset <= 0) { + message("[WARNING] Invalid offset. Using default: 7 days") + offset <- 7 + } + } else { + offset <- 7 + message("[INFO] No offset provided. Using default: 7 days") + } + + # ==== Load Configuration ==== + + # Set working directory if needed + if (basename(getwd()) == "r_app") { + setwd("..") + } + + tryCatch({ + source("r_app/parameters_project.R") + message("[INFO] ✓ Loaded parameters_project.R") + }, error = function(e) { + stop("[ERROR] Failed to load parameters_project.R: ", e$message) + }) + + tryCatch({ + source("r_app/40_mosaic_creation_per_field_utils.R") + message("[INFO] ✓ Loaded 40_mosaic_creation_per_field_utils.R") + }, error = function(e) { + stop("[ERROR] Failed to load utilities: ", e$message) + }) + + # ==== Get Project Directories ==== + + setup <- setup_project_directories(project_dir) + + # Determine input/output directories + # Input: field_tiles_CI/ (from Script 20) + field_tiles_ci_dir <- setup$field_tiles_ci_dir + + # Output: weekly_mosaic/ (for Scripts 90/91) + weekly_mosaic_output_dir <- file.path(setup$laravel_storage_dir, "weekly_mosaic") + + message(paste("[INFO] Input directory:", field_tiles_ci_dir)) + message(paste("[INFO] Output directory:", weekly_mosaic_output_dir)) + + # ==== Validate Input Directory ==== + + if (!dir.exists(field_tiles_ci_dir)) { + stop(paste("[ERROR] Input directory not found:", field_tiles_ci_dir, + "\nScript 20 (CI extraction) must be run first to create per-field TIFFs.")) + } + + # Check if directory has any TIFFs + field_dirs <- list.dirs(field_tiles_ci_dir, full.names = FALSE, recursive = FALSE) + if (length(field_dirs) == 0) { + stop(paste("[ERROR] No field subdirectories found in:", field_tiles_ci_dir)) + } + + message(paste("[INFO] Found", length(field_dirs), "field directories")) + + # ==== Generate Date Range ==== + + dates <- date_list(end_date, offset) + + # ==== Create Per-Field Weekly Mosaics ==== + + created_files <- create_all_field_weekly_mosaics( + dates = dates, + field_tiles_ci_dir = field_tiles_ci_dir, + output_dir = weekly_mosaic_output_dir + ) + + # ==== Summary ==== + + message("\n") + message("========================================================") + message(paste(" COMPLETED")) + message(paste(" Created:", length(created_files), "weekly field mosaics")) + message("========================================================\n") + + if (length(created_files) > 0) { + message("[SUCCESS] Weekly mosaics ready for reporting (Scripts 90/91)") + } else { + message("[WARNING] No mosaics created - check input data") + } + + return(invisible(created_files)) +} + +# Execute main if script is run directly +if (sys.nframe() == 0) { + tryCatch({ + created <- main() + quit(save = "no", status = 0) + }, error = function(e) { + message(paste("\n[FATAL ERROR]", e$message)) + quit(save = "no", status = 1) + }) +} diff --git a/r_app/40_mosaic_creation_per_field_utils.R b/r_app/40_mosaic_creation_per_field_utils.R new file mode 100644 index 0000000..bf49773 --- /dev/null +++ b/r_app/40_mosaic_creation_per_field_utils.R @@ -0,0 +1,275 @@ +# MOSAIC_CREATION_PER_FIELD_UTILS.R +# ================================== +# Utility functions for creating per-field weekly mosaics from per-field daily TIFFs. +# +# This module aggregates daily per-field 5-band TIFFs (R,G,B,NIR,CI from Script 20) +# into weekly per-field mosaics using MAX compositing. +# +# DATA FLOW: +# Script 20: field_tiles_CI/{FIELD}/{DATE}.tif (5-band daily, per-field) +# ↓ +# Script 40 NEW (this module): +# For each field: +# For each week: +# - Find all daily TIFFs for that week +# - Stack & create MAX composite +# - Save: weekly_mosaic/{FIELD}/week_WW_YYYY.tif +# ↓ +# Scripts 90/91: Read weekly_mosaic/{FIELD}/week_WW_YYYY.tif (unchanged interface) + +#' Safe logging function +#' @param message The message to log +#' @param level The log level (default: "INFO") +#' @return NULL (used for side effects) +#' +safe_log <- function(message, level = "INFO") { + if (exists("log_message")) { + log_message(message, level) + } else { + if (level %in% c("ERROR", "WARNING")) { + warning(message) + } else { + message(message) + } + } +} + +#' Generate date range for processing (ISO week-based) +#' +#' @param end_date The end date (Date object or YYYY-MM-DD string) +#' @param offset Number of days to look back from end_date (typically 7 for one week) +#' @return List with week, year, start_date, end_date, days_filter (vector of YYYY-MM-DD strings) +#' +date_list <- function(end_date, offset) { + if (!lubridate::is.Date(end_date)) { + end_date <- as.Date(end_date) + if (is.na(end_date)) { + stop("Invalid end_date. Expected Date object or YYYY-MM-DD string.") + } + } + + offset <- as.numeric(offset) + if (is.na(offset) || offset < 1) { + stop("Invalid offset. Expected positive number.") + } + + offset <- offset - 1 # Adjust to include end_date + start_date <- end_date - lubridate::days(offset) + + week <- lubridate::isoweek(end_date) + year <- lubridate::isoyear(end_date) + + days_filter <- seq(from = start_date, to = end_date, by = "day") + days_filter <- format(days_filter, "%Y-%m-%d") + + safe_log(paste("Date range:", start_date, "to", end_date, + "(week", week, "of", year, ")")) + + return(list( + week = week, + year = year, + start_date = start_date, + end_date = end_date, + days_filter = days_filter + )) +} + +#' Find all per-field daily TIFFs for a specific week +#' +#' @param field_tiles_ci_dir Base directory containing per-field TIFFs +#' (e.g., field_tiles_CI/) +#' @param days_filter Vector of YYYY-MM-DD dates to match +#' @return List with field names and their matching TIFF files for the week +#' +find_per_field_tiffs_for_week <- function(field_tiles_ci_dir, days_filter) { + + if (!dir.exists(field_tiles_ci_dir)) { + safe_log(paste("Field TIFFs directory not found:", field_tiles_ci_dir), "WARNING") + return(list()) + } + + # List all field subdirectories + field_dirs <- list.dirs(field_tiles_ci_dir, full.names = FALSE, recursive = FALSE) + + if (length(field_dirs) == 0) { + safe_log("No field subdirectories found in field_tiles_CI/", "WARNING") + return(list()) + } + + # For each field, find TIFF files matching the week's dates + field_tiffs <- list() + + for (field in field_dirs) { + field_path <- file.path(field_tiles_ci_dir, field) + + # Find all TIFF files in this field directory + tiff_files <- list.files(field_path, pattern = "\\.tif$", full.names = TRUE) + + if (length(tiff_files) == 0) { + next # Skip fields with no TIFFs + } + + # Filter to only those matching week's dates + matching_files <- tiff_files[grepl(paste(days_filter, collapse = "|"), tiff_files)] + + if (length(matching_files) > 0) { + field_tiffs[[field]] <- sort(matching_files) + } + } + + safe_log(paste("Found TIFFs for", length(field_tiffs), "fields in week")) + + return(field_tiffs) +} + +#' Create weekly MAX composite for a single field +#' +#' Loads all daily TIFFs for a field+week combination and creates a MAX composite +#' (per-band maximum across all days). +#' +#' @param tiff_files Vector of TIFF file paths for this field+week +#' @param field_name Name of the field (for logging) +#' @return SpatRaster with 5 bands (R,G,B,NIR,CI), or NULL if fails +#' +create_field_weekly_composite <- function(tiff_files, field_name) { + + if (length(tiff_files) == 0) { + safe_log(paste("No TIFF files for field:", field_name), "WARNING") + return(NULL) + } + + tryCatch({ + # Load all TIFFs + rasters <- list() + for (file in tiff_files) { + tryCatch({ + r <- terra::rast(file) + rasters[[length(rasters) + 1]] <- r + }, error = function(e) { + safe_log(paste("Warning: Could not load", basename(file), "for field", field_name), "WARNING") + }) + } + + if (length(rasters) == 0) { + safe_log(paste("Failed to load any rasters for field:", field_name), "ERROR") + return(NULL) + } + + # Create MAX composite + if (length(rasters) == 1) { + composite <- rasters[[1]] + safe_log(paste(" Field", field_name, "- single day (no compositing needed)")) + } else { + # Stack all rasters and apply MAX per pixel per band + collection <- terra::sprc(rasters) + composite <- terra::mosaic(collection, fun = "max") + safe_log(paste(" Field", field_name, "- MAX composite from", length(rasters), "days")) + } + + # Ensure 5 bands with expected names + if (terra::nlyr(composite) >= 5) { + composite <- terra::subset(composite, 1:5) + names(composite) <- c("Red", "Green", "Blue", "NIR", "CI") + } else { + safe_log(paste("Warning: Field", field_name, "has", terra::nlyr(composite), + "bands (expected 5)"), "WARNING") + } + + return(composite) + + }, error = function(e) { + safe_log(paste("Error creating composite for field", field_name, ":", e$message), "ERROR") + return(NULL) + }) +} + +#' Save per-field weekly mosaic +#' +#' @param raster SpatRaster to save +#' @param output_dir Base output directory (e.g., laravel_app/storage/app/{project}/weekly_mosaic/) +#' @param field_name Name of the field +#' @param week Week number (ISO week) +#' @param year Year (ISO year) +#' @return File path of saved TIFF, or NULL if fails +#' +save_field_weekly_mosaic <- function(raster, output_dir, field_name, week, year) { + + if (is.null(raster)) { + return(NULL) + } + + tryCatch({ + # Create field-specific output directory + field_output_dir <- file.path(output_dir, field_name) + dir.create(field_output_dir, recursive = TRUE, showWarnings = FALSE) + + # Generate filename: week_WW_YYYY.tif + filename <- sprintf("week_%02d_%04d.tif", week, year) + file_path <- file.path(field_output_dir, filename) + + # Save raster + terra::writeRaster(raster, file_path, overwrite = TRUE) + + safe_log(paste(" Saved:", basename(field_output_dir), "/", filename)) + + return(file_path) + + }, error = function(e) { + safe_log(paste("Error saving mosaic for field", field_name, ":", e$message), "ERROR") + return(NULL) + }) +} + +#' Create all weekly mosaics for all fields in a week +#' +#' Main orchestration function. Loops over all fields and creates weekly mosaics. +#' +#' @param dates List from date_list() - contains week, year, days_filter +#' @param field_tiles_ci_dir Input: field_tiles_CI/ directory +#' @param output_dir Output: weekly_mosaic/ directory +#' @return Vector of successfully created file paths +#' +create_all_field_weekly_mosaics <- function(dates, field_tiles_ci_dir, output_dir) { + + safe_log(paste("=== Creating Per-Field Weekly Mosaics ===")) + safe_log(paste("Week:", dates$week, "Year:", dates$year)) + + # Find all per-field TIFFs for this week + field_tiffs <- find_per_field_tiffs_for_week(field_tiles_ci_dir, dates$days_filter) + + if (length(field_tiffs) == 0) { + safe_log("No per-field TIFFs found for this week - returning empty list", "WARNING") + return(character()) + } + + safe_log(paste("Processing", length(field_tiffs), "fields...")) + + created_files <- character() + + # Process each field + for (field_name in names(field_tiffs)) { + tiff_files <- field_tiffs[[field_name]] + + # Create composite + composite <- create_field_weekly_composite(tiff_files, field_name) + + if (!is.null(composite)) { + # Save + saved_path <- save_field_weekly_mosaic( + composite, + output_dir, + field_name, + dates$week, + dates$year + ) + + if (!is.null(saved_path)) { + created_files <- c(created_files, saved_path) + } + } + } + + safe_log(paste("✓ Completed: Created", length(created_files), "weekly field mosaics")) + + return(created_files) +}