491 lines
21 KiB
R
491 lines
21 KiB
R
# REPORT_UTILS.R
|
|
# =============
|
|
# Utility functions for generating SmartCane reports with visualizations.
|
|
# These functions support the creation of maps, charts and report elements
|
|
# for the CI_report_dashboard_planet.Rmd document.
|
|
|
|
#' 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)
|
|
}
|
|
}
|
|
}
|
|
|
|
#' Creates a sub-chunk for use within RMarkdown documents
|
|
#'
|
|
#' @param g A ggplot object to render in the sub-chunk
|
|
#' @param fig_height Height of the figure in inches
|
|
#' @param fig_width Width of the figure in inches
|
|
#' @return NULL (writes chunk directly to output)
|
|
#'
|
|
subchunkify <- function(g, fig_height=7, fig_width=5) {
|
|
g_deparsed <- paste0(deparse(
|
|
function() {g}
|
|
), collapse = '')
|
|
|
|
sub_chunk <- paste0("
|
|
`","``{r sub_chunk_", floor(runif(1) * 10000), ", fig.height=", fig_height, ", fig.width=", fig_width, ", echo=FALSE}",
|
|
"\n(",
|
|
g_deparsed
|
|
, ")()",
|
|
"\n`","``
|
|
")
|
|
|
|
cat(knitr::knit(text = knitr::knit_expand(text = sub_chunk), quiet = TRUE))
|
|
}
|
|
|
|
#' Creates a Chlorophyll Index map for a pivot
|
|
#'
|
|
#' @param pivot_raster The raster data containing CI values
|
|
#' @param pivot_shape The shape of the pivot field
|
|
#' @param pivot_spans Additional boundary data for the field
|
|
#' @param show_legend Whether to show the legend (default: FALSE)
|
|
#' @param legend_is_portrait Whether to show the legend in portrait orientation (default: FALSE)
|
|
#' @param week Week number to display in the title
|
|
#' @param age Age of the crop in weeks
|
|
#' @param borders Whether to display field borders (default: FALSE)
|
|
#' @return A tmap object with the CI map
|
|
#'
|
|
create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week, age, borders = FALSE){
|
|
# Input validation
|
|
if (missing(pivot_raster) || is.null(pivot_raster)) {
|
|
stop("pivot_raster is required")
|
|
}
|
|
if (missing(pivot_shape) || is.null(pivot_shape)) {
|
|
stop("pivot_shape is required")
|
|
}
|
|
if (missing(pivot_spans) || is.null(pivot_spans)) {
|
|
stop("pivot_spans is required")
|
|
}
|
|
if (missing(week) || is.null(week)) {
|
|
stop("week parameter is required")
|
|
}
|
|
if (missing(age) || is.null(age)) {
|
|
stop("age parameter is required")
|
|
} # Create the base map
|
|
map <- tm_shape(pivot_raster, unit = "m") # Add raster with continuous spectrum (fixed scale 1-8 for consistent comparison)
|
|
map <- map + tm_raster(col.scale = tm_scale_continuous(values = "brewer.rd_yl_gn",
|
|
limits = c(1, 8)),
|
|
col.legend = tm_legend(title = "CI",
|
|
orientation = if(legend_is_portrait) "portrait" else "landscape",
|
|
show = show_legend,
|
|
position = if(show_legend) tm_pos_out("left", "center") else c("left", "bottom")
|
|
))
|
|
# Add layout elements
|
|
map <- map + tm_layout(main.title = paste0("Max CI week ", week,"\n", age, " weeks old"),
|
|
main.title.size = 0.7)
|
|
|
|
# Add borders if requested
|
|
if (borders) {
|
|
map <- map +
|
|
tm_shape(pivot_shape) +
|
|
tm_borders(lwd = 3) +
|
|
tm_text("sub_field", size = 1/2) +
|
|
tm_shape(pivot_spans) +
|
|
tm_borders(lwd = 0.5, alpha = 0.5)
|
|
}
|
|
|
|
return(map)
|
|
}
|
|
|
|
#' Creates a Chlorophyll Index difference map between two weeks
|
|
#'
|
|
#' @param pivot_raster The raster data containing CI difference values
|
|
#' @param pivot_shape The shape of the pivot field
|
|
#' @param pivot_spans Additional boundary data for the field
|
|
#' @param show_legend Whether to show the legend (default: FALSE)
|
|
#' @param legend_is_portrait Whether to show the legend in portrait orientation (default: FALSE)
|
|
#' @param week_1 First week number for comparison
|
|
#' @param week_2 Second week number for comparison
|
|
#' @param age Age of the crop in weeks
|
|
#' @param borders Whether to display field borders (default: TRUE)
|
|
#' @return A tmap object with the CI difference map
|
|
#'
|
|
create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week_1, week_2, age, borders = TRUE){
|
|
# Input validation
|
|
if (missing(pivot_raster) || is.null(pivot_raster)) {
|
|
stop("pivot_raster is required")
|
|
}
|
|
if (missing(pivot_shape) || is.null(pivot_shape)) {
|
|
stop("pivot_shape is required")
|
|
}
|
|
if (missing(pivot_spans) || is.null(pivot_spans)) {
|
|
stop("pivot_spans is required")
|
|
}
|
|
if (missing(week_1) || is.null(week_1) || missing(week_2) || is.null(week_2)) {
|
|
stop("week_1 and week_2 parameters are required")
|
|
}
|
|
if (missing(age) || is.null(age)) {
|
|
stop("age parameter is required")
|
|
} # Create the base map
|
|
map <- tm_shape(pivot_raster, unit = "m") # Add raster with continuous spectrum (centered at 0 for difference maps, fixed scale)
|
|
map <- map + tm_raster(col.scale = tm_scale_continuous(values = "brewer.rd_yl_gn",
|
|
midpoint = 0,
|
|
limits = c(-3, 3)),
|
|
col.legend = tm_legend(title = "CI diff.",
|
|
orientation = if(legend_is_portrait) "portrait" else "landscape",
|
|
show = show_legend,
|
|
position = if(show_legend) tm_pos_out("right", "center") else c("left", "bottom")
|
|
))
|
|
# Add layout elements
|
|
map <- map + tm_layout(main.title = paste0("CI change week ", week_1, " - week ", week_2, "\n", age, " weeks old"),
|
|
main.title.size = 0.7)
|
|
|
|
# Add borders if requested
|
|
if (borders) {
|
|
map <- map +
|
|
tm_shape(pivot_shape) +
|
|
tm_borders(lwd = 3) +
|
|
tm_text("sub_field", size = 1/2) +
|
|
tm_shape(pivot_spans) +
|
|
tm_borders(lwd = 0.5, alpha = 0.5)
|
|
}
|
|
|
|
return(map)
|
|
}
|
|
|
|
#' Creates a visualization of CI data for a specific pivot field
|
|
#'
|
|
#' @param pivotName The name or ID of the pivot field to visualize
|
|
#' @param field_boundaries Field boundaries spatial data (sf object)
|
|
#' @param current_ci Current week's Chlorophyll Index raster
|
|
#' @param ci_minus_1 Previous week's Chlorophyll Index raster
|
|
#' @param ci_minus_2 Two weeks ago Chlorophyll Index raster
|
|
#' @param last_week_diff Difference raster between current and last week
|
|
#' @param three_week_diff Difference raster between current and three weeks ago
|
|
#' @param harvesting_data Data frame containing field harvesting/planting information
|
|
#' @param week Current week number
|
|
#' @param week_minus_1 Previous week number
|
|
#' @param week_minus_2 Two weeks ago week number
|
|
#' @param week_minus_3 Three weeks ago week number
|
|
#' @param borders Whether to display field borders (default: TRUE)
|
|
#' @return NULL (adds output directly to R Markdown document)
|
|
#'
|
|
ci_plot <- function(pivotName,
|
|
field_boundaries = AllPivots0,
|
|
current_ci = CI,
|
|
ci_minus_1 = CI_m1,
|
|
ci_minus_2 = CI_m2,
|
|
last_week_diff = last_week_dif_raster_abs,
|
|
three_week_diff = three_week_dif_raster_abs,
|
|
harvesting_data = harvesting_data,
|
|
week = week,
|
|
week_minus_1 = week_minus_1,
|
|
week_minus_2 = week_minus_2,
|
|
week_minus_3 = week_minus_3,
|
|
borders = TRUE){
|
|
# Input validation
|
|
if (missing(pivotName) || is.null(pivotName) || pivotName == "") {
|
|
stop("pivotName is required")
|
|
}
|
|
if (missing(field_boundaries) || is.null(field_boundaries)) {
|
|
stop("field_boundaries is required")
|
|
}
|
|
if (missing(current_ci) || is.null(current_ci)) {
|
|
stop("current_ci is required")
|
|
}
|
|
if (missing(ci_minus_1) || is.null(ci_minus_1)) {
|
|
stop("ci_minus_1 is required")
|
|
}
|
|
if (missing(ci_minus_2) || is.null(ci_minus_2)) {
|
|
stop("ci_minus_2 is required")
|
|
}
|
|
if (missing(last_week_diff) || is.null(last_week_diff)) {
|
|
stop("last_week_diff is required")
|
|
}
|
|
if (missing(three_week_diff) || is.null(three_week_diff)) {
|
|
stop("three_week_diff is required")
|
|
}
|
|
if (missing(harvesting_data) || is.null(harvesting_data)) {
|
|
stop("harvesting_data is required")
|
|
}
|
|
|
|
# Extract pivot shape and age data
|
|
tryCatch({
|
|
pivotShape <- field_boundaries %>% terra::subset(field %in% pivotName) %>% sf::st_transform(terra::crs(current_ci))
|
|
age <- harvesting_data %>%
|
|
dplyr::filter(field %in% pivotName) %>%
|
|
sort("year") %>%
|
|
tail(., 1) %>%
|
|
dplyr::select(age) %>%
|
|
unique() %>%
|
|
pull() %>%
|
|
round()
|
|
|
|
# Filter for the specific pivot
|
|
AllPivots2 <- field_boundaries %>% dplyr::filter(field %in% pivotName)
|
|
|
|
# Create crop masks for different timepoints using terra functions
|
|
singlePivot <- terra::crop(current_ci, pivotShape) %>% terra::mask(., pivotShape)
|
|
singlePivot_m1 <- terra::crop(ci_minus_1, pivotShape) %>% terra::mask(., pivotShape)
|
|
singlePivot_m2 <- terra::crop(ci_minus_2, pivotShape) %>% terra::mask(., pivotShape)
|
|
|
|
# Create difference maps
|
|
abs_CI_last_week <- terra::crop(last_week_diff, pivotShape) %>% terra::mask(., pivotShape)
|
|
abs_CI_three_week <- terra::crop(three_week_diff, pivotShape) %>% terra::mask(., pivotShape)
|
|
|
|
# Get planting date
|
|
planting_date <- harvesting_data %>%
|
|
dplyr::filter(field %in% pivotName) %>%
|
|
ungroup() %>%
|
|
dplyr::select(season_start) %>%
|
|
unique()
|
|
|
|
# Create spans for borders
|
|
joined_spans2 <- field_boundaries %>%
|
|
sf::st_transform(sf::st_crs(pivotShape)) %>% dplyr::filter(field %in% pivotName)
|
|
|
|
# Create the maps for different timepoints
|
|
CImap_m2 <- create_CI_map(singlePivot_m2, AllPivots2, joined_spans2,
|
|
show_legend = TRUE, legend_is_portrait = TRUE,
|
|
week = week_minus_2, age = age - 2, borders = borders)
|
|
|
|
CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2, joined_spans2,
|
|
show_legend = FALSE, legend_is_portrait = FALSE,
|
|
week = week_minus_1, age = age - 1, borders = borders)
|
|
|
|
CImap <- create_CI_map(singlePivot, AllPivots2, joined_spans2,
|
|
show_legend = FALSE, legend_is_portrait = FALSE,
|
|
week = week, age = age, borders = borders)
|
|
# Create difference maps - only show legend on the second one to avoid redundancy
|
|
CI_max_abs_last_week <- create_CI_diff_map(abs_CI_last_week, AllPivots2, joined_spans2,
|
|
show_legend = FALSE, legend_is_portrait = FALSE,
|
|
week_1 = week, week_2 = week_minus_1, age = age, borders = borders)
|
|
|
|
CI_max_abs_three_week <- create_CI_diff_map(abs_CI_three_week, AllPivots2, joined_spans2,
|
|
show_legend = TRUE, legend_is_portrait = TRUE,
|
|
week_1 = week, week_2 = week_minus_3, age = age, borders = borders)
|
|
# Arrange the maps with equal widths
|
|
tst <- tmap_arrange(CImap_m2, CImap_m1, CImap, CI_max_abs_last_week, CI_max_abs_three_week,
|
|
nrow = 1, widths = c(0.23, 0.18, 0.18, 0.18, 0.23))
|
|
|
|
# Output heading and map to R Markdown
|
|
cat(paste("## Field", pivotName, "-", age, "weeks after planting/harvest", "\n"))
|
|
print(tst)
|
|
|
|
}, error = function(e) {
|
|
safe_log(paste("Error creating CI plot for pivot", pivotName, ":", e$message), "ERROR")
|
|
cat(paste("## Field", pivotName, "- Error creating visualization", "\n"))
|
|
cat(paste("Error:", e$message, "\n"))
|
|
})
|
|
}
|
|
|
|
#' Creates a plot showing Chlorophyll Index data over time for a pivot field
|
|
#'
|
|
#' @param pivotName The name or ID of the pivot field to visualize
|
|
#' @param ci_quadrant_data Data frame containing CI quadrant data with field, sub_field, Date, DOY, cumulative_CI, value and season columns
|
|
#' @param plot_type Type of plot to generate: "value", "CI_rate", or "cumulative_CI"
|
|
#' @param facet_on Whether to facet the plot by season (TRUE) or overlay all seasons (FALSE)
|
|
#' @param x_unit Unit for x-axis: "days" for DOY or "weeks" for week number (default: "days")
|
|
#' @return NULL (adds output directly to R Markdown document)
|
|
#'
|
|
cum_ci_plot <- function(pivotName, ci_quadrant_data = CI_quadrant, plot_type = "value", facet_on = FALSE, x_unit = "days") {
|
|
# Input validation
|
|
if (missing(pivotName) || is.null(pivotName) || pivotName == "") {
|
|
stop("pivotName is required")
|
|
}
|
|
if (missing(ci_quadrant_data) || is.null(ci_quadrant_data)) {
|
|
stop("ci_quadrant_data is required")
|
|
}
|
|
if (!plot_type %in% c("value", "CI_rate", "cumulative_CI")) {
|
|
stop("plot_type must be one of: 'value', 'CI_rate', or 'cumulative_CI'")
|
|
}
|
|
if (!x_unit %in% c("days", "weeks")) {
|
|
stop("x_unit must be either 'days' or 'weeks'")
|
|
}
|
|
|
|
# Filter data for the specified pivot
|
|
tryCatch({
|
|
data_ci <- ci_quadrant_data %>% dplyr::filter(field == pivotName)
|
|
|
|
if (nrow(data_ci) == 0) {
|
|
safe_log(paste("No CI data found for pivot", pivotName), "WARNING")
|
|
return(cum_ci_plot2(pivotName)) # Use fallback function when no data is available
|
|
}
|
|
|
|
# Process data
|
|
data_ci2 <- data_ci %>%
|
|
dplyr::mutate(CI_rate = cumulative_CI / DOY,
|
|
week = lubridate::week(Date)) %>%
|
|
dplyr::group_by(field) %>%
|
|
dplyr::mutate(mean_CIrate_rolling_10_days = zoo::rollapplyr(CI_rate, width = 10, FUN = mean, partial = TRUE),
|
|
mean_rolling_10_days = zoo::rollapplyr(value, width = 10, FUN = mean, partial = TRUE))
|
|
|
|
data_ci2 <- data_ci2 %>% dplyr::mutate(season = as.factor(season))
|
|
|
|
# Prepare date information by season
|
|
date_preparation_perfect_pivot <- data_ci2 %>%
|
|
dplyr::group_by(season) %>%
|
|
dplyr::summarise(min_date = min(Date),
|
|
max_date = max(Date),
|
|
days = max_date - min_date)
|
|
|
|
# Get the 3 most recent seasons
|
|
unique_seasons <- sort(unique(date_preparation_perfect_pivot$season), decreasing = TRUE)[1:3]
|
|
|
|
# Determine the y aesthetic based on the plot type
|
|
y_aesthetic <- switch(plot_type,
|
|
"CI_rate" = "mean_CIrate_rolling_10_days",
|
|
"cumulative_CI" = "cumulative_CI",
|
|
"value" = "mean_rolling_10_days")
|
|
|
|
y_label <- switch(plot_type,
|
|
"CI_rate" = "10-Day Rolling Mean CI Rate (cumulative CI / age)",
|
|
"cumulative_CI" = "Cumulative CI",
|
|
"value" = "10-Day Rolling Mean CI")
|
|
|
|
# Determine x-axis variable based on x_unit parameter
|
|
x_var <- if (x_unit == "days") {
|
|
if (facet_on) "Date" else "DOY"
|
|
} else {
|
|
"week"
|
|
}
|
|
|
|
x_label <- switch(x_unit,
|
|
"days" = if (facet_on) "Date" else "Age of Crop (Days)",
|
|
"weeks" = "Week Number")
|
|
|
|
# Create plot with either facets by season or overlay by DOY/week
|
|
if (facet_on) {
|
|
g <- ggplot2::ggplot(data = data_ci2 %>% dplyr::filter(season %in% unique_seasons)) +
|
|
ggplot2::facet_wrap(~season, scales = "free_x") +
|
|
ggplot2::geom_line(ggplot2::aes_string(x = x_var, y = y_aesthetic, col = "sub_field", group = "sub_field")) +
|
|
ggplot2::labs(title = paste("Plot of", y_label, "for Pivot", pivotName),
|
|
color = "Field Name",
|
|
y = y_label,
|
|
x = x_label) +
|
|
ggplot2::scale_x_date(date_breaks = "1 month", date_labels = "%m-%Y") +
|
|
ggplot2::theme_minimal() +
|
|
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 60, hjust = 1),
|
|
legend.justification = c(1, 0), legend.position = c(1, 0),
|
|
legend.title = ggplot2::element_text(size = 8),
|
|
legend.text = ggplot2::element_text(size = 8)) +
|
|
ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = TRUE))
|
|
} else {
|
|
g <- ggplot2::ggplot(data = data_ci2 %>% dplyr::filter(season %in% unique_seasons)) +
|
|
ggplot2::geom_line(ggplot2::aes_string(x = x_var, y = y_aesthetic, col = "season", group = "season")) +
|
|
ggplot2::labs(title = paste("Plot of", y_label, "for Pivot", pivotName),
|
|
color = "Season",
|
|
y = y_label,
|
|
x = x_label) +
|
|
ggplot2::theme_minimal() +
|
|
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 60, hjust = 1),
|
|
legend.justification = c(1, 0), legend.position = c(1, 0),
|
|
legend.title = ggplot2::element_text(size = 8),
|
|
legend.text = ggplot2::element_text(size = 8)) +
|
|
ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = TRUE))
|
|
}
|
|
|
|
# Output plot to R Markdown with reduced height
|
|
subchunkify(g, 3.2, 10) # Reduced from 3.2 to 2.8
|
|
|
|
}, error = function(e) {
|
|
safe_log(paste("Error creating CI trend plot for pivot", pivotName, ":", e$message), "ERROR")
|
|
cum_ci_plot2(pivotName) # Use fallback function in case of error
|
|
})
|
|
}
|
|
|
|
#' Fallback function for creating CI visualization when data is missing
|
|
#'
|
|
#' @param pivotName The name or ID of the pivot field to visualize
|
|
#' @return NULL (adds output directly to R Markdown document)
|
|
#'
|
|
cum_ci_plot2 <- function(pivotName){
|
|
# Input validation
|
|
if (missing(pivotName) || is.null(pivotName) || pivotName == "") {
|
|
stop("pivotName is required")
|
|
}
|
|
|
|
# Create a simple plot showing "No data available"
|
|
tryCatch({
|
|
end_date <- Sys.Date()
|
|
start_date <- end_date %m-% months(11) # 11 months ago from end_date
|
|
date_seq <- seq.Date(from = start_date, to = end_date, by = "month")
|
|
midpoint_date <- start_date + (end_date - start_date) / 2
|
|
|
|
g <- ggplot() +
|
|
scale_x_date(limits = c(start_date, end_date), date_breaks = "1 month", date_labels = "%m-%Y") +
|
|
scale_y_continuous(limits = c(0, 4)) +
|
|
labs(title = paste("14 day rolling MEAN CI rate - Field ", pivotName),
|
|
x = "Date", y = "CI Rate") +
|
|
theme_minimal() +
|
|
theme(axis.text.x = element_text(angle = 60, hjust = 1),
|
|
legend.justification = c(1, 0), legend.position = c(1, 0),
|
|
legend.title = element_text(size = 8),
|
|
legend.text = element_text(size = 8)) +
|
|
annotate("text", x = midpoint_date, y = 2, label = "No data available", size = 6, hjust = 0.5)
|
|
|
|
subchunkify(g, 3.2, 10)
|
|
|
|
}, error = function(e) {
|
|
safe_log(paste("Error creating fallback CI plot for pivot", pivotName, ":", e$message), "ERROR")
|
|
cat(paste("No data available for field", pivotName, "\n"))
|
|
})
|
|
}
|
|
|
|
#' Gets the file path for a specific week's mosaic
|
|
#'
|
|
#' @param mosaic_path Base directory containing mosaic files
|
|
#' @param input_date Reference date to calculate from
|
|
#' @param week_offset Number of weeks to offset from input date (positive or negative)
|
|
#' @return File path to the requested week's mosaic TIF file
|
|
#'
|
|
get_week_path <- function(mosaic_path, input_date, week_offset) {
|
|
# Input validation
|
|
if (missing(mosaic_path) || is.null(mosaic_path) || mosaic_path == "") {
|
|
stop("mosaic_path is required")
|
|
}
|
|
if (missing(input_date)) {
|
|
stop("input_date is required")
|
|
}
|
|
|
|
tryCatch({
|
|
# Convert input_date to Date object (in case it's a string)
|
|
input_date <- as.Date(input_date)
|
|
if (is.na(input_date)) {
|
|
stop("Invalid input_date. Expected a Date object or a string convertible to Date.")
|
|
}
|
|
|
|
# Validate week_offset
|
|
week_offset <- as.integer(week_offset)
|
|
if (is.na(week_offset)) {
|
|
stop("Invalid week_offset. Expected an integer value.")
|
|
}
|
|
|
|
# Get the start of the week for the input date (adjust to Monday as the start of the week)
|
|
start_of_week <- lubridate::floor_date(input_date, unit = "week", week_start = 1)
|
|
|
|
# Calculate the new date after applying the week offset
|
|
target_date <- start_of_week + lubridate::weeks(week_offset)
|
|
|
|
# Get the week number and year of the target date
|
|
target_week <- sprintf("%02d", lubridate::isoweek(target_date)) # Left-pad week number with a zero if needed
|
|
target_year <- lubridate::isoyear(target_date)
|
|
|
|
# Generate the file path for the target week
|
|
path_to_week <- here::here(mosaic_path, paste0("week_", target_week, "_", target_year, ".tif"))
|
|
|
|
# Log the path calculation
|
|
safe_log(paste("Calculated path for week", target_week, "of year", target_year, ":", path_to_week), "INFO")
|
|
|
|
# Return the path
|
|
return(path_to_week)
|
|
|
|
}, error = function(e) {
|
|
safe_log(paste("Error calculating week path:", e$message), "ERROR")
|
|
stop(e$message)
|
|
})
|
|
}
|
|
|