369 lines
12 KiB
R
369 lines
12 KiB
R
# SAR Visualization Script
|
|
# ========================
|
|
#
|
|
# Functions for creating maps and plots of Sentinel-1 SAR data
|
|
# Integrates with existing tmap and ggplot2 workflow
|
|
#
|
|
# Author: Timon
|
|
# Date: August 2025
|
|
|
|
# Load required libraries
|
|
# source("sar_analysis_functions.R") # Don't source here, will be sourced in main script
|
|
|
|
# Load additional plotting libraries
|
|
if (!require(RColorBrewer)) install.packages("RColorBrewer")
|
|
if (!require(scales)) install.packages("scales")
|
|
if (!require(patchwork)) install.packages("patchwork")
|
|
|
|
library(RColorBrewer)
|
|
library(scales)
|
|
library(patchwork)
|
|
|
|
# SAR Mapping Functions
|
|
# ====================
|
|
|
|
#' Create SAR backscatter map
|
|
#'
|
|
#' @param sar_data SAR raster data
|
|
#' @param field_boundaries Field boundary polygons
|
|
#' @param title Map title
|
|
#' @param band SAR band name for legend
|
|
#' @param interactive Create interactive map (default TRUE)
|
|
#' @return tmap object
|
|
create_sar_map <- function(sar_data, field_boundaries = NULL, title = "SAR Backscatter",
|
|
band = "VV_dB", interactive = TRUE) {
|
|
|
|
# Create base map
|
|
map_plot <- tm_shape(sar_data) +
|
|
tm_raster(
|
|
title = paste(band, "(dB)"),
|
|
palette = "viridis",
|
|
style = "cont",
|
|
n = 10,
|
|
alpha = 0.8
|
|
) +
|
|
tm_layout(
|
|
title = title,
|
|
title.size = 1.2
|
|
)
|
|
|
|
# Add field boundaries if provided
|
|
if (!is.null(field_boundaries)) {
|
|
map_plot <- map_plot +
|
|
tm_shape(field_boundaries) +
|
|
tm_borders(col = "white", lwd = 2, alpha = 0.9) +
|
|
tm_fill(alpha = 0.1, col = "lightblue") +
|
|
tm_text("Name", size = 0.7, col = "black", bg.color = "white", bg.alpha = 0.7)
|
|
}
|
|
|
|
# Add interactive features if in view mode
|
|
if (interactive && tmap_mode() == "view") {
|
|
if (!is.null(field_boundaries)) {
|
|
map_plot <- map_plot +
|
|
tm_shape(field_boundaries) +
|
|
tm_polygons(
|
|
alpha = 0.2,
|
|
col = "lightblue",
|
|
border.col = "white",
|
|
border.lwd = 2,
|
|
popup.vars = c("Field Name" = "Name", "Area (ha)" = "Area_ha"),
|
|
id = "field_label"
|
|
)
|
|
}
|
|
}
|
|
|
|
return(map_plot)
|
|
}
|
|
|
|
#' Create multi-temporal RGB SAR visualization
|
|
#' Your idea: Week1=Red, Week2=Green, Week3=Blue for change detection
|
|
#'
|
|
#' @param week1_data SAR data for week 1 (Red channel)
|
|
#' @param week2_data SAR data for week 2 (Green channel)
|
|
#' @param week3_data SAR data for week 3 (Blue channel)
|
|
#' @param field_boundaries Field boundaries
|
|
#' @return tmap object
|
|
create_sar_rgb_change_map <- function(week1_data, week2_data, week3_data,
|
|
field_boundaries = NULL) {
|
|
|
|
cat("Creating SAR RGB change detection map...\n")
|
|
|
|
# Normalize each week to 0-255 range for RGB
|
|
week1_norm <- (week1_data - global(week1_data, min, na.rm = TRUE)[1,1]) /
|
|
(global(week1_data, max, na.rm = TRUE)[1,1] - global(week1_data, min, na.rm = TRUE)[1,1]) * 255
|
|
|
|
week2_norm <- (week2_data - global(week2_data, min, na.rm = TRUE)[1,1]) /
|
|
(global(week2_data, max, na.rm = TRUE)[1,1] - global(week2_data, min, na.rm = TRUE)[1,1]) * 255
|
|
|
|
week3_norm <- (week3_data - global(week3_data, min, na.rm = TRUE)[1,1]) /
|
|
(global(week3_data, max, na.rm = TRUE)[1,1] - global(week3_data, min, na.rm = TRUE)[1,1]) * 255
|
|
|
|
# Stack for RGB
|
|
rgb_stack <- c(week1_norm, week2_norm, week3_norm)
|
|
names(rgb_stack) <- c("Red", "Green", "Blue")
|
|
|
|
# Create RGB map
|
|
map_plot <- tm_shape(rgb_stack) +
|
|
tm_rgb(
|
|
r = 1, g = 2, b = 3,
|
|
alpha = 0.8
|
|
) +
|
|
tm_layout(
|
|
title = "SAR Change Detection (RGB)\nRed=Week1, Green=Week2, Blue=Week3",
|
|
title.size = 1.0
|
|
)
|
|
|
|
# Add field boundaries
|
|
if (!is.null(field_boundaries)) {
|
|
map_plot <- map_plot +
|
|
tm_shape(field_boundaries) +
|
|
tm_borders(col = "white", lwd = 0.5, alpha = 0.7)
|
|
}
|
|
|
|
return(map_plot)
|
|
}
|
|
|
|
#' Create SAR comparison map (side-by-side weeks)
|
|
#'
|
|
#' @param sar_data1 First SAR dataset
|
|
#' @param sar_data2 Second SAR dataset
|
|
#' @param title1 Title for first map
|
|
#' @param title2 Title for second map
|
|
#' @param field_boundaries Field boundaries
|
|
#' @return Combined tmap
|
|
create_sar_comparison_map <- function(sar_data1, sar_data2, title1 = "Week 1",
|
|
title2 = "Week 2", field_boundaries = NULL) {
|
|
|
|
# Create individual maps
|
|
map1 <- create_sar_map(sar_data1, field_boundaries, title1)
|
|
map2 <- create_sar_map(sar_data2, field_boundaries, title2)
|
|
|
|
# Combine maps
|
|
combined_map <- tmap_arrange(map1, map2, ncol = 2)
|
|
|
|
return(combined_map)
|
|
}
|
|
|
|
#' Create SAR metrics map (RVI, Cross-pol ratio, etc.)
|
|
#'
|
|
#' @param metric_data Calculated SAR metric raster
|
|
#' @param metric_name Name of the metric
|
|
#' @param field_boundaries Field boundaries
|
|
#' @param color_palette Color palette to use
|
|
#' @param interactive Create interactive map
|
|
#' @return tmap object
|
|
create_sar_metrics_map <- function(metric_data, metric_name = "SAR Metric",
|
|
field_boundaries = NULL, color_palette = "RdYlGn",
|
|
interactive = TRUE) {
|
|
|
|
# Create map
|
|
map_plot <- tm_shape(metric_data) +
|
|
tm_raster(
|
|
title = metric_name,
|
|
palette = color_palette,
|
|
style = "quantile",
|
|
n = 8,
|
|
alpha = 0.8
|
|
) +
|
|
tm_layout(
|
|
title = paste("SAR", metric_name),
|
|
title.size = 1.2
|
|
)
|
|
|
|
# Add field boundaries with enhanced interactivity
|
|
if (!is.null(field_boundaries)) {
|
|
if (interactive && tmap_mode() == "view") {
|
|
map_plot <- map_plot +
|
|
tm_shape(field_boundaries) +
|
|
tm_polygons(
|
|
alpha = 0.2,
|
|
col = "lightblue",
|
|
border.col = "white",
|
|
border.lwd = 2,
|
|
popup.vars = c("Field Name" = "Name"),
|
|
id = "field_label"
|
|
)
|
|
} else {
|
|
map_plot <- map_plot +
|
|
tm_shape(field_boundaries) +
|
|
tm_borders(col = "white", lwd = 2, alpha = 0.9) +
|
|
tm_text("Name", size = 0.6, col = "black", bg.color = "white", bg.alpha = 0.7)
|
|
}
|
|
}
|
|
|
|
return(map_plot)
|
|
}
|
|
|
|
# SAR Time Series Plotting Functions
|
|
# ==================================
|
|
|
|
#' Plot SAR time series for multiple fields
|
|
#'
|
|
#' @param field_ids Vector of field IDs to plot
|
|
#' @param field_boundaries Field boundary polygons
|
|
#' @param week_range Range of weeks
|
|
#' @param band SAR band to plot
|
|
#' @param facet_fields Create separate panels for each field
|
|
#' @return ggplot object
|
|
plot_sar_timeseries <- function(field_ids, field_boundaries, week_range = 26:33,
|
|
band = "VV_dB_filtered", facet_fields = TRUE) {
|
|
|
|
cat("Creating SAR time series plot for", length(field_ids), "fields...\n")
|
|
|
|
# Collect data for all fields
|
|
all_data <- data.frame()
|
|
|
|
for (field_id in field_ids) {
|
|
field_data <- analyze_sar_field_trends(field_id, field_boundaries, week_range, band)
|
|
if (!is.null(field_data)) {
|
|
all_data <- rbind(all_data, field_data$timeseries)
|
|
}
|
|
}
|
|
|
|
if (nrow(all_data) == 0) {
|
|
stop("No data found for specified fields")
|
|
}
|
|
|
|
# Create plot
|
|
p <- ggplot(all_data, aes(x = week, y = value, color = field_id)) +
|
|
geom_line(size = 1.2, alpha = 0.8) +
|
|
geom_point(size = 2.5, alpha = 0.9) +
|
|
labs(
|
|
title = paste("SAR Time Series -", band),
|
|
subtitle = paste("Weeks", min(week_range), "to", max(week_range), "2025"),
|
|
x = "Week Number",
|
|
y = paste(band, "Backscatter (dB)"),
|
|
color = "Field ID"
|
|
) +
|
|
theme_minimal() +
|
|
theme(
|
|
plot.title = element_text(size = 14, face = "bold"),
|
|
plot.subtitle = element_text(size = 12),
|
|
axis.text = element_text(size = 10),
|
|
axis.title = element_text(size = 12),
|
|
legend.title = element_text(size = 11),
|
|
legend.text = element_text(size = 10)
|
|
) +
|
|
scale_color_viridis_d(option = "plasma") +
|
|
scale_x_continuous(breaks = week_range)
|
|
|
|
# Add faceting if requested
|
|
if (facet_fields && length(field_ids) > 1) {
|
|
p <- p + facet_wrap(~field_id, scales = "free_y", ncol = 2)
|
|
}
|
|
|
|
return(p)
|
|
}
|
|
|
|
#' Plot SAR vs Optical CI comparison
|
|
#'
|
|
#' @param field_id Field to analyze
|
|
#' @param field_boundaries Field boundaries
|
|
#' @param week_range Week range
|
|
#' @param ci_data_dir Directory with CI mosaics
|
|
#' @return ggplot object
|
|
plot_sar_vs_ci_comparison <- function(field_id, field_boundaries, week_range = 26:33,
|
|
ci_data_dir = "data/aura/weekly_CI_mosaic") {
|
|
|
|
cat("Creating SAR vs CI comparison for field:", field_id, "\n")
|
|
|
|
# Get SAR data
|
|
sar_data <- analyze_sar_field_trends(field_id, field_boundaries, week_range, "VV_dB_filtered")
|
|
|
|
# Get CI data (adapt this to your CI loading function)
|
|
ci_data <- data.frame()
|
|
field_poly <- field_boundaries[field_boundaries$Name == field_id, ]
|
|
|
|
for (week in week_range) {
|
|
ci_file <- file.path(ci_data_dir, paste0("week_", sprintf("%02d", week), "_2025.tif"))
|
|
if (file.exists(ci_file)) {
|
|
ci_raster <- rast(ci_file)
|
|
ci_value <- extract(ci_raster, field_poly, fun = mean, na.rm = TRUE)[1, 2]
|
|
ci_data <- rbind(ci_data, data.frame(week = week, ci_value = ci_value))
|
|
}
|
|
}
|
|
|
|
# Combine datasets
|
|
if (nrow(ci_data) > 0 && !is.null(sar_data)) {
|
|
combined_data <- merge(sar_data$timeseries, ci_data, by = "week", all = TRUE)
|
|
|
|
# Create dual-axis plot
|
|
p1 <- ggplot(combined_data, aes(x = week)) +
|
|
geom_line(aes(y = value, color = "SAR VV_dB"), size = 1.2) +
|
|
geom_point(aes(y = value, color = "SAR VV_dB"), size = 2.5) +
|
|
labs(
|
|
title = paste("SAR vs CI Comparison - Field", field_id),
|
|
x = "Week Number",
|
|
y = "SAR Backscatter (dB)",
|
|
color = "Data Type"
|
|
) +
|
|
theme_minimal() +
|
|
scale_color_manual(values = c("SAR VV_dB" = "#440154")) +
|
|
scale_x_continuous(breaks = week_range)
|
|
|
|
# Add CI on secondary axis (approximation)
|
|
if (!all(is.na(combined_data$ci_value))) {
|
|
p1 <- p1 +
|
|
geom_line(aes(y = scales::rescale(ci_value, to = range(value, na.rm = TRUE)),
|
|
color = "CI"), size = 1.2) +
|
|
geom_point(aes(y = scales::rescale(ci_value, to = range(value, na.rm = TRUE)),
|
|
color = "CI"), size = 2.5) +
|
|
scale_color_manual(values = c("SAR VV_dB" = "#440154", "CI" = "#35B779"))
|
|
}
|
|
|
|
return(p1)
|
|
} else {
|
|
warning("Insufficient data for comparison")
|
|
return(NULL)
|
|
}
|
|
}
|
|
|
|
#' Create SAR metrics dashboard
|
|
#'
|
|
#' @param week_num Week to analyze
|
|
#' @param field_boundaries Field boundaries
|
|
#' @param interactive Create interactive dashboard
|
|
#' @return List of individual maps for interactive display
|
|
create_sar_dashboard <- function(week_num = 33, field_boundaries, interactive = TRUE) {
|
|
|
|
cat("Creating SAR dashboard for week", week_num, "...\n")
|
|
|
|
# Load different SAR bands
|
|
vv_linear <- load_sar_mosaic(week_num, band = "VV")
|
|
vh_linear <- load_sar_mosaic(week_num, band = "VH")
|
|
vv_db <- load_sar_mosaic(week_num, band = "VV_dB_filtered")
|
|
vh_db <- load_sar_mosaic(week_num, band = "VH_dB_filtered")
|
|
|
|
# Calculate metrics
|
|
rvi <- calculate_rvi(vv_linear, vh_linear)
|
|
cross_ratio <- calculate_cross_pol_ratio(vv_db, vh_db)
|
|
csi <- calculate_crop_structure_index(vv_linear, vh_linear)
|
|
|
|
# Create individual maps for interactive display
|
|
map1 <- create_sar_map(vv_db, field_boundaries, "VV Backscatter", "VV_dB", interactive)
|
|
map2 <- create_sar_map(vh_db, field_boundaries, "VH Backscatter", "VH_dB", interactive)
|
|
map3 <- create_sar_metrics_map(rvi, "RVI", field_boundaries, "RdYlGn", interactive)
|
|
map4 <- create_sar_metrics_map(cross_ratio, "Cross-Pol Ratio", field_boundaries, "RdBu", interactive)
|
|
|
|
# Return list of maps for individual display in RMarkdown
|
|
return(list(
|
|
vv_map = map1,
|
|
vh_map = map2,
|
|
rvi_map = map3,
|
|
cross_ratio_map = map4,
|
|
rvi_data = rvi,
|
|
cross_ratio_data = cross_ratio
|
|
))
|
|
}
|
|
}
|
|
|
|
cat("SAR visualization functions loaded successfully!\n")
|
|
cat("Available functions:\n")
|
|
cat("- create_sar_map(): Basic SAR backscatter map\n")
|
|
cat("- create_sar_rgb_change_map(): Multi-temporal RGB change detection\n")
|
|
cat("- create_sar_comparison_map(): Side-by-side week comparison\n")
|
|
cat("- create_sar_metrics_map(): SAR indices visualization\n")
|
|
cat("- plot_sar_timeseries(): Time series plots for fields\n")
|
|
cat("- plot_sar_vs_ci_comparison(): SAR vs optical CI comparison\n")
|
|
cat("- create_sar_dashboard(): Complete SAR metrics dashboard\n")
|