156 lines
4.4 KiB
R
156 lines
4.4 KiB
R
# CI Analysis Functions - Clean Version
|
|
# =====================================
|
|
#
|
|
# Functions for interactive CI visualization of Aura fields
|
|
# Author: Timon, Date: August 2025
|
|
|
|
# Load required libraries
|
|
library(terra)
|
|
library(sf)
|
|
library(tmap)
|
|
library(dplyr)
|
|
library(ggplot2)
|
|
|
|
# Data Loading Functions
|
|
# ======================
|
|
|
|
#' Get available CI weeks
|
|
get_available_ci_weeks <- function(data_dir = "../../../laravel_app/storage/app/aura/weekly_mosaic") {
|
|
files <- list.files(data_dir, pattern = "week_.*_2025\\.tif$", full.names = FALSE)
|
|
week_nums <- as.numeric(gsub("week_(\\d+)_2025\\.tif", "\\1", files))
|
|
return(sort(week_nums))
|
|
}
|
|
|
|
#' Get latest CI week
|
|
get_latest_ci_week <- function(data_dir = "../../../laravel_app/storage/app/aura/weekly_mosaic") {
|
|
weeks <- get_available_ci_weeks(data_dir)
|
|
if (length(weeks) == 0) stop("No CI weekly mosaics found")
|
|
return(max(weeks))
|
|
}
|
|
|
|
#' Load CI weekly mosaic
|
|
load_ci_weekly <- function(week_num, year = 2025, data_dir = "../../../laravel_app/storage/app/aura/weekly_mosaic") {
|
|
filename <- file.path(data_dir, paste0("week_", week_num, "_", year, ".tif"))
|
|
|
|
if (!file.exists(filename)) {
|
|
warning(paste("CI file not found:", filename))
|
|
return(NULL)
|
|
}
|
|
|
|
ci_raster <- rast(filename)
|
|
|
|
cat("✓ Loaded CI week", week_num, "- layers:", nlyr(ci_raster), "\n")
|
|
cat(" Layer names:", paste(names(ci_raster), collapse = ", "), "\n")
|
|
|
|
# Extract only the CI band (usually the 5th band: R, G, B, NIR, CI)
|
|
if (nlyr(ci_raster) >= 5) {
|
|
ci_band <- ci_raster[[5]] # Select the 5th layer (CI)
|
|
names(ci_band) <- paste0("CI_week_", week_num, "_", year)
|
|
cat(" Extracted CI band (layer 5)\n")
|
|
return(ci_band)
|
|
} else if (nlyr(ci_raster) == 1) {
|
|
# If it's already a single band, assume it's CI
|
|
names(ci_raster) <- paste0("CI_week_", week_num, "_", year)
|
|
cat(" Using single band as CI\n")
|
|
return(ci_raster)
|
|
} else {
|
|
warning("Unexpected number of bands. Using first band as CI.")
|
|
ci_band <- ci_raster[[1]]
|
|
names(ci_band) <- paste0("CI_week_", week_num, "_", year)
|
|
return(ci_band)
|
|
}
|
|
}
|
|
|
|
#' Load current and previous week CI data
|
|
load_ci_comparison_data <- function(data_dir = "../../../laravel_app/storage/app/aura/weekly_mosaic") {
|
|
latest_week <- get_latest_ci_week(data_dir)
|
|
previous_week <- latest_week - 1
|
|
|
|
current_ci <- load_ci_weekly(latest_week, data_dir = data_dir)
|
|
previous_ci <- load_ci_weekly(previous_week, data_dir = data_dir)
|
|
|
|
return(list(
|
|
current_ci = current_ci,
|
|
previous_ci = previous_ci,
|
|
current_week = latest_week,
|
|
previous_week = previous_week
|
|
))
|
|
}
|
|
|
|
# Visualization Functions
|
|
# =======================
|
|
|
|
#' Create interactive CI map
|
|
create_interactive_ci_map <- function(ci_raster, field_boundaries = NULL, title = "CI Map", color_palette = "viridis") {
|
|
tmap_mode("view")
|
|
|
|
ci_map <- tm_shape(ci_raster) +
|
|
tm_raster(
|
|
title = "CI Value",
|
|
palette = color_palette,
|
|
style = "cont",
|
|
alpha = 0.8
|
|
) +
|
|
tm_layout(title = title, title.size = 1.2)
|
|
|
|
if (!is.null(field_boundaries)) {
|
|
ci_map <- ci_map +
|
|
tm_shape(field_boundaries) +
|
|
tm_polygons(
|
|
alpha = 0.1,
|
|
col = NA,
|
|
border.col = "white",
|
|
border.lwd = 2,
|
|
popup.vars = c("Field" = "field", "Sub Field" = "sub_field"),
|
|
id = "field"
|
|
)
|
|
}
|
|
|
|
return(ci_map)
|
|
}
|
|
|
|
#' Create CI change detection map
|
|
create_ci_change_map <- function(current_ci, previous_ci, field_boundaries = NULL, title = "CI Change") {
|
|
|
|
# Ensure both rasters are single-band
|
|
if (nlyr(current_ci) > 1) {
|
|
current_ci <- current_ci[[1]]
|
|
}
|
|
if (nlyr(previous_ci) > 1) {
|
|
previous_ci <- previous_ci[[1]]
|
|
}
|
|
|
|
# Calculate change
|
|
ci_change <- current_ci - previous_ci
|
|
|
|
# Set name for the change raster
|
|
names(ci_change) <- "CI_Change"
|
|
|
|
change_map <- tm_shape(ci_change) +
|
|
tm_raster(
|
|
title = "CI Change",
|
|
palette = "RdBu",
|
|
style = "cont",
|
|
alpha = 0.8,
|
|
midpoint = 0
|
|
) +
|
|
tm_layout(title = title, title.size = 1.2)
|
|
|
|
if (!is.null(field_boundaries)) {
|
|
change_map <- change_map +
|
|
tm_shape(field_boundaries) +
|
|
tm_polygons(
|
|
alpha = 0.1,
|
|
col = NA,
|
|
border.col = "yellow",
|
|
border.lwd = 2,
|
|
popup.vars = c("Field" = "field", "Sub Field" = "sub_field"),
|
|
id = "field"
|
|
)
|
|
}
|
|
|
|
return(change_map)
|
|
}
|
|
|
|
cat("✓ CI analysis functions loaded successfully\n")
|