410 lines
11 KiB
R
410 lines
11 KiB
R
# CI (Canopy Index) Analysis Functions
|
|
# ====================================
|
|
#
|
|
# Functions for loading and analyzing CI data from Planet imagery
|
|
# Creates interactive visualizations for this week, last week, and change detection
|
|
#
|
|
# Author: Timon
|
|
# Date: August 2025
|
|
|
|
# Load required libraries
|
|
if (!require(terra)) install.packages("terra")
|
|
if (!require(dplyr)) install.packages("dplyr")
|
|
if (!require(sf)) install.packages("sf")
|
|
if (!require(ggplot2)) install.packages("ggplot2")
|
|
if (!require(tmap)) install.packages("tmap")
|
|
if (!require(RColorBrewer)) install.packages("RColorBrewer")
|
|
if (!require(viridis)) install.packages("viridis")
|
|
|
|
library(terra)
|
|
library(dplyr)
|
|
library(sf)
|
|
library(ggplot2)
|
|
library(tmap)
|
|
library(RColorBrewer)
|
|
library(viridis)
|
|
|
|
# CI Data Loading Functions
|
|
# =========================
|
|
|
|
#' Load CI weekly mosaic for a specific week
|
|
#' @param week_num Week number (27-34 available)
|
|
#' @param year Year (2025)
|
|
#' @param data_dir Directory containing weekly CI mosaics
|
|
#' @return SpatRaster object
|
|
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 weekly mosaic not found:", filename))
|
|
return(NULL)
|
|
}
|
|
|
|
# Load raster
|
|
ci_raster <- rast(filename)
|
|
|
|
# Set appropriate names
|
|
names(ci_raster) <- paste0("CI_week_", week_num, "_", year)
|
|
|
|
cat("✓ Loaded CI for week", week_num, ":", filename, "\n")
|
|
cat(" Dimensions:", dim(ci_raster), "\n")
|
|
cat(" CRS:", as.character(crs(ci_raster)), "\n")
|
|
|
|
# Get value range
|
|
ci_range <- global(ci_raster, range, na.rm = TRUE)
|
|
cat(" CI value range:", round(ci_range[1,1], 3), "to", round(ci_range[2,1], 3), "\n")
|
|
|
|
return(ci_raster)
|
|
}
|
|
|
|
#' Get the most recent available CI week
|
|
#' @param data_dir Directory containing weekly CI mosaics
|
|
#' @return Week number of most recent data
|
|
get_latest_ci_week <- function(data_dir = "../../../laravel_app/storage/app/aura/weekly_mosaic") {
|
|
|
|
# List all available weeks
|
|
files <- list.files(data_dir, pattern = "week_.*_2025\\.tif$", full.names = FALSE)
|
|
|
|
if (length(files) == 0) {
|
|
stop("No CI weekly mosaics found in ", data_dir)
|
|
}
|
|
|
|
# Extract week numbers
|
|
week_nums <- as.numeric(gsub("week_(\\d+)_2025\\.tif", "\\1", files))
|
|
latest_week <- max(week_nums)
|
|
|
|
cat("Available CI weeks:", paste(sort(week_nums), collapse = ", "), "\n")
|
|
cat("Latest week:", latest_week, "\n")
|
|
|
|
return(latest_week)
|
|
}
|
|
|
|
#' Load current and previous week CI data
|
|
#' @param data_dir Directory containing weekly CI mosaics
|
|
#' @return Named list with current_week, previous_week CI rasters and week numbers
|
|
load_ci_comparison_data <- function(data_dir = "../../../laravel_app/storage/app/aura/weekly_mosaic") {
|
|
|
|
# Get latest week
|
|
latest_week <- get_latest_ci_week(data_dir)
|
|
previous_week <- latest_week - 1
|
|
|
|
# Load both weeks
|
|
current_ci <- load_ci_weekly(latest_week, data_dir = data_dir)
|
|
previous_ci <- load_ci_weekly(previous_week, data_dir = data_dir)
|
|
|
|
if (is.null(current_ci)) {
|
|
stop("Could not load current week CI data")
|
|
}
|
|
|
|
if (is.null(previous_ci)) {
|
|
warning("Could not load previous week CI data, using available week")
|
|
# Try to find any available previous week
|
|
available_weeks <- get_available_ci_weeks(data_dir)
|
|
available_weeks <- available_weeks[available_weeks < latest_week]
|
|
if (length(available_weeks) > 0) {
|
|
previous_week <- max(available_weeks)
|
|
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
|
|
))
|
|
}
|
|
|
|
#' Get all available CI weeks
|
|
#' @param data_dir Directory containing weekly CI mosaics
|
|
#' @return Vector of available week numbers
|
|
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))
|
|
}
|
|
|
|
# CI Interactive Visualization Functions
|
|
# ======================================
|
|
|
|
#' Create interactive CI map
|
|
#' @param ci_raster CI SpatRaster
|
|
#' @param field_boundaries Field boundary polygons (optional)
|
|
#' @param title Map title
|
|
#' @param color_palette Color palette for CI values
|
|
#' @return tmap object
|
|
create_interactive_ci_map <- function(ci_raster, field_boundaries = NULL,
|
|
title = "CI Map",
|
|
color_palette = "viridis") {
|
|
|
|
# Set tmap to interactive mode
|
|
tmap_mode("view")
|
|
|
|
# Create base CI map
|
|
ci_map <- tm_shape(ci_raster) +
|
|
tm_raster(
|
|
title = "CI Value",
|
|
palette = color_palette,
|
|
style = "cont",
|
|
alpha = 0.8,
|
|
breaks = seq(0, 1, 0.1)
|
|
) +
|
|
tm_layout(
|
|
title = title,
|
|
title.size = 1.2
|
|
)
|
|
|
|
# Add field boundaries if provided
|
|
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
|
|
#' @param current_ci Current week CI raster
|
|
#' @param previous_ci Previous week CI raster
|
|
#' @param field_boundaries Field boundary polygons (optional)
|
|
#' @param title Map title
|
|
#' @return tmap object
|
|
create_ci_change_map <- function(current_ci, previous_ci, field_boundaries = NULL,
|
|
title = "CI Change Detection") {
|
|
|
|
# Calculate change (current - previous)
|
|
ci_change <- current_ci - previous_ci
|
|
|
|
# Set names for clarity
|
|
names(ci_change) <- "CI_Change"
|
|
|
|
# Create change map
|
|
change_map <- tm_shape(ci_change) +
|
|
tm_raster(
|
|
title = "CI Change",
|
|
palette = "RdBu",
|
|
style = "cont",
|
|
alpha = 0.8,
|
|
midpoint = 0,
|
|
breaks = seq(-0.5, 0.5, 0.1)
|
|
) +
|
|
tm_layout(
|
|
title = title,
|
|
title.size = 1.2
|
|
)
|
|
|
|
# Add field boundaries if provided
|
|
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)
|
|
}
|
|
|
|
# Get the second most recent file
|
|
file_info <- file.info(rds_files)
|
|
sorted_files <- rds_files[order(file_info$mtime, decreasing = TRUE)]
|
|
previous_file <- sorted_files[2]
|
|
|
|
cat("Loading previous CI from:", basename(previous_file), "\n")
|
|
|
|
ci_data <- readRDS(previous_file)
|
|
|
|
if (inherits(ci_data, "SpatRaster")) {
|
|
return(ci_data)
|
|
} else if (is.list(ci_data) && "ci_map" %in% names(ci_data)) {
|
|
return(ci_data$ci_map)
|
|
} else {
|
|
warning("Unexpected data format in previous RDS file")
|
|
return(NULL)
|
|
}
|
|
}
|
|
|
|
#' Calculate CI change between current and previous week
|
|
#' @param current_ci Current week CI raster
|
|
#' @param previous_ci Previous week CI raster
|
|
#' @return SpatRaster showing CI change
|
|
calculate_ci_change <- function(current_ci, previous_ci) {
|
|
|
|
if (is.null(current_ci) || is.null(previous_ci)) {
|
|
warning("Both current and previous CI data required for change calculation")
|
|
return(NULL)
|
|
}
|
|
|
|
cat("Calculating CI change...\n")
|
|
|
|
# Ensure both rasters have the same extent and resolution
|
|
if (!compareGeom(current_ci, previous_ci)) {
|
|
cat("Resampling rasters to match...\n")
|
|
previous_ci <- resample(previous_ci, current_ci)
|
|
}
|
|
|
|
# Calculate change
|
|
ci_change <- current_ci - previous_ci
|
|
names(ci_change) <- "CI_Change"
|
|
|
|
return(ci_change)
|
|
}
|
|
|
|
# CI Visualization Functions
|
|
# ==========================
|
|
|
|
#' Create interactive CI map
|
|
#' @param ci_raster CI raster data
|
|
#' @param field_boundaries Field boundary polygons
|
|
#' @param title Map title
|
|
#' @param palette Color palette to use
|
|
#' @return tmap object
|
|
create_interactive_ci_map <- function(ci_raster, field_boundaries = NULL,
|
|
title = "Canopy Index",
|
|
palette = "RdYlGn") {
|
|
|
|
if (is.null(ci_raster)) {
|
|
cat("No CI data available for mapping\n")
|
|
return(NULL)
|
|
}
|
|
|
|
# Set tmap to interactive mode
|
|
tmap_mode("view")
|
|
|
|
# Create base map
|
|
ci_map <- tm_shape(ci_raster) +
|
|
tm_raster(
|
|
title = "CI",
|
|
palette = palette,
|
|
style = "cont",
|
|
alpha = 0.8,
|
|
n = 10
|
|
) +
|
|
tm_layout(
|
|
title = title,
|
|
title.size = 1.2
|
|
)
|
|
|
|
# Add field boundaries if available
|
|
if (!is.null(field_boundaries)) {
|
|
ci_map <- ci_map +
|
|
tm_shape(field_boundaries) +
|
|
tm_polygons(
|
|
alpha = 0.2,
|
|
col = "white",
|
|
border.col = "black",
|
|
border.lwd = 1,
|
|
popup.vars = c("Field" = "field", "Sub Field" = "sub_field"),
|
|
id = "field"
|
|
)
|
|
}
|
|
|
|
return(ci_map)
|
|
}
|
|
|
|
#' Create CI change map
|
|
#' @param change_raster CI change raster
|
|
#' @param field_boundaries Field boundary polygons
|
|
#' @return tmap object
|
|
create_ci_change_map <- function(change_raster, field_boundaries = NULL) {
|
|
|
|
if (is.null(change_raster)) {
|
|
cat("No CI change data available for mapping\n")
|
|
return(NULL)
|
|
}
|
|
|
|
tmap_mode("view")
|
|
|
|
change_map <- tm_shape(change_raster) +
|
|
tm_raster(
|
|
title = "CI Change",
|
|
palette = "RdBu",
|
|
style = "cont",
|
|
alpha = 0.8,
|
|
midpoint = 0,
|
|
n = 11
|
|
) +
|
|
tm_layout(
|
|
title = "CI Change: This Week vs Last Week",
|
|
title.size = 1.2
|
|
)
|
|
|
|
# Add field boundaries
|
|
if (!is.null(field_boundaries)) {
|
|
change_map <- change_map +
|
|
tm_shape(field_boundaries) +
|
|
tm_polygons(
|
|
alpha = 0.2,
|
|
col = "white",
|
|
border.col = "yellow",
|
|
border.lwd = 2,
|
|
popup.vars = c("Field" = "field", "Sub Field" = "sub_field"),
|
|
id = "field"
|
|
)
|
|
}
|
|
|
|
return(change_map)
|
|
}
|
|
|
|
#' Load field boundaries
|
|
#' @param geojson_path Path to field boundaries geojson
|
|
#' @return sf object
|
|
load_field_boundaries <- function(geojson_path = "../../../pivot.geojson") {
|
|
|
|
if (!file.exists(geojson_path)) {
|
|
warning("Field boundaries file not found:", geojson_path)
|
|
return(NULL)
|
|
}
|
|
|
|
cat("Loading field boundaries from:", basename(geojson_path), "\n")
|
|
boundaries <- st_read(geojson_path, quiet = TRUE)
|
|
|
|
return(boundaries)
|
|
}
|
|
|
|
#' Calculate CI statistics by field
|
|
#' @param ci_raster CI raster data
|
|
#' @param field_boundaries Field boundary polygons
|
|
#' @return data.frame with field statistics
|
|
calculate_field_ci_stats <- function(ci_raster, field_boundaries) {
|
|
|
|
if (is.null(ci_raster) || is.null(field_boundaries)) {
|
|
return(NULL)
|
|
}
|
|
|
|
cat("Calculating CI statistics by field...\n")
|
|
|
|
# Extract CI values for each field
|
|
field_stats <- extract(ci_raster, field_boundaries, fun = c(mean, sd, min, max), na.rm = TRUE)
|
|
|
|
# Combine with field information
|
|
stats_df <- data.frame(
|
|
field = field_boundaries$field,
|
|
sub_field = field_boundaries$sub_field,
|
|
ci_mean = field_stats[,2], # mean
|
|
ci_sd = field_stats[,3], # sd
|
|
ci_min = field_stats[,4], # min
|
|
ci_max = field_stats[,5] # max
|
|
)
|
|
|
|
return(stats_df)
|
|
}
|
|
|
|
cat("CI analysis functions loaded successfully!\n")
|