Timon
07aee7bed1
- Replace raster package with terra throughout the codebase - Update map visualizations with better layout and legends - Add descriptive headers to report sections - Improve map legend positioning and sizing - Enhance error handling for missing data - Remove redundant legends in field-specific visualizations - Optimize figure dimensions to prevent page overflow - Expand documentation of CI index and report components - Update package dependencies in packages.
84 lines
2.8 KiB
R
84 lines
2.8 KiB
R
library(ggplot2)
|
|
library(dplyr)
|
|
library(sf)
|
|
library(terra)
|
|
|
|
# Define the file paths
|
|
raster_dir <- "C:/Users/timon/Resilience BV/4020 SCane ESA DEMO - Documenten/General/4020 SCDEMO Team/4020 TechnicalData/WP3/smartcane/laravel_app/storage/app/chemba/merged_final_tif"
|
|
polygon_path <- "C:/Users/timon/Resilience BV/4020 SCane ESA DEMO - Documenten/General/4020 SCDEMO Team/4020 TechnicalData/WP3/smartcane/r_app/experiments/pivot.geojson"
|
|
|
|
# Load the polygon
|
|
polygon <- st_read(polygon_path)
|
|
|
|
# Filter the polygon for the specific field
|
|
polygon_filtered <- subset(polygon, field == "1.1")
|
|
|
|
# List all raster files in the directory
|
|
raster_files <- list.files(raster_dir, pattern = "\\.tif$", full.names = FALSE)
|
|
raster_file <- raster_files[1]
|
|
|
|
# Initialize an empty list to store CI values
|
|
ci_values <- list()
|
|
|
|
# Define a function to process each raster file
|
|
process_raster_file <- function(raster_file) {
|
|
# Extract the date from the file name (assuming the date is in YYYYMMDD format)
|
|
date <- as.Date(substr(raster_file, 1, 10), format = "%Y-%m-%d")
|
|
|
|
# Load the raster
|
|
raster_data <- terra::rast(file.path(raster_dir, raster_file))
|
|
|
|
# Crop the raster using the filtered polygon
|
|
cropped_raster <- terra::crop(raster_data, polygon_filtered)
|
|
|
|
# Extract the CI band (assuming the CI band is the first band)
|
|
ci_band <- cropped_raster$CI
|
|
|
|
# Extract the CI values from the CI band
|
|
return(data.frame(Date = date, CI = terra::values(ci_band)) %>%
|
|
filter(!is.na(CI)) %>%
|
|
mutate(Polygon = unique(polygon_filtered$field)))
|
|
}
|
|
|
|
# Use walk to apply the function to each raster file
|
|
ci_values <- purrr::map(raster_files, process_raster_file)
|
|
|
|
# Combine all CI values into a single data frame
|
|
ci_values_df <- do.call(rbind, ci_values)
|
|
|
|
head(ci_values_df)
|
|
|
|
# Plot the first raster
|
|
plot(terra::rast(file.path(raster_dir, raster_files[1]))[[1]])
|
|
|
|
# Calculate the mean CI value
|
|
ci_stats <- ci_values_df %>%
|
|
group_by(Date, Polygon) %>%
|
|
summarize(
|
|
mean_ci = mean(CI, na.rm = TRUE),
|
|
q5 = quantile(CI, 0.05, na.rm = TRUE),
|
|
q25 = quantile(CI, 0.25, na.rm = TRUE),
|
|
q50 = quantile(CI, 0.50, na.rm = TRUE),
|
|
q75 = quantile(CI, 0.75, na.rm = TRUE),
|
|
q95 = quantile(CI, 0.95, na.rm = TRUE),
|
|
cv = sd(CI, na.rm = TRUE) / mean(CI, na.rm = TRUE) * 100,
|
|
.groups = "drop"
|
|
)
|
|
|
|
# Plot the mean CI value over time
|
|
ggplot(ci_stats, aes(x = Date, y = mean_ci)) +
|
|
geom_line() +
|
|
geom_ribbon(aes(ymin = q25, ymax = q75), alpha = 0.2) +
|
|
labs(title = "Mean CI value over time",
|
|
x = "Date",
|
|
y = "Mean CI") +
|
|
theme_minimal()
|
|
|
|
# Plot the coefficient of variation over time
|
|
ggplot(ci_stats, aes(x = Date, y = cv)) +
|
|
geom_line() +
|
|
labs(title = "Coefficient of variation over time",
|
|
x = "Date",
|
|
y = "CV (%)") +
|
|
theme_minimal()
|