708 lines
26 KiB
R
708 lines
26 KiB
R
#
|
|
# YOUNG_FIELD_ANALYSIS.R
|
|
# ======================
|
|
# Specialized analysis for young sugarcane fields (0-12 months)
|
|
# Focuses on:
|
|
# 1. Germination issues (gap detection, uneven emergence)
|
|
# 2. Weed pressure detection (both patchy and uniform weeds)
|
|
# 3. Age-specific intervention recommendations
|
|
#
|
|
# Usage: Rscript young_field_analysis.R [current_week] [previous_week] [project_dir]
|
|
#
|
|
|
|
# 1. Load required packages
|
|
# -----------------------
|
|
suppressPackageStartupMessages({
|
|
library(sf)
|
|
library(terra)
|
|
library(tidyverse)
|
|
library(lubridate)
|
|
library(here)
|
|
library(readxl) # For reading harvest_data.xlsx
|
|
library(spdep) # For spatial statistics
|
|
})
|
|
|
|
# 2. Young field analysis configuration
|
|
# -----------------------------------
|
|
# Age thresholds (months)
|
|
GERMINATION_PHASE_MAX <- 4 # 0-4 months: critical germination period
|
|
WEED_CRITICAL_PHASE_MAX <- 8 # 4-8 months: weeds most competitive
|
|
YOUNG_FIELD_MAX <- 12 # 0-12 months: all young field analysis
|
|
|
|
# Detection thresholds
|
|
WEED_CI_CHANGE_THRESHOLD <- 1.5 # Weekly CI increase indicating possible weeds
|
|
SEVERE_WEED_CI_CHANGE <- 2.0 # Severe weed pressure
|
|
GERMINATION_GAP_CV_THRESHOLD <- 0.30 # High CV in young fields = poor germination
|
|
LOW_CI_GERMINATION_THRESHOLD <- 0.5 # Very low CI = poor emergence
|
|
|
|
# 3. Enhanced spectral indices for young crop detection
|
|
# ---------------------------------------------------
|
|
|
|
#' Calculate enhanced vegetation indices from RGBNIR bands
|
|
#' @param red Red band raster
|
|
#' @param green Green band raster
|
|
#' @param blue Blue band raster
|
|
#' @param nir Near-infrared band raster
|
|
#' @return List of vegetation index rasters
|
|
calculate_enhanced_indices <- function(red, green, blue, nir) {
|
|
|
|
cat("Calculating enhanced vegetation indices...\n")
|
|
|
|
# 1. Standard NDVI (baseline)
|
|
ndvi <- (nir - red) / (nir + red)
|
|
names(ndvi) <- "NDVI"
|
|
|
|
# 2. Green NDVI (more sensitive to early vegetation)
|
|
gndvi <- (nir - green) / (nir + green)
|
|
names(gndvi) <- "GNDVI"
|
|
|
|
# 3. Excess Green Index (early vegetation detection)
|
|
exg <- 2*green - red - blue
|
|
names(exg) <- "ExG"
|
|
|
|
# 4. Visible Atmospherically Resistant Index
|
|
vari <- (green - red) / (green + red + blue)
|
|
names(vari) <- "VARI"
|
|
|
|
# 5. Green Ratio Vegetation Index
|
|
grvi <- green / red
|
|
names(grvi) <- "GRVI"
|
|
|
|
# 6. Chlorophyll Index (CI = NIR / Green - 1, NOT NIR/Red)
|
|
# *** CRITICAL: Correct formula uses GREEN band, not RED ***
|
|
ci <- nir / green - 1
|
|
names(ci) <- "CI"
|
|
|
|
return(list(
|
|
NDVI = ndvi,
|
|
GNDVI = gndvi,
|
|
ExG = exg,
|
|
VARI = vari,
|
|
GRVI = grvi,
|
|
CI = ci
|
|
))
|
|
}
|
|
|
|
#' Create visual maps of different indices for comparison
|
|
#' @param indices_list List of index rasters
|
|
#' @param field_boundaries SF object with field boundaries
|
|
#' @param output_dir Directory to save maps
|
|
#' @param week_num Week number for filename
|
|
create_index_comparison_maps <- function(indices_list, field_boundaries, output_dir, week_num) {
|
|
|
|
if (!dir.exists(output_dir)) {
|
|
dir.create(output_dir, recursive = TRUE)
|
|
}
|
|
|
|
cat("Creating visual comparison maps and TIF exports...\n")
|
|
|
|
# Create individual maps for each index
|
|
for (index_name in names(indices_list)) {
|
|
|
|
tryCatch({
|
|
index_raster <- indices_list[[index_name]]
|
|
|
|
# Create filenames
|
|
map_filename <- file.path(output_dir, paste0("week_", sprintf("%02d", week_num), "_", index_name, "_map.png"))
|
|
tif_filename <- file.path(output_dir, paste0("week_", sprintf("%02d", week_num), "_", index_name, ".tif"))
|
|
|
|
# Export TIF file for QGIS
|
|
terra::writeRaster(index_raster, tif_filename, overwrite = TRUE)
|
|
cat("- Exported TIF:", tif_filename, "\n")
|
|
|
|
# Create PNG map without field boundaries
|
|
png(map_filename, width = 1200, height = 800, res = 150)
|
|
|
|
# Plot raster only (no field boundaries)
|
|
plot(index_raster, main = paste("Week", week_num, "-", index_name))
|
|
|
|
dev.off()
|
|
|
|
cat("- Saved map:", map_filename, "\n")
|
|
|
|
}, error = function(e) {
|
|
warning(paste("Error creating outputs for", index_name, ":", e$message))
|
|
})
|
|
}
|
|
|
|
cat("✓ Index comparison maps and TIF files created\n")
|
|
}
|
|
|
|
#' Load field age data from existing harvest_data object (loaded in parameters_project.R)
|
|
#' @param harvest_data Data frame with harvest/field information
|
|
#' @param field_boundaries SF object with field boundaries
|
|
#' @return Field boundaries with age information added
|
|
#' # harvest_data = harvesting_data
|
|
#' # field_boundaries = field_boundaries_sf
|
|
load_field_ages_from_existing <- function(harvest_data, field_boundaries) {
|
|
|
|
cat("Processing field age data from existing harvest_data object...\n")
|
|
|
|
tryCatch({
|
|
# Display structure of harvest data
|
|
cat("- Harvest data columns:", paste(colnames(harvest_data), collapse = ", "), "\n")
|
|
cat("- Harvest data rows:", nrow(harvest_data), "\n")
|
|
|
|
# Display field boundaries structure
|
|
cat("- Field boundaries columns:", paste(colnames(field_boundaries), collapse = ", "), "\n")
|
|
cat("- Field boundaries rows:", nrow(field_boundaries), "\n")
|
|
|
|
# Join harvest data with field boundaries on field and sub_field
|
|
# Convert field and sub_field to character to ensure consistent matching
|
|
# Filter for current season only (where season_end is today's date)
|
|
current_date <- Sys.Date()
|
|
|
|
harvest_data_clean <- harvest_data %>%
|
|
mutate(
|
|
field = as.character(field),
|
|
sub_field = as.character(sub_field)
|
|
) %>%
|
|
# Filter for current season only
|
|
filter(season_end == current_date | (season_end >= current_date - 7 & season_end <= current_date)) %>%
|
|
select(field, sub_field, age, season_start, season_end, tonnage_ha)
|
|
|
|
cat("- Harvest data after filtering for current season:", nrow(harvest_data_clean), "rows\n")
|
|
|
|
if (nrow(harvest_data_clean) == 0) {
|
|
cat("⚠️ No current season data found in harvest data\n")
|
|
cat("- Looking for season_end near:", current_date, "\n")
|
|
|
|
# Show available season_end dates to help debug
|
|
available_dates <- harvest_data %>%
|
|
select(season_end) %>%
|
|
distinct() %>%
|
|
arrange(season_end)
|
|
|
|
cat("- Available season_end dates in harvest data:\n")
|
|
print(available_dates)
|
|
|
|
return(field_boundaries)
|
|
}
|
|
|
|
field_boundaries_clean <- field_boundaries %>%
|
|
mutate(
|
|
field = as.character(field),
|
|
sub_field = as.character(sub_field)
|
|
)
|
|
|
|
# Perform the join
|
|
field_boundaries_with_age <- field_boundaries_clean %>%
|
|
left_join(harvest_data_clean, by = c("field", "sub_field"))
|
|
|
|
# Check join success
|
|
matched_fields <- sum(!is.na(field_boundaries_with_age$age))
|
|
total_fields <- nrow(field_boundaries_with_age)
|
|
|
|
cat("✓ Successfully joined harvest data\n")
|
|
cat("- Fields with age data:", matched_fields, "out of", total_fields, "\n")
|
|
|
|
if (matched_fields > 0) {
|
|
age_summary <- field_boundaries_with_age %>%
|
|
filter(!is.na(age)) %>%
|
|
pull(age) %>%
|
|
summary()
|
|
|
|
cat("- Age range (weeks):", paste(names(age_summary), "=", round(age_summary, 1), collapse = ", "), "\n")
|
|
|
|
# Convert age from weeks to months for analysis
|
|
field_boundaries_with_age <- field_boundaries_with_age %>%
|
|
mutate(age_months = round(age / 4.33, 1)) # 4.33 weeks per month average
|
|
|
|
cat("- Age range (months):", paste(round(range(field_boundaries_with_age$age_months, na.rm = TRUE), 1), collapse = " to "), "\n")
|
|
}
|
|
|
|
return(field_boundaries_with_age)
|
|
|
|
}, error = function(e) {
|
|
warning(paste("Error processing harvest data:", e$message))
|
|
cat("Returning original field boundaries without age information\n")
|
|
return(field_boundaries)
|
|
})
|
|
}
|
|
|
|
#' Detect germination issues in young fields
|
|
#' @param indices_list List of vegetation index rasters
|
|
#' @param field_boundaries SF object with field boundaries (with age info)
|
|
#' @param young_fields_only Logical: analyze only young fields?
|
|
#' @return List with germination analysis results
|
|
detect_germination_issues <- function(indices_list, field_boundaries, young_fields_only = TRUE) {
|
|
|
|
cat("=== GERMINATION ANALYSIS ===\n")
|
|
|
|
germination_results <- list()
|
|
field_boundaries_vect <- terra::vect(field_boundaries)
|
|
|
|
for (i in seq_len(nrow(field_boundaries))) {
|
|
field_name <- field_boundaries$field[i]
|
|
sub_field_name <- field_boundaries$sub_field[i]
|
|
field_id <- paste0(field_name, "_", sub_field_name)
|
|
|
|
# Get field age if available
|
|
field_age_months <- if ("age_months" %in% colnames(field_boundaries) && !is.na(field_boundaries$age_months[i])) {
|
|
field_boundaries$age_months[i]
|
|
} else {
|
|
NA
|
|
}
|
|
|
|
cat("Analyzing germination for field:", field_id)
|
|
if (!is.na(field_age_months)) {
|
|
cat(" (Age:", field_age_months, "months)")
|
|
}
|
|
cat("\n")
|
|
|
|
# Extract field boundary
|
|
field_vect <- field_boundaries_vect[i]
|
|
|
|
# Analyze each index for germination patterns
|
|
field_analysis <- list()
|
|
|
|
for (index_name in names(indices_list)) {
|
|
index_raster <- indices_list[[index_name]]
|
|
|
|
# Extract values for this field
|
|
field_values <- terra::extract(index_raster, field_vect, fun = NULL)
|
|
valid_values <- unlist(field_values)
|
|
valid_values <- valid_values[!is.na(valid_values) & is.finite(valid_values)]
|
|
|
|
if (length(valid_values) > 10) {
|
|
# Calculate germination-relevant metrics
|
|
mean_val <- mean(valid_values)
|
|
cv_val <- sd(valid_values) / mean_val
|
|
min_val <- min(valid_values)
|
|
|
|
# Only check for germination issues in fields ≤ 4 months old
|
|
check_germination <- is.na(field_age_months) || field_age_months <= GERMINATION_PHASE_MAX
|
|
|
|
if (check_germination) {
|
|
# Detect potential germination issues
|
|
high_variation <- cv_val > GERMINATION_GAP_CV_THRESHOLD
|
|
low_values <- mean_val < LOW_CI_GERMINATION_THRESHOLD
|
|
very_low_areas <- sum(valid_values < (mean_val - 2*sd(valid_values))) / length(valid_values) * 100
|
|
} else {
|
|
# Field too old for germination analysis
|
|
high_variation <- FALSE
|
|
low_values <- FALSE
|
|
very_low_areas <- 0
|
|
}
|
|
|
|
field_analysis[[index_name]] <- list(
|
|
mean = mean_val,
|
|
cv = cv_val,
|
|
min = min_val,
|
|
high_variation = high_variation,
|
|
low_values = low_values,
|
|
very_low_areas_pct = very_low_areas,
|
|
n_pixels = length(valid_values),
|
|
age_months = field_age_months,
|
|
germination_analysis_applied = check_germination
|
|
)
|
|
}
|
|
}
|
|
|
|
germination_results[[field_id]] <- field_analysis
|
|
}
|
|
|
|
return(germination_results)
|
|
}
|
|
|
|
#' Detect weed pressure using change analysis
|
|
#' @param current_indices List of current week indices
|
|
#' @param previous_indices List of previous week indices
|
|
#' @param field_boundaries SF object with field boundaries
|
|
#' @return List with weed detection results
|
|
detect_weed_pressure <- function(current_indices, previous_indices, field_boundaries) {
|
|
|
|
cat("=== WEED PRESSURE DETECTION ===\n")
|
|
|
|
weed_results <- list()
|
|
field_boundaries_vect <- terra::vect(field_boundaries)
|
|
|
|
for (i in seq_len(nrow(field_boundaries))) {
|
|
field_name <- field_boundaries$field[i]
|
|
sub_field_name <- field_boundaries$sub_field[i]
|
|
field_id <- paste0(field_name, "_", sub_field_name)
|
|
|
|
cat("Analyzing weed pressure for field:", field_id, "\n")
|
|
|
|
field_vect <- field_boundaries_vect[i]
|
|
field_weed_analysis <- list()
|
|
|
|
# Analyze change in each index
|
|
for (index_name in names(current_indices)) {
|
|
|
|
if (index_name %in% names(previous_indices)) {
|
|
|
|
current_raster <- current_indices[[index_name]]
|
|
previous_raster <- previous_indices[[index_name]]
|
|
|
|
# Extract values for both weeks
|
|
current_values <- unlist(terra::extract(current_raster, field_vect, fun = NULL))
|
|
previous_values <- unlist(terra::extract(previous_raster, field_vect, fun = NULL))
|
|
|
|
# Clean values
|
|
valid_idx <- !is.na(current_values) & !is.na(previous_values) &
|
|
is.finite(current_values) & is.finite(previous_values)
|
|
current_clean <- current_values[valid_idx]
|
|
previous_clean <- previous_values[valid_idx]
|
|
|
|
if (length(current_clean) > 10) {
|
|
|
|
# Calculate change metrics
|
|
change_values <- current_clean - previous_clean
|
|
mean_change <- mean(change_values)
|
|
change_cv <- sd(change_values) / abs(mean(current_clean))
|
|
|
|
# Weed detection criteria
|
|
# 1. Significant positive change (weeds growing faster than cane)
|
|
rapid_increase <- mean_change >= WEED_CI_CHANGE_THRESHOLD
|
|
severe_increase <- mean_change >= SEVERE_WEED_CI_CHANGE
|
|
|
|
# 2. Percentage of field with rapid increase
|
|
rapid_increase_pct <- sum(change_values >= WEED_CI_CHANGE_THRESHOLD) / length(change_values) * 100
|
|
|
|
# 3. Patchy vs uniform weed patterns
|
|
patchy_weeds <- change_cv > 0.5 && rapid_increase_pct > 10 && rapid_increase_pct < 80
|
|
uniform_weeds <- change_cv < 0.3 && rapid_increase_pct > 60 # Whole field weeds
|
|
|
|
# 4. Current vegetation level (high CI + rapid change = likely weeds)
|
|
current_mean <- mean(current_clean)
|
|
high_current_ci <- current_mean > 2.0 # Adjust threshold as needed
|
|
|
|
field_weed_analysis[[index_name]] <- list(
|
|
mean_change = mean_change,
|
|
change_cv = change_cv,
|
|
current_mean = current_mean,
|
|
rapid_increase = rapid_increase,
|
|
severe_increase = severe_increase,
|
|
rapid_increase_pct = rapid_increase_pct,
|
|
patchy_weeds = patchy_weeds,
|
|
uniform_weeds = uniform_weeds,
|
|
high_current_ci = high_current_ci,
|
|
n_pixels = length(current_clean)
|
|
)
|
|
}
|
|
}
|
|
}
|
|
|
|
weed_results[[field_id]] <- field_weed_analysis
|
|
}
|
|
|
|
return(weed_results)
|
|
}
|
|
|
|
#' Generate intervention recommendations based on detected issues
|
|
#' @param field_id Field identifier
|
|
#' @param germination_analysis Germination analysis results
|
|
#' @param weed_analysis Weed detection results
|
|
#' @param field_age_months Field age in months (if available)
|
|
#' @return List with recommendations
|
|
generate_young_field_recommendations <- function(field_id, germination_analysis, weed_analysis, field_age_months = NA) {
|
|
|
|
recommendations <- list()
|
|
|
|
# Priority level: 1=urgent, 2=high, 3=moderate, 4=monitoring
|
|
priority <- 4
|
|
messages <- c()
|
|
interventions <- c()
|
|
|
|
# Check for germination issues (if CI analysis available and field is young enough)
|
|
if ("CI" %in% names(germination_analysis)) {
|
|
ci_analysis <- germination_analysis[["CI"]]
|
|
|
|
# Only apply germination analysis for fields ≤ 4 months old
|
|
if (!is.null(ci_analysis$germination_analysis_applied) && ci_analysis$germination_analysis_applied) {
|
|
|
|
if (ci_analysis$high_variation && ci_analysis$low_values) {
|
|
priority <- min(priority, 1)
|
|
messages <- c(messages, "🚨 URGENT: Poor germination detected - high variation with low CI values")
|
|
interventions <- c(interventions, "Immediate replanting in gap areas", "Check seed quality and planting conditions")
|
|
} else if (ci_analysis$high_variation) {
|
|
priority <- min(priority, 2)
|
|
messages <- c(messages, "⚠️ Alert: Uneven emergence detected - moderate intervention needed")
|
|
interventions <- c(interventions, "Monitor gap areas for potential replanting")
|
|
}
|
|
|
|
if (ci_analysis$very_low_areas_pct > 15) {
|
|
priority <- min(priority, 2)
|
|
messages <- c(messages, paste0("⚠️ Alert: ", round(ci_analysis$very_low_areas_pct, 1), "% of field has very low vegetation"))
|
|
interventions <- c(interventions, "Investigate cause of low-performing areas")
|
|
}
|
|
|
|
} else {
|
|
# Field too old for germination analysis - check for other issues
|
|
if (ci_analysis$cv > 0.5) {
|
|
priority <- min(priority, 3)
|
|
messages <- c(messages, "📊 Info: Uneven crop growth detected (field too old for replanting)")
|
|
interventions <- c(interventions, "Monitor for potential management issues")
|
|
}
|
|
}
|
|
}
|
|
|
|
# Check for weed pressure (if CI analysis available)
|
|
if ("CI" %in% names(weed_analysis)) {
|
|
ci_weed <- weed_analysis[["CI"]]
|
|
|
|
if (ci_weed$severe_increase) {
|
|
priority <- min(priority, 1)
|
|
messages <- c(messages, paste0("🚨 CRITICAL: Severe weed pressure detected - CI increased by ", round(ci_weed$mean_change, 2)))
|
|
interventions <- c(interventions, "Immediate weed control required", "Consider herbicide application")
|
|
} else if (ci_weed$rapid_increase) {
|
|
priority <- min(priority, 2)
|
|
messages <- c(messages, paste0("⚠️ Alert: Weed pressure detected - CI increased by ", round(ci_weed$mean_change, 2)))
|
|
interventions <- c(interventions, "Schedule weed control within 1-2 weeks")
|
|
}
|
|
|
|
if (ci_weed$uniform_weeds) {
|
|
priority <- min(priority, 2)
|
|
messages <- c(messages, paste0("🔶 Pattern: Uniform weed emergence across ", round(ci_weed$rapid_increase_pct, 1), "% of field"))
|
|
interventions <- c(interventions, "Broad-scale weed management needed")
|
|
} else if (ci_weed$patchy_weeds) {
|
|
priority <- min(priority, 3)
|
|
messages <- c(messages, paste0("🔶 Pattern: Patchy weed emergence in ", round(ci_weed$rapid_increase_pct, 1), "% of field"))
|
|
interventions <- c(interventions, "Targeted spot treatment recommended")
|
|
}
|
|
}
|
|
|
|
# Age-specific recommendations
|
|
if (!is.na(field_age_months)) {
|
|
if (field_age_months <= GERMINATION_PHASE_MAX) {
|
|
interventions <- c(interventions, "Critical germination phase - maintain optimal moisture", "Monitor for pest damage")
|
|
} else if (field_age_months <= WEED_CRITICAL_PHASE_MAX) {
|
|
interventions <- c(interventions, "Peak weed competition period - prioritize weed control")
|
|
}
|
|
}
|
|
|
|
# Default monitoring if no issues
|
|
if (length(messages) == 0) {
|
|
messages <- c("✅ No significant issues detected")
|
|
interventions <- c("Continue routine monitoring")
|
|
}
|
|
|
|
return(list(
|
|
field_id = field_id,
|
|
priority = priority,
|
|
messages = messages,
|
|
interventions = interventions,
|
|
requires_action = priority <= 3
|
|
))
|
|
}
|
|
|
|
# 4. Main analysis function
|
|
# -----------------------
|
|
main_young_field_analysis <- function() {
|
|
|
|
# Get command line arguments
|
|
args <- commandArgs(trailingOnly = TRUE)
|
|
|
|
current_week <- if (length(args) >= 1) as.numeric(args[1]) else 30
|
|
previous_week <- if (length(args) >= 2) as.numeric(args[2]) else 29
|
|
project_dir <- if (length(args) >= 3) args[3] else "simba"
|
|
|
|
cat("=== YOUNG FIELD ANALYSIS SYSTEM ===\n")
|
|
cat("Project:", project_dir, "\n")
|
|
cat("Analyzing weeks:", previous_week, "→", current_week, "\n\n")
|
|
|
|
# Load project configuration
|
|
assign("project_dir", project_dir, envir = .GlobalEnv)
|
|
|
|
tryCatch({
|
|
source("parameters_project.R")
|
|
cat("✓ Project configuration loaded\n")
|
|
}, error = function(e) {
|
|
tryCatch({
|
|
source(here::here("r_app", "parameters_project.R"))
|
|
cat("✓ Project configuration loaded from r_app directory\n")
|
|
}, error = function(e) {
|
|
stop("Failed to load project configuration")
|
|
})
|
|
})
|
|
|
|
# Verify required variables
|
|
if (!exists("weekly_CI_mosaic") || !exists("field_boundaries_sf")) {
|
|
stop("Required project variables not found")
|
|
}
|
|
|
|
# Check if harvest data is available from parameters_project.R
|
|
if (exists("harvesting_data") && !is.null(harvesting_data)) {
|
|
cat("✓ Harvest data loaded from parameters_project.R\n")
|
|
field_boundaries_with_age <- load_field_ages_from_existing(harvesting_data, field_boundaries_sf)
|
|
} else {
|
|
cat("⚠️ No harvest data found in parameters_project.R\n")
|
|
field_boundaries_with_age <- field_boundaries_sf
|
|
}
|
|
|
|
# Construct mosaic file paths
|
|
year <- 2025
|
|
current_week_file <- sprintf("week_%02d_%d.tif", current_week, year)
|
|
previous_week_file <- sprintf("week_%02d_%d.tif", previous_week, year)
|
|
|
|
current_mosaic_path <- file.path(weekly_CI_mosaic, current_week_file)
|
|
previous_mosaic_path <- file.path(weekly_CI_mosaic, previous_week_file)
|
|
|
|
cat("Current week mosaic:", current_mosaic_path, "\n")
|
|
cat("Previous week mosaic:", previous_mosaic_path, "\n")
|
|
|
|
# Check if files exist
|
|
if (!file.exists(current_mosaic_path)) {
|
|
stop("Current week mosaic not found: ", current_mosaic_path)
|
|
}
|
|
if (!file.exists(previous_mosaic_path)) {
|
|
stop("Previous week mosaic not found: ", previous_mosaic_path)
|
|
}
|
|
|
|
# Load mosaics and extract bands
|
|
cat("\nLoading mosaics and extracting RGBNIR bands...\n")
|
|
|
|
current_mosaic <- terra::rast(current_mosaic_path)
|
|
previous_mosaic <- terra::rast(previous_mosaic_path)
|
|
|
|
cat("Current mosaic bands:", nlyr(current_mosaic), "\n")
|
|
cat("Previous mosaic bands:", nlyr(previous_mosaic), "\n")
|
|
|
|
# Extract RGBNIR bands (assuming standard order: R=1, G=2, B=3, NIR=4)
|
|
current_red <- current_mosaic[[1]]
|
|
current_green <- current_mosaic[[2]]
|
|
current_blue <- current_mosaic[[3]]
|
|
current_nir <- current_mosaic[[4]]
|
|
|
|
previous_red <- previous_mosaic[[1]]
|
|
previous_green <- previous_mosaic[[2]]
|
|
previous_blue <- previous_mosaic[[3]]
|
|
previous_nir <- previous_mosaic[[4]]
|
|
|
|
# Calculate enhanced indices for both weeks
|
|
cat("\nCalculating vegetation indices...\n")
|
|
current_indices <- calculate_enhanced_indices(current_red, current_green, current_blue, current_nir)
|
|
previous_indices <- calculate_enhanced_indices(previous_red, previous_green, previous_blue, previous_nir)
|
|
|
|
# Create output directory for maps
|
|
output_dir <- file.path(dirname(weekly_CI_mosaic), "young_field_analysis")
|
|
|
|
# Create visual comparison maps
|
|
cat("\nCreating visual maps...\n")
|
|
# Check if maps already exist before creating them
|
|
current_week_maps_exist <- all(sapply(names(current_indices), function(index_name) {
|
|
map_file <- file.path(output_dir, paste0("week_", sprintf("%02d", current_week), "_", index_name, "_map.png"))
|
|
tif_file <- file.path(output_dir, paste0("week_", sprintf("%02d", current_week), "_", index_name, ".tif"))
|
|
file.exists(map_file) && file.exists(tif_file)
|
|
}))
|
|
|
|
previous_week_maps_exist <- all(sapply(names(previous_indices), function(index_name) {
|
|
map_file <- file.path(output_dir, paste0("week_", sprintf("%02d", previous_week), "_", index_name, "_map.png"))
|
|
tif_file <- file.path(output_dir, paste0("week_", sprintf("%02d", previous_week), "_", index_name, ".tif"))
|
|
file.exists(map_file) && file.exists(tif_file)
|
|
}))
|
|
|
|
if (!current_week_maps_exist) {
|
|
cat("Creating current week maps (week", current_week, ")...\n")
|
|
create_index_comparison_maps(current_indices, field_boundaries_sf, output_dir, current_week)
|
|
} else {
|
|
cat("Current week maps already exist, skipping creation...\n")
|
|
}
|
|
|
|
if (!previous_week_maps_exist) {
|
|
cat("Creating previous week maps (week", previous_week, ")...\n")
|
|
create_index_comparison_maps(previous_indices, field_boundaries_sf, output_dir, previous_week)
|
|
} else {
|
|
cat("Previous week maps already exist, skipping creation...\n")
|
|
}
|
|
|
|
# Perform germination analysis
|
|
cat("\nPerforming germination analysis...\n")
|
|
germination_results <- detect_germination_issues(current_indices, field_boundaries_with_age)
|
|
|
|
# Perform weed detection analysis
|
|
cat("\nPerforming weed pressure analysis...\n")
|
|
weed_results <- detect_weed_pressure(current_indices, previous_indices, field_boundaries_with_age)
|
|
|
|
# Generate recommendations for each field
|
|
cat("\n=== FIELD RECOMMENDATIONS ===\n")
|
|
|
|
field_recommendations <- list()
|
|
action_needed_fields <- 0
|
|
|
|
for (field_id in names(germination_results)) {
|
|
|
|
germination_analysis <- germination_results[[field_id]]
|
|
weed_analysis <- if (field_id %in% names(weed_results)) weed_results[[field_id]] else list()
|
|
|
|
# Get field age if available
|
|
field_row <- field_boundaries_with_age[field_boundaries_with_age$field == strsplit(field_id, "_")[[1]][1] &
|
|
field_boundaries_with_age$sub_field == strsplit(field_id, "_")[[1]][2], ]
|
|
field_age_months <- if (nrow(field_row) > 0 && !is.na(field_row$age_months[1])) field_row$age_months[1] else NA
|
|
|
|
recommendations <- generate_young_field_recommendations(
|
|
field_id,
|
|
germination_analysis,
|
|
weed_analysis,
|
|
field_age_months
|
|
)
|
|
|
|
field_recommendations[[field_id]] <- recommendations
|
|
|
|
# Print recommendations
|
|
cat("\nFIELD:", field_id, "\n")
|
|
if (!is.na(field_age_months)) {
|
|
cat("Age:", field_age_months, "months\n")
|
|
}
|
|
cat("Priority Level:", recommendations$priority, "\n")
|
|
|
|
for (message in recommendations$messages) {
|
|
cat("-", message, "\n")
|
|
}
|
|
|
|
if (length(recommendations$interventions) > 0) {
|
|
cat("Recommended Actions:\n")
|
|
for (intervention in recommendations$interventions) {
|
|
cat(" •", intervention, "\n")
|
|
}
|
|
}
|
|
|
|
if (recommendations$requires_action) {
|
|
action_needed_fields <- action_needed_fields + 1
|
|
}
|
|
|
|
# Show index comparison for this field
|
|
if ("CI" %in% names(germination_analysis)) {
|
|
ci_stats <- germination_analysis[["CI"]]
|
|
cat("CI Stats: Mean =", round(ci_stats$mean, 3),
|
|
", CV =", round(ci_stats$cv, 3),
|
|
", Low areas =", round(ci_stats$very_low_areas_pct, 1), "%\n")
|
|
}
|
|
|
|
if ("CI" %in% names(weed_analysis)) {
|
|
weed_stats <- weed_analysis[["CI"]]
|
|
cat("Change Stats: CI Δ =", round(weed_stats$mean_change, 3),
|
|
", Rapid increase =", round(weed_stats$rapid_increase_pct, 1), "%\n")
|
|
}
|
|
}
|
|
|
|
# Summary
|
|
cat("\n=== ANALYSIS SUMMARY ===\n")
|
|
cat("Total fields analyzed:", length(field_recommendations), "\n")
|
|
cat("Fields requiring action:", action_needed_fields, "\n")
|
|
cat("Maps saved to:", output_dir, "\n")
|
|
|
|
cat("\nIndex comparison maps created for visual inspection:\n")
|
|
cat("- NDVI: Standard vegetation index\n")
|
|
cat("- GNDVI: Green NDVI (sensitive to early vegetation)\n")
|
|
cat("- ExG: Excess Green (early detection)\n")
|
|
cat("- VARI: Atmospherically resistant\n")
|
|
cat("- GRVI: Green ratio\n")
|
|
cat("- CI: Chlorophyll Index (current standard)\n")
|
|
|
|
cat("\n✓ Young field analysis complete\n")
|
|
|
|
return(list(
|
|
germination_results = germination_results,
|
|
weed_results = weed_results,
|
|
recommendations = field_recommendations,
|
|
indices_calculated = names(current_indices),
|
|
maps_created = file.path(output_dir)
|
|
))
|
|
}
|
|
|
|
# Run analysis if script called directly
|
|
if (sys.nframe() == 0) {
|
|
main_young_field_analysis()
|
|
} |