script 40 per-field weekly mosaics - tested with aura
This commit is contained in:
parent
9b8c971902
commit
1c08070965
173
r_app/40_mosaic_creation_per_field.R
Normal file
173
r_app/40_mosaic_creation_per_field.R
Normal file
|
|
@ -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)
|
||||
})
|
||||
}
|
||||
275
r_app/40_mosaic_creation_per_field_utils.R
Normal file
275
r_app/40_mosaic_creation_per_field_utils.R
Normal file
|
|
@ -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)
|
||||
}
|
||||
Loading…
Reference in a new issue