Timon
07aee7bed1
- Replace raster package with terra throughout the codebase - Update map visualizations with better layout and legends - Add descriptive headers to report sections - Improve map legend positioning and sizing - Enhance error handling for missing data - Remove redundant legends in field-specific visualizations - Optimize figure dimensions to prevent page overflow - Expand documentation of CI index and report components - Update package dependencies in packages.
310 lines
10 KiB
R
310 lines
10 KiB
R
# 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.
|
|
|
|
#' Safe logging function that works whether log_message exists or not
|
|
#'
|
|
#' @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 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
|
|
#'
|
|
date_list <- function(end_date, offset) {
|
|
# Input validation
|
|
if (!lubridate::is.Date(end_date)) {
|
|
end_date <- as.Date(end_date)
|
|
if (is.na(end_date)) {
|
|
stop("Invalid end_date provided. Expected a Date object or a string convertible to Date.")
|
|
}
|
|
}
|
|
|
|
offset <- as.numeric(offset)
|
|
if (is.na(offset) || offset < 1) {
|
|
stop("Invalid offset provided. Expected a positive number.")
|
|
}
|
|
|
|
# Calculate date range
|
|
offset <- offset - 1 # Adjust offset to include end_date
|
|
start_date <- end_date - lubridate::days(offset)
|
|
|
|
# Extract week and year information
|
|
week <- lubridate::week(start_date)
|
|
year <- lubridate::year(start_date)
|
|
|
|
# Generate sequence of dates
|
|
days_filter <- seq(from = start_date, to = end_date, by = "day")
|
|
days_filter <- format(days_filter, "%Y-%m-%d") # Format for consistent filtering
|
|
|
|
# Log the date range
|
|
safe_log(paste("Date range generated from", start_date, "to", end_date))
|
|
|
|
return(list(
|
|
"week" = week,
|
|
"year" = year,
|
|
"days_filter" = days_filter,
|
|
"start_date" = start_date,
|
|
"end_date" = end_date
|
|
))
|
|
}
|
|
|
|
#' 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 = TRUE) {
|
|
# Find VRT files for the specified date range
|
|
vrt_list <- find_vrt_files(daily_vrt_dir, dates)
|
|
|
|
# Find final raster files for fallback
|
|
raster_files_final <- list.files(merged_final_dir, full.names = TRUE, pattern = "\\.tif$")
|
|
|
|
# Process the mosaic if VRT files are available
|
|
if (length(vrt_list) > 0) {
|
|
safe_log("VRT list created, assessing cloud cover for mosaic creation")
|
|
|
|
# Calculate cloud cover statistics
|
|
missing_pixels_count <- count_cloud_coverage(vrt_list, field_boundaries)
|
|
|
|
# Create mosaic based on cloud cover assessment
|
|
mosaic <- create_mosaic(vrt_list, missing_pixels_count, field_boundaries, raster_files_final)
|
|
|
|
} else {
|
|
safe_log("No VRT files available for the date range, creating empty mosaic", "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]) %>%
|
|
terra::setValues(0) %>%
|
|
terra::crop(field_boundaries, mask = TRUE)
|
|
|
|
names(mosaic) <- c("Red", "Green", "Blue", "NIR", "CI")
|
|
}
|
|
|
|
# Save the mosaic
|
|
file_path <- save_mosaic(mosaic, output_dir, file_name_tif, create_plots)
|
|
|
|
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
|
|
vrt_files <- list.files(here::here(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
|
|
#'
|
|
#' @param vrt_list List of VRT files to analyze
|
|
#' @param field_boundaries Field boundaries vector for masking
|
|
#' @return Data frame with cloud coverage statistics
|
|
#'
|
|
count_cloud_coverage <- function(vrt_list, field_boundaries) {
|
|
if (length(vrt_list) == 0) {
|
|
warning("No VRT files provided for cloud coverage calculation")
|
|
return(NULL)
|
|
}
|
|
|
|
tryCatch({
|
|
# Calculate total pixel area using the first VRT file
|
|
total_pix_area <- terra::rast(vrt_list[1]) %>%
|
|
terra::subset(1) %>%
|
|
terra::setValues(1) %>%
|
|
terra::crop(field_boundaries, mask = TRUE) %>%
|
|
terra::global(fun = "notNA")
|
|
|
|
# Extract layer 1 from all rasters (for cloud detection)
|
|
layer_5_list <- purrr::map(vrt_list, function(file) {
|
|
terra::rast(file) %>% terra::subset(1)
|
|
}) %>% terra::rast()
|
|
|
|
# Calculate percentage of missing pixels (clouds)
|
|
missing_pixels_count <- terra::global(layer_5_list, fun = "notNA") %>%
|
|
dplyr::mutate(
|
|
total_pixels = total_pix_area$notNA,
|
|
missing_pixels_percentage = round(100 - ((notNA / total_pix_area$notNA) * 100)),
|
|
thres_5perc = as.integer(missing_pixels_percentage < 5),
|
|
thres_40perc = as.integer(missing_pixels_percentage < 45)
|
|
)
|
|
|
|
# Log results
|
|
safe_log(paste(
|
|
"Cloud cover assessment completed for", length(vrt_list), "files.",
|
|
sum(missing_pixels_count$thres_5perc), "files with <5% cloud cover,",
|
|
sum(missing_pixels_count$thres_40perc), "files with <45% cloud cover"
|
|
))
|
|
|
|
return(missing_pixels_count)
|
|
}, error = function(e) {
|
|
warning("Error in cloud coverage calculation: ", e$message)
|
|
return(NULL)
|
|
})
|
|
}
|
|
|
|
#' Create a mosaic from VRT files based on cloud coverage
|
|
#'
|
|
#' @param vrt_list List of VRT files to create mosaic from
|
|
#' @param missing_pixels_count Cloud coverage statistics from count_cloud_coverage()
|
|
#' @param field_boundaries Field boundaries vector for masking (optional)
|
|
#' @param raster_files_final List of processed raster files to use as fallback
|
|
#' @return A SpatRaster object with the mosaic
|
|
#'
|
|
create_mosaic <- function(vrt_list, missing_pixels_count, field_boundaries = NULL, raster_files_final = NULL) {
|
|
# If no VRT files, create an empty mosaic
|
|
if (length(vrt_list) == 0) {
|
|
if (length(raster_files_final) == 0 || is.null(field_boundaries)) {
|
|
stop("No VRT files available and no fallback raster files or field boundaries provided")
|
|
}
|
|
|
|
safe_log("No images available for this period, creating empty mosaic", "WARNING")
|
|
|
|
x <- terra::rast(raster_files_final[1]) %>%
|
|
terra::setValues(0) %>%
|
|
terra::crop(field_boundaries, mask = TRUE)
|
|
|
|
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
|
return(x)
|
|
}
|
|
|
|
# If missing pixel count was not calculated, use all files
|
|
if (is.null(missing_pixels_count)) {
|
|
safe_log("No cloud coverage data available, using all images", "WARNING")
|
|
|
|
rsrc <- terra::sprc(vrt_list)
|
|
x <- terra::mosaic(rsrc, fun = "max")
|
|
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
|
return(x)
|
|
}
|
|
|
|
# Determine best rasters to use based on cloud coverage
|
|
index_5perc <- which(missing_pixels_count$thres_5perc == max(missing_pixels_count$thres_5perc))
|
|
index_40perc <- which(missing_pixels_count$thres_40perc == max(missing_pixels_count$thres_40perc))
|
|
|
|
# Create mosaic based on available cloud-free images
|
|
if (sum(missing_pixels_count$thres_5perc) > 1) {
|
|
safe_log("Creating max composite from multiple cloud-free images (<5% clouds)")
|
|
|
|
cloudy_rasters_list <- vrt_list[index_5perc]
|
|
rsrc <- terra::sprc(cloudy_rasters_list)
|
|
x <- terra::mosaic(rsrc, fun = "max")
|
|
|
|
} else if (sum(missing_pixels_count$thres_5perc) == 1) {
|
|
safe_log("Using single cloud-free image (<5% clouds)")
|
|
|
|
x <- terra::rast(vrt_list[index_5perc[1]])
|
|
|
|
} else if (sum(missing_pixels_count$thres_40perc) > 1) {
|
|
safe_log("Creating max composite from partially cloudy images (<40% clouds)", "WARNING")
|
|
|
|
cloudy_rasters_list <- vrt_list[index_40perc]
|
|
rsrc <- terra::sprc(cloudy_rasters_list)
|
|
x <- terra::mosaic(rsrc, fun = "max")
|
|
|
|
} else if (sum(missing_pixels_count$thres_40perc) == 1) {
|
|
safe_log("Using single partially cloudy image (<40% clouds)", "WARNING")
|
|
|
|
x <- terra::rast(vrt_list[index_40perc[1]])
|
|
|
|
} else {
|
|
safe_log("No cloud-free images available, using all images", "WARNING")
|
|
|
|
rsrc <- terra::sprc(vrt_list)
|
|
x <- terra::mosaic(rsrc, fun = "max")
|
|
}
|
|
|
|
# Set consistent layer names
|
|
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
|
return(x)
|
|
}
|
|
|
|
#' 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)
|
|
#' @return The file path of the saved raster
|
|
#'
|
|
save_mosaic <- function(mosaic_raster, output_dir, file_name, plot_result = 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
|
|
file_path <- here::here(output_dir, file_name)
|
|
|
|
# Save raster
|
|
terra::writeRaster(mosaic_raster, file_path, overwrite = TRUE)
|
|
|
|
# Create plots if requested
|
|
if (plot_result) {
|
|
if ("CI" %in% names(mosaic_raster)) {
|
|
terra::plot(mosaic_raster$CI, main = paste("CI map", file_name))
|
|
}
|
|
|
|
if (all(c("Red", "Green", "Blue") %in% names(mosaic_raster))) {
|
|
terra::plotRGB(mosaic_raster, main = paste("RGB map", file_name))
|
|
}
|
|
}
|
|
|
|
# Log save completion
|
|
safe_log(paste("Mosaic saved to:", file_path))
|
|
|
|
return(file_path)
|
|
}
|