88 lines
3.1 KiB
R
88 lines
3.1 KiB
R
library(raster)
|
|
library(rgdal)
|
|
library(ggplot2)
|
|
library(dplyr)
|
|
|
|
# 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
|
|
library(sf)
|
|
|
|
# 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
|
|
library(terra)
|
|
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 <- rast(file.path(raster_dir, raster_file))
|
|
|
|
# Crop the raster using the filtered polygon
|
|
cropped_raster <- 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 = 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(rast(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
|
|
)
|
|
|
|
# Plot the coefficient of variation over time
|
|
ggplot(ci_stats, aes(x = Date, y = cv, group = Polygon)) +
|
|
geom_line(color = "green") +
|
|
labs(title = "Coefficient of Variation (CV) Over Time",
|
|
x = "Date",
|
|
y = "Coefficient of Variation (%)") +
|
|
theme_minimal()
|
|
|
|
# Plot the CI statistics over time with different quantiles and coefficient of variation
|
|
ggplot(ci_stats, aes(x = Date, group = Polygon)) +
|
|
geom_line(aes(y = q50), color = "red") +
|
|
geom_ribbon(aes(ymin = q25, ymax = q75), fill = "blue", alpha = 0.2) +
|
|
geom_ribbon(aes(ymin = q5, ymax = q95), fill = "blue", alpha = 0.1) +
|
|
labs(title = "CI Statistics Over Time",
|
|
x = "Date",
|
|
y = "CI Value / Coefficient of Variation (%)") +
|
|
theme_minimal()
|